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