Line data Source code
1 : module modal_aero_calcsize
2 :
3 : ! RCE 07.04.13: Adapted from MIRAGE2 code
4 :
5 : use shr_kind_mod, only: r8 => shr_kind_r8
6 : use spmd_utils, only: masterproc
7 : use physconst, only: pi, rhoh2o, gravit
8 :
9 : use ppgrid, only: pcols, pver
10 : use physics_types, only: physics_state, physics_ptend
11 : use physics_buffer, only: physics_buffer_desc, pbuf_get_index, pbuf_old_tim_idx, pbuf_get_field
12 :
13 : use phys_control, only: phys_getopts
14 : use rad_constituents, only: rad_cnst_get_info, rad_cnst_get_aer_mmr, rad_cnst_get_aer_props, &
15 : rad_cnst_get_mode_props, rad_cnst_get_mode_num
16 :
17 : use cam_logfile, only: iulog
18 : use cam_abortutils, only: endrun
19 : use cam_history, only: addfld, add_default, fieldname_len, horiz_only, outfld
20 : use constituents, only: pcnst, cnst_name
21 : use modal_aero_wateruptake, only: modal_strat_sulfate
22 :
23 : use ref_pres, only: top_lev => clim_modal_aero_top_lev
24 :
25 : #ifdef MODAL_AERO
26 :
27 : ! these are the variables needed for the diagnostic calculation of dry radius
28 : use modal_aero_data, only: ntot_amode, nspec_amode, nspec_max, &
29 : numptr_amode, &
30 : alnsg_amode, &
31 : voltonumbhi_amode, voltonumblo_amode, &
32 : dgnum_amode, dgnumhi_amode, dgnumlo_amode
33 :
34 :
35 : ! these variables are needed for the prognostic calculations to exchange mass
36 : ! between modes
37 : use modal_aero_data, only: numptrcw_amode, mprognum_amode, qqcw_get_field, lmassptrcw_amode, &
38 : lmassptr_amode, modeptr_accum, modeptr_aitken, &
39 : specmw_amode, specdens_amode, voltonumb_amode, &
40 : cnst_name_cw
41 :
42 : use modal_aero_rename, only: lspectooa_renamexf, lspecfrma_renamexf, lspectooc_renamexf, lspecfrmc_renamexf, &
43 : modetoo_renamexf, nspecfrm_renamexf, npair_renamexf, modefrm_renamexf
44 :
45 :
46 : #endif
47 :
48 :
49 : implicit none
50 : private
51 : save
52 :
53 : public modal_aero_calcsize_init, modal_aero_calcsize_sub, modal_aero_calcsize_diag
54 : public :: modal_aero_calcsize_reg
55 :
56 : logical :: do_adjust_default
57 : logical :: do_aitacc_transfer_default
58 :
59 : integer :: dgnum_idx = -1
60 : integer :: hygro_idx = -1
61 : integer :: dryvol_idx = -1
62 : integer :: dryrad_idx = -1
63 : integer :: drymass_idx = -1
64 : integer :: so4dryvol_idx = -1
65 : integer :: naer_idx = -1
66 : integer :: sulfeq_idx = -1
67 :
68 :
69 : !===============================================================================
70 : contains
71 : !===============================================================================
72 :
73 1536 : subroutine modal_aero_calcsize_reg()
74 : use physics_buffer, only: pbuf_add_field, dtype_r8
75 : use rad_constituents, only: rad_cnst_get_info
76 :
77 : integer :: nmodes
78 :
79 1536 : call rad_cnst_get_info(0, nmodes=nmodes)
80 :
81 6144 : call pbuf_add_field('DGNUM', 'global', dtype_r8, (/pcols, pver, nmodes/), dgnum_idx)
82 :
83 6144 : call pbuf_add_field('HYGRO', 'phys_pkg', dtype_r8, (/pcols,pver,nmodes/), hygro_idx)
84 6144 : call pbuf_add_field('DRYVOL', 'phys_pkg', dtype_r8, (/pcols,pver,nmodes/), dryvol_idx)
85 6144 : call pbuf_add_field('DRYRAD', 'phys_pkg', dtype_r8, (/pcols,pver,nmodes/), dryrad_idx)
86 6144 : call pbuf_add_field('DRYMASS', 'phys_pkg', dtype_r8, (/pcols,pver,nmodes/), drymass_idx)
87 6144 : call pbuf_add_field('SO4DRYVOL', 'phys_pkg', dtype_r8, (/pcols,pver,nmodes/), so4dryvol_idx)
88 6144 : call pbuf_add_field('NAER', 'phys_pkg', dtype_r8, (/pcols,pver,nmodes/), naer_idx)
89 :
90 1536 : end subroutine modal_aero_calcsize_reg
91 :
92 : !===============================================================================
93 : !===============================================================================
94 :
95 1536 : subroutine modal_aero_calcsize_init(pbuf2d)
96 1536 : use time_manager, only: is_first_step
97 : use physics_buffer,only: pbuf_set_field
98 :
99 : !-----------------------------------------------------------------------
100 : !
101 : ! Purpose:
102 : ! set do_adjust_default and do_aitacc_transfer_default flags
103 : ! create history fields for column tendencies associated with
104 : ! modal_aero_calcsize
105 : !
106 : ! Author: R. Easter
107 : !
108 : !-----------------------------------------------------------------------
109 :
110 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
111 :
112 : ! local
113 : integer :: ipair, iq
114 : integer :: jac
115 : integer :: lsfrm, lstoo
116 : integer :: n, nacc, nait
117 : logical :: history_aerosol
118 :
119 : character(len=fieldname_len) :: tmpnamea, tmpnameb
120 : character(len=fieldname_len+3) :: fieldname
121 : character(128) :: long_name
122 : character(8) :: unit
123 : !-----------------------------------------------------------------------
124 :
125 1536 : call phys_getopts(history_aerosol_out=history_aerosol)
126 :
127 : ! init entities required for both prescribed and prognostic modes
128 :
129 1536 : if (is_first_step()) then
130 : ! initialize fields in physics buffer
131 768 : call pbuf_set_field(pbuf2d, dgnum_idx, 0.0_r8)
132 : endif
133 :
134 : #ifndef MODAL_AERO
135 : do_adjust_default = .false.
136 : do_aitacc_transfer_default = .false.
137 : #else
138 : ! do_adjust_default allows adjustment to be turned on/off
139 1536 : do_adjust_default = .true.
140 :
141 : ! do_aitacc_transfer_default allows aitken <--> accum mode transfer to be turned on/off
142 : ! *** it can only be true when aitken & accum modes are both present
143 : ! and have prognosed number and diagnosed surface/sigmag
144 1536 : nait = modeptr_aitken
145 1536 : nacc = modeptr_accum
146 1536 : do_aitacc_transfer_default = .false.
147 : if ((modeptr_aitken > 0) .and. &
148 1536 : (modeptr_accum > 0) .and. &
149 : (modeptr_aitken /= modeptr_accum)) then
150 1536 : do_aitacc_transfer_default = .true.
151 1536 : if (mprognum_amode(nait) <= 0) do_aitacc_transfer_default = .false.
152 1536 : if (mprognum_amode(nacc) <= 0) do_aitacc_transfer_default = .false.
153 : end if
154 :
155 : if ( .not. do_adjust_default ) return
156 :
157 : ! define history fields for number-adjust source-sink for all modes
158 7680 : do n = 1, ntot_amode
159 6144 : if (mprognum_amode(n) <= 0) cycle
160 :
161 19968 : do jac = 1, 2
162 12288 : if (jac == 1) then
163 6144 : tmpnamea = cnst_name(numptr_amode(n))
164 : else
165 6144 : tmpnamea = cnst_name_cw(numptrcw_amode(n))
166 : end if
167 12288 : unit = '#/m2/s'
168 12288 : fieldname = trim(tmpnamea) // '_sfcsiz1'
169 12288 : long_name = trim(tmpnamea) // ' calcsize number-adjust column source'
170 12288 : call addfld( fieldname, horiz_only, 'A', unit, long_name )
171 12288 : if (history_aerosol) then
172 0 : call add_default(fieldname, 1, ' ')
173 : end if
174 12288 : if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname
175 :
176 12288 : fieldname = trim(tmpnamea) // '_sfcsiz2'
177 12288 : long_name = trim(tmpnamea) // ' calcsize number-adjust column sink'
178 12288 : call addfld( fieldname, horiz_only, 'A', unit, long_name )
179 12288 : if (history_aerosol) then
180 0 : call add_default(fieldname, 1, ' ')
181 : end if
182 18432 : if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname
183 : end do ! jac = ...
184 : end do ! n = ...
185 :
186 1536 : if ( .not. do_aitacc_transfer_default ) return
187 :
188 : ! check that renaming ipair=1 is aitken-->accum
189 1536 : ipair = 1
190 1536 : if ((modefrm_renamexf(ipair) .ne. nait) .or. &
191 : (modetoo_renamexf(ipair) .ne. nacc)) then
192 : write( 6, '(//2a//)' ) &
193 0 : '*** modal_aero_calcaersize_init error -- ', &
194 0 : 'modefrm/too_renamexf(1) are wrong'
195 0 : call endrun( 'modal_aero_calcaersize_init error' )
196 : end if
197 :
198 : ! define history fields for aitken-accum transfer
199 9216 : do iq = 1, nspecfrm_renamexf(ipair)
200 :
201 : ! jac=1 does interstitial ("_a"); jac=2 does activated ("_c");
202 24576 : do jac = 1, 2
203 :
204 : ! the lspecfrma_renamexf (and lspecfrmc_renamexf) are aitken species
205 : ! the lspectooa_renamexf (and lspectooc_renamexf) are accum species
206 15360 : if (jac .eq. 1) then
207 7680 : lsfrm = lspecfrma_renamexf(iq,ipair)
208 7680 : lstoo = lspectooa_renamexf(iq,ipair)
209 : else
210 7680 : lsfrm = lspecfrmc_renamexf(iq,ipair)
211 7680 : lstoo = lspectooc_renamexf(iq,ipair)
212 : end if
213 15360 : if ((lsfrm <= 0) .or. (lstoo <= 0)) cycle
214 :
215 15360 : if (jac .eq. 1) then
216 7680 : tmpnamea = cnst_name(lsfrm)
217 7680 : tmpnameb = cnst_name(lstoo)
218 : else
219 7680 : tmpnamea = cnst_name_cw(lsfrm)
220 7680 : tmpnameb = cnst_name_cw(lstoo)
221 : end if
222 :
223 15360 : unit = 'kg/m2/s'
224 15360 : if ((tmpnamea(1:3) == 'num') .or. &
225 3072 : (tmpnamea(1:3) == 'NUM')) unit = '#/m2/s'
226 15360 : fieldname = trim(tmpnamea) // '_sfcsiz3'
227 15360 : long_name = trim(tmpnamea) // ' calcsize aitken-to-accum adjust column tendency'
228 15360 : call addfld( fieldname, horiz_only, 'A', unit, long_name )
229 15360 : if (history_aerosol) then
230 0 : call add_default(fieldname, 1, ' ')
231 : end if
232 15360 : if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname
233 :
234 15360 : fieldname = trim(tmpnameb) // '_sfcsiz3'
235 15360 : long_name = trim(tmpnameb) // ' calcsize aitken-to-accum adjust column tendency'
236 15360 : call addfld( fieldname, horiz_only, 'A', unit, long_name )
237 15360 : if (history_aerosol) then
238 0 : call add_default(fieldname, 1, ' ')
239 : end if
240 15360 : if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname
241 :
242 15360 : fieldname = trim(tmpnamea) // '_sfcsiz4'
243 15360 : long_name = trim(tmpnamea) // ' calcsize accum-to-aitken adjust column tendency'
244 15360 : call addfld( fieldname, horiz_only, 'A', unit, long_name )
245 15360 : if (history_aerosol) then
246 0 : call add_default(fieldname, 1, ' ')
247 : end if
248 15360 : if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname
249 :
250 15360 : fieldname = trim(tmpnameb) // '_sfcsiz4'
251 15360 : long_name = trim(tmpnameb) // ' calcsize accum-to-aitken adjust column tendency'
252 15360 : call addfld( fieldname, horiz_only, 'A', unit, long_name )
253 15360 : if (history_aerosol) then
254 0 : call add_default(fieldname, 1, ' ')
255 : end if
256 23040 : if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname
257 :
258 : end do ! jac = ...
259 : end do ! iq = ...
260 :
261 : #endif
262 :
263 1536 : end subroutine modal_aero_calcsize_init
264 :
265 : !===============================================================================
266 :
267 58824 : subroutine modal_aero_calcsize_sub(state, ptend, deltat, pbuf, do_adjust_in, &
268 : do_aitacc_transfer_in)
269 :
270 : !-----------------------------------------------------------------------
271 : !
272 : ! Calculates aerosol size distribution parameters
273 : ! mprognum_amode > 0
274 : ! calculate Dgnum from mass, number, and fixed sigmag
275 : ! mprognum_amode <= 0
276 : ! calculate number from mass, fixed Dgnum, and fixed sigmag
277 : !
278 : ! Also (optionally) adjusts prognostic number to
279 : ! be within bounds determined by mass, Dgnum bounds, and sigma bounds
280 : !
281 : ! Author: R. Easter
282 : !
283 : !-----------------------------------------------------------------------
284 :
285 : ! arguments
286 : type(physics_state), target, intent(in) :: state ! Physics state variables
287 :
288 : type(physics_ptend), target, intent(inout) :: ptend ! indivdual parameterization tendencies
289 :
290 : real(r8), intent(in) :: deltat ! model time-step size (s)
291 : type(physics_buffer_desc), pointer :: pbuf(:) ! physics buffer
292 :
293 : logical, optional :: do_adjust_in
294 : logical, optional :: do_aitacc_transfer_in
295 :
296 : #ifdef MODAL_AERO
297 :
298 : ! local
299 :
300 : logical :: do_adjust
301 : logical :: do_aitacc_transfer
302 :
303 : integer :: lchnk ! chunk identifier
304 : integer :: ncol ! number of columns
305 :
306 58824 : real(r8), pointer :: t(:,:) ! Temperature in Kelvin
307 58824 : real(r8), pointer :: pmid(:,:) ! pressure at model levels (Pa)
308 58824 : real(r8), pointer :: pdel(:,:) ! pressure thickness of levels
309 58824 : real(r8), pointer :: q(:,:,:) ! Tracer MR array
310 :
311 58824 : logical, pointer :: dotend(:) ! flag for doing tendency
312 58824 : real(r8), pointer :: dqdt(:,:,:) ! TMR tendency array
313 :
314 58824 : real(r8), pointer :: dgncur_a(:,:,:)
315 :
316 : integer :: i, icol_diag, iduma, ipair, iq
317 : integer :: ixfer_acc2ait, ixfer_ait2acc
318 : integer :: ixfer_acc2ait_sv(pcols,pver), ixfer_ait2acc_sv(pcols,pver)
319 : integer :: j, jac, jsrflx, k
320 : integer :: l, l1, la, lc, lna, lnc, lsfrm, lstoo
321 : integer :: n, nacc, nait
322 :
323 : integer, save :: idiagaa = 1
324 :
325 : logical :: dotendqqcw(pcnst)
326 117648 : logical :: noxf_acc2ait(nspec_max)
327 :
328 : character(len=fieldname_len) :: tmpnamea, tmpnameb
329 : character(len=fieldname_len+3) :: fieldname
330 :
331 : real(r8), parameter :: third = 1.0_r8/3.0_r8
332 58824 : real(r8), pointer :: fldcw(:,:)
333 : real(r8) :: delnum_a2, delnum_c2 ! work variables
334 : real(r8) :: delnum_a3, delnum_c3, delnum_t3 ! work variables
335 : real(r8) :: deltatinv ! 1/deltat
336 117648 : real(r8) :: dgncur_c(pcols,pver,ntot_amode)
337 : real(r8) :: dgnyy, dgnxx ! dgnumlo/hi of current mode
338 : real(r8) :: dqqcwdt(pcols,pver,pcnst) ! cloudborne TMR tendency array
339 : real(r8) :: drv_a, drv_c, drv_t ! dry volume (cm3/mol_air)
340 : real(r8) :: drv_t0
341 : real(r8) :: drv_a_noxf, drv_c_noxf, drv_t_noxf
342 : real(r8) :: drv_a_acc, drv_c_acc
343 : real(r8) :: drv_a_accsv(pcols,pver), drv_c_accsv(pcols,pver)
344 : real(r8) :: drv_a_aitsv(pcols,pver), drv_c_aitsv(pcols,pver)
345 117648 : real(r8) :: drv_a_sv(pcols,pver,ntot_amode), drv_c_sv(pcols,pver,ntot_amode)
346 : real(r8) :: dryvol_a(pcols,pver) ! interstital aerosol dry
347 : ! volume (cm^3/mol_air)
348 : real(r8) :: dryvol_c(pcols,pver) ! activated aerosol dry volume
349 : real(r8) :: duma, dumb, dumc, dumd ! work variables
350 : real(r8) :: dumfac, dummwdens ! work variables
351 : real(r8) :: frelaxadj ! relaxation factor applied
352 : ! to size bounds
353 : real(r8) :: fracadj ! deltat/tadj
354 : real(r8) :: num_a0, num_c0, num_t0 ! initial number (#/mol_air)
355 : real(r8) :: num_a1, num_c1 ! working number (#/mol_air)
356 : real(r8) :: num_a2, num_c2, num_t2 ! working number (#/mol_air)
357 : real(r8) :: num_a, num_c, num_t ! final number (#/mol_air)
358 : real(r8) :: num_t_noxf
359 : real(r8) :: numbnd ! bounded number
360 : real(r8) :: num_a_acc, num_c_acc
361 : real(r8) :: num_a_accsv(pcols,pver), num_c_accsv(pcols,pver)
362 : real(r8) :: num_a_aitsv(pcols,pver), num_c_aitsv(pcols,pver)
363 117648 : real(r8) :: num_a_sv(pcols,pver,ntot_amode), num_c_sv(pcols,pver,ntot_amode)
364 : real(r8) :: pdel_fac !
365 : real(r8) :: tadj ! adjustment time scale
366 : real(r8) :: tadjinv ! 1/tadj
367 117648 : real(r8) :: v2ncur_a(pcols,pver,ntot_amode)
368 58824 : real(r8) :: v2ncur_c(pcols,pver,ntot_amode)
369 : real(r8) :: v2nyy, v2nxx, v2nzz ! voltonumblo/hi of current mode
370 : real(r8) :: v2nyyrl, v2nxxrl ! relaxed voltonumblo/hi
371 : real(r8) :: xfercoef
372 : real(r8) :: xfercoef_num_acc2ait, xfercoef_vol_acc2ait
373 : real(r8) :: xfercoef_num_ait2acc, xfercoef_vol_ait2acc
374 : real(r8) :: xferfrac_num_acc2ait, xferfrac_vol_acc2ait
375 : real(r8) :: xferfrac_num_ait2acc, xferfrac_vol_ait2acc
376 : real(r8) :: xfertend, xfertend_num(2,2)
377 :
378 : integer, parameter :: nsrflx = 4 ! last dimension of qsrflx
379 : real(r8) :: qsrflx(pcols,pcnst,nsrflx,2)
380 : ! process-specific column tracer tendencies
381 : ! 3rd index --
382 : ! 1="standard" number adjust gain;
383 : ! 2="standard" number adjust loss;
384 : ! 3=aitken-->accum renaming; 4=accum-->aitken)
385 : ! 4th index --
386 : ! 1="a" species; 2="c" species
387 : !-----------------------------------------------------------------------
388 :
389 58824 : if (present(do_adjust_in)) then
390 0 : do_adjust = do_adjust_in
391 : else
392 58824 : do_adjust = do_adjust_default
393 : end if
394 :
395 58824 : if (present(do_aitacc_transfer_in)) then
396 0 : do_aitacc_transfer = do_aitacc_transfer_in
397 : else
398 58824 : do_aitacc_transfer = do_aitacc_transfer_default
399 : end if
400 :
401 58824 : lchnk = state%lchnk
402 58824 : ncol = state%ncol
403 :
404 58824 : t => state%t
405 58824 : pmid => state%pmid
406 58824 : pdel => state%pdel
407 58824 : q => state%q
408 :
409 58824 : dotend => ptend%lq
410 58824 : dqdt => ptend%q
411 :
412 58824 : call pbuf_get_field(pbuf, dgnum_idx, dgncur_a)
413 :
414 58824 : dotendqqcw(:) = .false.
415 58824 : dqqcwdt(:,:,:) = 0.0_r8
416 58824 : qsrflx(:,:,:,:) = 0.0_r8
417 :
418 58824 : nait = modeptr_aitken
419 58824 : nacc = modeptr_accum
420 :
421 58824 : deltatinv = 1.0_r8/(deltat*(1.0_r8 + 1.0e-15_r8))
422 : ! tadj = adjustment time scale for number, surface when they are prognosed
423 : ! currently set to deltat
424 : tadj = deltat
425 58824 : tadj = 86400
426 58824 : tadj = max( tadj, deltat )
427 58824 : tadjinv = 1.0_r8/(tadj*(1.0_r8 + 1.0e-15_r8))
428 58824 : fracadj = deltat*tadjinv
429 58824 : fracadj = max( 0.0_r8, min( 1.0_r8, fracadj ) )
430 :
431 :
432 : !
433 : !
434 : ! the "do 40000" loop does the original (pre jan-2006)
435 : ! number adjustment, one mode at a time
436 : ! this artificially adjusts number when mean particle size is too large
437 : ! or too small
438 : !
439 : !
440 294120 : do n = 1, ntot_amode
441 :
442 :
443 : ! initialize all parameters to the default values for the mode
444 22117824 : do k=top_lev,pver
445 365622624 : do i=1,ncol
446 : ! sgcur_a(i,k,n) = sigmag_amode(n)
447 : ! sgcur_c(i,k,n) = sigmag_amode(n)
448 343504800 : dgncur_a(i,k,n) = dgnum_amode(n)
449 343504800 : dgncur_c(i,k,n) = dgnum_amode(n)
450 343504800 : v2ncur_a(i,k,n) = voltonumb_amode(n)
451 343504800 : v2ncur_c(i,k,n) = voltonumb_amode(n)
452 343504800 : dryvol_a(i,k) = 0.0_r8
453 365387328 : dryvol_c(i,k) = 0.0_r8
454 : end do
455 : end do
456 :
457 : ! compute dry volume mixrats =
458 : ! sum_over_components{ component_mass mixrat / density }
459 1117656 : do l1 = 1, nspec_amode(n)
460 : ! need qmass*dummwdens = (kg/kg-air) * [1/(kg/m3)] = m3/kg-air
461 882360 : dummwdens = 1.0_r8 / specdens_amode(l1,n)
462 882360 : la = lmassptr_amode(l1,n)
463 82941840 : do k=top_lev,pver
464 1371084840 : do i=1,ncol
465 2576286000 : dryvol_a(i,k) = dryvol_a(i,k) &
466 3946488480 : + max(0.0_r8,q(i,k,la))*dummwdens
467 : end do
468 : end do
469 :
470 882360 : fldcw => qqcw_get_field(pbuf,lmassptrcw_amode(l1,n),lchnk)
471 83177136 : do k=top_lev,pver
472 1371084840 : do i=1,ncol
473 2576286000 : dryvol_c(i,k) = dryvol_c(i,k) &
474 3946488480 : + max(0.0_r8,fldcw(i,k))*dummwdens
475 : end do
476 : end do
477 : end do
478 :
479 : ! set "short-hand" number pointers
480 235296 : lna = numptr_amode(n)
481 235296 : lnc = numptrcw_amode(n)
482 235296 : fldcw => qqcw_get_field(pbuf,numptrcw_amode(n),lchnk,.true.)
483 :
484 :
485 : ! go to section for appropriate number/surface diagnosed/prognosed options
486 235296 : if (mprognum_amode(n) <= 0) then
487 :
488 : ! option 1 -- number diagnosed (fixed dgnum and sigmag)
489 : ! compute number tendencies that will bring numbers to their
490 : ! current diagnosed values
491 : !
492 0 : if (lna > 0) then
493 0 : dotend(lna) = .true.
494 0 : do k=top_lev,pver
495 0 : do i=1,ncol
496 0 : dqdt(i,k,lna) = (dryvol_a(i,k)*voltonumb_amode(n) &
497 0 : - q(i,k,lna)) * deltatinv
498 : end do
499 : end do
500 : end if
501 0 : if (lnc > 0) then
502 0 : dotendqqcw(lnc) = .true.
503 0 : do k=top_lev,pver
504 0 : do i=1,ncol
505 0 : dqqcwdt(i,k,lnc) = (dryvol_c(i,k)*voltonumb_amode(n) &
506 0 : - fldcw(i,k)) * deltatinv
507 : end do
508 : end do
509 : end if
510 : else
511 :
512 :
513 : !
514 : ! option 2 -- number prognosed (variable dgnum, fixed sigmag)
515 : ! Compute number tendencies to adjust numbers if they are outside
516 : ! the limits determined by current volume and dgnumlo/hi
517 : ! The interstitial and activated aerosol fractions can, at times,
518 : ! be the lower or upper tail of the "total" distribution. Thus they
519 : ! can be expected to have a greater range of size parameters than
520 : ! what is specified for the total distribution (via dgnumlo/hi)
521 : ! When both the interstitial and activated dry volumes are positive,
522 : ! the adjustment strategy is to (1) adjust the interstitial and activated
523 : ! numbers towards relaxed bounds, then (2) adjust the total/combined
524 : ! number towards the primary bounds.
525 : !
526 : ! note
527 : ! v2nyy = voltonumblo_amode is proportional to dgnumlo**(-3),
528 : ! and produces the maximum allowed number for a given volume
529 : ! v2nxx = voltonumbhi_amode is proportional to dgnumhi**(-3),
530 : ! and produces the minimum allowed number for a given volume
531 : ! v2nxxrl and v2nyyrl are their "relaxed" equivalents.
532 : ! Setting frelaxadj=27=3**3 means that
533 : ! dgnumlo_relaxed = dgnumlo/3 and dgnumhi_relaxed = dgnumhi*3
534 : !
535 : ! if do_aitacc_transfer is .true., then
536 : ! for n=nacc, multiply v2nyy by 1.0e6 to effectively turn off the
537 : ! adjustment when number is too big (size is too small)
538 : ! for n=nait, divide v2nxx by 1.0e6 to effectively turn off the
539 : ! adjustment when number is too small (size is too big)
540 : !OLD however, do not change the v2nyyrl/v2nxxrl so that
541 : !OLD the interstitial<-->activated adjustment is not changed
542 : !NEW also change the v2nyyrl/v2nxxrl so that
543 : !NEW the interstitial<-->activated adjustment is turned off
544 : !
545 : end if
546 235296 : frelaxadj = 27.0_r8
547 235296 : dumfac = exp(4.5_r8*alnsg_amode(n)**2)*pi/6.0_r8
548 235296 : v2nxx = voltonumbhi_amode(n)
549 235296 : v2nyy = voltonumblo_amode(n)
550 235296 : v2nxxrl = v2nxx/frelaxadj
551 235296 : v2nyyrl = v2nyy*frelaxadj
552 235296 : dgnxx = dgnumhi_amode(n)
553 235296 : dgnyy = dgnumlo_amode(n)
554 235296 : if ( do_aitacc_transfer ) then
555 235296 : if (n == nait) v2nxx = v2nxx/1.0e6_r8
556 235296 : if (n == nacc) v2nyy = v2nyy*1.0e6_r8
557 235296 : v2nxxrl = v2nxx/frelaxadj ! NEW
558 235296 : v2nyyrl = v2nyy*frelaxadj ! NEW
559 : end if
560 :
561 235296 : if (do_adjust) then
562 235296 : dotend(lna) = .true.
563 235296 : dotendqqcw(lnc) = .true.
564 : end if
565 :
566 22176648 : do k = top_lev, pver
567 365622624 : do i = 1, ncol
568 :
569 343504800 : drv_a = dryvol_a(i,k)
570 343504800 : num_a0 = q(i,k,lna)
571 343504800 : num_a = max( 0.0_r8, num_a0 )
572 343504800 : drv_c = dryvol_c(i,k)
573 343504800 : num_c0 = fldcw(i,k)
574 343504800 : num_c = max( 0.0_r8, num_c0 )
575 :
576 343504800 : if ( do_adjust) then
577 :
578 : !
579 : ! do number adjustment for interstitial and activated particles
580 : ! adjustments that (1) make numbers non-negative or (2) make numbers
581 : ! zero when volume is zero are applied over time-scale deltat
582 : ! adjustments that bring numbers to within specified bounds are
583 : ! applied over time-scale tadj
584 : !
585 343504800 : if ((drv_a <= 0.0_r8) .and. (drv_c <= 0.0_r8)) then
586 : ! both interstitial and activated volumes are zero
587 : ! adjust both numbers to zero
588 0 : num_a = 0.0_r8
589 0 : dqdt(i,k,lna) = -num_a0*deltatinv
590 0 : num_c = 0.0_r8
591 0 : dqqcwdt(i,k,lnc) = -num_c0*deltatinv
592 343504800 : else if (drv_c <= 0.0_r8) then
593 : ! activated volume is zero, so interstitial number/volume == total/combined
594 : ! apply step 1 and 3, but skip the relaxed adjustment (step 2, see below)
595 286574185 : num_c = 0.0_r8
596 286574185 : dqqcwdt(i,k,lnc) = -num_c0*deltatinv
597 286574185 : num_a1 = num_a
598 286574185 : numbnd = max( drv_a*v2nxx, min( drv_a*v2nyy, num_a1 ) )
599 286574185 : num_a = num_a1 + (numbnd - num_a1)*fracadj
600 286574185 : dqdt(i,k,lna) = (num_a - num_a0)*deltatinv
601 :
602 56930615 : else if (drv_a <= 0.0_r8) then
603 : ! interstitial volume is zero, treat similar to above
604 0 : num_a = 0.0_r8
605 0 : dqdt(i,k,lna) = -num_a0*deltatinv
606 0 : num_c1 = num_c
607 0 : numbnd = max( drv_c*v2nxx, min( drv_c*v2nyy, num_c1 ) )
608 0 : num_c = num_c1 + (numbnd - num_c1)*fracadj
609 0 : dqqcwdt(i,k,lnc) = (num_c - num_c0)*deltatinv
610 : else
611 : ! both volumes are positive
612 : ! apply 3 adjustment steps
613 : ! step1: num_a,c0 --> num_a,c1 forces non-negative values
614 56930615 : num_a1 = num_a
615 56930615 : num_c1 = num_c
616 : ! step2: num_a,c1 --> num_a,c2 applies relaxed bounds to the interstitial
617 : ! and activated number (individually)
618 : ! if only only a or c changes, adjust the other in the opposite direction
619 : ! as much as possible to conserve a+c
620 56930615 : numbnd = max( drv_a*v2nxxrl, min( drv_a*v2nyyrl, num_a1 ) )
621 56930615 : delnum_a2 = (numbnd - num_a1)*fracadj
622 56930615 : num_a2 = num_a1 + delnum_a2
623 56930615 : numbnd = max( drv_c*v2nxxrl, min( drv_c*v2nyyrl, num_c1 ) )
624 56930615 : delnum_c2 = (numbnd - num_c1)*fracadj
625 56930615 : num_c2 = num_c1 + delnum_c2
626 56930615 : if ((delnum_a2 == 0.0_r8) .and. (delnum_c2 /= 0.0_r8)) then
627 : num_a2 = max( drv_a*v2nxxrl, min( drv_a*v2nyyrl, &
628 26788789 : num_a1-delnum_c2 ) )
629 30141826 : else if ((delnum_a2 /= 0.0_r8) .and. (delnum_c2 == 0.0_r8)) then
630 : num_c2 = max( drv_c*v2nxxrl, min( drv_c*v2nyyrl, &
631 4 : num_c1-delnum_a2 ) )
632 : end if
633 : ! step3: num_a,c2 --> num_a,c3 applies stricter bounds to the
634 : ! combined/total number
635 56930615 : drv_t = drv_a + drv_c
636 56930615 : num_t2 = num_a2 + num_c2
637 56930615 : delnum_a3 = 0.0_r8
638 56930615 : delnum_c3 = 0.0_r8
639 56930615 : if (num_t2 < drv_t*v2nxx) then
640 5831203 : delnum_t3 = (drv_t*v2nxx - num_t2)*fracadj
641 : ! if you are here then (num_a2 < drv_a*v2nxx) and/or
642 : ! (num_c2 < drv_c*v2nxx) must be true
643 5831203 : if ((num_a2 < drv_a*v2nxx) .and. (num_c2 < drv_c*v2nxx)) then
644 5830737 : delnum_a3 = delnum_t3*(num_a2/num_t2)
645 5830737 : delnum_c3 = delnum_t3*(num_c2/num_t2)
646 466 : else if (num_c2 < drv_c*v2nxx) then
647 : delnum_c3 = delnum_t3
648 33 : else if (num_a2 < drv_a*v2nxx) then
649 33 : delnum_a3 = delnum_t3
650 : end if
651 51099412 : else if (num_t2 > drv_t*v2nyy) then
652 2995320 : delnum_t3 = (drv_t*v2nyy - num_t2)*fracadj
653 : ! if you are here then (num_a2 > drv_a*v2nyy) and/or
654 : ! (num_c2 > drv_c*v2nyy) must be true
655 2995320 : if ((num_a2 > drv_a*v2nyy) .and. (num_c2 > drv_c*v2nyy)) then
656 1937520 : delnum_a3 = delnum_t3*(num_a2/num_t2)
657 1937520 : delnum_c3 = delnum_t3*(num_c2/num_t2)
658 1057800 : else if (num_c2 > drv_c*v2nyy) then
659 : delnum_c3 = delnum_t3
660 1025092 : else if (num_a2 > drv_a*v2nyy) then
661 1025092 : delnum_a3 = delnum_t3
662 : end if
663 : end if
664 56930615 : num_a = num_a2 + delnum_a3
665 56930615 : dqdt(i,k,lna) = (num_a - num_a0)*deltatinv
666 56930615 : num_c = num_c2 + delnum_c3
667 56930615 : dqqcwdt(i,k,lnc) = (num_c - num_c0)*deltatinv
668 : end if
669 :
670 : end if ! do_adjust
671 :
672 : !
673 : ! now compute current dgn and v2n
674 : !
675 343504800 : if (drv_a > 0.0_r8) then
676 343504800 : if (num_a <= drv_a*v2nxx) then
677 37732270 : dgncur_a(i,k,n) = dgnxx
678 37732270 : v2ncur_a(i,k,n) = v2nxx
679 305772530 : else if (num_a >= drv_a*v2nyy) then
680 25259708 : dgncur_a(i,k,n) = dgnyy
681 25259708 : v2ncur_a(i,k,n) = v2nyy
682 : else
683 280512822 : dgncur_a(i,k,n) = (drv_a/(dumfac*num_a))**third
684 280512822 : v2ncur_a(i,k,n) = num_a/drv_a
685 : end if
686 : end if
687 343504800 : pdel_fac = pdel(i,k)/gravit ! = rho*dz
688 343504800 : jac = 1
689 343504800 : qsrflx(i,lna,1,jac) = qsrflx(i,lna,1,jac) + max(0.0_r8,dqdt(i,k,lna))*pdel_fac
690 343504800 : qsrflx(i,lna,2,jac) = qsrflx(i,lna,2,jac) + min(0.0_r8,dqdt(i,k,lna))*pdel_fac
691 :
692 343504800 : if (drv_c > 0.0_r8) then
693 56930615 : if (num_c <= drv_c*v2nxx) then
694 27195320 : dgncur_c(i,k,n) = dgnumhi_amode(n)
695 27195320 : v2ncur_c(i,k,n) = v2nxx
696 29735295 : else if (num_c >= drv_c*v2nyy) then
697 3207004 : dgncur_c(i,k,n) = dgnumlo_amode(n)
698 3207004 : v2ncur_c(i,k,n) = v2nyy
699 : else
700 26528291 : dgncur_c(i,k,n) = (drv_c/(dumfac*num_c))**third
701 26528291 : v2ncur_c(i,k,n) = num_c/drv_c
702 : end if
703 : end if
704 343504800 : jac = 2
705 343504800 : qsrflx(i,lnc,1,jac) = qsrflx(i,lnc,1,jac) + max(0.0_r8,dqqcwdt(i,k,lnc))*pdel_fac
706 343504800 : qsrflx(i,lnc,2,jac) = qsrflx(i,lnc,2,jac) + min(0.0_r8,dqqcwdt(i,k,lnc))*pdel_fac
707 :
708 :
709 : ! save number and dryvol for aitken <--> accum renaming
710 343504800 : if ( do_aitacc_transfer ) then
711 343504800 : if (n == nait) then
712 85876200 : drv_a_aitsv(i,k) = drv_a
713 85876200 : num_a_aitsv(i,k) = num_a
714 85876200 : drv_c_aitsv(i,k) = drv_c
715 85876200 : num_c_aitsv(i,k) = num_c
716 257628600 : else if (n == nacc) then
717 85876200 : drv_a_accsv(i,k) = drv_a
718 85876200 : num_a_accsv(i,k) = num_a
719 85876200 : drv_c_accsv(i,k) = drv_c
720 85876200 : num_c_accsv(i,k) = num_c
721 : end if
722 : end if
723 343504800 : drv_a_sv(i,k,n) = drv_a
724 343504800 : num_a_sv(i,k,n) = num_a
725 343504800 : drv_c_sv(i,k,n) = drv_c
726 365387328 : num_c_sv(i,k,n) = num_c
727 :
728 : end do
729 : end do
730 :
731 :
732 : !
733 : ! option 3 -- number and surface prognosed (variable dgnum and sigmag)
734 : ! this is not implemented
735 : !
736 : end do ! do n = 1, ntot_amode
737 :
738 :
739 : !
740 : !
741 : ! the following section (from here to label 49000)
742 : ! does aitken <--> accum mode transfer
743 : !
744 : ! when the aitken mode mean size is too big, the largest
745 : ! aitken particles are transferred into the accum mode
746 : ! to reduce the aitken mode mean size
747 : ! when the accum mode mean size is too small, the smallest
748 : ! accum particles are transferred into the aitken mode
749 : ! to increase the accum mode mean size
750 : !
751 : !
752 58824 : ixfer_ait2acc_sv(:,:) = 0
753 58824 : ixfer_acc2ait_sv(:,:) = 0
754 58824 : if ( do_aitacc_transfer ) then
755 :
756 : ! old - on time first step, npair_renamexf will be <= 0,
757 : ! in which case need to do modal_aero_rename_init
758 : ! new - init is now done through chem_init and things below it
759 58824 : if (npair_renamexf .le. 0) then
760 0 : npair_renamexf = 0
761 : ! call modal_aero_rename_init
762 : if (npair_renamexf .le. 0) then
763 : write( 6, '(//a//)' ) &
764 0 : '*** modal_aero_calcaersize_sub error -- npair_renamexf <= 0'
765 0 : call endrun( 'modal_aero_calcaersize_sub error' )
766 : end if
767 : end if
768 :
769 : ! check that renaming ipair=1 is aitken-->accum
770 58824 : ipair = 1
771 58824 : if ((modefrm_renamexf(ipair) .ne. nait) .or. &
772 : (modetoo_renamexf(ipair) .ne. nacc)) then
773 : write( 6, '(//2a//)' ) &
774 0 : '*** modal_aero_calcaersize_sub error -- ', &
775 0 : 'modefrm/too_renamexf(1) are wrong'
776 0 : call endrun( 'modal_aero_calcaersize_sub error' )
777 : end if
778 :
779 : ! set dotend() for species that will be transferred
780 352944 : do iq = 1, nspecfrm_renamexf(ipair)
781 294120 : lsfrm = lspecfrma_renamexf(iq,ipair)
782 294120 : lstoo = lspectooa_renamexf(iq,ipair)
783 294120 : if ((lsfrm > 0) .and. (lstoo > 0)) then
784 294120 : dotend(lsfrm) = .true.
785 294120 : dotend(lstoo) = .true.
786 : end if
787 294120 : lsfrm = lspecfrmc_renamexf(iq,ipair)
788 294120 : lstoo = lspectooc_renamexf(iq,ipair)
789 352944 : if ((lsfrm > 0) .and. (lstoo > 0)) then
790 294120 : dotendqqcw(lsfrm) = .true.
791 294120 : dotendqqcw(lstoo) = .true.
792 : end if
793 : end do
794 :
795 : ! identify accum species cannot be transferred to aitken mode
796 411768 : noxf_acc2ait(:) = .true.
797 411768 : do l1 = 1, nspec_amode(nacc)
798 352944 : la = lmassptr_amode(l1,nacc)
799 2176488 : do iq = 1, nspecfrm_renamexf(ipair)
800 2117664 : if (lspectooa_renamexf(iq,ipair) == la) then
801 235296 : noxf_acc2ait(l1) = .false.
802 : end if
803 : end do
804 : end do
805 :
806 : ! v2nzz is voltonumb at the "geometrically-defined" mid-point
807 : ! between the aitken and accum modes
808 58824 : v2nzz = sqrt(voltonumb_amode(nait)*voltonumb_amode(nacc))
809 :
810 : ! loop over columns and levels
811 5529456 : do k = top_lev, pver
812 91405656 : do i = 1, ncol
813 :
814 85876200 : pdel_fac = pdel(i,k)/gravit ! = rho*dz
815 85876200 : xfertend_num(:,:) = 0.0_r8
816 :
817 : ! compute aitken --> accum transfer rates
818 85876200 : ixfer_ait2acc = 0
819 85876200 : xfercoef_num_ait2acc = 0.0_r8
820 85876200 : xfercoef_vol_ait2acc = 0.0_r8
821 :
822 85876200 : drv_t = drv_a_aitsv(i,k) + drv_c_aitsv(i,k)
823 85876200 : num_t = num_a_aitsv(i,k) + num_c_aitsv(i,k)
824 85876200 : if (drv_t > 0.0_r8) then
825 85876200 : if (num_t < drv_t*v2nzz) then
826 4936248 : ixfer_ait2acc = 1
827 4936248 : if (num_t < drv_t*voltonumb_amode(nacc)) then
828 : xferfrac_num_ait2acc = 1.0_r8
829 : xferfrac_vol_ait2acc = 1.0_r8
830 : else
831 : xferfrac_vol_ait2acc = ((num_t/drv_t) - v2nzz)/ &
832 4936185 : (voltonumb_amode(nacc) - v2nzz)
833 : xferfrac_num_ait2acc = xferfrac_vol_ait2acc* &
834 4936185 : (drv_t*voltonumb_amode(nacc)/num_t)
835 4936185 : if ((xferfrac_num_ait2acc <= 0.0_r8) .or. &
836 : (xferfrac_vol_ait2acc <= 0.0_r8)) then
837 : xferfrac_num_ait2acc = 0.0_r8
838 : xferfrac_vol_ait2acc = 0.0_r8
839 4936185 : else if ((xferfrac_num_ait2acc >= 1.0_r8) .or. &
840 : (xferfrac_vol_ait2acc >= 1.0_r8)) then
841 0 : xferfrac_num_ait2acc = 1.0_r8
842 0 : xferfrac_vol_ait2acc = 1.0_r8
843 : end if
844 : end if
845 4936248 : xfercoef_num_ait2acc = xferfrac_num_ait2acc*tadjinv
846 4936248 : xfercoef_vol_ait2acc = xferfrac_vol_ait2acc*tadjinv
847 4936248 : xfertend_num(1,1) = num_a_aitsv(i,k)*xfercoef_num_ait2acc
848 4936248 : xfertend_num(1,2) = num_c_aitsv(i,k)*xfercoef_num_ait2acc
849 : end if
850 : end if
851 :
852 : ! compute accum --> aitken transfer rates
853 : ! accum may have some species (seasalt, dust, poa, lll) that are
854 : ! not in aitken mode
855 : ! so first divide the accum drv & num into not-transferred (noxf) species
856 : ! and transferred species, and use the transferred-species
857 : ! portion in what follows
858 85876200 : ixfer_acc2ait = 0
859 85876200 : xfercoef_num_acc2ait = 0.0_r8
860 85876200 : xfercoef_vol_acc2ait = 0.0_r8
861 :
862 85876200 : drv_t = drv_a_accsv(i,k) + drv_c_accsv(i,k)
863 85876200 : num_t = num_a_accsv(i,k) + num_c_accsv(i,k)
864 85876200 : drv_a_noxf = 0.0_r8
865 85876200 : drv_c_noxf = 0.0_r8
866 85876200 : if (drv_t > 0.0_r8) then
867 85876200 : if (num_t > drv_t*v2nzz) then
868 56661220 : do l1 = 1, nspec_amode(nacc)
869 :
870 56661220 : if ( noxf_acc2ait(l1) ) then
871 : ! need qmass*dummwdens = (kg/kg-air) * [1/(kg/m3)] = m3/kg-air
872 16188920 : dummwdens = 1.0_r8 / specdens_amode(l1,nacc)
873 16188920 : la = lmassptr_amode(l1,nacc)
874 : drv_a_noxf = drv_a_noxf &
875 16188920 : + max(0.0_r8,q(i,k,la))*dummwdens
876 16188920 : lc = lmassptrcw_amode(l1,nacc)
877 :
878 16188920 : fldcw => qqcw_get_field(pbuf,lmassptrcw_amode(l1,nacc),lchnk)
879 : drv_c_noxf = drv_c_noxf &
880 16188920 : + max(0.0_r8,fldcw(i,k))*dummwdens
881 : end if
882 : end do
883 8094460 : drv_t_noxf = drv_a_noxf + drv_c_noxf
884 8094460 : num_t_noxf = drv_t_noxf*voltonumblo_amode(nacc)
885 8094460 : num_t0 = num_t
886 8094460 : drv_t0 = drv_t
887 8094460 : num_t = max( 0.0_r8, num_t - num_t_noxf )
888 8094460 : drv_t = max( 0.0_r8, drv_t - drv_t_noxf )
889 : end if
890 : end if
891 :
892 85876200 : if (drv_t > 0.0_r8) then
893 85876200 : if (num_t > drv_t*v2nzz) then
894 8094460 : ixfer_acc2ait = 1
895 8094460 : if (num_t > drv_t*voltonumb_amode(nait)) then
896 : xferfrac_num_acc2ait = 1.0_r8
897 : xferfrac_vol_acc2ait = 1.0_r8
898 : else
899 : xferfrac_vol_acc2ait = ((num_t/drv_t) - v2nzz)/ &
900 7344162 : (voltonumb_amode(nait) - v2nzz)
901 : xferfrac_num_acc2ait = xferfrac_vol_acc2ait* &
902 7344162 : (drv_t*voltonumb_amode(nait)/num_t)
903 7344162 : if ((xferfrac_num_acc2ait <= 0.0_r8) .or. &
904 : (xferfrac_vol_acc2ait <= 0.0_r8)) then
905 : xferfrac_num_acc2ait = 0.0_r8
906 : xferfrac_vol_acc2ait = 0.0_r8
907 7344162 : else if ((xferfrac_num_acc2ait >= 1.0_r8) .or. &
908 : (xferfrac_vol_acc2ait >= 1.0_r8)) then
909 0 : xferfrac_num_acc2ait = 1.0_r8
910 0 : xferfrac_vol_acc2ait = 1.0_r8
911 : end if
912 : end if
913 8094460 : duma = 1.0e-37_r8
914 : xferfrac_num_acc2ait = xferfrac_num_acc2ait* &
915 8094460 : num_t/max( duma, num_t0 )
916 8094460 : xfercoef_num_acc2ait = xferfrac_num_acc2ait*tadjinv
917 8094460 : xfercoef_vol_acc2ait = xferfrac_vol_acc2ait*tadjinv
918 8094460 : xfertend_num(2,1) = num_a_accsv(i,k)*xfercoef_num_acc2ait
919 8094460 : xfertend_num(2,2) = num_c_accsv(i,k)*xfercoef_num_acc2ait
920 : end if
921 : end if
922 :
923 : ! jump to end-of-loop if no transfer is needed at current i,k
924 91346832 : if (ixfer_ait2acc+ixfer_acc2ait > 0) then
925 13030708 : ixfer_ait2acc_sv(i,k) = ixfer_ait2acc
926 13030708 : ixfer_acc2ait_sv(i,k) = ixfer_acc2ait
927 :
928 : !
929 : ! compute new dgncur & v2ncur for aitken & accum modes
930 : !
931 : ! currently inactive
932 39092124 : do n = nait, nacc, (nacc-nait)
933 26061416 : if (n .eq. nait) then
934 13030708 : duma = (xfertend_num(1,1) - xfertend_num(2,1))*deltat
935 13030708 : num_a = max( 0.0_r8, num_a_aitsv(i,k) - duma )
936 13030708 : num_a_acc = max( 0.0_r8, num_a_accsv(i,k) + duma )
937 : duma = (drv_a_aitsv(i,k)*xfercoef_vol_ait2acc - &
938 13030708 : (drv_a_accsv(i,k)-drv_a_noxf)*xfercoef_vol_acc2ait)*deltat
939 13030708 : drv_a = max( 0.0_r8, drv_a_aitsv(i,k) - duma )
940 13030708 : drv_a_acc = max( 0.0_r8, drv_a_accsv(i,k) + duma )
941 13030708 : duma = (xfertend_num(1,2) - xfertend_num(2,2))*deltat
942 13030708 : num_c = max( 0.0_r8, num_c_aitsv(i,k) - duma )
943 13030708 : num_c_acc = max( 0.0_r8, num_c_accsv(i,k) + duma )
944 : duma = (drv_c_aitsv(i,k)*xfercoef_vol_ait2acc - &
945 13030708 : (drv_c_accsv(i,k)-drv_c_noxf)*xfercoef_vol_acc2ait)*deltat
946 13030708 : drv_c = max( 0.0_r8, drv_c_aitsv(i,k) - duma )
947 13030708 : drv_c_acc = max( 0.0_r8, drv_c_accsv(i,k) + duma )
948 : else
949 : num_a = num_a_acc
950 : drv_a = drv_a_acc
951 : num_c = num_c_acc
952 : drv_c = drv_c_acc
953 : end if
954 :
955 26061416 : if (drv_a > 0.0_r8) then
956 26061416 : if (num_a <= drv_a*voltonumbhi_amode(n)) then
957 4926141 : dgncur_a(i,k,n) = dgnumhi_amode(n)
958 4926141 : v2ncur_a(i,k,n) = voltonumbhi_amode(n)
959 21135275 : else if (num_a >= drv_a*voltonumblo_amode(n)) then
960 9165983 : dgncur_a(i,k,n) = dgnumlo_amode(n)
961 9165983 : v2ncur_a(i,k,n) = voltonumblo_amode(n)
962 : else
963 11969292 : dgncur_a(i,k,n) = (drv_a/(dumfac*num_a))**third
964 11969292 : v2ncur_a(i,k,n) = num_a/drv_a
965 : end if
966 : else
967 0 : dgncur_a(i,k,n) = dgnum_amode(n)
968 0 : v2ncur_a(i,k,n) = voltonumb_amode(n)
969 : end if
970 :
971 39092124 : if (drv_c > 0.0_r8) then
972 264318 : if (num_c <= drv_c*voltonumbhi_amode(n)) then
973 168528 : dgncur_c(i,k,n) = dgnumhi_amode(n)
974 168528 : v2ncur_c(i,k,n) = voltonumbhi_amode(n)
975 95790 : else if (num_c >= drv_c*voltonumblo_amode(n)) then
976 0 : dgncur_c(i,k,n) = dgnumlo_amode(n)
977 0 : v2ncur_c(i,k,n) = voltonumblo_amode(n)
978 : else
979 95790 : dgncur_c(i,k,n) = (drv_c/(dumfac*num_c))**third
980 95790 : v2ncur_c(i,k,n) = num_c/drv_c
981 : end if
982 : else
983 25797098 : dgncur_c(i,k,n) = dgnum_amode(n)
984 25797098 : v2ncur_c(i,k,n) = voltonumb_amode(n)
985 : end if
986 :
987 : end do
988 :
989 :
990 : !
991 : ! compute tendency amounts for aitken <--> accum transfer
992 : !
993 :
994 13030708 : if ( masterproc ) then
995 16338 : if (idiagaa > 0) then
996 6 : do j = 1, 2
997 26 : do iq = 1, nspecfrm_renamexf(ipair)
998 64 : do jac = 1, 2
999 40 : if (j .eq. 1) then
1000 20 : if (jac .eq. 1) then
1001 10 : lsfrm = lspecfrma_renamexf(iq,ipair)
1002 10 : lstoo = lspectooa_renamexf(iq,ipair)
1003 : else
1004 10 : lsfrm = lspecfrmc_renamexf(iq,ipair)
1005 10 : lstoo = lspectooc_renamexf(iq,ipair)
1006 : end if
1007 : else
1008 20 : if (jac .eq. 1) then
1009 10 : lsfrm = lspectooa_renamexf(iq,ipair)
1010 10 : lstoo = lspecfrma_renamexf(iq,ipair)
1011 : else
1012 10 : lsfrm = lspectooc_renamexf(iq,ipair)
1013 10 : lstoo = lspecfrmc_renamexf(iq,ipair)
1014 : end if
1015 : end if
1016 40 : write( 6, '(a,3i3,2i4)' ) 'calcsize j,iq,jac, lsfrm,lstoo', &
1017 100 : j,iq,jac, lsfrm,lstoo
1018 : end do
1019 : end do
1020 : end do
1021 : end if
1022 : end if
1023 13030708 : idiagaa = -1
1024 :
1025 :
1026 : ! j=1 does aitken-->accum; j=2 does accum-->aitken
1027 39092124 : do j = 1, 2
1028 :
1029 26061416 : if ((j .eq. 1 .and. ixfer_ait2acc > 0) .or. &
1030 13030708 : (j .eq. 2 .and. ixfer_acc2ait > 0)) then
1031 :
1032 13030708 : jsrflx = j+2
1033 13030708 : if (j .eq. 1) then
1034 : xfercoef = xfercoef_vol_ait2acc
1035 : else
1036 8094460 : xfercoef = xfercoef_vol_acc2ait
1037 : end if
1038 :
1039 78184248 : do iq = 1, nspecfrm_renamexf(ipair)
1040 :
1041 : ! jac=1 does interstitial ("_a"); jac=2 does activated ("_c");
1042 208491328 : do jac = 1, 2
1043 :
1044 : ! the lspecfrma_renamexf (and lspecfrmc_renamexf) are aitken species
1045 : ! the lspectooa_renamexf (and lspectooc_renamexf) are accum species
1046 : ! for j=1, want lsfrm=aitken species, lstoo=accum species
1047 : ! for j=2, want lsfrm=accum species, lstoo=aitken species
1048 130307080 : if (j .eq. 1) then
1049 49362480 : if (jac .eq. 1) then
1050 24681240 : lsfrm = lspecfrma_renamexf(iq,ipair)
1051 24681240 : lstoo = lspectooa_renamexf(iq,ipair)
1052 : else
1053 24681240 : lsfrm = lspecfrmc_renamexf(iq,ipair)
1054 24681240 : lstoo = lspectooc_renamexf(iq,ipair)
1055 : end if
1056 : else
1057 80944600 : if (jac .eq. 1) then
1058 40472300 : lsfrm = lspectooa_renamexf(iq,ipair)
1059 40472300 : lstoo = lspecfrma_renamexf(iq,ipair)
1060 : else
1061 40472300 : lsfrm = lspectooc_renamexf(iq,ipair)
1062 40472300 : lstoo = lspecfrmc_renamexf(iq,ipair)
1063 : end if
1064 : end if
1065 :
1066 195460620 : if ((lsfrm > 0) .and. (lstoo > 0)) then
1067 130307080 : if (jac .eq. 1) then
1068 65153540 : if (iq .eq. 1) then
1069 13030708 : xfertend = xfertend_num(j,jac)
1070 : else
1071 52122832 : xfertend = max(0.0_r8,q(i,k,lsfrm))*xfercoef
1072 : end if
1073 65153540 : dqdt(i,k,lsfrm) = dqdt(i,k,lsfrm) - xfertend
1074 65153540 : dqdt(i,k,lstoo) = dqdt(i,k,lstoo) + xfertend
1075 : else
1076 65153540 : if (iq .eq. 1) then
1077 13030708 : xfertend = xfertend_num(j,jac)
1078 : else
1079 52122832 : fldcw => qqcw_get_field(pbuf,lsfrm,lchnk)
1080 52122832 : xfertend = max(0.0_r8,fldcw(i,k))*xfercoef
1081 : end if
1082 65153540 : dqqcwdt(i,k,lsfrm) = dqqcwdt(i,k,lsfrm) - xfertend
1083 65153540 : dqqcwdt(i,k,lstoo) = dqqcwdt(i,k,lstoo) + xfertend
1084 : end if
1085 130307080 : qsrflx(i,lsfrm,jsrflx,jac) = qsrflx(i,lsfrm,jsrflx,jac) - xfertend*pdel_fac
1086 130307080 : qsrflx(i,lstoo,jsrflx,jac) = qsrflx(i,lstoo,jsrflx,jac) + xfertend*pdel_fac
1087 : end if
1088 :
1089 : end do
1090 : end do
1091 : end if
1092 : end do
1093 :
1094 : end if
1095 : end do
1096 : end do
1097 :
1098 :
1099 : end if ! do_aitacc_transfer
1100 58824 : lsfrm = -123456789 ! executable statement for debugging
1101 :
1102 :
1103 : !
1104 : ! apply tendencies to cloud-borne species MRs
1105 : !
1106 2470608 : do l = 1, pcnst
1107 2411784 : lc = l
1108 2470608 : if ( lc>0 .and. dotendqqcw(lc) ) then
1109 705888 : fldcw=> qqcw_get_field(pbuf,l,lchnk)
1110 66353472 : do k = top_lev, pver
1111 1096867872 : do i = 1, ncol
1112 0 : fldcw(i,k) = max( 0.0_r8, &
1113 1096161984 : (fldcw(i,k) + dqqcwdt(i,k,lc)*deltat) )
1114 : end do
1115 : end do
1116 : end if
1117 : end do
1118 :
1119 : !
1120 : ! do outfld calls
1121 : !
1122 :
1123 : ! history fields for number-adjust source-sink for all modes
1124 58824 : if ( .not. do_adjust ) return
1125 :
1126 294120 : do n = 1, ntot_amode
1127 235296 : if (mprognum_amode(n) <= 0) cycle
1128 :
1129 764712 : do jac = 1, 2
1130 470592 : if (jac == 1) then
1131 235296 : l = numptr_amode(n)
1132 235296 : tmpnamea = cnst_name(l)
1133 : else
1134 235296 : l = numptrcw_amode(n)
1135 235296 : tmpnamea = cnst_name_cw(l)
1136 : end if
1137 470592 : fieldname = trim(tmpnamea) // '_sfcsiz1'
1138 470592 : call outfld( fieldname, qsrflx(:,l,1,jac), pcols, lchnk)
1139 :
1140 470592 : fieldname = trim(tmpnamea) // '_sfcsiz2'
1141 705888 : call outfld( fieldname, qsrflx(:,l,2,jac), pcols, lchnk)
1142 : end do ! jac = ...
1143 :
1144 : end do ! n = ...
1145 :
1146 :
1147 : ! history fields for aitken-accum transfer
1148 58824 : if ( .not. do_aitacc_transfer ) return
1149 :
1150 352944 : do iq = 1, nspecfrm_renamexf(ipair)
1151 :
1152 : ! jac=1 does interstitial ("_a"); jac=2 does activated ("_c");
1153 941184 : do jac = 1, 2
1154 :
1155 : ! the lspecfrma_renamexf (and lspecfrmc_renamexf) are aitken species
1156 : ! the lspectooa_renamexf (and lspectooc_renamexf) are accum species
1157 588240 : if (jac .eq. 1) then
1158 294120 : lsfrm = lspecfrma_renamexf(iq,ipair)
1159 294120 : lstoo = lspectooa_renamexf(iq,ipair)
1160 : else
1161 294120 : lsfrm = lspecfrmc_renamexf(iq,ipair)
1162 294120 : lstoo = lspectooc_renamexf(iq,ipair)
1163 : end if
1164 588240 : if ((lsfrm <= 0) .or. (lstoo <= 0)) cycle
1165 :
1166 588240 : if (jac .eq. 1) then
1167 294120 : tmpnamea = cnst_name(lsfrm)
1168 294120 : tmpnameb = cnst_name(lstoo)
1169 : else
1170 294120 : tmpnamea = cnst_name_cw(lsfrm)
1171 294120 : tmpnameb = cnst_name_cw(lstoo)
1172 : end if
1173 : if ((lsfrm <= 0) .or. (lstoo <= 0)) cycle
1174 :
1175 588240 : fieldname = trim(tmpnamea) // '_sfcsiz3'
1176 588240 : call outfld( fieldname, qsrflx(:,lsfrm,3,jac), pcols, lchnk)
1177 :
1178 588240 : fieldname = trim(tmpnameb) // '_sfcsiz3'
1179 588240 : call outfld( fieldname, qsrflx(:,lstoo,3,jac), pcols, lchnk)
1180 :
1181 588240 : fieldname = trim(tmpnamea) // '_sfcsiz4'
1182 588240 : call outfld( fieldname, qsrflx(:,lsfrm,4,jac), pcols, lchnk)
1183 :
1184 588240 : fieldname = trim(tmpnameb) // '_sfcsiz4'
1185 882360 : call outfld( fieldname, qsrflx(:,lstoo,4,jac), pcols, lchnk)
1186 :
1187 : end do ! jac = ...
1188 : end do ! iq = ...
1189 :
1190 58824 : call modal_aero_calcdry(state, pbuf)
1191 :
1192 : #endif
1193 :
1194 60360 : end subroutine modal_aero_calcsize_sub
1195 :
1196 :
1197 : !----------------------------------------------------------------------
1198 :
1199 :
1200 6192 : subroutine modal_aero_calcsize_diag(state, pbuf, list_idx_in, dgnum_m, &
1201 : hygro_m, dryvol_m, dryrad_m, drymass_m, so4dryvol_m, naer_m)
1202 :
1203 : !-----------------------------------------------------------------------
1204 : !
1205 : ! Calculate aerosol size distribution parameters
1206 : !
1207 : ! ***N.B.*** DGNUM for the modes in the climate list are put directly into
1208 : ! the physics buffer. For diagnostic list calculations use the
1209 : ! optional list_idx and dgnum args.
1210 : !-----------------------------------------------------------------------
1211 :
1212 : ! arguments
1213 : type(physics_state), intent(in), target :: state ! Physics state variables
1214 : type(physics_buffer_desc), pointer :: pbuf(:) ! physics buffer
1215 :
1216 : integer, optional, intent(in) :: list_idx_in ! diagnostic list index
1217 : real(r8), optional, pointer :: dgnum_m(:,:,:) ! interstital aerosol dry number mode radius (m)
1218 : real(r8), optional, pointer :: hygro_m(:,:,:)
1219 : real(r8), optional, pointer :: dryvol_m(:,:,:)
1220 : real(r8), optional, pointer :: dryrad_m(:,:,:)
1221 : real(r8), optional, pointer :: drymass_m(:,:,:)
1222 : real(r8), optional, pointer :: so4dryvol_m(:,:,:)
1223 : real(r8), optional, pointer :: naer_m(:,:,:)
1224 :
1225 :
1226 : ! local
1227 : integer :: i, k, l1, n
1228 : integer :: lchnk, ncol
1229 : integer :: list_idx, stat
1230 : integer :: nmodes
1231 : integer :: nspec
1232 :
1233 3096 : real(r8), pointer :: dgncur_a(:,:) ! (pcols,pver)
1234 :
1235 :
1236 : real(r8), parameter :: third = 1.0_r8/3.0_r8
1237 :
1238 3096 : real(r8), pointer :: mode_num(:,:) ! mode number mixing ratio
1239 3096 : real(r8), pointer :: specmmr(:,:) ! specie mmr
1240 : real(r8) :: specdens ! specie density
1241 :
1242 : real(r8) :: dryvol_a(pcols,pver) ! interstital aerosol dry volume (cm^3/mol_air)
1243 :
1244 : real(r8) :: dgnum, dgnumhi, dgnumlo
1245 : real(r8) :: dgnyy, dgnxx ! dgnumlo/hi of current mode
1246 : real(r8) :: drv_a ! dry volume (cm3/mol_air)
1247 : real(r8) :: dumfac, dummwdens ! work variables
1248 : real(r8) :: num_a0 ! initial number (#/mol_air)
1249 : real(r8) :: num_a ! final number (#/mol_air)
1250 : real(r8) :: voltonumbhi, voltonumblo
1251 : real(r8) :: v2nyy, v2nxx ! voltonumblo/hi of current mode
1252 : real(r8) :: sigmag, alnsg
1253 : !-----------------------------------------------------------------------
1254 :
1255 3096 : lchnk = state%lchnk
1256 3096 : ncol = state%ncol
1257 :
1258 3096 : list_idx = 0 ! climate list by default
1259 0 : if (present(list_idx_in)) list_idx = list_idx_in
1260 :
1261 3096 : call rad_cnst_get_info(list_idx, nmodes=nmodes)
1262 :
1263 3096 : if (list_idx /= 0) then
1264 0 : if (.not. present(dgnum_m)) then
1265 : call endrun('modal_aero_calcsize_diag called for'// &
1266 0 : 'diagnostic list but dgnum_m pointer not present')
1267 : end if
1268 0 : if (.not. associated(dgnum_m)) then
1269 : call endrun('modal_aero_calcsize_diag called for'// &
1270 0 : 'diagnostic list but dgnum_m not associated')
1271 : end if
1272 :
1273 0 : if (.not. present(hygro_m)) then
1274 : call endrun('modal_aero_calcsize_diag called for'// &
1275 0 : 'diagnostic list but hygro_m pointer not present')
1276 : end if
1277 0 : if (.not. associated(hygro_m)) then
1278 : call endrun('modal_aero_calcsize_diag called for'// &
1279 0 : 'diagnostic list but hygro_m not associated')
1280 : end if
1281 :
1282 0 : if (.not. present(dryvol_m)) then
1283 : call endrun('modal_aero_calcsize_diag called for'// &
1284 0 : 'diagnostic list but dryvol_m pointer not present')
1285 : end if
1286 0 : if (.not. associated(dryvol_m)) then
1287 : call endrun('modal_aero_calcsize_diag called for'// &
1288 0 : 'diagnostic list but dryvol_m not associated')
1289 : end if
1290 :
1291 0 : if (.not. present(dryrad_m)) then
1292 : call endrun('modal_aero_calcsize_diag called for'// &
1293 0 : 'diagnostic list but dryrad_m pointer not present')
1294 : end if
1295 0 : if (.not. associated(dryrad_m)) then
1296 : call endrun('modal_aero_calcsize_diag called for'// &
1297 0 : 'diagnostic list but dryrad_m not associated')
1298 : end if
1299 :
1300 0 : if (.not. present(drymass_m)) then
1301 : call endrun('modal_aero_calcsize_diag called for'// &
1302 0 : 'diagnostic list but drymass_m pointer not present')
1303 : end if
1304 0 : if (.not. associated(drymass_m)) then
1305 : call endrun('modal_aero_calcsize_diag called for'// &
1306 0 : 'diagnostic list but drymass_m not associated')
1307 : end if
1308 :
1309 0 : if (.not. present(so4dryvol_m)) then
1310 : call endrun('modal_aero_calcsize_diag called for'// &
1311 0 : 'diagnostic list but so4dryvol_m pointer not present')
1312 : end if
1313 0 : if (.not. associated(so4dryvol_m)) then
1314 : call endrun('modal_aero_calcsize_diag called for'// &
1315 0 : 'diagnostic list but so4dryvol_m not associated')
1316 : end if
1317 :
1318 0 : if (.not. present(naer_m)) then
1319 : call endrun('modal_aero_calcsize_diag called for'// &
1320 0 : 'diagnostic list but naer_m pointer not present')
1321 : end if
1322 0 : if (.not. associated(naer_m)) then
1323 : call endrun('modal_aero_calcsize_diag called for'// &
1324 0 : 'diagnostic list but naer_m not associated')
1325 : end if
1326 :
1327 : end if
1328 :
1329 15480 : do n = 1, nmodes
1330 :
1331 12384 : if (list_idx == 0) then
1332 49536 : call pbuf_get_field(pbuf, dgnum_idx, dgncur_a, start=(/1,1,n/), kount=(/pcols,pver,1/))
1333 : else
1334 0 : dgncur_a => dgnum_m(:,:,n)
1335 : end if
1336 :
1337 : ! get mode properties
1338 : call rad_cnst_get_mode_props(list_idx, n, dgnum=dgnum, dgnumhi=dgnumhi, dgnumlo=dgnumlo, &
1339 12384 : sigmag=sigmag)
1340 :
1341 : ! get mode number mixing ratio
1342 12384 : call rad_cnst_get_mode_num(list_idx, n, 'a', state, pbuf, mode_num)
1343 :
1344 19591488 : dgncur_a(:,:) = dgnum
1345 12384 : dryvol_a(:,:) = 0.0_r8
1346 :
1347 : ! compute dry volume mixrats =
1348 : ! sum_over_components{ component_mass mixrat / density }
1349 12384 : call rad_cnst_get_info(list_idx, n, nspec=nspec)
1350 58824 : do l1 = 1, nspec
1351 :
1352 46440 : call rad_cnst_get_aer_mmr(list_idx, n, l1, 'a', state, pbuf, specmmr)
1353 46440 : call rad_cnst_get_aer_props(list_idx, n, l1, density_aer=specdens)
1354 :
1355 : ! need qmass*dummwdens = (kg/kg-air) * [1/(kg/m3)] = m3/kg-air
1356 46440 : dummwdens = 1.0_r8 / specdens
1357 :
1358 4377744 : do k=top_lev,pver
1359 72162360 : do i=1,ncol
1360 135594000 : dryvol_a(i,k) = dryvol_a(i,k) &
1361 207709920 : + max(0.0_r8, specmmr(i,k))*dummwdens
1362 : end do
1363 : end do
1364 : end do
1365 :
1366 12384 : alnsg = log( sigmag )
1367 12384 : dumfac = exp(4.5_r8*alnsg**2)*pi/6.0_r8
1368 12384 : voltonumblo = 1._r8 / ( (pi/6._r8)*(dgnumlo**3)*exp(4.5_r8*alnsg**2) )
1369 12384 : voltonumbhi = 1._r8 / ( (pi/6._r8)*(dgnumhi**3)*exp(4.5_r8*alnsg**2) )
1370 12384 : v2nxx = voltonumbhi
1371 12384 : v2nyy = voltonumblo
1372 12384 : dgnxx = dgnumhi
1373 12384 : dgnyy = dgnumlo
1374 :
1375 1191960 : do k = top_lev, pver
1376 19243296 : do i = 1, ncol
1377 :
1378 18079200 : drv_a = dryvol_a(i,k)
1379 18079200 : num_a0 = mode_num(i,k)
1380 18079200 : num_a = max( 0.0_r8, num_a0 )
1381 :
1382 19230912 : if (drv_a > 0.0_r8) then
1383 18079200 : if (num_a <= drv_a*v2nxx) then
1384 2716438 : dgncur_a(i,k) = dgnxx
1385 15362762 : else if (num_a >= drv_a*v2nyy) then
1386 1760257 : dgncur_a(i,k) = dgnyy
1387 : else
1388 13602505 : dgncur_a(i,k) = (drv_a/(dumfac*num_a))**third
1389 : end if
1390 : end if
1391 :
1392 : end do
1393 : end do
1394 :
1395 : end do ! nmodes
1396 :
1397 24768 : call modal_aero_calcdry(state, pbuf, list_idx_in, dgnum_m, hygro_m, dryvol_m, dryrad_m, drymass_m, so4dryvol_m, naer_m)
1398 :
1399 3096 : end subroutine modal_aero_calcsize_diag
1400 :
1401 61920 : subroutine modal_aero_calcdry(state, pbuf, list_idx_in, dgnumdry_m, hygro_m, dryvol_m, dryrad_m, drymass_m, so4dryvol_m, naer_m)
1402 :
1403 : type(physics_state), target, intent(in) :: state ! Physics state variables
1404 : type(physics_buffer_desc), pointer :: pbuf(:) ! physics buffer
1405 : integer, optional, intent(in) :: list_idx_in ! diagnostic list index
1406 : real(r8), optional, pointer :: dgnumdry_m(:,:,:)
1407 : real(r8), optional, pointer :: hygro_m(:,:,:)
1408 : real(r8), optional, pointer :: dryvol_m(:,:,:)
1409 : real(r8), optional, pointer :: dryrad_m(:,:,:)
1410 : real(r8), optional, pointer :: drymass_m(:,:,:)
1411 : real(r8), optional, pointer :: so4dryvol_m(:,:,:)
1412 : real(r8), optional, pointer :: naer_m(:,:,:)
1413 :
1414 : real(r8), parameter :: third = 1._r8/3._r8
1415 : real(r8), parameter :: pi43 = pi*4.0_r8/3.0_r8
1416 :
1417 61920 : real(r8), pointer :: maer(:,:) ! aerosol wet mass MR (including water) (kg/kg-air)
1418 61920 : real(r8), pointer :: hygro(:,:,:) ! volume-weighted mean hygroscopicity (--)
1419 61920 : real(r8), pointer :: dryvol(:,:,:) ! single-particle-mean dry volume (m3)
1420 61920 : real(r8), pointer :: dryrad(:,:,:) ! dry volume mean radius of aerosol (m)
1421 61920 : real(r8), pointer :: drymass(:,:,:) ! single-particle-mean dry mass (kg)
1422 61920 : real(r8), pointer :: so4dryvol(:,:,:) ! single-particle-mean so4 dry volume (m3)
1423 61920 : real(r8), pointer :: naer(:,:,:) ! aerosol number MR (bounded!) (#/kg-air)
1424 :
1425 61920 : real(r8), pointer :: dgncur_a(:,:,:)
1426 61920 : real(r8), pointer :: raer(:,:) ! aerosol species MRs (kg/kg and #/kg)
1427 :
1428 : real(r8), pointer :: sulfeq(:,:,:) ! H2SO4 equilibrium mixing ratios over particles (mol/mol)
1429 :
1430 : real(r8) :: dryvolmr(pcols,pver) ! volume MR for aerosol mode (m3/kg)
1431 : real(r8) :: so4dryvolmr(pcols,pver) ! volume MR for sulfate aerosol in mode (m3/kg)
1432 :
1433 : real(r8) :: specdens
1434 : real(r8) :: spechygro, spechygro_1
1435 : real(r8) :: sigmag
1436 : real(r8) :: duma, dumb
1437 : real(r8) :: alnsg
1438 :
1439 : real(r8) :: v2ncur_a
1440 : real(r8) :: drydens ! dry particle density (kg/m^3)
1441 :
1442 : character(len=fieldname_len+3) :: fieldname
1443 : character(len=32) :: spectype
1444 :
1445 : integer :: nmodes, lchnk, ncol, list_idx, i, k, l, m
1446 : integer :: nspec
1447 :
1448 :
1449 61920 : lchnk = state%lchnk
1450 61920 : ncol = state%ncol
1451 :
1452 61920 : list_idx = 0
1453 61920 : if (present(list_idx_in)) then
1454 0 : list_idx = list_idx_in
1455 :
1456 : ! check that all optional args are present
1457 0 : if (.not. present(dgnumdry_m)) then
1458 : call endrun('modal_aero_calcdry called for'// &
1459 0 : 'diagnostic list but required args not present')
1460 : end if
1461 :
1462 : ! arrays for diagnostic calculations must be associated
1463 0 : if (.not. associated(dgnumdry_m)) then
1464 : call endrun('modal_aero_calcdry called for'// &
1465 0 : 'diagnostic list but required args not associated')
1466 : end if
1467 : end if
1468 :
1469 : ! loop over all aerosol modes
1470 61920 : call rad_cnst_get_info(list_idx, nmodes=nmodes)
1471 :
1472 61920 : allocate( maer(pcols,pver))
1473 :
1474 61920 : if (list_idx == 0) then
1475 61920 : call pbuf_get_field(pbuf, dgnum_idx, dgncur_a )
1476 61920 : call pbuf_get_field(pbuf, hygro_idx, hygro)
1477 61920 : call pbuf_get_field(pbuf, dryvol_idx, dryvol)
1478 61920 : call pbuf_get_field(pbuf, dryrad_idx, dryrad)
1479 61920 : call pbuf_get_field(pbuf, drymass_idx, drymass)
1480 61920 : call pbuf_get_field(pbuf, so4dryvol_idx, so4dryvol)
1481 61920 : call pbuf_get_field(pbuf, naer_idx, naer)
1482 : else
1483 0 : dgncur_a => dgnumdry_m
1484 0 : hygro => hygro_m
1485 0 : dryvol => dryvol_m
1486 0 : dryrad => dryrad_m
1487 0 : drymass => drymass_m
1488 0 : so4dryvol => so4dryvol_m
1489 0 : naer => naer_m
1490 : end if
1491 :
1492 391891680 : hygro(:,:,:) = 0._r8
1493 391891680 : so4dryvol(:,:,:) = 0._r8
1494 :
1495 309600 : do m = 1, nmodes
1496 :
1497 391829760 : maer(:,:) = 0._r8
1498 247680 : dryvolmr(:,:) = 0._r8
1499 247680 : so4dryvolmr(:,:) = 0._r8
1500 :
1501 : ! get mode properties
1502 247680 : call rad_cnst_get_mode_props(list_idx, m, sigmag=sigmag)
1503 :
1504 : ! get mode info
1505 247680 : call rad_cnst_get_info(list_idx, m, nspec=nspec)
1506 :
1507 1176480 : do l = 1, nspec
1508 :
1509 : ! get species interstitial mixing ratio ('a')
1510 928800 : call rad_cnst_get_aer_mmr(list_idx, m, l, 'a', state, pbuf, raer)
1511 : call rad_cnst_get_aer_props(list_idx, m, l, density_aer=specdens, &
1512 928800 : hygro_aer=spechygro, spectype=spectype)
1513 :
1514 928800 : if (l == 1) then
1515 : ! save off these values to be used as defaults
1516 247680 : spechygro_1 = spechygro
1517 : end if
1518 :
1519 88483680 : do k = top_lev, pver
1520 1443247200 : do i = 1, ncol
1521 1355940000 : duma = raer(i,k) ! kg/kg air
1522 1355940000 : maer(i,k) = maer(i,k) + duma
1523 1355940000 : dumb = duma/specdens ! m3/kg air
1524 1355940000 : dryvolmr(i,k) = dryvolmr(i,k) + dumb
1525 1355940000 : if (modal_strat_sulfate .and. (trim(spectype).eq.'sulfate')) then
1526 0 : so4dryvolmr(i,k) = so4dryvolmr(i,k) + dumb
1527 : end if
1528 1442318400 : hygro(i,k,m) = hygro(i,k,m) + dumb*spechygro
1529 : end do
1530 : end do
1531 : end do
1532 :
1533 247680 : alnsg = log(sigmag)
1534 :
1535 23281920 : do k = top_lev, pver
1536 384865920 : do i = 1, ncol
1537 :
1538 361584000 : if (dryvolmr(i,k) > 1.0e-30_r8) then
1539 357814963 : hygro(i,k,m) = hygro(i,k,m)/dryvolmr(i,k)
1540 : else
1541 3769037 : hygro(i,k,m) = spechygro_1
1542 : end if
1543 :
1544 : ! dry aerosol properties
1545 :
1546 361584000 : v2ncur_a = 1._r8 / ( (pi/6._r8)*(dgncur_a(i,k,m)**3._r8)*exp(4.5_r8*alnsg**2._r8) )
1547 : ! naer = aerosol number (#/kg)
1548 361584000 : naer(i,k,m) = dryvolmr(i,k)*v2ncur_a
1549 :
1550 : ! compute mean (1 particle) dry volume and mass for each mode
1551 : ! old coding is replaced because the new (1/v2ncur_a) is equal to
1552 : ! the mean particle volume
1553 : ! also moletomass forces maer >= 1.0e-30, so (maer/dryvolmr)
1554 : ! should never cause problems (but check for maer < 1.0e-31 anyway)
1555 361584000 : if (maer(i,k) .gt. 1.0e-31_r8) then
1556 359353540 : drydens = maer(i,k)/dryvolmr(i,k) ! kg/m3 aerosol
1557 : else
1558 : drydens = 1.0_r8
1559 : end if
1560 361584000 : dryvol(i,k,m) = 1.0_r8/v2ncur_a ! m3/particle
1561 361584000 : drymass(i,k,m) = drydens*dryvol(i,k,m) ! kg/particle
1562 384618240 : dryrad(i,k,m) = (dryvol(i,k,m)/pi43)**third ! m
1563 : end do ! i = 1, ncol
1564 : end do ! k = top_lev, pver
1565 :
1566 :
1567 557280 : if (modal_strat_sulfate) then
1568 0 : do k = top_lev, pver
1569 0 : do i = 1, ncol
1570 0 : if (so4dryvolmr(i,k) .gt. 1.0e-31_r8) then
1571 0 : so4dryvol(i,k,m) = dryvol(i,k,m)*so4dryvolmr(i,k)/dryvolmr(i,k)
1572 : else
1573 0 : so4dryvol(i,k,m) = 0.0_r8
1574 : end if
1575 :
1576 : end do ! i = 1, ncol
1577 : end do ! k = top_lev, pver
1578 :
1579 : end if
1580 :
1581 : end do ! m = 1, nmodes
1582 :
1583 61920 : deallocate( maer)
1584 :
1585 :
1586 61920 : end subroutine modal_aero_calcdry
1587 : !----------------------------------------------------------------------
1588 :
1589 : end module modal_aero_calcsize
|
