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 0 : 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 0 : call rad_cnst_get_info(0, nmodes=nmodes)
80 :
81 0 : call pbuf_add_field('DGNUM', 'global', dtype_r8, (/pcols, pver, nmodes/), dgnum_idx)
82 :
83 0 : call pbuf_add_field('HYGRO', 'phys_pkg', dtype_r8, (/pcols,pver,nmodes/), hygro_idx)
84 0 : call pbuf_add_field('DRYVOL', 'phys_pkg', dtype_r8, (/pcols,pver,nmodes/), dryvol_idx)
85 0 : call pbuf_add_field('DRYRAD', 'phys_pkg', dtype_r8, (/pcols,pver,nmodes/), dryrad_idx)
86 0 : call pbuf_add_field('DRYMASS', 'phys_pkg', dtype_r8, (/pcols,pver,nmodes/), drymass_idx)
87 0 : call pbuf_add_field('SO4DRYVOL', 'phys_pkg', dtype_r8, (/pcols,pver,nmodes/), so4dryvol_idx)
88 0 : call pbuf_add_field('NAER', 'phys_pkg', dtype_r8, (/pcols,pver,nmodes/), naer_idx)
89 :
90 0 : end subroutine modal_aero_calcsize_reg
91 :
92 : !===============================================================================
93 : !===============================================================================
94 :
95 0 : subroutine modal_aero_calcsize_init(pbuf2d)
96 0 : 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 0 : call phys_getopts(history_aerosol_out=history_aerosol)
126 :
127 : ! init entities required for both prescribed and prognostic modes
128 :
129 0 : if (is_first_step()) then
130 : ! initialize fields in physics buffer
131 0 : call pbuf_set_field(pbuf2d, dgnum_idx, 0.0_r8)
132 : endif
133 :
134 : #ifndef MODAL_AERO
135 0 : do_adjust_default = .false.
136 0 : do_aitacc_transfer_default = .false.
137 : #else
138 : ! do_adjust_default allows adjustment to be turned on/off
139 : 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 : nait = modeptr_aitken
145 : nacc = modeptr_accum
146 : do_aitacc_transfer_default = .false.
147 : if ((modeptr_aitken > 0) .and. &
148 : (modeptr_accum > 0) .and. &
149 : (modeptr_aitken /= modeptr_accum)) then
150 : do_aitacc_transfer_default = .true.
151 : if (mprognum_amode(nait) <= 0) do_aitacc_transfer_default = .false.
152 : 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 : do n = 1, ntot_amode
159 : if (mprognum_amode(n) <= 0) cycle
160 :
161 : do jac = 1, 2
162 : if (jac == 1) then
163 : tmpnamea = cnst_name(numptr_amode(n))
164 : else
165 : tmpnamea = cnst_name_cw(numptrcw_amode(n))
166 : end if
167 : unit = '#/m2/s'
168 : fieldname = trim(tmpnamea) // '_sfcsiz1'
169 : long_name = trim(tmpnamea) // ' calcsize number-adjust column source'
170 : call addfld( fieldname, horiz_only, 'A', unit, long_name )
171 : if (history_aerosol) then
172 : call add_default(fieldname, 1, ' ')
173 : end if
174 : if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname
175 :
176 : fieldname = trim(tmpnamea) // '_sfcsiz2'
177 : long_name = trim(tmpnamea) // ' calcsize number-adjust column sink'
178 : call addfld( fieldname, horiz_only, 'A', unit, long_name )
179 : if (history_aerosol) then
180 : call add_default(fieldname, 1, ' ')
181 : end if
182 : if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname
183 : end do ! jac = ...
184 : end do ! n = ...
185 :
186 : if ( .not. do_aitacc_transfer_default ) return
187 :
188 : ! check that renaming ipair=1 is aitken-->accum
189 : ipair = 1
190 : if ((modefrm_renamexf(ipair) .ne. nait) .or. &
191 : (modetoo_renamexf(ipair) .ne. nacc)) then
192 : write( 6, '(//2a//)' ) &
193 : '*** modal_aero_calcaersize_init error -- ', &
194 : 'modefrm/too_renamexf(1) are wrong'
195 : call endrun( 'modal_aero_calcaersize_init error' )
196 : end if
197 :
198 : ! define history fields for aitken-accum transfer
199 : do iq = 1, nspecfrm_renamexf(ipair)
200 :
201 : ! jac=1 does interstitial ("_a"); jac=2 does activated ("_c");
202 : 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 : if (jac .eq. 1) then
207 : lsfrm = lspecfrma_renamexf(iq,ipair)
208 : lstoo = lspectooa_renamexf(iq,ipair)
209 : else
210 : lsfrm = lspecfrmc_renamexf(iq,ipair)
211 : lstoo = lspectooc_renamexf(iq,ipair)
212 : end if
213 : if ((lsfrm <= 0) .or. (lstoo <= 0)) cycle
214 :
215 : if (jac .eq. 1) then
216 : tmpnamea = cnst_name(lsfrm)
217 : tmpnameb = cnst_name(lstoo)
218 : else
219 : tmpnamea = cnst_name_cw(lsfrm)
220 : tmpnameb = cnst_name_cw(lstoo)
221 : end if
222 :
223 : unit = 'kg/m2/s'
224 : if ((tmpnamea(1:3) == 'num') .or. &
225 : (tmpnamea(1:3) == 'NUM')) unit = '#/m2/s'
226 : fieldname = trim(tmpnamea) // '_sfcsiz3'
227 : long_name = trim(tmpnamea) // ' calcsize aitken-to-accum adjust column tendency'
228 : call addfld( fieldname, horiz_only, 'A', unit, long_name )
229 : if (history_aerosol) then
230 : call add_default(fieldname, 1, ' ')
231 : end if
232 : if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname
233 :
234 : fieldname = trim(tmpnameb) // '_sfcsiz3'
235 : long_name = trim(tmpnameb) // ' calcsize aitken-to-accum adjust column tendency'
236 : call addfld( fieldname, horiz_only, 'A', unit, long_name )
237 : if (history_aerosol) then
238 : call add_default(fieldname, 1, ' ')
239 : end if
240 : if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname
241 :
242 : fieldname = trim(tmpnamea) // '_sfcsiz4'
243 : long_name = trim(tmpnamea) // ' calcsize accum-to-aitken adjust column tendency'
244 : call addfld( fieldname, horiz_only, 'A', unit, long_name )
245 : if (history_aerosol) then
246 : call add_default(fieldname, 1, ' ')
247 : end if
248 : if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname
249 :
250 : fieldname = trim(tmpnameb) // '_sfcsiz4'
251 : long_name = trim(tmpnameb) // ' calcsize accum-to-aitken adjust column tendency'
252 : call addfld( fieldname, horiz_only, 'A', unit, long_name )
253 : if (history_aerosol) then
254 : call add_default(fieldname, 1, ' ')
255 : end if
256 : if ( masterproc ) write(*,'(2a)') 'calcsize addfld - ', fieldname
257 :
258 : end do ! jac = ...
259 : end do ! iq = ...
260 :
261 : #endif
262 :
263 0 : end subroutine modal_aero_calcsize_init
264 :
265 : !===============================================================================
266 :
267 0 : 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 : real(r8), pointer :: t(:,:) ! Temperature in Kelvin
307 : real(r8), pointer :: pmid(:,:) ! pressure at model levels (Pa)
308 : real(r8), pointer :: pdel(:,:) ! pressure thickness of levels
309 : real(r8), pointer :: q(:,:,:) ! Tracer MR array
310 :
311 : logical, pointer :: dotend(:) ! flag for doing tendency
312 : real(r8), pointer :: dqdt(:,:,:) ! TMR tendency array
313 :
314 : 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 : 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 : 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 : 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 : 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 : 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 : real(r8) :: v2ncur_a(pcols,pver,ntot_amode)
368 : 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 : if (present(do_adjust_in)) then
390 : do_adjust = do_adjust_in
391 : else
392 : do_adjust = do_adjust_default
393 : end if
394 :
395 : if (present(do_aitacc_transfer_in)) then
396 : do_aitacc_transfer = do_aitacc_transfer_in
397 : else
398 : do_aitacc_transfer = do_aitacc_transfer_default
399 : end if
400 :
401 : lchnk = state%lchnk
402 : ncol = state%ncol
403 :
404 : t => state%t
405 : pmid => state%pmid
406 : pdel => state%pdel
407 : q => state%q
408 :
409 : dotend => ptend%lq
410 : dqdt => ptend%q
411 :
412 : call pbuf_get_field(pbuf, dgnum_idx, dgncur_a)
413 :
414 : dotendqqcw(:) = .false.
415 : dqqcwdt(:,:,:) = 0.0_r8
416 : qsrflx(:,:,:,:) = 0.0_r8
417 :
418 : nait = modeptr_aitken
419 : nacc = modeptr_accum
420 :
421 : 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 : tadj = 86400
426 : tadj = max( tadj, deltat )
427 : tadjinv = 1.0_r8/(tadj*(1.0_r8 + 1.0e-15_r8))
428 : fracadj = deltat*tadjinv
429 : 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 : do n = 1, ntot_amode
441 :
442 :
443 : ! initialize all parameters to the default values for the mode
444 : do k=top_lev,pver
445 : do i=1,ncol
446 : ! sgcur_a(i,k,n) = sigmag_amode(n)
447 : ! sgcur_c(i,k,n) = sigmag_amode(n)
448 : dgncur_a(i,k,n) = dgnum_amode(n)
449 : dgncur_c(i,k,n) = dgnum_amode(n)
450 : v2ncur_a(i,k,n) = voltonumb_amode(n)
451 : v2ncur_c(i,k,n) = voltonumb_amode(n)
452 : dryvol_a(i,k) = 0.0_r8
453 : 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 : do l1 = 1, nspec_amode(n)
460 : ! need qmass*dummwdens = (kg/kg-air) * [1/(kg/m3)] = m3/kg-air
461 : dummwdens = 1.0_r8 / specdens_amode(l1,n)
462 : la = lmassptr_amode(l1,n)
463 : do k=top_lev,pver
464 : do i=1,ncol
465 : dryvol_a(i,k) = dryvol_a(i,k) &
466 : + max(0.0_r8,q(i,k,la))*dummwdens
467 : end do
468 : end do
469 :
470 : fldcw => qqcw_get_field(pbuf,lmassptrcw_amode(l1,n),lchnk)
471 : do k=top_lev,pver
472 : do i=1,ncol
473 : dryvol_c(i,k) = dryvol_c(i,k) &
474 : + 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 : lna = numptr_amode(n)
481 : lnc = numptrcw_amode(n)
482 : fldcw => qqcw_get_field(pbuf,numptrcw_amode(n),lchnk,.true.)
483 :
484 :
485 : ! go to section for appropriate number/surface diagnosed/prognosed options
486 : 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 : if (lna > 0) then
493 : dotend(lna) = .true.
494 : do k=top_lev,pver
495 : do i=1,ncol
496 : dqdt(i,k,lna) = (dryvol_a(i,k)*voltonumb_amode(n) &
497 : - q(i,k,lna)) * deltatinv
498 : end do
499 : end do
500 : end if
501 : if (lnc > 0) then
502 : dotendqqcw(lnc) = .true.
503 : do k=top_lev,pver
504 : do i=1,ncol
505 : dqqcwdt(i,k,lnc) = (dryvol_c(i,k)*voltonumb_amode(n) &
506 : - 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 : frelaxadj = 27.0_r8
547 : dumfac = exp(4.5_r8*alnsg_amode(n)**2)*pi/6.0_r8
548 : v2nxx = voltonumbhi_amode(n)
549 : v2nyy = voltonumblo_amode(n)
550 : v2nxxrl = v2nxx/frelaxadj
551 : v2nyyrl = v2nyy*frelaxadj
552 : dgnxx = dgnumhi_amode(n)
553 : dgnyy = dgnumlo_amode(n)
554 : if ( do_aitacc_transfer ) then
555 : if (n == nait) v2nxx = v2nxx/1.0e6_r8
556 : if (n == nacc) v2nyy = v2nyy*1.0e6_r8
557 : v2nxxrl = v2nxx/frelaxadj ! NEW
558 : v2nyyrl = v2nyy*frelaxadj ! NEW
559 : end if
560 :
561 : if (do_adjust) then
562 : dotend(lna) = .true.
563 : dotendqqcw(lnc) = .true.
564 : end if
565 :
566 : do k = top_lev, pver
567 : do i = 1, ncol
568 :
569 : drv_a = dryvol_a(i,k)
570 : num_a0 = q(i,k,lna)
571 : num_a = max( 0.0_r8, num_a0 )
572 : drv_c = dryvol_c(i,k)
573 : num_c0 = fldcw(i,k)
574 : num_c = max( 0.0_r8, num_c0 )
575 :
576 : 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 : 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 : num_a = 0.0_r8
589 : dqdt(i,k,lna) = -num_a0*deltatinv
590 : num_c = 0.0_r8
591 : dqqcwdt(i,k,lnc) = -num_c0*deltatinv
592 : 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 : num_c = 0.0_r8
596 : dqqcwdt(i,k,lnc) = -num_c0*deltatinv
597 : num_a1 = num_a
598 : numbnd = max( drv_a*v2nxx, min( drv_a*v2nyy, num_a1 ) )
599 : num_a = num_a1 + (numbnd - num_a1)*fracadj
600 : dqdt(i,k,lna) = (num_a - num_a0)*deltatinv
601 :
602 : else if (drv_a <= 0.0_r8) then
603 : ! interstitial volume is zero, treat similar to above
604 : num_a = 0.0_r8
605 : dqdt(i,k,lna) = -num_a0*deltatinv
606 : num_c1 = num_c
607 : numbnd = max( drv_c*v2nxx, min( drv_c*v2nyy, num_c1 ) )
608 : num_c = num_c1 + (numbnd - num_c1)*fracadj
609 : 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 : num_a1 = num_a
615 : 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 : numbnd = max( drv_a*v2nxxrl, min( drv_a*v2nyyrl, num_a1 ) )
621 : delnum_a2 = (numbnd - num_a1)*fracadj
622 : num_a2 = num_a1 + delnum_a2
623 : numbnd = max( drv_c*v2nxxrl, min( drv_c*v2nyyrl, num_c1 ) )
624 : delnum_c2 = (numbnd - num_c1)*fracadj
625 : num_c2 = num_c1 + delnum_c2
626 : 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 : num_a1-delnum_c2 ) )
629 : 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 : 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 : drv_t = drv_a + drv_c
636 : num_t2 = num_a2 + num_c2
637 : delnum_a3 = 0.0_r8
638 : delnum_c3 = 0.0_r8
639 : if (num_t2 < drv_t*v2nxx) then
640 : 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 : if ((num_a2 < drv_a*v2nxx) .and. (num_c2 < drv_c*v2nxx)) then
644 : delnum_a3 = delnum_t3*(num_a2/num_t2)
645 : delnum_c3 = delnum_t3*(num_c2/num_t2)
646 : else if (num_c2 < drv_c*v2nxx) then
647 : delnum_c3 = delnum_t3
648 : else if (num_a2 < drv_a*v2nxx) then
649 : delnum_a3 = delnum_t3
650 : end if
651 : else if (num_t2 > drv_t*v2nyy) then
652 : 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 : if ((num_a2 > drv_a*v2nyy) .and. (num_c2 > drv_c*v2nyy)) then
656 : delnum_a3 = delnum_t3*(num_a2/num_t2)
657 : delnum_c3 = delnum_t3*(num_c2/num_t2)
658 : else if (num_c2 > drv_c*v2nyy) then
659 : delnum_c3 = delnum_t3
660 : else if (num_a2 > drv_a*v2nyy) then
661 : delnum_a3 = delnum_t3
662 : end if
663 : end if
664 : num_a = num_a2 + delnum_a3
665 : dqdt(i,k,lna) = (num_a - num_a0)*deltatinv
666 : num_c = num_c2 + delnum_c3
667 : 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 : if (drv_a > 0.0_r8) then
676 : if (num_a <= drv_a*v2nxx) then
677 : dgncur_a(i,k,n) = dgnxx
678 : v2ncur_a(i,k,n) = v2nxx
679 : else if (num_a >= drv_a*v2nyy) then
680 : dgncur_a(i,k,n) = dgnyy
681 : v2ncur_a(i,k,n) = v2nyy
682 : else
683 : dgncur_a(i,k,n) = (drv_a/(dumfac*num_a))**third
684 : v2ncur_a(i,k,n) = num_a/drv_a
685 : end if
686 : end if
687 : pdel_fac = pdel(i,k)/gravit ! = rho*dz
688 : jac = 1
689 : qsrflx(i,lna,1,jac) = qsrflx(i,lna,1,jac) + max(0.0_r8,dqdt(i,k,lna))*pdel_fac
690 : qsrflx(i,lna,2,jac) = qsrflx(i,lna,2,jac) + min(0.0_r8,dqdt(i,k,lna))*pdel_fac
691 :
692 : if (drv_c > 0.0_r8) then
693 : if (num_c <= drv_c*v2nxx) then
694 : dgncur_c(i,k,n) = dgnumhi_amode(n)
695 : v2ncur_c(i,k,n) = v2nxx
696 : else if (num_c >= drv_c*v2nyy) then
697 : dgncur_c(i,k,n) = dgnumlo_amode(n)
698 : v2ncur_c(i,k,n) = v2nyy
699 : else
700 : dgncur_c(i,k,n) = (drv_c/(dumfac*num_c))**third
701 : v2ncur_c(i,k,n) = num_c/drv_c
702 : end if
703 : end if
704 : jac = 2
705 : qsrflx(i,lnc,1,jac) = qsrflx(i,lnc,1,jac) + max(0.0_r8,dqqcwdt(i,k,lnc))*pdel_fac
706 : 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 : if ( do_aitacc_transfer ) then
711 : if (n == nait) then
712 : drv_a_aitsv(i,k) = drv_a
713 : num_a_aitsv(i,k) = num_a
714 : drv_c_aitsv(i,k) = drv_c
715 : num_c_aitsv(i,k) = num_c
716 : else if (n == nacc) then
717 : drv_a_accsv(i,k) = drv_a
718 : num_a_accsv(i,k) = num_a
719 : drv_c_accsv(i,k) = drv_c
720 : num_c_accsv(i,k) = num_c
721 : end if
722 : end if
723 : drv_a_sv(i,k,n) = drv_a
724 : num_a_sv(i,k,n) = num_a
725 : drv_c_sv(i,k,n) = drv_c
726 : 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 : ixfer_ait2acc_sv(:,:) = 0
753 : ixfer_acc2ait_sv(:,:) = 0
754 : 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 : if (npair_renamexf .le. 0) then
760 : npair_renamexf = 0
761 : ! call modal_aero_rename_init
762 : if (npair_renamexf .le. 0) then
763 : write( 6, '(//a//)' ) &
764 : '*** modal_aero_calcaersize_sub error -- npair_renamexf <= 0'
765 : call endrun( 'modal_aero_calcaersize_sub error' )
766 : end if
767 : end if
768 :
769 : ! check that renaming ipair=1 is aitken-->accum
770 : ipair = 1
771 : if ((modefrm_renamexf(ipair) .ne. nait) .or. &
772 : (modetoo_renamexf(ipair) .ne. nacc)) then
773 : write( 6, '(//2a//)' ) &
774 : '*** modal_aero_calcaersize_sub error -- ', &
775 : 'modefrm/too_renamexf(1) are wrong'
776 : call endrun( 'modal_aero_calcaersize_sub error' )
777 : end if
778 :
779 : ! set dotend() for species that will be transferred
780 : do iq = 1, nspecfrm_renamexf(ipair)
781 : lsfrm = lspecfrma_renamexf(iq,ipair)
782 : lstoo = lspectooa_renamexf(iq,ipair)
783 : if ((lsfrm > 0) .and. (lstoo > 0)) then
784 : dotend(lsfrm) = .true.
785 : dotend(lstoo) = .true.
786 : end if
787 : lsfrm = lspecfrmc_renamexf(iq,ipair)
788 : lstoo = lspectooc_renamexf(iq,ipair)
789 : if ((lsfrm > 0) .and. (lstoo > 0)) then
790 : dotendqqcw(lsfrm) = .true.
791 : dotendqqcw(lstoo) = .true.
792 : end if
793 : end do
794 :
795 : ! identify accum species cannot be transferred to aitken mode
796 : noxf_acc2ait(:) = .true.
797 : do l1 = 1, nspec_amode(nacc)
798 : la = lmassptr_amode(l1,nacc)
799 : do iq = 1, nspecfrm_renamexf(ipair)
800 : if (lspectooa_renamexf(iq,ipair) == la) then
801 : 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 : v2nzz = sqrt(voltonumb_amode(nait)*voltonumb_amode(nacc))
809 :
810 : ! loop over columns and levels
811 : do k = top_lev, pver
812 : do i = 1, ncol
813 :
814 : pdel_fac = pdel(i,k)/gravit ! = rho*dz
815 : xfertend_num(:,:) = 0.0_r8
816 :
817 : ! compute aitken --> accum transfer rates
818 : ixfer_ait2acc = 0
819 : xfercoef_num_ait2acc = 0.0_r8
820 : xfercoef_vol_ait2acc = 0.0_r8
821 :
822 : drv_t = drv_a_aitsv(i,k) + drv_c_aitsv(i,k)
823 : num_t = num_a_aitsv(i,k) + num_c_aitsv(i,k)
824 : if (drv_t > 0.0_r8) then
825 : if (num_t < drv_t*v2nzz) then
826 : ixfer_ait2acc = 1
827 : 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 : (voltonumb_amode(nacc) - v2nzz)
833 : xferfrac_num_ait2acc = xferfrac_vol_ait2acc* &
834 : (drv_t*voltonumb_amode(nacc)/num_t)
835 : 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 : else if ((xferfrac_num_ait2acc >= 1.0_r8) .or. &
840 : (xferfrac_vol_ait2acc >= 1.0_r8)) then
841 : xferfrac_num_ait2acc = 1.0_r8
842 : xferfrac_vol_ait2acc = 1.0_r8
843 : end if
844 : end if
845 : xfercoef_num_ait2acc = xferfrac_num_ait2acc*tadjinv
846 : xfercoef_vol_ait2acc = xferfrac_vol_ait2acc*tadjinv
847 : xfertend_num(1,1) = num_a_aitsv(i,k)*xfercoef_num_ait2acc
848 : 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 : ixfer_acc2ait = 0
859 : xfercoef_num_acc2ait = 0.0_r8
860 : xfercoef_vol_acc2ait = 0.0_r8
861 :
862 : drv_t = drv_a_accsv(i,k) + drv_c_accsv(i,k)
863 : num_t = num_a_accsv(i,k) + num_c_accsv(i,k)
864 : drv_a_noxf = 0.0_r8
865 : drv_c_noxf = 0.0_r8
866 : if (drv_t > 0.0_r8) then
867 : if (num_t > drv_t*v2nzz) then
868 : do l1 = 1, nspec_amode(nacc)
869 :
870 : if ( noxf_acc2ait(l1) ) then
871 : ! need qmass*dummwdens = (kg/kg-air) * [1/(kg/m3)] = m3/kg-air
872 : dummwdens = 1.0_r8 / specdens_amode(l1,nacc)
873 : la = lmassptr_amode(l1,nacc)
874 : drv_a_noxf = drv_a_noxf &
875 : + max(0.0_r8,q(i,k,la))*dummwdens
876 : lc = lmassptrcw_amode(l1,nacc)
877 :
878 : fldcw => qqcw_get_field(pbuf,lmassptrcw_amode(l1,nacc),lchnk)
879 : drv_c_noxf = drv_c_noxf &
880 : + max(0.0_r8,fldcw(i,k))*dummwdens
881 : end if
882 : end do
883 : drv_t_noxf = drv_a_noxf + drv_c_noxf
884 : num_t_noxf = drv_t_noxf*voltonumblo_amode(nacc)
885 : num_t0 = num_t
886 : drv_t0 = drv_t
887 : num_t = max( 0.0_r8, num_t - num_t_noxf )
888 : drv_t = max( 0.0_r8, drv_t - drv_t_noxf )
889 : end if
890 : end if
891 :
892 : if (drv_t > 0.0_r8) then
893 : if (num_t > drv_t*v2nzz) then
894 : ixfer_acc2ait = 1
895 : 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 : (voltonumb_amode(nait) - v2nzz)
901 : xferfrac_num_acc2ait = xferfrac_vol_acc2ait* &
902 : (drv_t*voltonumb_amode(nait)/num_t)
903 : 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 : else if ((xferfrac_num_acc2ait >= 1.0_r8) .or. &
908 : (xferfrac_vol_acc2ait >= 1.0_r8)) then
909 : xferfrac_num_acc2ait = 1.0_r8
910 : xferfrac_vol_acc2ait = 1.0_r8
911 : end if
912 : end if
913 : duma = 1.0e-37_r8
914 : xferfrac_num_acc2ait = xferfrac_num_acc2ait* &
915 : num_t/max( duma, num_t0 )
916 : xfercoef_num_acc2ait = xferfrac_num_acc2ait*tadjinv
917 : xfercoef_vol_acc2ait = xferfrac_vol_acc2ait*tadjinv
918 : xfertend_num(2,1) = num_a_accsv(i,k)*xfercoef_num_acc2ait
919 : 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 : if (ixfer_ait2acc+ixfer_acc2ait > 0) then
925 : ixfer_ait2acc_sv(i,k) = ixfer_ait2acc
926 : ixfer_acc2ait_sv(i,k) = ixfer_acc2ait
927 :
928 : !
929 : ! compute new dgncur & v2ncur for aitken & accum modes
930 : !
931 : ! currently inactive
932 : do n = nait, nacc, (nacc-nait)
933 : if (n .eq. nait) then
934 : duma = (xfertend_num(1,1) - xfertend_num(2,1))*deltat
935 : num_a = max( 0.0_r8, num_a_aitsv(i,k) - duma )
936 : num_a_acc = max( 0.0_r8, num_a_accsv(i,k) + duma )
937 : duma = (drv_a_aitsv(i,k)*xfercoef_vol_ait2acc - &
938 : (drv_a_accsv(i,k)-drv_a_noxf)*xfercoef_vol_acc2ait)*deltat
939 : drv_a = max( 0.0_r8, drv_a_aitsv(i,k) - duma )
940 : drv_a_acc = max( 0.0_r8, drv_a_accsv(i,k) + duma )
941 : duma = (xfertend_num(1,2) - xfertend_num(2,2))*deltat
942 : num_c = max( 0.0_r8, num_c_aitsv(i,k) - duma )
943 : num_c_acc = max( 0.0_r8, num_c_accsv(i,k) + duma )
944 : duma = (drv_c_aitsv(i,k)*xfercoef_vol_ait2acc - &
945 : (drv_c_accsv(i,k)-drv_c_noxf)*xfercoef_vol_acc2ait)*deltat
946 : drv_c = max( 0.0_r8, drv_c_aitsv(i,k) - duma )
947 : 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 : if (drv_a > 0.0_r8) then
956 : if (num_a <= drv_a*voltonumbhi_amode(n)) then
957 : dgncur_a(i,k,n) = dgnumhi_amode(n)
958 : v2ncur_a(i,k,n) = voltonumbhi_amode(n)
959 : else if (num_a >= drv_a*voltonumblo_amode(n)) then
960 : dgncur_a(i,k,n) = dgnumlo_amode(n)
961 : v2ncur_a(i,k,n) = voltonumblo_amode(n)
962 : else
963 : dgncur_a(i,k,n) = (drv_a/(dumfac*num_a))**third
964 : v2ncur_a(i,k,n) = num_a/drv_a
965 : end if
966 : else
967 : dgncur_a(i,k,n) = dgnum_amode(n)
968 : v2ncur_a(i,k,n) = voltonumb_amode(n)
969 : end if
970 :
971 : if (drv_c > 0.0_r8) then
972 : if (num_c <= drv_c*voltonumbhi_amode(n)) then
973 : dgncur_c(i,k,n) = dgnumhi_amode(n)
974 : v2ncur_c(i,k,n) = voltonumbhi_amode(n)
975 : else if (num_c >= drv_c*voltonumblo_amode(n)) then
976 : dgncur_c(i,k,n) = dgnumlo_amode(n)
977 : v2ncur_c(i,k,n) = voltonumblo_amode(n)
978 : else
979 : dgncur_c(i,k,n) = (drv_c/(dumfac*num_c))**third
980 : v2ncur_c(i,k,n) = num_c/drv_c
981 : end if
982 : else
983 : dgncur_c(i,k,n) = dgnum_amode(n)
984 : 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 : if ( masterproc ) then
995 : if (idiagaa > 0) then
996 : do j = 1, 2
997 : do iq = 1, nspecfrm_renamexf(ipair)
998 : do jac = 1, 2
999 : if (j .eq. 1) then
1000 : if (jac .eq. 1) then
1001 : lsfrm = lspecfrma_renamexf(iq,ipair)
1002 : lstoo = lspectooa_renamexf(iq,ipair)
1003 : else
1004 : lsfrm = lspecfrmc_renamexf(iq,ipair)
1005 : lstoo = lspectooc_renamexf(iq,ipair)
1006 : end if
1007 : else
1008 : if (jac .eq. 1) then
1009 : lsfrm = lspectooa_renamexf(iq,ipair)
1010 : lstoo = lspecfrma_renamexf(iq,ipair)
1011 : else
1012 : lsfrm = lspectooc_renamexf(iq,ipair)
1013 : lstoo = lspecfrmc_renamexf(iq,ipair)
1014 : end if
1015 : end if
1016 : write( 6, '(a,3i3,2i4)' ) 'calcsize j,iq,jac, lsfrm,lstoo', &
1017 : j,iq,jac, lsfrm,lstoo
1018 : end do
1019 : end do
1020 : end do
1021 : end if
1022 : end if
1023 : idiagaa = -1
1024 :
1025 :
1026 : ! j=1 does aitken-->accum; j=2 does accum-->aitken
1027 : do j = 1, 2
1028 :
1029 : if ((j .eq. 1 .and. ixfer_ait2acc > 0) .or. &
1030 : (j .eq. 2 .and. ixfer_acc2ait > 0)) then
1031 :
1032 : jsrflx = j+2
1033 : if (j .eq. 1) then
1034 : xfercoef = xfercoef_vol_ait2acc
1035 : else
1036 : xfercoef = xfercoef_vol_acc2ait
1037 : end if
1038 :
1039 : do iq = 1, nspecfrm_renamexf(ipair)
1040 :
1041 : ! jac=1 does interstitial ("_a"); jac=2 does activated ("_c");
1042 : 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 : if (j .eq. 1) then
1049 : if (jac .eq. 1) then
1050 : lsfrm = lspecfrma_renamexf(iq,ipair)
1051 : lstoo = lspectooa_renamexf(iq,ipair)
1052 : else
1053 : lsfrm = lspecfrmc_renamexf(iq,ipair)
1054 : lstoo = lspectooc_renamexf(iq,ipair)
1055 : end if
1056 : else
1057 : if (jac .eq. 1) then
1058 : lsfrm = lspectooa_renamexf(iq,ipair)
1059 : lstoo = lspecfrma_renamexf(iq,ipair)
1060 : else
1061 : lsfrm = lspectooc_renamexf(iq,ipair)
1062 : lstoo = lspecfrmc_renamexf(iq,ipair)
1063 : end if
1064 : end if
1065 :
1066 : if ((lsfrm > 0) .and. (lstoo > 0)) then
1067 : if (jac .eq. 1) then
1068 : if (iq .eq. 1) then
1069 : xfertend = xfertend_num(j,jac)
1070 : else
1071 : xfertend = max(0.0_r8,q(i,k,lsfrm))*xfercoef
1072 : end if
1073 : dqdt(i,k,lsfrm) = dqdt(i,k,lsfrm) - xfertend
1074 : dqdt(i,k,lstoo) = dqdt(i,k,lstoo) + xfertend
1075 : else
1076 : if (iq .eq. 1) then
1077 : xfertend = xfertend_num(j,jac)
1078 : else
1079 : fldcw => qqcw_get_field(pbuf,lsfrm,lchnk)
1080 : xfertend = max(0.0_r8,fldcw(i,k))*xfercoef
1081 : end if
1082 : dqqcwdt(i,k,lsfrm) = dqqcwdt(i,k,lsfrm) - xfertend
1083 : dqqcwdt(i,k,lstoo) = dqqcwdt(i,k,lstoo) + xfertend
1084 : end if
1085 : qsrflx(i,lsfrm,jsrflx,jac) = qsrflx(i,lsfrm,jsrflx,jac) - xfertend*pdel_fac
1086 : 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 : lsfrm = -123456789 ! executable statement for debugging
1101 :
1102 :
1103 : !
1104 : ! apply tendencies to cloud-borne species MRs
1105 : !
1106 : do l = 1, pcnst
1107 : lc = l
1108 : if ( lc>0 .and. dotendqqcw(lc) ) then
1109 : fldcw=> qqcw_get_field(pbuf,l,lchnk)
1110 : do k = top_lev, pver
1111 : do i = 1, ncol
1112 : fldcw(i,k) = max( 0.0_r8, &
1113 : (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 : if ( .not. do_adjust ) return
1125 :
1126 : do n = 1, ntot_amode
1127 : if (mprognum_amode(n) <= 0) cycle
1128 :
1129 : do jac = 1, 2
1130 : if (jac == 1) then
1131 : l = numptr_amode(n)
1132 : tmpnamea = cnst_name(l)
1133 : else
1134 : l = numptrcw_amode(n)
1135 : tmpnamea = cnst_name_cw(l)
1136 : end if
1137 : fieldname = trim(tmpnamea) // '_sfcsiz1'
1138 : call outfld( fieldname, qsrflx(:,l,1,jac), pcols, lchnk)
1139 :
1140 : fieldname = trim(tmpnamea) // '_sfcsiz2'
1141 : 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 : if ( .not. do_aitacc_transfer ) return
1149 :
1150 : do iq = 1, nspecfrm_renamexf(ipair)
1151 :
1152 : ! jac=1 does interstitial ("_a"); jac=2 does activated ("_c");
1153 : 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 : if (jac .eq. 1) then
1158 : lsfrm = lspecfrma_renamexf(iq,ipair)
1159 : lstoo = lspectooa_renamexf(iq,ipair)
1160 : else
1161 : lsfrm = lspecfrmc_renamexf(iq,ipair)
1162 : lstoo = lspectooc_renamexf(iq,ipair)
1163 : end if
1164 : if ((lsfrm <= 0) .or. (lstoo <= 0)) cycle
1165 :
1166 : if (jac .eq. 1) then
1167 : tmpnamea = cnst_name(lsfrm)
1168 : tmpnameb = cnst_name(lstoo)
1169 : else
1170 : tmpnamea = cnst_name_cw(lsfrm)
1171 : tmpnameb = cnst_name_cw(lstoo)
1172 : end if
1173 : if ((lsfrm <= 0) .or. (lstoo <= 0)) cycle
1174 :
1175 : fieldname = trim(tmpnamea) // '_sfcsiz3'
1176 : call outfld( fieldname, qsrflx(:,lsfrm,3,jac), pcols, lchnk)
1177 :
1178 : fieldname = trim(tmpnameb) // '_sfcsiz3'
1179 : call outfld( fieldname, qsrflx(:,lstoo,3,jac), pcols, lchnk)
1180 :
1181 : fieldname = trim(tmpnamea) // '_sfcsiz4'
1182 : call outfld( fieldname, qsrflx(:,lsfrm,4,jac), pcols, lchnk)
1183 :
1184 : fieldname = trim(tmpnameb) // '_sfcsiz4'
1185 : call outfld( fieldname, qsrflx(:,lstoo,4,jac), pcols, lchnk)
1186 :
1187 : end do ! jac = ...
1188 : end do ! iq = ...
1189 :
1190 : call modal_aero_calcdry(state, pbuf)
1191 :
1192 : #endif
1193 :
1194 0 : end subroutine modal_aero_calcsize_sub
1195 :
1196 :
1197 : !----------------------------------------------------------------------
1198 :
1199 :
1200 0 : 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 0 : real(r8), pointer :: dgncur_a(:,:) ! (pcols,pver)
1234 :
1235 :
1236 : real(r8), parameter :: third = 1.0_r8/3.0_r8
1237 :
1238 0 : real(r8), pointer :: mode_num(:,:) ! mode number mixing ratio
1239 0 : 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 0 : lchnk = state%lchnk
1256 0 : ncol = state%ncol
1257 :
1258 0 : list_idx = 0 ! climate list by default
1259 0 : if (present(list_idx_in)) list_idx = list_idx_in
1260 :
1261 0 : call rad_cnst_get_info(list_idx, nmodes=nmodes)
1262 :
1263 0 : 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 0 : do n = 1, nmodes
1330 :
1331 0 : if (list_idx == 0) then
1332 0 : 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 0 : sigmag=sigmag)
1340 :
1341 : ! get mode number mixing ratio
1342 0 : call rad_cnst_get_mode_num(list_idx, n, 'a', state, pbuf, mode_num)
1343 :
1344 0 : dgncur_a(:,:) = dgnum
1345 0 : dryvol_a(:,:) = 0.0_r8
1346 :
1347 : ! compute dry volume mixrats =
1348 : ! sum_over_components{ component_mass mixrat / density }
1349 0 : call rad_cnst_get_info(list_idx, n, nspec=nspec)
1350 0 : do l1 = 1, nspec
1351 :
1352 0 : call rad_cnst_get_aer_mmr(list_idx, n, l1, 'a', state, pbuf, specmmr)
1353 0 : 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 0 : dummwdens = 1.0_r8 / specdens
1357 :
1358 0 : do k=top_lev,pver
1359 0 : do i=1,ncol
1360 0 : dryvol_a(i,k) = dryvol_a(i,k) &
1361 0 : + max(0.0_r8, specmmr(i,k))*dummwdens
1362 : end do
1363 : end do
1364 : end do
1365 :
1366 0 : alnsg = log( sigmag )
1367 0 : dumfac = exp(4.5_r8*alnsg**2)*pi/6.0_r8
1368 0 : voltonumblo = 1._r8 / ( (pi/6._r8)*(dgnumlo**3)*exp(4.5_r8*alnsg**2) )
1369 0 : voltonumbhi = 1._r8 / ( (pi/6._r8)*(dgnumhi**3)*exp(4.5_r8*alnsg**2) )
1370 0 : v2nxx = voltonumbhi
1371 0 : v2nyy = voltonumblo
1372 0 : dgnxx = dgnumhi
1373 0 : dgnyy = dgnumlo
1374 :
1375 0 : do k = top_lev, pver
1376 0 : do i = 1, ncol
1377 :
1378 0 : drv_a = dryvol_a(i,k)
1379 0 : num_a0 = mode_num(i,k)
1380 0 : num_a = max( 0.0_r8, num_a0 )
1381 :
1382 0 : if (drv_a > 0.0_r8) then
1383 0 : if (num_a <= drv_a*v2nxx) then
1384 0 : dgncur_a(i,k) = dgnxx
1385 0 : else if (num_a >= drv_a*v2nyy) then
1386 0 : dgncur_a(i,k) = dgnyy
1387 : else
1388 0 : 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 0 : 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 0 : end subroutine modal_aero_calcsize_diag
1400 :
1401 0 : 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 0 : real(r8), pointer :: maer(:,:) ! aerosol wet mass MR (including water) (kg/kg-air)
1418 0 : real(r8), pointer :: hygro(:,:,:) ! volume-weighted mean hygroscopicity (--)
1419 0 : real(r8), pointer :: dryvol(:,:,:) ! single-particle-mean dry volume (m3)
1420 0 : real(r8), pointer :: dryrad(:,:,:) ! dry volume mean radius of aerosol (m)
1421 0 : real(r8), pointer :: drymass(:,:,:) ! single-particle-mean dry mass (kg)
1422 0 : real(r8), pointer :: so4dryvol(:,:,:) ! single-particle-mean so4 dry volume (m3)
1423 0 : real(r8), pointer :: naer(:,:,:) ! aerosol number MR (bounded!) (#/kg-air)
1424 :
1425 0 : real(r8), pointer :: dgncur_a(:,:,:)
1426 0 : 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 0 : lchnk = state%lchnk
1450 0 : ncol = state%ncol
1451 :
1452 0 : list_idx = 0
1453 0 : 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 0 : call rad_cnst_get_info(list_idx, nmodes=nmodes)
1471 :
1472 0 : allocate( maer(pcols,pver))
1473 :
1474 0 : if (list_idx == 0) then
1475 0 : call pbuf_get_field(pbuf, dgnum_idx, dgncur_a )
1476 0 : call pbuf_get_field(pbuf, hygro_idx, hygro)
1477 0 : call pbuf_get_field(pbuf, dryvol_idx, dryvol)
1478 0 : call pbuf_get_field(pbuf, dryrad_idx, dryrad)
1479 0 : call pbuf_get_field(pbuf, drymass_idx, drymass)
1480 0 : call pbuf_get_field(pbuf, so4dryvol_idx, so4dryvol)
1481 0 : 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 0 : hygro(:,:,:) = 0._r8
1493 0 : so4dryvol(:,:,:) = 0._r8
1494 :
1495 0 : do m = 1, nmodes
1496 :
1497 0 : maer(:,:) = 0._r8
1498 0 : dryvolmr(:,:) = 0._r8
1499 0 : so4dryvolmr(:,:) = 0._r8
1500 :
1501 : ! get mode properties
1502 0 : call rad_cnst_get_mode_props(list_idx, m, sigmag=sigmag)
1503 :
1504 : ! get mode info
1505 0 : call rad_cnst_get_info(list_idx, m, nspec=nspec)
1506 :
1507 0 : do l = 1, nspec
1508 :
1509 : ! get species interstitial mixing ratio ('a')
1510 0 : 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 0 : hygro_aer=spechygro, spectype=spectype)
1513 :
1514 0 : if (l == 1) then
1515 : ! save off these values to be used as defaults
1516 0 : spechygro_1 = spechygro
1517 : end if
1518 :
1519 0 : do k = top_lev, pver
1520 0 : do i = 1, ncol
1521 0 : duma = raer(i,k) ! kg/kg air
1522 0 : maer(i,k) = maer(i,k) + duma
1523 0 : dumb = duma/specdens ! m3/kg air
1524 0 : dryvolmr(i,k) = dryvolmr(i,k) + dumb
1525 0 : if (modal_strat_sulfate .and. (trim(spectype).eq.'sulfate')) then
1526 0 : so4dryvolmr(i,k) = so4dryvolmr(i,k) + dumb
1527 : end if
1528 0 : hygro(i,k,m) = hygro(i,k,m) + dumb*spechygro
1529 : end do
1530 : end do
1531 : end do
1532 :
1533 0 : alnsg = log(sigmag)
1534 :
1535 0 : do k = top_lev, pver
1536 0 : do i = 1, ncol
1537 :
1538 0 : if (dryvolmr(i,k) > 1.0e-30_r8) then
1539 0 : hygro(i,k,m) = hygro(i,k,m)/dryvolmr(i,k)
1540 : else
1541 0 : hygro(i,k,m) = spechygro_1
1542 : end if
1543 :
1544 : ! dry aerosol properties
1545 :
1546 0 : 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 0 : 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 0 : if (maer(i,k) .gt. 1.0e-31_r8) then
1556 0 : drydens = maer(i,k)/dryvolmr(i,k) ! kg/m3 aerosol
1557 : else
1558 : drydens = 1.0_r8
1559 : end if
1560 0 : dryvol(i,k,m) = 1.0_r8/v2ncur_a ! m3/particle
1561 0 : drymass(i,k,m) = drydens*dryvol(i,k,m) ! kg/particle
1562 0 : 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 0 : 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 0 : deallocate( maer)
1584 :
1585 :
1586 0 : end subroutine modal_aero_calcdry
1587 : !----------------------------------------------------------------------
1588 :
1589 : end module modal_aero_calcsize
|
