Line data Source code
1 : !=======================================================================
2 :
3 : ! Aerosol tracer within sea ice
4 : !
5 : ! authors Marika Holland, NCAR
6 : ! David Bailey, NCAR
7 :
8 : module icepack_aerosol
9 :
10 : use icepack_kinds
11 : use icepack_parameters, only: c0, c1, c2, puny, rhoi, rhos, hs_min
12 : use icepack_parameters, only: hi_ssl, hs_ssl
13 : use icepack_tracers, only: max_aero
14 : use icepack_warnings, only: warnstr, icepack_warnings_add
15 : use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
16 :
17 : use icepack_zbgc_shared, only: kscavz
18 :
19 : implicit none
20 :
21 : private
22 : public :: update_aerosol, update_snow_bgc
23 :
24 : !=======================================================================
25 :
26 : contains
27 :
28 : !=======================================================================
29 :
30 : ! Increase aerosol in ice or snow surface due to deposition
31 : ! and vertical cycling
32 :
33 52109425 : subroutine update_aerosol(dt, &
34 : nilyr, nslyr, &
35 : n_aero, &
36 : meltt, melts, &
37 : meltb, congel, &
38 : snoice, &
39 : fsnow, &
40 52109425 : aerosno, aeroice, &
41 : aice_old, &
42 : vice_old, vsno_old, &
43 : vicen, vsnon, aicen, &
44 104218850 : faero_atm, faero_ocn)
45 :
46 : integer (kind=int_kind), intent(in) :: &
47 : nilyr, nslyr, n_aero
48 :
49 : real (kind=dbl_kind), intent(in) :: &
50 : dt, & ! time step
51 : meltt, & ! thermodynamic melt/growth rates
52 : melts, &
53 : meltb, &
54 : congel, &
55 : snoice, &
56 : fsnow, &
57 : vicen, & ! ice volume (m)
58 : vsnon, & ! snow volume (m)
59 : aicen, & ! ice area fraction
60 : aice_old, & ! values prior to thermodynamic changes
61 : vice_old, &
62 : vsno_old
63 :
64 : real (kind=dbl_kind), dimension(:), &
65 : intent(in) :: &
66 : faero_atm ! aerosol deposition rate (W/m^2 s)
67 :
68 : real (kind=dbl_kind), dimension(:), &
69 : intent(inout) :: &
70 : faero_ocn ! aerosol flux to ocean (W/m^2 s)
71 :
72 : real (kind=dbl_kind), dimension(:,:), intent(inout) :: &
73 : aerosno, aeroice ! kg/m^2
74 :
75 : ! local variables
76 : integer (kind=int_kind) :: k
77 :
78 : real (kind=dbl_kind) :: &
79 3952149 : dzssl, dzssl_new, & ! snow ssl thickness
80 3952149 : dzint, & ! snow interior thickness
81 3952149 : dzssli, dzssli_new, & ! ice ssl thickness
82 3952149 : dzinti, & ! ice interior thickness
83 3952149 : dznew, & ! tracks thickness changes
84 7904298 : hs, hi, & ! snow/ice thickness (m)
85 7904298 : dhs_evap, dhi_evap, & ! snow/ice thickness change due to evap
86 7904298 : dhs_melts, dhi_meltt, & ! ... due to surface melt
87 7904298 : dhs_snoice, dhi_snoice, & ! ... due to snow-ice formation
88 7904298 : dhi_congel, dhi_meltb, & ! ... due to bottom growth, melt
89 7904298 : hslyr, hilyr, & ! snow, ice layer thickness (m)
90 7904298 : hslyr_old, hilyr_old, & ! old snow, ice layer thickness (m)
91 7904298 : hs_old, hi_old, & ! old snow, ice thickness (m)
92 3952149 : sloss1, sloss2, & ! aerosol mass loss (kg/m^2)
93 3952149 : ar ! 1/aicen(i,j)
94 :
95 : real (kind=dbl_kind), dimension(max_aero) :: &
96 : kscav, kscavsi ! scavenging by melt water
97 :
98 : real (kind=dbl_kind), dimension(n_aero) :: &
99 108170999 : aerotot, aerotot0, & ! for conservation check
100 104218850 : focn_old ! for conservation check
101 :
102 : ! echmod: this assumes max_aero=6
103 : data kscav / .03_dbl_kind, .20_dbl_kind, .02_dbl_kind, &
104 : .02_dbl_kind, .01_dbl_kind, .01_dbl_kind /
105 : data kscavsi / .03_dbl_kind, .20_dbl_kind, .02_dbl_kind, &
106 : .02_dbl_kind, .01_dbl_kind, .01_dbl_kind /
107 :
108 : character(len=*),parameter :: subname='(update_aerosol)'
109 :
110 : !-------------------------------------------------------------------
111 : ! initialize
112 : !-------------------------------------------------------------------
113 104218850 : focn_old(:) = faero_ocn(:)
114 :
115 52109425 : hs_old = vsno_old/aice_old
116 52109425 : hi_old = vice_old/aice_old
117 52109425 : hslyr_old = hs_old/real(nslyr,kind=dbl_kind)
118 52109425 : hilyr_old = hi_old/real(nilyr,kind=dbl_kind)
119 :
120 52109425 : dzssl = min(hslyr_old/c2, hs_ssl)
121 52109425 : dzssli = min(hilyr_old/c2, hi_ssl)
122 52109425 : dzint = hs_old - dzssl
123 52109425 : dzinti = hi_old - dzssli
124 :
125 52109425 : if (aicen > c0) then
126 52109425 : ar = c1/aicen
127 52109425 : hs = vsnon*ar
128 52109425 : hi = vicen*ar
129 52109425 : dhs_melts = -melts*ar
130 52109425 : dhi_snoice = snoice*ar
131 52109425 : dhs_snoice = dhi_snoice*rhoi/rhos
132 52109425 : dhi_meltt = -meltt*ar
133 52109425 : dhi_meltb = -meltb*ar
134 52109425 : dhi_congel = congel*ar
135 : else ! ice disappeared during time step
136 0 : hs = vsnon/aice_old
137 0 : hi = vicen/aice_old
138 0 : dhs_melts = -melts/aice_old
139 0 : dhi_snoice = snoice/aice_old
140 0 : dhs_snoice = dhi_snoice*rhoi/rhos
141 0 : dhi_meltt = -meltt/aice_old
142 0 : dhi_meltb = -meltb/aice_old
143 0 : dhi_congel = congel/aice_old
144 : endif
145 :
146 : dhs_evap = hs - (hs_old + dhs_melts - dhs_snoice &
147 52109425 : + fsnow/rhos*dt)
148 : dhi_evap = hi - (hi_old + dhi_meltt + dhi_meltb &
149 52109425 : + dhi_congel + dhi_snoice)
150 :
151 104218850 : do k = 1, n_aero
152 3952149 : aerotot0(k) = aerosno(k,2) + aerosno(k,1) &
153 104218850 : + aeroice(k,2) + aeroice(k,1)
154 : enddo
155 :
156 : !-------------------------------------------------------------------
157 : ! evaporation
158 : !-------------------------------------------------------------------
159 52109425 : dzint = dzint + min(dzssl + dhs_evap, c0)
160 52109425 : dzinti = dzinti + min(dzssli + dhi_evap, c0)
161 52109425 : dzssl = max(dzssl + dhs_evap, c0)
162 52109425 : dzssli = max(dzssli + dhi_evap, c0)
163 :
164 : !-------------------------------------------------------------------
165 : ! basal ice growth
166 : !-------------------------------------------------------------------
167 52109425 : dzinti = dzinti + dhi_congel
168 :
169 : !-------------------------------------------------------------------
170 : ! surface snow melt
171 : !-------------------------------------------------------------------
172 52109425 : if (-dhs_melts > puny) then
173 7316724 : do k = 1, n_aero
174 3658362 : sloss1 = c0
175 3658362 : sloss2 = c0
176 3658362 : if (dzssl > puny) &
177 14836 : sloss1 = kscav(k)*aerosno(k,1) &
178 3657867 : *min(-dhs_melts,dzssl)/dzssl
179 3658362 : aerosno(k,1) = aerosno(k,1) - sloss1
180 3658362 : if (dzint > puny) &
181 14838 : sloss2 = kscav(k)*aerosno(k,2) &
182 3658058 : *max(-dhs_melts-dzssl,c0)/dzint
183 3658362 : aerosno(k,2) = aerosno(k,2) - sloss2
184 7316724 : faero_ocn(k) = faero_ocn(k) + (sloss1+sloss2)/dt
185 : enddo ! n_aero
186 :
187 : ! update snow thickness
188 3658362 : dzint=dzint+min(dzssl+dhs_melts, c0)
189 3658362 : dzssl=max(dzssl+dhs_melts, c0)
190 :
191 3658362 : if ( dzssl <= puny ) then ! ssl melts away
192 1100956 : aerosno(:,2) = aerosno(:,1) + aerosno(:,2)
193 1100956 : aerosno(:,1) = c0
194 550478 : dzssl = max(dzssl, c0)
195 : endif
196 3658362 : if (dzint <= puny ) then ! all snow melts away
197 0 : aeroice(:,1) = aeroice(:,1) &
198 1017598 : + aerosno(:,1) + aerosno(:,2)
199 2543995 : aerosno(:,:) = c0
200 508799 : dzint = max(dzint, c0)
201 : endif
202 : endif
203 :
204 : !-------------------------------------------------------------------
205 : ! surface ice melt
206 : !-------------------------------------------------------------------
207 52109425 : if (-dhi_meltt > puny) then
208 3424344 : do k = 1, n_aero
209 1712172 : sloss1 = c0
210 1712172 : sloss2 = c0
211 1712172 : if (dzssli > puny) &
212 359 : sloss1 = kscav(k)*aeroice(k,1) &
213 1711803 : *min(-dhi_meltt,dzssli)/dzssli
214 1712172 : aeroice(k,1) = aeroice(k,1) - sloss1
215 1712172 : if (dzinti > puny) &
216 378 : sloss2 = kscav(k)*aeroice(k,2) &
217 1712172 : *max(-dhi_meltt-dzssli,c0)/dzinti
218 1712172 : aeroice(k,2) = aeroice(k,2) - sloss2
219 3424344 : faero_ocn(k) = faero_ocn(k) + (sloss1+sloss2)/dt
220 : enddo
221 :
222 1712172 : dzinti = dzinti + min(dzssli+dhi_meltt, c0)
223 1712172 : dzssli = max(dzssli+dhi_meltt, c0)
224 1712172 : if (dzssli <= puny) then ! ssl ice melts away
225 1962 : do k = 1, n_aero
226 981 : aeroice(k,2) = aeroice(k,1) + aeroice(k,2)
227 1962 : aeroice(k,1) = c0
228 : enddo
229 981 : dzssli = max(dzssli, c0)
230 : endif
231 1712172 : if (dzinti <= puny) then ! all ice melts away
232 0 : do k = 1, n_aero
233 0 : faero_ocn(k) = faero_ocn(k) &
234 0 : + (aeroice(k,1)+aeroice(k,2))/dt
235 0 : aeroice(k,:)=c0
236 : enddo
237 0 : dzinti = max(dzinti, c0)
238 : endif
239 : endif
240 :
241 : !-------------------------------------------------------------------
242 : ! basal ice melt. Assume all aero lost in basal melt
243 : !-------------------------------------------------------------------
244 52109425 : if (-dhi_meltb > puny) then
245 28633562 : do k=1,n_aero
246 14316781 : sloss1=c0
247 14316781 : sloss2=c0
248 14316781 : if (dzssli > puny) &
249 : sloss1 = max(-dhi_meltb-dzinti, c0) &
250 14315664 : *aeroice(k,1)/dzssli
251 14316781 : aeroice(k,1) = aeroice(k,1) - sloss1
252 14316781 : if (dzinti > puny) &
253 : sloss2 = min(-dhi_meltb, dzinti) &
254 14315918 : *aeroice(k,2)/dzinti
255 14316781 : aeroice(k,2) = aeroice(k,2) - sloss2
256 28633562 : faero_ocn(k) = faero_ocn(k) + (sloss1+sloss2)/dt
257 : enddo
258 :
259 14316781 : dzssli = dzssli + min(dzinti+dhi_meltb, c0)
260 14316781 : dzinti = max(dzinti+dhi_meltb, c0)
261 : endif
262 :
263 : !-------------------------------------------------------------------
264 : ! snowfall
265 : !-------------------------------------------------------------------
266 52109425 : if (fsnow > c0) dzssl = dzssl + fsnow/rhos*dt
267 :
268 : !-------------------------------------------------------------------
269 : ! snow-ice formation
270 : !-------------------------------------------------------------------
271 52109425 : if (dhs_snoice > puny) then
272 2862578 : do k = 1, n_aero
273 1431289 : sloss1 = c0
274 1431289 : sloss2 = c0
275 1431289 : if (dzint > puny) &
276 : sloss2 = min(dhs_snoice, dzint) &
277 1431289 : *aerosno(k,2)/dzint
278 1431289 : aerosno(k,2) = aerosno(k,2) - sloss2
279 1431289 : if (dzssl > puny) &
280 : sloss1 = max(dhs_snoice-dzint, c0) &
281 1431289 : *aerosno(k,1)/dzssl
282 1431289 : aerosno(k,1) = aerosno(k,1) - sloss1
283 178248 : aeroice(k,1) = aeroice(k,1) &
284 1431289 : + (c1-kscavsi(k))*(sloss2+sloss1)
285 178248 : faero_ocn(k) = faero_ocn(k) &
286 2862578 : + kscavsi(k)*(sloss2+sloss1)/dt
287 : enddo
288 1431289 : dzssl = dzssl - max(dhs_snoice-dzint, c0)
289 1431289 : dzint = max(dzint-dhs_snoice, c0)
290 1431289 : dzssli = dzssli + dhi_snoice
291 : endif
292 :
293 : !-------------------------------------------------------------------
294 : ! aerosol deposition
295 : !-------------------------------------------------------------------
296 52109425 : if (aicen > c0) then
297 52109425 : hs = vsnon * ar
298 : else
299 0 : hs = c0
300 : endif
301 52109425 : if (hs > hs_min) then ! should this really be hs_min or 0?
302 : ! should use same hs_min value as in radiation
303 100464474 : do k=1,n_aero
304 3938636 : aerosno(k,1) = aerosno(k,1) &
305 100464474 : + faero_atm(k)*dt*aicen
306 : enddo
307 : else
308 3754376 : do k=1,n_aero
309 13513 : aeroice(k,1) = aeroice(k,1) &
310 3754376 : + faero_atm(k)*dt*aicen
311 : enddo
312 : endif
313 :
314 : !-------------------------------------------------------------------
315 : ! redistribute aerosol within vertical layers
316 : !-------------------------------------------------------------------
317 52109425 : if (aicen > c0) then
318 52109425 : hs = vsnon * ar ! new snow thickness
319 52109425 : hi = vicen * ar ! new ice thickness
320 : else
321 0 : hs = c0
322 0 : hi = c0
323 : endif
324 52109425 : if (dzssl <= puny) then ! nothing in SSL
325 3624234 : do k=1,n_aero
326 1812117 : aerosno(k,2) = aerosno(k,2) + aerosno(k,1)
327 3624234 : aerosno(k,1) = c0
328 : enddo
329 : endif
330 52109425 : if (dzint <= puny) then ! nothing in Snow Int
331 5825724 : do k = 1, n_aero
332 2912862 : aeroice(k,1) = aeroice(k,1) + aerosno(k,2)
333 5825724 : aerosno(k,2) = c0
334 : enddo
335 : endif
336 52109425 : if (dzssli <= puny) then ! nothing in Ice SSL
337 14225800 : do k = 1, n_aero
338 7112900 : aeroice(k,2) = aeroice(k,2) + aeroice(k,1)
339 14225800 : aeroice(k,1) = c0
340 : enddo
341 : endif
342 :
343 52109425 : if (dzinti <= puny) then ! nothing in Ice INT
344 12314730 : do k = 1, n_aero
345 466902 : faero_ocn(k) = faero_ocn(k) &
346 6157365 : + (aeroice(k,1)+aeroice(k,2))/dt
347 24629460 : aeroice(k,:)=c0
348 : enddo
349 : endif
350 :
351 52109425 : hslyr = hs/real(nslyr,kind=dbl_kind)
352 52109425 : hilyr = hi/real(nilyr,kind=dbl_kind)
353 52109425 : dzssl_new = min(hslyr/c2, hs_ssl)
354 52109425 : dzssli_new = min(hilyr/c2, hi_ssl)
355 :
356 52109425 : if (hs > hs_min) then
357 100464474 : do k = 1, n_aero
358 50232237 : dznew = min(dzssl_new-dzssl, c0)
359 50232237 : sloss1 = c0
360 50232237 : if (dzssl > puny) &
361 50201138 : sloss1 = dznew*aerosno(k,1)/dzssl ! not neccesarily a loss
362 50232237 : dznew = max(dzssl_new-dzssl, c0)
363 50232237 : if (dzint > puny) &
364 49110445 : sloss1 = sloss1 + aerosno(k,2)*dznew/dzint
365 50232237 : aerosno(k,1) = aerosno(k,1) + sloss1
366 100464474 : aerosno(k,2) = aerosno(k,2) - sloss1
367 : enddo
368 : else
369 0 : aeroice(:,1) = aeroice(:,1) &
370 3754376 : + aerosno(:,1) + aerosno(:,2)
371 9385940 : aerosno(:,:) = c0
372 : endif
373 :
374 52109425 : if (vicen > puny) then ! may want a limit on hi instead?
375 104167618 : do k = 1, n_aero
376 52083809 : sloss2 = c0
377 52083809 : dznew = min(dzssli_new-dzssli, c0)
378 52083809 : if (dzssli > puny) &
379 44978425 : sloss2 = dznew*aeroice(k,1)/dzssli
380 52083809 : dznew = max(dzssli_new-dzssli, c0)
381 52083809 : if (dzinti > puny) &
382 45944544 : sloss2 = sloss2 + aeroice(k,2)*dznew/dzinti
383 52083809 : aeroice(k,1) = aeroice(k,1) + sloss2
384 104167618 : aeroice(k,2) = aeroice(k,2) - sloss2
385 : enddo
386 : else
387 51232 : faero_ocn(:) = faero_ocn(:) + (aeroice(:,1)+aeroice(:,2))/dt
388 128080 : aeroice(:,:) = c0
389 : endif
390 :
391 : !-------------------------------------------------------------------
392 : ! check conservation
393 : !-------------------------------------------------------------------
394 104218850 : do k = 1, n_aero
395 3952149 : aerotot(k) = aerosno(k,2) + aerosno(k,1) &
396 52109425 : + aeroice(k,2) + aeroice(k,1)
397 56061574 : if ((aerotot(k)-aerotot0(k)) &
398 3952149 : - ( faero_atm(k)*aicen &
399 56061574 : - (faero_ocn(k)-focn_old(k)) )*dt > puny) then
400 :
401 0 : write(warnstr,*) subname, 'aerosol tracer: ',k
402 0 : call icepack_warnings_add(warnstr)
403 0 : write(warnstr,*) subname, 'aerotot-aerotot0 ',aerotot(k)-aerotot0(k)
404 0 : call icepack_warnings_add(warnstr)
405 0 : write(warnstr,*) subname, 'faero_atm-faero_ocn ', &
406 0 : (faero_atm(k)*aicen-(faero_ocn(k)-focn_old(k)))*dt
407 0 : call icepack_warnings_add(warnstr)
408 : endif
409 : enddo
410 :
411 : !-------------------------------------------------------------------
412 : ! check for negative values
413 : !-------------------------------------------------------------------
414 :
415 : !echmod: note that this does not test or fix all aero tracers
416 3952149 : if (aeroice(1,1) < -puny .or. &
417 3952149 : aeroice(1,2) < -puny .or. &
418 52109425 : aerosno(1,1) < -puny .or. &
419 : aerosno(1,2) < -puny) then
420 :
421 0 : write(warnstr,*) subname, 'aerosol negative in aerosol code'
422 0 : call icepack_warnings_add(warnstr)
423 :
424 0 : aeroice(1,1) = max(aeroice(1,1), c0)
425 0 : aeroice(1,2) = max(aeroice(1,2), c0)
426 0 : aerosno(1,1) = max(aerosno(1,1), c0)
427 0 : aerosno(1,2) = max(aerosno(1,2), c0)
428 :
429 : endif
430 :
431 52109425 : end subroutine update_aerosol
432 :
433 : !=======================================================================
434 :
435 : ! Increase aerosol in snow surface due to deposition
436 : ! and vertical cycling : after update_aerosol
437 :
438 29415679 : subroutine update_snow_bgc (dt, nblyr, &
439 : nslyr, &
440 : meltt, melts, &
441 : meltb, congel, &
442 : snoice, nbtrcr, &
443 : fsnow, ntrcr, &
444 29415679 : trcrn, bio_index, &
445 29415679 : aice_old, zbgc_snow, &
446 : vice_old, vsno_old, &
447 : vicen, vsnon, &
448 29415679 : aicen, flux_bio_atm,&
449 29415679 : zbgc_atm, flux_bio)
450 :
451 : integer (kind=int_kind), intent(in) :: &
452 : nbtrcr, & ! number of distinct snow tracers
453 : nblyr, & ! number of bio layers
454 : nslyr, & ! number of snow layers
455 : ntrcr ! number of tracers
456 :
457 : integer (kind=int_kind), dimension (nbtrcr), intent(in) :: &
458 : bio_index
459 :
460 : real (kind=dbl_kind), intent(in) :: &
461 : dt ! time step
462 :
463 : real (kind=dbl_kind), intent(in) :: &
464 : meltt, & ! thermodynamic melt/growth rates
465 : melts, &
466 : meltb, &
467 : congel, &
468 : snoice, &
469 : fsnow, &
470 : vicen, & ! ice volume (m)
471 : vsnon, & ! snow volume (m)
472 : aicen, & ! ice area fraction
473 : aice_old, & ! values prior to thermodynamic changes
474 : vice_old, &
475 : vsno_old
476 :
477 : real (kind=dbl_kind),dimension(nbtrcr), intent(inout) :: &
478 : zbgc_snow, & ! aerosol contribution from snow to ice
479 : zbgc_atm, & ! and atm to ice concentration * volume (kg or mmol/m^3*m)
480 : flux_bio ! total ocean tracer flux (mmol/m^2/s)
481 :
482 : real (kind=dbl_kind), dimension(nbtrcr), &
483 : intent(in) :: &
484 : flux_bio_atm ! aerosol deposition rate (kg or mmol/m^2 s)
485 :
486 : real (kind=dbl_kind), dimension(ntrcr), &
487 : intent(inout) :: &
488 : trcrn ! ice/snow tracer array
489 :
490 : ! local variables
491 :
492 : integer (kind=int_kind) :: k, n
493 :
494 : real (kind=dbl_kind) :: &
495 80291 : dzssl, dzssl_new, & ! snow ssl thickness
496 80291 : dzint, dzint_new, & ! snow interior thickness
497 80291 : hs, & ! snow thickness (m)
498 80291 : dhs_evap, & ! snow thickness change due to evap
499 80291 : dhs_melts, & ! ... due to surface melt
500 80291 : dhs_snoice, & ! ... due to snow-ice formation
501 80291 : hslyr, & ! snow layer thickness (m)
502 80291 : hslyr_old, & ! old snow layer thickness (m)
503 80291 : hs_old, & ! old snow thickness (m)
504 80291 : dznew, & ! change in the snow sl (m)
505 80291 : sloss1, sloss2, & ! aerosol mass loss (kg/m^2)
506 80291 : ar ! 1/aicen(i,j)
507 :
508 : real (kind=dbl_kind), dimension(nbtrcr) :: &
509 61802125 : aerotot, aerotot0, & ! for conservation check (mmol/m^3)
510 60276596 : aero_cons , & ! for conservation check (mmol/m^2)
511 60276596 : flux_bio_o ! initial ocean tracer flux (mmol/m^2/s)
512 :
513 : real (kind=dbl_kind), dimension(nbtrcr,2) :: &
514 61962707 : aerosno, & ! kg/m^2
515 62042998 : aerosno0 ! for diagnostic prints
516 :
517 : character(len=*),parameter :: subname='(update_snow_bgc)'
518 :
519 : !-------------------------------------------------------------------
520 : ! initialize
521 : !-------------------------------------------------------------------
522 1206042839 : aerosno (:,:) = c0
523 1206042839 : aerosno0(:,:) = c0
524 588313580 : aero_cons(:) = c0
525 588313580 : zbgc_snow(:) = c0
526 588313580 : zbgc_atm(:) = c0
527 :
528 29415679 : hs_old = vsno_old/aice_old
529 29415679 : hslyr_old = hs_old/real(nslyr,kind=dbl_kind)
530 :
531 29415679 : dzssl = min(hslyr_old/c2, hs_ssl)
532 29415679 : dzint = hs_old - dzssl
533 :
534 29415679 : if (aicen > c0) then
535 29415679 : ar = c1/aicen
536 29415679 : hs = vsnon*ar
537 29415679 : dhs_melts = -melts
538 29415679 : dhs_snoice = snoice*rhoi/rhos
539 : else ! ice disappeared during time step
540 0 : ar = c1
541 0 : hs = vsnon/aice_old
542 0 : dhs_melts = -melts
543 0 : dhs_snoice = snoice*rhoi/rhos
544 : endif
545 :
546 : dhs_evap = hs - (hs_old + dhs_melts - dhs_snoice &
547 29415679 : + fsnow/rhos*dt)
548 :
549 : ! trcrn() has units kg/m^3
550 :
551 29415679 : if ((vsno_old .le. puny) .or. (vsnon .le. puny)) then
552 :
553 24968601 : do k=1,nbtrcr
554 152817 : flux_bio(k) = flux_bio(k) + &
555 305634 : (trcrn(bio_index(k)+ nblyr+1)*dzssl+ &
556 22743456 : trcrn(bio_index(k)+ nblyr+2)*dzint)/dt
557 22590639 : trcrn(bio_index(k) + nblyr+1) = c0
558 22590639 : trcrn(bio_index(k) + nblyr+2) = c0
559 152817 : zbgc_atm(k) = zbgc_atm(k) &
560 23779620 : + flux_bio_atm(k)*dt
561 : enddo
562 :
563 : else
564 :
565 564533960 : do k=1,nbtrcr
566 536307262 : flux_bio_o(k) = flux_bio(k)
567 536307262 : aerosno (k,1) = trcrn(bio_index(k)+ nblyr+1) * dzssl
568 536307262 : aerosno (k,2) = trcrn(bio_index(k)+ nblyr+2) * dzint
569 1608921786 : aerosno0(k,:) = aerosno(k,:)
570 564533960 : aerotot0(k) = aerosno(k,2) + aerosno(k,1)
571 : enddo
572 :
573 : !-------------------------------------------------------------------
574 : ! evaporation
575 : !-------------------------------------------------------------------
576 28226698 : dzint = dzint + min(dzssl + dhs_evap, c0)
577 28226698 : dzssl = max(dzssl + dhs_evap, c0)
578 :
579 : !-------------------------------------------------------------------
580 : ! surface snow melt
581 : !-------------------------------------------------------------------
582 28226698 : if (-dhs_melts > puny) then
583 15376560 : do k = 1, nbtrcr
584 14607732 : sloss1 = c0
585 14607732 : sloss2 = c0
586 14607732 : if (dzssl > puny) &
587 2299 : sloss1 = kscavz(k)*aerosno(k,1) &
588 14607428 : *min(-dhs_melts,dzssl)/dzssl
589 14607732 : aerosno(k,1) = aerosno(k,1) - sloss1
590 14607732 : if (dzint > puny) &
591 2318 : sloss2 = kscavz(k)*aerosno(k,2) &
592 14607732 : *max(-dhs_melts-dzssl,c0)/dzint
593 14607732 : aerosno(k,2) = aerosno(k,2) - sloss2
594 15376560 : zbgc_snow(k) = zbgc_snow(k) + (sloss1+sloss2)
595 : enddo !
596 :
597 : ! update snow thickness
598 768828 : dzint=dzint+min(dzssl+dhs_melts, c0)
599 768828 : dzssl=max(dzssl+dhs_melts, c0)
600 :
601 768828 : if ( dzssl <= puny ) then ! ssl melts away
602 306120 : aerosno(:,2) = aerosno(:,1) + aerosno(:,2)
603 306120 : aerosno(:,1) = c0
604 15306 : dzssl = max(dzssl, c0)
605 : endif
606 768828 : if (dzint <= puny ) then ! all snow melts away
607 0 : zbgc_snow(:) = zbgc_snow(:) &
608 460 : + max(c0,aerosno(:,1) + aerosno(:,2))
609 943 : aerosno(:,:) = c0
610 23 : dzint = max(dzint, c0)
611 : endif
612 : endif
613 :
614 : !-------------------------------------------------------------------
615 : ! snowfall
616 : !-------------------------------------------------------------------
617 28226698 : if (fsnow > c0) dzssl = dzssl + fsnow/rhos*dt
618 :
619 : !-------------------------------------------------------------------
620 : ! snow-ice formation
621 : !-------------------------------------------------------------------
622 28226698 : if (dhs_snoice > puny) then
623 20403880 : do k = 1, nbtrcr
624 19383686 : sloss1 = c0
625 19383686 : sloss2 = c0
626 19383686 : if (dzint > puny) &
627 : sloss2 = min(dhs_snoice, dzint) &
628 19383686 : *aerosno(k,2)/dzint
629 19383686 : aerosno(k,2) = aerosno(k,2) - sloss2
630 19383686 : if (dzssl > puny) &
631 : sloss1 = max(dhs_snoice-dzint, c0) &
632 19377663 : *aerosno(k,1)/dzssl
633 19383686 : aerosno(k,1) = aerosno(k,1) - sloss1
634 4351 : zbgc_snow(k) = zbgc_snow(k) &
635 20403880 : + (sloss2+sloss1)
636 : enddo
637 1020194 : dzssl = dzssl - max(dhs_snoice-dzint, c0)
638 1020194 : dzint = max(dzint-dhs_snoice, c0)
639 : endif
640 :
641 : !-------------------------------------------------------------------
642 : ! aerosol deposition
643 : !-------------------------------------------------------------------
644 28226698 : if (aicen > c0) then
645 28226698 : hs = vsnon * ar
646 : else
647 0 : hs = c0
648 : endif
649 28226698 : if (hs >= hs_min) then !should this really be hs_min or 0?
650 : ! should use same hs_min value as in radiation
651 550551700 : do k=1,nbtrcr
652 1173421 : aerosno(k,1) = aerosno(k,1) &
653 550551700 : + flux_bio_atm(k)*dt
654 : enddo
655 : else
656 13982260 : do k=1,nbtrcr
657 199291 : zbgc_atm(k) = zbgc_atm(k) &
658 13982260 : + flux_bio_atm(k)*dt
659 : enddo
660 : endif
661 :
662 : !-------------------------------------------------------------------
663 : ! redistribute aerosol within vertical layers
664 : !-------------------------------------------------------------------
665 28226698 : if (aicen > c0) then
666 28226698 : hs = vsnon * ar ! new snow thickness
667 : else
668 0 : hs = c0
669 : endif
670 28226698 : if (dzssl <= puny) then ! nothing in SSL
671 356860 : do k=1,nbtrcr
672 339017 : aerosno(k,2) = aerosno(k,2) + aerosno(k,1)
673 356860 : aerosno(k,1) = c0
674 : enddo
675 : endif
676 28226698 : if (dzint <= puny) then ! nothing in Snow Int
677 14944340 : do k = 1, nbtrcr
678 14197123 : zbgc_snow(k) = zbgc_snow(k) + max(c0,aerosno(k,2))
679 14944340 : aerosno(k,2) = c0
680 : enddo
681 : endif
682 :
683 28226698 : hslyr = hs/real(nslyr,kind=dbl_kind)
684 28226698 : dzssl_new = min(hslyr/c2, hs_ssl)
685 28226698 : dzint_new = hs - dzssl_new
686 :
687 28226698 : if (hs > hs_min) then !should this really be hs_min or 0?
688 550551700 : do k = 1, nbtrcr
689 523024115 : dznew = min(dzssl_new-dzssl, c0)
690 523024115 : sloss1 = c0
691 523024115 : if (dzssl > puny) &
692 522763074 : sloss1 = dznew*aerosno(k,1)/dzssl ! not neccesarily a loss
693 523024115 : dznew = max(dzssl_new-dzssl, c0)
694 523024115 : if (dzint > puny) &
695 518820004 : sloss1 = sloss1 + aerosno(k,2)*dznew/dzint
696 523024115 : aerosno(k,1) = aerosno(k,1) + sloss1
697 550551700 : aerosno(k,2) = aerosno(k,2) - sloss1
698 : enddo
699 : else
700 0 : zbgc_snow(:) = zbgc_snow(:) &
701 13982260 : + max(c0,aerosno(:,1) + aerosno(:,2))
702 28663633 : aerosno(:,:) = c0
703 : endif
704 :
705 : !-------------------------------------------------------------------
706 : ! check conservation
707 : !-------------------------------------------------------------------
708 564533960 : do k = 1, nbtrcr
709 1372712 : aerotot(k) = aerosno(k,2) + aerosno(k,1) &
710 536307262 : + zbgc_snow(k) + zbgc_atm(k)
711 1372712 : aero_cons(k) = aerotot(k)-aerotot0(k) &
712 1372712 : - ( flux_bio_atm(k) &
713 536307262 : - (flux_bio(k)-flux_bio_o(k))) * dt
714 564533960 : if (aero_cons(k) > puny .or. zbgc_snow(k) + zbgc_atm(k) < c0) then
715 0 : write(warnstr,*) subname, 'Conservation failure: aerosols in snow'
716 0 : call icepack_warnings_add(warnstr)
717 0 : write(warnstr,*) subname, 'test aerosol 1'
718 0 : call icepack_warnings_add(warnstr)
719 0 : write(warnstr,*) subname, 'aerosol tracer: ',k
720 0 : call icepack_warnings_add(warnstr)
721 0 : write(warnstr,*) subname, 'aero_cons(k),puny:', aero_cons(k),puny
722 0 : call icepack_warnings_add(warnstr)
723 0 : write(warnstr,*) subname, 'aerotot,aerotot0 ',aerotot(k),aerotot0(k)
724 0 : call icepack_warnings_add(warnstr)
725 0 : write(warnstr,*) subname, ' aerosno(k,2),aerosno(k,1) ', aerosno(k,2),aerosno(k,1)
726 0 : call icepack_warnings_add(warnstr)
727 0 : write(warnstr,*) subname, 'flux_bio_atm(k)*aicen*dt', &
728 0 : flux_bio_atm(k)*aicen*dt
729 0 : call icepack_warnings_add(warnstr)
730 0 : write(warnstr,*) subname, 'zbgc_snow(k)', &
731 0 : zbgc_snow(k)
732 0 : call icepack_warnings_add(warnstr)
733 0 : write(warnstr,*) subname, 'zbgc_atm(k)', &
734 0 : zbgc_atm(k)
735 0 : call icepack_warnings_add(warnstr)
736 : endif
737 : enddo
738 :
739 : !-------------------------------------------------------------------
740 : ! reload tracers
741 : !-------------------------------------------------------------------
742 28226698 : if (vsnon > puny) then
743 564533960 : do k = 1,nbtrcr
744 536307262 : trcrn(bio_index(k)+nblyr+1)=aerosno(k,1)/dzssl_new
745 564533960 : trcrn(bio_index(k)+nblyr+2)=aerosno(k,2)/dzint_new
746 : enddo
747 : else
748 0 : do k = 1,nbtrcr
749 0 : zbgc_snow(k) = (zbgc_snow(k) + aerosno(k,1) + aerosno(k,2))
750 0 : trcrn(bio_index(k)+nblyr+1)= c0
751 0 : trcrn(bio_index(k)+nblyr+2)= c0
752 : enddo
753 : endif
754 : !-------------------------------------------------------------------
755 : ! check for negative values
756 : !-------------------------------------------------------------------
757 1100913470 : if (minval(aerosno(:,1)) < -puny .or. &
758 72248 : minval(aerosno(:,2)) < -puny) then
759 :
760 0 : write(warnstr,*) subname, 'Snow aerosol negative in update_snow_bgc'
761 0 : call icepack_warnings_add(warnstr)
762 0 : write(warnstr,*) subname, 'aicen= ' ,aicen
763 0 : call icepack_warnings_add(warnstr)
764 0 : write(warnstr,*) subname, 'vicen= ' ,vicen
765 0 : call icepack_warnings_add(warnstr)
766 0 : write(warnstr,*) subname, 'vsnon= ' ,vsnon
767 0 : call icepack_warnings_add(warnstr)
768 0 : write(warnstr,*) subname, 'viceold= ' ,vice_old
769 0 : call icepack_warnings_add(warnstr)
770 0 : write(warnstr,*) subname, 'vsnoold= ' ,vsno_old
771 0 : call icepack_warnings_add(warnstr)
772 0 : write(warnstr,*) subname, 'melts= ' ,melts
773 0 : call icepack_warnings_add(warnstr)
774 0 : write(warnstr,*) subname, 'meltt= ' ,meltt
775 0 : call icepack_warnings_add(warnstr)
776 0 : write(warnstr,*) subname, 'meltb= ' ,meltb
777 0 : call icepack_warnings_add(warnstr)
778 0 : write(warnstr,*) subname, 'congel= ' ,congel
779 0 : call icepack_warnings_add(warnstr)
780 0 : write(warnstr,*) subname, 'snoice= ' ,snoice
781 0 : call icepack_warnings_add(warnstr)
782 0 : write(warnstr,*) subname, 'aero evap snow= ' ,dhs_evap
783 0 : call icepack_warnings_add(warnstr)
784 0 : write(warnstr,*) subname, 'fsnow= ' ,fsnow
785 0 : call icepack_warnings_add(warnstr)
786 0 : do k = 1, nbtrcr
787 0 : write(warnstr,*) subname, 'NBTRCR value k = ', k
788 0 : call icepack_warnings_add(warnstr)
789 0 : write(warnstr,*) subname, 'aero snowssl (k)= ' ,aerosno0(k,1)
790 0 : call icepack_warnings_add(warnstr)
791 0 : write(warnstr,*) subname, 'aero new snowssl (k)= ',aerosno(k,1)
792 0 : call icepack_warnings_add(warnstr)
793 0 : write(warnstr,*) subname, 'aero snowint (k)= ' ,aerosno0(k,2)
794 0 : call icepack_warnings_add(warnstr)
795 0 : write(warnstr,*) subname, 'aero new snowint(k)= ',aerosno(k,2)
796 0 : call icepack_warnings_add(warnstr)
797 0 : write(warnstr,*) subname, 'flux_bio_atm(k)= ' , flux_bio_atm(k)
798 0 : call icepack_warnings_add(warnstr)
799 0 : write(warnstr,*) subname, 'zbgc_snow(k)= ' ,zbgc_snow(k)
800 0 : call icepack_warnings_add(warnstr)
801 0 : write(warnstr,*) subname, 'zbgc_atm(k)= ' ,zbgc_atm(k)
802 0 : call icepack_warnings_add(warnstr)
803 :
804 28226698 : do n = 1,2
805 0 : trcrn(bio_index(k)+nblyr+n)=max(trcrn(bio_index(k)+nblyr+n), c0)
806 : enddo
807 : enddo
808 : endif
809 : endif
810 :
811 29415679 : end subroutine update_snow_bgc
812 :
813 : !=======================================================================
814 :
815 : end module icepack_aerosol
816 :
817 : !=======================================================================
|
