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, p5, puny, rhoi, rhos, hs_min
12 : use icepack_parameters, only: hi_ssl, hs_ssl, hs_ssl_min
13 : use icepack_tracers, only: max_aero, nilyr, nslyr, nblyr, ntrcr, nbtrcr, n_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 : ! Deposition and vertical cycling of aerosol in ice or snow
31 : ! Called from icepack_step_therm1 when tr_aero=T (not used for zbgc tracers)
32 :
33 43852465 : subroutine update_aerosol(dt, &
34 : meltt, melts, & ! LCOV_EXCL_LINE
35 : meltb, congel, & ! LCOV_EXCL_LINE
36 : snoice, & ! LCOV_EXCL_LINE
37 : fsnow, & ! LCOV_EXCL_LINE
38 43852465 : aerosno, aeroice, & ! LCOV_EXCL_LINE
39 : aice_old, & ! LCOV_EXCL_LINE
40 : vice_old, vsno_old, & ! LCOV_EXCL_LINE
41 : vicen, vsnon, aicen, & ! LCOV_EXCL_LINE
42 87704930 : faero_atm, faero_ocn)
43 :
44 : real (kind=dbl_kind), intent(in) :: &
45 : dt, & ! time step ! LCOV_EXCL_LINE
46 : meltt, & ! thermodynamic melt/growth rates ! LCOV_EXCL_LINE
47 : melts, & ! LCOV_EXCL_LINE
48 : meltb, & ! LCOV_EXCL_LINE
49 : congel, & ! LCOV_EXCL_LINE
50 : snoice, & ! LCOV_EXCL_LINE
51 : fsnow, & ! LCOV_EXCL_LINE
52 : vicen, & ! ice volume (m) ! LCOV_EXCL_LINE
53 : vsnon, & ! snow volume (m) ! LCOV_EXCL_LINE
54 : aicen, & ! ice area fraction ! LCOV_EXCL_LINE
55 : aice_old, & ! values prior to thermodynamic changes ! LCOV_EXCL_LINE
56 : vice_old, & ! LCOV_EXCL_LINE
57 : vsno_old
58 :
59 : real (kind=dbl_kind), dimension(:), intent(in) :: &
60 : faero_atm ! aerosol deposition rate (W/m^2 s)
61 :
62 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
63 : faero_ocn ! aerosol flux to ocean (W/m^2 s)
64 :
65 : real (kind=dbl_kind), dimension(:,:), intent(inout) :: &
66 : aerosno, aeroice ! kg/m^2
67 :
68 : ! local variables
69 : integer (kind=int_kind) :: k
70 :
71 : real (kind=dbl_kind) :: &
72 : dzssl, dzssl_new, & ! snow ssl thickness ! LCOV_EXCL_LINE
73 : dzint, & ! snow interior thickness ! LCOV_EXCL_LINE
74 : dzssli, dzssli_new, & ! ice ssl thickness ! LCOV_EXCL_LINE
75 : dzinti, & ! ice interior thickness ! LCOV_EXCL_LINE
76 : dznew, & ! tracks thickness changes ! LCOV_EXCL_LINE
77 : hs, hi, & ! snow/ice thickness (m) ! LCOV_EXCL_LINE
78 : dhs_evap, dhi_evap, & ! snow/ice thickness change due to evap ! LCOV_EXCL_LINE
79 : dhs_melts, dhi_meltt, & ! ... due to surface melt ! LCOV_EXCL_LINE
80 : dhs_snoice, dhi_snoice, & ! ... due to snow-ice formation ! LCOV_EXCL_LINE
81 : dhi_congel, dhi_meltb, & ! ... due to bottom growth, melt ! LCOV_EXCL_LINE
82 : hslyr, hilyr, & ! snow, ice layer thickness (m) ! LCOV_EXCL_LINE
83 : hslyr_old, hilyr_old, & ! old snow, ice layer thickness (m) ! LCOV_EXCL_LINE
84 : hs_old, hi_old, & ! old snow, ice thickness (m) ! LCOV_EXCL_LINE
85 : sloss1, sloss2, & ! aerosol mass loss (kg/m^2) ! LCOV_EXCL_LINE
86 : ar ! 1/aicen(i,j)
87 :
88 : real (kind=dbl_kind), dimension(max_aero) :: &
89 : kscav, kscavsi ! scavenging by melt water
90 :
91 : real (kind=dbl_kind), dimension(n_aero) :: &
92 87704930 : aerotot, aerotot0, & ! for conservation check ! LCOV_EXCL_LINE
93 87704930 : focn_old ! for conservation check
94 :
95 : ! echmod: this assumes max_aero=6
96 : data kscav / .03_dbl_kind, .20_dbl_kind, .02_dbl_kind, &
97 : .02_dbl_kind, .01_dbl_kind, .01_dbl_kind /
98 : data kscavsi / .03_dbl_kind, .20_dbl_kind, .02_dbl_kind, &
99 : .02_dbl_kind, .01_dbl_kind, .01_dbl_kind /
100 :
101 : character(len=*),parameter :: subname='(update_aerosol)'
102 :
103 : !-------------------------------------------------------------------
104 : ! initialize
105 : !-------------------------------------------------------------------
106 87704930 : focn_old(:) = faero_ocn(:)
107 :
108 43852465 : hs_old = vsno_old/aice_old
109 43852465 : hi_old = vice_old/aice_old
110 43852465 : hslyr_old = hs_old/real(nslyr,kind=dbl_kind)
111 43852465 : hilyr_old = hi_old/real(nilyr,kind=dbl_kind)
112 :
113 43852465 : dzssl = min(hslyr_old/c2, hs_ssl)
114 43852465 : dzssli = min(hilyr_old/c2, hi_ssl)
115 43852465 : dzint = hs_old - dzssl
116 43852465 : dzinti = hi_old - dzssli
117 :
118 43852465 : if (aicen > c0) then
119 43852162 : ar = c1/aicen
120 : else ! ice disappeared during time step
121 303 : ar = c1/aice_old
122 : endif
123 :
124 43852465 : hs = vsnon*ar
125 43852465 : hi = vicen*ar
126 : ! fluxes were divided by aice for thermo, not yet multiplied by aice
127 43852465 : dhs_melts = -melts
128 43852465 : dhi_snoice = snoice
129 43852465 : dhs_snoice = dhi_snoice*rhoi/rhos
130 43852465 : dhi_meltt = -meltt
131 43852465 : dhi_meltb = -meltb
132 43852465 : dhi_congel = congel
133 :
134 : dhs_evap = hs - (hs_old + dhs_melts - dhs_snoice &
135 43852465 : + fsnow/rhos*dt)
136 : dhi_evap = hi - (hi_old + dhi_meltt + dhi_meltb &
137 43852465 : + dhi_congel + dhi_snoice)
138 :
139 87704930 : do k = 1, n_aero
140 : aerotot0(k) = aerosno(k,2) + aerosno(k,1) &
141 87704930 : + aeroice(k,2) + aeroice(k,1)
142 : enddo
143 :
144 : !-------------------------------------------------------------------
145 : ! evaporation
146 : !-------------------------------------------------------------------
147 43852465 : dzint = dzint + min(dzssl + dhs_evap, c0)
148 43852465 : dzinti = dzinti + min(dzssli + dhi_evap, c0)
149 43852465 : dzssl = max(dzssl + dhs_evap, c0)
150 43852465 : dzssli = max(dzssli + dhi_evap, c0)
151 :
152 : !-------------------------------------------------------------------
153 : ! basal ice growth
154 : !-------------------------------------------------------------------
155 43852465 : dzinti = dzinti + dhi_congel
156 :
157 : !-------------------------------------------------------------------
158 : ! surface snow melt
159 : !-------------------------------------------------------------------
160 43852465 : if (-dhs_melts > puny) then
161 8117280 : do k = 1, n_aero
162 4058640 : sloss1 = c0
163 4058640 : sloss2 = c0
164 4058640 : if (dzssl > puny) &
165 : sloss1 = kscav(k)*aerosno(k,1) & ! LCOV_EXCL_LINE
166 4046919 : *min(-dhs_melts,dzssl)/dzssl
167 4058640 : aerosno(k,1) = aerosno(k,1) - sloss1
168 4058640 : if (dzint > puny) &
169 : sloss2 = kscav(k)*aerosno(k,2) & ! LCOV_EXCL_LINE
170 4058062 : *max(-dhs_melts-dzssl,c0)/dzint
171 4058640 : aerosno(k,2) = aerosno(k,2) - sloss2
172 8117280 : faero_ocn(k) = faero_ocn(k) + (sloss1+sloss2)/dt
173 : enddo ! n_aero
174 :
175 : ! update snow thickness
176 4058640 : dzint=dzint+min(dzssl+dhs_melts, c0)
177 4058640 : dzssl=max(dzssl+dhs_melts, c0)
178 :
179 4058640 : if ( dzssl <= puny ) then ! ssl melts away
180 1041102 : aerosno(:,2) = aerosno(:,1) + aerosno(:,2)
181 1041102 : aerosno(:,1) = c0
182 520551 : dzssl = max(dzssl, c0)
183 : endif
184 4058640 : if (dzint <= puny ) then ! all snow melts away
185 : aeroice(:,1) = aeroice(:,1) &
186 967382 : + aerosno(:,1) + aerosno(:,2)
187 2418455 : aerosno(:,:) = c0
188 483691 : dzint = max(dzint, c0)
189 : endif
190 : endif
191 :
192 : !-------------------------------------------------------------------
193 : ! surface ice melt
194 : !-------------------------------------------------------------------
195 43852465 : if (-dhi_meltt > puny) then
196 3649702 : do k = 1, n_aero
197 1824851 : sloss1 = c0
198 1824851 : sloss2 = c0
199 1824851 : if (dzssli > puny) &
200 : sloss1 = kscav(k)*aeroice(k,1) & ! LCOV_EXCL_LINE
201 1824851 : *min(-dhi_meltt,dzssli)/dzssli
202 1824851 : aeroice(k,1) = aeroice(k,1) - sloss1
203 1824851 : if (dzinti > puny) &
204 : sloss2 = kscav(k)*aeroice(k,2) & ! LCOV_EXCL_LINE
205 1824851 : *max(-dhi_meltt-dzssli,c0)/dzinti
206 1824851 : aeroice(k,2) = aeroice(k,2) - sloss2
207 3649702 : faero_ocn(k) = faero_ocn(k) + (sloss1+sloss2)/dt
208 : enddo
209 :
210 1824851 : dzinti = dzinti + min(dzssli+dhi_meltt, c0)
211 1824851 : dzssli = max(dzssli+dhi_meltt, c0)
212 1824851 : if (dzssli <= puny) then ! ssl ice melts away
213 15568 : do k = 1, n_aero
214 7784 : aeroice(k,2) = aeroice(k,1) + aeroice(k,2)
215 15568 : aeroice(k,1) = c0
216 : enddo
217 7784 : dzssli = max(dzssli, c0)
218 : endif
219 1824851 : if (dzinti <= puny) then ! all ice melts away
220 0 : do k = 1, n_aero
221 : faero_ocn(k) = faero_ocn(k) &
222 0 : + (aeroice(k,1)+aeroice(k,2))/dt
223 0 : aeroice(k,:)=c0
224 : enddo
225 0 : dzinti = max(dzinti, c0)
226 : endif
227 : endif
228 :
229 : !-------------------------------------------------------------------
230 : ! basal ice melt. Assume all aero lost in basal melt
231 : !-------------------------------------------------------------------
232 43852465 : if (-dhi_meltb > puny) then
233 31717478 : do k=1,n_aero
234 15858739 : sloss1=c0
235 15858739 : sloss2=c0
236 15858739 : if (dzssli > puny) &
237 : sloss1 = max(-dhi_meltb-dzinti, c0) & ! LCOV_EXCL_LINE
238 15850955 : *aeroice(k,1)/dzssli
239 15858739 : aeroice(k,1) = aeroice(k,1) - sloss1
240 15858739 : if (dzinti > puny) &
241 : sloss2 = min(-dhi_meltb, dzinti) & ! LCOV_EXCL_LINE
242 15858739 : *aeroice(k,2)/dzinti
243 15858739 : aeroice(k,2) = aeroice(k,2) - sloss2
244 31717478 : faero_ocn(k) = faero_ocn(k) + (sloss1+sloss2)/dt
245 : enddo
246 :
247 15858739 : dzssli = dzssli + min(dzinti+dhi_meltb, c0)
248 15858739 : dzinti = max(dzinti+dhi_meltb, c0)
249 : endif
250 :
251 : !-------------------------------------------------------------------
252 : ! snowfall
253 : !-------------------------------------------------------------------
254 43852465 : if (fsnow > c0) dzssl = dzssl + fsnow/rhos*dt
255 :
256 : !-------------------------------------------------------------------
257 : ! snow-ice formation
258 : !-------------------------------------------------------------------
259 43852465 : if (dhs_snoice > puny) then
260 3062722 : do k = 1, n_aero
261 1531361 : sloss1 = c0
262 1531361 : sloss2 = c0
263 1531361 : if (dzint > puny) &
264 : sloss2 = min(dhs_snoice, dzint) & ! LCOV_EXCL_LINE
265 1531360 : *aerosno(k,2)/dzint
266 1531361 : aerosno(k,2) = aerosno(k,2) - sloss2
267 1531361 : if (dzssl > puny) &
268 : sloss1 = max(dhs_snoice-dzint, c0) & ! LCOV_EXCL_LINE
269 1531352 : *aerosno(k,1)/dzssl
270 1531361 : aerosno(k,1) = aerosno(k,1) - sloss1
271 : aeroice(k,1) = aeroice(k,1) &
272 1531361 : + (c1-kscavsi(k))*(sloss2+sloss1)
273 : faero_ocn(k) = faero_ocn(k) &
274 3062722 : + kscavsi(k)*(sloss2+sloss1)/dt
275 : enddo
276 1531361 : dzssl = dzssl - max(dhs_snoice-dzint, c0)
277 1531361 : dzint = max(dzint-dhs_snoice, c0)
278 1531361 : dzssli = dzssli + dhi_snoice
279 : endif
280 :
281 : !-------------------------------------------------------------------
282 : ! aerosol deposition
283 : !-------------------------------------------------------------------
284 43852465 : if (aicen > c0) then
285 43852162 : hs = vsnon * ar
286 : else
287 303 : hs = c0
288 : endif
289 43852465 : if (hs > hs_min) then ! should this really be hs_min or 0?
290 : ! should use same hs_min value as in radiation
291 82251748 : do k=1,n_aero
292 : aerosno(k,1) = aerosno(k,1) &
293 82251748 : + faero_atm(k)*dt*aicen
294 : enddo
295 : else
296 5453182 : do k=1,n_aero
297 : aeroice(k,1) = aeroice(k,1) &
298 5453182 : + faero_atm(k)*dt*aicen
299 : enddo
300 : endif
301 :
302 : !-------------------------------------------------------------------
303 : ! redistribute aerosol within vertical layers
304 : !-------------------------------------------------------------------
305 43852465 : if (aicen > c0) then
306 43852162 : hs = vsnon * ar ! new snow thickness
307 43852162 : hi = vicen * ar ! new ice thickness
308 : else
309 303 : hs = c0
310 303 : hi = c0
311 : endif
312 43852465 : if (dzssl <= puny) then ! nothing in SSL
313 4616826 : do k=1,n_aero
314 2308413 : aerosno(k,2) = aerosno(k,2) + aerosno(k,1)
315 4616826 : aerosno(k,1) = c0
316 : enddo
317 : endif
318 43852465 : if (dzint <= puny) then ! nothing in Snow Int
319 4794024 : do k = 1, n_aero
320 2397012 : aeroice(k,1) = aeroice(k,1) + aerosno(k,2)
321 4794024 : aerosno(k,2) = c0
322 : enddo
323 : endif
324 43852465 : if (dzssli <= puny) then ! nothing in Ice SSL
325 15946 : do k = 1, n_aero
326 7973 : aeroice(k,2) = aeroice(k,2) + aeroice(k,1)
327 15946 : aeroice(k,1) = c0
328 : enddo
329 : endif
330 :
331 43852465 : if (dzinti <= puny) then ! nothing in Ice INT
332 1770 : do k = 1, n_aero
333 : faero_ocn(k) = faero_ocn(k) &
334 885 : + (aeroice(k,1)+aeroice(k,2))/dt
335 3540 : aeroice(k,:)=c0
336 : enddo
337 : endif
338 :
339 43852465 : hslyr = hs/real(nslyr,kind=dbl_kind)
340 43852465 : hilyr = hi/real(nilyr,kind=dbl_kind)
341 43852465 : dzssl_new = min(hslyr/c2, hs_ssl)
342 43852465 : dzssli_new = min(hilyr/c2, hi_ssl)
343 :
344 43852465 : if (hs > hs_min) then
345 82251748 : do k = 1, n_aero
346 41125874 : dznew = min(dzssl_new-dzssl, c0)
347 41125874 : sloss1 = c0
348 41125874 : if (dzssl > puny) &
349 41095046 : sloss1 = dznew*aerosno(k,1)/dzssl ! not neccesarily a loss
350 41125874 : dznew = max(dzssl_new-dzssl, c0)
351 41125874 : if (dzint > puny) &
352 41110178 : sloss1 = sloss1 + aerosno(k,2)*dznew/dzint
353 41125874 : aerosno(k,1) = aerosno(k,1) + sloss1
354 82251748 : aerosno(k,2) = aerosno(k,2) - sloss1
355 : enddo
356 : else
357 : aeroice(:,1) = aeroice(:,1) &
358 5453182 : + aerosno(:,1) + aerosno(:,2)
359 13632955 : aerosno(:,:) = c0
360 : endif
361 :
362 43852465 : if (vicen > puny) then ! may want a limit on hi instead?
363 87655774 : do k = 1, n_aero
364 43827887 : sloss2 = c0
365 43827887 : dznew = min(dzssli_new-dzssli, c0)
366 43827887 : if (dzssli > puny) &
367 43820446 : sloss2 = dznew*aeroice(k,1)/dzssli
368 43827887 : dznew = max(dzssli_new-dzssli, c0)
369 43827887 : if (dzinti > puny) &
370 43827334 : sloss2 = sloss2 + aeroice(k,2)*dznew/dzinti
371 43827887 : aeroice(k,1) = aeroice(k,1) + sloss2
372 87655774 : aeroice(k,2) = aeroice(k,2) - sloss2
373 : enddo
374 : else
375 49156 : faero_ocn(:) = faero_ocn(:) + (aeroice(:,1)+aeroice(:,2))/dt
376 122890 : aeroice(:,:) = c0
377 : endif
378 :
379 : !-------------------------------------------------------------------
380 : ! check conservation
381 : !-------------------------------------------------------------------
382 87704930 : do k = 1, n_aero
383 : aerotot(k) = aerosno(k,2) + aerosno(k,1) &
384 43852465 : + aeroice(k,2) + aeroice(k,1)
385 43852465 : if ((aerotot(k)-aerotot0(k)) &
386 : - ( faero_atm(k)*aicen & ! LCOV_EXCL_LINE
387 43852465 : - (faero_ocn(k)-focn_old(k)) )*dt > puny) then
388 :
389 0 : write(warnstr,*) subname, 'aerosol tracer: ',k
390 0 : call icepack_warnings_add(warnstr)
391 0 : write(warnstr,*) subname, 'aerotot-aerotot0 ',aerotot(k)-aerotot0(k)
392 0 : call icepack_warnings_add(warnstr)
393 0 : write(warnstr,*) subname, 'faero_atm-faero_ocn ', &
394 0 : (faero_atm(k)*aicen-(faero_ocn(k)-focn_old(k)))*dt
395 0 : call icepack_warnings_add(warnstr)
396 : endif
397 : enddo
398 :
399 : !-------------------------------------------------------------------
400 : ! check for negative values
401 : !-------------------------------------------------------------------
402 :
403 : !echmod: note that this does not test or fix all aero tracers
404 : if (aeroice(1,1) < -puny .or. &
405 : aeroice(1,2) < -puny .or. & ! LCOV_EXCL_LINE
406 43852465 : aerosno(1,1) < -puny .or. & ! LCOV_EXCL_LINE
407 : aerosno(1,2) < -puny) then
408 :
409 0 : write(warnstr,*) subname, 'aerosol negative in aerosol code'
410 0 : call icepack_warnings_add(warnstr)
411 :
412 0 : aeroice(1,1) = max(aeroice(1,1), c0)
413 0 : aeroice(1,2) = max(aeroice(1,2), c0)
414 0 : aerosno(1,1) = max(aerosno(1,1), c0)
415 0 : aerosno(1,2) = max(aerosno(1,2), c0)
416 :
417 : endif
418 :
419 43852465 : end subroutine update_aerosol
420 :
421 : !=======================================================================
422 :
423 : ! Aerosol in snow for vertical biogeochemistry with mushy thermodynamics
424 : ! Called from icepack_algae.F90 when z_tracers=T (replaces update_aerosol)
425 :
426 33835340 : subroutine update_snow_bgc(dt, &
427 : meltt, melts, & ! LCOV_EXCL_LINE
428 : meltb, congel, & ! LCOV_EXCL_LINE
429 : snoice, fsnow, & ! LCOV_EXCL_LINE
430 33835340 : trcrn, bio_index, & ! LCOV_EXCL_LINE
431 33835340 : aice_old, zbgc_snow, & ! LCOV_EXCL_LINE
432 : vice_old, vsno_old, & ! LCOV_EXCL_LINE
433 : vicen, vsnon, & ! LCOV_EXCL_LINE
434 33835340 : aicen, flux_bio_atm,& ! LCOV_EXCL_LINE
435 33835340 : zbgc_atm, flux_bio, & ! LCOV_EXCL_LINE
436 33835340 : bio_index_o)
437 :
438 : integer (kind=int_kind), dimension (nbtrcr), intent(in) :: &
439 : bio_index, & ! LCOV_EXCL_LINE
440 : bio_index_o ! provides index of scavenging (kscavz) data array
441 :
442 : real (kind=dbl_kind), intent(in) :: &
443 : dt ! time step
444 :
445 : real (kind=dbl_kind), intent(in) :: &
446 : meltt, & ! thermodynamic melt/growth rates ! LCOV_EXCL_LINE
447 : melts, & ! LCOV_EXCL_LINE
448 : meltb, & ! LCOV_EXCL_LINE
449 : congel, & ! LCOV_EXCL_LINE
450 : snoice, & ! LCOV_EXCL_LINE
451 : fsnow, & ! LCOV_EXCL_LINE
452 : vicen, & ! ice volume (m) ! LCOV_EXCL_LINE
453 : vsnon, & ! snow volume (m) ! LCOV_EXCL_LINE
454 : aicen, & ! ice area fraction ! LCOV_EXCL_LINE
455 : aice_old, & ! values prior to thermodynamic changes ! LCOV_EXCL_LINE
456 : vice_old, & ! LCOV_EXCL_LINE
457 : vsno_old
458 :
459 : real (kind=dbl_kind), dimension(nbtrcr), intent(out) :: &
460 : zbgc_snow, & ! aerosol contribution from snow to ice ! LCOV_EXCL_LINE
461 : zbgc_atm ! and atm to ice concentration * volume (kg or mmol/m^3*m)
462 :
463 : real (kind=dbl_kind),dimension(nbtrcr), intent(inout) :: &
464 : flux_bio ! total ocean tracer flux (mmol/m^2/s)
465 :
466 : real (kind=dbl_kind), dimension(nbtrcr), intent(in) :: &
467 : flux_bio_atm ! aerosol deposition rate (kg or mmol/m^2 s)
468 :
469 : real (kind=dbl_kind), dimension(ntrcr), intent(inout) :: &
470 : trcrn ! ice/snow tracer array
471 :
472 : ! local variables
473 :
474 : integer (kind=int_kind) :: k, n
475 :
476 : real (kind=dbl_kind) :: &
477 : dzssl, dzssl_new, & ! snow ssl thickness ! LCOV_EXCL_LINE
478 : dzint, dzint_new, & ! snow interior thickness ! LCOV_EXCL_LINE
479 : dz, & ! ! LCOV_EXCL_LINE
480 : hi, & ! ice thickness (m) ! LCOV_EXCL_LINE
481 : hilyr, & ! ice layer thickness (m) ! LCOV_EXCL_LINE
482 : hs, & ! snow thickness (m) ! LCOV_EXCL_LINE
483 : dhs_evap, & ! snow thickness change due to evap ! LCOV_EXCL_LINE
484 : dhs_melts, & ! ... due to surface melt ! LCOV_EXCL_LINE
485 : dhs_snoice, & ! ... due to snow-ice formation ! LCOV_EXCL_LINE
486 : hslyr, & ! snow layer thickness (m) ! LCOV_EXCL_LINE
487 : hslyr_old, & ! old snow layer thickness (m) ! LCOV_EXCL_LINE
488 : hs_old, & ! old snow thickness (m) ! LCOV_EXCL_LINE
489 : dznew, & ! change in the snow sl (m) ! LCOV_EXCL_LINE
490 : sloss1, sloss2, & ! aerosol mass loss (kg/m^2) ! LCOV_EXCL_LINE
491 : ar ! 1/aicen(i,j)
492 :
493 : real (kind=dbl_kind), dimension(nbtrcr) :: &
494 67670680 : aerotot, aerotot0, & ! for conservation check (mmol/m^3) ! LCOV_EXCL_LINE
495 67670680 : aero_cons , & ! for conservation check (mmol/m^2) ! LCOV_EXCL_LINE
496 67670680 : flux_bio_o ! initial ocean tracer flux (mmol/m^2/s)
497 :
498 : real (kind=dbl_kind), dimension(nbtrcr,2) :: &
499 67670680 : aerosno, & ! kg/m^2 ! LCOV_EXCL_LINE
500 67670680 : aerosno0 ! for diagnostic prints
501 :
502 : character(len=*),parameter :: subname='(update_snow_bgc)'
503 :
504 : !-------------------------------------------------------------------
505 : ! initialize
506 : !-------------------------------------------------------------------
507 1225590436 : aerosno (:,:) = c0
508 1225590436 : aerosno0(:,:) = c0
509 595877548 : aero_cons(:) = c0
510 595877548 : zbgc_snow(:) = c0
511 595877548 : zbgc_atm(:) = c0
512 :
513 33835340 : hs_old = vsno_old/aice_old
514 33835340 : if (aice_old .gt. puny) then
515 33835340 : hs_old = vsno_old/aice_old
516 : else
517 0 : hs_old = c0
518 : end if
519 33835340 : hslyr_old = hs_old/real(nslyr,kind=dbl_kind)
520 :
521 33835340 : dzssl = hslyr_old/c2
522 33835340 : dzint = hs_old - dzssl
523 :
524 33835340 : if (aicen > c0) then
525 33835340 : ar = c1/aicen
526 33835340 : hs = vsnon*ar
527 33835340 : hi = vicen*ar
528 : else ! ice disappeared during time step
529 0 : ar = c1
530 0 : hs = c0
531 0 : hi = c0
532 0 : if (aice_old > c0) hs = vsnon/aice_old
533 : endif
534 33835340 : hilyr = hi/real(nblyr,kind=dbl_kind)
535 33835340 : hslyr = hs/real(nslyr,kind=dbl_kind)
536 33835340 : dzssl_new = hslyr/c2
537 33835340 : dhs_melts = -melts
538 33835340 : dhs_snoice = snoice*rhoi/rhos
539 : dhs_evap = hs - (hs_old + dhs_melts - dhs_snoice &
540 33835340 : + fsnow/rhos*dt)
541 :
542 : ! trcrn() has units kg/m^3
543 :
544 33835340 : if (dzssl_new .lt. hs_ssl_min) then ! Put atm BC/dust flux directly into the sea ice
545 274209482 : do k=1,nbtrcr
546 259882173 : flux_bio_o(k) = flux_bio(k)
547 259882173 : if (hilyr .lt. hs_ssl_min) then
548 : flux_bio(k) = flux_bio(k) + &
549 : (trcrn(bio_index(k)+ nblyr+1)*dzssl+ & ! LCOV_EXCL_LINE
550 273758 : trcrn(bio_index(k)+ nblyr+2)*dzint)/dt
551 273758 : flux_bio(k) = flux_bio(k) + flux_bio_atm(k)
552 : else
553 : zbgc_snow(k) = zbgc_snow(k) + &
554 : (trcrn(bio_index(k)+ nblyr+1)*dzssl+ & ! LCOV_EXCL_LINE
555 259608415 : trcrn(bio_index(k)+ nblyr+2)*dzint)
556 : zbgc_atm(k) = zbgc_atm(k) &
557 259608415 : + flux_bio_atm(k)*dt
558 : end if
559 259882173 : trcrn(bio_index(k) + nblyr+1) = c0
560 274209482 : trcrn(bio_index(k) + nblyr+2) = c0
561 : enddo
562 :
563 : else
564 :
565 321668066 : do k=1,nbtrcr
566 302160035 : flux_bio_o(k) = flux_bio(k)
567 302160035 : aerosno (k,1) = trcrn(bio_index(k)+ nblyr+1) * dzssl
568 302160035 : aerosno (k,2) = trcrn(bio_index(k)+ nblyr+2) * dzint
569 906480105 : aerosno0(k,:) = aerosno(k,:)
570 321668066 : aerotot0(k) = aerosno(k,2) + aerosno(k,1)
571 : enddo
572 :
573 : !-------------------------------------------------------------------
574 : ! evaporation
575 : !-------------------------------------------------------------------
576 19508031 : dzint = dzint + min(dzssl + dhs_evap, c0)
577 19508031 : dzssl = max(dzssl + dhs_evap, c0)
578 19508031 : if (dzssl <= puny) then
579 1046 : do k = 1,nbtrcr
580 963 : aerosno(k,2) = aerosno(k,2) + aerosno(k,1)
581 1046 : aerosno(k,1) = c0
582 : end do
583 : end if
584 19508031 : if (dzint <= puny) then
585 84 : do k = 1,nbtrcr
586 80 : zbgc_snow(k) = zbgc_snow(k) + (aerosno(k,2) + aerosno(k,1))
587 80 : aerosno(k,2) = c0
588 84 : aerosno(k,1) = c0
589 : end do
590 : end if
591 :
592 : !------------------------------------------------------------------
593 : ! snowfall
594 : !-------------------------------------------------------------------
595 19508031 : if (fsnow > c0) then
596 15644626 : sloss1 = c0
597 15644626 : dz = min(fsnow/rhos*dt,dzssl)
598 268575511 : do k = 1, nbtrcr
599 252930885 : if (dzssl > puny) &
600 252929922 : sloss1 = aerosno(k,1)*dz/dzssl
601 252930885 : aerosno(k,1) = max(c0,aerosno(k,1) - sloss1)
602 268575511 : aerosno(k,2) = aerosno(k,2) + sloss1
603 : end do
604 15644626 : dzssl = dzssl - dz + fsnow/rhos*dt
605 15644626 : dzint = dzint + dz
606 : end if
607 19508031 : if (dzssl <= puny) then
608 0 : do k = 1,nbtrcr
609 0 : aerosno(k,2) = aerosno(k,2) + aerosno(k,1)
610 0 : aerosno(k,1) = c0
611 : end do
612 : end if
613 19508031 : if (dzint <= puny) then
614 84 : do k = 1,nbtrcr
615 80 : zbgc_snow(k) = zbgc_snow(k) + (aerosno(k,2) + aerosno(k,1))
616 80 : aerosno(k,2) = c0
617 84 : aerosno(k,1) = c0
618 : end do
619 : end if
620 :
621 : !-------------------------------------------------------------------
622 : ! surface snow melt
623 : !-------------------------------------------------------------------
624 19508031 : if (-dhs_melts > puny) then
625 9534906 : do k = 1, nbtrcr
626 8493865 : sloss1 = c0
627 8493865 : sloss2 = c0
628 8493865 : if (dzssl > puny) &
629 : sloss1 = kscavz(bio_index_o(k))*aerosno(k,1) & ! LCOV_EXCL_LINE
630 8493865 : *min(-dhs_melts,dzssl)/dzssl
631 8493865 : aerosno(k,1) = max(c0,aerosno(k,1) - sloss1)
632 8493865 : if (dzint > puny) &
633 : sloss2 = kscavz(bio_index_o(k))*aerosno(k,2) & ! LCOV_EXCL_LINE
634 8493865 : *max(-dhs_melts-dzssl,c0)/dzint
635 8493865 : aerosno(k,2) = max(c0,aerosno(k,2) - sloss2)
636 9534906 : zbgc_snow(k) = zbgc_snow(k) + (sloss1+sloss2) ! all not scavenged ends in ice
637 : enddo
638 :
639 : ! update snow thickness
640 1041041 : dzint=dzint+min(dzssl+dhs_melts, c0)
641 1041041 : dzssl=max(dzssl+dhs_melts, c0)
642 :
643 1041041 : if ( dzssl .le. puny ) then ! ssl melts away
644 908948 : do k = 1,nbtrcr
645 860160 : aerosno(k,2) = aerosno(k,1) + aerosno(k,2)
646 908948 : aerosno(k,1) = c0
647 : end do
648 48788 : dzssl = max(dzssl, c0)
649 : endif
650 1041041 : if (dzint .le. puny ) then ! all snow melts away
651 0 : do k = 1,nbtrcr
652 : zbgc_snow(k) = zbgc_snow(k) &
653 0 : + aerosno(k,1) + aerosno(k,2)
654 0 : aerosno(k,:) = c0
655 : enddo
656 0 : dzint = max(dzint, c0)
657 : endif
658 : endif ! -dhs_melts > puny
659 :
660 : !-------------------------------------------------------------------
661 : ! snow-ice formation
662 : !-------------------------------------------------------------------
663 19508031 : if (dhs_snoice > puny) then
664 5865147 : do k = 1, nbtrcr
665 5514566 : sloss1 = c0
666 5514566 : sloss2 = c0
667 5514566 : if (dzint > puny .and. aerosno(k,2) > c0) &
668 : sloss2 = min(dhs_snoice, dzint) & ! LCOV_EXCL_LINE
669 1033485 : *aerosno(k,2)/dzint
670 5514566 : aerosno(k,2) = max(c0,aerosno(k,2) - sloss2)
671 5514566 : if (dzssl > puny .and. aerosno(k,1) > c0) &
672 : sloss1 = max(dhs_snoice-dzint, c0) & ! LCOV_EXCL_LINE
673 1027232 : *aerosno(k,1)/dzssl
674 :
675 5514566 : aerosno(k,1) = max(c0,aerosno(k,1) - sloss1)
676 : flux_bio(k) = flux_bio(k) &
677 5514566 : + kscavz(bio_index_o(k)) * (sloss2+sloss1)/dt
678 : zbgc_snow(k) = zbgc_snow(k) &
679 5865147 : + (c1-kscavz(bio_index_o(k)))*(sloss2+sloss1)
680 : enddo
681 350581 : dzssl = max(c0,dzssl - max(dhs_snoice-dzint, c0))
682 350581 : dzint = max(dzint-dhs_snoice, c0)
683 : endif ! dhs_snowice > puny
684 :
685 : !-------------------------------------------------------------------
686 : ! aerosol deposition
687 : !-------------------------------------------------------------------
688 : ! Spread out the atm dust flux in the snow interior for small snow surface layers
689 19508031 : if (dzssl .ge. hs_ssl*p5) then
690 :
691 42657420 : do k=1,nbtrcr
692 : aerosno(k,1) = aerosno(k,1) &
693 42657420 : + flux_bio_atm(k)*dt
694 : enddo
695 : else
696 14650587 : dz = (hs_ssl*p5 - dzssl)/(hs_ssl*p5)
697 279010646 : do k=1,nbtrcr
698 : aerosno(k,1) = aerosno(k,1) &
699 264360059 : + flux_bio_atm(k)*dt*(c1-dz)
700 : aerosno(k,2) = aerosno(k,2) &
701 279010646 : + flux_bio_atm(k)*dt*dz
702 : enddo
703 : endif
704 :
705 : !-------------------------------------------------------------------
706 : ! redistribute aerosol within vertical layers
707 : !-------------------------------------------------------------------
708 19508031 : if (aicen > c0) then
709 19508031 : hs = vsnon * ar ! new snow thickness
710 : else
711 0 : hs = c0
712 : endif
713 19508031 : if (dzssl <= puny) then ! nothing in SSL
714 908948 : do k=1,nbtrcr
715 860160 : aerosno(k,2) = aerosno(k,2) + aerosno(k,1)
716 908948 : aerosno(k,1) = c0
717 : enddo
718 : endif
719 19508031 : if (dzint <= puny) then ! nothing in Snow Int
720 62636 : do k = 1, nbtrcr
721 54593 : zbgc_snow(k) = zbgc_snow(k) + max(c0,aerosno(k,2)+aerosno(k,1))
722 54593 : aerosno(k,1) = c0
723 62636 : aerosno(k,2) = c0
724 : enddo
725 : endif
726 :
727 19508031 : hslyr = hs/real(nslyr,kind=dbl_kind)
728 19508031 : dzssl_new = hslyr/c2
729 19508031 : dzint_new = max(c0,hs - dzssl_new)
730 :
731 19508031 : if (hs > hs_min) then
732 321668066 : do k = 1, nbtrcr
733 302160035 : dznew = min(dzssl_new-dzssl, c0)
734 302160035 : sloss1 = c0
735 302160035 : if (dzssl > puny .and. aerosno(k,1) > c0) &
736 58353600 : sloss1 = dznew*aerosno(k,1)/dzssl ! not neccesarily a loss
737 302160035 : dznew = max(dzssl_new-dzssl, c0)
738 302160035 : if (dzint > puny .and. aerosno(k,2) > c0) &
739 58445904 : sloss1 = aerosno(k,2)*dznew/dzint
740 302160035 : aerosno(k,1) = max(c0,aerosno(k,1) + sloss1)
741 321668066 : aerosno(k,2) = max(c0,aerosno(k,2) - sloss1)
742 : enddo
743 : else
744 : zbgc_snow(:) = zbgc_snow(:) &
745 0 : + aerosno(:,1) + aerosno(:,2)
746 0 : aerosno(:,:) = c0
747 : endif
748 :
749 : !-------------------------------------------------------------------
750 : ! check conservation
751 : !-------------------------------------------------------------------
752 321668066 : do k = 1, nbtrcr
753 : aerotot(k) = aerosno(k,2) + aerosno(k,1) &
754 302160035 : + zbgc_snow(k) + zbgc_atm(k)
755 : aero_cons(k) = aerotot(k)-aerotot0(k) &
756 : - ( flux_bio_atm(k) & ! LCOV_EXCL_LINE
757 302160035 : - (flux_bio(k)-flux_bio_o(k))) * dt
758 302160035 : if (aerotot0(k) > aerotot(k) .and. aerotot0(k) > c0) then
759 12848 : aero_cons(k) = aero_cons(k)/aerotot0(k)
760 302147187 : else if (aerotot(k) > c0) then
761 58511245 : aero_cons(k) = aero_cons(k)/aerotot(k)
762 : end if
763 321668066 : if (aero_cons(k) > puny .or. zbgc_snow(k) + zbgc_atm(k) < c0) then
764 0 : write(warnstr,*) subname, 'Conservation failure: aerosols in snow'
765 0 : call icepack_warnings_add(warnstr)
766 0 : write(warnstr,*) subname, 'test aerosol 1'
767 0 : call icepack_warnings_add(warnstr)
768 0 : write(warnstr,*) subname, 'aerosol tracer: ',k
769 0 : call icepack_warnings_add(warnstr)
770 0 : write(warnstr,*) subname, 'aero_cons(k),puny:', aero_cons(k),puny
771 0 : call icepack_warnings_add(warnstr)
772 0 : write(warnstr,*) subname, 'aerotot,aerotot0 ',aerotot(k),aerotot0(k)
773 0 : call icepack_warnings_add(warnstr)
774 0 : write(warnstr,*) subname, ' aerosno(k,2),aerosno(k,1) ', aerosno(k,2),aerosno(k,1)
775 0 : call icepack_warnings_add(warnstr)
776 0 : write(warnstr,*) subname, 'flux_bio_atm(k)*aicen*dt', &
777 0 : flux_bio_atm(k)*aicen*dt
778 0 : call icepack_warnings_add(warnstr)
779 0 : write(warnstr,*) subname, 'zbgc_snow(k)', &
780 0 : zbgc_snow(k)
781 0 : call icepack_warnings_add(warnstr)
782 0 : write(warnstr,*) subname, 'zbgc_atm(k)', &
783 0 : zbgc_atm(k)
784 0 : call icepack_warnings_add(warnstr)
785 : endif
786 : enddo
787 :
788 : !-------------------------------------------------------------------
789 : ! reload tracers
790 : !-------------------------------------------------------------------
791 19508031 : if (dzssl_new > puny .and. dzint_new > puny .and. vsnon > puny) then
792 319079117 : do k = 1,nbtrcr
793 299704103 : trcrn(bio_index(k)+nblyr+1)=aerosno(k,1)/dzssl_new
794 319079117 : trcrn(bio_index(k)+nblyr+2)=aerosno(k,2)/dzint_new
795 : enddo
796 : else
797 2588949 : do k = 1,nbtrcr
798 2455932 : trcrn(bio_index(k)+nblyr+1)= c0
799 2588949 : trcrn(bio_index(k)+nblyr+2)= c0
800 : enddo
801 : endif
802 :
803 : !-------------------------------------------------------------------
804 : ! check for negative values
805 : !-------------------------------------------------------------------
806 623828101 : if (minval(aerosno(:,1)) < -puny .or. &
807 : minval(aerosno(:,2)) < -puny) then
808 :
809 0 : write(warnstr,*) subname, 'Snow aerosol negative in update_snow_bgc'
810 0 : call icepack_warnings_add(warnstr)
811 0 : write(warnstr,*) subname, 'aicen= ' ,aicen
812 0 : call icepack_warnings_add(warnstr)
813 0 : write(warnstr,*) subname, 'vicen= ' ,vicen
814 0 : call icepack_warnings_add(warnstr)
815 0 : write(warnstr,*) subname, 'vsnon= ' ,vsnon
816 0 : call icepack_warnings_add(warnstr)
817 0 : write(warnstr,*) subname, 'viceold= ' ,vice_old
818 0 : call icepack_warnings_add(warnstr)
819 0 : write(warnstr,*) subname, 'vsnoold= ' ,vsno_old
820 0 : call icepack_warnings_add(warnstr)
821 0 : write(warnstr,*) subname, 'melts= ' ,melts
822 0 : call icepack_warnings_add(warnstr)
823 0 : write(warnstr,*) subname, 'meltt= ' ,meltt
824 0 : call icepack_warnings_add(warnstr)
825 0 : write(warnstr,*) subname, 'meltb= ' ,meltb
826 0 : call icepack_warnings_add(warnstr)
827 0 : write(warnstr,*) subname, 'congel= ' ,congel
828 0 : call icepack_warnings_add(warnstr)
829 0 : write(warnstr,*) subname, 'snoice= ' ,snoice
830 0 : call icepack_warnings_add(warnstr)
831 0 : write(warnstr,*) subname, 'aero evap snow= ' ,dhs_evap
832 0 : call icepack_warnings_add(warnstr)
833 0 : write(warnstr,*) subname, 'fsnow= ' ,fsnow
834 0 : call icepack_warnings_add(warnstr)
835 0 : do k = 1, nbtrcr
836 0 : write(warnstr,*) subname, 'NBTRCR value k = ', k
837 0 : call icepack_warnings_add(warnstr)
838 0 : write(warnstr,*) subname, 'aero snowssl (k)= ' ,aerosno0(k,1)
839 0 : call icepack_warnings_add(warnstr)
840 0 : write(warnstr,*) subname, 'aero new snowssl (k)= ',aerosno(k,1)
841 0 : call icepack_warnings_add(warnstr)
842 0 : write(warnstr,*) subname, 'aero snowint (k)= ' ,aerosno0(k,2)
843 0 : call icepack_warnings_add(warnstr)
844 0 : write(warnstr,*) subname, 'aero new snowint(k)= ',aerosno(k,2)
845 0 : call icepack_warnings_add(warnstr)
846 0 : write(warnstr,*) subname, 'flux_bio_atm(k)= ' , flux_bio_atm(k)
847 0 : call icepack_warnings_add(warnstr)
848 0 : write(warnstr,*) subname, 'zbgc_snow(k)= ' ,zbgc_snow(k)
849 0 : call icepack_warnings_add(warnstr)
850 0 : write(warnstr,*) subname, 'zbgc_atm(k)= ' ,zbgc_atm(k)
851 0 : call icepack_warnings_add(warnstr)
852 :
853 0 : do n = 1,2
854 0 : trcrn(bio_index(k)+nblyr+n)=max(trcrn(bio_index(k)+nblyr+n), c0)
855 : enddo
856 : enddo
857 : endif
858 : endif
859 :
860 33835340 : end subroutine update_snow_bgc
861 :
862 : !=======================================================================
863 :
864 : end module icepack_aerosol
865 :
866 : !=======================================================================
|
