Line data Source code
1 :
2 : module modal_aero_convproc
3 : !---------------------------------------------------------------------------------
4 : ! Purpose:
5 : !
6 : ! CAM interface to aerosol/trace-gas convective cloud processing scheme
7 : !
8 : ! currently these routines assume stratiform and convective clouds only interact
9 : ! through the detrainment of convective cloudborne material into stratiform clouds
10 : !
11 : ! thus the stratiform-cloudborne aerosols (in the qqcw array) are not processed
12 : ! by the convective up/downdrafts, but are affected by the detrainment
13 : !
14 : ! Author: R. C. Easter
15 : !
16 : !---------------------------------------------------------------------------------
17 :
18 : use shr_kind_mod, only: r8=>shr_kind_r8
19 : use spmd_utils, only: masterproc
20 : use physconst, only: gravit, rair
21 : use ppgrid, only: pver, pcols, pverp
22 : use constituents, only: pcnst, cnst_name
23 : use constituents, only: cnst_species_class, cnst_spec_class_aerosol, cnst_spec_class_gas
24 : use phys_control, only: phys_getopts
25 :
26 : use physics_types, only: physics_state, physics_ptend
27 : use physics_buffer, only: physics_buffer_desc, pbuf_get_index, pbuf_get_field
28 :
29 : use time_manager, only: get_nstep
30 : use cam_history, only: outfld, addfld, add_default, horiz_only
31 : use cam_logfile, only: iulog
32 : use cam_abortutils, only: endrun
33 :
34 : use modal_aero_data, only: lmassptr_amode, nspec_amode, ntot_amode, numptr_amode
35 : use modal_aero_data, only: lptr_so4_a_amode, lptr_dust_a_amode, lptr_nacl_a_amode, mode_size_order
36 : use modal_aero_data, only: lptr2_pom_a_amode, lptr2_soa_a_amode, lptr2_bc_a_amode, nsoa, npoa, nbc
37 : use modal_aero_data, only: lptr_msa_a_amode, lptr_nh4_a_amode, lptr_no3_a_amode
38 :
39 : use modal_aerosol_properties_mod, only: modal_aerosol_properties
40 :
41 : implicit none
42 : private
43 : save
44 :
45 : public :: &
46 : ma_convproc_register, &
47 : ma_convproc_init, &
48 : ma_convproc_intr, &
49 : ma_convproc_readnl
50 :
51 : ! namelist options
52 : ! NOTE: These are the defaults for CAM6.
53 : logical, protected, public :: convproc_do_gas = .false.
54 : logical, protected, public :: deepconv_wetdep_history = .true.
55 : logical, protected, public :: convproc_do_deep = .true.
56 : ! NOTE: Shallow convection processing does not currently work with CLUBB.
57 : logical, protected, public :: convproc_do_shallow = .false.
58 : ! NOTE: These are the defaults for the Eaton/Wang parameterization.
59 : logical, protected, public :: convproc_do_evaprain_atonce = .false.
60 : real(r8), protected, public :: convproc_pom_spechygro = -1._r8
61 : real(r8), protected, public :: convproc_wup_max = 4.0_r8
62 :
63 : logical, parameter :: use_cwaer_for_activate_maxsat = .false.
64 : logical, parameter :: apply_convproc_tend_to_ptend = .true.
65 :
66 : real(r8) :: hund_ovr_g ! = 100.0_r8/gravit
67 : ! used with zm_conv mass fluxes and delta-p
68 : ! for mu = [mbar/s], mu*hund_ovr_g = [kg/m2/s]
69 : ! for dp = [mbar] and q = [kg/kg], q*dp*hund_ovr_g = [kg/m2]
70 :
71 : ! method1_activate_nlayers = number of layers (including cloud base) where activation is applied
72 : integer, parameter :: method1_activate_nlayers = 2
73 : ! method2_activate_smaxmax = the uniform or peak supersat value (as 0-1 fraction = percent*0.01)
74 : real(r8), parameter :: method2_activate_smaxmax = 0.003_r8
75 :
76 : ! method_reduce_actfrac = 1 -- multiply activation fractions by factor_reduce_actfrac
77 : ! (this works ok with convproc_method_activate = 1 but not for ... = 2)
78 : ! = 2 -- do 2 iterations to get an overall reduction by factor_reduce_actfrac
79 : ! (this works ok with convproc_method_activate = 1 or 2)
80 : ! = other -- do nothing involving reduce_actfrac
81 : integer, parameter :: method_reduce_actfrac = 0
82 : real(r8), parameter :: factor_reduce_actfrac = 0.5_r8
83 :
84 : ! convproc_method_activate - 1=apply abdulrazzak-ghan to entrained aerosols for lowest nlayers
85 : ! 2=do secondary activation with prescribed supersat
86 : integer, parameter :: convproc_method_activate = 2
87 :
88 : logical :: convproc_do_aer
89 :
90 : ! physics buffer indices
91 : integer :: fracis_idx = 0
92 :
93 : integer :: rprddp_idx = 0
94 : integer :: rprdsh_idx = 0
95 : integer :: nevapr_shcu_idx = 0
96 : integer :: nevapr_dpcu_idx = 0
97 :
98 : integer :: icwmrdp_idx = 0
99 : integer :: icwmrsh_idx = 0
100 : integer :: sh_frac_idx = 0
101 : integer :: dp_frac_idx = 0
102 :
103 : integer :: zm_mu_idx = 0
104 : integer :: zm_eu_idx = 0
105 : integer :: zm_du_idx = 0
106 : integer :: zm_md_idx = 0
107 : integer :: zm_ed_idx = 0
108 : integer :: zm_dp_idx = 0
109 : integer :: zm_dsubcld_idx = 0
110 : integer :: zm_jt_idx = 0
111 : integer :: zm_maxg_idx = 0
112 : integer :: zm_ideep_idx = 0
113 :
114 : integer :: cmfmc_sh_idx = 0
115 : integer :: sh_e_ed_ratio_idx = 0
116 :
117 : integer :: istat
118 :
119 : logical, parameter :: debug=.false.
120 :
121 : type(modal_aerosol_properties), pointer :: aero_props_obj => null()
122 :
123 : !=========================================================================================
124 : contains
125 : !=========================================================================================
126 :
127 167511 : subroutine ma_convproc_register
128 :
129 0 : end subroutine ma_convproc_register
130 :
131 : !=========================================================================================
132 1536 : subroutine ma_convproc_readnl(nlfile)
133 :
134 : use namelist_utils, only: find_group_name
135 : use spmd_utils, only: mpicom, masterprocid, mpi_real8, mpi_logical
136 :
137 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
138 :
139 : ! Local variables
140 : integer :: unitn, ierr
141 : character(len=*), parameter :: subname = 'ma_convproc_readnl'
142 :
143 : namelist /aerosol_convproc_opts/ convproc_do_gas, deepconv_wetdep_history, convproc_do_deep, &
144 : convproc_do_shallow, convproc_do_evaprain_atonce, convproc_pom_spechygro, convproc_wup_max
145 :
146 : ! Read namelist
147 1536 : if (masterproc) then
148 2 : open( newunit=unitn, file=trim(nlfile), status='old' )
149 2 : call find_group_name(unitn, 'aerosol_convproc_opts', status=ierr)
150 2 : if (ierr == 0) then
151 2 : read(unitn, aerosol_convproc_opts, iostat=ierr)
152 2 : if (ierr /= 0) then
153 0 : call endrun(subname // ':: ERROR reading namelist')
154 : end if
155 : end if
156 2 : close(unitn)
157 : end if
158 :
159 : ! Broadcast namelist variables
160 1536 : call mpi_bcast( convproc_do_gas, 1, mpi_logical, masterprocid, mpicom, ierr)
161 1536 : call mpi_bcast( deepconv_wetdep_history, 1, mpi_logical, masterprocid, mpicom, ierr)
162 1536 : call mpi_bcast( convproc_do_deep, 1, mpi_logical, masterprocid, mpicom, ierr)
163 1536 : call mpi_bcast( convproc_do_shallow, 1, mpi_logical, masterprocid, mpicom, ierr)
164 1536 : call mpi_bcast( convproc_do_evaprain_atonce, 1, mpi_logical, masterprocid, mpicom, ierr)
165 1536 : call mpi_bcast( convproc_pom_spechygro, 1, mpi_real8, masterprocid, mpicom, ierr)
166 1536 : call mpi_bcast( convproc_wup_max, 1, mpi_real8, masterprocid, mpicom, ierr)
167 :
168 1536 : if (masterproc) then
169 2 : write(iulog,*) subname//': convproc_do_gas = ', convproc_do_gas
170 2 : write(iulog,*) subname//': deepconv_wetdep_history = ',deepconv_wetdep_history
171 2 : write(iulog,*) subname//': convproc_do_deep = ',convproc_do_deep
172 2 : write(iulog,*) subname//': convproc_do_shallow = ',convproc_do_shallow
173 2 : write(iulog,*) subname//': convproc_do_evaprain_atonce = ',convproc_do_evaprain_atonce
174 2 : write(iulog,*) subname//': convproc_pom_spechygro = ',convproc_pom_spechygro
175 2 : write(iulog,*) subname//': convproc_wup_max = ', convproc_wup_max
176 : end if
177 :
178 1536 : end subroutine ma_convproc_readnl
179 :
180 : !=========================================================================================
181 :
182 1536 : subroutine ma_convproc_init
183 :
184 : integer :: n, l, ll
185 : integer :: npass_calc_updraft
186 : logical :: history_aerosol
187 :
188 : call phys_getopts( history_aerosol_out=history_aerosol, &
189 1536 : convproc_do_aer_out = convproc_do_aer )
190 :
191 : call addfld('SH_MFUP_MAX', horiz_only, 'A', 'kg/m2', &
192 1536 : 'Shallow conv. column-max updraft mass flux' )
193 : call addfld('SH_WCLDBASE', horiz_only, 'A', 'm/s', &
194 1536 : 'Shallow conv. cloudbase vertical velocity' )
195 : call addfld('SH_KCLDBASE', horiz_only, 'A', '1', &
196 1536 : 'Shallow conv. cloudbase level index' )
197 :
198 : call addfld('DP_MFUP_MAX', horiz_only, 'A', 'kg/m2', &
199 1536 : 'Deep conv. column-max updraft mass flux' )
200 : call addfld('DP_WCLDBASE', horiz_only, 'A', 'm/s', &
201 1536 : 'Deep conv. cloudbase vertical velocity' )
202 : call addfld('DP_KCLDBASE', horiz_only, 'A', '1', &
203 1536 : 'Deep conv. cloudbase level index' )
204 :
205 : ! output wet deposition fields to history
206 : ! I = in-cloud removal; E = precip-evap resuspension
207 : ! C = convective (total); D = deep convective
208 : ! note that the precip-evap resuspension includes that resulting from
209 : ! below-cloud removal, calculated in mz_aero_wet_intr
210 1536 : if (convproc_do_aer .and. apply_convproc_tend_to_ptend ) then
211 7680 : do n = 1, ntot_amode
212 36864 : do ll = 0, nspec_amode(n)
213 29184 : if (ll == 0) then
214 6144 : l = numptr_amode(n)
215 : else
216 23040 : l = lmassptr_amode(ll,n)
217 : end if
218 :
219 29184 : call addfld (trim(cnst_name(l))//'SFSEC', &
220 58368 : horiz_only, 'A','kg/m2/s','Wet deposition flux (precip evap, convective) at surface')
221 29184 : if (history_aerosol) then
222 0 : call add_default(trim(cnst_name(l))//'SFSEC', 1, ' ')
223 : end if
224 :
225 35328 : if ( deepconv_wetdep_history ) then
226 : call addfld (trim(cnst_name(l))//'SFSID', &
227 29184 : horiz_only, 'A','kg/m2/s','Wet deposition flux (incloud, deep convective) at surface')
228 : call addfld (trim(cnst_name(l))//'SFSED', &
229 29184 : horiz_only, 'A','kg/m2/s','Wet deposition flux (precip evap, deep convective) at surface')
230 29184 : if (history_aerosol) then
231 0 : call add_default(trim(cnst_name(l))//'SFSID', 1, ' ')
232 0 : call add_default(trim(cnst_name(l))//'SFSED', 1, ' ')
233 : end if
234 : end if
235 : end do
236 : end do
237 : end if
238 :
239 1536 : if ( history_aerosol .and. &
240 : ( convproc_do_aer .or. convproc_do_gas) ) then
241 0 : call add_default( 'SH_MFUP_MAX', 1, ' ' )
242 0 : call add_default( 'SH_WCLDBASE', 1, ' ' )
243 0 : call add_default( 'SH_KCLDBASE', 1, ' ' )
244 0 : call add_default( 'DP_MFUP_MAX', 1, ' ' )
245 0 : call add_default( 'DP_WCLDBASE', 1, ' ' )
246 0 : call add_default( 'DP_KCLDBASE', 1, ' ' )
247 : end if
248 :
249 1536 : fracis_idx = pbuf_get_index('FRACIS')
250 :
251 1536 : rprddp_idx = pbuf_get_index('RPRDDP')
252 1536 : rprdsh_idx = pbuf_get_index('RPRDSH')
253 1536 : nevapr_dpcu_idx = pbuf_get_index('NEVAPR_DPCU')
254 1536 : nevapr_shcu_idx = pbuf_get_index('NEVAPR_SHCU')
255 :
256 1536 : icwmrdp_idx = pbuf_get_index('ICWMRDP')
257 1536 : icwmrsh_idx = pbuf_get_index('ICWMRSH')
258 1536 : dp_frac_idx = pbuf_get_index('DP_FRAC')
259 1536 : sh_frac_idx = pbuf_get_index('SH_FRAC')
260 :
261 1536 : zm_mu_idx = pbuf_get_index('ZM_MU')
262 1536 : zm_eu_idx = pbuf_get_index('ZM_EU')
263 1536 : zm_du_idx = pbuf_get_index('ZM_DU')
264 1536 : zm_md_idx = pbuf_get_index('ZM_MD')
265 1536 : zm_ed_idx = pbuf_get_index('ZM_ED')
266 1536 : zm_dp_idx = pbuf_get_index('ZM_DP')
267 1536 : zm_dsubcld_idx = pbuf_get_index('ZM_DSUBCLD')
268 1536 : zm_jt_idx = pbuf_get_index('ZM_JT')
269 1536 : zm_maxg_idx = pbuf_get_index('ZM_MAXG')
270 1536 : zm_ideep_idx = pbuf_get_index('ZM_IDEEP')
271 :
272 1536 : cmfmc_sh_idx = pbuf_get_index('CMFMC_SH')
273 1536 : sh_e_ed_ratio_idx = pbuf_get_index('SH_E_ED_RATIO', istat)
274 :
275 1536 : if (masterproc ) then
276 :
277 2 : write(iulog,'(a,l12)') 'ma_convproc_init - convproc_do_aer = ', &
278 4 : convproc_do_aer
279 2 : write(iulog,'(a,l12)') 'ma_convproc_init - convproc_do_gas = ', &
280 4 : convproc_do_gas
281 2 : write(iulog,'(a,l12)') 'ma_convproc_init - use_cwaer_for_activate_maxsat = ', &
282 4 : use_cwaer_for_activate_maxsat
283 2 : write(iulog,'(a,l12)') 'ma_convproc_init - apply_convproc_tend_to_ptend = ', &
284 4 : apply_convproc_tend_to_ptend
285 2 : write(iulog,'(a,i12)') 'ma_convproc_init - convproc_method_activate = ', &
286 4 : convproc_method_activate
287 2 : write(iulog,'(a,i12)') 'ma_convproc_init - method1_activate_nlayers = ', &
288 4 : method1_activate_nlayers
289 2 : write(iulog,'(a,1pe12.4)') 'ma_convproc_init - method2_activate_smaxmax = ', &
290 4 : method2_activate_smaxmax
291 2 : write(iulog,'(a,i12)') 'ma_convproc_init - method_reduce_actfrac = ', &
292 4 : method_reduce_actfrac
293 2 : write(iulog,'(a,1pe12.4)') 'ma_convproc_init - factor_reduce_actfrac = ', &
294 4 : factor_reduce_actfrac
295 :
296 2 : npass_calc_updraft = 1
297 : if ( (method_reduce_actfrac == 2) .and. &
298 : (factor_reduce_actfrac >= 0.0_r8) .and. &
299 : (factor_reduce_actfrac <= 1.0_r8) ) npass_calc_updraft = 2
300 2 : write(iulog,'(a,i12)') 'ma_convproc_init - npass_calc_updraft = ', &
301 4 : npass_calc_updraft
302 :
303 : end if
304 :
305 1536 : aero_props_obj => modal_aerosol_properties()
306 1536 : if (.not.associated(aero_props_obj)) then
307 0 : call endrun('ma_convproc_init: modal_aerosol_properties constructor failed')
308 : end if
309 :
310 1536 : end subroutine ma_convproc_init
311 :
312 : !=========================================================================================
313 :
314 58824 : subroutine ma_convproc_intr( state, ptend, pbuf, ztodt, &
315 58824 : nsrflx_mzaer2cnvpr, qsrflx_mzaer2cnvpr, &
316 : aerdepwetis, dcondt_resusp3d )
317 : !-----------------------------------------------------------------------
318 : !
319 : ! Convective cloud processing (transport, activation/resuspension,
320 : ! wet removal) of aerosols and trace gases.
321 : ! (Currently no aqueous chemistry and no trace-gas wet removal)
322 : ! Does aerosols when convproc_do_aer is .true.
323 : ! Does trace gases when convproc_do_gas is .true.
324 : !
325 : ! Does deep and shallow convection
326 : ! Uses mass fluxes, cloud water, precip production from the
327 : ! convective cloud routines
328 : !
329 : ! Author: R. Easter
330 : !
331 : !-----------------------------------------------------------------------
332 :
333 :
334 : ! Arguments
335 : type(physics_state), intent(in ) :: state ! Physics state variables
336 : type(physics_ptend), intent(inout) :: ptend ! %lq set in aero_model_wetdep
337 : type(physics_buffer_desc), pointer :: pbuf(:)
338 : real(r8), intent(in) :: ztodt ! 2 delta t (model time increment)
339 :
340 : integer, intent(in) :: nsrflx_mzaer2cnvpr
341 : real(r8), intent(in) :: qsrflx_mzaer2cnvpr(pcols,pcnst,nsrflx_mzaer2cnvpr)
342 : real(r8), intent(inout) :: aerdepwetis(pcols,pcnst) ! aerosol wet deposition (interstitial)
343 : real(r8), intent(inout) :: dcondt_resusp3d(2*pcnst,pcols,pver)
344 :
345 : ! Local variables
346 : integer, parameter :: nsrflx = 5 ! last dimension of qsrflx
347 : integer :: l, ll, lchnk
348 : integer :: n, ncol
349 :
350 : real(r8) :: dqdt(pcols,pver,pcnst)
351 : real(r8) :: dt
352 : real(r8) :: qa(pcols,pver,pcnst), qb(pcols,pver,pcnst)
353 : real(r8) :: qsrflx(pcols,pcnst,nsrflx)
354 : real(r8) :: sflxic(pcols,pcnst)
355 : real(r8) :: sflxid(pcols,pcnst)
356 : real(r8) :: sflxec(pcols,pcnst)
357 : real(r8) :: sflxed(pcols,pcnst)
358 :
359 : logical :: dotend(pcnst)
360 : !-------------------------------------------------------------------------------------------------
361 :
362 : ! Initialize
363 58824 : lchnk = state%lchnk
364 58824 : ncol = state%ncol
365 58824 : dt = ztodt
366 :
367 58824 : hund_ovr_g = 100.0_r8/gravit
368 : ! used with zm_conv mass fluxes and delta-p
369 : ! for mu = [mbar/s], mu*hund_ovr_g = [kg/m2/s]
370 : ! for dp = [mbar] and q = [kg/kg], q*dp*hund_ovr_g = [kg/m2]
371 :
372 58824 : sflxic(:,:) = 0.0_r8
373 58824 : sflxid(:,:) = 0.0_r8
374 58824 : sflxec(:,:) = 0.0_r8
375 58824 : sflxed(:,:) = 0.0_r8
376 2470608 : do l = 1, pcnst
377 2470608 : if ( (cnst_species_class(l) == cnst_spec_class_aerosol) .and. ptend%lq(l) ) then
378 18662256 : sflxec(1:ncol,l) = qsrflx_mzaer2cnvpr(1:ncol,l,1)
379 18662256 : sflxed(1:ncol,l) = qsrflx_mzaer2cnvpr(1:ncol,l,2)
380 : end if
381 : end do
382 :
383 : ! prepare for deep conv processing
384 2470608 : do l = 1, pcnst
385 2470608 : if ( ptend%lq(l) ) then
386 : ! calc new q (after calcaersize and mz_aero_wet_intr)
387 1736707464 : qa(1:ncol,:,l) = state%q(1:ncol,:,l) + dt*ptend%q(1:ncol,:,l)
388 1736707464 : qb(1:ncol,:,l) = max( 0.0_r8, qa(1:ncol,:,l) )
389 : else
390 : ! use old q
391 2010924432 : qb(1:ncol,:,l) = state%q(1:ncol,:,l)
392 : end if
393 : end do
394 58824 : dqdt(:,:,:) = 0.0_r8
395 58824 : qsrflx(:,:,:) = 0.0_r8
396 :
397 58824 : if (convproc_do_aer .or. convproc_do_gas) then
398 :
399 : ! do deep conv processing
400 58824 : if (convproc_do_deep) then
401 : call ma_convproc_dp_intr( &
402 : state, pbuf, dt, &
403 58824 : qb, dqdt, dotend, nsrflx, qsrflx, dcondt_resusp3d )
404 :
405 :
406 : ! apply deep conv processing tendency and prepare for shallow conv processing
407 2470608 : do l = 1, pcnst
408 2411784 : if ( .not. dotend(l) ) cycle
409 :
410 : ! calc new q (after ma_convproc_dp_intr)
411 1736707464 : qa(1:ncol,:,l) = qb(1:ncol,:,l) + dt*dqdt(1:ncol,:,l)
412 1736707464 : qb(1:ncol,:,l) = max( 0.0_r8, qa(1:ncol,:,l) )
413 :
414 : if ( apply_convproc_tend_to_ptend ) then
415 : ! add dqdt onto ptend%q and set ptend%lq
416 1736707464 : ptend%q(1:ncol,:,l) = ptend%q(1:ncol,:,l) + dqdt(1:ncol,:,l)
417 1117656 : ptend%lq(l) = .true.
418 : end if
419 :
420 1117656 : if ((cnst_species_class(l) == cnst_spec_class_aerosol) .or. &
421 : (cnst_species_class(l) == cnst_spec_class_gas )) then
422 : ! these used for history file wetdep diagnostics
423 18662256 : sflxic(1:ncol,l) = sflxic(1:ncol,l) + qsrflx(1:ncol,l,4)
424 18662256 : sflxid(1:ncol,l) = sflxid(1:ncol,l) + qsrflx(1:ncol,l,4)
425 18662256 : sflxec(1:ncol,l) = sflxec(1:ncol,l) + qsrflx(1:ncol,l,5)
426 18662256 : sflxed(1:ncol,l) = sflxed(1:ncol,l) + qsrflx(1:ncol,l,5)
427 : end if
428 :
429 1176480 : if (cnst_species_class(l) == cnst_spec_class_aerosol) then
430 : ! this used for surface coupling
431 : aerdepwetis(1:ncol,l) = aerdepwetis(1:ncol,l) &
432 18662256 : + qsrflx(1:ncol,l,4) + qsrflx(1:ncol,l,5)
433 : end if
434 : end do
435 : end if
436 :
437 58824 : dqdt(:,:,:) = 0.0_r8
438 58824 : qsrflx(:,:,:) = 0.0_r8
439 58824 : if (convproc_do_shallow) then
440 : call ma_convproc_sh_intr( &
441 : state, pbuf, dt, &
442 0 : qb, dqdt, dotend, nsrflx, qsrflx, dcondt_resusp3d )
443 :
444 : ! apply shallow conv processing tendency
445 0 : do l = 1, pcnst
446 0 : if ( .not. dotend(l) ) cycle
447 :
448 : ! calc new q (after ma_convproc_sh_intr)
449 0 : qa(1:ncol,:,l) = qb(1:ncol,:,l) + dt*dqdt(1:ncol,:,l)
450 0 : qb(1:ncol,:,l) = max( 0.0_r8, qa(1:ncol,:,l) )
451 :
452 : if ( apply_convproc_tend_to_ptend ) then
453 : ! add dqdt onto ptend%q and set ptend%lq
454 0 : ptend%q(1:ncol,:,l) = ptend%q(1:ncol,:,l) + dqdt(1:ncol,:,l)
455 0 : ptend%lq(l) = .true.
456 : end if
457 :
458 0 : if ((cnst_species_class(l) == cnst_spec_class_aerosol) .or. &
459 : (cnst_species_class(l) == cnst_spec_class_gas )) then
460 0 : sflxic(1:ncol,l) = sflxic(1:ncol,l) + qsrflx(1:ncol,l,4)
461 0 : sflxec(1:ncol,l) = sflxec(1:ncol,l) + qsrflx(1:ncol,l,5)
462 : end if
463 :
464 0 : if (cnst_species_class(l) == cnst_spec_class_aerosol) then
465 : aerdepwetis(1:ncol,l) = aerdepwetis(1:ncol,l) &
466 0 : + qsrflx(1:ncol,l,4) + qsrflx(1:ncol,l,5)
467 : end if
468 :
469 : end do
470 : end if
471 :
472 : end if ! (convproc_do_aer .or. convproc_do_gas) then
473 :
474 :
475 58824 : if (convproc_do_aer .and. apply_convproc_tend_to_ptend ) then
476 294120 : do n = 1, ntot_amode
477 1411776 : do ll = 0, nspec_amode(n)
478 1117656 : if (ll == 0) then
479 235296 : l = numptr_amode(n)
480 : else
481 882360 : l = lmassptr_amode(ll,n)
482 : end if
483 :
484 1117656 : call outfld( trim(cnst_name(l))//'SFWET', aerdepwetis(:,l), pcols, lchnk )
485 1117656 : call outfld( trim(cnst_name(l))//'SFSIC', sflxic(:,l), pcols, lchnk )
486 1117656 : call outfld( trim(cnst_name(l))//'SFSEC', sflxec(:,l), pcols, lchnk )
487 :
488 1352952 : if ( deepconv_wetdep_history ) then
489 1117656 : call outfld( trim(cnst_name(l))//'SFSID', sflxid(:,l), pcols, lchnk )
490 1117656 : call outfld( trim(cnst_name(l))//'SFSED', sflxed(:,l), pcols, lchnk )
491 : end if
492 : end do
493 : end do
494 : end if
495 :
496 58824 : end subroutine ma_convproc_intr
497 :
498 : !=========================================================================================
499 :
500 58824 : subroutine ma_convproc_dp_intr( &
501 : state, pbuf, dt, &
502 58824 : q, dqdt, dotend, nsrflx, qsrflx, dcondt_resusp3d)
503 : !-----------------------------------------------------------------------
504 : !
505 : ! Convective cloud processing (transport, activation/resuspension,
506 : ! wet removal) of aerosols and trace gases.
507 : ! (Currently no aqueous chemistry and no trace-gas wet removal)
508 : ! Does aerosols when convproc_do_aer is .true.
509 : ! Does trace gases when convproc_do_gas is .true.
510 : !
511 : ! This routine does deep convection
512 : ! Uses mass fluxes, cloud water, precip production from the
513 : ! convective cloud routines
514 : !
515 : ! Author: R. Easter
516 : !
517 : !-----------------------------------------------------------------------
518 :
519 :
520 : ! Arguments
521 : type(physics_state), intent(in ) :: state ! Physics state variables
522 : type(physics_buffer_desc), pointer :: pbuf(:)
523 :
524 : real(r8), intent(in) :: dt ! delta t (model time increment)
525 :
526 : real(r8), intent(in) :: q(pcols,pver,pcnst)
527 : real(r8), intent(inout) :: dqdt(pcols,pver,pcnst)
528 : logical, intent(out) :: dotend(pcnst)
529 : integer, intent(in) :: nsrflx
530 : real(r8), intent(inout) :: qsrflx(pcols,pcnst,nsrflx)
531 : real(r8), intent(inout) :: dcondt_resusp3d(pcnst*2,pcols,pver)
532 :
533 : integer :: i
534 : integer :: itmpveca(pcols)
535 : integer :: l, lchnk, lun, ncol
536 : integer :: nstep
537 :
538 : real(r8) :: dpdry(pcols,pver) ! layer delta-p-dry (mb)
539 : real(r8) :: fracice(pcols,pver) ! Ice fraction of cloud droplets
540 : real(r8) :: qaa(pcols,pver,pcnst)
541 : real(r8) :: xx_mfup_max(pcols), xx_wcldbase(pcols), xx_kcldbase(pcols)
542 :
543 : ! physics buffer fields
544 58824 : real(r8), pointer :: fracis(:,:,:) ! fraction of transported species that are insoluble
545 58824 : real(r8), pointer :: rprddp(:,:) ! Deep conv precip production (kg/kg/s - grid avg)
546 58824 : real(r8), pointer :: evapcdp(:,:) ! Deep conv precip evaporation (kg/kg/s - grid avg)
547 58824 : real(r8), pointer :: icwmrdp(:,:) ! Deep conv cloud condensate (kg/kg - in cloud)
548 58824 : real(r8), pointer :: dp_frac(:,:) ! Deep conv cloud frac (0-1)
549 : ! mu, md, ..., ideep, lengath are all deep conv variables
550 58824 : real(r8), pointer :: mu(:,:) ! Updraft mass flux (positive) (pcols,pver)
551 58824 : real(r8), pointer :: md(:,:) ! Downdraft mass flux (negative) (pcols,pver)
552 58824 : real(r8), pointer :: du(:,:) ! Mass detrain rate from updraft (pcols,pver)
553 58824 : real(r8), pointer :: eu(:,:) ! Mass entrain rate into updraft (pcols,pver)
554 58824 : real(r8), pointer :: ed(:,:) ! Mass entrain rate into downdraft (pcols,pver)
555 : ! eu, ed, du are "d(massflux)/dp" and are all positive
556 58824 : real(r8), pointer :: dp(:,:) ! Delta pressure between interfaces (pcols,pver)
557 58824 : real(r8), pointer :: dsubcld(:) ! Delta pressure from cloud base to sfc (pcols)
558 :
559 58824 : integer, pointer :: jt(:) ! Index of cloud top for each column (pcols)
560 58824 : integer, pointer :: maxg(:) ! Index of cloud top for each column (pcols)
561 58824 : integer, pointer :: ideep(:) ! Gathering array (pcols)
562 : integer :: lengath ! Gathered min lon indices over which to operate
563 :
564 :
565 : ! Initialize
566 :
567 58824 : lchnk = state%lchnk
568 58824 : ncol = state%ncol
569 117648 : nstep = get_nstep()
570 58824 : lun = iulog
571 :
572 : ! Associate pointers with physics buffer fields
573 58824 : call pbuf_get_field(pbuf, fracis_idx, fracis)
574 58824 : call pbuf_get_field(pbuf, rprddp_idx, rprddp)
575 58824 : call pbuf_get_field(pbuf, nevapr_dpcu_idx, evapcdp)
576 58824 : call pbuf_get_field(pbuf, icwmrdp_idx, icwmrdp)
577 58824 : call pbuf_get_field(pbuf, dp_frac_idx, dp_frac)
578 58824 : call pbuf_get_field(pbuf, fracis_idx, fracis)
579 58824 : call pbuf_get_field(pbuf, zm_mu_idx, mu)
580 58824 : call pbuf_get_field(pbuf, zm_eu_idx, eu)
581 58824 : call pbuf_get_field(pbuf, zm_du_idx, du)
582 58824 : call pbuf_get_field(pbuf, zm_md_idx, md)
583 58824 : call pbuf_get_field(pbuf, zm_ed_idx, ed)
584 58824 : call pbuf_get_field(pbuf, zm_dp_idx, dp)
585 58824 : call pbuf_get_field(pbuf, zm_dsubcld_idx, dsubcld)
586 58824 : call pbuf_get_field(pbuf, zm_jt_idx, jt)
587 58824 : call pbuf_get_field(pbuf, zm_maxg_idx, maxg)
588 58824 : call pbuf_get_field(pbuf, zm_ideep_idx, ideep)
589 :
590 1000008 : lengath = count(ideep > 0)
591 58824 : if (lengath > ncol) lengath = ncol ! should not happen, but force it to not be larger than ncol for safety sake
592 :
593 58824 : fracice(:,:) = 0.0_r8
594 :
595 : ! initialize dpdry (units=mb), which is used for tracers of dry mixing ratio type
596 58824 : dpdry = 0._r8
597 223698 : do i = 1, lengath
598 15556980 : dpdry(i,:) = state%pdeldry(ideep(i),:)/100._r8
599 : end do
600 :
601 58824 : qaa = q
602 :
603 : ! turn on/off calculations for aerosols and trace gases
604 2470608 : do l = 1, pcnst
605 2411784 : dotend(l) = .false.
606 2470608 : if (cnst_species_class(l) == cnst_spec_class_aerosol) then
607 1117656 : if (convproc_do_aer) dotend(l) = .true.
608 1294128 : else if (cnst_species_class(l) == cnst_spec_class_gas) then
609 647064 : if (convproc_do_gas) dotend(l) = .true.
610 : end if
611 : end do
612 :
613 1000008 : itmpveca(:) = -1
614 :
615 : call ma_convproc_tend( &
616 : 'deep', &
617 : lchnk, pcnst, nstep, dt, &
618 : state%t, state%pmid, state%pdel, qaa, &
619 : mu, md, du, eu, &
620 : ed, dp, dpdry, jt, &
621 : maxg, ideep, 1, lengath, &
622 : dp_frac, icwmrdp, rprddp, evapcdp, &
623 : fracice, &
624 : dqdt, dotend, nsrflx, qsrflx, &
625 : xx_mfup_max, xx_wcldbase, xx_kcldbase, &
626 58824 : lun, itmpveca, dcondt_resusp3d )
627 :
628 58824 : call outfld( 'DP_MFUP_MAX', xx_mfup_max, pcols, lchnk )
629 58824 : call outfld( 'DP_WCLDBASE', xx_wcldbase, pcols, lchnk )
630 58824 : call outfld( 'DP_KCLDBASE', xx_kcldbase, pcols, lchnk )
631 :
632 58824 : end subroutine ma_convproc_dp_intr
633 :
634 :
635 :
636 : !=========================================================================================
637 0 : subroutine ma_convproc_sh_intr( &
638 : state, pbuf, dt, &
639 0 : q, dqdt, dotend, nsrflx, qsrflx, dcondt_resusp3d )
640 : !-----------------------------------------------------------------------
641 : !
642 : ! Purpose:
643 : ! Convective cloud processing (transport, activation/resuspension,
644 : ! wet removal) of aerosols and trace gases.
645 : ! (Currently no aqueous chemistry and no trace-gas wet removal)
646 : ! Does aerosols when convproc_do_aer is .true.
647 : ! Does trace gases when convproc_do_gas is .true.
648 : !
649 : ! This routine does shallow convection
650 : ! Uses mass fluxes, cloud water, precip production from the
651 : ! convective cloud routines
652 : !
653 : ! Author: R. Easter
654 : !
655 : !-----------------------------------------------------------------------
656 :
657 : ! Arguments
658 : type(physics_state), intent(in ) :: state ! Physics state variables
659 : type(physics_buffer_desc), pointer :: pbuf(:)
660 :
661 : real(r8), intent(in) :: dt ! delta t (model time increment)
662 :
663 : real(r8), intent(in) :: q(pcols,pver,pcnst)
664 : real(r8), intent(inout) :: dqdt(pcols,pver,pcnst)
665 : logical, intent(out) :: dotend(pcnst)
666 : integer, intent(in) :: nsrflx
667 : real(r8), intent(inout) :: qsrflx(pcols,pcnst,nsrflx)
668 : real(r8), intent(inout) :: dcondt_resusp3d(pcnst*2,pcols,pver)
669 :
670 : integer :: i
671 : integer :: itmpveca(pcols)
672 : integer :: k, kaa, kbb, kk
673 : integer :: l, lchnk, lun
674 : integer :: maxg_minval
675 : integer :: ncol, nstep
676 :
677 : real(r8) :: dpdry(pcols,pver) ! layer delta-p-dry (mb)
678 : real(r8) :: fracice(pcols,pver) ! Ice fraction of cloud droplets
679 : real(r8) :: qaa(pcols,pver,pcnst)
680 : real(r8) :: tmpa, tmpb
681 : real(r8) :: xx_mfup_max(pcols), xx_wcldbase(pcols), xx_kcldbase(pcols)
682 :
683 : ! variables that mimic the zm-deep counterparts
684 : real(r8) :: mu(pcols,pver) ! Updraft mass flux (positive)
685 : real(r8) :: md(pcols,pver) ! Downdraft mass flux (negative)
686 : real(r8) :: du(pcols,pver) ! Mass detrain rate from updraft
687 : real(r8) :: eu(pcols,pver) ! Mass entrain rate into updraft
688 : real(r8) :: ed(pcols,pver) ! Mass entrain rate into downdraft
689 : ! eu, ed, du are "d(massflux)/dp" and are all positive
690 : real(r8) :: dp(pcols,pver) ! Delta pressure between interfaces
691 :
692 : integer :: jt(pcols) ! Index of cloud top for each column
693 : integer :: maxg(pcols) ! Index of cloud bot for each column
694 : integer :: ideep(pcols) ! Gathering array
695 : integer :: lengath ! Gathered min lon indices over which to operate
696 :
697 : ! physics buffer fields
698 0 : real(r8), pointer :: rprdsh(:,:) ! Shallow conv precip production (kg/kg/s - grid avg)
699 0 : real(r8), pointer :: evapcsh(:,:) ! Shal conv precip evaporation (kg/kg/s - grid avg)
700 0 : real(r8), pointer :: icwmrsh(:,:) ! Shal conv cloud condensate (kg/kg - in cloud)
701 0 : real(r8), pointer :: sh_frac(:,:) ! Shal conv cloud frac (0-1)
702 :
703 0 : real(r8), pointer :: cmfmcsh(:,:) ! Shallow conv mass flux (pcols,pverp) (kg/m2/s)
704 0 : real(r8), pointer :: sh_e_ed_ratio(:,:) ! shallow conv [ent/(ent+det)] ratio (pcols,pver)
705 :
706 : ! Initialize
707 :
708 0 : lchnk = state%lchnk
709 0 : ncol = state%ncol
710 0 : nstep = get_nstep()
711 0 : lun = iulog
712 :
713 : ! Associate pointers with physics buffer fields
714 0 : call pbuf_get_field(pbuf, rprdsh_idx, rprdsh)
715 0 : call pbuf_get_field(pbuf, nevapr_shcu_idx, evapcsh)
716 0 : call pbuf_get_field(pbuf, icwmrsh_idx, icwmrsh)
717 0 : call pbuf_get_field(pbuf, sh_frac_idx, sh_frac)
718 0 : call pbuf_get_field(pbuf, cmfmc_sh_idx, cmfmcsh)
719 0 : if (sh_e_ed_ratio_idx .gt. 0) then
720 0 : call pbuf_get_field(pbuf, sh_e_ed_ratio_idx, sh_e_ed_ratio)
721 : end if
722 :
723 0 : fracice(:,:) = 0.0_r8
724 :
725 : ! create mass flux, entrainment, detrainment, and delta-p arrays
726 : ! with same units as the zm-deep
727 0 : mu(:,:) = 0.0_r8
728 0 : md(:,:) = 0.0_r8
729 0 : du(:,:) = 0.0_r8
730 0 : eu(:,:) = 0.0_r8
731 0 : ed(:,:) = 0.0_r8
732 0 : jt(:) = -1
733 0 : maxg(:) = -1
734 0 : ideep(:) = -1
735 0 : lengath = ncol
736 0 : maxg_minval = pver*2
737 :
738 : ! these dp and dpdry have units of mb
739 0 : dpdry(1:ncol,:) = state%pdeldry(1:ncol,:)/100._r8
740 0 : dp( 1:ncol,:) = state%pdel( 1:ncol,:)/100._r8
741 :
742 0 : do i = 1, ncol
743 0 : ideep(i) = i
744 :
745 : ! load updraft mass flux from cmfmcsh
746 0 : kk = 0
747 0 : do k = 2, pver
748 : ! if mass-flux < 1e-7 kg/m2/s ~= 1e-7 m/s ~= 1 cm/day, treat as zero
749 0 : if (cmfmcsh(i,k) >= 1.0e-7_r8) then
750 : ! mu has units of mb/s
751 0 : mu(i,k) = cmfmcsh(i,k) / hund_ovr_g
752 0 : kk = kk + 1
753 0 : if (kk == 1) jt(i) = k - 1
754 0 : maxg(i) = k
755 : end if
756 : end do
757 0 : if (kk <= 0) cycle ! current column has no convection
758 :
759 : ! extend below-cloud source region downwards (how far?)
760 0 : maxg_minval = min( maxg_minval, maxg(i) )
761 0 : kaa = maxg(i)
762 0 : kbb = min( kaa+4, pver )
763 : ! kbb = pver
764 0 : if (kbb > kaa) then
765 0 : tmpa = sum( dpdry(i,kaa:kbb) )
766 0 : do k = kaa+1, kbb
767 0 : mu(i,k) = mu(i,kaa)*sum( dpdry(i,k:kbb) )/tmpa
768 : end do
769 0 : maxg(i) = kbb
770 : end if
771 :
772 : ! calc ent / detrainment, using the [ent/(ent+det)] ratio from uw scheme
773 : ! which is equal to [fer_out/(fer_out+fdr_out)] (see uwshcu.F90)
774 : !
775 : ! note that the ratio is set to -1.0 (invalid) when both fer and fdr are very small
776 : ! and the ratio values are often strange (??) at topmost layer
777 : !
778 : ! for initial testing, impose a limit of
779 : ! entrainment <= 4 * (net entrainment), OR
780 : ! detrainment <= 4 * (net detrainment)
781 0 : do k = jt(i), maxg(i)
782 0 : if (k < pver) then
783 0 : tmpa = (mu(i,k) - mu(i,k+1))/dpdry(i,k)
784 : else
785 0 : tmpa = mu(i,k)/dpdry(i,k)
786 : end if
787 0 : if (sh_e_ed_ratio_idx .gt. 0) then
788 0 : tmpb = sh_e_ed_ratio(i,k)
789 : else
790 : tmpb = -1.0_r8 ! force ent only or det only
791 : end if
792 0 : if (tmpb < -1.0e-5_r8) then
793 : ! do ent only or det only
794 0 : if (tmpa >= 0.0_r8) then
795 : ! net entrainment
796 0 : eu(i,k) = tmpa
797 : else
798 : ! net detrainment
799 0 : du(i,k) = -tmpa
800 : end if
801 : else
802 0 : if (tmpa >= 0.0_r8) then
803 : ! net entrainment
804 0 : if (k >= kaa .or. tmpb < 0.0_r8) then
805 : ! layers at/below initial maxg, or sh_e_ed_ratio is invalid
806 0 : eu(i,k) = tmpa
807 : else
808 0 : tmpb = max( tmpb, 0.571_r8 )
809 0 : eu(i,k) = tmpa*(tmpb/(2.0_r8*tmpb - 1.0_r8))
810 0 : du(i,k) = eu(i,k) - tmpa
811 : end if
812 : else
813 : ! net detrainment
814 0 : tmpa = -tmpa
815 0 : if (k <= jt(i) .or. tmpb < 0.0_r8) then
816 : ! layers at/above jt (where ratio is strange??), or sh_e_ed_ratio is invalid
817 0 : du(i,k) = tmpa
818 : else
819 0 : tmpb = min( tmpb, 0.429_r8 )
820 0 : du(i,k) = tmpa*(1.0_r8 - tmpb)/(1.0_r8 - 2.0_r8*tmpb)
821 0 : eu(i,k) = du(i,k) - tmpa
822 : end if
823 : end if
824 : end if
825 : end do ! k
826 :
827 : end do ! i
828 :
829 0 : qaa = q
830 :
831 : ! turn on/off calculations for aerosols and trace gases
832 0 : do l = 1, pcnst
833 0 : dotend(l) = .false.
834 0 : if (cnst_species_class(l) == cnst_spec_class_aerosol) then
835 0 : if (convproc_do_aer) dotend(l) = .true.
836 0 : else if (cnst_species_class(l) == cnst_spec_class_gas) then
837 0 : if (convproc_do_gas) dotend(l) = .true.
838 : end if
839 : end do
840 :
841 :
842 0 : itmpveca(:) = -1
843 :
844 : call ma_convproc_tend( &
845 : 'uwsh', &
846 : lchnk, pcnst, nstep, dt, &
847 : state%t, state%pmid, state%pdel, qaa, &
848 : mu, md, du, eu, &
849 : ed, dp, dpdry, jt, &
850 : maxg, ideep, 1, lengath, &
851 : sh_frac, icwmrsh, rprdsh, evapcsh, &
852 : fracice, &
853 : dqdt, dotend, nsrflx, qsrflx, &
854 : xx_mfup_max, xx_wcldbase, xx_kcldbase, &
855 0 : lun, itmpveca, dcondt_resusp3d)
856 :
857 0 : call outfld( 'SH_MFUP_MAX', xx_mfup_max, pcols, lchnk )
858 0 : call outfld( 'SH_WCLDBASE', xx_wcldbase, pcols, lchnk )
859 0 : call outfld( 'SH_KCLDBASE', xx_kcldbase, pcols, lchnk )
860 :
861 0 : end subroutine ma_convproc_sh_intr
862 :
863 : !=========================================================================================
864 :
865 58824 : subroutine ma_convproc_tend( &
866 : convtype, &
867 : lchnk, ncnst, nstep, dt, &
868 58824 : t, pmid, pdel, q, &
869 : mu, md, du, eu, &
870 : ed, dp, dpdry, jt, &
871 : mx, ideep, il1g, il2g, &
872 : cldfrac, icwmr, rprd, evapc, &
873 : fracice, &
874 58824 : dqdt, doconvproc, nsrflx, qsrflx, &
875 : xx_mfup_max, xx_wcldbase, xx_kcldbase, &
876 : lun, idiag_in, dcondt_resusp3d )
877 :
878 : !-----------------------------------------------------------------------
879 : !
880 : ! Purpose:
881 : ! Convective transport of trace species.
882 : ! The trace species need not be conservative, and source/sink terms for
883 : ! activation, resuspension, aqueous chemistry and gas uptake, and
884 : ! wet removal are all applied.
885 : ! Currently this works with the ZM deep convection, but we should be able
886 : ! to adapt it for both Hack and McCaa shallow convection
887 : !
888 : !
889 : ! Compare to subr convproc which does conservative trace species.
890 : !
891 : ! A distinction between "moist" and "dry" mixing ratios is not currently made.
892 : ! (P. Rasch comment: Note that we are still assuming that the tracers are
893 : ! in a moist mixing ratio this will change soon)
894 :
895 : !
896 : ! Method:
897 : ! Computes tracer mixing ratios in updraft and downdraft "cells" in a
898 : ! Lagrangian manner, with source/sinks applied in the updraft other.
899 : ! Then computes grid-cell-mean tendencies by considering
900 : ! updraft and downdraft fluxes across layer boundaries
901 : ! environment subsidence/lifting fluxes across layer boundaries
902 : ! sources and sinks in the updraft
903 : ! resuspension of activated species in the grid-cell as a whole
904 : !
905 : ! Note1: A better estimate or calculation of either the updraft velocity
906 : ! or fractional area is needed.
907 : ! Note2: If updraft area is a small fraction of over cloud area,
908 : ! then aqueous chemistry is underestimated. These are both
909 : ! research areas.
910 : !
911 : ! Authors: O. Seland and R. Easter, based on convtran by P. Rasch
912 : !
913 : !-----------------------------------------------------------------------
914 :
915 : use modal_aero_data, only: cnst_name_cw, &
916 : lmassptr_amode, lmassptrcw_amode, &
917 : ntot_amode, ntot_amode, &
918 : nspec_amode, numptr_amode, numptrcw_amode
919 :
920 : implicit none
921 :
922 : !-----------------------------------------------------------------------
923 : !
924 : ! Input arguments
925 : !
926 : character(len=*), intent(in) :: convtype ! identifies the type of
927 : ! convection ("deep", "shcu")
928 : integer, intent(in) :: lchnk ! chunk identifier
929 : integer, intent(in) :: ncnst ! number of tracers to transport
930 : integer, intent(in) :: nstep ! Time step index
931 : real(r8), intent(in) :: dt ! Model timestep
932 : real(r8), intent(in) :: t(pcols,pver) ! Temperature
933 : real(r8), intent(in) :: pmid(pcols,pver) ! Pressure at model levels
934 : real(r8), intent(in) :: pdel(pcols,pver) ! Pressure thickness of levels
935 : real(r8), intent(in) :: q(pcols,pver,ncnst) ! Tracer array including moisture
936 :
937 : real(r8), intent(in) :: mu(pcols,pver) ! Updraft mass flux (positive)
938 : real(r8), intent(in) :: md(pcols,pver) ! Downdraft mass flux (negative)
939 : real(r8), intent(in) :: du(pcols,pver) ! Mass detrain rate from updraft
940 : real(r8), intent(in) :: eu(pcols,pver) ! Mass entrain rate into updraft
941 : real(r8), intent(in) :: ed(pcols,pver) ! Mass entrain rate into downdraft
942 : ! *** note1 - mu, md, eu, ed, du, dp, dpdry are GATHERED ARRAYS ***
943 : ! *** note2 - mu and md units are (mb/s), which is used in the zm_conv code
944 : ! - eventually these should be changed to (kg/m2/s)
945 : ! *** note3 - eu, ed, du are "d(massflux)/dp" (with dp units = mb), and are all >= 0
946 :
947 : real(r8), intent(in) :: dp(pcols,pver) ! Delta pressure between interfaces (mb)
948 : real(r8), intent(in) :: dpdry(pcols,pver) ! Delta dry-pressure (mb)
949 : ! real(r8), intent(in) :: dsubcld(pcols) ! Delta pressure from cloud base to sfc
950 : integer, intent(in) :: jt(pcols) ! Index of cloud top for each column
951 : integer, intent(in) :: mx(pcols) ! Index of cloud top for each column
952 : integer, intent(in) :: ideep(pcols) ! Gathering array indices
953 : integer, intent(in) :: il1g ! Gathered min lon indices over which to operate
954 : integer, intent(in) :: il2g ! Gathered max lon indices over which to operate
955 : ! *** note4 -- for il1g <= i <= il2g, icol = ideep(i) is the "normal" chunk column index
956 :
957 : real(r8), intent(in) :: cldfrac(pcols,pver) ! Convective cloud fractional area
958 : real(r8), intent(in) :: icwmr(pcols,pver) ! Convective cloud water from zhang
959 : real(r8), intent(in) :: rprd(pcols,pver) ! Convective precipitation formation rate
960 : real(r8), intent(in) :: evapc(pcols,pver) ! Convective precipitation evaporation rate
961 : real(r8), intent(in) :: fracice(pcols,pver) ! Ice fraction of cloud droplets
962 :
963 : real(r8), intent(out):: dqdt(pcols,pver,ncnst) ! Tracer tendency array
964 : logical, intent(in) :: doconvproc(ncnst) ! flag for doing convective transport
965 : integer, intent(in) :: nsrflx ! last dimension of qsrflx
966 : real(r8), intent(out):: qsrflx(pcols,pcnst,nsrflx)
967 : ! process-specific column tracer tendencies
968 : ! (1=activation, 2=resuspension, 3=aqueous rxn,
969 : ! 4=wet removal, 5=renaming)
970 : real(r8), intent(out) :: xx_mfup_max(pcols)
971 : real(r8), intent(out) :: xx_wcldbase(pcols)
972 : real(r8), intent(out) :: xx_kcldbase(pcols)
973 : integer, intent(in) :: lun ! unit number for diagnostic output
974 : integer, intent(in) :: idiag_in(pcols) ! flag for diagnostic output
975 : real(r8), intent(inout) :: dcondt_resusp3d(pcnst*2,pcols,pver)
976 :
977 : !--------------------------Local Variables------------------------------
978 :
979 : ! cloudborne aerosol, so the arrays are dimensioned with pcnst_extd = pcnst*2
980 : integer, parameter :: pcnst_extd = pcnst*2
981 :
982 : integer :: i, icol ! Work index
983 : integer :: iconvtype ! 1=deep, 2=uw shallow
984 : integer :: idiag_act ! Work index
985 : integer :: iflux_method ! 1=as in convtran (deep), 2=simpler
986 : integer :: ipass_calc_updraft
987 : integer :: itmpa, itmpb ! Work variable
988 : integer :: j, jtsub ! Work index
989 : integer :: k ! Work index
990 : integer :: kactcnt ! Counter for no. of levels having activation
991 : integer :: kactcntb ! Counter for activation diagnostic output
992 : integer :: kactfirst ! Lowest layer with activation (= cloudbase)
993 : integer :: kbot ! Cloud-flux bottom layer for current i (=mx(i))
994 : integer :: kbot_prevap ! Lowest layer for doing resuspension from evaporating precip
995 : integer :: ktop ! Cloud-flux top layer for current i (=jt(i))
996 : ! Layers between kbot,ktop have mass fluxes
997 : ! but not all have cloud water, because the
998 : ! updraft starts below the cloud base
999 : integer :: km1, km1x ! Work index
1000 : integer :: kp1, kp1x ! Work index
1001 : integer :: l, ll, la, lc ! Work index
1002 : integer :: m, n ! Work index
1003 : integer :: merr ! number of errors (i.e., failed diagnostics)
1004 : ! for current column
1005 : integer :: nerr ! number of errors for entire run
1006 : integer :: nerrmax ! maximum number of errors to report
1007 : integer :: ncnst_extd
1008 : integer :: npass_calc_updraft
1009 : integer :: ntsub !
1010 :
1011 : logical do_act_this_lev ! flag for doing activation at current level
1012 : logical doconvproc_extd(pcnst_extd) ! flag for doing convective transport
1013 :
1014 : real(r8) aqfrac(pcnst_extd) ! aqueous fraction of constituent in updraft
1015 : real(r8) cldfrac_i(pver) ! cldfrac at current i (with adjustments)
1016 :
1017 : real(r8) chat(pcnst_extd,pverp) ! mix ratio in env at interfaces
1018 : real(r8) cond(pcnst_extd,pverp) ! mix ratio in downdraft at interfaces
1019 : real(r8) const(pcnst_extd,pver) ! gathered tracer array
1020 : real(r8) conu(pcnst_extd,pverp) ! mix ratio in updraft at interfaces
1021 :
1022 : real(r8) dcondt(pcnst_extd,pver) ! grid-average TMR tendency for current column
1023 : real(r8) dcondt_prevap(pcnst_extd,pver) ! portion of dcondt from precip evaporation
1024 : real(r8) dcondt_resusp(pcnst_extd,pver) ! portion of dcondt from resuspension
1025 :
1026 : real(r8) dcondt_wetdep(pcnst_extd,pver) ! portion of dcondt from wet deposition
1027 : real(r8) dconudt_activa(pcnst_extd,pverp) ! d(conu)/dt by activation
1028 : real(r8) dconudt_aqchem(pcnst_extd,pverp) ! d(conu)/dt by aqueous chem
1029 : real(r8) dconudt_wetdep(pcnst_extd,pverp) ! d(conu)/dt by wet removal
1030 :
1031 : real(r8) maxflux(pcnst_extd) ! maximum (over layers) of fluxin and fluxout
1032 : real(r8) maxflux2(pcnst_extd) ! ditto but computed using method-2 fluxes
1033 : real(r8) maxprevap(pcnst_extd) ! maximum (over layers) of dcondt_prevap*dp
1034 : real(r8) maxresusp(pcnst_extd) ! maximum (over layers) of dcondt_resusp*dp
1035 : real(r8) maxsrce(pcnst_extd) ! maximum (over layers) of netsrce
1036 :
1037 : real(r8) sumflux(pcnst_extd) ! sum (over layers) of netflux
1038 : real(r8) sumflux2(pcnst_extd) ! ditto but computed using method-2 fluxes
1039 : real(r8) sumsrce(pcnst_extd) ! sum (over layers) of dp*netsrce
1040 : real(r8) sumchng(pcnst_extd) ! sum (over layers) of dp*dcondt
1041 : real(r8) sumchng3(pcnst_extd) ! ditto but after call to resusp_conv
1042 : real(r8) sumactiva(pcnst_extd) ! sum (over layers) of dp*dconudt_activa
1043 : real(r8) sumaqchem(pcnst_extd) ! sum (over layers) of dp*dconudt_aqchem
1044 : real(r8) sumprevap(pcnst_extd) ! sum (over layers) of dp*dcondt_prevap
1045 : real(r8) sumresusp(pcnst_extd) ! sum (over layers) of dp*dcondt_resusp
1046 : real(r8) sumwetdep(pcnst_extd) ! sum (over layers) of dp*dconudt_wetdep
1047 :
1048 : real(r8) cabv ! mix ratio of constituent above
1049 : real(r8) cbel ! mix ratio of constituent below
1050 : real(r8) cdifr ! normalized diff between cabv and cbel
1051 : real(r8) cdt(pver) ! (in-updraft first order wet removal rate) * dt
1052 : real(r8) clw_cut ! threshold clw value for doing updraft
1053 : ! transformation and removal
1054 : real(r8) courantmax ! maximum courant no.
1055 : real(r8) dddp(pver) ! dd(i,k)*dp(i,k) at current i
1056 : real(r8) dp_i(pver) ! dp(i,k) at current i
1057 : real(r8) dt_u(pver) ! lagrangian transport time in the updraft
1058 : real(r8) dudp(pver) ! du(i,k)*dp(i,k) at current i
1059 : real(r8) dqdt_i(pver,pcnst) ! dqdt(i,k,m) at current i
1060 : real(r8) dtsub ! dt/ntsub
1061 : real(r8) dz ! working layer thickness (m)
1062 : real(r8) eddp(pver) ! ed(i,k)*dp(i,k) at current i
1063 : real(r8) eudp(pver) ! eu(i,k)*dp(i,k) at current i
1064 : real(r8) expcdtm1 ! a work variable
1065 : real(r8) fa_u(pver) ! fractional area of in the updraft
1066 : real(r8) fa_u_dp ! current fa_u(k)*dp_i(k)
1067 : real(r8) f_ent ! fraction of the "before-detrainment" updraft
1068 : ! massflux at k/k-1 interface resulting from
1069 : ! entrainment of level k air
1070 : real(r8) fluxin ! a work variable
1071 : real(r8) fluxout ! a work variable
1072 : real(r8) maxc ! a work variable
1073 : real(r8) mbsth ! Threshold for mass fluxes
1074 : real(r8) minc ! a work variable
1075 : real(r8) md_m_eddp ! a work variable
1076 : real(r8) md_i(pverp) ! md(i,k) at current i (note pverp dimension)
1077 : real(r8) md_x(pverp) ! md(i,k) at current i (note pverp dimension)
1078 : real(r8) mu_i(pverp) ! mu(i,k) at current i (note pverp dimension)
1079 : real(r8) mu_x(pverp) ! mu(i,k) at current i (note pverp dimension)
1080 : ! md_i, md_x, mu_i, mu_x are all "dry" mass fluxes
1081 : ! the mu_x/md_x are initially calculated from the incoming mu/md by applying dp/dpdry
1082 : ! the mu_i/md_i are next calculated by applying the mbsth threshold
1083 : real(r8) mu_p_eudp(pver) ! = mu_i(kp1) + eudp(k)
1084 : real(r8) netflux ! a work variable
1085 : real(r8) netsrce ! a work variable
1086 : real(r8) q_i(pver,pcnst) ! q(i,k,m) at current i
1087 117648 : real(r8) qsrflx_i(pcnst,nsrflx) ! qsrflx(i,m,n) at current i
1088 : real(r8) relerr_cut ! relative error criterion for diagnostics
1089 : real(r8) rhoair_i(pver) ! air density at current i
1090 : real(r8) small ! a small number
1091 : real(r8) tmpa, tmpb ! work variables
1092 : real(r8) tmpf ! work variables
1093 : real(r8) tmpveca(pcnst_extd) ! work variables
1094 : real(r8) tmpmata(pcnst_extd,3) ! work variables
1095 : real(r8) xinv_ntsub ! 1.0/ntsub
1096 : real(r8) wup(pver) ! working updraft velocity (m/s)
1097 : real(r8) zmagl(pver) ! working height above surface (m)
1098 : real(r8) zkm ! working height above surface (km)
1099 :
1100 : character(len=16) :: cnst_name_extd(pcnst_extd)
1101 :
1102 : !Fractional area of ensemble mean updrafts in ZM scheme set to 0.01
1103 : !Chosen to reproduce vertical vecocities in GATEIII GIGALES (Khairoutdinov etal 2009, JAMES)
1104 : real(r8), parameter :: zm_areafrac = 0.01_r8
1105 : !-----------------------------------------------------------------------
1106 : !
1107 :
1108 : ! if (nstep > 1) call endrun()
1109 :
1110 58824 : if (convtype == 'deep') then
1111 : iconvtype = 1
1112 : iflux_method = 1
1113 0 : else if (convtype == 'uwsh') then
1114 : iconvtype = 2
1115 : iflux_method = 2
1116 : else
1117 0 : call endrun( '*** ma_convproc_tend -- convtype is not |deep| or |uwsh|' )
1118 : end if
1119 :
1120 58824 : nerr = 0
1121 58824 : nerrmax = 99
1122 :
1123 58824 : ncnst_extd = pcnst_extd
1124 58824 : dcondt_resusp3d(:,:,:) = 0._r8
1125 :
1126 58824 : small = 1.e-36_r8
1127 : ! mbsth is the threshold below which we treat the mass fluxes as zero (in mb/s)
1128 58824 : mbsth = 1.e-15_r8
1129 :
1130 205354584 : qsrflx(:,:,:) = 0.0_r8
1131 3815501112 : dqdt(:,:,:) = 0.0_r8
1132 58824 : xx_mfup_max(:) = 0.0_r8
1133 58824 : xx_wcldbase(:) = 0.0_r8
1134 58824 : xx_kcldbase(:) = 0.0_r8
1135 :
1136 58824 : wup(:) = 0.0_r8
1137 :
1138 : ! set doconvproc_extd (extended array) values
1139 : ! inititialize aqfrac to 1.0 for activated aerosol species, 0.0 otherwise
1140 58824 : doconvproc_extd(:) = .false.
1141 2411784 : doconvproc_extd(2:ncnst) = doconvproc(2:ncnst)
1142 58824 : aqfrac(:) = 0.0_r8
1143 294120 : do n = 1, ntot_amode
1144 1411776 : do ll = 0, nspec_amode(n)
1145 1117656 : if (ll == 0) then
1146 235296 : la = numptr_amode(n)
1147 235296 : lc = numptrcw_amode(n) + pcnst
1148 : else
1149 882360 : la = lmassptr_amode(ll,n)
1150 882360 : lc = lmassptrcw_amode(ll,n) + pcnst
1151 : end if
1152 1352952 : if ( doconvproc(la) ) then
1153 1117656 : doconvproc_extd(lc) = .true.
1154 1117656 : aqfrac(lc) = 1.0_r8
1155 : end if
1156 : enddo
1157 : enddo ! n
1158 :
1159 4882392 : do l = 1, pcnst_extd
1160 4882392 : if (l <= pcnst) then
1161 2411784 : cnst_name_extd(l) = cnst_name(l)
1162 : else
1163 2411784 : cnst_name_extd(l) = trim(cnst_name(l-pcnst)) // '_cw'
1164 : end if
1165 : end do
1166 :
1167 :
1168 : ! Loop ever each column that has convection
1169 : ! *** i is index to gathered arrays; ideep(i) is index to "normal" chunk arrays
1170 : i_loop_main_aa: &
1171 223698 : do i = il1g, il2g
1172 164874 : icol = ideep(i)
1173 :
1174 :
1175 164874 : if ( (jt(i) <= 0) .and. (mx(i) <= 0) .and. (iconvtype /= 1) ) then
1176 : ! shallow conv case with jt,mx <= 0, which means there is no shallow conv
1177 : ! in this column -- skip this column
1178 : cycle i_loop_main_aa
1179 :
1180 164874 : else if ( (jt(i) < 1) .or. (mx(i) > pver) .or. (jt(i) > mx(i)) ) then
1181 : ! invalid cloudtop and cloudbase indices -- skip this column
1182 0 : write(lun,9010) 'illegal jt, mx', convtype, lchnk, icol, i, &
1183 0 : jt(i), mx(i)
1184 : 9010 format( '*** ma_convproc_tend error -- ', a, 5x, 'convtype = ', a / &
1185 : '*** lchnk, icol, il, jt, mx = ', 5(1x,i10) )
1186 0 : cycle i_loop_main_aa
1187 :
1188 164874 : else if (jt(i) == mx(i)) then
1189 : ! cloudtop = cloudbase (1 layer cloud) -- skip this column
1190 0 : write(lun,9010) 'jt == mx', convtype, lchnk, icol, i, jt(i), mx(i)
1191 0 : cycle i_loop_main_aa
1192 :
1193 : end if
1194 :
1195 :
1196 : !
1197 : ! cloudtop and cloudbase indices are valid so proceed with calculations
1198 : !
1199 :
1200 : ! Load dp_i and cldfrac_i, and calc rhoair_i
1201 15498156 : do k = 1, pver
1202 15333282 : dp_i(k) = dpdry(i,k)
1203 15333282 : cldfrac_i(k) = cldfrac(icol,k)
1204 15498156 : rhoair_i(k) = pmid(icol,k)/(rair*t(icol,k))
1205 : end do
1206 :
1207 : ! Calc dry mass fluxes
1208 : ! This is approximate because the updraft air is has different temp and qv than
1209 : ! the grid mean, but the whole convective parameterization is highly approximate
1210 164874 : mu_x(:) = 0.0_r8
1211 164874 : md_x(:) = 0.0_r8
1212 : ! (eu-du) = d(mu)/dp -- integrate upwards, multiplying by dpdry
1213 15498156 : do k = pver, 1, -1
1214 15333282 : mu_x(k) = mu_x(k+1) + (eu(i,k)-du(i,k))*dp_i(k)
1215 15498156 : xx_mfup_max(icol) = max( xx_mfup_max(icol), mu_x(k) )
1216 : end do
1217 : ! (ed) = d(md)/dp -- integrate downwards, multiplying by dpdry
1218 15333282 : do k = 2, pver
1219 15333282 : md_x(k) = md_x(k-1) - ed(i,k-1)*dp_i(k-1)
1220 : end do
1221 :
1222 : ! Load mass fluxes over cloud layers
1223 : ! (Note - use of arrays dimensioned k=1,pver+1 simplifies later coding)
1224 : ! Zero out values below threshold
1225 : ! Zero out values at "top of cloudtop", "base of cloudbase"
1226 164874 : ktop = jt(i)
1227 164874 : kbot = mx(i)
1228 : ! usually the updraft ( & downdraft) start ( & end ) at kbot=pver, but sometimes kbot < pver
1229 : ! transport, activation, resuspension, and wet removal only occur between kbot >= k >= ktop
1230 : ! resuspension from evaporating precip can occur at k > kbot when kbot < pver
1231 164874 : kbot_prevap = pver
1232 164874 : mu_i(:) = 0.0_r8
1233 164874 : md_i(:) = 0.0_r8
1234 2539434 : do k = ktop+1, kbot
1235 2374560 : mu_i(k) = mu_x(k)
1236 2374560 : if (mu_i(k) <= mbsth) mu_i(k) = 0.0_r8
1237 2374560 : md_i(k) = md_x(k)
1238 2539434 : if (md_i(k) >= -mbsth) md_i(k) = 0.0_r8
1239 : end do
1240 164874 : mu_i(ktop) = 0.0_r8
1241 164874 : md_i(ktop) = 0.0_r8
1242 164874 : mu_i(kbot+1) = 0.0_r8
1243 164874 : md_i(kbot+1) = 0.0_r8
1244 :
1245 : ! Compute updraft and downdraft "entrainment*dp" from eu and ed
1246 : ! Compute "detrainment*dp" from mass conservation
1247 164874 : eudp(:) = 0.0_r8
1248 164874 : dudp(:) = 0.0_r8
1249 164874 : eddp(:) = 0.0_r8
1250 164874 : dddp(:) = 0.0_r8
1251 164874 : courantmax = 0.0_r8
1252 2704308 : do k = ktop, kbot
1253 2539434 : if ((mu_i(k) > 0) .or. (mu_i(k+1) > 0)) then
1254 2533070 : if (du(i,k) <= 0.0_r8) then
1255 2287521 : eudp(k) = mu_i(k) - mu_i(k+1)
1256 : else
1257 245549 : eudp(k) = max( eu(i,k)*dp_i(k), 0.0_r8 )
1258 245549 : dudp(k) = (mu_i(k+1) + eudp(k)) - mu_i(k)
1259 245549 : if (dudp(k) < 1.0e-12_r8*eudp(k)) then
1260 0 : eudp(k) = mu_i(k) - mu_i(k+1)
1261 0 : dudp(k) = 0.0_r8
1262 : end if
1263 : end if
1264 : end if
1265 2539434 : if ((md_i(k) < 0) .or. (md_i(k+1) < 0)) then
1266 2102215 : eddp(k) = max( ed(i,k)*dp_i(k), 0.0_r8 )
1267 2102215 : dddp(k) = (md_i(k+1) + eddp(k)) - md_i(k)
1268 2102215 : if (dddp(k) < 1.0e-12_r8*eddp(k)) then
1269 1942861 : eddp(k) = md_i(k) - md_i(k+1)
1270 1942861 : dddp(k) = 0.0_r8
1271 : end if
1272 : end if
1273 : ! courantmax = max( courantmax, (eudp(k)+eddp(k))*dt/dp_i(k) ) ! old version - incorrect
1274 2704308 : courantmax = max( courantmax, ( mu_i(k+1)+eudp(k)-md_i(k)+eddp(k) )*dt/dp_i(k) )
1275 : end do ! k
1276 :
1277 : ! number of time substeps needed to maintain "courant number" <= 1
1278 164874 : ntsub = 1
1279 164874 : if (courantmax > (1.0_r8 + 1.0e-6_r8)) then
1280 2637 : ntsub = 1 + int( courantmax )
1281 : end if
1282 164874 : xinv_ntsub = 1.0_r8/ntsub
1283 164874 : dtsub = dt*xinv_ntsub
1284 164874 : courantmax = courantmax*xinv_ntsub
1285 :
1286 : ! zmagl(k) = height above surface for middle of level k
1287 164874 : zmagl(pver) = 0.0_r8
1288 15498156 : do k = pver, 1, -1
1289 15333282 : if (k < pver) then
1290 15168408 : zmagl(k) = zmagl(k+1) + 0.5_r8*dz
1291 : end if
1292 15333282 : dz = dp_i(k)*hund_ovr_g/rhoair_i(k)
1293 15498156 : zmagl(k) = zmagl(k) + 0.5_r8*dz
1294 : end do
1295 :
1296 : ! load tracer mixing ratio array, which will be updated at the end of each jtsub interation
1297 635589270 : q_i(1:pver,1:pcnst) = q(icol,1:pver,1:pcnst)
1298 :
1299 : !
1300 : ! when method_reduce_actfrac = 2, need to do the updraft calc twice
1301 : ! (1st to get non-adjusted activation amount, 2nd to apply reduction factor)
1302 164874 : npass_calc_updraft = 1
1303 : if ( (method_reduce_actfrac == 2) .and. &
1304 : (factor_reduce_actfrac >= 0.0_r8) .and. &
1305 : (factor_reduce_actfrac <= 1.0_r8) ) npass_calc_updraft = 2
1306 :
1307 :
1308 : jtsub_loop_main_aa: &
1309 391209 : do jtsub = 1, ntsub
1310 :
1311 :
1312 : ipass_calc_updraft_loop: &
1313 499896 : do ipass_calc_updraft = 1, npass_calc_updraft
1314 :
1315 :
1316 167511 : if (idiag_in(icol) > 0) &
1317 0 : write(lun,'(/a,3x,a,1x,i9,5i5)') 'qakr - convtype,lchnk,i,jt,mx,jtsub,ipass=', &
1318 0 : trim(convtype), lchnk, icol, jt(i), mx(i), jtsub, ipass_calc_updraft
1319 :
1320 35344821 : qsrflx_i(:,:) = 0.0_r8
1321 167511 : dqdt_i(:,:) = 0.0_r8
1322 :
1323 167511 : const(:,:) = 0.0_r8 ! zero cloud-phase species
1324 167511 : chat(:,:) = 0.0_r8 ! zero cloud-phase species
1325 167511 : conu(:,:) = 0.0_r8
1326 167511 : cond(:,:) = 0.0_r8
1327 :
1328 167511 : dcondt(:,:) = 0.0_r8
1329 167511 : dcondt_resusp(:,:) = 0.0_r8
1330 167511 : dcondt_wetdep(:,:) = 0.0_r8
1331 167511 : dcondt_prevap(:,:) = 0.0_r8
1332 167511 : dconudt_aqchem(:,:) = 0.0_r8
1333 167511 : dconudt_wetdep(:,:) = 0.0_r8
1334 : ! only initialize the activation tendency on ipass=1
1335 167511 : if (ipass_calc_updraft == 1) dconudt_activa(:,:) = 0.0_r8
1336 :
1337 : ! initialize mixing ratio arrays (chat, const, conu, cond)
1338 6867951 : do m = 2, ncnst
1339 6867951 : if ( doconvproc_extd(m) ) then
1340 :
1341 : ! Gather up the constituent
1342 299174646 : do k = 1,pver
1343 299174646 : const(m,k) = q_i(k,m)
1344 : end do
1345 :
1346 : ! From now on work only with gathered data
1347 : ! Interpolate environment tracer values to interfaces
1348 299174646 : do k = 1,pver
1349 295991937 : km1 = max(1,k-1)
1350 295991937 : minc = min(const(m,km1),const(m,k))
1351 295991937 : maxc = max(const(m,km1),const(m,k))
1352 295991937 : if (minc < 0) then
1353 : cdifr = 0._r8
1354 : else
1355 295991937 : cdifr = abs(const(m,k)-const(m,km1))/max(maxc,small)
1356 : endif
1357 :
1358 : ! If the two layers differ significantly use a geometric averaging procedure
1359 : ! But only do that for deep convection. For shallow, use the simple
1360 : ! averaging which is used in subr cmfmca
1361 295991937 : if (iconvtype /= 1) then
1362 0 : chat(m,k) = 0.5_r8* (const(m,k)+const(m,km1))
1363 295991937 : else if (cdifr > 1.E-6_r8) then
1364 : ! if (cdifr > 1.E-6) then
1365 292379023 : cabv = max(const(m,km1),maxc*1.e-12_r8)
1366 292379023 : cbel = max(const(m,k),maxc*1.e-12_r8)
1367 292379023 : chat(m,k) = log(cabv/cbel)/(cabv-cbel)*cabv*cbel
1368 : else ! Small diff, so just arithmetic mean
1369 3612914 : chat(m,k) = 0.5_r8* (const(m,k)+const(m,km1))
1370 : end if
1371 :
1372 : ! Set provisional up and down draft values, and tendencies
1373 295991937 : conu(m,k) = chat(m,k)
1374 299174646 : cond(m,k) = chat(m,k)
1375 : end do ! k
1376 :
1377 : ! Values at surface inferface == values in lowest layer
1378 3182709 : chat(m,pver+1) = const(m,pver)
1379 3182709 : conu(m,pver+1) = const(m,pver)
1380 3182709 : cond(m,pver+1) = const(m,pver)
1381 : end if
1382 : end do ! m
1383 :
1384 :
1385 :
1386 :
1387 : ! Compute updraft mixing ratios from cloudbase to cloudtop
1388 : ! No special treatment is needed at k=pver because arrays
1389 : ! are dimensioned 1:pver+1
1390 : ! A time-split approach is used. First, entrainment is applied to produce
1391 : ! an initial conu(m,k) from conu(m,k+1). Next, chemistry/physics are
1392 : ! applied to the initial conu(m,k) to produce a final conu(m,k).
1393 : ! Detrainment from the updraft uses this final conu(m,k).
1394 : ! Note that different time-split approaches would give somewhat different
1395 : ! results
1396 167511 : kactcnt = 0 ; kactcntb = 0 ; kactfirst = 1
1397 : k_loop_main_bb: &
1398 2741894 : do k = kbot, ktop, -1
1399 2574383 : kp1 = k+1
1400 :
1401 : ! cldfrac = conv cloud fractional area. This could represent anvil cirrus area,
1402 : ! and may not useful for aqueous chem and wet removal calculations
1403 2574383 : cldfrac_i(k) = max( cldfrac_i(k), 0.005_r8 )
1404 : ! mu_p_eudp(k) = updraft massflux at k, without detrainment between kp1,k
1405 2574383 : mu_p_eudp(k) = mu_i(kp1) + eudp(k)
1406 :
1407 2574383 : fa_u(k) = 0.0_r8 !BSINGH(10/15/2014): Initialized so that it has a value if the following "if" check yeilds .false.
1408 2741894 : if (mu_p_eudp(k) > mbsth) then
1409 : ! if (mu_p_eudp(k) <= mbsth) the updraft mass flux is negligible at base and top
1410 : ! of current layer,
1411 : ! so current layer is a "gap" between two unconnected updrafts,
1412 : ! so essentially skip all the updraft calculations for this layer
1413 :
1414 : ! First apply changes from entrainment
1415 2568019 : f_ent = eudp(k)/mu_p_eudp(k)
1416 2568019 : f_ent = max( 0.0_r8, min( 1.0_r8, f_ent ) )
1417 2568019 : tmpa = 1.0_r8 - f_ent
1418 210577558 : do m = 2, ncnst_extd
1419 210577558 : if (doconvproc_extd(m)) then
1420 97584722 : conu(m,k) = tmpa*conu(m,kp1) + f_ent*const(m,k)
1421 : end if
1422 : end do
1423 :
1424 : ! estimate updraft velocity (wup)
1425 2568019 : if (iconvtype /= 1) then
1426 : ! shallow - wup = (mup in kg/m2/s) / [rhoair * (updraft area)]
1427 0 : wup(k) = (mu_i(kp1) + mu_i(k))*0.5_r8*hund_ovr_g &
1428 0 : / (rhoair_i(k) * (cldfrac_i(k)*0.5_r8))
1429 : else
1430 : ! deep - as in shallow, but assumed constant updraft_area with height zm_areafrac
1431 2568019 : wup(k) = (mu_i(kp1) + mu_i(k))*0.5_r8*hund_ovr_g &
1432 2568019 : / (rhoair_i(k) * zm_areafrac)
1433 : end if
1434 :
1435 : ! compute lagrangian transport time (dt_u) and updraft fractional area (fa_u)
1436 : ! *** these must obey dt_u(k)*mu_p_eudp(k) = dp_i(k)*fa_u(k)
1437 2568019 : dt_u(k) = dz/wup(k)
1438 2568019 : dt_u(k) = min( dt_u(k), dt )
1439 2568019 : fa_u(k) = dt_u(k)*(mu_p_eudp(k)/dp_i(k))
1440 :
1441 :
1442 : ! Now apply transformation and removal changes
1443 : ! Skip levels where icwmr(icol,k) <= clw_cut (= 1.0e-6) to eliminate
1444 : ! occasional very small icwmr values from the ZM module
1445 2568019 : clw_cut = 1.0e-6_r8
1446 :
1447 :
1448 : if (convproc_method_activate <= 1) then
1449 : ! aerosol activation - method 1
1450 : ! skip levels that are completely glaciated (fracice(icol,k) == 1.0)
1451 : ! when kactcnt=1 (first/lowest layer with cloud water) apply
1452 : ! activatation to the entire updraft
1453 : ! when kactcnt>1 apply activatation to the amount entrained at this level
1454 : if ((icwmr(icol,k) > clw_cut) .and. (fracice(icol,k) < 1.0_r8)) then
1455 : kactcnt = kactcnt + 1
1456 :
1457 : idiag_act = idiag_in(icol)
1458 : if ((kactcnt == 1) .or. (f_ent > 0.0_r8)) then
1459 : kactcntb = kactcntb + 1
1460 : if ((kactcntb == 1) .and. (idiag_act > 0)) then
1461 : write(lun,'(/a,i9,2i4)') &
1462 : 'qaku act_conv lchnk,i,jtsub', lchnk, icol, jtsub
1463 : end if
1464 : end if
1465 :
1466 : if (kactcnt == 1) then
1467 : ! diagnostic fields
1468 : ! xx_wcldbase = w at first cloudy layer, estimated from mu and cldfrac
1469 : xx_wcldbase(icol) = (mu_i(kp1) + mu_i(k))*0.5_r8*hund_ovr_g &
1470 : / (rhoair_i(k) * (cldfrac_i(k)*0.5_r8))
1471 : xx_kcldbase(icol) = k
1472 :
1473 : kactfirst = k
1474 : tmpa = 1.0_r8
1475 : call ma_activate_convproc( &
1476 : conu(:,k), dconudt_activa(:,k), conu(:,k), &
1477 : tmpa, dt_u(k), wup(k), &
1478 : t(icol,k), rhoair_i(k), fracice(icol,k), &
1479 : pcnst_extd, lun, idiag_act, &
1480 : lchnk, icol, k, &
1481 : ipass_calc_updraft )
1482 : else if (f_ent > 0.0_r8) then
1483 : ! current layer is above cloud base (=first layer with activation)
1484 : ! only allow activation at k = kactfirst thru kactfirst-(method1_activate_nlayers-1)
1485 : if (k >= kactfirst-(method1_activate_nlayers-1)) then
1486 : call ma_activate_convproc( &
1487 : conu(:,k), dconudt_activa(:,k), const(:,k), &
1488 : f_ent, dt_u(k), wup(k), &
1489 : t(icol,k), rhoair_i(k), fracice(icol,k), &
1490 : pcnst_extd, lun, idiag_act, &
1491 : lchnk, icol, k, &
1492 : ipass_calc_updraft )
1493 : end if
1494 : end if
1495 : ! the following was for cam2 shallow convection (hack),
1496 : ! but is not appropriate for cam5 (uwshcu)
1497 : ! else if ((kactcnt > 0) .and. (iconvtype /= 1)) then
1498 : ! ! for shallow conv, when you move from activation occuring to
1499 : ! ! not occuring, reset kactcnt=0, because the hack scheme can
1500 : ! ! produce multiple "1.5 layer clouds" separated by clear air
1501 : ! kactcnt = 0
1502 : ! end if
1503 : end if ! ((icwmr(icol,k) > clw_cut) .and. (fracice(icol,k) < 1.0)) then
1504 :
1505 : else ! (convproc_method_activate >= 2)
1506 : ! aerosol activation - method 2
1507 : ! skip levels that are completely glaciated (fracice(icol,k) == 1.0)
1508 : ! when kactcnt=1 (first/lowest layer with cloud water)
1509 : ! apply "primary" activatation to the entire updraft
1510 : ! when kactcnt>1
1511 : ! apply secondary activatation to the entire updraft
1512 : ! do this for all levels above cloud base (even if completely glaciated)
1513 : ! (this is something for sensitivity testing)
1514 2568019 : do_act_this_lev = .false.
1515 2568019 : if (kactcnt <= 0) then
1516 581935 : if (icwmr(icol,k) > clw_cut) then
1517 161991 : do_act_this_lev = .true.
1518 161991 : kactcnt = 1
1519 161991 : kactfirst = k
1520 : ! diagnostic fields
1521 : ! xx_wcldbase = w at first cloudy layer, estimated from mu and cldfrac
1522 161991 : xx_wcldbase(icol) = (mu_i(kp1) + mu_i(k))*0.5_r8*hund_ovr_g &
1523 161991 : / (rhoair_i(k) * (cldfrac_i(k)*0.5_r8))
1524 161991 : xx_kcldbase(icol) = k
1525 : end if
1526 : else
1527 : ! if ((icwmr(icol,k) > clw_cut) .and. (fracice(icol,k) < 1.0)) then
1528 1986084 : do_act_this_lev = .true.
1529 1986084 : kactcnt = kactcnt + 1
1530 : ! end if
1531 : end if
1532 :
1533 2568019 : idiag_act = idiag_in(icol)
1534 2568019 : if ( do_act_this_lev ) then
1535 2148075 : kactcntb = kactcntb + 1
1536 2148075 : if ((kactcntb == 1) .and. (idiag_act > 0)) then
1537 : write(lun,'(/a,i9,2i4)') &
1538 0 : 'qaku act_conv lchnk,i,jtsub', lchnk, icol, jtsub
1539 : end if
1540 :
1541 : call ma_activate_convproc_method2( &
1542 : conu(:,k), dconudt_activa(:,k), &
1543 2148075 : f_ent, dt_u(k), wup(k), &
1544 2148075 : t(icol,k), rhoair_i(k), fracice(icol,k), &
1545 : pcnst_extd, lun, idiag_act, &
1546 : lchnk, icol, k, &
1547 6444225 : kactfirst, ipass_calc_updraft )
1548 : end if
1549 :
1550 : end if ! (convproc_method_activate <= 1)
1551 :
1552 : ! aqueous chemistry
1553 : ! do glaciated levels as aqchem_conv will eventually do acid vapor uptake
1554 : ! to ice, and aqchem_conv module checks fracice before doing liquid wtr stuff
1555 2568019 : if (icwmr(icol,k) > clw_cut) then
1556 : ! call aqchem_conv( conu(1,k), dconudt_aqchem(1,k), aqfrac, &
1557 : ! t(icol,k), fracice(icol,k), icwmr(icol,k), rhoair_i(k), &
1558 : ! lh2o2(icol,k), lo3(icol,k), dt_u(k) )
1559 : end if
1560 :
1561 : ! wet removal
1562 : !
1563 : ! mirage2
1564 : ! rprd = precip formation as a grid-cell average (kgW/kgA/s)
1565 : ! icwmr = cloud water MR within updraft area (kgW/kgA)
1566 : ! fupdr = updraft fractional area (--)
1567 : ! A = rprd/fupdr = precip formation rate within updraft area (kgW/kgA/s)
1568 : ! B = A/icwmr = rprd/(icwmr*fupdr)
1569 : ! = first-order removal rate (1/s)
1570 : ! C = dp/(mup/fupdr) = updraft air residence time in the layer (s)
1571 : !
1572 : ! fraction removed = (1.0 - exp(-cdt)) where
1573 : ! cdt = B*C = (dp/mup)*rprd/icwmr
1574 : !
1575 : ! Note1: fupdr cancels out in cdt, so need not be specified
1576 : ! Note2: dp & mup units need only be consistent (e.g., mb & mb/s)
1577 : ! Note3: for shallow conv, cdt = 1-beta (beta defined in Hack scheme)
1578 : ! Note4: the "dp" in C above and code below should be the moist dp
1579 : !
1580 : ! cam5
1581 : ! clw_preloss = cloud water MR before loss to precip
1582 : ! = icwmr + dt*(rprd/fupdr)
1583 : ! B = A/clw_preloss = (rprd/fupdr)/(icwmr + dt*rprd/fupdr)
1584 : ! = rprd/(fupdr*icwmr + dt*rprd)
1585 : ! = first-order removal rate (1/s)
1586 : !
1587 : ! fraction removed = (1.0 - exp(-cdt)) where
1588 : ! cdt = B*C = (fupdr*dp/mup)*[rprd/(fupdr*icwmr + dt*rprd)]
1589 : !
1590 : ! Note1: *** cdt is now sensitive to fupdr, which we do not really know,
1591 : ! and is not the same as the convective cloud fraction
1592 : ! Note2: dt is appropriate in the above cdt expression, not dtsub
1593 : !
1594 : ! Apply wet removal at levels where
1595 : ! icwmr(icol,k) > clw_cut AND rprd(icol,k) > 0.0
1596 : ! as wet removal occurs in both liquid and ice clouds
1597 : !
1598 2568019 : cdt(k) = 0.0_r8
1599 2568019 : if ((icwmr(icol,k) > clw_cut) .and. (rprd(icol,k) > 0.0_r8)) then
1600 : ! if (iconvtype == 1) then
1601 1963222 : tmpf = 0.5_r8*cldfrac_i(k)
1602 1963222 : cdt(k) = (tmpf*dp(i,k)/mu_p_eudp(k)) * rprd(icol,k) / &
1603 3926444 : (tmpf*icwmr(icol,k) + dt*rprd(icol,k))
1604 : ! else if (k < pver) then
1605 : ! if (eudp(k+1) > 0) cdt(k) = &
1606 : ! rprd(icol,k)*dp(i,k)/(icwmr(icol,k)*eudp(k+1))
1607 : ! end if
1608 : end if
1609 2568019 : if (cdt(k) > 0.0_r8) then
1610 1963222 : expcdtm1 = exp(-cdt(k)) - 1.0_r8
1611 160984204 : do m = 2, ncnst_extd
1612 160984204 : if (doconvproc_extd(m)) then
1613 74602436 : dconudt_wetdep(m,k) = conu(m,k)*aqfrac(m)*expcdtm1
1614 74602436 : conu(m,k) = conu(m,k) + dconudt_wetdep(m,k)
1615 74602436 : dconudt_wetdep(m,k) = dconudt_wetdep(m,k) / dt_u(k)
1616 : end if
1617 : enddo
1618 : end if
1619 :
1620 : end if ! "(mu_p_eudp(k) > mbsth)"
1621 : end do k_loop_main_bb ! "k = kbot, ktop, -1"
1622 :
1623 : ! when doing updraft calcs twice, only need to go this far on the first pass
1624 : if ( (ipass_calc_updraft == 1) .and. &
1625 : (npass_calc_updraft == 2) ) cycle ipass_calc_updraft_loop
1626 :
1627 167511 : if (idiag_in(icol) > 0) then
1628 : ! do wet removal diagnostics here
1629 0 : do k = kbot, ktop, -1
1630 0 : if (mu_p_eudp(k) > mbsth) &
1631 : write(lun,'(a,i9,3i4,1p,6e10.3)') &
1632 0 : 'qakr - l,i,k,jt; cdt, cldfrac, icwmr, rprd, ...', lchnk, icol, k, jtsub, &
1633 0 : cdt(k), cldfrac_i(k), icwmr(icol,k), rprd(icol,k), dp(i,k), mu_p_eudp(k)
1634 : end do
1635 : end if
1636 :
1637 :
1638 : ! Compute downdraft mixing ratios from cloudtop to cloudbase
1639 : ! No special treatment is needed at k=2
1640 : ! No transformation or removal is applied in the downdraft
1641 2741894 : do k = ktop, kbot
1642 2574383 : kp1 = k + 1
1643 : ! md_m_eddp = downdraft massflux at kp1, without detrainment between k,kp1
1644 2574383 : md_m_eddp = md_i(k) - eddp(k)
1645 2741894 : if (md_m_eddp < -mbsth) then
1646 174562748 : do m = 2, ncnst_extd
1647 174562748 : if (doconvproc_extd(m)) then
1648 80894932 : cond(m,kp1) = ( md_i(k)*cond(m,k) &
1649 161789864 : - eddp(k)*const(m,k) ) / md_m_eddp
1650 : endif
1651 : end do
1652 : end if
1653 : end do ! k
1654 :
1655 :
1656 : ! Now computes fluxes and tendencies
1657 : ! NOTE: The approach used in convtran applies to inert tracers and
1658 : ! must be modified to include source and sink terms
1659 167511 : sumflux(:) = 0.0_r8
1660 167511 : sumflux2(:) = 0.0_r8
1661 167511 : sumsrce(:) = 0.0_r8
1662 167511 : sumchng(:) = 0.0_r8
1663 167511 : sumchng3(:) = 0.0_r8
1664 167511 : sumactiva(:) = 0.0_r8
1665 167511 : sumaqchem(:) = 0.0_r8
1666 167511 : sumwetdep(:) = 0.0_r8
1667 167511 : sumresusp(:) = 0.0_r8
1668 167511 : sumprevap(:) = 0.0_r8
1669 :
1670 167511 : maxflux(:) = 0.0_r8
1671 167511 : maxflux2(:) = 0.0_r8
1672 167511 : maxresusp(:) = 0.0_r8
1673 167511 : maxsrce(:) = 0.0_r8
1674 167511 : maxprevap(:) = 0.0_r8
1675 :
1676 : k_loop_main_cc: &
1677 2741894 : do k = ktop, kbot
1678 2574383 : kp1 = k+1
1679 2574383 : km1 = k-1
1680 2574383 : kp1x = min( kp1, pver )
1681 2574383 : km1x = max( km1, 1 )
1682 2574383 : fa_u_dp = fa_u(k)*dp_i(k)
1683 211266917 : do m = 2, ncnst_extd
1684 211099406 : if (doconvproc_extd(m)) then
1685 :
1686 : ! First compute fluxes using environment subsidence/lifting and
1687 : ! entrainment/detrainment into up/downdrafts,
1688 : ! to provide an additional mass balance check
1689 : ! (this could be deleted after the code is well tested)
1690 195653108 : fluxin = mu_i(k)*min(chat(m,k),const(m,km1x)) &
1691 195653108 : - md_i(kp1)*min(chat(m,kp1),const(m,kp1x)) &
1692 489132770 : + dudp(k)*conu(m,k) + dddp(k)*cond(m,kp1)
1693 : fluxout = mu_i(kp1)*min(chat(m,kp1),const(m,k)) &
1694 : - md_i(k)*min(chat(m,k),const(m,k)) &
1695 97826554 : + (eudp(k) + eddp(k))*const(m,k)
1696 :
1697 97826554 : netflux = fluxin - fluxout
1698 :
1699 97826554 : sumflux2(m) = sumflux2(m) + netflux
1700 97826554 : maxflux2(m) = max( maxflux2(m), abs(fluxin), abs(fluxout) )
1701 :
1702 : ! Now compute fluxes as in convtran, and also source/sink terms
1703 : ! (version 3 limit fluxes outside convection to mass in appropriate layer
1704 : ! (these limiters are probably only safe for positive definite quantitities
1705 : ! (it assumes that mu and md already satify a courant number limit of 1)
1706 97826554 : if (iflux_method /= 2) then
1707 : fluxin = mu_i(kp1)*conu(m,kp1) &
1708 : + mu_i(k )*min(chat(m,k ),const(m,km1x)) &
1709 : - ( md_i(k )*cond(m,k) &
1710 97826554 : + md_i(kp1)*min(chat(m,kp1),const(m,kp1x)) )
1711 : fluxout = mu_i(k )*conu(m,k) &
1712 : + mu_i(kp1)*min(chat(m,kp1),const(m,k )) &
1713 : - ( md_i(kp1)*cond(m,kp1) &
1714 97826554 : + md_i(k )*min(chat(m,k ),const(m,k )) )
1715 : else
1716 : fluxin = mu_i(kp1)*conu(m,kp1) &
1717 0 : - ( md_i(k )*cond(m,k) )
1718 : fluxout = mu_i(k )*conu(m,k) &
1719 0 : - ( md_i(kp1)*cond(m,kp1) )
1720 0 : tmpveca(1) = fluxin ; tmpveca(4) = -fluxout
1721 :
1722 : ! new method -- simple upstream method for the env subsidence
1723 : ! tmpa = net env mass flux (positive up) at top of layer k
1724 0 : tmpa = -( mu_i(k ) + md_i(k ) )
1725 0 : if (tmpa <= 0.0_r8) then
1726 0 : fluxin = fluxin - tmpa*const(m,km1x)
1727 : else
1728 0 : fluxout = fluxout + tmpa*const(m,k )
1729 : end if
1730 0 : tmpveca(2) = fluxin ; tmpveca(5) = -fluxout
1731 : ! tmpa = net env mass flux (positive up) at base of layer k
1732 0 : tmpa = -( mu_i(kp1) + md_i(kp1) )
1733 0 : if (tmpa >= 0.0_r8) then
1734 0 : fluxin = fluxin + tmpa*const(m,kp1x)
1735 : else
1736 0 : fluxout = fluxout - tmpa*const(m,k )
1737 : end if
1738 0 : tmpveca(3) = fluxin ; tmpveca(6) = -fluxout
1739 : end if
1740 :
1741 97826554 : netflux = fluxin - fluxout
1742 : netsrce = fa_u_dp*(dconudt_aqchem(m,k) + &
1743 97826554 : dconudt_activa(m,k) + dconudt_wetdep(m,k))
1744 97826554 : dcondt(m,k) = (netflux+netsrce)/dp_i(k)
1745 :
1746 97826554 : dcondt_wetdep(m,k) = fa_u_dp*dconudt_wetdep(m,k)/dp_i(k)
1747 :
1748 97826554 : sumflux(m) = sumflux(m) + netflux
1749 97826554 : maxflux(m) = max( maxflux(m), abs(fluxin), abs(fluxout) )
1750 97826554 : sumsrce(m) = sumsrce(m) + netsrce
1751 : maxsrce(m) = max( maxsrce(m), &
1752 : fa_u_dp*max( abs(dconudt_aqchem(m,k)), &
1753 97826554 : abs(dconudt_activa(m,k)), abs(dconudt_wetdep(m,k)) ) )
1754 97826554 : sumchng(m) = sumchng(m) + dcondt(m,k)*dp_i(k)
1755 97826554 : sumactiva(m) = sumactiva(m) + fa_u_dp*dconudt_activa(m,k)
1756 97826554 : sumaqchem(m) = sumaqchem(m) + fa_u_dp*dconudt_aqchem(m,k)
1757 97826554 : sumwetdep(m) = sumwetdep(m) + fa_u_dp*dconudt_wetdep(m,k)
1758 :
1759 97826554 : if ( idiag_in(icol)>0 .and. k==26 .and. &
1760 : (m==16 .or. m==23 .or. m==16+pcnst .or. m==23+pcnst) ) then
1761 0 : if (m==16) &
1762 : write(lun,'(a,i9,4i4,1p,22x, 2x,11x, 2x,6e11.3)') &
1763 0 : 'qakww0-'//convtype(1:4), lchnk, icol, k, -1, jtsub, &
1764 0 : dtsub*mu_i(k+1)/dp_i(k), dtsub*mu_i(k)/dp_i(k), dtsub*eudp(k)/dp_i(k), &
1765 0 : dtsub*md_i(k+1)/dp_i(k), dtsub*md_i(k)/dp_i(k), dtsub*eddp(k)/dp_i(k)
1766 :
1767 : write(lun,'(a,i9,4i4,1p,2e11.3,2x,e11.3,2x,6e11.3)') &
1768 0 : 'qakww1-'//convtype(1:4), lchnk, icol, k, m, jtsub, &
1769 0 : const(m,k), const(m,k)+dtsub*dcondt(m,k), dtsub*dcondt(m,k), &
1770 0 : dtsub*fluxin/dp_i(k), -dtsub*fluxout/dp_i(k), &
1771 0 : dtsub*fa_u_dp*dconudt_aqchem(m,k)/dp_i(k), &
1772 0 : dtsub*fa_u_dp*dconudt_activa(m,k)/dp_i(k), &
1773 0 : dtsub*fa_u_dp*dconudt_wetdep(m,k)/dp_i(k)
1774 : write(lun,'(a,i9,4i4,1p,22x, 2x,11x, 2x,6e11.3)') &
1775 0 : 'qakww1-'//convtype(1:4), lchnk, icol, k, m, jtsub, &
1776 0 : dtsub*tmpveca(1:6)/dp_i(k)
1777 : end if
1778 :
1779 : end if ! "(doconvproc_extd(m))"
1780 : end do ! "m = 2,ncnst_extd"
1781 : end do k_loop_main_cc ! "k = ktop, kbot"
1782 :
1783 :
1784 : ! calculate effects of precipitation evaporation
1785 : call ma_precpevap_convproc( dcondt, dcondt_wetdep, dcondt_prevap, &
1786 : rprd, evapc, dp_i, &
1787 : icol, ktop, pcnst_extd, &
1788 167511 : lun, idiag_in(icol), lchnk, &
1789 335022 : doconvproc_extd )
1790 167511 : if ( idiag_in(icol)>0 ) then
1791 0 : k = 26
1792 0 : do m = 16, 23, 7
1793 : write(lun,'(a,i9,4i4,1p,2e11.3,2x,e11.3,2x,5e11.3)') &
1794 0 : 'qakww2-'//convtype(1:4), lchnk, icol, k, m, jtsub, &
1795 0 : const(m,k), const(m,k)+dtsub*dcondt(m,k), dtsub*dcondt(m,k)
1796 : end do
1797 0 : do m = 16+pcnst, 23+pcnst, 7
1798 : write(lun,'(a,i9,4i4,1p,2e11.3,2x,e11.3,2x,5e11.3)') &
1799 0 : 'qakww2-'//convtype(1:4), lchnk, icol, k, m, jtsub, &
1800 0 : const(m,k), const(m,k)+dtsub*dcondt(m,k), dtsub*dcondt(m,k)
1801 : end do
1802 : end if
1803 :
1804 :
1805 :
1806 : ! make adjustments to dcondt for activated & unactivated aerosol species
1807 : ! pairs to account any (or total) resuspension of convective-cloudborne aerosol
1808 : call ma_resuspend_convproc( dcondt, dcondt_resusp, &
1809 167511 : const, dp_i, ktop, kbot_prevap, pcnst_extd )
1810 :
1811 : ! Do resuspension of aerosols from rain only when the rain has
1812 : ! totally evaporated.
1813 167511 : if (convproc_do_evaprain_atonce) then
1814 654465477 : dcondt_resusp3d(pcnst+1:pcnst_extd,icol,:) = dcondt_resusp(pcnst+1:pcnst_extd,:)
1815 654465477 : dcondt_resusp(pcnst+1:pcnst_extd,:) = 0._r8
1816 : end if
1817 :
1818 167511 : if ( idiag_in(icol)>0 ) then
1819 0 : k = 26
1820 0 : do m = 16, 23, 7
1821 : write(lun,'(a,i9,4i4,1p,2e11.3,2x,e11.3,2x,5e11.3)') &
1822 0 : 'qakww3-'//convtype(1:4), lchnk, icol, k, m, jtsub, &
1823 0 : const(m,k), const(m,k)+dtsub*dcondt(m,k), dtsub*dcondt(m,k)
1824 : end do
1825 0 : do m = 16+pcnst, 23+pcnst, 7
1826 : write(lun,'(a,i9,4i4,1p,2e11.3,2x,e11.3,2x,5e11.3)') &
1827 0 : 'qakww3-'//convtype(1:4), lchnk, icol, k, m, jtsub, &
1828 0 : const(m,k), const(m,k)+dtsub*dcondt(m,k), dtsub*dcondt(m,k)
1829 : end do
1830 : end if
1831 :
1832 :
1833 : ! calculate new column-tendency variables
1834 13735902 : do m = 2, ncnst_extd
1835 13735902 : if (doconvproc_extd(m)) then
1836 132731872 : do k = ktop, kbot_prevap
1837 126366454 : sumchng3(m) = sumchng3(m) + dcondt(m,k)*dp_i(k)
1838 126366454 : sumresusp(m) = sumresusp(m) + dcondt_resusp(m,k)*dp_i(k)
1839 : maxresusp(m) = max( maxresusp(m), &
1840 126366454 : abs(dcondt_resusp(m,k)*dp_i(k)) )
1841 126366454 : sumprevap(m) = sumprevap(m) + dcondt_prevap(m,k)*dp_i(k)
1842 : maxprevap(m) = max( maxprevap(m), &
1843 132731872 : abs(dcondt_prevap(m,k)*dp_i(k)) )
1844 : end do
1845 : end if
1846 : end do ! m
1847 :
1848 :
1849 : ! do checks for mass conservation
1850 : ! do not expect errors > 1.0e-14, but use a conservative 1.0e-10 here,
1851 : ! as an error of this size is still not a big concern
1852 167511 : relerr_cut = 1.0e-10_r8
1853 167511 : if (nerr < nerrmax) then
1854 167511 : merr = 0
1855 167511 : if (courantmax > (1.0_r8 + 1.0e-6_r8)) then
1856 0 : write(lun,9161) '-', trim(convtype), courantmax
1857 0 : merr = merr + 1
1858 : end if
1859 13735902 : do m = 2, ncnst_extd
1860 13735902 : if (doconvproc_extd(m)) then
1861 6365418 : itmpa = 0
1862 : ! sumflux should be ~=0.0 because fluxout of one layer cancels
1863 : ! fluxin to adjacent layer
1864 6365418 : tmpa = sumflux(m)
1865 6365418 : tmpb = max( maxflux(m), small )
1866 6365418 : if (abs(tmpa) > relerr_cut*tmpb) then
1867 0 : write(lun,9151) '1', m, cnst_name_extd(m), tmpb, tmpa, (tmpa/tmpb)
1868 0 : itmpa = itmpa + 1
1869 : end if
1870 : ! sumflux2 involve environment fluxes and entrainment/detrainment
1871 : ! to up/downdrafts, and it should be equal to sumchng,
1872 : ! and so (sumflux2 - sumsrce) should be ~=0.0
1873 6365418 : tmpa = sumflux2(m) - sumsrce(m)
1874 6365418 : tmpb = max( maxflux2(m), maxsrce(m), small )
1875 6365418 : if (abs(tmpa) > relerr_cut*tmpb) then
1876 0 : write(lun,9151) '2', m, cnst_name_extd(m), tmpb, tmpa, (tmpa/tmpb)
1877 0 : itmpa = itmpa + 10
1878 : end if
1879 : ! sunchng = sumflux + sumsrce, so (sumchng - sumsrc) should be ~=0.0
1880 6365418 : tmpa = sumchng(m) - sumsrce(m)
1881 6365418 : tmpb = max( maxflux(m), maxsrce(m), small )
1882 6365418 : if (abs(tmpa) > relerr_cut*tmpb) then
1883 0 : write(lun,9151) '3', m, cnst_name_extd(m), tmpb, tmpa, (tmpa/tmpb)
1884 0 : itmpa = itmpa + 100
1885 : end if
1886 : ! sumchng3 = sumchng + sumresusp + sumprevap,
1887 : ! so tmpa (below) should be ~=0.0
1888 : ! NOTE: This check needs to be redone if the rain is being
1889 : ! evaporated all at once. Until then, skip this check for that case.
1890 6365418 : if (.not. convproc_do_evaprain_atonce) then
1891 0 : tmpa = sumchng3(m) - (sumsrce(m) + sumresusp(m) + sumprevap(m))
1892 0 : tmpb = max( maxflux(m), maxsrce(m), maxresusp(m), maxprevap(m), small )
1893 :
1894 0 : if (abs(tmpa) > relerr_cut*tmpb) then
1895 0 : write(lun,9151) '4', m, cnst_name_extd(m), tmpb, tmpa, (tmpa/tmpb)
1896 0 : itmpa = itmpa + 1000
1897 : end if
1898 : end if
1899 :
1900 6365418 : if (itmpa > 0) merr = merr + 1
1901 : end if
1902 : end do ! m
1903 167511 : if (merr > 0) write(lun,9181) convtype, lchnk, icol, i, jt(i), mx(i)
1904 167511 : nerr = nerr + merr
1905 167511 : if (nerr >= nerrmax) write(lun,9171) nerr
1906 : end if ! (nerr < nerrmax) then
1907 :
1908 : 9151 format( '*** ma_convproc_tend error, massbal', a, 1x, i5,1x,a, &
1909 : ' -- maxflux, sumflux, relerr =', 3(1pe14.6) )
1910 : 9161 format( '*** ma_convproc_tend error, courantmax', 2a, 3x, 1pe14.6 )
1911 : 9171 format( '*** ma_convproc_tend error, stopping messages after nerr =', i10 )
1912 :
1913 : 9181 format( '*** ma_convproc_tend error -- convtype, lchnk, icol, il, jt, mx = ', a,2x,5(1x,i10) )
1914 :
1915 :
1916 : !
1917 : ! note again the ma_convproc_tend does not apply convective cloud processing
1918 : ! to the stratiform-cloudborne aerosol
1919 : ! within this routine, cloudborne aerosols are convective-cloudborne
1920 : !
1921 : ! before tendencies (dcondt, which is loaded into dqdt) are returned,
1922 : ! the convective-cloudborne aerosol tendencies must be combined
1923 : ! with the interstitial tendencies
1924 : ! ma_resuspend_convproc has already done this for the dcondt
1925 : !
1926 : ! the individual process column tendencies (sumwetdep, sumprevap, ...)
1927 : ! are just diagnostic fields that can be written to history
1928 : ! tendencies for interstitial and convective-cloudborne aerosol could
1929 : ! both be passed back and output, if desired
1930 : ! currently, however, the interstitial and convective-cloudborne tendencies
1931 : ! are combined (in the next code block) before being passed back (in qsrflx)
1932 : !
1933 837555 : do n = 1, ntot_amode
1934 4020264 : do ll = 0, nspec_amode(n)
1935 3182709 : if (ll == 0) then
1936 670044 : la = numptr_amode(n)
1937 670044 : lc = numptrcw_amode(n) + pcnst
1938 : else
1939 2512665 : la = lmassptr_amode(ll,n)
1940 2512665 : lc = lmassptrcw_amode(ll,n) + pcnst
1941 : end if
1942 3852753 : if (doconvproc(la)) then
1943 3182709 : sumactiva(la) = sumactiva(la) + sumactiva(lc)
1944 3182709 : sumresusp(la) = sumresusp(la) + sumresusp(lc)
1945 3182709 : sumaqchem(la) = sumaqchem(la) + sumaqchem(lc)
1946 3182709 : sumwetdep(la) = sumwetdep(la) + sumwetdep(lc)
1947 3182709 : sumprevap(la) = sumprevap(la) + sumprevap(lc)
1948 : ! if (n==1 .and. ll==1) then
1949 : ! write(lun,*) 'la, sumaqchem(la) =', la, sumaqchem(la)
1950 : ! endif
1951 : end if
1952 : enddo ! ll
1953 : enddo ! n
1954 :
1955 : !
1956 : ! scatter overall tendency back to full array
1957 : !
1958 6867951 : do m = 2, ncnst
1959 6867951 : if (doconvproc(m)) then
1960 66365936 : do k = ktop, kbot_prevap
1961 63183227 : dqdt_i(k,m) = dcondt(m,k)
1962 66365936 : dqdt(icol,k,m) = dqdt(icol,k,m) + dqdt_i(k,m)*xinv_ntsub
1963 : end do
1964 : ! dqdt_i(:,m) = 0.
1965 : end if
1966 : end do ! m
1967 :
1968 : ! scatter column burden tendencies for various processes to qsrflx
1969 6867951 : do m = 2, ncnst
1970 6867951 : if (doconvproc(m)) then
1971 3182709 : qsrflx_i(m,1) = sumactiva(m)*hund_ovr_g
1972 3182709 : qsrflx_i(m,2) = sumresusp(m)*hund_ovr_g
1973 3182709 : qsrflx_i(m,3) = sumaqchem(m)*hund_ovr_g
1974 3182709 : qsrflx_i(m,4) = sumwetdep(m)*hund_ovr_g
1975 3182709 : qsrflx_i(m,5) = sumprevap(m)*hund_ovr_g
1976 : ! qsrflx_i(m,1:4) = 0.
1977 19096254 : qsrflx(icol,m,1:5) = qsrflx(icol,m,1:5) + qsrflx_i(m,1:5)*xinv_ntsub
1978 : end if
1979 : end do ! m
1980 :
1981 :
1982 : ! diagnostic output of profiles before
1983 167511 : if (idiag_in(icol) > 0) then
1984 0 : write(lun, '(/3a,i9,2i4)' ) 'qakr-', trim(convtype), ' - lchnk,i,jtsub', lchnk, icol, jtsub
1985 0 : n = 1
1986 :
1987 0 : do j = 1, 2
1988 0 : if (j == 1) then
1989 : write(lun, '(4a,i4)' ) &
1990 0 : 'qakr-', trim(convtype), ' - k, mu,md; then mode-1 ', &
1991 0 : 'numb & numbcw for q, const, conu, cond, delq(a/c/ac noresu)', jtsub
1992 : else
1993 : write(lun, '(/4a,i4)' ) &
1994 0 : 'qakr-', trim(convtype), ' - k, mu,md; then mode-1 ', &
1995 0 : 'mass & masscw for q, const, conu, cond, delq(a/c/ac noresu)', jtsub
1996 : end if
1997 :
1998 0 : do k = 10, pver
1999 0 : tmpveca(:) = 0.0_r8
2000 0 : do ll = 1, nspec_amode(n)
2001 0 : if (j == 1) then
2002 0 : la = numptr_amode(n)
2003 0 : lc = numptr_amode(n) + pcnst
2004 : else
2005 0 : la = lmassptr_amode(ll,n)
2006 0 : lc = lmassptr_amode(ll,n) + pcnst
2007 : end if
2008 0 : tmpveca(1) = tmpveca(1) + q_i(k,la)
2009 0 : tmpveca(2) = tmpveca(2) + const(la,k)
2010 0 : tmpveca(3) = tmpveca(3) + const(lc,k)
2011 0 : tmpveca(4) = tmpveca(4) + conu( la,k)
2012 0 : tmpveca(5) = tmpveca(5) + conu( lc,k)
2013 0 : tmpveca(6) = tmpveca(6) + cond( la,k)
2014 0 : tmpveca(7) = tmpveca(7) + cond( lc,k)
2015 0 : tmpveca(8) = tmpveca(8) + (dcondt(la,k)-dcondt_resusp(la,k))*dtsub
2016 0 : tmpveca(9) = tmpveca(9) + (dcondt(lc,k)-dcondt_resusp(lc,k))*dtsub
2017 0 : tmpveca(10) = tmpveca(8) + tmpveca(9)
2018 0 : if (j == 1) exit
2019 : end do ! ll
2020 0 : if ((k > 15) .and. (mod(k,5) == 1)) write(lun,'(a)')
2021 0 : write(lun, '(a,i3,1p,2e10.2, e11.2, 3(2x,2e9.2), 2x,3e10.2 )' ) 'qakr', k, &
2022 0 : mu_i(k), md_i(k), tmpveca(1:10)
2023 : end do ! k
2024 : end do ! j
2025 :
2026 : if (pcnst < 0) then
2027 : write(lun, '(/a,i4)' ) &
2028 : 'qakr - name; burden; qsrflx tot, activa,resusp,aqchem,wetdep,resid', jtsub
2029 : do m = 2, ncnst
2030 : if ( .not. doconvproc(m) ) cycle
2031 : tmpveca(1) = sum( q_i(:,m)*dp_i(:) ) * hund_ovr_g
2032 : tmpveca(2) = sum( dqdt_i(:,m)*dp_i(:) ) * hund_ovr_g
2033 : tmpveca(3:6) = qsrflx_i(m,1:4)
2034 : tmpveca(7) = tmpveca(2) - sum( tmpveca(3:6) )
2035 : write(lun, '(2a,1p,2(2x,e11.3),2x,4e11.3,2x,e11.3)' ) &
2036 : 'qakr ', cnst_name_extd(m)(1:10), tmpveca(1:7)
2037 : end do ! m
2038 : end if ! (pcnst < 0) then
2039 :
2040 0 : write(lun, '(/3a,i4)' ) 'qakr-', trim(convtype), &
2041 0 : ' - name; burden; sumchng3, sumactiva,resusp,aqchem,wetdep, resid,resid*dt/burden', jtsub
2042 : ! write(lun, '(/2a)' ) &
2043 : ! 'qakr - name; burden; sumchng3; ', &
2044 : ! 'sumactiva,resusp,aqchem,wetdep,prevap; resid,resid*dtsub/burden'
2045 0 : tmpb = 0.0_r8
2046 0 : itmpb = 0
2047 0 : do m = 2, pcnst
2048 0 : if ( .not. doconvproc_extd(m) ) cycle
2049 :
2050 0 : tmpmata(:,:) = 0.0_r8
2051 0 : do j = 1, 3
2052 0 : l = m
2053 0 : if (j == 3) l = m + pcnst
2054 0 : if ( .not. doconvproc_extd(l) ) cycle
2055 :
2056 0 : if (j == 1) then
2057 0 : tmpmata(1,j) = sum( q_i(:,l)*dp_i(:) ) * hund_ovr_g
2058 0 : tmpmata(2,j) = sum( dqdt_i(:,l)*dp_i(:) ) * hund_ovr_g
2059 0 : tmpmata(3:7,j) = qsrflx_i(l,1:5)
2060 : else
2061 0 : tmpmata(1,j) = sum( const(l,1:pver)*dp_i(1:pver) ) * hund_ovr_g
2062 0 : tmpmata(2,j) = sumchng3( l) * hund_ovr_g
2063 0 : tmpmata(3,j) = sumactiva(l) * hund_ovr_g
2064 0 : tmpmata(4,j) = sumresusp(l) * hund_ovr_g
2065 0 : tmpmata(5,j) = sumaqchem(l) * hund_ovr_g
2066 0 : tmpmata(6,j) = sumwetdep(l) * hund_ovr_g
2067 0 : tmpmata(7,j) = sumprevap(l) * hund_ovr_g
2068 : end if
2069 : end do ! j
2070 :
2071 0 : tmpmata(3:7,2) = tmpmata(3:7,2) - tmpmata(3:7,3) ! because lc values were added onto la
2072 0 : do j = 1, 3
2073 0 : tmpmata(8,j) = tmpmata(2,j) - sum( tmpmata(3:7,j) ) ! residual
2074 0 : tmpa = max( tmpmata(1,min(j,2)), 1.0e-20_r8 )
2075 0 : tmpmata(9,j) = tmpmata(8,j) * dtsub / tmpa
2076 0 : if (abs(tmpmata(9,j)) > tmpb) then
2077 0 : tmpb = abs(tmpmata(9,j))
2078 0 : itmpb = m
2079 : end if
2080 : end do
2081 :
2082 : ! write(lun, '(/2a,1p,2(2x,e11.3),2x,4e11.3,2x,2e11.3)' ) &
2083 : ! 'qakr1 ', cnst_name_extd(m)(1:10), tmpmata(1:6,1), tmpmata(8:9,1)
2084 : write(lun, '(/2a,1p,2(2x,e11.3),2x,5e11.3,2x,2e11.3)' ) &
2085 0 : 'qakr1 ', cnst_name_extd(m)(1:10), tmpmata(1:9,1)
2086 : ! write(lun, '( 2a,1p,2(2x,e11.3),2x,4e11.3,2x,2e11.3)' ) &
2087 : ! 'qakr2 ', cnst_name_extd(m)(1:10), tmpmata(1:6,2), tmpmata(8:9,2)
2088 : write(lun, '( 2a,1p,2(2x,e11.3),2x,5e11.3,2x,2e11.3)' ) &
2089 0 : 'qakr2 ', cnst_name_extd(m)(1:10), tmpmata(1:9,2)
2090 0 : if ( .not. doconvproc_extd(l) ) cycle
2091 : ! write(lun, '( 2a,1p,2(2x,e11.3),2x,4e11.3,2x,2e11.3)' ) &
2092 : ! 'qakr3 ', cnst_name_cw(m)(1:10), tmpmata(1:6,3), tmpmata(8:9,3)
2093 : write(lun, '( 2a,1p,2(2x,e11.3),2x,5e11.3,2x,2e11.3)' ) &
2094 0 : 'qakr3 ', cnst_name_cw(m)(1:10), tmpmata(1:9,3)
2095 : end do ! m
2096 0 : write(lun, '(/3a,2i4,1p,e11.2)' ) 'qakr-', trim(convtype), &
2097 0 : ' - max(resid*dt/burden)', jtsub, itmpb, tmpb
2098 :
2099 : end if ! (idiag_in(icol) > 0) then
2100 :
2101 :
2102 335022 : if (jtsub < ntsub) then
2103 : ! update the q_i for the next interation of the jtsub loop
2104 108117 : do m = 2, ncnst
2105 108117 : if (doconvproc(m)) then
2106 1034721 : do k = ktop, kbot_prevap
2107 1034721 : q_i(k,m) = max( (q_i(k,m) + dqdt_i(k,m)*dtsub), 0.0_r8 )
2108 : end do
2109 : end if
2110 : end do ! m
2111 : end if
2112 :
2113 : end do ipass_calc_updraft_loop
2114 :
2115 : end do jtsub_loop_main_aa ! of the main "do jtsub = 1, ntsub" loop
2116 :
2117 :
2118 : end do i_loop_main_aa ! of the main "do i = il1g, il2g" loop
2119 :
2120 58824 : return
2121 58824 : end subroutine ma_convproc_tend
2122 :
2123 :
2124 :
2125 : !=========================================================================================
2126 167511 : subroutine ma_precpevap_convproc( &
2127 167511 : dcondt, dcondt_wetdep, dcondt_prevap, &
2128 : rprd, evapc, dp_i, &
2129 : icol, ktop, pcnst_extd, &
2130 : lun, idiag_prevap, lchnk, &
2131 167511 : doconvproc_extd )
2132 : !-----------------------------------------------------------------------
2133 : !
2134 : ! Purpose:
2135 : ! Calculate resuspension of wet-removed aerosol species resulting
2136 : ! precip evaporation
2137 : !
2138 : ! Author: R. Easter
2139 : !
2140 : !-----------------------------------------------------------------------
2141 :
2142 : use modal_aero_data, only: &
2143 58824 : lmassptrcw_amode, nspec_amode, numptrcw_amode
2144 :
2145 : implicit none
2146 :
2147 : !-----------------------------------------------------------------------
2148 : ! arguments
2149 : ! (note: TMR = tracer mixing ratio)
2150 : integer, intent(in) :: pcnst_extd
2151 :
2152 : real(r8), intent(inout) :: dcondt(pcnst_extd,pver)
2153 : ! overall TMR tendency from convection
2154 : real(r8), intent(in) :: dcondt_wetdep(pcnst_extd,pver)
2155 : ! portion of TMR tendency due to wet removal
2156 : real(r8), intent(inout) :: dcondt_prevap(pcnst_extd,pver)
2157 : ! portion of TMR tendency due to precip evaporation
2158 : ! (actually, due to the adjustments made here)
2159 : ! (on entry, this is 0.0)
2160 :
2161 : real(r8), intent(in) :: rprd(pcols,pver) ! conv precip production rate (gathered)
2162 : real(r8), intent(in) :: evapc(pcols,pver) ! conv precip evaporation rate (gathered)
2163 : real(r8), intent(in) :: dp_i(pver) ! pressure thickness of level (in mb)
2164 :
2165 : integer, intent(in) :: icol ! normal (ungathered) i index for current column
2166 : integer, intent(in) :: ktop ! index of top cloud level for current column
2167 : integer, intent(in) :: lun ! logical unit for diagnostic output
2168 : integer, intent(in) :: idiag_prevap ! flag for diagnostic output
2169 : integer, intent(in) :: lchnk ! chunk index
2170 :
2171 : logical, intent(in) :: doconvproc_extd(pcnst_extd) ! indicates which species to process
2172 :
2173 : !-----------------------------------------------------------------------
2174 : ! local variables
2175 : integer :: k, l, ll, m, n
2176 : real(r8) :: del_pr_flux_prod ! change to precip flux from production [(kg/kg/s)*mb]
2177 : real(r8) :: del_pr_flux_evap ! change to precip flux from evaporation [(kg/kg/s)*mb]
2178 : real(r8) :: del_wd_flux_evap ! change to wet deposition flux from evaporation [(kg/kg/s)*mb]
2179 : real(r8) :: fdel_pr_flux_evap ! fractional change to precip flux from evaporation
2180 : real(r8) :: pr_flux ! precip flux at base of current layer [(kg/kg/s)*mb]
2181 : real(r8) :: pr_flux_old
2182 : real(r8) :: tmpa, tmpb, tmpc, tmpd
2183 : real(r8) :: tmpdp ! delta-pressure (mb)
2184 335022 : real(r8) :: wd_flux(pcnst_extd) ! tracer wet deposition flux at base of current layer [(kg/kg/s)*mb]
2185 : integer :: i
2186 : character(len=4) :: spcstr
2187 : !-----------------------------------------------------------------------
2188 :
2189 :
2190 167511 : pr_flux = 0.0_r8
2191 13903413 : wd_flux(:) = 0.0_r8
2192 :
2193 167511 : if (idiag_prevap > 0) then
2194 0 : write(lun,'(a,i9,i4,5x,a)') 'qakx - lchnk,i', lchnk, icol, &
2195 0 : '// k; pr_flux old,new; delprod,devap; mode-1 numb wetdep,prevap; mass ...'
2196 : end if
2197 :
2198 3492944 : do k = ktop, pver
2199 3325433 : tmpdp = dp_i(k)
2200 :
2201 3325433 : pr_flux_old = pr_flux
2202 3325433 : del_pr_flux_prod = tmpdp*max(0.0_r8, rprd(icol,k))
2203 3325433 : pr_flux = pr_flux_old + del_pr_flux_prod
2204 :
2205 3325433 : del_pr_flux_evap = min( pr_flux, tmpdp*max(0.0_r8, evapc(icol,k)) )
2206 :
2207 : ! Do resuspension of aerosols from rain only when the rain has
2208 : ! totally evaporated in one layer.
2209 3325433 : if (convproc_do_evaprain_atonce .and. &
2210 3141096 : (del_pr_flux_evap.ne.pr_flux)) del_pr_flux_evap = 0._r8
2211 :
2212 3325433 : fdel_pr_flux_evap = del_pr_flux_evap / max(pr_flux, 1.0e-35_r8)
2213 :
2214 272685506 : do m = 2, pcnst_extd
2215 272685506 : if ( doconvproc_extd(m) ) then
2216 : ! use -dcondt_wetdep(m,k) as it is negative (or zero)
2217 126366454 : wd_flux(m) = wd_flux(m) + tmpdp*max(0.0_r8, -dcondt_wetdep(m,k))
2218 126366454 : del_wd_flux_evap = wd_flux(m)*fdel_pr_flux_evap
2219 126366454 : wd_flux(m) = max( 0.0_r8, wd_flux(m)-del_wd_flux_evap )
2220 :
2221 126366454 : dcondt_prevap(m,k) = del_wd_flux_evap/tmpdp
2222 126366454 : dcondt(m,k) = dcondt(m,k) + dcondt_prevap(m,k)
2223 : end if
2224 : end do
2225 :
2226 : ! Do resuspension of aerosol species from rain to coarse mode (large particle) rather
2227 : ! than to individual modes.
2228 3325433 : if (convproc_do_evaprain_atonce) then
2229 :
2230 3325433 : call accumulate_to_larger_mode( 'SO4', lptr_so4_a_amode, dcondt_prevap(:,k) )
2231 3325433 : call accumulate_to_larger_mode( 'DUST',lptr_dust_a_amode,dcondt_prevap(:,k) )
2232 3325433 : call accumulate_to_larger_mode( 'NACL',lptr_nacl_a_amode,dcondt_prevap(:,k) )
2233 3325433 : call accumulate_to_larger_mode( 'MSA', lptr_msa_a_amode, dcondt_prevap(:,k) )
2234 3325433 : call accumulate_to_larger_mode( 'NH4', lptr_nh4_a_amode, dcondt_prevap(:,k) )
2235 3325433 : call accumulate_to_larger_mode( 'NO3', lptr_no3_a_amode, dcondt_prevap(:,k) )
2236 :
2237 3325433 : spcstr = ' '
2238 6650866 : do i = 1,nsoa
2239 3325433 : if (nsoa>1) write(spcstr,'(i4)') i
2240 6650866 : call accumulate_to_larger_mode( 'SOA'//adjustl(spcstr), lptr2_soa_a_amode(:,i), dcondt_prevap(:,k) )
2241 : enddo
2242 3325433 : spcstr = ' '
2243 6650866 : do i = 1,npoa
2244 3325433 : if (npoa>1) write(spcstr,'(i4)') i
2245 6650866 : call accumulate_to_larger_mode( 'POM'//adjustl(spcstr), lptr2_pom_a_amode(:,i), dcondt_prevap(:,k) )
2246 : enddo
2247 3325433 : spcstr = ' '
2248 6650866 : do i = 1,nbc
2249 3325433 : if (nbc>1) write(spcstr,'(i4)') i
2250 6650866 : call accumulate_to_larger_mode( 'BC'//adjustl(spcstr), lptr2_bc_a_amode(:,i), dcondt_prevap(:,k) )
2251 : enddo
2252 :
2253 : end if
2254 :
2255 3325433 : pr_flux = max( 0.0_r8, pr_flux-del_pr_flux_evap )
2256 :
2257 3492944 : if (idiag_prevap > 0) then
2258 0 : n = 1
2259 0 : l = numptrcw_amode(n) + pcnst
2260 0 : tmpa = dcondt_wetdep(l,k)
2261 0 : tmpb = dcondt_prevap(l,k)
2262 0 : tmpc = 0.0_r8
2263 0 : tmpd = 0.0_r8
2264 0 : do ll = 1, nspec_amode(n)
2265 0 : l = lmassptrcw_amode(ll,n) + pcnst
2266 0 : tmpc = tmpc + dcondt_wetdep(l,k)
2267 0 : tmpd = tmpd + dcondt_prevap(l,k)
2268 : end do
2269 0 : write(lun,'(a,i4,1p,4(2x,2e10.2))') 'qakx', k, &
2270 0 : pr_flux_old, pr_flux, del_pr_flux_prod, -del_pr_flux_evap, &
2271 0 : -tmpa, tmpb, -tmpc, tmpd
2272 : end if
2273 : end do ! k
2274 :
2275 167511 : return
2276 167511 : end subroutine ma_precpevap_convproc
2277 :
2278 : !=========================================================================================
2279 29928897 : subroutine accumulate_to_larger_mode( spc_name, lptr, prevap )
2280 :
2281 : character(len=*), intent(in) :: spc_name
2282 : integer, intent(in) :: lptr(:)
2283 : real(r8), intent(inout) :: prevap(:)
2284 :
2285 : integer :: m,n, nl,ns
2286 :
2287 : ! find constituent index of the largest mode for the species
2288 69834093 : loop1: do m = 1,ntot_amode-1
2289 59857794 : nl = lptr(mode_size_order(m))
2290 69834093 : if (nl>0) exit loop1
2291 : end do loop1
2292 :
2293 29928897 : if (.not. nl>0) return
2294 :
2295 : ! accumulate the smaller modes into the largest mode
2296 69834093 : do n = m+1,ntot_amode
2297 49881495 : ns = lptr(mode_size_order(n))
2298 69834093 : if (ns>0) then
2299 29928897 : prevap(nl) = prevap(nl) + prevap(ns)
2300 29928897 : prevap(ns) = 0._r8
2301 : if (masterproc .and. debug) then
2302 : write(iulog,'(a,i3,a,i3)') trim(spc_name)//' mode number accumulate ',ns,'->',nl
2303 : endif
2304 : endif
2305 : end do
2306 :
2307 167511 : end subroutine accumulate_to_larger_mode
2308 :
2309 : !=========================================================================================
2310 : subroutine ma_activate_convproc( &
2311 : conu, dconudt, conent, &
2312 : f_ent, dt_u, wup, &
2313 : tair, rhoair, fracice, &
2314 : pcnst_extd, lun, idiag_act, &
2315 : lchnk, i, k, &
2316 : ipass_calc_updraft )
2317 : !-----------------------------------------------------------------------
2318 : !
2319 : ! Purpose:
2320 : ! Calculate activation of aerosol species in convective updraft
2321 : ! for a single column and level
2322 : !
2323 : ! Method:
2324 : ! conu(l) = Updraft TMR (tracer mixing ratio) at k/k-1 interface
2325 : ! conent(l) = TMR of air that is entrained into the updraft from level k
2326 : ! f_ent = Fraction of the "before-detrainment" updraft massflux at
2327 : ! k/k-1 interface" resulting from entrainment of level k air
2328 : ! (where k is the current level in subr ma_convproc_tend)
2329 : !
2330 : ! On entry to this routine, the conu(l) represents the updraft TMR
2331 : ! after entrainment, but before chemistry/physics and detrainment,
2332 : ! and is equal to
2333 : ! conu(l) = f_ent*conent(l) + (1.0-f_ent)*conu_below(l)
2334 : ! where
2335 : ! conu_below(l) = updraft TMR at the k+1/k interface, and
2336 : ! f_ent = (eudp/mu_p_eudp) is the fraction of the updraft massflux
2337 : ! from level k entrainment
2338 : !
2339 : ! This routine applies aerosol activation to the entrained tracer,
2340 : ! then adjusts the conu so that on exit,
2341 : ! conu(la) = conu_incoming(la) - f_ent*conent(la)*f_act(la)
2342 : ! conu(lc) = conu_incoming(lc) + f_ent*conent(la)*f_act(la)
2343 : ! where
2344 : ! la, lc = indices for an unactivated/activated aerosol component pair
2345 : ! f_act = fraction of conent(la) that is activated. The f_act are
2346 : ! calculated with the Razzak-Ghan activation parameterization.
2347 : ! The f_act differ for each mode, and for number/surface/mass.
2348 : !
2349 : ! Note: At the lowest layer with cloud water, subr convproc calls this
2350 : ! routine with conent==conu and f_ent==1.0, with the result that
2351 : ! activation is applied to the entire updraft tracer flux
2352 : !
2353 : ! *** The updraft velocity used for activation calculations is rather
2354 : ! uncertain and needs more work. However, an updraft of 1-3 m/s
2355 : ! will activate essentially all of accumulation and coarse mode particles.
2356 : !
2357 : ! Author: R. Easter
2358 : !
2359 : !-----------------------------------------------------------------------
2360 :
2361 : use ndrop, only: activate_aerosol
2362 :
2363 : use modal_aero_data, only: lmassptr_amode, lmassptrcw_amode, &
2364 : ntot_amode, &
2365 : nspec_amode, ntot_amode, numptr_amode, numptrcw_amode, &
2366 : specdens_amode, spechygro, &
2367 : voltonumblo_amode, voltonumbhi_amode
2368 :
2369 : implicit none
2370 :
2371 : !-----------------------------------------------------------------------
2372 : ! arguments (note: TMR = tracer mixing ratio)
2373 : integer, intent(in) :: pcnst_extd
2374 : ! conu = tracer mixing ratios in updraft at top of this (current) level
2375 : ! The conu are changed by activation
2376 : real(r8), intent(inout) :: conu(pcnst_extd)
2377 : ! conent = TMRs in the entrained air at this level
2378 : real(r8), intent(in) :: conent(pcnst_extd)
2379 : real(r8), intent(inout) :: dconudt(pcnst_extd) ! TMR tendencies due to activation
2380 :
2381 : real(r8), intent(in) :: f_ent ! fraction of updraft massflux that was
2382 : ! entrained across this layer == eudp/mu_p_eudp
2383 : real(r8), intent(in) :: dt_u ! lagrangian transport time (s) in the
2384 : ! updraft at current level
2385 : real(r8), intent(in) :: wup ! mean updraft vertical velocity (m/s)
2386 : ! at current level updraft
2387 :
2388 : real(r8), intent(in) :: tair ! Temperature in Kelvin
2389 : real(r8), intent(in) :: rhoair ! air density (kg/m3)
2390 :
2391 : real(r8), intent(in) :: fracice ! Fraction of ice within the cloud
2392 : ! used as in-cloud wet removal rate
2393 : integer, intent(in) :: lun ! logical unit for diagnostic output
2394 : integer, intent(in) :: idiag_act ! flag for diagnostic output
2395 : integer, intent(in) :: lchnk ! chunk index
2396 : integer, intent(in) :: i ! column index
2397 : integer, intent(in) :: k ! level index
2398 : integer, intent(in) :: ipass_calc_updraft
2399 :
2400 : !-----------------------------------------------------------------------
2401 : ! local variables
2402 : integer :: ll, la, lc, n
2403 :
2404 : real(r8) :: delact ! working variable
2405 : real(r8) :: dt_u_inv ! 1.0/dt_u
2406 : real(r8) :: fluxm(ntot_amode) ! to understand this, see subr activate_aerosol
2407 : real(r8) :: fluxn(ntot_amode) ! to understand this, see subr activate_aerosol
2408 : real(r8) :: flux_fullact ! to understand this, see subr activate_aerosol
2409 : real(r8) :: fm(ntot_amode) ! mass fraction of aerosols activated
2410 : real(r8) :: fn(ntot_amode) ! number fraction of aerosols activated
2411 : real(r8) :: hygro(ntot_amode) ! current hygroscopicity for int+act
2412 : real(r8) :: naerosol(ntot_amode) ! interstitial+activated number conc (#/m3)
2413 : real(r8) :: sigw ! standard deviation of updraft velocity (cm/s)
2414 : real(r8) :: tmpa, tmpb, tmpc ! working variable
2415 : real(r8) :: tmp_fact ! working variable
2416 : real(r8) :: vaerosol(ntot_amode) ! int+act volume (m3/m3)
2417 : real(r8) :: wbar ! mean updraft velocity (cm/s)
2418 : real(r8) :: wdiab ! diabatic vertical velocity (cm/s)
2419 : real(r8) :: wminf, wmaxf ! limits for integration over updraft spectrum (cm/s)
2420 :
2421 :
2422 : !-----------------------------------------------------------------------
2423 :
2424 :
2425 : ! when ipass_calc_updraft == 2, apply the activation tendencies
2426 : ! from pass 1, but multiplied by factor_reduce_actfrac
2427 : ! (can only have ipass_calc_updraft == 2 when method_reduce_actfrac = 2)
2428 : if (ipass_calc_updraft == 2) then
2429 :
2430 : dt_u_inv = 1.0_r8/dt_u
2431 : do n = 1, ntot_amode
2432 : do ll = 0, nspec_amode(n)
2433 : if (ll == 0) then
2434 : la = numptr_amode(n)
2435 : lc = numptrcw_amode(n) + pcnst
2436 : else
2437 : la = lmassptr_amode(ll,n)
2438 : lc = lmassptrcw_amode(ll,n) + pcnst
2439 : end if
2440 :
2441 : delact = dconudt(lc)*dt_u * factor_reduce_actfrac
2442 : delact = min( delact, conu(la) )
2443 : delact = max( delact, 0.0_r8 )
2444 : conu(la) = conu(la) - delact
2445 : conu(lc) = conu(lc) + delact
2446 : dconudt(la) = -delact*dt_u_inv
2447 : dconudt(lc) = delact*dt_u_inv
2448 : end do
2449 : end do ! "n = 1, ntot_amode"
2450 : return
2451 :
2452 : end if ! (ipass_calc_updraft == 2)
2453 :
2454 :
2455 : ! check f_ent > 0
2456 : if (f_ent <= 0.0_r8) return
2457 :
2458 :
2459 : do n = 1, ntot_amode
2460 : ! compute a (or a+cw) volume and hygroscopicity
2461 : tmpa = 0.0_r8
2462 : tmpb = 0.0_r8
2463 : do ll = 1, nspec_amode(n)
2464 : tmpc = max( conent(lmassptr_amode(ll,n)), 0.0_r8 )
2465 : if ( use_cwaer_for_activate_maxsat ) &
2466 : tmpc = tmpc + max( conent(lmassptrcw_amode(ll,n)+pcnst), 0.0_r8 )
2467 : tmpc = tmpc / specdens_amode(ll,n)
2468 : tmpa = tmpa + tmpc
2469 : tmpb = tmpb + tmpc * spechygro(ll,n)
2470 : end do
2471 : vaerosol(n) = tmpa * rhoair
2472 : if (tmpa < 1.0e-35_r8) then
2473 : hygro(n) = 0.2_r8
2474 : else
2475 : hygro(n) = tmpb/tmpa
2476 : end if
2477 :
2478 : ! load a (or a+cw) number and bound it
2479 : tmpa = max( conent(numptr_amode(n)), 0.0_r8 )
2480 : if ( use_cwaer_for_activate_maxsat ) &
2481 : tmpa = tmpa + max( conent(numptrcw_amode(n)+pcnst), 0.0_r8 )
2482 : naerosol(n) = tmpa * rhoair
2483 : naerosol(n) = max( naerosol(n), &
2484 : vaerosol(n)*voltonumbhi_amode(n) )
2485 : naerosol(n) = min( naerosol(n), &
2486 : vaerosol(n)*voltonumblo_amode(n) )
2487 :
2488 : ! diagnostic output for testing/development
2489 : ! if (lun > 0) then
2490 : ! if (n == 1) then
2491 : ! write(lun,9500)
2492 : ! write(lun,9510) (cnst_name(l), conu(l), l=1,pcnst_extd)
2493 : ! write(lun,9520) tair, rhoaircgs, airconcgs
2494 : ! end if
2495 : ! write(lun,9530) n, ntype(n), vaerosol
2496 : ! write(lun,9540) naerosol(n), tmp*airconcgs, &
2497 : ! voltonumbhi_amode(n), voltonumblo_amode(n)
2498 : ! write(lun,9550) (maerosol(l,n), l=1,ntype(n))
2499 : !9500 format( / 'activate_conv output -- conu values' )
2500 : !9510 format( 3( a, 1pe11.3, 4x ) )
2501 : !9520 format( 'ta, rhoa, acon ', 3(1pe11.3) )
2502 : !9530 format( 'n, ntype, sg, vol ', i6, i5, 2(1pe11.3) )
2503 : !9540 format( 'num, num0, v2nhi&lo', 4(1pe11.3) )
2504 : !9550 format( 'masses ', 6(1pe11.3) )
2505 : ! end if
2506 :
2507 : end do
2508 :
2509 :
2510 : ! call Razzak-Ghan activation routine with single updraft
2511 : wbar = max( wup, 0.5_r8 ) ! force wbar >= 0.5 m/s for now
2512 : sigw = 0.0_r8
2513 : wdiab = 0.0_r8
2514 : wminf = wbar
2515 : wmaxf = wbar
2516 :
2517 : call activate_aerosol( &
2518 : wbar, sigw, wdiab, wminf, wmaxf, tair, rhoair, &
2519 : naerosol, ntot_amode, vaerosol, hygro, aero_props_obj, &
2520 : fn, fm, fluxn, fluxm, flux_fullact )
2521 :
2522 :
2523 :
2524 : ! diagnostic output for testing/development
2525 : if (idiag_act > 0) then
2526 : n = min( ntot_amode, 3 )
2527 : write(lun, '(a,i3,2f6.3, 1p,2(2x,3e10.2), 0p,3(2x,3f6.3) )' ) &
2528 : 'qaku k,w,qn,qm,hy,fn,fm', k, wup, wbar, &
2529 : naerosol(1:n)/rhoair, vaerosol(1:n)*1.8e3_r8/rhoair, &
2530 : hygro(1:n), fn(1:n), fm(1:n)
2531 : ! convert naer, vaer to number and (approx) mass TMRs
2532 : end if
2533 : ! if (lun > 0) then
2534 : ! write(lun,9560) (fn(n), n=1,ntot_amode)
2535 : ! write(lun,9570) (fm(n), n=1,ntot_amode)
2536 : !9560 format( 'fnact values ', 6(1pe11.3) )
2537 : !9570 format( 'fmact values ', 6(1pe11.3) )
2538 : ! end if
2539 :
2540 :
2541 : ! apply the activation fractions to the updraft aerosol mixing ratios
2542 : dt_u_inv = 1.0_r8/dt_u
2543 :
2544 : do n = 1, ntot_amode
2545 : do ll = 0, nspec_amode(n)
2546 : if (ll == 0) then
2547 : la = numptr_amode(n)
2548 : lc = numptrcw_amode(n) + pcnst
2549 : tmp_fact = fn(n)
2550 : else
2551 : la = lmassptr_amode(ll,n)
2552 : lc = lmassptrcw_amode(ll,n) + pcnst
2553 : tmp_fact = fm(n)
2554 : end if
2555 :
2556 : if ( (method_reduce_actfrac == 1) .and. &
2557 : (factor_reduce_actfrac >= 0.0_r8) .and. &
2558 : (factor_reduce_actfrac < 1.0_r8) ) &
2559 : tmp_fact = tmp_fact * factor_reduce_actfrac
2560 :
2561 : delact = min( conent(la)*tmp_fact*f_ent, conu(la) )
2562 : delact = max( delact, 0.0_r8 )
2563 : conu(la) = conu(la) - delact
2564 : conu(lc) = conu(lc) + delact
2565 : dconudt(la) = -delact*dt_u_inv
2566 : dconudt(lc) = delact*dt_u_inv
2567 : end do
2568 : end do ! "n = 1, ntot_amode"
2569 :
2570 : return
2571 : end subroutine ma_activate_convproc
2572 :
2573 :
2574 :
2575 : !=========================================================================================
2576 2148075 : subroutine ma_activate_convproc_method2( &
2577 2148075 : conu, dconudt, &
2578 : f_ent, dt_u, wup, &
2579 : tair, rhoair, fracice, &
2580 : pcnst_extd, lun, idiag_act, &
2581 : lchnk, i, k, &
2582 : kactfirst, ipass_calc_updraft )
2583 : !-----------------------------------------------------------------------
2584 : !
2585 : ! Purpose:
2586 : ! Calculate activation of aerosol species in convective updraft
2587 : ! for a single column and level
2588 : !
2589 : ! Method:
2590 : ! conu(l) = Updraft TMR (tracer mixing ratio) at k/k-1 interface
2591 : ! f_ent = Fraction of the "before-detrainment" updraft massflux at
2592 : ! k/k-1 interface" resulting from entrainment of level k air
2593 : ! (where k is the current level in subr ma_convproc_tend)
2594 : !
2595 : ! On entry to this routine, the conu(l) represents the updraft TMR
2596 : ! after entrainment, but before chemistry/physics and detrainment.
2597 : !
2598 : ! This routine applies aerosol activation to the conu tracer mixing ratios,
2599 : ! then adjusts the conu so that on exit,
2600 : ! conu(la) = conu_incoming(la) - conu(la)*f_act(la)
2601 : ! conu(lc) = conu_incoming(lc) + conu(la)*f_act(la)
2602 : ! where
2603 : ! la, lc = indices for an unactivated/activated aerosol component pair
2604 : ! f_act = fraction of conu(la) that is activated. The f_act are
2605 : ! calculated with the Razzak-Ghan activation parameterization.
2606 : ! The f_act differ for each mode, and for number/surface/mass.
2607 : !
2608 : ! At cloud base (k==kactfirst), primary activation is done using the
2609 : ! "standard" code in subr activate do diagnose maximum supersaturation.
2610 : ! Above cloud base, secondary activation is done using a
2611 : ! prescribed supersaturation.
2612 : !
2613 : ! *** The updraft velocity used for activation calculations is rather
2614 : ! uncertain and needs more work. However, an updraft of 1-3 m/s
2615 : ! will activate essentially all of accumulation and coarse mode particles.
2616 : !
2617 : ! Author: R. Easter
2618 : !
2619 : !-----------------------------------------------------------------------
2620 :
2621 : use ndrop, only: activate_aerosol
2622 :
2623 : use modal_aero_data, only: lmassptr_amode, lmassptrcw_amode, &
2624 : ntot_amode, &
2625 : nspec_amode, ntot_amode, numptr_amode, numptrcw_amode, &
2626 : specdens_amode, spechygro, &
2627 : voltonumblo_amode, voltonumbhi_amode
2628 :
2629 : use rad_constituents,only: rad_cnst_get_info
2630 :
2631 : implicit none
2632 :
2633 : !-----------------------------------------------------------------------
2634 : ! arguments (note: TMR = tracer mixing ratio)
2635 : integer, intent(in) :: pcnst_extd
2636 : ! conu = tracer mixing ratios in updraft at top of this (current) level
2637 : ! The conu are changed by activation
2638 : real(r8), intent(inout) :: conu(pcnst_extd)
2639 : real(r8), intent(inout) :: dconudt(pcnst_extd) ! TMR tendencies due to activation
2640 :
2641 : real(r8), intent(in) :: f_ent ! fraction of updraft massflux that was
2642 : ! entrained across this layer == eudp/mu_p_eudp
2643 : real(r8), intent(in) :: dt_u ! lagrangian transport time (s) in the
2644 : ! updraft at current level
2645 : real(r8), intent(in) :: wup ! mean updraft vertical velocity (m/s)
2646 : ! at current level updraft
2647 :
2648 : real(r8), intent(in) :: tair ! Temperature in Kelvin
2649 : real(r8), intent(in) :: rhoair ! air density (kg/m3)
2650 :
2651 : real(r8), intent(in) :: fracice ! Fraction of ice within the cloud
2652 : ! used as in-cloud wet removal rate
2653 : integer, intent(in) :: lun ! logical unit for diagnostic output
2654 : integer, intent(in) :: idiag_act ! flag for diagnostic output
2655 : integer, intent(in) :: lchnk ! chunk index
2656 : integer, intent(in) :: i ! column index
2657 : integer, intent(in) :: k ! level index
2658 : integer, intent(in) :: kactfirst ! k at cloud base
2659 : integer, intent(in) :: ipass_calc_updraft
2660 :
2661 : !-----------------------------------------------------------------------
2662 : ! local variables
2663 : integer :: ll, la, lc, n
2664 :
2665 : real(r8) :: delact ! working variable
2666 : real(r8) :: dt_u_inv ! 1.0/dt_u
2667 4296150 : real(r8) :: fluxm(ntot_amode) ! to understand this, see subr activate_aerosol
2668 4296150 : real(r8) :: fluxn(ntot_amode) ! to understand this, see subr activate_aerosol
2669 : real(r8) :: flux_fullact ! to understand this, see subr activate_aerosol
2670 4296150 : real(r8) :: fm(ntot_amode) ! mass fraction of aerosols activated
2671 4296150 : real(r8) :: fn(ntot_amode) ! number fraction of aerosols activated
2672 4296150 : real(r8) :: hygro(ntot_amode) ! current hygroscopicity for int+act
2673 4296150 : real(r8) :: naerosol(ntot_amode) ! interstitial+activated number conc (#/m3)
2674 : real(r8) :: sigw ! standard deviation of updraft velocity (cm/s)
2675 : real(r8) :: smax_prescribed ! prescribed supersaturation for secondary activation (0-1 fraction)
2676 : real(r8) :: tmpa, tmpb, tmpc ! working variable
2677 : real(r8) :: tmp_fact ! working variable
2678 4296150 : real(r8) :: vaerosol(ntot_amode) ! int+act volume (m3/m3)
2679 : real(r8) :: wbar ! mean updraft velocity (cm/s)
2680 : real(r8) :: wdiab ! diabatic vertical velocity (cm/s)
2681 : real(r8) :: wminf, wmaxf ! limits for integration over updraft spectrum (cm/s)
2682 :
2683 : character(len=32) :: spec_type
2684 :
2685 : !-----------------------------------------------------------------------
2686 :
2687 :
2688 : ! when ipass_calc_updraft == 2, apply the activation tendencies
2689 : ! from pass 1, but multiplied by factor_reduce_actfrac
2690 : ! (can only have ipass_calc_updraft == 2 when method_reduce_actfrac = 2)
2691 2148075 : if (ipass_calc_updraft == 2) then
2692 :
2693 0 : dt_u_inv = 1.0_r8/dt_u
2694 0 : do n = 1, ntot_amode
2695 0 : do ll = 0, nspec_amode(n)
2696 0 : if (ll == 0) then
2697 0 : la = numptr_amode(n)
2698 0 : lc = numptrcw_amode(n) + pcnst
2699 : else
2700 0 : la = lmassptr_amode(ll,n)
2701 0 : lc = lmassptrcw_amode(ll,n) + pcnst
2702 : end if
2703 :
2704 0 : delact = dconudt(lc)*dt_u * factor_reduce_actfrac
2705 0 : delact = min( delact, conu(la) )
2706 0 : delact = max( delact, 0.0_r8 )
2707 0 : conu(la) = conu(la) - delact
2708 0 : conu(lc) = conu(lc) + delact
2709 0 : dconudt(la) = -delact*dt_u_inv
2710 0 : dconudt(lc) = delact*dt_u_inv
2711 : end do
2712 : end do ! "n = 1, ntot_amode"
2713 : return
2714 :
2715 : end if ! (ipass_calc_updraft == 2)
2716 :
2717 :
2718 : ! check f_ent > 0
2719 2148075 : if (f_ent <= 0.0_r8) return
2720 :
2721 :
2722 9930420 : do n = 1, ntot_amode
2723 : ! compute a (or a+cw) volume and hygroscopicity
2724 7944336 : tmpa = 0.0_r8
2725 7944336 : tmpb = 0.0_r8
2726 37735596 : do ll = 1, nspec_amode(n)
2727 29791260 : tmpc = max( conu(lmassptr_amode(ll,n)), 0.0_r8 )
2728 : if ( use_cwaer_for_activate_maxsat ) &
2729 : tmpc = tmpc + max( conu(lmassptrcw_amode(ll,n)+pcnst), 0.0_r8 )
2730 29791260 : tmpc = tmpc / specdens_amode(ll,n)
2731 29791260 : tmpa = tmpa + tmpc
2732 :
2733 : ! Change the hygroscopicity of POM based on the discussion with Prof.
2734 : ! Xiaohong Liu. Some observational studies found that the primary organic
2735 : ! material from biomass burning emission shows very high hygroscopicity.
2736 : ! Also, found that BC mass will be overestimated if all the aerosols in
2737 : ! the primary mode are free to be removed. Therefore, set the hygroscopicity
2738 : ! of POM here as 0.2 to enhance the wet scavenge of primary BC and POM.
2739 :
2740 29791260 : call rad_cnst_get_info(0, n, ll, spec_type=spec_type)
2741 37735596 : if (spec_type=='p-organic' .and. convproc_pom_spechygro>0._r8) then
2742 3972168 : tmpb = tmpb + tmpc * convproc_pom_spechygro
2743 : else
2744 25819092 : tmpb = tmpb + tmpc * spechygro(ll,n)
2745 : end if
2746 : end do
2747 7944336 : vaerosol(n) = tmpa * rhoair
2748 7944336 : if (tmpa < 1.0e-35_r8) then
2749 0 : hygro(n) = 0.2_r8
2750 : else
2751 7944336 : hygro(n) = tmpb/tmpa
2752 : end if
2753 :
2754 : ! load a (or a+cw) number and bound it
2755 7944336 : tmpa = max( conu(numptr_amode(n)), 0.0_r8 )
2756 : if ( use_cwaer_for_activate_maxsat ) &
2757 : tmpa = tmpa + max( conu(numptrcw_amode(n)+pcnst), 0.0_r8 )
2758 7944336 : naerosol(n) = tmpa * rhoair
2759 : naerosol(n) = max( naerosol(n), &
2760 7944336 : vaerosol(n)*voltonumbhi_amode(n) )
2761 : naerosol(n) = min( naerosol(n), &
2762 9930420 : vaerosol(n)*voltonumblo_amode(n) )
2763 :
2764 : ! diagnostic output for testing/development
2765 : ! if (lun > 0) then
2766 : ! if (n == 1) then
2767 : ! write(lun,9500)
2768 : ! write(lun,9510) (cnst_name(l), conu(l), l=1,pcnst_extd)
2769 : ! write(lun,9520) tair, rhoaircgs, airconcgs
2770 : ! end if
2771 : ! write(lun,9530) n, ntype(n), vaerosol
2772 : ! write(lun,9540) naerosol(n), tmp*airconcgs, &
2773 : ! voltonumbhi_amode(n), voltonumblo_amode(n)
2774 : ! write(lun,9550) (maerosol(l,n), l=1,ntype(n))
2775 : !9500 format( / 'activate_conv output -- conu values' )
2776 : !9510 format( 3( a, 1pe11.3, 4x ) )
2777 : !9520 format( 'ta, rhoa, acon ', 3(1pe11.3) )
2778 : !9530 format( 'n, ntype, sg, vol ', i6, i5, 2(1pe11.3) )
2779 : !9540 format( 'num, num0, v2nhi&lo', 4(1pe11.3) )
2780 : !9550 format( 'masses ', 6(1pe11.3) )
2781 : ! end if
2782 :
2783 : end do
2784 :
2785 :
2786 : ! call Razzak-Ghan activation routine with single updraft
2787 1986084 : wbar = max( wup, 0.5_r8 ) ! force wbar >= 0.5 m/s for now
2788 1986084 : sigw = 0.0_r8
2789 1986084 : wdiab = 0.0_r8
2790 1986084 : wminf = wbar
2791 1986084 : wmaxf = wbar
2792 :
2793 1986084 : if (k == kactfirst) then
2794 :
2795 : call activate_aerosol( &
2796 : wbar, sigw, wdiab, wminf, wmaxf, tair, rhoair, &
2797 : naerosol, ntot_amode, vaerosol, hygro, aero_props_obj, &
2798 161991 : fn, fm, fluxn, fluxm, flux_fullact )
2799 :
2800 :
2801 : else
2802 : ! above cloud base - do secondary activation with prescribed supersat
2803 : ! that is constant with height
2804 1824093 : smax_prescribed = method2_activate_smaxmax
2805 : call activate_aerosol( &
2806 : wbar, sigw, wdiab, wminf, wmaxf, tair, rhoair, &
2807 : naerosol, ntot_amode, vaerosol, hygro, aero_props_obj, &
2808 1824093 : fn, fm, fluxn, fluxm, flux_fullact, smax_prescribed )
2809 : end if
2810 :
2811 :
2812 : ! diagnostic output for testing/development
2813 1986084 : if (idiag_act > 0) then
2814 0 : n = min( ntot_amode, 3 )
2815 : write(lun, '(a,i3,2f6.3, 1p,2(2x,3e10.2), 0p,3(2x,3f6.3) )' ) &
2816 0 : 'qaku k,w,qn,qm,hy,fn,fm', k, wup, wbar, &
2817 0 : naerosol(1:n)/rhoair, vaerosol(1:n)*1.8e3_r8/rhoair, &
2818 0 : hygro(1:n), fn(1:n), fm(1:n)
2819 : ! convert naer, vaer to number and (approx) mass TMRs
2820 : end if
2821 : ! if (lun > 0) then
2822 : ! write(lun,9560) (fn(n), n=1,ntot_amode)
2823 : ! write(lun,9570) (fm(n), n=1,ntot_amode)
2824 : !9560 format( 'fnact values ', 6(1pe11.3) )
2825 : !9570 format( 'fmact values ', 6(1pe11.3) )
2826 : ! end if
2827 :
2828 :
2829 : ! apply the activation fractions to the updraft aerosol mixing ratios
2830 1986084 : dt_u_inv = 1.0_r8/dt_u
2831 :
2832 9930420 : do n = 1, ntot_amode
2833 47666016 : do ll = 0, nspec_amode(n)
2834 37735596 : if (ll == 0) then
2835 7944336 : la = numptr_amode(n)
2836 7944336 : lc = numptrcw_amode(n) + pcnst
2837 7944336 : tmp_fact = fn(n)
2838 : else
2839 29791260 : la = lmassptr_amode(ll,n)
2840 29791260 : lc = lmassptrcw_amode(ll,n) + pcnst
2841 29791260 : tmp_fact = fm(n)
2842 : end if
2843 :
2844 : if ( (method_reduce_actfrac == 1) .and. &
2845 : (factor_reduce_actfrac >= 0.0_r8) .and. &
2846 : (factor_reduce_actfrac < 1.0_r8) ) &
2847 : tmp_fact = tmp_fact * factor_reduce_actfrac
2848 :
2849 37735596 : delact = min( conu(la)*tmp_fact, conu(la) )
2850 37735596 : delact = max( delact, 0.0_r8 )
2851 37735596 : conu(la) = conu(la) - delact
2852 37735596 : conu(lc) = conu(lc) + delact
2853 37735596 : dconudt(la) = -delact*dt_u_inv
2854 45679932 : dconudt(lc) = delact*dt_u_inv
2855 : end do
2856 : end do ! "n = 1, ntot_amode"
2857 :
2858 : return
2859 2148075 : end subroutine ma_activate_convproc_method2
2860 :
2861 :
2862 :
2863 : !=========================================================================================
2864 167511 : subroutine ma_resuspend_convproc( &
2865 167511 : dcondt, dcondt_resusp, &
2866 : const, dp_i, ktop, kbot_prevap, pcnst_extd )
2867 : !-----------------------------------------------------------------------
2868 : !
2869 : ! Purpose:
2870 : ! Calculate resuspension of activated aerosol species resulting from both
2871 : ! detrainment from updraft and downdraft into environment
2872 : ! subsidence and lifting of environment, which may move air from
2873 : ! levels with large-scale cloud to levels with no large-scale cloud
2874 : !
2875 : ! Method:
2876 : ! Three possible approaches were considered:
2877 : !
2878 : ! 1. Ad-hoc #1 approach. At each level, adjust dcondt for the activated
2879 : ! and unactivated portions of a particular aerosol species so that the
2880 : ! ratio of dcondt (activated/unactivate) is equal to the ratio of the
2881 : ! mixing ratios before convection.
2882 : ! THIS WAS IMPLEMENTED IN MIRAGE2
2883 : !
2884 : ! 2. Ad-hoc #2 approach. At each level, adjust dcondt for the activated
2885 : ! and unactivated portions of a particular aerosol species so that the
2886 : ! change to the activated portion is minimized (zero if possible). The
2887 : ! would minimize effects of convection on the large-scale cloud.
2888 : ! THIS IS CURRENTLY IMPLEMENTED IN CAM5 where we assume that convective
2889 : ! clouds have no impact on the stratiform-cloudborne aerosol
2890 : !
2891 : ! 3. Mechanistic approach that treats the details of interactions between
2892 : ! the large-scale and convective clouds. (Something for the future.)
2893 : !
2894 : ! Author: R. Easter
2895 : !
2896 : !-----------------------------------------------------------------------
2897 :
2898 : use modal_aero_data, only: lmassptr_amode, lmassptrcw_amode, &
2899 2148075 : nspec_amode, ntot_amode, numptr_amode, numptrcw_amode
2900 :
2901 : implicit none
2902 :
2903 : !-----------------------------------------------------------------------
2904 : ! arguments
2905 : ! (note: TMR = tracer mixing ratio)
2906 : integer, intent(in) :: pcnst_extd
2907 : real(r8), intent(inout) :: dcondt(pcnst_extd,pver)
2908 : ! overall TMR tendency from convection
2909 : real(r8), intent(inout) :: dcondt_resusp(pcnst_extd,pver)
2910 : ! portion of TMR tendency due to resuspension
2911 : ! (actually, due to the adjustments made here)
2912 : real(r8), intent(in) :: const(pcnst_extd,pver) ! TMRs before convection
2913 :
2914 : real(r8), intent(in) :: dp_i(pver) ! pressure thickness of level (in mb)
2915 : integer, intent(in) :: ktop, kbot_prevap ! indices of top and bottom cloud levels
2916 :
2917 : !-----------------------------------------------------------------------
2918 : ! local variables
2919 : integer :: k, ll, la, lc, n
2920 : real(r8) :: qa, qc, qac ! working variables (mixing ratios)
2921 : real(r8) :: qdota, qdotc, qdotac ! working variables (MR tendencies)
2922 : !-----------------------------------------------------------------------
2923 :
2924 :
2925 837555 : do n = 1, ntot_amode
2926 :
2927 4020264 : do ll = 0, nspec_amode(n)
2928 3182709 : if (ll == 0) then
2929 670044 : la = numptr_amode(n)
2930 670044 : lc = numptrcw_amode(n) + pcnst
2931 : else
2932 2512665 : la = lmassptr_amode(ll,n)
2933 2512665 : lc = lmassptrcw_amode(ll,n) + pcnst
2934 : end if
2935 :
2936 : ! apply adjustments to dcondt for pairs of unactivated (la) and
2937 : ! activated (lc) aerosol species
2938 3182709 : if ( (la <= 0) .or. (la > pcnst_extd) ) cycle
2939 3182709 : if ( (lc <= 0) .or. (lc > pcnst_extd) ) cycle
2940 :
2941 67035980 : do k = ktop, kbot_prevap
2942 63183227 : qdota = dcondt(la,k)
2943 63183227 : qdotc = dcondt(lc,k)
2944 63183227 : qdotac = qdota + qdotc
2945 :
2946 : ! mirage2 approach
2947 : ! qa = max( const(la,k), 0.0_r8 )
2948 : ! qc = max( const(lc,k), 0.0_r8 )
2949 : ! qac = qa + qc
2950 : ! if (qac <= 0.0) then
2951 : ! dcondt(la,k) = qdotac
2952 : ! dcondt(lc,k) = 0.0
2953 : ! else
2954 : ! dcondt(la,k) = qdotac*(qa/qac)
2955 : ! dcondt(lc,k) = qdotac*(qc/qac)
2956 : ! end if
2957 :
2958 : ! cam5 approach
2959 66365936 : if (convproc_do_evaprain_atonce) then
2960 : dcondt(la,k) = qdota
2961 : dcondt(lc,k) = qdotc
2962 :
2963 63183227 : dcondt_resusp(la,k) = dcondt(la,k)
2964 63183227 : dcondt_resusp(lc,k) = dcondt(lc,k)
2965 : else
2966 0 : dcondt(la,k) = qdotac
2967 0 : dcondt(lc,k) = 0.0_r8
2968 :
2969 0 : dcondt_resusp(la,k) = (dcondt(la,k) - qdota)
2970 0 : dcondt_resusp(lc,k) = (dcondt(lc,k) - qdotc)
2971 : end if
2972 : end do
2973 :
2974 : end do ! "ll = -1, nspec_amode(n)"
2975 : end do ! "n = 1, ntot_amode"
2976 :
2977 167511 : return
2978 167511 : end subroutine ma_resuspend_convproc
2979 :
2980 :
2981 :
2982 : !=========================================================================================
2983 :
2984 :
2985 :
2986 : end module modal_aero_convproc
|
