Line data Source code
1 : module wetdep
2 :
3 : !-----------------------------------------------------------------------
4 : !
5 : ! Wet deposition routines for both aerosols and gas phase constituents.
6 : !
7 : !-----------------------------------------------------------------------
8 :
9 : use shr_kind_mod, only: r8 => shr_kind_r8
10 : use ppgrid, only: pcols, pver
11 : use physconst, only: gravit, rair, tmelt
12 : use phys_control, only: cam_physpkg_is
13 : use cam_logfile, only: iulog
14 : use cam_abortutils, only: endrun
15 :
16 : implicit none
17 : save
18 : private
19 :
20 : public :: wetdepa_v1 ! scavenging codes for very soluble aerosols -- CAM4 version
21 : public :: wetdepa_v2 ! scavenging codes for very soluble aerosols -- CAM5 version
22 : public :: wetdepg ! scavenging of gas phase constituents by henry's law
23 : public :: clddiag ! calc of cloudy volume and rain mixing ratio
24 :
25 : public :: wetdep_inputs_t
26 : public :: wetdep_init
27 : public :: wetdep_inputs_set
28 :
29 : real(r8), parameter :: cmftau = 3600._r8
30 : real(r8), parameter :: rhoh2o = 1000._r8 ! density of water
31 : real(r8), parameter :: molwta = 28.97_r8 ! molecular weight dry air gm/mole
32 : real(r8), parameter :: omsm = 1._r8-2*epsilon(1._r8) ! used to prevent roundoff errors below zero
33 :
34 : type wetdep_inputs_t
35 : real(r8), pointer :: cldt(:,:) => null() ! cloud fraction
36 : real(r8), pointer :: qme(:,:) => null()
37 : real(r8), pointer :: prain(:,:) => null()
38 : real(r8), pointer :: bergso(:,:) => null()
39 : real(r8), pointer :: evapr(:,:) => null()
40 : real(r8) :: cldcu(pcols,pver) ! convective cloud fraction, currently empty
41 : real(r8) :: evapc(pcols,pver) ! Evaporation rate of convective precipitation
42 : real(r8) :: cmfdqr(pcols,pver) ! convective production of rain
43 : real(r8) :: conicw(pcols,pver) ! convective in-cloud water
44 : real(r8) :: totcond(pcols, pver) ! total condensate
45 : real(r8) :: cldv(pcols,pver) ! cloudy volume undergoing wet chem and scavenging
46 : real(r8) :: cldvcu(pcols,pver) ! Convective precipitation area at the top interface of current layer
47 : real(r8) :: cldvst(pcols,pver) ! Stratiform precipitation area at the top interface of current layer
48 : end type wetdep_inputs_t
49 :
50 : integer :: cld_idx = 0
51 : integer :: qme_idx = 0
52 : integer :: prain_idx = 0
53 : integer :: bergso_idx = 0
54 : integer :: nevapr_idx = 0
55 :
56 : integer :: icwmrdp_idx = 0
57 : integer :: icwmrsh_idx = 0
58 : integer :: rprddp_idx = 0
59 : integer :: rprdsh_idx = 0
60 : integer :: sh_frac_idx = 0
61 : integer :: dp_frac_idx = 0
62 : integer :: nevapr_shcu_idx = 0
63 : integer :: nevapr_dpcu_idx = 0
64 : integer :: ixcldice, ixcldliq
65 :
66 : !==============================================================================
67 : contains
68 : !==============================================================================
69 :
70 : !==============================================================================
71 : !==============================================================================
72 1536 : subroutine wetdep_init()
73 : use physics_buffer, only: pbuf_get_index
74 : use constituents, only: cnst_get_ind
75 :
76 : integer :: ierr
77 :
78 1536 : cld_idx = pbuf_get_index('CLD')
79 1536 : qme_idx = pbuf_get_index('QME')
80 1536 : prain_idx = pbuf_get_index('PRAIN')
81 1536 : bergso_idx = pbuf_get_index('BERGSO', errcode=ierr )
82 1536 : nevapr_idx = pbuf_get_index('NEVAPR')
83 :
84 1536 : icwmrdp_idx = pbuf_get_index('ICWMRDP')
85 1536 : rprddp_idx = pbuf_get_index('RPRDDP')
86 1536 : icwmrsh_idx = pbuf_get_index('ICWMRSH')
87 1536 : rprdsh_idx = pbuf_get_index('RPRDSH')
88 1536 : sh_frac_idx = pbuf_get_index('SH_FRAC' )
89 1536 : dp_frac_idx = pbuf_get_index('DP_FRAC')
90 1536 : nevapr_shcu_idx = pbuf_get_index('NEVAPR_SHCU')
91 1536 : nevapr_dpcu_idx = pbuf_get_index('NEVAPR_DPCU')
92 :
93 1536 : call cnst_get_ind('CLDICE', ixcldice)
94 1536 : call cnst_get_ind('CLDLIQ', ixcldliq)
95 :
96 1536 : endsubroutine wetdep_init
97 :
98 : !==============================================================================
99 : ! gathers up the inputs needed for the wetdepa routines
100 : !==============================================================================
101 0 : subroutine wetdep_inputs_set( state, pbuf, inputs )
102 1536 : use physics_types, only: physics_state
103 : use physics_buffer, only: physics_buffer_desc, pbuf_get_field, pbuf_old_tim_idx
104 :
105 : ! args
106 :
107 : type(physics_state), intent(in ) :: state !! physics state
108 : type(physics_buffer_desc), pointer :: pbuf(:) !! physics buffer
109 : type(wetdep_inputs_t), intent(out) :: inputs !! collection of wetdepa inputs
110 :
111 : ! local vars
112 :
113 0 : real(r8), pointer :: icwmrdp(:,:) ! in cloud water mixing ratio, deep convection
114 0 : real(r8), pointer :: rprddp(:,:) ! rain production, deep convection
115 0 : real(r8), pointer :: icwmrsh(:,:) ! in cloud water mixing ratio, deep convection
116 0 : real(r8), pointer :: rprdsh(:,:) ! rain production, deep convection
117 0 : real(r8), pointer :: sh_frac(:,:) ! Shallow convective cloud fraction
118 0 : real(r8), pointer :: dp_frac(:,:) ! Deep convective cloud fraction
119 0 : real(r8), pointer :: evapcsh(:,:) ! Evaporation rate of shallow convective precipitation >=0.
120 0 : real(r8), pointer :: evapcdp(:,:) ! Evaporation rate of deep convective precipitation >=0.
121 :
122 : real(r8) :: rainmr(pcols,pver) ! mixing ratio of rain within cloud volume
123 : real(r8) :: cldst(pcols,pver) ! Stratiform cloud fraction
124 :
125 : integer :: itim, ncol
126 :
127 0 : ncol = state%ncol
128 0 : itim = pbuf_old_tim_idx()
129 :
130 0 : call pbuf_get_field(pbuf, cld_idx, inputs%cldt, start=(/1,1,itim/), kount=(/pcols,pver,1/) )
131 0 : call pbuf_get_field(pbuf, qme_idx, inputs%qme )
132 0 : call pbuf_get_field(pbuf, prain_idx, inputs%prain )
133 0 : call pbuf_get_field(pbuf, nevapr_idx, inputs%evapr )
134 0 : call pbuf_get_field(pbuf, icwmrdp_idx, icwmrdp )
135 0 : call pbuf_get_field(pbuf, icwmrsh_idx, icwmrsh )
136 0 : call pbuf_get_field(pbuf, rprddp_idx, rprddp )
137 0 : call pbuf_get_field(pbuf, rprdsh_idx, rprdsh )
138 0 : call pbuf_get_field(pbuf, sh_frac_idx, sh_frac )
139 0 : call pbuf_get_field(pbuf, dp_frac_idx, dp_frac )
140 0 : call pbuf_get_field(pbuf, nevapr_shcu_idx, evapcsh )
141 0 : call pbuf_get_field(pbuf, nevapr_dpcu_idx, evapcdp )
142 :
143 0 : if (bergso_idx>0) then
144 0 : call pbuf_get_field(pbuf, bergso_idx, inputs%bergso )
145 : else
146 0 : if (.not. associated(inputs%bergso)) then
147 0 : allocate(inputs%bergso(pcols,pver))
148 0 : inputs%bergso(:,:) = 0.0_r8
149 : endif
150 : endif
151 :
152 0 : inputs%cldcu(:ncol,:) = dp_frac(:ncol,:) + sh_frac(:ncol,:)
153 0 : cldst(:ncol,:) = inputs%cldt(:ncol,:) - inputs%cldcu(:ncol,:) ! Stratiform cloud fraction
154 0 : inputs%evapc(:ncol,:) = evapcsh(:ncol,:) + evapcdp(:ncol,:)
155 0 : inputs%cmfdqr(:ncol,:) = rprddp(:ncol,:) + rprdsh(:ncol,:)
156 :
157 : ! sum deep and shallow convection contributions
158 0 : if (cam_physpkg_is('cam5') .or. cam_physpkg_is('cam6')) then
159 : ! Dec.29.2009. Sungsu
160 0 : inputs%conicw(:ncol,:) = (icwmrdp(:ncol,:)*dp_frac(:ncol,:) + icwmrsh(:ncol,:)*sh_frac(:ncol,:))/ &
161 0 : max(0.01_r8, sh_frac(:ncol,:) + dp_frac(:ncol,:))
162 : else
163 0 : inputs%conicw(:ncol,:) = icwmrdp(:ncol,:) + icwmrsh(:ncol,:)
164 : end if
165 :
166 0 : inputs%totcond(:ncol,:) = state%q(:ncol,:,ixcldliq) + state%q(:ncol,:,ixcldice)
167 :
168 : call clddiag( state%t, state%pmid, state%pdel, inputs%cmfdqr, inputs%evapc, &
169 : inputs%cldt, inputs%cldcu, cldst, inputs%qme, inputs%evapr, &
170 : inputs%prain, inputs%cldv, inputs%cldvcu, inputs%cldvst, rainmr, &
171 0 : state%ncol )
172 :
173 0 : end subroutine wetdep_inputs_set
174 :
175 0 : subroutine clddiag(t, pmid, pdel, cmfdqr, evapc, &
176 : cldt, cldcu, cldst, cme, evapr, &
177 : prain, cldv, cldvcu, cldvst, rain, &
178 : ncol)
179 :
180 : ! ------------------------------------------------------------------------------------
181 : ! Estimate the cloudy volume which is occupied by rain or cloud water as
182 : ! the max between the local cloud amount or the
183 : ! sum above of (cloud*positive precip production) sum total precip from above
184 : ! ---------------------------------- x ------------------------
185 : ! sum above of (positive precip ) sum positive precip from above
186 : ! Author: P. Rasch
187 : ! Sungsu Park. Mar.2010
188 : ! ------------------------------------------------------------------------------------
189 :
190 : ! Input arguments:
191 : real(r8), intent(in) :: t(pcols,pver) ! temperature (K)
192 : real(r8), intent(in) :: pmid(pcols,pver) ! pressure at layer midpoints
193 : real(r8), intent(in) :: pdel(pcols,pver) ! pressure difference across layers
194 : real(r8), intent(in) :: cmfdqr(pcols,pver) ! dq/dt due to convective rainout
195 : real(r8), intent(in) :: evapc(pcols,pver) ! Evaporation rate of convective precipitation ( >= 0 )
196 : real(r8), intent(in) :: cldt(pcols,pver) ! total cloud fraction
197 : real(r8), intent(in) :: cldcu(pcols,pver) ! Cumulus cloud fraction
198 : real(r8), intent(in) :: cldst(pcols,pver) ! Stratus cloud fraction
199 : real(r8), intent(in) :: cme(pcols,pver) ! rate of cond-evap within the cloud
200 : real(r8), intent(in) :: evapr(pcols,pver) ! rate of evaporation of falling precipitation (kg/kg/s)
201 : real(r8), intent(in) :: prain(pcols,pver) ! rate of conversion of condensate to precipitation (kg/kg/s)
202 : integer, intent(in) :: ncol
203 :
204 : ! Output arguments:
205 : real(r8), intent(out) :: cldv(pcols,pver) ! fraction occupied by rain or cloud water
206 : real(r8), intent(out) :: cldvcu(pcols,pver) ! Convective precipitation volume
207 : real(r8), intent(out) :: cldvst(pcols,pver) ! Stratiform precipitation volume
208 : real(r8), intent(out) :: rain(pcols,pver) ! mixing ratio of rain (kg/kg)
209 :
210 : ! Local variables:
211 : integer i, k
212 : real(r8) convfw ! used in fallspeed calculation; taken from findmcnew
213 : real(r8) sumppr(pcols) ! precipitation rate (kg/m2-s)
214 : real(r8) sumpppr(pcols) ! sum of positive precips from above
215 : real(r8) cldv1(pcols) ! precip weighted cloud fraction from above
216 : real(r8) lprec ! local production rate of precip (kg/m2/s)
217 : real(r8) lprecp ! local production rate of precip (kg/m2/s) if positive
218 : real(r8) rho ! air density
219 : real(r8) vfall
220 : real(r8) sumppr_cu(pcols) ! Convective precipitation rate (kg/m2-s)
221 : real(r8) sumpppr_cu(pcols) ! Sum of positive convective precips from above
222 : real(r8) cldv1_cu(pcols) ! Convective precip weighted convective cloud fraction from above
223 : real(r8) lprec_cu ! Local production rate of convective precip (kg/m2/s)
224 : real(r8) lprecp_cu ! Local production rate of convective precip (kg/m2/s) if positive
225 : real(r8) sumppr_st(pcols) ! Stratiform precipitation rate (kg/m2-s)
226 : real(r8) sumpppr_st(pcols) ! Sum of positive stratiform precips from above
227 : real(r8) cldv1_st(pcols) ! Stratiform precip weighted stratiform cloud fraction from above
228 : real(r8) lprec_st ! Local production rate of stratiform precip (kg/m2/s)
229 : real(r8) lprecp_st ! Local production rate of stratiform precip (kg/m2/s) if positive
230 : ! -----------------------------------------------------------------------
231 :
232 0 : convfw = 1.94_r8*2.13_r8*sqrt(rhoh2o*gravit*2.7e-4_r8)
233 0 : do i=1,ncol
234 0 : sumppr(i) = 0._r8
235 0 : cldv1(i) = 0._r8
236 0 : sumpppr(i) = 1.e-36_r8
237 0 : sumppr_cu(i) = 0._r8
238 0 : cldv1_cu(i) = 0._r8
239 0 : sumpppr_cu(i) = 1.e-36_r8
240 0 : sumppr_st(i) = 0._r8
241 0 : cldv1_st(i) = 0._r8
242 0 : sumpppr_st(i) = 1.e-36_r8
243 : end do
244 :
245 0 : do k = 1,pver
246 0 : do i = 1,ncol
247 0 : cldv(i,k) = &
248 : max(min(1._r8, &
249 : cldv1(i)/sumpppr(i) &
250 : )*sumppr(i)/sumpppr(i), &
251 : cldt(i,k) &
252 0 : )
253 : lprec = pdel(i,k)/gravit &
254 0 : *(prain(i,k)+cmfdqr(i,k)-evapr(i,k))
255 0 : lprecp = max(lprec,1.e-30_r8)
256 0 : cldv1(i) = cldv1(i) + cldt(i,k)*lprecp
257 0 : sumppr(i) = sumppr(i) + lprec
258 0 : sumpppr(i) = sumpppr(i) + lprecp
259 :
260 : ! For convective precipitation volume at the top interface of each layer. Neglect the current layer.
261 0 : cldvcu(i,k) = max(min(1._r8,cldv1_cu(i)/sumpppr_cu(i))*(sumppr_cu(i)/sumpppr_cu(i)),0._r8)
262 0 : lprec_cu = (pdel(i,k)/gravit)*(cmfdqr(i,k)-evapc(i,k))
263 0 : lprecp_cu = max(lprec_cu,1.e-30_r8)
264 0 : cldv1_cu(i) = cldv1_cu(i) + cldcu(i,k)*lprecp_cu
265 0 : sumppr_cu(i) = sumppr_cu(i) + lprec_cu
266 0 : sumpppr_cu(i) = sumpppr_cu(i) + lprecp_cu
267 :
268 : ! For stratiform precipitation volume at the top interface of each layer. Neglect the current layer.
269 0 : cldvst(i,k) = max(min(1._r8,cldv1_st(i)/sumpppr_st(i))*(sumppr_st(i)/sumpppr_st(i)),0._r8)
270 0 : lprec_st = (pdel(i,k)/gravit)*(prain(i,k)-evapr(i,k))
271 0 : lprecp_st = max(lprec_st,1.e-30_r8)
272 0 : cldv1_st(i) = cldv1_st(i) + cldst(i,k)*lprecp_st
273 0 : sumppr_st(i) = sumppr_st(i) + lprec_st
274 0 : sumpppr_st(i) = sumpppr_st(i) + lprecp_st
275 :
276 0 : rain(i,k) = 0._r8
277 0 : if(t(i,k) > tmelt) then
278 0 : rho = pmid(i,k)/(rair*t(i,k))
279 0 : vfall = convfw/sqrt(rho)
280 0 : rain(i,k) = sumppr(i)/(rho*vfall)
281 0 : if (rain(i,k) < 1.e-14_r8) rain(i,k) = 0._r8
282 : endif
283 : end do
284 : end do
285 :
286 0 : end subroutine clddiag
287 :
288 : !==============================================================================
289 :
290 : ! This is the CAM5 version of wetdepa.
291 :
292 0 : subroutine wetdepa_v2( &
293 : p, q, pdel, cldt, cldc, &
294 : cmfdqr, evapc, conicw, precs, conds, &
295 : evaps, cwat, tracer, deltat, scavt, &
296 : iscavt, cldvcu, cldvst, dlf, fracis, &
297 : sol_fact, ncol, scavcoef, is_strat_cloudborne, qqcw, &
298 : f_act_conv, icscavt, isscavt, bcscavt, bsscavt, &
299 : convproc_do_aer, rcscavt, rsscavt, &
300 : sol_facti_in, sol_factic_in, convproc_do_evaprain_atonce_in, bergso_in )
301 :
302 : !-----------------------------------------------------------------------
303 : !
304 : ! scavenging code for very soluble aerosols
305 : !
306 : !-----------------------------------------------------------------------
307 :
308 : real(r8), intent(in) ::&
309 : p(pcols,pver), &! pressure
310 : q(pcols,pver), &! moisture
311 : pdel(pcols,pver), &! pressure thikness
312 : cldt(pcols,pver), &! total cloud fraction
313 : cldc(pcols,pver), &! convective cloud fraction
314 : cmfdqr(pcols,pver), &! rate of production of convective precip
315 : evapc(pcols,pver), &! Evaporation rate of convective precipitation
316 : conicw(pcols,pver), &! convective cloud water
317 : cwat(pcols,pver), &! cloud water amount
318 : precs(pcols,pver), &! rate of production of stratiform precip
319 : conds(pcols,pver), &! rate of production of condensate
320 : evaps(pcols,pver), &! rate of evaporation of precip
321 : cldvcu(pcols,pver), &! Convective precipitation area at the top interface of each layer
322 : cldvst(pcols,pver), &! Stratiform precipitation area at the top interface of each layer
323 : dlf(pcols,pver), &! Detrainment of convective condensate [kg/kg/s]
324 : deltat, &! time step
325 : tracer(pcols,pver) ! trace species
326 :
327 : ! If subroutine is called with just sol_fact:
328 : ! sol_fact is used for both in- and below-cloud scavenging
329 : ! If subroutine is called with optional argument sol_facti_in:
330 : ! sol_fact is used for below cloud scavenging
331 : ! sol_facti is used for in cloud scavenging
332 :
333 : real(r8), intent(in) :: sol_fact(pcols,pver)
334 : integer, intent(in) :: ncol
335 : real(r8), intent(in) :: scavcoef(pcols,pver) ! Dana and Hales coefficient (/mm) (0.1 if not MODAL_AERO)
336 : real(r8), intent(out) ::&
337 : scavt(pcols,pver), &! scavenging tend
338 : iscavt(pcols,pver), &! incloud scavenging tends
339 : fracis(pcols,pver) ! fraction of species not scavenged
340 :
341 : ! Setting is_strat_cloudborne=.true. indicates that tracer is stratiform-cloudborne aerosol.
342 : ! This is only used by MAM code. The optional args qqcw and f_act_conv are not referenced
343 : ! in this case.
344 : ! Setting is_strat_cloudborne=.false. is being used to indicate that the tracers are the
345 : ! interstitial modal aerosols. In this case the optional qqcw (the cloud borne mixing ratio
346 : ! corresponding to the interstitial aerosol) must be provided, as well as the optional f_act_conv.
347 : logical, intent(in), optional :: is_strat_cloudborne
348 : real(r8), intent(in), optional :: qqcw(pcols,pver)
349 : real(r8), intent(in), optional :: f_act_conv(pcols,pver)
350 :
351 : real(r8), intent(in), optional :: sol_facti_in(pcols,pver) ! solubility factor (frac of aerosol scavenged in cloud)
352 : real(r8), intent(in), optional :: sol_factic_in(pcols,pver) ! sol_facti_in for convective clouds
353 :
354 :
355 : real(r8), intent(out), optional :: icscavt(pcols,pver) ! incloud, convective
356 : real(r8), intent(out), optional :: isscavt(pcols,pver) ! incloud, stratiform
357 : real(r8), intent(out), optional :: bcscavt(pcols,pver) ! below cloud, convective
358 : real(r8), intent(out), optional :: bsscavt(pcols,pver) ! below cloud, stratiform
359 :
360 : ! Setting convproc_do_aer=.true. removes the resuspension term from bcscavt and
361 : ! bsscavt and returns those terms as rcscavt and rsscavt respectively.
362 : logical, intent(in), optional :: convproc_do_aer
363 : real(r8), intent(out), optional :: rcscavt(pcols,pver) ! resuspension, convective
364 : real(r8), intent(out), optional :: rsscavt(pcols,pver) ! resuspension, stratiform
365 : logical, intent(in), optional :: convproc_do_evaprain_atonce_in
366 : real(r8), intent(in), optional :: bergso_in(pcols,pver)
367 :
368 : ! local variables
369 :
370 : integer :: i, k
371 : logical :: out_resuspension
372 :
373 : real(r8) :: clds(pcols) ! stratiform cloud fraction
374 : real(r8) :: fracev(pcols) ! fraction of precip from above that is evaporating
375 : real(r8) :: fracev_cu(pcols) ! Fraction of convective precip from above that is evaporating
376 : real(r8) :: fracp(pcols) ! fraction of cloud water converted to precip
377 : real(r8) :: pdog(pcols) ! work variable (pdel/gravit)
378 : real(r8) :: rpdog(pcols) ! work variable (gravit/pdel)
379 : real(r8) :: precabc(pcols) ! conv precip from above (work array)
380 : real(r8) :: precabs(pcols) ! strat precip from above (work array)
381 : real(r8) :: rat(pcols) ! ratio of amount available to amount removed
382 : real(r8) :: scavab(pcols) ! scavenged tracer flux from above (work array)
383 : real(r8) :: scavabc(pcols) ! scavenged tracer flux from above (work array)
384 : real(r8) :: srcc(pcols) ! tend for convective rain
385 : real(r8) :: srcs(pcols) ! tend for stratiform rain
386 : real(r8) :: srct(pcols) ! work variable
387 :
388 : real(r8) :: fins(pcols) ! fraction of rem. rate by strat rain
389 : real(r8) :: finc(pcols) ! fraction of rem. rate by conv. rain
390 : real(r8) :: conv_scav_ic(pcols) ! convective scavenging incloud
391 : real(r8) :: conv_scav_bc(pcols) ! convective scavenging below cloud
392 : real(r8) :: st_scav_ic(pcols) ! stratiform scavenging incloud
393 : real(r8) :: st_scav_bc(pcols) ! stratiform scavenging below cloud
394 :
395 : real(r8) :: odds(pcols) ! limit on removal rate (proportional to prec)
396 : real(r8) :: dblchek(pcols)
397 : logical :: found
398 :
399 : real(r8) :: trac_qqcw(pcols)
400 : real(r8) :: tracer_incu(pcols)
401 : real(r8) :: tracer_mean(pcols)
402 :
403 : ! For stratiform cloud, cloudborne aerosol is treated explicitly,
404 : ! and sol_facti is 1.0 for cloudborne, 0.0 for interstitial.
405 : ! For convective cloud, cloudborne aerosol is not treated explicitly,
406 : ! and sol_factic is 1.0 for both cloudborne and interstitial.
407 :
408 : real(r8) :: sol_facti(pcols,pver) ! in cloud fraction of aerosol scavenged
409 : real(r8) :: sol_factb(pcols,pver) ! below cloud fraction of aerosol scavenged
410 : real(r8) :: sol_factic(pcols,pver) ! in cloud fraction of aerosol scavenged for convective clouds
411 :
412 : real(r8) :: rdeltat
413 : logical :: convproc_do_evaprain_atonce
414 :
415 : ! ------------------------------------------------------------------------
416 :
417 0 : if (present(convproc_do_evaprain_atonce_in)) then
418 0 : convproc_do_evaprain_atonce = convproc_do_evaprain_atonce_in
419 : else
420 : convproc_do_evaprain_atonce = .false.
421 : endif
422 :
423 : ! default (if other sol_facts aren't in call, set all to required sol_fact)
424 0 : sol_facti = sol_fact
425 0 : sol_factb = sol_fact
426 :
427 0 : if ( present(sol_facti_in) ) sol_facti = sol_facti_in
428 :
429 0 : sol_factic = sol_facti
430 0 : if ( present(sol_factic_in ) ) sol_factic = sol_factic_in
431 :
432 : ! Determine whether resuspension fields are output.
433 0 : out_resuspension = .false.
434 0 : if (present(convproc_do_aer)) then
435 0 : if (convproc_do_aer) then
436 : if (present(bcscavt) .and. present(bsscavt) .and. &
437 0 : present(rcscavt) .and. present(rsscavt) ) then
438 : out_resuspension = .true.
439 : else
440 : call endrun('wetdepa_v2: bcscavt, bsscavt, rcscavt, rsscavt'// &
441 0 : ' must be present when convproc_do_aero true')
442 : end if
443 : end if
444 : end if
445 :
446 : ! this section of code is for highly soluble aerosols,
447 : ! the assumption is that within the cloud that
448 : ! all the tracer is in the cloud water
449 : !
450 : ! for both convective and stratiform clouds,
451 : ! the fraction of cloud water converted to precip defines
452 : ! the amount of tracer which is pulled out.
453 :
454 0 : precabs(:ncol) = 0.0_r8
455 0 : precabc(:ncol) = 0.0_r8
456 0 : scavab(:ncol) = 0.0_r8
457 0 : scavabc(:ncol) = 0.0_r8
458 :
459 0 : do k = 1, pver
460 0 : do i = 1, ncol
461 :
462 0 : clds(i) = cldt(i,k) - cldc(i,k)
463 0 : pdog(i) = pdel(i,k)/gravit
464 0 : rpdog(i) = gravit/pdel(i,k)
465 0 : rdeltat = 1.0_r8/deltat
466 :
467 : ! ****************** Evaporation **************************
468 : ! calculate the fraction of strat precip from above
469 : ! which evaporates within this layer
470 : fracev(i) = evaps(i,k)*pdog(i) &
471 0 : /max(1.e-12_r8,precabs(i))
472 :
473 : ! If resuspending aerosol only when all the rain has totally
474 : ! evaporated then zero out any aerosol tendency for partial
475 : ! evaporation.
476 0 : if (convproc_do_evaprain_atonce .and. (evaps(i,k)*pdog(i)/=precabs(i))) then
477 0 : fracev(i) = 0._r8
478 : end if
479 :
480 : ! trap to ensure reasonable ratio bounds
481 0 : fracev(i) = max(0._r8,min(1._r8,fracev(i)))
482 :
483 : ! Same as above but convective precipitation part
484 0 : fracev_cu(i) = evapc(i,k)*pdog(i)/max(1.e-12_r8,precabc(i))
485 0 : fracev_cu(i) = max(0._r8,min(1._r8,fracev_cu(i)))
486 :
487 : ! ****************** Convection ***************************
488 : !
489 : ! set odds proportional to fraction of the grid box that is swept by the
490 : ! precipitation =precabc/rhoh20*(area of sphere projected on plane
491 : ! /volume of sphere)*deltat
492 : ! assume the radius of a raindrop is 1 e-3 m from Rogers and Yau,
493 : ! unless the fraction of the area that is cloud is less than odds, in which
494 : ! case use the cloud fraction (assumes precabs is in kg/m2/s)
495 : ! is really: precabs*3/4/1000./1e-3*deltat
496 : ! here I use .1 from Balkanski
497 : !
498 : ! use a local rate of convective rain production for incloud scav
499 : !
500 : ! Fraction of convective cloud water converted to rain. This version is used
501 : ! in 2 of the 3 branches below before fracp is reused in the stratiform calc.
502 : ! NB: In below formula for fracp conicw is a LWC/IWC that has already
503 : ! precipitated out, i.e., conicw does not contain precipitation
504 :
505 : fracp(i) = cmfdqr(i,k)*deltat / &
506 0 : max( 1.e-12_r8, cldc(i,k)*conicw(i,k) + (cmfdqr(i,k)+dlf(i,k))*deltat )
507 0 : fracp(i) = max( min( 1._r8, fracp(i)), 0._r8 )
508 :
509 0 : if ( present(is_strat_cloudborne) ) then
510 :
511 0 : if ( is_strat_cloudborne ) then
512 :
513 : ! convective scavenging
514 :
515 0 : conv_scav_ic(i) = 0._r8
516 :
517 0 : conv_scav_bc(i) = 0._r8
518 :
519 : ! stratiform scavenging
520 0 : if (convproc_do_evaprain_atonce .and. present(bergso_in)) then
521 : fracp(i) = (precs(i,k)-bergso_in(i,k))*deltat / &
522 0 : max( 1.e-12_r8, cwat(i,k) + precs(i,k)*deltat )
523 : else
524 : fracp(i) = precs(i,k)*deltat / &
525 0 : max( 1.e-12_r8, cwat(i,k) + precs(i,k)*deltat )
526 : endif
527 :
528 0 : fracp(i) = max( 0._r8, min(1._r8, fracp(i)) )
529 :
530 0 : st_scav_ic(i) = sol_facti(i,k) *fracp(i)*tracer(i,k)*rdeltat
531 :
532 0 : st_scav_bc(i) = 0._r8
533 :
534 : else
535 :
536 : ! convective scavenging
537 :
538 : trac_qqcw(i) = min(qqcw(i,k), &
539 0 : tracer(i,k)*( clds(i)/max( 0.01_r8, 1._r8-clds(i) ) ) )
540 :
541 0 : tracer_incu(i) = f_act_conv(i,k)*(tracer(i,k) + trac_qqcw(i))
542 :
543 0 : conv_scav_ic(i) = sol_factic(i,k)*cldc(i,k)*fracp(i)*tracer_incu(i)*rdeltat
544 :
545 : tracer_mean(i) = tracer(i,k)*(1._r8 - cldc(i,k)*f_act_conv(i,k)) - &
546 0 : cldc(i,k)*f_act_conv(i,k)*trac_qqcw(i)
547 0 : tracer_mean(i) = max(0._r8,tracer_mean(i))
548 :
549 0 : odds(i) = precabc(i)/max(cldvcu(i,k),1.e-5_r8)*scavcoef(i,k)*deltat
550 0 : odds(i) = max(min(1._r8,odds(i)),0._r8)
551 0 : conv_scav_bc(i) = sol_factb(i,k) *cldvcu(i,k)*odds(i)*tracer_mean(i)*rdeltat
552 :
553 :
554 : ! stratiform scavenging
555 :
556 0 : st_scav_ic(i) = 0._r8
557 :
558 0 : odds(i) = precabs(i)/max(cldvst(i,k),1.e-5_r8)*scavcoef(i,k)*deltat
559 0 : odds(i) = max(min(1._r8,odds(i)),0._r8)
560 0 : st_scav_bc(i) = sol_factb(i,k) *cldvst(i,k)*odds(i)*tracer_mean(i)*rdeltat
561 :
562 : end if
563 :
564 : else
565 :
566 : ! convective scavenging
567 :
568 0 : conv_scav_ic(i) = sol_factic(i,k)*cldc(i,k)*fracp(i)*tracer(i,k)*rdeltat
569 :
570 0 : odds(i) = precabc(i)/max(cldvcu(i,k), 1.e-5_r8)*scavcoef(i,k)*deltat
571 0 : odds(i) = max( min(1._r8, odds(i)), 0._r8)
572 0 : conv_scav_bc(i) = sol_factb(i,k)*cldvcu(i,k)*odds(i)*tracer(i,k)*rdeltat
573 :
574 : ! stratiform scavenging
575 :
576 : ! fracp is the fraction of cloud water converted to precip
577 : ! NB: In below formula for fracp cwat is a LWC/IWC that has already
578 : ! precipitated out, i.e., cwat does not contain precipitation
579 : fracp(i) = precs(i,k)*deltat / &
580 0 : max( 1.e-12_r8, cwat(i,k) + precs(i,k)*deltat )
581 0 : fracp(i) = max( 0._r8, min( 1._r8, fracp(i) ) )
582 :
583 : ! assume the corresponding amnt of tracer is removed
584 0 : st_scav_ic(i) = sol_facti(i,k)*clds(i)*fracp(i)*tracer(i,k)*rdeltat
585 :
586 0 : odds(i) = precabs(i)/max(cldvst(i,k),1.e-5_r8)*scavcoef(i,k)*deltat
587 0 : odds(i) = max(min(1._r8,odds(i)),0._r8)
588 0 : st_scav_bc(i) =sol_factb(i,k)*(cldvst(i,k)*odds(i)) *tracer(i,k)*rdeltat
589 :
590 : end if
591 :
592 : ! total convective scavenging
593 0 : srcc(i) = conv_scav_ic(i) + conv_scav_bc(i)
594 0 : finc(i) = conv_scav_ic(i)/(srcc(i) + 1.e-36_r8)
595 :
596 : ! total stratiform scavenging
597 0 : srcs(i) = st_scav_ic(i) + st_scav_bc(i)
598 0 : fins(i) = st_scav_ic(i)/(srcs(i) + 1.e-36_r8)
599 :
600 : ! make sure we dont take out more than is there
601 : ! ratio of amount available to amount removed
602 0 : rat(i) = tracer(i,k)/max(deltat*(srcc(i)+srcs(i)),1.e-36_r8)
603 0 : if (rat(i)<1._r8) then
604 0 : srcs(i) = srcs(i)*rat(i)
605 0 : srcc(i) = srcc(i)*rat(i)
606 : endif
607 0 : srct(i) = (srcc(i)+srcs(i))*omsm
608 :
609 :
610 : ! fraction that is not removed within the cloud
611 : ! (assumed to be interstitial, and subject to convective transport)
612 0 : fracp(i) = deltat*srct(i)/max(cldvst(i,k)*tracer(i,k),1.e-36_r8) ! amount removed
613 0 : fracp(i) = max(0._r8,min(1._r8,fracp(i)))
614 0 : fracis(i,k) = 1._r8 - fracp(i)
615 :
616 : ! tend is all tracer removed by scavenging, plus all re-appearing from evaporation above
617 : ! Sungsu added cumulus contribution in the below 3 blocks
618 0 : scavt(i,k) = -srct(i) + (fracev(i)*scavab(i)+fracev_cu(i)*scavabc(i))*rpdog(i)
619 0 : iscavt(i,k) = -(srcc(i)*finc(i) + srcs(i)*fins(i))*omsm
620 :
621 0 : if ( present(icscavt) ) icscavt(i,k) = -(srcc(i)*finc(i)) * omsm
622 0 : if ( present(isscavt) ) isscavt(i,k) = -(srcs(i)*fins(i)) * omsm
623 :
624 0 : if (.not. out_resuspension) then
625 0 : if (present(bcscavt)) bcscavt(i,k) = -(srcc(i) * (1-finc(i))) * omsm + &
626 0 : fracev_cu(i)*scavabc(i)*rpdog(i)
627 :
628 0 : if (present(bsscavt)) bsscavt(i,k) = -(srcs(i) * (1-fins(i))) * omsm + &
629 0 : fracev(i)*scavab(i)*rpdog(i)
630 : else
631 0 : bcscavt(i,k) = -(srcc(i) * (1-finc(i))) * omsm
632 0 : rcscavt(i,k) = fracev_cu(i)*scavabc(i)*rpdog(i)
633 :
634 0 : bsscavt(i,k) = -(srcs(i) * (1-fins(i))) * omsm
635 0 : rsscavt(i,k) = fracev(i)*scavab(i)*rpdog(i)
636 : end if
637 :
638 0 : dblchek(i) = tracer(i,k) + deltat*scavt(i,k)
639 :
640 : ! now keep track of scavenged mass and precip
641 0 : scavab(i) = scavab(i)*(1-fracev(i)) + srcs(i)*pdog(i)
642 0 : precabs(i) = precabs(i) + (precs(i,k) - evaps(i,k))*pdog(i)
643 0 : scavabc(i) = scavabc(i)*(1-fracev_cu(i)) + srcc(i)*pdog(i)
644 0 : precabc(i) = precabc(i) + (cmfdqr(i,k) - evapc(i,k))*pdog(i)
645 :
646 : end do ! End of i = 1, ncol
647 :
648 0 : found = .false.
649 0 : do i = 1,ncol
650 0 : if ( dblchek(i) < 0._r8 ) then
651 : found = .true.
652 : exit
653 : end if
654 : end do
655 :
656 0 : if ( found ) then
657 0 : do i = 1,ncol
658 0 : if (dblchek(i) < 0._r8) then
659 0 : write(iulog,*) ' wetdapa: negative value ', i, k, tracer(i,k), &
660 0 : dblchek(i), scavt(i,k), srct(i), rat(i), fracev(i)
661 : endif
662 : end do
663 : endif
664 :
665 : end do ! End of k = 1, pver
666 :
667 0 : end subroutine wetdepa_v2
668 :
669 :
670 : !==============================================================================
671 :
672 : ! This is the frozen CAM4 version of wetdepa.
673 :
674 :
675 0 : subroutine wetdepa_v1( t, p, q, pdel, &
676 : cldt, cldc, cmfdqr, conicw, precs, conds, &
677 : evaps, cwat, tracer, deltat, &
678 : scavt, iscavt, cldv, fracis, sol_fact, ncol, &
679 : scavcoef,icscavt, isscavt, bcscavt, bsscavt, &
680 : sol_facti_in, sol_factbi_in, sol_factii_in, &
681 : sol_factic_in, sol_factiic_in )
682 :
683 : !-----------------------------------------------------------------------
684 : ! Purpose:
685 : ! scavenging code for very soluble aerosols
686 : !
687 : ! Author: P. Rasch
688 : ! Modified by T. Bond 3/2003 to track different removals
689 : !-----------------------------------------------------------------------
690 :
691 : implicit none
692 :
693 : real(r8), intent(in) ::&
694 : t(pcols,pver), &! temperature
695 : p(pcols,pver), &! pressure
696 : q(pcols,pver), &! moisture
697 : pdel(pcols,pver), &! pressure thikness
698 : cldt(pcols,pver), &! total cloud fraction
699 : cldc(pcols,pver), &! convective cloud fraction
700 : cmfdqr(pcols,pver), &! rate of production of convective precip
701 : conicw(pcols,pver), &! convective cloud water
702 : cwat(pcols,pver), &! cloud water amount
703 : precs(pcols,pver), &! rate of production of stratiform precip
704 : conds(pcols,pver), &! rate of production of condensate
705 : evaps(pcols,pver), &! rate of evaporation of precip
706 : cldv(pcols,pver), &! total cloud fraction
707 : deltat, &! time step
708 : tracer(pcols,pver) ! trace species
709 : ! If subroutine is called with just sol_fact:
710 : ! sol_fact is used for both in- and below-cloud scavenging
711 : ! If subroutine is called with optional argument sol_facti_in:
712 : ! sol_fact is used for below cloud scavenging
713 : ! sol_facti is used for in cloud scavenging
714 : real(r8), intent(in) :: sol_fact ! solubility factor (fraction of aer scavenged below & in, or just below or sol_facti_in is provided)
715 : real(r8), intent(in), optional :: sol_facti_in ! solubility factor (frac of aerosol scavenged in cloud)
716 : real(r8), intent(in), optional :: sol_factbi_in ! solubility factor (frac of aerosol scavenged below cloud by ice)
717 : real(r8), intent(in), optional :: sol_factii_in ! solubility factor (frac of aerosol scavenged in cloud by ice)
718 : real(r8), intent(in), optional :: sol_factic_in(pcols,pver) ! sol_facti_in for convective clouds
719 : real(r8), intent(in), optional :: sol_factiic_in ! sol_factii_in for convective clouds
720 : real(r8), intent(in) :: scavcoef(pcols,pver) ! Dana and Hales coefficient (/mm) (0.1 if not MODAL_AERO)
721 :
722 : integer, intent(in) :: ncol
723 :
724 : real(r8), intent(out) ::&
725 : scavt(pcols,pver), &! scavenging tend
726 : iscavt(pcols,pver), &! incloud scavenging tends
727 : fracis(pcols,pver) ! fraction of species not scavenged
728 :
729 : real(r8), intent(out), optional :: icscavt(pcols,pver) ! incloud, convective
730 : real(r8), intent(out), optional :: isscavt(pcols,pver) ! incloud, stratiform
731 : real(r8), intent(out), optional :: bcscavt(pcols,pver) ! below cloud, convective
732 : real(r8), intent(out), optional :: bsscavt(pcols,pver) ! below cloud, stratiform
733 :
734 : ! local variables
735 :
736 : integer i ! x index
737 : integer k ! z index
738 :
739 : real(r8) adjfac ! factor stolen from cmfmca
740 : real(r8) aqfrac ! fraction of tracer in aqueous phase
741 : real(r8) cwatc ! local convective total water amount
742 : real(r8) cwats ! local stratiform total water amount
743 : real(r8) cwatp ! local water amount falling from above precip
744 : real(r8) fracev(pcols) ! fraction of precip from above that is evaporating
745 : real(r8) fracp ! fraction of cloud water converted to precip
746 : real(r8) gafrac ! fraction of tracer in gas phasea
747 : real(r8) hconst ! henry's law solubility constant when equation is expressed
748 : ! in terms of mixing ratios
749 : real(r8) mpla ! moles / liter H2O entering the layer from above
750 : real(r8) mplb ! moles / liter H2O leaving the layer below
751 : real(r8) part ! partial pressure of tracer in atmospheres
752 : real(r8) patm ! total pressure in atmospheres
753 : real(r8) pdog ! work variable (pdel/gravit)
754 : real(r8) precabc(pcols) ! conv precip from above (work array)
755 : real(r8) precabs(pcols) ! strat precip from above (work array)
756 : real(r8) precbl ! precip falling out of level (work array)
757 : real(r8) precmin ! minimum convective precip causing scavenging
758 : real(r8) rat(pcols) ! ratio of amount available to amount removed
759 : real(r8) scavab(pcols) ! scavenged tracer flux from above (work array)
760 : real(r8) scavabc(pcols) ! scavenged tracer flux from above (work array)
761 : real(r8) srcc ! tend for convective rain
762 : real(r8) srcs ! tend for stratiform rain
763 : real(r8) srct(pcols) ! work variable
764 : real(r8) tracab(pcols) ! column integrated tracer amount
765 :
766 : real(r8) fins ! fraction of rem. rate by strat rain
767 : real(r8) finc ! fraction of rem. rate by conv. rain
768 : real(r8) srcs1 ! work variable
769 : real(r8) srcs2 ! work variable
770 : real(r8) tc ! temp in celcius
771 : real(r8) weight ! fraction of condensate which is ice
772 : real(r8) cldmabs(pcols) ! maximum cloud at or above this level
773 : real(r8) cldmabc(pcols) ! maximum cloud at or above this level
774 : real(r8) odds ! limit on removal rate (proportional to prec)
775 : real(r8) dblchek(pcols)
776 : logical :: found
777 :
778 : real(r8) sol_facti, sol_factb ! in cloud and below cloud fraction of aerosol scavenged
779 : real(r8) sol_factii, sol_factbi ! in cloud and below cloud fraction of aerosol scavenged by ice
780 : real(r8) sol_factic(pcols,pver) ! sol_facti for convective clouds
781 : real(r8) sol_factiic ! sol_factii for convective clouds
782 : ! sol_factic & solfact_iic added for MODAL_AERO.
783 : ! For stratiform cloud, cloudborne aerosol is treated explicitly,
784 : ! and sol_facti is 1.0 for cloudborne, 0.0 for interstitial.
785 : ! For convective cloud, cloudborne aerosol is not treated explicitly,
786 : ! and sol_factic is 1.0 for both cloudborne and interstitial.
787 :
788 : ! ------------------------------------------------------------------------
789 0 : precmin = 0.1_r8/8.64e4_r8 ! set critical value to 0.1 mm/day in kg/m2/s
790 :
791 0 : adjfac = deltat/(max(deltat,cmftau)) ! adjustment factor from hack scheme
792 :
793 : ! default (if other sol_facts aren't in call, set all to required sol_fact
794 0 : sol_facti = sol_fact
795 0 : sol_factb = sol_fact
796 0 : sol_factii = sol_fact
797 0 : sol_factbi = sol_fact
798 :
799 0 : if ( present(sol_facti_in) ) sol_facti = sol_facti_in
800 0 : if ( present(sol_factii_in) ) sol_factii = sol_factii_in
801 0 : if ( present(sol_factbi_in) ) sol_factbi = sol_factbi_in
802 :
803 0 : sol_factic = sol_facti
804 0 : sol_factiic = sol_factii
805 0 : if ( present(sol_factic_in ) ) sol_factic = sol_factic_in
806 0 : if ( present(sol_factiic_in) ) sol_factiic = sol_factiic_in
807 :
808 : ! this section of code is for highly soluble aerosols,
809 : ! the assumption is that within the cloud that
810 : ! all the tracer is in the cloud water
811 : !
812 : ! for both convective and stratiform clouds,
813 : ! the fraction of cloud water converted to precip defines
814 : ! the amount of tracer which is pulled out.
815 : !
816 :
817 0 : do i = 1,pcols
818 0 : precabs(i) = 0
819 0 : precabc(i) = 0
820 0 : scavab(i) = 0
821 0 : scavabc(i) = 0
822 0 : tracab(i) = 0
823 0 : cldmabs(i) = 0
824 0 : cldmabc(i) = 0
825 : end do
826 :
827 0 : do k = 1,pver
828 0 : do i = 1,ncol
829 0 : tc = t(i,k) - tmelt
830 : weight = max(0._r8,min(-tc*0.05_r8,1.0_r8)) ! fraction of condensate that is ice
831 0 : weight = 0._r8 ! assume no ice
832 :
833 0 : pdog = pdel(i,k)/gravit
834 :
835 : ! ****************** Evaporation **************************
836 : ! calculate the fraction of strat precip from above
837 : ! which evaporates within this layer
838 : fracev(i) = evaps(i,k)*pdel(i,k)/gravit &
839 0 : /max(1.e-12_r8,precabs(i))
840 :
841 : ! trap to ensure reasonable ratio bounds
842 0 : fracev(i) = max(0._r8,min(1._r8,fracev(i)))
843 :
844 : ! ****************** Convection ***************************
845 : ! now do the convective scavenging
846 :
847 : ! set odds proportional to fraction of the grid box that is swept by the
848 : ! precipitation =precabc/rhoh20*(area of sphere projected on plane
849 : ! /volume of sphere)*deltat
850 : ! assume the radius of a raindrop is 1 e-3 m from Rogers and Yau,
851 : ! unless the fraction of the area that is cloud is less than odds, in which
852 : ! case use the cloud fraction (assumes precabs is in kg/m2/s)
853 : ! is really: precabs*3/4/1000./1e-3*deltat
854 : ! here I use .1 from Balkanski
855 : !
856 : ! use a local rate of convective rain production for incloud scav
857 : !odds=max(min(1._r8, &
858 : ! cmfdqr(i,k)*pdel(i,k)/gravit*0.1_r8*deltat),0._r8)
859 : !++mcb -- change cldc to cldt; change cldt to cldv (9/17/96)
860 : ! srcs1 = cldt(i,k)*odds*tracer(i,k)*(1.-weight) &
861 : ! srcs1 = cldv(i,k)*odds*tracer(i,k)*(1.-weight) &
862 : !srcs1 = cldc(i,k)*odds*tracer(i,k)*(1.-weight) &
863 : ! /deltat
864 :
865 : ! fraction of convective cloud water converted to rain
866 0 : fracp = cmfdqr(i,k)*deltat/max(1.e-8_r8,conicw(i,k))
867 : ! note cmfdrq can be negative from evap of rain, so constrain it
868 0 : fracp = max(min(1._r8,fracp),0._r8)
869 : ! remove that amount from within the convective area
870 : ! srcs1 = cldc(i,k)*fracp*tracer(i,k)*(1._r8-weight)/deltat ! liquid only
871 : ! srcs1 = cldc(i,k)*fracp*tracer(i,k)/deltat ! any condensation
872 : ! srcs1 = 0.
873 : srcs1 = sol_factic(i,k)*cldt(i,k)*fracp*tracer(i,k)*(1._r8-weight)/deltat & ! liquid
874 0 : + sol_factiic*cldt(i,k)*fracp*tracer(i,k)*(weight)/deltat ! ice
875 :
876 :
877 : !--mcb
878 :
879 : ! scavenge below cloud
880 :
881 : ! cldmabc(i) = max(cldc(i,k),cldmabc(i))
882 : ! cldmabc(i) = max(cldt(i,k),cldmabc(i))
883 0 : cldmabc(i) = max(cldv(i,k),cldmabc(i))
884 0 : cldmabc(i) = cldv(i,k)
885 :
886 : odds=max( &
887 : min(1._r8,precabc(i)/max(cldmabc(i),1.e-5_r8) &
888 0 : *scavcoef(i,k)*deltat),0._r8) ! Dana and Hales coefficient (/mm)
889 : srcs2 = sol_factb*cldmabc(i)*odds*tracer(i,k)*(1._r8-weight)/deltat & ! liquid
890 0 : + sol_factbi*cldmabc(i)*odds*tracer(i,k)*(weight)/deltat !ice
891 : !Note that using the temperature-determined weight doesn't make much sense here
892 :
893 :
894 0 : srcc = srcs1 + srcs2 ! convective tend by both processes
895 0 : finc = srcs1/(srcc + 1.e-36_r8) ! fraction in-cloud
896 :
897 : ! ****************** Stratiform ***********************
898 : ! now do the stratiform scavenging
899 :
900 : ! incloud scavenging
901 :
902 : ! fracp is the fraction of cloud water converted to precip
903 0 : fracp = precs(i,k)*deltat/max(cwat(i,k),1.e-12_r8)
904 0 : fracp = max(0._r8,min(1._r8,fracp))
905 : ! fracp = 0. ! for debug
906 :
907 : ! assume the corresponding amnt of tracer is removed
908 : !++mcb -- remove cldc; change cldt to cldv
909 : ! srcs1 = (cldt(i,k)-cldc(i,k))*fracp*tracer(i,k)/deltat
910 : ! srcs1 = cldv(i,k)*fracp*tracer(i,k)/deltat &
911 : ! srcs1 = cldt(i,k)*fracp*tracer(i,k)/deltat ! all condensate
912 : srcs1 = sol_facti*cldt(i,k)*fracp*tracer(i,k)/deltat*(1._r8-weight) & ! liquid
913 0 : + sol_factii*cldt(i,k)*fracp*tracer(i,k)/deltat*(weight) ! ice
914 :
915 :
916 : ! below cloud scavenging
917 :
918 : ! volume undergoing below cloud scavenging
919 0 : cldmabs(i) = cldv(i,k) ! precipitating volume
920 : ! cldmabs(i) = cldt(i,k) ! local cloud volume
921 :
922 0 : odds = precabs(i)/max(cldmabs(i),1.e-5_r8)*scavcoef(i,k)*deltat
923 0 : odds = max(min(1._r8,odds),0._r8)
924 : srcs2 =sol_factb*(cldmabs(i)*odds) *tracer(i,k)*(1._r8-weight)/deltat & ! liquid
925 0 : + sol_factbi*(cldmabs(i)*odds) *tracer(i,k)*(weight)/deltat ! ice
926 : !Note that using the temperature-determined weight doesn't make much sense here
927 :
928 :
929 0 : srcs = srcs1 + srcs2 ! total stratiform scavenging
930 0 : fins=srcs1/(srcs + 1.e-36_r8) ! fraction taken by incloud processes
931 :
932 : ! make sure we dont take out more than is there
933 : ! ratio of amount available to amount removed
934 0 : rat(i) = tracer(i,k)/max(deltat*(srcc+srcs),1.e-36_r8)
935 0 : if (rat(i) < 1._r8) then
936 0 : srcs = srcs*rat(i)
937 0 : srcc = srcc*rat(i)
938 : endif
939 0 : srct(i) = (srcc+srcs)*omsm
940 :
941 :
942 : ! fraction that is not removed within the cloud
943 : ! (assumed to be interstitial, and subject to convective transport)
944 0 : fracp = deltat*srct(i)/max(cldmabs(i)*tracer(i,k),1.e-36_r8) ! amount removed
945 0 : fracp = max(0._r8,min(1._r8,fracp))
946 0 : fracis(i,k) = 1._r8 - fracp
947 :
948 : ! tend is all tracer removed by scavenging, plus all re-appearing from evaporation above
949 0 : scavt(i,k) = -srct(i) + fracev(i)*scavab(i)*gravit/pdel(i,k)
950 0 : iscavt(i,k) = -(srcc*finc + srcs*fins)*omsm
951 :
952 0 : if ( present(icscavt) ) icscavt(i,k) = -(srcc*finc) * omsm
953 0 : if ( present(isscavt) ) isscavt(i,k) = -(srcs*fins) * omsm
954 0 : if ( present(bcscavt) ) bcscavt(i,k) = -(srcc * (1-finc)) * omsm
955 0 : if ( present(bsscavt) ) bsscavt(i,k) = -(srcs * (1-fins)) * omsm + &
956 0 : fracev(i)*scavab(i)*gravit/pdel(i,k)
957 :
958 0 : dblchek(i) = tracer(i,k) + deltat*scavt(i,k)
959 :
960 : ! now keep track of scavenged mass and precip
961 0 : scavab(i) = scavab(i)*(1-fracev(i)) + srcs*pdel(i,k)/gravit
962 0 : precabs(i) = precabs(i) + (precs(i,k) - evaps(i,k))*pdel(i,k)/gravit
963 0 : scavabc(i) = scavabc(i) + srcc*pdel(i,k)/gravit
964 0 : precabc(i) = precabc(i) + (cmfdqr(i,k))*pdel(i,k)/gravit
965 0 : tracab(i) = tracab(i) + tracer(i,k)*pdel(i,k)/gravit
966 :
967 : end do
968 :
969 0 : found = .false.
970 0 : do i = 1,ncol
971 0 : if ( dblchek(i) < 0._r8 ) then
972 : found = .true.
973 : exit
974 : end if
975 : end do
976 :
977 0 : if ( found ) then
978 0 : do i = 1,ncol
979 0 : if (dblchek(i) < 0._r8) then
980 0 : write(iulog,*) ' wetdapa: negative value ', i, k, tracer(i,k), &
981 0 : dblchek(i), scavt(i,k), srct(i), rat(i), fracev(i)
982 : endif
983 : end do
984 : endif
985 :
986 : end do
987 :
988 0 : end subroutine wetdepa_v1
989 :
990 : !==============================================================================
991 :
992 : ! wetdepg is currently being used for both CAM4 and CAM5 by making use of the
993 : ! cam_physpkg_is method.
994 :
995 0 : subroutine wetdepg( t, p, q, pdel, &
996 : cldt, cldc, cmfdqr, evapc, precs, evaps, &
997 : rain, cwat, tracer, deltat, molwt, &
998 : solconst, scavt, iscavt, cldv, icwmr1, &
999 : icwmr2, fracis, ncol )
1000 :
1001 : !-----------------------------------------------------------------------
1002 : ! Purpose:
1003 : ! scavenging of gas phase constituents by henry's law
1004 : !
1005 : ! Author: P. Rasch
1006 : !-----------------------------------------------------------------------
1007 :
1008 : real(r8), intent(in) ::&
1009 : t(pcols,pver), &! temperature
1010 : p(pcols,pver), &! pressure
1011 : q(pcols,pver), &! moisture
1012 : pdel(pcols,pver), &! pressure thikness
1013 : cldt(pcols,pver), &! total cloud fraction
1014 : cldc(pcols,pver), &! convective cloud fraction
1015 : cmfdqr(pcols,pver), &! rate of production of convective precip
1016 : rain (pcols,pver), &! total rainwater mixing ratio
1017 : cwat(pcols,pver), &! cloud water amount
1018 : precs(pcols,pver), &! rate of production of stratiform precip
1019 : evaps(pcols,pver), &! rate of evaporation of precip
1020 : ! Sungsu
1021 : evapc(pcols,pver), &! Rate of evaporation of convective precipitation
1022 : ! Sungsu
1023 : cldv(pcols,pver), &! estimate of local volume occupied by clouds
1024 : icwmr1 (pcols,pver), &! in cloud water mixing ration for zhang scheme
1025 : icwmr2 (pcols,pver), &! in cloud water mixing ration for hack scheme
1026 : deltat, &! time step
1027 : tracer(pcols,pver), &! trace species
1028 : molwt ! molecular weights
1029 :
1030 : integer, intent(in) :: ncol
1031 :
1032 : real(r8) &
1033 : solconst(pcols,pver) ! Henry's law coefficient
1034 :
1035 : real(r8), intent(out) ::&
1036 : scavt(pcols,pver), &! scavenging tend
1037 : iscavt(pcols,pver), &! incloud scavenging tends
1038 : fracis(pcols, pver) ! fraction of constituent that is insoluble
1039 :
1040 : ! local variables
1041 :
1042 : integer i ! x index
1043 : integer k ! z index
1044 :
1045 : real(r8) adjfac ! factor stolen from cmfmca
1046 : real(r8) aqfrac ! fraction of tracer in aqueous phase
1047 : real(r8) cwatc ! local convective total water amount
1048 : real(r8) cwats ! local stratiform total water amount
1049 : real(r8) cwatl ! local cloud liq water amount
1050 : real(r8) cwatp ! local water amount falling from above precip
1051 : real(r8) cwatpl ! local water amount falling from above precip (liq)
1052 : real(r8) cwatt ! local sum of strat + conv total water amount
1053 : real(r8) cwatti ! cwatt/cldv = cloudy grid volume mixing ratio
1054 : real(r8) fracev ! fraction of precip from above that is evaporating
1055 : real(r8) fracp ! fraction of cloud water converted to precip
1056 : real(r8) gafrac ! fraction of tracer in gas phasea
1057 : real(r8) hconst ! henry's law solubility constant when equation is expressed
1058 : ! in terms of mixing ratios
1059 : real(r8) mpla ! moles / liter H2O entering the layer from above
1060 : real(r8) mplb ! moles / liter H2O leaving the layer below
1061 : real(r8) part ! partial pressure of tracer in atmospheres
1062 : real(r8) patm ! total pressure in atmospheres
1063 : real(r8) pdog ! work variable (pdel/gravit)
1064 : real(r8) precab(pcols) ! precip from above (work array)
1065 : real(r8) precbl ! precip work variable
1066 : real(r8) precxx ! precip work variable
1067 : real(r8) precxx2 !
1068 : real(r8) precic ! precip work variable
1069 : real(r8) rat ! ratio of amount available to amount removed
1070 : real(r8) scavab(pcols) ! scavenged tracer flux from above (work array)
1071 : real(r8) scavabc(pcols) ! scavenged tracer flux from above (work array)
1072 :
1073 : real(r8) scavmax ! an estimate of the max tracer avail for removal
1074 : real(r8) scavbl ! flux removed at bottom of layer
1075 : real(r8) fins ! in cloud fraction removed by strat rain
1076 : real(r8) finc ! in cloud fraction removed by conv rain
1077 : real(r8) rate ! max removal rate estimate
1078 : real(r8) scavlimt ! limiting value 1
1079 : real(r8) scavt1 ! limiting value 2
1080 : real(r8) scavin ! scavenging by incloud processes
1081 : real(r8) scavbc ! scavenging by below cloud processes
1082 : real(r8) tc
1083 : real(r8) weight ! ice fraction
1084 : real(r8) wtpl ! work variable
1085 : real(r8) cldmabs(pcols) ! maximum cloud at or above this level
1086 : real(r8) cldmabc(pcols) ! maximum cloud at or above this level
1087 : !-----------------------------------------------------------
1088 :
1089 0 : adjfac = deltat/(max(deltat,cmftau)) ! adjustment factor from hack scheme
1090 :
1091 : ! zero accumulators
1092 0 : do i = 1,pcols
1093 0 : precab(i) = 1.e-36_r8
1094 0 : scavab(i) = 0._r8
1095 0 : cldmabs(i) = 0._r8
1096 : end do
1097 :
1098 0 : do k = 1,pver
1099 0 : do i = 1,ncol
1100 :
1101 0 : tc = t(i,k) - tmelt
1102 0 : weight = max(0._r8,min(-tc*0.05_r8,1.0_r8)) ! fraction of condensate that is ice
1103 :
1104 0 : cldmabs(i) = max(cldmabs(i),cldt(i,k))
1105 :
1106 : ! partitioning coefs for gas and aqueous phase
1107 : ! take as a cloud water amount, the sum of the stratiform amount
1108 : ! plus the convective rain water amount
1109 :
1110 : ! convective amnt is just the local precip rate from the hack scheme
1111 : ! since there is no storage of water, this ignores that falling from above
1112 : ! cwatc = cmfdqr(i,k)*deltat/adjfac
1113 : !++mcb -- test cwatc
1114 0 : cwatc = (icwmr1(i,k) + icwmr2(i,k)) * (1._r8-weight)
1115 : !--mcb
1116 :
1117 : ! strat cloud water amount and also ignore the part falling from above
1118 0 : cwats = cwat(i,k)
1119 :
1120 : ! cloud water as liq
1121 : !++mcb -- add cwatc later (in cwatti)
1122 : ! cwatl = (1.-weight)*(cwatc+cwats)
1123 0 : cwatl = (1._r8-weight)*cwats
1124 : ! cloud water as ice
1125 : !*not used cwati = weight*(cwatc+cwats)
1126 :
1127 : ! total suspended condensate as liquid
1128 0 : cwatt = cwatl + rain(i,k)
1129 :
1130 : ! incloud version
1131 : !++mcb -- add cwatc here
1132 0 : cwatti = cwatt/max(cldv(i,k), 0.00001_r8) + cwatc
1133 :
1134 : ! partitioning terms
1135 0 : patm = p(i,k)/1.013e5_r8 ! pressure in atmospheres
1136 0 : hconst = molwta*patm*solconst(i,k)*cwatti/rhoh2o
1137 0 : aqfrac = hconst/(1._r8+hconst)
1138 0 : gafrac = 1/(1._r8+hconst)
1139 0 : fracis(i,k) = gafrac
1140 :
1141 :
1142 : ! partial pressure of the tracer in the gridbox in atmospheres
1143 0 : part = patm*gafrac*tracer(i,k)*molwta/molwt
1144 :
1145 : ! use henrys law to give moles tracer /liter of water
1146 : ! in this volume
1147 : ! then convert to kg tracer /liter of water (kg tracer / kg water)
1148 0 : mplb = solconst(i,k)*part*molwt/1000._r8
1149 :
1150 :
1151 0 : pdog = pdel(i,k)/gravit
1152 :
1153 : ! this part of precip will be carried downward but at a new molarity of mpl
1154 0 : precic = pdog*(precs(i,k) + cmfdqr(i,k))
1155 :
1156 : ! we cant take out more than entered, plus that available in the cloud
1157 : ! scavmax = scavab(i)+tracer(i,k)*cldt(i,k)/deltat*pdog
1158 0 : scavmax = scavab(i)+tracer(i,k)*cldv(i,k)/deltat*pdog
1159 :
1160 : ! flux of tracer by incloud processes
1161 0 : scavin = precic*(1._r8-weight)*mplb
1162 :
1163 : ! fraction of precip which entered above that leaves below
1164 0 : if (cam_physpkg_is('cam5') .or. cam_physpkg_is('cam6')) then
1165 : ! Sungsu added evaporation of convective precipitation below.
1166 0 : precxx = precab(i)-pdog*(evaps(i,k)+evapc(i,k))
1167 : else
1168 0 : precxx = precab(i)-pdog*evaps(i,k)
1169 : end if
1170 0 : precxx = max (precxx,0.0_r8)
1171 :
1172 : ! flux of tracer by below cloud processes
1173 : !++mcb -- removed wtpl because it is now not assigned and previously
1174 : ! when it was assigned it was unnecessary: if(tc.gt.0)wtpl=1
1175 0 : if (tc>0.0_r8) then
1176 : ! scavbc = precxx*wtpl*mplb ! if liquid
1177 0 : scavbc = precxx*mplb ! if liquid
1178 : else
1179 0 : precxx2=max(precxx,1.e-36_r8)
1180 0 : scavbc = scavab(i)*precxx2/(precab(i)) ! if ice
1181 : endif
1182 :
1183 0 : scavbl = min(scavbc + scavin, scavmax)
1184 :
1185 : ! first guess assuming that henries law works
1186 0 : scavt1 = (scavab(i)-scavbl)/pdog*omsm
1187 :
1188 : ! pjr this should not be required, but we put it in to make sure we cant remove too much
1189 : ! remember, scavt1 is generally negative (indicating removal)
1190 0 : scavt1 = max(scavt1,-tracer(i,k)*cldv(i,k)/deltat)
1191 :
1192 : !++mcb -- remove this limitation for gas species
1193 : !c use the dana and hales or balkanski limit on scavenging
1194 : !c rate = precab(i)*0.1
1195 : ! rate = (precic + precxx)*0.1
1196 : ! scavlimt = -tracer(i,k)*cldv(i,k)
1197 : ! $ *rate/(1.+rate*deltat)
1198 :
1199 : ! scavt(i,k) = max(scavt1, scavlimt)
1200 :
1201 : ! instead just set scavt to scavt1
1202 0 : scavt(i,k) = scavt1
1203 : !--mcb
1204 :
1205 : ! now update the amount leaving the layer
1206 0 : scavbl = scavab(i) - scavt(i,k)*pdog
1207 :
1208 : ! in cloud amount is that formed locally over the total flux out bottom
1209 0 : fins = scavin/(scavin + scavbc + 1.e-36_r8)
1210 0 : iscavt(i,k) = scavt(i,k)*fins
1211 :
1212 0 : scavab(i) = scavbl
1213 0 : precab(i) = max(precxx + precic,1.e-36_r8)
1214 :
1215 :
1216 :
1217 : end do
1218 : end do
1219 :
1220 0 : end subroutine wetdepg
1221 :
1222 : !##############################################################################
1223 :
1224 0 : end module wetdep
|
