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