updated top-level README and version_decl for V4.5 (#1847)
[WRF.git] / chem / module_cam_mam_wetdep.F
blob6a02d2ffa4d1dd7c8faedd940dc91c10851ac421
1 ! Updated to CESM1.0.3 (CAM5.1.01) by Balwinder.Singh@pnnl.gov
2 #define WRF_PORT
3 #define MODAL_AERO
4 module wetdep
6 !----------------------------------------------------------------------- 
8 ! Wet deposition routines for both aerosols and gas phase constituents.
9
10 !-----------------------------------------------------------------------
12 use shr_kind_mod, only: r8 => shr_kind_r8
13 #ifndef WRF_PORT
14 use ppgrid,       only: pcols, pver
15 #endif
16 use physconst,    only: gravit, rair, tmelt
17 #ifndef WRF_PORT
18 use phys_control, only: cam_physpkg_is
19 use cam_logfile,  only: iulog
20 #else
21 use module_cam_support, only: pcols, pver, iulog
22 #endif
24 implicit none
25 save
26 private
28 public :: wetdepa_v1  ! scavenging codes for very soluble aerosols -- CAM4 version
29 public :: wetdepa_v2  ! scavenging codes for very soluble aerosols -- CAM5 version
30 #ifndef WRF_PORT
31 public :: wetdepg     ! scavenging of gas phase constituents by henry's law
32 #endif
33 public :: clddiag     ! calc of cloudy volume and rain mixing ratio
35 real(r8), parameter :: cmftau = 3600._r8
36 real(r8), parameter :: rhoh2o = 1000._r8            ! density of water
37 real(r8), parameter :: molwta = 28.97_r8            ! molecular weight dry air gm/mole
39 !==============================================================================
40 contains
41 !==============================================================================
43 subroutine clddiag(t, pmid, pdel, cmfdqr, evapc, &
44                    cldt, cldcu, cldst, cme, evapr, &
45                    prain, cldv, cldvcu, cldvst, rain, &
46                    ncol)
48    ! ------------------------------------------------------------------------------------ 
49    ! Estimate the cloudy volume which is occupied by rain or cloud water as
50    ! the max between the local cloud amount or the
51    ! sum above of (cloud*positive precip production)      sum total precip from above
52    !              ----------------------------------   x ------------------------
53    ! sum above of     (positive precip           )        sum positive precip from above
54    ! Author: P. Rasch
55    !         Sungsu Park. Mar.2010 
56    ! ------------------------------------------------------------------------------------
58    ! Input arguments:
59    real(r8), intent(in) :: t(pcols,pver)        ! temperature (K)
60    real(r8), intent(in) :: pmid(pcols,pver)     ! pressure at layer midpoints
61    real(r8), intent(in) :: pdel(pcols,pver)     ! pressure difference across layers
62    real(r8), intent(in) :: cmfdqr(pcols,pver)   ! dq/dt due to convective rainout 
63    real(r8), intent(in) :: evapc(pcols,pver)    ! Evaporation rate of convective precipitation ( >= 0 ) 
64    real(r8), intent(in) :: cldt(pcols,pver)    ! total cloud fraction
65    real(r8), intent(in) :: cldcu(pcols,pver)    ! Cumulus cloud fraction
66    real(r8), intent(in) :: cldst(pcols,pver)    ! Stratus cloud fraction
67    real(r8), intent(in) :: cme(pcols,pver)      ! rate of cond-evap within the cloud
68    real(r8), intent(in) :: evapr(pcols,pver)    ! rate of evaporation of falling precipitation (kg/kg/s)
69    real(r8), intent(in) :: prain(pcols,pver)    ! rate of conversion of condensate to precipitation (kg/kg/s)
70    integer, intent(in) :: ncol
72    ! Output arguments:
73    real(r8), intent(out) :: cldv(pcols,pver)     ! fraction occupied by rain or cloud water 
74    real(r8), intent(out) :: cldvcu(pcols,pver)   ! Convective precipitation volume
75    real(r8), intent(out) :: cldvst(pcols,pver)   ! Stratiform precipitation volume
76    real(r8), intent(out) :: rain(pcols,pver)     ! mixing ratio of rain (kg/kg)
78    ! Local variables:
79    integer  i, k
80    real(r8) convfw         ! used in fallspeed calculation; taken from findmcnew
81    real(r8) sumppr(pcols)        ! precipitation rate (kg/m2-s)
82    real(r8) sumpppr(pcols)       ! sum of positive precips from above
83    real(r8) cldv1(pcols)         ! precip weighted cloud fraction from above
84    real(r8) lprec                ! local production rate of precip (kg/m2/s)
85    real(r8) lprecp               ! local production rate of precip (kg/m2/s) if positive
86    real(r8) rho                  ! air density
87    real(r8) vfall
88    real(r8) sumppr_cu(pcols)     ! Convective precipitation rate (kg/m2-s)
89    real(r8) sumpppr_cu(pcols)    ! Sum of positive convective precips from above
90    real(r8) cldv1_cu(pcols)      ! Convective precip weighted convective cloud fraction from above
91    real(r8) lprec_cu             ! Local production rate of convective precip (kg/m2/s)
92    real(r8) lprecp_cu            ! Local production rate of convective precip (kg/m2/s) if positive
93    real(r8) sumppr_st(pcols)     ! Stratiform precipitation rate (kg/m2-s)
94    real(r8) sumpppr_st(pcols)    ! Sum of positive stratiform precips from above
95    real(r8) cldv1_st(pcols)      ! Stratiform precip weighted stratiform cloud fraction from above
96    real(r8) lprec_st             ! Local production rate of stratiform precip (kg/m2/s)
97    real(r8) lprecp_st            ! Local production rate of stratiform precip (kg/m2/s) if positive
98    ! -----------------------------------------------------------------------
100    convfw = 1.94_r8*2.13_r8*sqrt(rhoh2o*gravit*2.7e-4_r8)
101    do i=1,ncol
102       sumppr(i) = 0._r8
103       cldv1(i) = 0._r8
104       sumpppr(i) = 1.e-36_r8
105       sumppr_cu(i)  = 0._r8
106       cldv1_cu(i)   = 0._r8
107       sumpppr_cu(i) = 1.e-36_r8
108       sumppr_st(i)  = 0._r8
109       cldv1_st(i)   = 0._r8
110       sumpppr_st(i) = 1.e-36_r8
111    end do
113    do k = 1,pver
114       do i = 1,ncol
115          cldv(i,k) = &
116             max(min(1._r8, &
117             cldv1(i)/sumpppr(i) &
118             )*sumppr(i)/sumpppr(i), &
119             cldt(i,k) &
120             )
121          lprec = pdel(i,k)/gravit &
122             *(prain(i,k)+cmfdqr(i,k)-evapr(i,k))
123          lprecp = max(lprec,1.e-30_r8)
124          cldv1(i) = cldv1(i)  + cldt(i,k)*lprecp
125          sumppr(i) = sumppr(i) + lprec
126          sumpppr(i) = sumpppr(i) + lprecp
128          ! For convective precipitation volume at the top interface of each layer. Neglect the current layer.
129          cldvcu(i,k)   = max(min(1._r8,cldv1_cu(i)/sumpppr_cu(i))*(sumppr_cu(i)/sumpppr_cu(i)),0._r8)
130          lprec_cu      = (pdel(i,k)/gravit)*(cmfdqr(i,k)-evapc(i,k))
131          lprecp_cu     = max(lprec_cu,1.e-30_r8)
132          cldv1_cu(i)   = cldv1_cu(i) + cldcu(i,k)*lprecp_cu
133          sumppr_cu(i)  = sumppr_cu(i) + lprec_cu
134          sumpppr_cu(i) = sumpppr_cu(i) + lprecp_cu
136          ! For stratiform precipitation volume at the top interface of each layer. Neglect the current layer.
137          cldvst(i,k)   = max(min(1._r8,cldv1_st(i)/sumpppr_st(i))*(sumppr_st(i)/sumpppr_st(i)),0._r8)
138          lprec_st      = (pdel(i,k)/gravit)*(prain(i,k)-evapr(i,k))
139          lprecp_st     = max(lprec_st,1.e-30_r8)
140          cldv1_st(i)   = cldv1_st(i) + cldst(i,k)*lprecp_st
141          sumppr_st(i)  = sumppr_st(i) + lprec_st
142          sumpppr_st(i) = sumpppr_st(i) + lprecp_st
144          rain(i,k) = 0._r8
145          if(t(i,k) .gt. tmelt) then
146             rho = pmid(i,k)/(rair*t(i,k))
147             vfall = convfw/sqrt(rho)
148             rain(i,k) = sumppr(i)/(rho*vfall)
149             if (rain(i,k).lt.1.e-14_r8) rain(i,k) = 0._r8
150          endif
151       end do
152    end do
154 end subroutine clddiag
156 !==============================================================================
158 ! This is the CAM5 version of wetdepa.
160 subroutine wetdepa_v2(t, p, q, pdel, &
161                    cldt, cldc, cmfdqr, evapc, conicw, precs, conds, &
162                        evaps, cwat, tracer, deltat, &
163                        scavt, iscavt, cldv, cldvcu, cldvst, dlf, fracis, sol_fact, ncol, &
164                        scavcoef,  &
165 #ifdef MODAL_AERO
166                        is_strat_cloudborne, rate1ord_cw2pr_st, qqcw, f_act_conv, &     ! rce 2010/05/01
167 #endif
168                        icscavt, isscavt, bcscavt, bsscavt, &
169                        sol_facti_in, sol_factbi_in, sol_factii_in, &
170                        sol_factic_in, sol_factiic_in )
172       !----------------------------------------------------------------------- 
173       ! Purpose: 
174       ! scavenging code for very soluble aerosols
175       ! 
176       ! Author: P. Rasch
177       ! Modified by T. Bond 3/2003 to track different removals
178       ! Sungsu Park. Mar.2010 : Impose consistencies with a few changes in physics.
179       !-----------------------------------------------------------------------
183       implicit none
185       real(r8), intent(in) ::&
186          t(pcols,pver),        &! temperature
187          p(pcols,pver),        &! pressure
188          q(pcols,pver),        &! moisture
189          pdel(pcols,pver),     &! pressure thikness
190          cldt(pcols,pver),    &! total cloud fraction
191          cldc(pcols,pver),     &! convective cloud fraction
192          cmfdqr(pcols,pver),   &! rate of production of convective precip
193 ! Sungsu
194          evapc(pcols,pver),    &! Evaporation rate of convective precipitation
195 ! Sungsu
196          conicw(pcols,pver),   &! convective cloud water
197          cwat(pcols,pver),     &! cloud water amount 
198          precs(pcols,pver),    &! rate of production of stratiform precip
199          conds(pcols,pver),    &! rate of production of condensate
200          evaps(pcols,pver),    &! rate of evaporation of precip
201          cldv(pcols,pver),     &! total cloud fraction
202 ! Sungsu
203          cldvcu(pcols,pver),   &! Convective precipitation area at the top interface of each layer
204          cldvst(pcols,pver),   &! Stratiform precipitation area at the top interface of each layer
205          dlf(pcols,pver),      &! Detrainment of convective condensate [kg/kg/s]
206 ! Sungsu
207          deltat,               &! time step
208          tracer(pcols,pver)     ! trace species
209       ! If subroutine is called with just sol_fact:
210             ! sol_fact is used for both in- and below-cloud scavenging
211       ! If subroutine is called with optional argument sol_facti_in:
212             ! sol_fact  is used for below cloud scavenging
213             ! sol_facti is used for in cloud scavenging
214          real(r8), intent(in) :: sol_fact ! solubility factor (fraction of aer scavenged below & in, or just below or sol_facti_in is provided)
215          integer, intent(in) :: ncol
216          real(r8), intent(in) :: scavcoef(pcols,pver) ! Dana and Hales coefficient (/mm) (0.1 if not MODAL_AERO)
217 #ifdef MODAL_AERO
218          ! rce 2010/05/01
219          ! is_strat_cloudborne = .true. if tracer is stratiform-cloudborne aerosol; else .false. 
220          logical, intent(in), optional :: is_strat_cloudborne   
221          ! rate1ord_cw2pr_st = 1st order rate for strat cw to precip (1/s) 
222          real(r8), intent(in), optional  :: rate1ord_cw2pr_st(pcols,pver)
223          ! qqcw = strat-cloudborne aerosol corresponding to tracer when is_strat_cloudborne==.false.; else 0.0 
224          real(r8), intent(in), optional  :: qqcw(pcols,pver)
225          ! f_act_conv = conv-cloud activation fraction when is_strat_cloudborne==.false.; else 0.0 
226          real(r8), intent(in), optional  :: f_act_conv(pcols,pver)
227          ! end rce 2010/05/01
228 #endif
230          real(r8), intent(in), optional :: sol_facti_in   ! solubility factor (frac of aerosol scavenged in cloud)
231          real(r8), intent(in), optional :: sol_factbi_in  ! solubility factor (frac of aerosol scavenged below cloud by ice)
232          real(r8), intent(in), optional :: sol_factii_in  ! solubility factor (frac of aerosol scavenged in cloud by ice)
233          real(r8), intent(in), optional :: sol_factic_in(pcols,pver)  ! sol_facti_in for convective clouds
234          real(r8), intent(in), optional :: sol_factiic_in ! sol_factii_in for convective clouds
235          
236       real(r8), intent(out) ::&
237          scavt(pcols,pver),    &! scavenging tend 
238          iscavt(pcols,pver),   &! incloud scavenging tends
239          fracis(pcols,pver)     ! fraction of species not scavenged
241       real(r8), intent(out), optional ::    icscavt(pcols,pver)     ! incloud, convective
242       real(r8), intent(out), optional ::    isscavt(pcols,pver)     ! incloud, stratiform
243       real(r8), intent(out), optional ::    bcscavt(pcols,pver)     ! below cloud, convective
244       real(r8), intent(out), optional ::    bsscavt(pcols,pver)     ! below cloud, stratiform
246       ! local variables
248       integer i                 ! x index
249       integer k                 ! z index
251       real(r8) adjfac               ! factor stolen from cmfmca
252       real(r8) aqfrac               ! fraction of tracer in aqueous phase
253       real(r8) cwatc                ! local convective total water amount 
254       real(r8) cwats                ! local stratiform total water amount 
255       real(r8) cwatp                ! local water amount falling from above precip
256       real(r8) fracev(pcols)        ! fraction of precip from above that is evaporating
257 ! Sungsu
258       real(r8) fracev_cu(pcols)     ! Fraction of convective precip from above that is evaporating
259 ! Sungsu
260       real(r8) fracp                ! fraction of cloud water converted to precip
261       real(r8) gafrac               ! fraction of tracer in gas phasea
262       real(r8) hconst               ! henry's law solubility constant when equation is expressed
263                                 ! in terms of mixing ratios
264       real(r8) mpla                 ! moles / liter H2O entering the layer from above
265       real(r8) mplb                 ! moles / liter H2O leaving the layer below
266       real(r8) omsm                 ! 1 - (a small number)
267       real(r8) part                 !  partial pressure of tracer in atmospheres
268       real(r8) patm                 ! total pressure in atmospheres
269       real(r8) pdog                 ! work variable (pdel/gravit)
270       real(r8) precabc(pcols)       ! conv precip from above (work array)
271       real(r8) precabs(pcols)       ! strat precip from above (work array)
272       real(r8) precbl               ! precip falling out of level (work array)
273       real(r8) precmin              ! minimum convective precip causing scavenging
274       real(r8) rat(pcols)           ! ratio of amount available to amount removed
275       real(r8) scavab(pcols)        ! scavenged tracer flux from above (work array)
276       real(r8) scavabc(pcols)       ! scavenged tracer flux from above (work array)
277       real(r8) srcc                 ! tend for convective rain
278       real(r8) srcs                 ! tend for stratiform rain
279       real(r8) srct(pcols)          ! work variable
280       real(r8) tracab(pcols)        ! column integrated tracer amount
281 !      real(r8) vfall                ! fall speed of precip
282       real(r8) fins                 ! fraction of rem. rate by strat rain
283       real(r8) finc                 ! fraction of rem. rate by conv. rain
284       real(r8) srcs1                ! work variable
285       real(r8) srcs2                ! work variable
286       real(r8) tc                   ! temp in celcius
287       real(r8) weight               ! fraction of condensate which is ice
288       real(r8) cldmabs(pcols)       ! maximum cloud at or above this level
289       real(r8) cldmabc(pcols)       ! maximum cloud at or above this level
290       real(r8) odds                 ! limit on removal rate (proportional to prec)
291       real(r8) dblchek(pcols)
292       logical :: found
294     ! Jan.16.2009. Sungsu for wet scavenging below clouds.
295     ! real(r8) cldovr_cu(pcols)     ! Convective precipitation area at the base of each layer
296     ! real(r8) cldovr_st(pcols)     ! Stratiform precipitation area at the base of each layer
298       real(r8) tracer_incu
299       real(r8) tracer_mean
301     ! End by Sungsu
303       real(r8) sol_facti,  sol_factb  ! in cloud and below cloud fraction of aerosol scavenged
304       real(r8) sol_factii, sol_factbi ! in cloud and below cloud fraction of aerosol scavenged by ice
305       real(r8) sol_factic(pcols,pver)             ! sol_facti for convective clouds
306       real(r8) sol_factiic            ! sol_factii for convective clouds
307       ! sol_factic & solfact_iic added for MODAL_AERO.  
308       ! For stratiform cloud, cloudborne aerosol is treated explicitly,
309       !    and sol_facti is 1.0 for cloudborne, 0.0 for interstitial.
310       ! For convective cloud, cloudborne aerosol is not treated explicitly,
311       !    and sol_factic is 1.0 for both cloudborne and interstitial.
313       ! ------------------------------------------------------------------------
314 !      omsm = 1.-1.e-10          ! used to prevent roundoff errors below zero
315       omsm = 1._r8-2*epsilon(1._r8) ! used to prevent roundoff errors below zero
316       precmin =  0.1_r8/8.64e4_r8      ! set critical value to 0.1 mm/day in kg/m2/s
318       adjfac = deltat/(max(deltat,cmftau)) ! adjustment factor from hack scheme
320       ! assume 4 m/s fall speed currently (should be improved)
321 !      vfall = 4.
322         
323       ! default (if other sol_facts aren't in call, set all to required sol_fact
324       sol_facti = sol_fact
325       sol_factb = sol_fact
326       sol_factii = sol_fact
327       sol_factbi = sol_fact
329       if ( present(sol_facti_in) )  sol_facti = sol_facti_in
330       if ( present(sol_factii_in) )  sol_factii = sol_factii_in
331       if ( present(sol_factbi_in) )  sol_factbi = sol_factbi_in
333       sol_factic  = sol_facti
334       sol_factiic = sol_factii
335       if ( present(sol_factic_in ) )  sol_factic  = sol_factic_in
336       if ( present(sol_factiic_in) )  sol_factiic = sol_factiic_in
338       ! this section of code is for highly soluble aerosols,
339       ! the assumption is that within the cloud that
340       ! all the tracer is in the cloud water
341       !
342       ! for both convective and stratiform clouds, 
343       ! the fraction of cloud water converted to precip defines
344       ! the amount of tracer which is pulled out.
345       !
347       do i = 1,pcols
348          precabs(i) = 0
349          precabc(i) = 0
350          scavab(i) = 0
351          scavabc(i) = 0
352          tracab(i) = 0
353          cldmabs(i) = 0
354          cldmabc(i) = 0
356        ! Jan.16. Sungsu 
357        ! I added below to compute vertically projected cumulus and stratus fractions from the top to the
358        ! current model layer by assuming a simple independent maximum overlapping assumption for 
359        ! each cloud.
360        ! cldovr_cu(i) = 0._r8
361        ! cldovr_st(i) = 0._r8
362        ! End by Sungsu
364       end do
366       do k = 1,pver
367          do i = 1,ncol
368             tc     = t(i,k) - tmelt
369             weight = max(0._r8,min(-tc*0.05_r8,1.0_r8)) ! fraction of condensate that is ice
370             weight = 0._r8                                 ! assume no ice
372             pdog = pdel(i,k)/gravit
374             ! ****************** Evaporation **************************
375             ! calculate the fraction of strat precip from above 
376             !                 which evaporates within this layer
377             fracev(i) = evaps(i,k)*pdel(i,k)/gravit &
378                      /max(1.e-12_r8,precabs(i))
380             ! trap to ensure reasonable ratio bounds
381             fracev(i) = max(0._r8,min(1._r8,fracev(i)))
383 ! Sungsu : Same as above but convective precipitation part
384             fracev_cu(i) = evapc(i,k)*pdel(i,k)/gravit/max(1.e-12_r8,precabc(i))
385             fracev_cu(i) = max(0._r8,min(1._r8,fracev_cu(i)))
386 ! Sungsu
387             ! ****************** Convection ***************************
388             ! now do the convective scavenging
390             ! set odds proportional to fraction of the grid box that is swept by the 
391             ! precipitation =precabc/rhoh20*(area of sphere projected on plane
392             !                                /volume of sphere)*deltat
393             ! assume the radius of a raindrop is 1 e-3 m from Rogers and Yau,
394             ! unless the fraction of the area that is cloud is less than odds, in which
395             ! case use the cloud fraction (assumes precabs is in kg/m2/s)
396             ! is really: precabs*3/4/1000./1e-3*deltat
397             ! here I use .1 from Balkanski
398             !
399             ! use a local rate of convective rain production for incloud scav
400             !odds=max(min(1._r8, &
401             !     cmfdqr(i,k)*pdel(i,k)/gravit*0.1_r8*deltat),0._r8)
402             !++mcb -- change cldc to cldt; change cldt to cldv (9/17/96)
403             !            srcs1 =  cldt(i,k)*odds*tracer(i,k)*(1.-weight) &
404             ! srcs1 =  cldv(i,k)*odds*tracer(i,k)*(1.-weight) &
405             !srcs1 =  cldc(i,k)*odds*tracer(i,k)*(1.-weight) &
406             !         /deltat 
408             ! fraction of convective cloud water converted to rain
409             ! Dec.29.2009 : Sungsu multiplied cldc(i,k) to conicw(i,k) below
410             ! fracp = cmfdqr(i,k)*deltat/max(1.e-8_r8,conicw(i,k))
411             ! fracp = cmfdqr(i,k)*deltat/max(1.e-8_r8,cldc(i,k)*conicw(i,k))
412             ! Sungsu: Below new formula of 'fracp' is necessary since 'conicw' is a LWC/IWC
413             !         that has already precipitated out, that is, 'conicw' does not contain
414             !         precipitation at all ! 
415               fracp = cmfdqr(i,k)*deltat/max(1.e-12_r8,cldc(i,k)*conicw(i,k)+(cmfdqr(i,k)+dlf(i,k))*deltat) ! Sungsu.Mar.19.2010.
416             ! Dec.29.2009
417             ! Note cmfdrq can be negative from evap of rain, so constrain it <-- This is wrong. cmfdqr does not
418             ! contain evaporation of precipitation.
419             fracp = max(min(1._r8,fracp),0._r8)
420             ! remove that amount from within the convective area
421 !           srcs1 = cldc(i,k)*fracp*tracer(i,k)*(1._r8-weight)/deltat ! liquid only
422 !           srcs1 = cldc(i,k)*fracp*tracer(i,k)/deltat             ! any condensation
423 !           srcs1 = 0.
424 !           Jan.02.2010. Sungsu : cldt --> cldc below.
425 #ifdef MODAL_AERO
426             ! rce 2010/05/01
427             if ( is_strat_cloudborne ) then
428                ! only strat in-cloud removal affects strat-cloudborne aerosol
429                srcs1 = 0._r8
430             else
431                tracer_incu = f_act_conv(i,k)*(tracer(i,k)+& 
432                              min(qqcw(i,k),tracer(i,k)*((cldt(i,k)-cldc(i,k))/max(0.01_r8,(1._r8-(cldt(i,k)-cldc(i,k)))))))              
433                srcs1 = sol_factic(i,k)*cldc(i,k)*fracp*tracer_incu*(1._r8-weight)/deltat &  ! Liquid
434                      + sol_factiic    *cldc(i,k)*fracp*tracer_incu*(weight)/deltat          ! Ice
435             end if
436 #else
437             srcs1 = sol_factic(i,k)*cldc(i,k)*fracp*tracer(i,k)*(1._r8-weight)/deltat &  ! liquid
438                  +  sol_factiic*cldc(i,k)*fracp*tracer(i,k)*(weight)/deltat      ! ice
439 #endif
442             !--mcb
444             ! scavenge below cloud
446             !            cldmabc(i) = max(cldc(i,k),cldmabc(i))
447             !            cldmabc(i) = max(cldt(i,k),cldmabc(i))
448             ! cldmabc(i) = max(cldv(i,k),cldmabc(i))
449             ! cldmabc(i) = cldv(i,k)
450             cldmabc(i) = cldvcu(i,k)
452             ! Jan. 16. 2010. Sungsu
453             ! cldmabc(i) = cldmabc(i) * cldovr_cu(i) / max( 0.01_r8, cldovr_cu(i) + cldovr_st(i) )
454             ! End by Sungsu
456 #ifdef MODAL_AERO
457             ! rce 2010/05/01
458             if ( is_strat_cloudborne ) then
459                ! only strat in-cloud removal affects strat-cloudborne aerosol
460                srcs2 = 0._r8
461             else
462                tracer_mean = tracer(i,k)*(1._r8-cldc(i,k)*f_act_conv(i,k))-cldc(i,k)*f_act_conv(i,k)*&
463                              min(qqcw(i,k),tracer(i,k)*((cldt(i,k)-cldc(i,k))/max(0.01_r8,(1._r8-(cldt(i,k)-cldc(i,k))))))
464                tracer_mean = max(0._r8,tracer_mean) 
465                odds  = max(min(1._r8,precabc(i)/max(cldmabc(i),1.e-5_r8)*scavcoef(i,k)*deltat),0._r8) ! Dana and Hales coefficient (/mm)
466                srcs2 = sol_factb *cldmabc(i)*odds*tracer_mean*(1._r8-weight)/deltat & ! Liquid
467                      + sol_factbi*cldmabc(i)*odds*tracer_mean*(weight)/deltat         ! Ice
468             end if
469 #else
470             odds=max( &
471                  min(1._r8,precabc(i)/max(cldmabc(i),1.e-5_r8) &
472                  *scavcoef(i,k)*deltat),0._r8) ! Dana and Hales coefficient (/mm)
473             srcs2 = sol_factb*cldmabc(i)*odds*tracer(i,k)*(1._r8-weight)/deltat & ! liquid
474                  +  sol_factbi*cldmabc(i)*odds*tracer(i,k)*(weight)/deltat    !ice
475 #endif
478             !Note that using the temperature-determined weight doesn't make much sense here
481             srcc = srcs1 + srcs2  ! convective tend by both processes
482             finc = srcs1/(srcc + 1.e-36_r8) ! fraction in-cloud
484             ! ****************** Stratiform ***********************
485             ! now do the stratiform scavenging
487             ! incloud scavenging
489 #ifdef MODAL_AERO
490             ! rce 2010/05/01
491             if ( is_strat_cloudborne ) then
492                ! new code for stratiform incloud scav of cloudborne (modal) aerosol 
493                ! >> use the 1st order cw to precip rate calculated in microphysics routine
494                ! >> cloudborne aerosol resides in cloudy portion of grid cell, so do not apply "cldt" factor
495              ! fracp = rate1ord_cw2pr_st(i,k)*deltat
496              ! fracp = max(0._r8,min(1._r8,fracp))
497                fracp = precs(i,k)*deltat/max(cwat(i,k)+precs(i,k)*deltat,1.e-12_r8) ! Sungsu. Mar.19.2010.
498                fracp = max(0._r8,min(1._r8,fracp))
499                srcs1 = sol_facti *fracp*tracer(i,k)/deltat*(1._r8-weight) &  ! Liquid
500                      + sol_factii*fracp*tracer(i,k)/deltat*(weight)          ! Ice
501             else
502                ! strat in-cloud removal only affects strat-cloudborne aerosol
503                srcs1 = 0._r8
504             end if
505             ! end rce 2010/05/01
506 #else
508             ! fracp is the fraction of cloud water converted to precip
509             ! Sungsu modified fracp as the convectiv case.
510             !        Below new formula by Sungsu of 'fracp' is necessary since 'cwat' is a LWC/IWC
511             !        that has already precipitated out, that is, 'cwat' does not contain
512             !        precipitation at all ! 
513 !           fracp =  precs(i,k)*deltat/max(cwat(i,k),1.e-12_r8)
514             fracp =  precs(i,k)*deltat/max(cwat(i,k)+precs(i,k)*deltat,1.e-12_r8) ! Sungsu. Mar.19.2010.
515             fracp = max(0._r8,min(1._r8,fracp))
516 !           fracp = 0.     ! for debug
518             ! assume the corresponding amnt of tracer is removed
519             !++mcb -- remove cldc; change cldt to cldv 
520             !            srcs1 = (cldt(i,k)-cldc(i,k))*fracp*tracer(i,k)/deltat
521             !            srcs1 = cldv(i,k)*fracp*tracer(i,k)/deltat &
522 !           srcs1 = cldt(i,k)*fracp*tracer(i,k)/deltat            ! all condensate
523 !           Jan.02.2010. Sungsu : cldt --> cldt - cldc below.
524             srcs1 = sol_facti*(cldt(i,k)-cldc(i,k))*fracp*tracer(i,k)/deltat*(1._r8-weight) &  ! liquid
525                  + sol_factii*(cldt(i,k)-cldc(i,k))*fracp*tracer(i,k)/deltat*(weight)       ! ice
526 #endif
529             ! below cloud scavenging
531 !           volume undergoing below cloud scavenging
532 !           cldmabs(i) = cldv(i,k)   ! precipitating volume
533 !           cldmabs(i) = cldt(i,k)   ! local cloud volume
534             cldmabs(i) = cldvst(i,k) ! Stratiform precipitation area at the top interface of current layer
536             ! Jan. 16. 2010. Sungsu
537             ! cldmabs(i) = cldmabs(i) * cldovr_st(i) / max( 0.01_r8, cldovr_cu(i) + cldovr_st(i) )
538             ! End by Sungsu
540 #ifdef MODAL_AERO
541             ! rce 2010/05/01
542             if ( is_strat_cloudborne ) then
543                ! only strat in-cloud removal affects strat-cloudborne aerosol
544                srcs2 = 0._r8
545             else
546                odds = precabs(i)/max(cldmabs(i),1.e-5_r8)*scavcoef(i,k)*deltat
547                odds = max(min(1._r8,odds),0._r8)
548                srcs2 = sol_factb *cldmabs(i)*odds*tracer_mean*(1._r8-weight)/deltat & ! Liquid
549                      + sol_factbi*cldmabs(i)*odds*tracer_mean*(weight)/deltat         ! Ice
550             end if
551 #else
552             odds = precabs(i)/max(cldmabs(i),1.e-5_r8)*scavcoef(i,k)*deltat
553             odds = max(min(1._r8,odds),0._r8)
554             srcs2 =sol_factb*(cldmabs(i)*odds) *tracer(i,k)*(1._r8-weight)/deltat & ! liquid
555                  + sol_factbi*(cldmabs(i)*odds) *tracer(i,k)*(weight)/deltat       ! ice
558             !Note that using the temperature-determined weight doesn't make much sense here
559 #endif
562             srcs = srcs1 + srcs2             ! total stratiform scavenging
563             fins=srcs1/(srcs + 1.e-36_r8)    ! fraction taken by incloud processes
565             ! make sure we dont take out more than is there
566             ! ratio of amount available to amount removed
567             rat(i) = tracer(i,k)/max(deltat*(srcc+srcs),1.e-36_r8)
568             if (rat(i).lt.1._r8) then
569                srcs = srcs*rat(i)
570                srcc = srcc*rat(i)
571             endif
572             srct(i) = (srcc+srcs)*omsm
574             
575             ! fraction that is not removed within the cloud
576             ! (assumed to be interstitial, and subject to convective transport)
577             fracp = deltat*srct(i)/max(cldmabs(i)*tracer(i,k),1.e-36_r8)  ! amount removed
578             fracp = max(0._r8,min(1._r8,fracp))
579             fracis(i,k) = 1._r8 - fracp
581             ! tend is all tracer removed by scavenging, plus all re-appearing from evaporation above
582             ! Sungsu added cumulus contribution in the below 3 blocks
583          
584             scavt(i,k) = -srct(i) + (fracev(i)*scavab(i)+fracev_cu(i)*scavabc(i))*gravit/pdel(i,k)
585             iscavt(i,k) = -(srcc*finc + srcs*fins)*omsm
587             if ( present(icscavt) ) icscavt(i,k) = -(srcc*finc) * omsm
588             if ( present(isscavt) ) isscavt(i,k) = -(srcs*fins) * omsm
589             if ( present(bcscavt) ) bcscavt(i,k) = -(srcc * (1-finc)) * omsm +  &
590                  fracev_cu(i)*scavabc(i)*gravit/pdel(i,k)
591             if ( present(bsscavt) ) bsscavt(i,k) = -(srcs * (1-fins)) * omsm +  &
592                  fracev(i)*scavab(i)*gravit/pdel(i,k)
594             dblchek(i) = tracer(i,k) + deltat*scavt(i,k)
596             ! now keep track of scavenged mass and precip
597             scavab(i) = scavab(i)*(1-fracev(i)) + srcs*pdel(i,k)/gravit
598             precabs(i) = precabs(i) + (precs(i,k) - evaps(i,k))*pdel(i,k)/gravit
599             scavabc(i) = scavabc(i)*(1-fracev_cu(i)) + srcc*pdel(i,k)/gravit
600             precabc(i) = precabc(i) + (cmfdqr(i,k) - evapc(i,k))*pdel(i,k)/gravit
601             tracab(i) = tracab(i) + tracer(i,k)*pdel(i,k)/gravit
603        ! Jan.16.2010. Sungsu
604        ! Compute convective and stratiform precipitation areas at the base interface
605        ! of current layer. These are for computing 'below cloud scavenging' in the 
606        ! next layer below.
608        ! cldovr_cu(i) = max( cldovr_cu(i), cldc(i,k) )
609        ! cldovr_st(i) = max( cldovr_st(i), max( 0._r8, cldt(i,k) - cldc(i,k) ) )
611        ! cldovr_cu(i) = max( 0._r8, min ( 1._r8, cldovr_cu(i) ) )
612        ! cldovr_st(i) = max( 0._r8, min ( 1._r8, cldovr_st(i) ) )
614        ! End by Sungsu
616          end do ! End of i = 1, ncol
618          found = .false.
619          do i = 1,ncol
620             if ( dblchek(i) < 0._r8 ) then
621                found = .true.
622                exit
623             end if
624          end do
626          if ( found ) then
627             do i = 1,ncol
628                if (dblchek(i) .lt. 0._r8) then
629                   write(iulog,*) ' wetdapa: negative value ', i, k, tracer(i,k), &
630                        dblchek(i), scavt(i,k), srct(i), rat(i), fracev(i)
631 #ifdef WRF_PORT
632                   call wrf_message(iulog)
633 #endif 
634                endif
635             end do
636          endif
638       end do ! End of k = 1, pver
640    end subroutine wetdepa_v2
643 !==============================================================================
645 ! This is the frozen CAM4 version of wetdepa.
648    subroutine wetdepa_v1( t, p, q, pdel, &
649                        cldt, cldc, cmfdqr, conicw, precs, conds, &
650                        evaps, cwat, tracer, deltat, &
651                        scavt, iscavt, cldv, fracis, sol_fact, ncol, &
652                        scavcoef,icscavt, isscavt, bcscavt, bsscavt, &
653                        sol_facti_in, sol_factbi_in, sol_factii_in, &
654                        sol_factic_in, sol_factiic_in )
656       !----------------------------------------------------------------------- 
657       ! Purpose: 
658       ! scavenging code for very soluble aerosols
659       ! 
660       ! Author: P. Rasch
661       ! Modified by T. Bond 3/2003 to track different removals
662       !-----------------------------------------------------------------------
664       implicit none
666       real(r8), intent(in) ::&
667          t(pcols,pver),        &! temperature
668          p(pcols,pver),        &! pressure
669          q(pcols,pver),        &! moisture
670          pdel(pcols,pver),     &! pressure thikness
671          cldt(pcols,pver),    &! total cloud fraction
672          cldc(pcols,pver),     &! convective cloud fraction
673          cmfdqr(pcols,pver),   &! rate of production of convective precip
674          conicw(pcols,pver),   &! convective cloud water
675          cwat(pcols,pver),     &! cloud water amount 
676          precs(pcols,pver),    &! rate of production of stratiform precip
677          conds(pcols,pver),    &! rate of production of condensate
678          evaps(pcols,pver),    &! rate of evaporation of precip
679          cldv(pcols,pver),     &! total cloud fraction
680          deltat,               &! time step
681          tracer(pcols,pver)     ! trace species
682       ! If subroutine is called with just sol_fact:
683             ! sol_fact is used for both in- and below-cloud scavenging
684       ! If subroutine is called with optional argument sol_facti_in:
685             ! sol_fact  is used for below cloud scavenging
686             ! sol_facti is used for in cloud scavenging
687          real(r8), intent(in) :: sol_fact ! solubility factor (fraction of aer scavenged below & in, or just below or sol_facti_in is provided)
688          real(r8), intent(in), optional :: sol_facti_in   ! solubility factor (frac of aerosol scavenged in cloud)
689          real(r8), intent(in), optional :: sol_factbi_in  ! solubility factor (frac of aerosol scavenged below cloud by ice)
690          real(r8), intent(in), optional :: sol_factii_in  ! solubility factor (frac of aerosol scavenged in cloud by ice)
691          real(r8), intent(in), optional :: sol_factic_in(pcols,pver)  ! sol_facti_in for convective clouds
692          real(r8), intent(in), optional :: sol_factiic_in ! sol_factii_in for convective clouds
693          real(r8), intent(in) :: scavcoef(pcols,pver) ! Dana and Hales coefficient (/mm) (0.1 if not MODAL_AERO)
694          
695       integer, intent(in) :: ncol
697       real(r8), intent(out) ::&
698          scavt(pcols,pver),    &! scavenging tend 
699          iscavt(pcols,pver),   &! incloud scavenging tends
700          fracis(pcols,pver)     ! fraction of species not scavenged
702       real(r8), intent(out), optional ::    icscavt(pcols,pver)     ! incloud, convective
703       real(r8), intent(out), optional ::    isscavt(pcols,pver)     ! incloud, stratiform
704       real(r8), intent(out), optional ::    bcscavt(pcols,pver)     ! below cloud, convective
705       real(r8), intent(out), optional ::    bsscavt(pcols,pver)     ! below cloud, stratiform
707       ! local variables
709       integer i                 ! x index
710       integer k                 ! z index
712       real(r8) adjfac               ! factor stolen from cmfmca
713       real(r8) aqfrac               ! fraction of tracer in aqueous phase
714       real(r8) cwatc                ! local convective total water amount 
715       real(r8) cwats                ! local stratiform total water amount 
716       real(r8) cwatp                ! local water amount falling from above precip
717       real(r8) fracev(pcols)        ! fraction of precip from above that is evaporating
718       real(r8) fracp                ! fraction of cloud water converted to precip
719       real(r8) gafrac               ! fraction of tracer in gas phasea
720       real(r8) hconst               ! henry's law solubility constant when equation is expressed
721                                 ! in terms of mixing ratios
722       real(r8) mpla                 ! moles / liter H2O entering the layer from above
723       real(r8) mplb                 ! moles / liter H2O leaving the layer below
724       real(r8) omsm                 ! 1 - (a small number)
725       real(r8) part                 !  partial pressure of tracer in atmospheres
726       real(r8) patm                 ! total pressure in atmospheres
727       real(r8) pdog                 ! work variable (pdel/gravit)
728       real(r8) precabc(pcols)       ! conv precip from above (work array)
729       real(r8) precabs(pcols)       ! strat precip from above (work array)
730       real(r8) precbl               ! precip falling out of level (work array)
731       real(r8) precmin              ! minimum convective precip causing scavenging
732       real(r8) rat(pcols)           ! ratio of amount available to amount removed
733       real(r8) scavab(pcols)        ! scavenged tracer flux from above (work array)
734       real(r8) scavabc(pcols)       ! scavenged tracer flux from above (work array)
735       real(r8) srcc                 ! tend for convective rain
736       real(r8) srcs                 ! tend for stratiform rain
737       real(r8) srct(pcols)          ! work variable
738       real(r8) tracab(pcols)        ! column integrated tracer amount
739 !      real(r8) vfall                ! fall speed of precip
740       real(r8) fins                 ! fraction of rem. rate by strat rain
741       real(r8) finc                 ! fraction of rem. rate by conv. rain
742       real(r8) srcs1                ! work variable
743       real(r8) srcs2                ! work variable
744       real(r8) tc                   ! temp in celcius
745       real(r8) weight               ! fraction of condensate which is ice
746       real(r8) cldmabs(pcols)       ! maximum cloud at or above this level
747       real(r8) cldmabc(pcols)       ! maximum cloud at or above this level
748       real(r8) odds                 ! limit on removal rate (proportional to prec)
749       real(r8) dblchek(pcols)
750       logical :: found
752       real(r8) sol_facti,  sol_factb  ! in cloud and below cloud fraction of aerosol scavenged
753       real(r8) sol_factii, sol_factbi ! in cloud and below cloud fraction of aerosol scavenged by ice
754       real(r8) sol_factic(pcols,pver)             ! sol_facti for convective clouds
755       real(r8) sol_factiic            ! sol_factii for convective clouds
756       ! sol_factic & solfact_iic added for MODAL_AERO.  
757       ! For stratiform cloud, cloudborne aerosol is treated explicitly,
758       !    and sol_facti is 1.0 for cloudborne, 0.0 for interstitial.
759       ! For convective cloud, cloudborne aerosol is not treated explicitly,
760       !    and sol_factic is 1.0 for both cloudborne and interstitial.
762       ! ------------------------------------------------------------------------
763 !      omsm = 1.-1.e-10          ! used to prevent roundoff errors below zero
764       omsm = 1._r8-2*epsilon(1._r8) ! used to prevent roundoff errors below zero
765       precmin =  0.1_r8/8.64e4_r8      ! set critical value to 0.1 mm/day in kg/m2/s
767       adjfac = deltat/(max(deltat,cmftau)) ! adjustment factor from hack scheme
769       ! assume 4 m/s fall speed currently (should be improved)
770 !      vfall = 4.
771         
772       ! default (if other sol_facts aren't in call, set all to required sol_fact
773       sol_facti = sol_fact
774       sol_factb = sol_fact
775       sol_factii = sol_fact
776       sol_factbi = sol_fact
778       if ( present(sol_facti_in) )  sol_facti = sol_facti_in
779       if ( present(sol_factii_in) )  sol_factii = sol_factii_in
780       if ( present(sol_factbi_in) )  sol_factbi = sol_factbi_in
782       sol_factic  = sol_facti
783       sol_factiic = sol_factii
784       if ( present(sol_factic_in ) )  sol_factic  = sol_factic_in
785       if ( present(sol_factiic_in) )  sol_factiic = sol_factiic_in
787       ! this section of code is for highly soluble aerosols,
788       ! the assumption is that within the cloud that
789       ! all the tracer is in the cloud water
790       !
791       ! for both convective and stratiform clouds, 
792       ! the fraction of cloud water converted to precip defines
793       ! the amount of tracer which is pulled out.
794       !
796       do i = 1,pcols
797          precabs(i) = 0
798          precabc(i) = 0
799          scavab(i) = 0
800          scavabc(i) = 0
801          tracab(i) = 0
802          cldmabs(i) = 0
803          cldmabc(i) = 0
804       end do
806       do k = 1,pver
807          do i = 1,ncol
808             tc     = t(i,k) - tmelt
809             weight = max(0._r8,min(-tc*0.05_r8,1.0_r8)) ! fraction of condensate that is ice
810             weight = 0._r8                                 ! assume no ice
812             pdog = pdel(i,k)/gravit
814             ! ****************** Evaporation **************************
815             ! calculate the fraction of strat precip from above 
816             !                 which evaporates within this layer
817             fracev(i) = evaps(i,k)*pdel(i,k)/gravit &
818                      /max(1.e-12_r8,precabs(i))
820             ! trap to ensure reasonable ratio bounds
821             fracev(i) = max(0._r8,min(1._r8,fracev(i)))
823             ! ****************** Convection ***************************
824             ! now do the convective scavenging
826             ! set odds proportional to fraction of the grid box that is swept by the 
827             ! precipitation =precabc/rhoh20*(area of sphere projected on plane
828             !                                /volume of sphere)*deltat
829             ! assume the radius of a raindrop is 1 e-3 m from Rogers and Yau,
830             ! unless the fraction of the area that is cloud is less than odds, in which
831             ! case use the cloud fraction (assumes precabs is in kg/m2/s)
832             ! is really: precabs*3/4/1000./1e-3*deltat
833             ! here I use .1 from Balkanski
834             !
835             ! use a local rate of convective rain production for incloud scav
836             !odds=max(min(1._r8, &
837             !     cmfdqr(i,k)*pdel(i,k)/gravit*0.1_r8*deltat),0._r8)
838             !++mcb -- change cldc to cldt; change cldt to cldv (9/17/96)
839             !            srcs1 =  cldt(i,k)*odds*tracer(i,k)*(1.-weight) &
840             ! srcs1 =  cldv(i,k)*odds*tracer(i,k)*(1.-weight) &
841             !srcs1 =  cldc(i,k)*odds*tracer(i,k)*(1.-weight) &
842             !         /deltat 
844             ! fraction of convective cloud water converted to rain
845             fracp = cmfdqr(i,k)*deltat/max(1.e-8_r8,conicw(i,k))
846             ! note cmfdrq can be negative from evap of rain, so constrain it
847             fracp = max(min(1._r8,fracp),0._r8)
848             ! remove that amount from within the convective area
849 !           srcs1 = cldc(i,k)*fracp*tracer(i,k)*(1._r8-weight)/deltat ! liquid only
850 !           srcs1 = cldc(i,k)*fracp*tracer(i,k)/deltat             ! any condensation
851 !           srcs1 = 0.
852             srcs1 = sol_factic(i,k)*cldt(i,k)*fracp*tracer(i,k)*(1._r8-weight)/deltat &  ! liquid
853                  +  sol_factiic*cldt(i,k)*fracp*tracer(i,k)*(weight)/deltat      ! ice
856             !--mcb
858             ! scavenge below cloud
860             !            cldmabc(i) = max(cldc(i,k),cldmabc(i))
861             !            cldmabc(i) = max(cldt(i,k),cldmabc(i))
862             cldmabc(i) = max(cldv(i,k),cldmabc(i))
863             cldmabc(i) = cldv(i,k)
865             odds=max( &
866                  min(1._r8,precabc(i)/max(cldmabc(i),1.e-5_r8) &
867                  *scavcoef(i,k)*deltat),0._r8) ! Dana and Hales coefficient (/mm)
868             srcs2 = sol_factb*cldmabc(i)*odds*tracer(i,k)*(1._r8-weight)/deltat & ! liquid
869                  +  sol_factbi*cldmabc(i)*odds*tracer(i,k)*(weight)/deltat    !ice
870             !Note that using the temperature-determined weight doesn't make much sense here
873             srcc = srcs1 + srcs2  ! convective tend by both processes
874             finc = srcs1/(srcc + 1.e-36_r8) ! fraction in-cloud
876             ! ****************** Stratiform ***********************
877             ! now do the stratiform scavenging
879             ! incloud scavenging
881             ! fracp is the fraction of cloud water converted to precip
882             fracp =  precs(i,k)*deltat/max(cwat(i,k),1.e-12_r8)
883             fracp = max(0._r8,min(1._r8,fracp))
884 !           fracp = 0.     ! for debug
886             ! assume the corresponding amnt of tracer is removed
887             !++mcb -- remove cldc; change cldt to cldv 
888             !            srcs1 = (cldt(i,k)-cldc(i,k))*fracp*tracer(i,k)/deltat
889             !            srcs1 = cldv(i,k)*fracp*tracer(i,k)/deltat &
890 !           srcs1 = cldt(i,k)*fracp*tracer(i,k)/deltat            ! all condensate
891             srcs1 = sol_facti*cldt(i,k)*fracp*tracer(i,k)/deltat*(1._r8-weight) &  ! liquid
892                  + sol_factii*cldt(i,k)*fracp*tracer(i,k)/deltat*(weight)       ! ice
895             ! below cloud scavenging
897 !           volume undergoing below cloud scavenging
898             cldmabs(i) = cldv(i,k)   ! precipitating volume
899 !           cldmabs(i) = cldt(i,k)   ! local cloud volume
901             odds = precabs(i)/max(cldmabs(i),1.e-5_r8)*scavcoef(i,k)*deltat
902             odds = max(min(1._r8,odds),0._r8)
903             srcs2 =sol_factb*(cldmabs(i)*odds) *tracer(i,k)*(1._r8-weight)/deltat & ! liquid
904                  + sol_factbi*(cldmabs(i)*odds) *tracer(i,k)*(weight)/deltat       ! ice
905             !Note that using the temperature-determined weight doesn't make much sense here
908             srcs = srcs1 + srcs2             ! total stratiform scavenging
909             fins=srcs1/(srcs + 1.e-36_r8)    ! fraction taken by incloud processes
911             ! make sure we dont take out more than is there
912             ! ratio of amount available to amount removed
913             rat(i) = tracer(i,k)/max(deltat*(srcc+srcs),1.e-36_r8)
914             if (rat(i).lt.1._r8) then
915                srcs = srcs*rat(i)
916                srcc = srcc*rat(i)
917             endif
918             srct(i) = (srcc+srcs)*omsm
920             
921             ! fraction that is not removed within the cloud
922             ! (assumed to be interstitial, and subject to convective transport)
923             fracp = deltat*srct(i)/max(cldmabs(i)*tracer(i,k),1.e-36_r8)  ! amount removed
924             fracp = max(0._r8,min(1._r8,fracp))
925             fracis(i,k) = 1._r8 - fracp
927             ! tend is all tracer removed by scavenging, plus all re-appearing from evaporation above
928             scavt(i,k) = -srct(i) + fracev(i)*scavab(i)*gravit/pdel(i,k)
929             iscavt(i,k) = -(srcc*finc + srcs*fins)*omsm
931             if ( present(icscavt) ) icscavt(i,k) = -(srcc*finc) * omsm
932             if ( present(isscavt) ) isscavt(i,k) = -(srcs*fins) * omsm
933             if ( present(bcscavt) ) bcscavt(i,k) = -(srcc * (1-finc)) * omsm
934             if ( present(bsscavt) ) bsscavt(i,k) = -(srcs * (1-fins)) * omsm +  &
935                  fracev(i)*scavab(i)*gravit/pdel(i,k)
937             dblchek(i) = tracer(i,k) + deltat*scavt(i,k)
939             ! now keep track of scavenged mass and precip
940             scavab(i) = scavab(i)*(1-fracev(i)) + srcs*pdel(i,k)/gravit
941             precabs(i) = precabs(i) + (precs(i,k) - evaps(i,k))*pdel(i,k)/gravit
942             scavabc(i) = scavabc(i) + srcc*pdel(i,k)/gravit
943             precabc(i) = precabc(i) + (cmfdqr(i,k))*pdel(i,k)/gravit
944             tracab(i) = tracab(i) + tracer(i,k)*pdel(i,k)/gravit
946          end do
948          found = .false.
949          do i = 1,ncol
950             if ( dblchek(i) < 0._r8 ) then
951                found = .true.
952                exit
953             end if
954          end do
956          if ( found ) then
957             do i = 1,ncol
958                if (dblchek(i) .lt. 0._r8) then
959                   write(iulog,*) ' wetdapa: negative value ', i, k, tracer(i,k), &
960                        dblchek(i), scavt(i,k), srct(i), rat(i), fracev(i)
961 #ifdef WRF_PORT
962                   call wrf_message(iulog)
963 #endif 
964                endif
965             end do
966          endif
968       end do
970    end subroutine wetdepa_v1
971 #ifndef WRF_PORT
972 !==============================================================================
974 ! wetdepg is currently being used for both CAM4 and CAM5 by making use of the
975 ! cam_physpkg_is method.
977    subroutine wetdepg( t, p, q, pdel, &
978                        cldt, cldc, cmfdqr, evapc, precs, evaps, &
979                        rain, cwat, tracer, deltat, molwt, &
980                        solconst, scavt, iscavt, cldv, icwmr1, &
981                        icwmr2, fracis, ncol )
983       !----------------------------------------------------------------------- 
984       ! Purpose: 
985       ! scavenging of gas phase constituents by henry's law
986       ! 
987       ! Author: P. Rasch
988       !-----------------------------------------------------------------------
990       real(r8), intent(in) ::&
991          t(pcols,pver),        &! temperature
992          p(pcols,pver),        &! pressure
993          q(pcols,pver),        &! moisture
994          pdel(pcols,pver),     &! pressure thikness
995          cldt(pcols,pver),     &! total cloud fraction
996          cldc(pcols,pver),     &! convective cloud fraction
997          cmfdqr(pcols,pver),   &! rate of production of convective precip
998          rain (pcols,pver),    &! total rainwater mixing ratio
999          cwat(pcols,pver),     &! cloud water amount 
1000          precs(pcols,pver),    &! rate of production of stratiform precip
1001          evaps(pcols,pver),    &! rate of evaporation of precip
1002 ! Sungsu
1003          evapc(pcols,pver),    &! Rate of evaporation of convective precipitation
1004 ! Sungsu 
1005          cldv(pcols,pver),     &! estimate of local volume occupied by clouds
1006          icwmr1 (pcols,pver),  &! in cloud water mixing ration for zhang scheme
1007          icwmr2 (pcols,pver),  &! in cloud water mixing ration for hack  scheme
1008          deltat,               &! time step
1009          tracer(pcols,pver),   &! trace species
1010          molwt                  ! molecular weights
1012       integer, intent(in) :: ncol
1014       real(r8) &
1015          solconst(pcols,pver)   ! Henry's law coefficient
1017       real(r8), intent(out) ::&
1018          scavt(pcols,pver),    &! scavenging tend 
1019          iscavt(pcols,pver),   &! incloud scavenging tends
1020          fracis(pcols, pver)    ! fraction of constituent that is insoluble
1022       ! local variables
1024       integer i                 ! x index
1025       integer k                 ! z index
1027       real(r8) adjfac               ! factor stolen from cmfmca
1028       real(r8) aqfrac               ! fraction of tracer in aqueous phase
1029       real(r8) cwatc                ! local convective total water amount 
1030       real(r8) cwats                ! local stratiform total water amount 
1031       real(r8) cwatl                ! local cloud liq water amount 
1032       real(r8) cwatp                ! local water amount falling from above precip
1033       real(r8) cwatpl               ! local water amount falling from above precip (liq)
1034       real(r8) cwatt                ! local sum of strat + conv total water amount 
1035       real(r8) cwatti               ! cwatt/cldv = cloudy grid volume mixing ratio
1036       real(r8) fracev               ! fraction of precip from above that is evaporating
1037       real(r8) fracp                ! fraction of cloud water converted to precip
1038       real(r8) gafrac               ! fraction of tracer in gas phasea
1039       real(r8) hconst               ! henry's law solubility constant when equation is expressed
1040                                 ! in terms of mixing ratios
1041       real(r8) mpla                 ! moles / liter H2O entering the layer from above
1042       real(r8) mplb                 ! moles / liter H2O leaving the layer below
1043       real(r8) omsm                 ! 1 - (a small number)
1044       real(r8) part                 !  partial pressure of tracer in atmospheres
1045       real(r8) patm                 ! total pressure in atmospheres
1046       real(r8) pdog                 ! work variable (pdel/gravit)
1047       real(r8) precab(pcols)        ! precip from above (work array)
1048       real(r8) precbl               ! precip work variable
1049       real(r8) precxx               ! precip work variable
1050       real(r8) precxx2               !
1051       real(r8) precic               ! precip work variable
1052       real(r8) rat                  ! ratio of amount available to amount removed
1053       real(r8) scavab(pcols)        ! scavenged tracer flux from above (work array)
1054       real(r8) scavabc(pcols)       ! scavenged tracer flux from above (work array)
1055       !      real(r8) vfall                ! fall speed of precip
1056       real(r8) scavmax              ! an estimate of the max tracer avail for removal
1057       real(r8) scavbl               ! flux removed at bottom of layer
1058       real(r8) fins                 ! in cloud fraction removed by strat rain
1059       real(r8) finc                 ! in cloud fraction removed by conv rain
1060       real(r8) rate                 ! max removal rate estimate
1061       real(r8) scavlimt             ! limiting value 1
1062       real(r8) scavt1               ! limiting value 2
1063       real(r8) scavin               ! scavenging by incloud processes
1064       real(r8) scavbc               ! scavenging by below cloud processes
1065       real(r8) tc
1066       real(r8) weight               ! ice fraction
1067       real(r8) wtpl                 ! work variable
1068       real(r8) cldmabs(pcols)       ! maximum cloud at or above this level
1069       real(r8) cldmabc(pcols)       ! maximum cloud at or above this level
1070       !-----------------------------------------------------------
1072       omsm = 1._r8-2*epsilon(1._r8)   ! used to prevent roundoff errors below zero
1074       adjfac = deltat/(max(deltat,cmftau)) ! adjustment factor from hack scheme
1076       ! assume 4 m/s fall speed currently (should be improved)
1077       !      vfall = 4.
1079       ! zero accumulators
1080       do i = 1,pcols
1081          precab(i) = 1.e-36_r8
1082          scavab(i) = 0._r8
1083          cldmabs(i) = 0._r8
1084       end do
1086       do k = 1,pver
1087          do i = 1,ncol
1089             tc     = t(i,k) - tmelt
1090             weight = max(0._r8,min(-tc*0.05_r8,1.0_r8)) ! fraction of condensate that is ice
1092             cldmabs(i) = max(cldmabs(i),cldt(i,k))
1094             ! partitioning coefs for gas and aqueous phase
1095             !              take as a cloud water amount, the sum of the stratiform amount
1096             !              plus the convective rain water amount 
1098             ! convective amnt is just the local precip rate from the hack scheme
1099             !              since there is no storage of water, this ignores that falling from above
1100             !            cwatc = cmfdqr(i,k)*deltat/adjfac
1101             !++mcb -- test cwatc
1102             cwatc = (icwmr1(i,k) + icwmr2(i,k)) * (1._r8-weight)
1103             !--mcb 
1105             ! strat cloud water amount and also ignore the part falling from above
1106             cwats = cwat(i,k) 
1108             ! cloud water as liq
1109             !++mcb -- add cwatc later (in cwatti)
1110             !            cwatl = (1.-weight)*(cwatc+cwats)
1111             cwatl = (1._r8-weight)*cwats
1112             ! cloud water as ice
1113             !*not used        cwati = weight*(cwatc+cwats)
1115             ! total suspended condensate as liquid
1116             cwatt = cwatl + rain(i,k)
1118             ! incloud version 
1119             !++mcb -- add cwatc here
1120             cwatti = cwatt/max(cldv(i,k), 0.00001_r8) + cwatc
1122             ! partitioning terms
1123             patm = p(i,k)/1.013e5_r8 ! pressure in atmospheres
1124             hconst = molwta*patm*solconst(i,k)*cwatti/rhoh2o
1125             aqfrac = hconst/(1._r8+hconst)
1126             gafrac = 1/(1._r8+hconst)
1127             fracis(i,k) = gafrac
1130             ! partial pressure of the tracer in the gridbox in atmospheres
1131             part = patm*gafrac*tracer(i,k)*molwta/molwt
1133             ! use henrys law to give moles tracer /liter of water
1134             ! in this volume 
1135             ! then convert to kg tracer /liter of water (kg tracer / kg water)
1136             mplb = solconst(i,k)*part*molwt/1000._r8
1139             pdog = pdel(i,k)/gravit
1141             ! this part of precip will be carried downward but at a new molarity of mpl 
1142             precic = pdog*(precs(i,k) + cmfdqr(i,k))
1144             ! we cant take out more than entered, plus that available in the cloud
1145             !                  scavmax = scavab(i)+tracer(i,k)*cldt(i,k)/deltat*pdog
1146             scavmax = scavab(i)+tracer(i,k)*cldv(i,k)/deltat*pdog
1148             ! flux of tracer by incloud processes
1149             scavin = precic*(1._r8-weight)*mplb
1151             ! fraction of precip which entered above that leaves below
1152             if (cam_physpkg_is('cam5')) then
1153                ! Sungsu added evaporation of convective precipitation below.
1154                precxx = precab(i)-pdog*(evaps(i,k)+evapc(i,k))
1155             else
1156                precxx = precab(i)-pdog*evaps(i,k)
1157             end if
1158             precxx = max (precxx,0.0_r8)
1160             ! flux of tracer by below cloud processes
1161             !++mcb -- removed wtpl because it is now not assigned and previously
1162             !          when it was assigned it was unnecessary:  if(tc.gt.0)wtpl=1
1163             if (tc.gt.0) then
1164                !               scavbc = precxx*wtpl*mplb ! if liquid
1165                scavbc = precxx*mplb ! if liquid
1166             else
1167                precxx2=max(precxx,1.e-36_r8)
1168                scavbc = scavab(i)*precxx2/(precab(i)) ! if ice
1169             endif
1171             scavbl = min(scavbc + scavin, scavmax)
1173             ! first guess assuming that henries law works
1174             scavt1 = (scavab(i)-scavbl)/pdog*omsm
1176             ! pjr this should not be required, but we put it in to make sure we cant remove too much
1177             ! remember, scavt1 is generally negative (indicating removal)
1178             scavt1 = max(scavt1,-tracer(i,k)*cldv(i,k)/deltat)
1180             !++mcb -- remove this limitation for gas species
1181             !c use the dana and hales or balkanski limit on scavenging
1182             !c            rate = precab(i)*0.1
1183             !            rate = (precic + precxx)*0.1
1184             !            scavlimt = -tracer(i,k)*cldv(i,k)
1185             !     $           *rate/(1.+rate*deltat)
1187             !            scavt(i,k) = max(scavt1, scavlimt)
1189             ! instead just set scavt to scavt1
1190             scavt(i,k) = scavt1
1191             !--mcb
1193             ! now update the amount leaving the layer
1194             scavbl = scavab(i) - scavt(i,k)*pdog 
1196             ! in cloud amount is that formed locally over the total flux out bottom
1197             fins = scavin/(scavin + scavbc + 1.e-36_r8)
1198             iscavt(i,k) = scavt(i,k)*fins
1200             scavab(i) = scavbl
1201             precab(i) = max(precxx + precic,1.e-36_r8)
1203         
1204             
1205          end do
1206       end do
1207       
1208    end subroutine wetdepg
1209 #endif
1210 !##############################################################################
1212 end module wetdep