ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions/Subroutines | Variables
mp_wsm6 Module Reference

Functions/Subroutines

subroutine, public mp_wsm6_init (den0, denr, dens, cl, cpv, hail_opt, errmsg, errflg)
 
subroutine, public mp_wsm6_finalize (errmsg, errflg)
 
subroutine, public mp_wsm6_run (t, q, qc, qi, qr, qs, qg, den, p, delz, delt, g, cpd, cpv, rd, rv, t0c, ep1, ep2, qmin, xls, xlv0, xlf0, den0, denr, cliq, cice, psat, rain, rainncv, sr, snow, snowncv, graupel, graupelncv, rainprod2d, evapprod2d, its, ite, kts, kte, microphysics_debug, diag_i_dbg, diag_j_dbg, diag_k_raw_base, errmsg, errflg)
 
real(kind=kind_phys) function rgmma (x)
 
real(kind=kind_phys) function fpvs (t, ice, rd, rv, cvap, cliq, cice, hvap, hsub, psat, t0c)
 
subroutine slope_wsm6 (qrs, den, denfac, t, rslope, rslopeb, rslope2, rslope3, vt, its, ite, kts, kte)
 
subroutine slope_rain (qrs, den, denfac, t, rslope, rslopeb, rslope2, rslope3, vt, its, ite, kts, kte)
 
subroutine slope_snow (qrs, den, denfac, t, rslope, rslopeb, rslope2, rslope3, vt, its, ite, kts, kte)
 
subroutine slope_graup (qrs, den, denfac, t, rslope, rslopeb, rslope2, rslope3, vt, its, ite, kts, kte)
 
subroutine nislfv_rain_plm (im, km, denl, denfacl, tkl, dzl, wwl, rql, precip, dt, id, iter)
 
subroutine nislfv_rain_plm6 (im, km, denl, denfacl, tkl, dzl, wwl, rql, rql2, precip1, precip2, dt, id, iter)
 
subroutine, public refl10cm_wsm6 (qv1d, qr1d, qs1d, qg1d, t1d, p1d, dBZ, kts, kte)
 

Variables

integer, parameter kind_phys = c_double
 
real(kind=kind_phys), parameter, private dtcldcr = 120.
 
real(kind=kind_phys), parameter, private n0r = 8.e6
 
real(kind=kind_phys), parameter, private avtr = 841.9
 
real(kind=kind_phys), parameter, private bvtr = 0.8
 
real(kind=kind_phys), parameter, private r0 = .8e-5
 
real(kind=kind_phys), parameter, private peaut = .55
 
real(kind=kind_phys), parameter, private xncr = 3.e8
 
real(kind=kind_phys), parameter, private xmyu = 1.718e-5
 
real(kind=kind_phys), parameter, private avts = 11.72
 
real(kind=kind_phys), parameter, private bvts = .41
 
real(kind=kind_phys), parameter, private lamdarmax = 8.e4
 
real(kind=kind_phys), parameter, private lamdasmax = 1.e5
 
real(kind=kind_phys), parameter, private dicon = 11.9
 
real(kind=kind_phys), parameter, private dimax = 500.e-6
 
real(kind=kind_phys), parameter, private pfrz1 = 100.
 
real(kind=kind_phys), parameter, private pfrz2 = 0.66
 
real(kind=kind_phys), parameter, private qcrmin = 1.e-9
 
real(kind=kind_phys), parameter, private eacrc = 1.0
 
real(kind=kind_phys), parameter, private dens = 100.0
 
real(kind=kind_phys), parameter, private qs0 = 6.e-4
 
real(kind=kind_phys), parameter, public n0smax = 1.e11
 
real(kind=kind_phys), parameter, public n0s = 2.e6
 
real(kind=kind_phys), parameter, public alpha = .12
 
real(kind=kind_phys), save qc0
 
real(kind=kind_phys), save qck1
 
real(kind=kind_phys), save bvtr1
 
real(kind=kind_phys), save bvtr2
 
real(kind=kind_phys), save bvtr3
 
real(kind=kind_phys), save bvtr4
 
real(kind=kind_phys), save g1pbr
 
real(kind=kind_phys), save g3pbr
 
real(kind=kind_phys), save g4pbr
 
real(kind=kind_phys), save g5pbro2
 
real(kind=kind_phys), save pvtr
 
real(kind=kind_phys), save eacrr
 
real(kind=kind_phys), save pacrr
 
real(kind=kind_phys), save bvtr6
 
real(kind=kind_phys), save g6pbr
 
real(kind=kind_phys), save precr1
 
real(kind=kind_phys), save precr2
 
real(kind=kind_phys), save roqimax
 
real(kind=kind_phys), save bvts1
 
real(kind=kind_phys), save bvts2
 
real(kind=kind_phys), save bvts3
 
real(kind=kind_phys), save bvts4
 
real(kind=kind_phys), save g1pbs
 
real(kind=kind_phys), save g3pbs
 
real(kind=kind_phys), save g4pbs
 
real(kind=kind_phys), save n0g
 
real(kind=kind_phys), save avtg
 
real(kind=kind_phys), save bvtg
 
real(kind=kind_phys), save deng
 
real(kind=kind_phys), save lamdagmax
 
real(kind=kind_phys), save g5pbso2
 
real(kind=kind_phys), save pvts
 
real(kind=kind_phys), save pacrs
 
real(kind=kind_phys), save precs1
 
real(kind=kind_phys), save precs2
 
real(kind=kind_phys), save pidn0r
 
real(kind=kind_phys), save xlv1
 
real(kind=kind_phys), save pacrc
 
real(kind=kind_phys), save pi
 
real(kind=kind_phys), save bvtg1
 
real(kind=kind_phys), save bvtg2
 
real(kind=kind_phys), save bvtg3
 
real(kind=kind_phys), save bvtg4
 
real(kind=kind_phys), save g1pbg
 
real(kind=kind_phys), save g3pbg
 
real(kind=kind_phys), save g4pbg
 
real(kind=kind_phys), save g5pbgo2
 
real(kind=kind_phys), save pvtg
 
real(kind=kind_phys), save pacrg
 
real(kind=kind_phys), save precg1
 
real(kind=kind_phys), save precg2
 
real(kind=kind_phys), save pidn0g
 
real(kind=kind_phys), save rslopermax
 
real(kind=kind_phys), save rslopesmax
 
real(kind=kind_phys), save rslopegmax
 
real(kind=kind_phys), save rsloperbmax
 
real(kind=kind_phys), save rslopesbmax
 
real(kind=kind_phys), save rslopegbmax
 
real(kind=kind_phys), save rsloper2max
 
real(kind=kind_phys), save rslopes2max
 
real(kind=kind_phys), save rslopeg2max
 
real(kind=kind_phys), save rsloper3max
 
real(kind=kind_phys), save rslopes3max
 
real(kind=kind_phys), save rslopeg3max
 
real(kind=kind_phys), save, public pidn0s
 
real(kind=kind_phys), save, public pidnc
 

Function/Subroutine Documentation

◆ fpvs()

real(kind=kind_phys) function mp_wsm6::fpvs ( real(kind=kind_phys), intent(in)  t,
integer, intent(in)  ice,
real(kind=kind_phys), intent(in)  rd,
real(kind=kind_phys), intent(in)  rv,
real(kind=kind_phys), intent(in)  cvap,
real(kind=kind_phys), intent(in)  cliq,
real(kind=kind_phys), intent(in)  cice,
real(kind=kind_phys), intent(in)  hvap,
real(kind=kind_phys), intent(in)  hsub,
real(kind=kind_phys), intent(in)  psat,
real(kind=kind_phys), intent(in)  t0c 
)
private
1610 !=================================================================================================================
1611 
1612  integer,intent(in):: ice
1613  real(kind=kind_phys),intent(in):: cice,cliq,cvap,hsub,hvap,psat,rd,rv,t0c
1614  real(kind=kind_phys),intent(in):: t
1615 
1616  real(kind=kind_phys):: tr,ttp,dldt,dldti,xa,xb,xai,xbi
1617 
1618 !-----------------------------------------------------------------------------------------------------------------
1619 
1620  ttp=t0c+0.01
1621  dldt=cvap-cliq
1622  xa=-dldt/rv
1623  xb=xa+hvap/(rv*ttp)
1624  dldti=cvap-cice
1625  xai=-dldti/rv
1626  xbi=xai+hsub/(rv*ttp)
1627  tr=ttp/t
1628  if(t.lt.ttp.and.ice.eq.1) then
1629  fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
1630  else
1631  fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
1632  endif
1633 

◆ mp_wsm6_finalize()

subroutine, public mp_wsm6::mp_wsm6_finalize ( character(len=*), intent(out)  errmsg,
integer, intent(out)  errflg 
)

arg_table_mp_wsm6_finalize

\html

195 !=================================================================================================================
196 
197 !--- output arguments:
198  character(len=*),intent(out):: errmsg
199  integer,intent(out):: errflg
200 
201 !-----------------------------------------------------------------------------------------------------------------
202 
203  errmsg = 'mp_wsm6_finalize OK'
204  errflg = 0
205 

◆ mp_wsm6_init()

subroutine, public mp_wsm6::mp_wsm6_init ( real(kind=kind_phys), intent(in)  den0,
real(kind=kind_phys), intent(in)  denr,
real(kind=kind_phys), intent(in)  dens,
real(kind=kind_phys), intent(in)  cl,
real(kind=kind_phys), intent(in)  cpv,
integer, intent(in)  hail_opt,
character(len=*), intent(out)  errmsg,
integer, intent(out)  errflg 
)

arg_table_mp_wsm6_init

\html

73 !=================================================================================================================
74 
75 !input arguments:
76  integer,intent(in):: hail_opt ! RAS
77  real(kind=kind_phys),intent(in):: den0,denr,dens,cl,cpv
78 
79 !output arguments:
80  character(len=*),intent(out):: errmsg
81  integer,intent(out):: errflg
82 
83 !-----------------------------------------------------------------------------------------------------------------
84 
85 ! RAS13.1 define graupel parameters as graupel-like or hail-like,
86 ! depending on namelist option
87  if(hail_opt .eq. 1) then !Hail!
88  n0g = 4.e4
89  deng = 700.
90  avtg = 285.0
91  bvtg = 0.8
92  lamdagmax = 2.e4
93  else !Graupel!
94  n0g = 4.e6
95  deng = 500
96  avtg = 330.0
97  bvtg = 0.8
98  lamdagmax = 6.e4
99  endif
100 !
101  pi = 4.*atan(1.)
102  xlv1 = cl-cpv
103 !
104  qc0 = 4./3.*pi*denr*r0**3*xncr/den0 ! 0.419e-3 -- .61e-3
105  qck1 = .104*9.8*peaut/(xncr*denr)**(1./3.)/xmyu*den0**(4./3.) ! 7.03
106  pidnc = pi*denr/6. ! syb
107 !
108  bvtr1 = 1.+bvtr
109  bvtr2 = 2.5+.5*bvtr
110  bvtr3 = 3.+bvtr
111  bvtr4 = 4.+bvtr
112  bvtr6 = 6.+bvtr
113  g1pbr = rgmma(bvtr1)
114  g3pbr = rgmma(bvtr3)
115  g4pbr = rgmma(bvtr4) ! 17.837825
116  g6pbr = rgmma(bvtr6)
117  g5pbro2 = rgmma(bvtr2) ! 1.8273
118  pvtr = avtr*g4pbr/6.
119  eacrr = 1.0
120  pacrr = pi*n0r*avtr*g3pbr*.25*eacrr
121  precr1 = 2.*pi*n0r*.78
122  precr2 = 2.*pi*n0r*.31*avtr**.5*g5pbro2
123  roqimax = 2.08e22*dimax**8
124 !
125  bvts1 = 1.+bvts
126  bvts2 = 2.5+.5*bvts
127  bvts3 = 3.+bvts
128  bvts4 = 4.+bvts
129  g1pbs = rgmma(bvts1) !.8875
130  g3pbs = rgmma(bvts3)
131  g4pbs = rgmma(bvts4) ! 12.0786
132  g5pbso2 = rgmma(bvts2)
133  pvts = avts*g4pbs/6.
134  pacrs = pi*n0s*avts*g3pbs*.25
135  precs1 = 4.*n0s*.65
136  precs2 = 4.*n0s*.44*avts**.5*g5pbso2
137  pidn0r = pi*denr*n0r
138  pidn0s = pi*dens*n0s
139 !
140  pacrc = pi*n0s*avts*g3pbs*.25*eacrc
141 !
142  bvtg1 = 1.+bvtg
143  bvtg2 = 2.5+.5*bvtg
144  bvtg3 = 3.+bvtg
145  bvtg4 = 4.+bvtg
146  g1pbg = rgmma(bvtg1)
147  g3pbg = rgmma(bvtg3)
148  g4pbg = rgmma(bvtg4)
149  pacrg = pi*n0g*avtg*g3pbg*.25
150  g5pbgo2 = rgmma(bvtg2)
151  pvtg = avtg*g4pbg/6.
152  precg1 = 2.*pi*n0g*.78
153  precg2 = 2.*pi*n0g*.31*avtg**.5*g5pbgo2
154  pidn0g = pi*deng*n0g
155 !
156  rslopermax = 1./lamdarmax
157  rslopesmax = 1./lamdasmax
158  rslopegmax = 1./lamdagmax
159  rsloperbmax = rslopermax ** bvtr
160  rslopesbmax = rslopesmax ** bvts
161  rslopegbmax = rslopegmax ** bvtg
162  rsloper2max = rslopermax * rslopermax
163  rslopes2max = rslopesmax * rslopesmax
164  rslopeg2max = rslopegmax * rslopegmax
165  rsloper3max = rsloper2max * rslopermax
166  rslopes3max = rslopes2max * rslopesmax
167  rslopeg3max = rslopeg2max * rslopegmax
168 
169 !+---+-----------------------------------------------------------------+
170 !.. Set these variables needed for computing radar reflectivity. These
171 !.. get used within radar_init to create other variables used in the
172 !.. radar module.
173  xam_r = pi*denr/6.
174  xbm_r = 3.
175  xmu_r = 0.
176  xam_s = pi*dens/6.
177  xbm_s = 3.
178  xmu_s = 0.
179  xam_g = pi*deng/6.
180  xbm_g = 3.
181  xmu_g = 0.
182 
183  call radar_init
184 
185  errmsg = 'mp_wsm6_init OK'
186  errflg = 0
187 

Referenced by mp_wsm6_isohelper::mp_wsm6_init_c().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ mp_wsm6_run()

subroutine, public mp_wsm6::mp_wsm6_run ( real(kind=kind_phys), dimension(its:,:), intent(inout)  t,
real(kind=kind_phys), dimension(its:,:), intent(inout)  q,
real(kind=kind_phys), dimension(its:,:), intent(inout)  qc,
real(kind=kind_phys), dimension(its:,:), intent(inout)  qi,
real(kind=kind_phys), dimension(its:,:), intent(inout)  qr,
real(kind=kind_phys), dimension(its:,:), intent(inout)  qs,
real(kind=kind_phys), dimension(its:,:), intent(inout)  qg,
real(kind=kind_phys), dimension(its:,:), intent(in)  den,
real(kind=kind_phys), dimension(its:,:), intent(in)  p,
real(kind=kind_phys), dimension(its:,:), intent(in)  delz,
real(kind=kind_phys), intent(in)  delt,
real(kind=kind_phys), intent(in)  g,
real(kind=kind_phys), intent(in)  cpd,
real(kind=kind_phys), intent(in)  cpv,
real(kind=kind_phys), intent(in)  rd,
real(kind=kind_phys), intent(in)  rv,
real(kind=kind_phys), intent(in)  t0c,
real(kind=kind_phys), intent(in)  ep1,
real(kind=kind_phys), intent(in)  ep2,
real(kind=kind_phys), intent(in)  qmin,
real(kind=kind_phys), intent(in)  xls,
real(kind=kind_phys), intent(in)  xlv0,
real(kind=kind_phys), intent(in)  xlf0,
real(kind=kind_phys), intent(in)  den0,
real(kind=kind_phys), intent(in)  denr,
real(kind=kind_phys), intent(in)  cliq,
real(kind=kind_phys), intent(in)  cice,
real(kind=kind_phys), intent(in)  psat,
real(kind=kind_phys), dimension(its:), intent(inout)  rain,
real(kind=kind_phys), dimension(its:), intent(inout)  rainncv,
real(kind=kind_phys), dimension(its:), intent(inout)  sr,
real(kind=kind_phys), dimension(its:), intent(inout), optional  snow,
real(kind=kind_phys), dimension(its:), intent(inout), optional  snowncv,
real(kind=kind_phys), dimension(its:), intent(inout), optional  graupel,
real(kind=kind_phys), dimension(its:), intent(inout), optional  graupelncv,
real(kind=kind_phys), dimension(its:,:), intent(inout), optional  rainprod2d,
real(kind=kind_phys), dimension(its:,:), intent(inout), optional  evapprod2d,
integer, intent(in)  its,
integer, intent(in)  ite,
integer, intent(in)  kts,
integer, intent(in)  kte,
integer, intent(in), optional  microphysics_debug,
integer, intent(in), optional  diag_i_dbg,
integer, intent(in), optional  diag_j_dbg,
integer, intent(in), optional  diag_k_raw_base,
character(len=*), intent(out)  errmsg,
integer, intent(out)  errflg 
)

arg_table_mp_wsm6_run

\html

220 !=================================================================================================================!
221 ! This code is a 6-class GRAUPEL phase microphyiscs scheme (WSM6) of the
222 ! Single-Moment MicroPhyiscs (WSMMP). The WSMMP assumes that ice nuclei
223 ! number concentration is a function of temperature, and seperate assumption
224 ! is developed, in which ice crystal number concentration is a function
225 ! of ice amount. A theoretical background of the ice-microphysics and related
226 ! processes in the WSMMPs are described in Hong et al. (2004).
227 ! All production terms in the WSM6 scheme are described in Hong and Lim (2006).
228 ! All units are in m.k.s. and source/sink terms in kgkg-1s-1.
229 !
230 ! WSM6 cloud scheme
231 !
232 ! Coded by Song-You Hong and Jeong-Ock Jade Lim (Yonsei Univ.)
233 ! Summer 2003
234 !
235 ! Implemented by Song-You Hong (Yonsei Univ.) and Jimy Dudhia (NCAR)
236 ! Summer 2004
237 !
238 ! further modifications :
239 ! semi-lagrangian sedimentation (JH,2010),hong, aug 2009
240 ! ==> higher accuracy and efficient at lower resolutions
241 ! reflectivity computation from greg thompson, lim, jun 2011
242 ! ==> only diagnostic, but with removal of too large drops
243 ! add hail option from afwa, aug 2014
244 ! ==> switch graupel or hail by changing no, den, fall vel.
245 ! effective radius of hydrometeors, bae from kiaps, jan 2015
246 ! ==> consistency in solar insolation of rrtmg radiation
247 ! bug fix in melting terms, bae from kiaps, nov 2015
248 ! ==> density of air is divided, which has not been
249 !
250 ! Reference) Hong, Dudhia, Chen (HDC, 2004) Mon. Wea. Rev.
251 ! Hong and Lim (HL, 2006) J. Korean Meteor. Soc.
252 ! Dudhia, Hong and Lim (DHL, 2008) J. Meteor. Soc. Japan
253 ! Lin, Farley, Orville (LFO, 1983) J. Appl. Meteor.
254 ! Rutledge, Hobbs (RH83, 1983) J. Atmos. Sci.
255 ! Rutledge, Hobbs (RH84, 1984) J. Atmos. Sci.
256 ! Juang and Hong (JH, 2010) Mon. Wea. Rev.
257 !
258 
259 !input arguments:
260  integer,intent(in):: its,ite,kts,kte
261  integer,intent(in),optional :: microphysics_debug
262  integer,intent(in),optional :: diag_i_dbg
263  integer,intent(in),optional :: diag_j_dbg
264  integer,intent(in),optional :: diag_k_raw_base
265 
266  real(kind=kind_phys),intent(in),dimension(its:,:):: &
267  den, &
268  p, &
269  delz
270  real(kind=kind_phys),intent(in):: &
271  delt, &
272  g, &
273  cpd, &
274  cpv, &
275  t0c, &
276  den0, &
277  rd, &
278  rv, &
279  ep1, &
280  ep2, &
281  qmin, &
282  xls, &
283  xlv0, &
284  xlf0, &
285  cliq, &
286  cice, &
287  psat, &
288  denr
289 
290 !inout arguments:
291  real(kind=kind_phys),intent(inout),dimension(its:,:):: &
292  t
293  real(kind=kind_phys),intent(inout),dimension(its:,:):: &
294  q, &
295  qc, &
296  qi, &
297  qr, &
298  qs, &
299  qg
300  real(kind=kind_phys),intent(inout),dimension(its:):: &
301  rain, &
302  rainncv, &
303  sr
304 
305  real(kind=kind_phys),intent(inout),dimension(its:),optional:: &
306  snow, &
307  snowncv
308 
309  real(kind=kind_phys),intent(inout),dimension(its:),optional:: &
310  graupel, &
311  graupelncv
312 
313  real(kind=kind_phys),intent(inout),dimension(its:,:),optional:: &
314  rainprod2d, &
315  evapprod2d
316 
317 !output arguments:
318  character(len=*),intent(out):: errmsg
319  integer,intent(out):: errflg
320 
321 !local variables and arrays:
322  real(kind=kind_phys),dimension(its:ite,kts:kte,3):: &
323  rh, &
324  qsat, &
325  rslope, &
326  rslope2, &
327  rslope3, &
328  rslopeb, &
329  qrs_tmp, &
330  falk, &
331  fall, &
332  work1
333  real(kind=kind_phys),dimension(its:ite,kts:kte):: &
334  fallc, &
335  falkc, &
336  work1c, &
337  work2c, &
338  workr, &
339  worka
340  real(kind=kind_phys),dimension(its:ite,kts:kte):: &
341  den_tmp, &
342  delz_tmp
343  real(kind=kind_phys),dimension(its:ite,kts:kte):: &
344  pigen, &
345  pidep, &
346  pcond, &
347  prevp, &
348  psevp, &
349  pgevp, &
350  psdep, &
351  pgdep, &
352  praut, &
353  psaut, &
354  pgaut, &
355  piacr, &
356  pracw, &
357  praci, &
358  pracs, &
359  psacw, &
360  psaci, &
361  psacr, &
362  pgacw, &
363  pgaci, &
364  pgacr, &
365  pgacs, &
366  paacw, &
367  psmlt, &
368  pgmlt, &
369  pseml, &
370  pgeml, &
371  pimlt, &
372  pihmf, &
373  pihtf, &
374  pgfrz
375  real(kind=kind_phys),dimension(its:ite,kts:kte):: &
376  qsum, &
377  xl, &
378  cpm, &
379  work2, &
380  denfac, &
381  xni, &
382  denqrs1, &
383  denqrs2, &
384  denqrs3, &
385  denqci, &
386  n0sfac
387  real(kind=kind_phys),dimension(its:ite):: &
388  delqrs1, &
389  delqrs2, &
390  delqrs3, &
391  delqi
392  real(kind=kind_phys),dimension(its:ite):: &
393  tstepsnow, &
394  tstepgraup
395  integer,dimension(its:ite):: &
396  mstep, &
397  numdt
398  logical,dimension(its:ite):: flgcld
399  real(kind=kind_phys):: &
400  cpmcal, xlcal, diffus, &
401  viscos, xka, venfac, conden, diffac, &
402  x, y, z, a, b, c, d, e, &
403  qdt, holdrr, holdrs, holdrg, supcol, supcolt, pvt, &
404  coeres, supsat, dtcld, xmi, eacrs, satdt, &
405  qimax, diameter, xni0, roqi0, &
406  fallsum, fallsum_qsi, fallsum_qg, &
407  vt2i,vt2r,vt2s,vt2g,acrfac,egs,egi, &
408  xlwork2, factor, source, value, &
409  xlf, pfrzdtc, pfrzdtr, supice, alpha2, delta2, delta3
410 real(kind=kind_phys):: vt2ave
411 real(kind=kind_phys):: holdc, holdci
412 integer:: i, j, k, mstepmax, &
413  iprt, latd, lond, loop, loops, ifsat, n, idim, kdim, idbg_col
414  integer :: mpdbg_level
415  integer :: i_dbg_local
416  integer :: j_dbg_local
417  integer :: k_raw_base_local
418 
419 !Temporaries used for inlining fpvs function
420  real(kind=kind_phys):: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
421 
422 ! variables for optimization
423  real(kind=kind_phys),dimension(its:ite):: dvec1,tvec1
424  real(kind=kind_phys):: temp
425 
426 !-----------------------------------------------------------------------------------------------------------------
427 
428 ! compute internal functions
429 !
430  cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
431  xlcal(x) = xlv0-xlv1*(x-t0c)
432 !----------------------------------------------------------------
433 ! diffus: diffusion coefficient of the water vapor
434 ! viscos: kinematic viscosity(m2s-1)
435 ! Optimizatin : A**B => exp(log(A)*(B))
436 !
437  diffus(x,y) = 8.794e-5 * exp(log(x)*(1.81)) / y ! 8.794e-5*x**1.81/y
438  viscos(x,y) = 1.496e-6 * (x*sqrt(x)) /(x+120.)/y ! 1.496e-6*x**1.5/(x+120.)/y
439  xka(x,y) = 1.414e3*viscos(x,y)*y
440  diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b))
441  venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333))) &
442  /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
443  conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+d*d/(rv*e)*c/(a*a))
444 !
445 !
446  idim = ite-its+1
447  kdim = kte-kts+1
448  mpdbg_level = 0
449  if (present(microphysics_debug)) mpdbg_level = microphysics_debug
450  i_dbg_local = its
451  if (present(diag_i_dbg)) i_dbg_local = max(its, min(ite, diag_i_dbg))
452  j_dbg_local = -1
453  if (present(diag_j_dbg)) j_dbg_local = diag_j_dbg
454  k_raw_base_local = 1
455  if (present(diag_k_raw_base)) k_raw_base_local = diag_k_raw_base
456 !
457 !----------------------------------------------------------------
458 ! paddint 0 for negative values generated by dynamics
459 !
460  do k = kts, kte
461  do i = its, ite
462  qc(i,k) = max(qc(i,k),0.0)
463  qr(i,k) = max(qr(i,k),0.0)
464  qi(i,k) = max(qi(i,k),0.0)
465  qs(i,k) = max(qs(i,k),0.0)
466  qg(i,k) = max(qg(i,k),0.0)
467  enddo
468  enddo
469 !
470 !----------------------------------------------------------------
471 ! latent heat for phase changes and heat capacity. neglect the
472 ! changes during microphysical process calculation emanuel(1994)
473 !
474  do k = kts, kte
475  do i = its, ite
476  cpm(i,k) = cpmcal(q(i,k))
477  xl(i,k) = xlcal(t(i,k))
478  enddo
479  enddo
480  do k = kts, kte
481  do i = its, ite
482  delz_tmp(i,k) = delz(i,k)
483  den_tmp(i,k) = den(i,k)
484  enddo
485  enddo
486 !
487 !----------------------------------------------------------------
488 ! initialize the surface rain, snow, graupel
489 !
490  do i = its, ite
491  rainncv(i) = 0.
492  if(present(snowncv) .and. present(snow)) snowncv(i) = 0.
493  if(present(graupelncv) .and. present(graupel)) graupelncv(i) = 0.
494  sr(i) = 0.
495 ! new local array to catch step snow and graupel
496  tstepsnow(i) = 0.
497  tstepgraup(i) = 0.
498  enddo
499 !
500 !----------------------------------------------------------------
501 ! compute the minor time steps.
502 !
503  loops = max(nint(delt/dtcldcr),1)
504  dtcld = delt/loops
505  if(delt.le.dtcldcr) dtcld = delt
506 !
507  do loop = 1,loops
508 !
509 !----------------------------------------------------------------
510 ! initialize the large scale variables
511 !
512  do i = its, ite
513  mstep(i) = 1
514  flgcld(i) = .true.
515  enddo
516 !
517 ! do k = kts, kte
518 ! do i = its, ite
519 ! denfac(i,k) = sqrt(den0/den(i,k))
520 ! enddo
521 ! enddo
522  do k = kts, kte
523  do i = its,ite
524  dvec1(i) = den(i,k)
525  enddo
526  call vrec(tvec1,dvec1,ite-its+1)
527  do i = its, ite
528  tvec1(i) = tvec1(i)*den0
529  enddo
530  call vsqrt(dvec1,tvec1,ite-its+1)
531  do i = its,ite
532  denfac(i,k) = dvec1(i)
533  enddo
534  enddo
535 !
536 ! Inline expansion for fpvs
537 ! qsat(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
538 ! qsat(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
539  hsub = xls
540  hvap = xlv0
541  cvap = cpv
542  ttp=t0c+0.01
543  dldt=cvap-cliq
544  xa=-dldt/rv
545  xb=xa+hvap/(rv*ttp)
546  dldti=cvap-cice
547  xai=-dldti/rv
548  xbi=xai+hsub/(rv*ttp)
549  do k = kts, kte
550  do i = its, ite
551  tr=ttp/t(i,k)
552  qsat(i,k,1)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
553  qsat(i,k,1) = min(qsat(i,k,1),0.99*p(i,k))
554  qsat(i,k,1) = ep2 * qsat(i,k,1) / (p(i,k) - qsat(i,k,1))
555  qsat(i,k,1) = max(qsat(i,k,1),qmin)
556  rh(i,k,1) = max(q(i,k) / qsat(i,k,1),qmin)
557  tr=ttp/t(i,k)
558  if(t(i,k).lt.ttp) then
559  qsat(i,k,2)=psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
560  else
561  qsat(i,k,2)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
562  endif
563  qsat(i,k,2) = min(qsat(i,k,2),0.99*p(i,k))
564  qsat(i,k,2) = ep2 * qsat(i,k,2) / (p(i,k) - qsat(i,k,2))
565  qsat(i,k,2) = max(qsat(i,k,2),qmin)
566  rh(i,k,2) = max(q(i,k) / qsat(i,k,2),qmin)
567  enddo
568  enddo
569 !
570 !----------------------------------------------------------------
571 ! initialize the variables for microphysical physics
572 !
573 !
574  ! WSM6-F90 TAG: RATES_ZERO
575  ! legacy_group: G2
576  ! process: Zero all process rates each sub-step
577  ! compare_vars: prevp, psdep, pgdep, praut, psaut, pgaut, pracw, praci, psaci, pracs, pidep, pcond, psmlt, pgmlt, pseml, psevp, fall_r, fall_s, fall_g, fallc
578  do k = kts, kte
579  do i = its, ite
580  prevp(i,k) = 0.
581  psdep(i,k) = 0.
582  pgdep(i,k) = 0.
583  praut(i,k) = 0.
584  psaut(i,k) = 0.
585  pgaut(i,k) = 0.
586  pracw(i,k) = 0.
587  praci(i,k) = 0.
588  piacr(i,k) = 0.
589  psaci(i,k) = 0.
590  psacw(i,k) = 0.
591  pracs(i,k) = 0.
592  psacr(i,k) = 0.
593  pgacw(i,k) = 0.
594  paacw(i,k) = 0.
595  pgaci(i,k) = 0.
596  pgacr(i,k) = 0.
597  pgacs(i,k) = 0.
598  pigen(i,k) = 0.
599  pidep(i,k) = 0.
600  pcond(i,k) = 0.
601  psmlt(i,k) = 0.
602  pgmlt(i,k) = 0.
603  pseml(i,k) = 0.
604  pgeml(i,k) = 0.
605  psevp(i,k) = 0.
606  pgevp(i,k) = 0.
607  falk(i,k,1) = 0.
608  falk(i,k,2) = 0.
609  falk(i,k,3) = 0.
610  fall(i,k,1) = 0.
611  fall(i,k,2) = 0.
612  fall(i,k,3) = 0.
613  fallc(i,k) = 0.
614  falkc(i,k) = 0.
615  xni(i,k) = 1.e3
616  enddo
617  enddo
618 ! WSM6-F90 TAG: XNI
619 ! legacy_group: G3
620 ! process: Ice crystal number concentration
621 ! compare_vars: xni, qi, den
622 !-------------------------------------------------------------
623 ! Ni: ice crystal number concentraiton [HDC 5c]
624 !-------------------------------------------------------------
625  do k = kts, kte
626  do i = its, ite
627  temp = (den(i,k)*max(qi(i,k),qmin))
628  temp = sqrt(sqrt(temp*temp*temp))
629  xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
630  enddo
631  enddo
632 !
633 !----------------------------------------------------------------
634 ! compute the fallout term:
635 ! first, vertical terminal velosity for minor loops
636 !----------------------------------------------------------------
637  ! WSM6-F90 TAG: SLOPE1
638  ! legacy_group: G4
639  ! process: First slope calculation
640  ! compare_vars: rslope, rslope2, rslope3, rslopeb, falk, fall, work1
641  do k = kts, kte
642  do i = its, ite
643  qrs_tmp(i,k,1) = qr(i,k)
644  qrs_tmp(i,k,2) = qs(i,k)
645  qrs_tmp(i,k,3) = qg(i,k)
646  enddo
647  enddo
648  call slope_wsm6(qrs_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2,rslope3, &
649  work1,its,ite,kts,kte)
650 !
651  do k = kte, kts, -1
652  do i = its, ite
653  workr(i,k) = work1(i,k,1)
654  qsum(i,k) = max( (qs(i,k)+qg(i,k)), 1.e-15)
655  if( qsum(i,k) .gt. 1.e-15 ) then
656  worka(i,k) = (work1(i,k,2)*qs(i,k) + work1(i,k,3)*qg(i,k)) &
657  / qsum(i,k)
658  else
659  worka(i,k) = 0.
660  endif
661  denqrs1(i,k) = den(i,k)*qr(i,k)
662  denqrs2(i,k) = den(i,k)*qs(i,k)
663  denqrs3(i,k) = den(i,k)*qg(i,k)
664  if(qr(i,k).le.0.0) workr(i,k) = 0.0
665  enddo
666  enddo
667  call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,workr,denqrs1, &
668  delqrs1,dtcld,1,1)
669  call nislfv_rain_plm6(idim,kdim,den_tmp,denfac,t,delz_tmp,worka, &
670  denqrs2,denqrs3,delqrs2,delqrs3,dtcld,1,1)
671  do k = kts, kte
672  do i = its, ite
673  qr(i,k) = max(denqrs1(i,k)/den(i,k),0.)
674  qs(i,k) = max(denqrs2(i,k)/den(i,k),0.)
675  qg(i,k) = max(denqrs3(i,k)/den(i,k),0.)
676  fall(i,k,1) = denqrs1(i,k)*workr(i,k)/delz(i,k)
677  fall(i,k,2) = denqrs2(i,k)*worka(i,k)/delz(i,k)
678  fall(i,k,3) = denqrs3(i,k)*worka(i,k)/delz(i,k)
679  enddo
680  enddo
681  do i = its, ite
682  fall(i,1,1) = delqrs1(i)/delz(i,1)/dtcld
683  fall(i,1,2) = delqrs2(i)/delz(i,1)/dtcld
684  fall(i,1,3) = delqrs3(i)/delz(i,1)/dtcld
685  enddo
686  ! WSM6-F90 TAG: SLOPE2
687  ! legacy_group: G6
688  ! process: Second slope calculation after sedimentation
689  ! compare_vars: rslope, rslope2, rslope3, rslopeb, falk, fall, work1
690  do k = kts, kte
691  do i = its, ite
692  qrs_tmp(i,k,1) = qr(i,k)
693  qrs_tmp(i,k,2) = qs(i,k)
694  qrs_tmp(i,k,3) = qg(i,k)
695  enddo
696  enddo
697  call slope_wsm6(qrs_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2,rslope3, &
698  work1,its,ite,kts,kte)
699  ! WSM6-F90 TAG: MELT
700  ! legacy_group: G7
701  ! process: Melting of snow/graupel and latent heating coupling
702  ! compare_vars: psmlt, pgmlt, t, qrs, qci, q
703 !
704  do k = kte, kts, -1
705  do i = its, ite
706  supcol = t0c-t(i,k)
707  n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
708  if(t(i,k).gt.t0c) then
709 !---------------------------------------------------------------
710 ! psmlt: melting of snow [HL A33] [RH83 A25]
711 ! (T>T0: S->R)
712 !---------------------------------------------------------------
713  xlf = xlf0
714  work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
715  if(qs(i,k).gt.0.) then
716  coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
717  psmlt(i,k) = xka(t(i,k),den(i,k))/xlf*(t0c-t(i,k))*pi/2. &
718  *n0sfac(i,k)*(precs1*rslope2(i,k,2) &
719  +precs2*work2(i,k)*coeres)/den(i,k)
720  psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep(i), &
721  -qs(i,k)/mstep(i)),0.)
722  qs(i,k) = qs(i,k) + psmlt(i,k)
723  qr(i,k) = qr(i,k) - psmlt(i,k)
724  t(i,k) = t(i,k) + xlf/cpm(i,k)*psmlt(i,k)
725  endif
726 !---------------------------------------------------------------
727 ! pgmlt: melting of graupel [HL A23] [LFO 47]
728 ! (T>T0: G->R)
729 !---------------------------------------------------------------
730  if(qg(i,k).gt.0.) then
731  coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3))
732  pgmlt(i,k) = xka(t(i,k),den(i,k))/xlf &
733  *(t0c-t(i,k))*(precg1*rslope2(i,k,3) &
734  +precg2*work2(i,k)*coeres)/den(i,k)
735  pgmlt(i,k) = min(max(pgmlt(i,k)*dtcld/mstep(i), &
736  -qg(i,k)/mstep(i)),0.)
737  qg(i,k) = qg(i,k) + pgmlt(i,k)
738  qr(i,k) = qr(i,k) - pgmlt(i,k)
739  t(i,k) = t(i,k) + xlf/cpm(i,k)*pgmlt(i,k)
740  endif
741  endif
742  enddo
743  enddo
744 ! WSM6-F90 TAG: VICE
745 ! legacy_group: G8
746 ! process: Cloud ice sedimentation/fallout
747 ! compare_vars: qci, fall, fallc, den, p
748 !---------------------------------------------------------------
749 ! Vice [ms-1] : fallout of ice crystal [HDC 5a]
750 !---------------------------------------------------------------
751  do k = kte, kts, -1
752  do i = its, ite
753  if(qi(i,k).le.0.) then
754  work1c(i,k) = 0.
755  else
756  xmi = den(i,k)*qi(i,k)/xni(i,k)
757  diameter = max(min(dicon * sqrt(xmi),dimax), 1.e-25)
758  work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31))
759  endif
760  enddo
761  enddo
762 !
763 ! forward semi-laglangian scheme (JH), PCM (piecewise constant), (linear)
764 !
765  do k = kte, kts, -1
766  do i = its, ite
767  denqci(i,k) = den(i,k)*qi(i,k)
768  enddo
769  enddo
770  call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,work1c,denqci, &
771  delqi,dtcld,1,0)
772  do k = kts, kte
773  do i = its, ite
774  qi(i,k) = max(denqci(i,k)/den(i,k),0.)
775  enddo
776  enddo
777  do i = its, ite
778  fallc(i,1) = delqi(i)/delz(i,1)/dtcld
779  enddo
780 !
781 !----------------------------------------------------------------
782 ! rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
783 !
784  ! WSM6-F90 TAG: PRECIP
785  ! legacy_group: G9
786  ! process: Surface precipitation accumulation
787  ! compare_vars: rainncv, snowncv, graupelncv, sr
788  do i = its, ite
789  fallsum = fall(i,kts,1)+fall(i,kts,2)+fall(i,kts,3)+fallc(i,kts)
790  fallsum_qsi = fall(i,kts,2)+fallc(i,kts)
791  fallsum_qg = fall(i,kts,3)
792  if(fallsum.gt.0.) then
793  rainncv(i) = fallsum*delz(i,kts)/denr*dtcld*1000. + rainncv(i)
794  rain(i) = fallsum*delz(i,kts)/denr*dtcld*1000. + rain(i)
795  endif
796  if(fallsum_qsi.gt.0.) then
797  tstepsnow(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. &
798  + tstepsnow(i)
799  if(present(snowncv) .and. present(snow)) then
800  snowncv(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. &
801  + snowncv(i)
802  snow(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + snow(i)
803  endif
804  endif
805  if(fallsum_qg.gt.0.) then
806  tstepgraup(i) = fallsum_qg*delz(i,kts)/denr*dtcld*1000. &
807  + tstepgraup(i)
808  if(present (graupelncv) .and. present (graupel)) then
809  graupelncv(i) = fallsum_qg*delz(i,kts)/denr*dtcld*1000. &
810  + graupelncv(i)
811  graupel(i) = fallsum_qg*delz(i,kts)/denr*dtcld*1000. + graupel(i)
812  endif
813  endif
814  if(present (snowncv)) then
815  if(fallsum.gt.0.)sr(i)=(snowncv(i) + graupelncv(i))/(rainncv(i)+1.e-12)
816  else
817  if(fallsum.gt.0.)sr(i)=(tstepsnow(i) + tstepgraup(i))/(rainncv(i)+1.e-12)
818  endif
819  enddo
820 !
821 ! WSM6-F90 TAG: PHASE
822 ! legacy_group: G10
823 ! process: Instantaneous phase changes
824 ! compare_vars: t, q, qci, qrs
825 !---------------------------------------------------------------
826 ! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
827 ! (T>T0: I->C)
828 !---------------------------------------------------------------
829  do k = kts, kte
830  do i = its, ite
831  pimlt(i,k) = 0.
832  pihmf(i,k) = 0.
833  pihtf(i,k) = 0.
834  pgfrz(i,k) = 0.
835  supcol = t0c-t(i,k)
836  xlf = xls-xl(i,k)
837  if(supcol.lt.0.) xlf = xlf0
838  if(supcol.lt.0.and.qi(i,k).gt.0.) then
839  pimlt(i,k) = qi(i,k)
840  qc(i,k) = qc(i,k) + qi(i,k)
841  t(i,k) = t(i,k) - xlf/cpm(i,k)*qi(i,k)
842  qi(i,k) = 0.
843  endif
844 !---------------------------------------------------------------
845 ! pihmf: homogeneous freezing of cloud water below -40c [HL A45]
846 ! (T<-40C: C->I)
847 !---------------------------------------------------------------
848  if(supcol.gt.40..and.qc(i,k).gt.0.) then
849  pihmf(i,k) = qc(i,k)
850  qi(i,k) = qi(i,k) + qc(i,k)
851  t(i,k) = t(i,k) + xlf/cpm(i,k)*qc(i,k)
852  qc(i,k) = 0.
853  endif
854 !---------------------------------------------------------------
855 ! pihtf: heterogeneous freezing of cloud water [HL A44]
856 ! (T0>T>-40C: C->I)
857 !---------------------------------------------------------------
858  if(supcol.gt.0..and.qc(i,k).gt.qmin) then
859 ! pfrzdtc = min(pfrz1*(exp(pfrz2*supcol)-1.) &
860 ! * den(i,k)/denr/xncr*qc(i,k)**2*dtcld,qc(i,k))
861  supcolt=min(supcol,50.)
862  pfrzdtc = min(pfrz1*(exp(pfrz2*supcolt)-1.) &
863  * den(i,k)/denr/xncr*qc(i,k)*qc(i,k)*dtcld,qc(i,k))
864  pihtf(i,k) = pfrzdtc
865  qi(i,k) = qi(i,k) + pfrzdtc
866  t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtc
867  qc(i,k) = qc(i,k)-pfrzdtc
868  endif
869 !---------------------------------------------------------------
870 ! pgfrz: freezing of rain water [HL A20] [LFO 45]
871 ! (T<T0, R->G)
872 !---------------------------------------------------------------
873  if(supcol.gt.0..and.qr(i,k).gt.0.) then
874 ! pfrzdtr = min(20.*pi**2*pfrz1*n0r*denr/den(i,k) &
875 ! * (exp(pfrz2*supcol)-1.)*rslope3(i,k,1)**2 &
876 ! * rslope(i,k,1)*dtcld,qr(i,k))
877  temp = rslope3(i,k,1)
878  temp = temp*temp*rslope(i,k,1)
879  supcolt=min(supcol,50.)
880  pfrzdtr = min(20.*(pi*pi)*pfrz1*n0r*denr/den(i,k) &
881  *(exp(pfrz2*supcolt)-1.)*temp*dtcld, &
882  qr(i,k))
883  pgfrz(i,k) = pfrzdtr
884  qg(i,k) = qg(i,k) + pfrzdtr
885  t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtr
886  qr(i,k) = qr(i,k)-pfrzdtr
887  endif
888  enddo
889  enddo
890 !
891 !
892 !----------------------------------------------------------------
893 ! update the slope parameters for microphysics computation
894 !
895  ! WSM6-F90 TAG: SLOPE3
896  ! legacy_group: G11
897  ! process: Third slope calculation
898  ! compare_vars: rslope, rslope2, rslope3, rslopeb, falk, fall, work1
899  do k = kts, kte
900  do i = its, ite
901  qrs_tmp(i,k,1) = qr(i,k)
902  qrs_tmp(i,k,2) = qs(i,k)
903  qrs_tmp(i,k,3) = qg(i,k)
904  enddo
905  enddo
906  call slope_wsm6(qrs_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2,rslope3, &
907  work1,its,ite,kts,kte)
908 ! WSM6-F90 TAG: DIFF_PREP
909 ! legacy_group: G12
910 ! process: Prepare diffusion/work terms
911 ! compare_vars: workdiffw, workdiffi, work2
912 !------------------------------------------------------------------
913 ! work1: the thermodynamic term in the denominator associated with
914 ! heat conduction and vapor diffusion
915 ! (ry88, y93, h85)
916 ! work2: parameter associated with the ventilation effects(y93)
917 !
918  do k = kts, kte
919  do i = its, ite
920  work1(i,k,1) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qsat(i,k,1))
921  work1(i,k,2) = diffac(xls,p(i,k),t(i,k),den(i,k),qsat(i,k,2))
922  work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
923  enddo
924  enddo
925 !
926 !===============================================================
927 !
928 ! warm rain processes
929 !
930 ! - follows the processes in RH83 and LFO except for autoconcersion
931 !
932 !===============================================================
933 !
934  do k = kts, kte
935  do i = its, ite
936  supsat = max(q(i,k),qmin)-qsat(i,k,1)
937  satdt = supsat/dtcld
938 !---------------------------------------------------------------
939 ! praut: auto conversion rate from cloud to rain [HDC 16]
940 ! (C->R)
941 !---------------------------------------------------------------
942  if(qc(i,k).gt.qc0) then
943  praut(i,k) = qck1*qc(i,k)**(7./3.)
944  praut(i,k) = min(praut(i,k),qc(i,k)/dtcld)
945  endif
946 !---------------------------------------------------------------
947 ! pracw: accretion of cloud water by rain [HL A40] [LFO 51]
948 ! (C->R)
949 !---------------------------------------------------------------
950  if(qr(i,k).gt.qcrmin.and.qc(i,k).gt.qmin) then
951  pracw(i,k) = min(pacrr*rslope3(i,k,1)*rslopeb(i,k,1) &
952  * qc(i,k)*denfac(i,k),qc(i,k)/dtcld)
953  endif
954 !---------------------------------------------------------------
955 ! prevp: evaporation/condensation rate of rain [HDC 14]
956 ! (V->R or R->V)
957 !---------------------------------------------------------------
958  if(qr(i,k).gt.0.) then
959  coeres = rslope2(i,k,1)*sqrt(rslope(i,k,1)*rslopeb(i,k,1))
960  prevp(i,k) = (rh(i,k,1)-1.)*(precr1*rslope2(i,k,1) &
961  + precr2*work2(i,k)*coeres)/work1(i,k,1)
962  if(prevp(i,k).lt.0.) then
963  prevp(i,k) = max(prevp(i,k),-qr(i,k)/dtcld)
964  prevp(i,k) = max(prevp(i,k),satdt/2)
965  else
966  prevp(i,k) = min(prevp(i,k),satdt/2)
967  endif
968  endif
969  enddo
970  enddo
971 !
972 !===============================================================
973 !
974 ! cold rain processes
975 !
976 ! - follows the revised ice microphysics processes in HDC
977 ! - the processes same as in RH83 and RH84 and LFO behave
978 ! following ice crystal hapits defined in HDC, inclduing
979 ! intercept parameter for snow (n0s), ice crystal number
980 ! concentration (ni), ice nuclei number concentration
981 ! (n0i), ice diameter (d)
982 !
983 !===============================================================
984 !
985  do k = kts, kte
986  do i = its, ite
987  supcol = t0c-t(i,k)
988  n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
989  supsat = max(q(i,k),qmin)-qsat(i,k,2)
990  satdt = supsat/dtcld
991  ifsat = 0
992 !-------------------------------------------------------------
993 ! Ni: ice crystal number concentraiton [HDC 5c]
994 !-------------------------------------------------------------
995 ! xni(i,k) = min(max(5.38e7*(den(i,k) &
996 ! * max(qi(i,k),qmin))**0.75,1.e3),1.e6)
997  temp = (den(i,k)*max(qi(i,k),qmin))
998  temp = sqrt(sqrt(temp*temp*temp))
999  xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
1000  eacrs = exp(0.07*(-supcol))
1001 !
1002  xmi = den(i,k)*qi(i,k)/xni(i,k)
1003  diameter = min(dicon * sqrt(xmi),dimax)
1004  vt2i = 1.49e4*diameter**1.31
1005  vt2r=pvtr*rslopeb(i,k,1)*denfac(i,k)
1006  vt2s=pvts*rslopeb(i,k,2)*denfac(i,k)
1007  vt2g=pvtg*rslopeb(i,k,3)*denfac(i,k)
1008  qsum(i,k) = max( (qs(i,k)+qg(i,k)), 1.e-15)
1009  if(qsum(i,k) .gt. 1.e-15) then
1010  vt2ave=(vt2s*qs(i,k)+vt2g*qg(i,k))/(qsum(i,k))
1011  else
1012  vt2ave=0.
1013  endif
1014  if(supcol.gt.0.and.qi(i,k).gt.qmin) then
1015  if(qr(i,k).gt.qcrmin) then
1016 ! WSM6-F90 TAG: PRACI
1017 ! legacy_group: G13b
1018 ! process: Accretion of cloud ice by rain
1019 ! compare_vars: praci, qi, qr, den
1020 !-------------------------------------------------------------
1021 ! praci: accretion of cloud ice by rain [HL A15] [LFO 25]
1022 ! (T<T0: I->R)
1023 !-------------------------------------------------------------
1024  acrfac = 2.*rslope3(i,k,1)+2.*diameter*rslope2(i,k,1) &
1025  + diameter**2*rslope(i,k,1)
1026  praci(i,k) = pi*qi(i,k)*n0r*abs(vt2r-vt2i)*acrfac/4.
1027 ! reduce collection efficiency (suggested by B. Wilt)
1028  praci(i,k) = praci(i,k)*min(max(0.0,qr(i,k)/qi(i,k)),1.)**2
1029  praci(i,k) = min(praci(i,k),qi(i,k)/dtcld)
1030 !-------------------------------------------------------------
1031 ! piacr: accretion of rain by cloud ice [HL A19] [LFO 26]
1032 ! (T<T0: R->S or R->G)
1033 !-------------------------------------------------------------
1034  piacr(i,k) = pi**2*avtr*n0r*denr*xni(i,k)*denfac(i,k) &
1035  * g6pbr*rslope3(i,k,1)*rslope3(i,k,1) &
1036  * rslopeb(i,k,1)/24./den(i,k)
1037 ! reduce collection efficiency (suggested by B. Wilt)
1038  piacr(i,k) = piacr(i,k)*min(max(0.0,qi(i,k)/qr(i,k)),1.)**2
1039  piacr(i,k) = min(piacr(i,k),qr(i,k)/dtcld)
1040  endif
1041 ! WSM6-F90 TAG: PSACI
1042 ! legacy_group: G13e
1043 ! process: Accretion of cloud ice by snow
1044 ! compare_vars: psaci, qi, qs, den
1045 !-------------------------------------------------------------
1046 ! psaci: accretion of cloud ice by snow [HDC 10]
1047 ! (T<T0: I->S)
1048 !-------------------------------------------------------------
1049  if(qs(i,k).gt.qcrmin) then
1050  acrfac = 2.*rslope3(i,k,2)+2.*diameter*rslope2(i,k,2) &
1051  + diameter**2*rslope(i,k,2)
1052  psaci(i,k) = pi*qi(i,k)*eacrs*n0s*n0sfac(i,k) &
1053  * abs(vt2ave-vt2i)*acrfac/4.
1054  psaci(i,k) = min(psaci(i,k),qi(i,k)/dtcld)
1055  endif
1056 !-------------------------------------------------------------
1057 ! pgaci: accretion of cloud ice by graupel [HL A17] [LFO 41]
1058 ! (T<T0: I->G)
1059 !-------------------------------------------------------------
1060  if(qg(i,k).gt.qcrmin) then
1061  egi = exp(0.07*(-supcol))
1062  acrfac = 2.*rslope3(i,k,3)+2.*diameter*rslope2(i,k,3) &
1063  + diameter**2*rslope(i,k,3)
1064  pgaci(i,k) = pi*egi*qi(i,k)*n0g*abs(vt2ave-vt2i)*acrfac/4.
1065  pgaci(i,k) = min(pgaci(i,k),qi(i,k)/dtcld)
1066  endif
1067  endif
1068 !-------------------------------------------------------------
1069 ! psacw: accretion of cloud water by snow [HL A7] [LFO 24]
1070 ! (T<T0: C->S, and T>=T0: C->R)
1071 !-------------------------------------------------------------
1072  if(qs(i,k).gt.qcrmin.and.qc(i,k).gt.qmin) then
1073  psacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2)*rslopeb(i,k,2) &
1074 ! reduce collection efficiency (suggested by B. Wilt)
1075  * min(max(0.0,qs(i,k)/qc(i,k)),1.)**2 &
1076  * qc(i,k)*denfac(i,k),qc(i,k)/dtcld)
1077  endif
1078 !-------------------------------------------------------------
1079 ! pgacw: accretion of cloud water by graupel [HL A6] [LFO 40]
1080 ! (T<T0: C->G, and T>=T0: C->R)
1081 !-------------------------------------------------------------
1082  if(qg(i,k).gt.qcrmin.and.qc(i,k).gt.qmin) then
1083  pgacw(i,k) = min(pacrg*rslope3(i,k,3)*rslopeb(i,k,3) &
1084 ! reduce collection efficiency (suggested by B. Wilt)
1085  * min(max(0.0,qg(i,k)/qc(i,k)),1.)**2 &
1086  * qc(i,k)*denfac(i,k),qc(i,k)/dtcld)
1087  endif
1088 !-------------------------------------------------------------
1089 ! paacw: accretion of cloud water by averaged snow/graupel
1090 ! (T<T0: C->G or S, and T>=T0: C->R)
1091 !-------------------------------------------------------------
1092  if(qsum(i,k) .gt. 1.e-15) then
1093  paacw(i,k) = (qs(i,k)*psacw(i,k)+qg(i,k)*pgacw(i,k)) &
1094  /(qsum(i,k))
1095  endif
1096 ! WSM6-F90 TAG: PRACS
1097 ! legacy_group: G13f
1098 ! process: Accretion of snow by rain / rain-snow interaction
1099 ! compare_vars: pracs, qr, qs, den
1100 !-------------------------------------------------------------
1101 ! pracs: accretion of snow by rain [HL A11] [LFO 27]
1102 ! (T<T0: S->G)
1103 !-------------------------------------------------------------
1104  if(qs(i,k).gt.qcrmin.and.qr(i,k).gt.qcrmin) then
1105  if(supcol.gt.0) then
1106  acrfac = 5.*rslope3(i,k,2)*rslope3(i,k,2)*rslope(i,k,1) &
1107  + 2.*rslope3(i,k,2)*rslope2(i,k,2)*rslope2(i,k,1) &
1108  + .5*rslope2(i,k,2)*rslope2(i,k,2)*rslope3(i,k,1)
1109  pracs(i,k) = pi**2*n0r*n0s*n0sfac(i,k)*abs(vt2r-vt2ave) &
1110  * (dens/den(i,k))*acrfac
1111 ! reduce collection efficiency (suggested by B. Wilt)
1112  pracs(i,k) = pracs(i,k)*min(max(0.0,qr(i,k)/qs(i,k)),1.)**2
1113  pracs(i,k) = min(pracs(i,k),qs(i,k)/dtcld)
1114  endif
1115 !-------------------------------------------------------------
1116 ! psacr: accretion of rain by snow [HL A10] [LFO 28]
1117 ! (T<T0:R->S or R->G) (T>=T0: enhance melting of snow)
1118 !-------------------------------------------------------------
1119  acrfac = 5.*rslope3(i,k,1)*rslope3(i,k,1)*rslope(i,k,2) &
1120  + 2.*rslope3(i,k,1)*rslope2(i,k,1)*rslope2(i,k,2) &
1121  +.5*rslope2(i,k,1)*rslope2(i,k,1)*rslope3(i,k,2)
1122  psacr(i,k) = pi**2*n0r*n0s*n0sfac(i,k)*abs(vt2ave-vt2r) &
1123  * (denr/den(i,k))*acrfac
1124 ! reduce collection efficiency (suggested by B. Wilt)
1125  psacr(i,k) = psacr(i,k)*min(max(0.0,qs(i,k)/qr(i,k)),1.)**2
1126  psacr(i,k) = min(psacr(i,k),qr(i,k)/dtcld)
1127  endif
1128 !-------------------------------------------------------------
1129 ! pgacr: accretion of rain by graupel [HL A12] [LFO 42]
1130 ! (T<T0: R->G) (T>=T0: enhance melting of graupel)
1131 !-------------------------------------------------------------
1132  if(qg(i,k).gt.qcrmin.and.qr(i,k).gt.qcrmin) then
1133  acrfac = 5.*rslope3(i,k,1)*rslope3(i,k,1)*rslope(i,k,3) &
1134  + 2.*rslope3(i,k,1)*rslope2(i,k,1)*rslope2(i,k,3) &
1135  + .5*rslope2(i,k,1)*rslope2(i,k,1)*rslope3(i,k,3)
1136  pgacr(i,k) = pi**2*n0r*n0g*abs(vt2ave-vt2r)*(denr/den(i,k)) &
1137  * acrfac
1138 ! reduce collection efficiency (suggested by B. Wilt)
1139  pgacr(i,k) = pgacr(i,k)*min(max(0.0,qg(i,k)/qr(i,k)),1.)**2
1140  pgacr(i,k) = min(pgacr(i,k),qr(i,k)/dtcld)
1141  endif
1142 !
1143 !-------------------------------------------------------------
1144 ! pgacs: accretion of snow by graupel [HL A13] [LFO 29]
1145 ! (S->G): This process is eliminated in V3.0 with the
1146 ! new combined snow/graupel fall speeds
1147 !-------------------------------------------------------------
1148  if(qg(i,k).gt.qcrmin.and.qs(i,k).gt.qcrmin) then
1149  pgacs(i,k) = 0.
1150  endif
1151  if(supcol.le.0) then
1152  xlf = xlf0
1153 ! WSM6-F90 TAG: PSEML
1154 ! legacy_group: G13g
1155 ! process: Snow evaporation/sublimation
1156 ! compare_vars: pseml, qs, qv, qsat, den
1157 !-------------------------------------------------------------
1158 ! pseml: enhanced melting of snow by accretion of water [HL A34]
1159 ! (T>=T0: S->R)
1160 !-------------------------------------------------------------
1161  if(qs(i,k).gt.0.) &
1162  pseml(i,k) = min(max(cliq*supcol*(paacw(i,k)+psacr(i,k)) &
1163  / xlf,-qs(i,k)/dtcld),0.)
1164 !-------------------------------------------------------------
1165 ! pgeml: enhanced melting of graupel by accretion of water [HL A24] [RH84 A21-A22]
1166 ! (T>=T0: G->R)
1167 !-------------------------------------------------------------
1168  if(qg(i,k).gt.0.) &
1169  pgeml(i,k) = min(max(cliq*supcol*(paacw(i,k)+pgacr(i,k)) &
1170  / xlf,-qg(i,k)/dtcld),0.)
1171  endif
1172  if(supcol.gt.0) then
1173 ! WSM6-F90 TAG: PIDEP
1174 ! legacy_group: G13h
1175 ! process: Ice deposition/sublimation
1176 ! compare_vars: pidep, qi, qv, qsat, den, t
1177 !-------------------------------------------------------------
1178 ! pidep: deposition/Sublimation rate of ice [HDC 9]
1179 ! (T<T0: V->I or I->V)
1180 !-------------------------------------------------------------
1181  if(qi(i,k).gt.0.and.ifsat.ne.1) then
1182  pidep(i,k) = 4.*diameter*xni(i,k)*(rh(i,k,2)-1.)/work1(i,k,2)
1183  supice = satdt-prevp(i,k)
1184  if(pidep(i,k).lt.0.) then
1185  pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
1186  pidep(i,k) = max(pidep(i,k),-qi(i,k)/dtcld)
1187  else
1188  pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
1189  endif
1190  if(abs(prevp(i,k)+pidep(i,k)).ge.abs(satdt)) ifsat = 1
1191  endif
1192 !-------------------------------------------------------------
1193 ! psdep: deposition/sublimation rate of snow [HDC 14]
1194 ! (T<T0: V->S or S->V)
1195 !-------------------------------------------------------------
1196  if(qs(i,k).gt.0..and.ifsat.ne.1) then
1197  coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
1198  psdep(i,k) = (rh(i,k,2)-1.)*n0sfac(i,k)*(precs1*rslope2(i,k,2) &
1199  + precs2*work2(i,k)*coeres)/work1(i,k,2)
1200  supice = satdt-prevp(i,k)-pidep(i,k)
1201  if(psdep(i,k).lt.0.) then
1202  psdep(i,k) = max(psdep(i,k),-qs(i,k)/dtcld)
1203  psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
1204  else
1205  psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
1206  endif
1207  if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt)) &
1208  ifsat = 1
1209  endif
1210 !-------------------------------------------------------------
1211 ! pgdep: deposition/sublimation rate of graupel [HL A21] [LFO 46]
1212 ! (T<T0: V->G or G->V)
1213 !-------------------------------------------------------------
1214  if(qg(i,k).gt.0..and.ifsat.ne.1) then
1215  coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3))
1216  pgdep(i,k) = (rh(i,k,2)-1.)*(precg1*rslope2(i,k,3) &
1217  + precg2*work2(i,k)*coeres)/work1(i,k,2)
1218  supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
1219  if(pgdep(i,k).lt.0.) then
1220  pgdep(i,k) = max(pgdep(i,k),-qg(i,k)/dtcld)
1221  pgdep(i,k) = max(max(pgdep(i,k),satdt/2),supice)
1222  else
1223  pgdep(i,k) = min(min(pgdep(i,k),satdt/2),supice)
1224  endif
1225  if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)+pgdep(i,k)).ge. &
1226  abs(satdt)) ifsat = 1
1227  endif
1228 !-------------------------------------------------------------
1229 ! pigen: generation(nucleation) of ice from vapor [HL 50] [HDC 7-8]
1230 ! (T<T0: V->I)
1231 !-------------------------------------------------------------
1232  if(supsat.gt.0.and.ifsat.ne.1) then
1233  supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)-pgdep(i,k)
1234  xni0 = 1.e3*exp(0.1*supcol)
1235  roqi0 = 4.92e-11*xni0**1.33
1236  pigen(i,k) = max(0.,(roqi0/den(i,k)-max(qi(i,k),0.))/dtcld)
1237  pigen(i,k) = min(min(pigen(i,k),satdt),supice)
1238  endif
1239 !
1240 ! WSM6-F90 TAG: PSAUT
1241 ! legacy_group: G13i
1242 ! process: Autoconversion to snow
1243 ! compare_vars: psaut, qi, qs, den
1244 !-------------------------------------------------------------
1245 ! psaut: conversion(aggregation) of ice to snow [HDC 12]
1246 ! (T<T0: I->S)
1247 !-------------------------------------------------------------
1248  if(qi(i,k).gt.0.) then
1249  qimax = roqimax/den(i,k)
1250  psaut(i,k) = max(0.,(qi(i,k)-qimax)/dtcld)
1251  endif
1252 !
1253 !-------------------------------------------------------------
1254 ! pgaut: conversion(aggregation) of snow to graupel [HL A4] [LFO 37]
1255 ! (T<T0: S->G)
1256 !-------------------------------------------------------------
1257  if(qs(i,k).gt.0.) then
1258  alpha2 = 1.e-3*exp(0.09*(-supcol))
1259  pgaut(i,k) = min(max(0.,alpha2*(qs(i,k)-qs0)),qs(i,k)/dtcld)
1260  endif
1261  endif
1262 !
1263 ! WSM6-F90 TAG: PSEVP
1264 ! legacy_group: G13j
1265 ! process: Graupel evaporation/sublimation
1266 ! compare_vars: psevp, qg, qv, qsat, den
1267 !-------------------------------------------------------------
1268 ! psevp: evaporation of melting snow [HL A35] [RH83 A27]
1269 ! (T>=T0: S->V)
1270 !-------------------------------------------------------------
1271  if(supcol.lt.0.) then
1272  if(qs(i,k).gt.0..and.rh(i,k,1).lt.1.) then
1273  coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
1274  psevp(i,k) = (rh(i,k,1)-1.)*n0sfac(i,k)*(precs1 &
1275  * rslope2(i,k,2)+precs2*work2(i,k) &
1276  * coeres)/work1(i,k,1)
1277  psevp(i,k) = min(max(psevp(i,k),-qs(i,k)/dtcld),0.)
1278  endif
1279 !-------------------------------------------------------------
1280 ! pgevp: evaporation of melting graupel [HL A25] [RH84 A19]
1281 ! (T>=T0: G->V)
1282 !-------------------------------------------------------------
1283  if(qg(i,k).gt.0..and.rh(i,k,1).lt.1.) then
1284  coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3))
1285  pgevp(i,k) = (rh(i,k,1)-1.)*(precg1*rslope2(i,k,3) &
1286  + precg2*work2(i,k)*coeres)/work1(i,k,1)
1287  pgevp(i,k) = min(max(pgevp(i,k),-qg(i,k)/dtcld),0.)
1288  endif
1289  endif
1290  enddo
1291  enddo
1292 !
1293 !
1294 !----------------------------------------------------------------
1295 ! check mass conservation of generation terms and feedback to the
1296 ! large scale
1297 !
1298  ! WSM6-F90 TAG: UPDATE
1299  ! legacy_group: G14
1300  ! process: Mass conservation and state update
1301  ! compare_vars: t, q, qci, qrs, qv
1302  do k = kts, kte
1303  do i = its, ite
1304 !
1305  delta2=0.
1306  delta3=0.
1307  if(qr(i,k).lt.1.e-4.and.qs(i,k).lt.1.e-4) delta2=1.
1308  if(qr(i,k).lt.1.e-4) delta3=1.
1309  if(t(i,k).le.t0c) then
1310 !
1311 ! cloud water
1312 !
1313  value = max(qmin,qc(i,k))
1314  source = (praut(i,k)+pracw(i,k)+paacw(i,k)+paacw(i,k))*dtcld
1315  if (source.gt.value) then
1316  factor = value/source
1317  praut(i,k) = praut(i,k)*factor
1318  pracw(i,k) = pracw(i,k)*factor
1319  paacw(i,k) = paacw(i,k)*factor
1320  endif
1321 !
1322 ! cloud ice
1323 !
1324  value = max(qmin,qi(i,k))
1325  source = (psaut(i,k)-pigen(i,k)-pidep(i,k)+praci(i,k)+psaci(i,k) &
1326  + pgaci(i,k))*dtcld
1327  if (source.gt.value) then
1328  factor = value/source
1329  psaut(i,k) = psaut(i,k)*factor
1330  pigen(i,k) = pigen(i,k)*factor
1331  pidep(i,k) = pidep(i,k)*factor
1332  praci(i,k) = praci(i,k)*factor
1333  psaci(i,k) = psaci(i,k)*factor
1334  pgaci(i,k) = pgaci(i,k)*factor
1335  endif
1336 !
1337 ! rain
1338 !
1339  value = max(qmin,qr(i,k))
1340  source = (-praut(i,k)-prevp(i,k)-pracw(i,k)+piacr(i,k)+psacr(i,k) &
1341  + pgacr(i,k))*dtcld
1342  if (source.gt.value) then
1343  factor = value/source
1344  praut(i,k) = praut(i,k)*factor
1345  prevp(i,k) = prevp(i,k)*factor
1346  pracw(i,k) = pracw(i,k)*factor
1347  piacr(i,k) = piacr(i,k)*factor
1348  psacr(i,k) = psacr(i,k)*factor
1349  pgacr(i,k) = pgacr(i,k)*factor
1350  endif
1351 !
1352 ! snow
1353 !
1354  value = max(qmin,qs(i,k))
1355  source = -(psdep(i,k)+psaut(i,k)-pgaut(i,k)+paacw(i,k)+piacr(i,k) &
1356  * delta3+praci(i,k)*delta3-pracs(i,k)*(1.-delta2) &
1357  + psacr(i,k)*delta2+psaci(i,k)-pgacs(i,k) )*dtcld
1358  if (source.gt.value) then
1359  factor = value/source
1360  psdep(i,k) = psdep(i,k)*factor
1361  psaut(i,k) = psaut(i,k)*factor
1362  pgaut(i,k) = pgaut(i,k)*factor
1363  paacw(i,k) = paacw(i,k)*factor
1364  piacr(i,k) = piacr(i,k)*factor
1365  praci(i,k) = praci(i,k)*factor
1366  psaci(i,k) = psaci(i,k)*factor
1367  pracs(i,k) = pracs(i,k)*factor
1368  psacr(i,k) = psacr(i,k)*factor
1369  pgacs(i,k) = pgacs(i,k)*factor
1370  endif
1371 !
1372 ! graupel
1373 !
1374  value = max(qmin,qg(i,k))
1375  source = -(pgdep(i,k)+pgaut(i,k) &
1376  + piacr(i,k)*(1.-delta3)+praci(i,k)*(1.-delta3) &
1377  + psacr(i,k)*(1.-delta2)+pracs(i,k)*(1.-delta2) &
1378  + pgaci(i,k)+paacw(i,k)+pgacr(i,k)+pgacs(i,k))*dtcld
1379  if (source.gt.value) then
1380  factor = value/source
1381  pgdep(i,k) = pgdep(i,k)*factor
1382  pgaut(i,k) = pgaut(i,k)*factor
1383  piacr(i,k) = piacr(i,k)*factor
1384  praci(i,k) = praci(i,k)*factor
1385  psacr(i,k) = psacr(i,k)*factor
1386  pracs(i,k) = pracs(i,k)*factor
1387  paacw(i,k) = paacw(i,k)*factor
1388  pgaci(i,k) = pgaci(i,k)*factor
1389  pgacr(i,k) = pgacr(i,k)*factor
1390  pgacs(i,k) = pgacs(i,k)*factor
1391  endif
1392 !
1393  work2(i,k)=-(prevp(i,k)+psdep(i,k)+pgdep(i,k)+pigen(i,k)+pidep(i,k))
1394 ! update
1395  q(i,k) = q(i,k)+work2(i,k)*dtcld
1396  qc(i,k) = max(qc(i,k)-(praut(i,k)+pracw(i,k) &
1397  + paacw(i,k)+paacw(i,k))*dtcld,0.)
1398  qr(i,k) = max(qr(i,k)+(praut(i,k)+pracw(i,k) &
1399  + prevp(i,k)-piacr(i,k)-pgacr(i,k) &
1400  - psacr(i,k))*dtcld,0.)
1401  qi(i,k) = max(qi(i,k)-(psaut(i,k)+praci(i,k) &
1402  + psaci(i,k)+pgaci(i,k)-pigen(i,k)-pidep(i,k)) &
1403  * dtcld,0.)
1404  qs(i,k) = max(qs(i,k)+(psdep(i,k)+psaut(i,k)+paacw(i,k) &
1405  - pgaut(i,k)+piacr(i,k)*delta3 &
1406  + praci(i,k)*delta3+psaci(i,k)-pgacs(i,k) &
1407  - pracs(i,k)*(1.-delta2)+psacr(i,k)*delta2) &
1408  * dtcld,0.)
1409  qg(i,k) = max(qg(i,k)+(pgdep(i,k)+pgaut(i,k) &
1410  + piacr(i,k)*(1.-delta3) &
1411  + praci(i,k)*(1.-delta3)+psacr(i,k)*(1.-delta2) &
1412  + pracs(i,k)*(1.-delta2)+pgaci(i,k)+paacw(i,k) &
1413  + pgacr(i,k)+pgacs(i,k))*dtcld,0.)
1414  xlf = xls-xl(i,k)
1415  xlwork2 = -xls*(psdep(i,k)+pgdep(i,k)+pidep(i,k)+pigen(i,k)) &
1416  -xl(i,k)*prevp(i,k)-xlf*(piacr(i,k)+paacw(i,k) &
1417  +paacw(i,k)+pgacr(i,k)+psacr(i,k))
1418  t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
1419  else
1420 !
1421 ! cloud water
1422 !
1423  value = max(qmin,qc(i,k))
1424  source=(praut(i,k)+pracw(i,k)+paacw(i,k)+paacw(i,k))*dtcld
1425  if (source.gt.value) then
1426  factor = value/source
1427  praut(i,k) = praut(i,k)*factor
1428  pracw(i,k) = pracw(i,k)*factor
1429  paacw(i,k) = paacw(i,k)*factor
1430  endif
1431 !
1432 ! rain
1433 !
1434  value = max(qmin,qr(i,k))
1435  source = (-paacw(i,k)-praut(i,k)+pseml(i,k)+pgeml(i,k)-pracw(i,k) &
1436  -paacw(i,k)-prevp(i,k))*dtcld
1437  if (source.gt.value) then
1438  factor = value/source
1439  praut(i,k) = praut(i,k)*factor
1440  prevp(i,k) = prevp(i,k)*factor
1441  pracw(i,k) = pracw(i,k)*factor
1442  paacw(i,k) = paacw(i,k)*factor
1443  pseml(i,k) = pseml(i,k)*factor
1444  pgeml(i,k) = pgeml(i,k)*factor
1445  endif
1446 !
1447 ! snow
1448 !
1449  value = max(qcrmin,qs(i,k))
1450  source=(pgacs(i,k)-pseml(i,k)-psevp(i,k))*dtcld
1451  if (source.gt.value) then
1452  factor = value/source
1453  pgacs(i,k) = pgacs(i,k)*factor
1454  psevp(i,k) = psevp(i,k)*factor
1455  pseml(i,k) = pseml(i,k)*factor
1456  endif
1457 !
1458 ! graupel
1459 !
1460  value = max(qcrmin,qg(i,k))
1461  source=-(pgacs(i,k)+pgevp(i,k)+pgeml(i,k))*dtcld
1462  if (source.gt.value) then
1463  factor = value/source
1464  pgacs(i,k) = pgacs(i,k)*factor
1465  pgevp(i,k) = pgevp(i,k)*factor
1466  pgeml(i,k) = pgeml(i,k)*factor
1467  endif
1468 !
1469  work2(i,k)=-(prevp(i,k)+psevp(i,k)+pgevp(i,k))
1470 ! update
1471  q(i,k) = q(i,k)+work2(i,k)*dtcld
1472  qc(i,k) = max(qc(i,k)-(praut(i,k)+pracw(i,k) &
1473  + paacw(i,k)+paacw(i,k))*dtcld,0.)
1474  qr(i,k) = max(qr(i,k)+(praut(i,k)+pracw(i,k) &
1475  + prevp(i,k)+paacw(i,k)+paacw(i,k)-pseml(i,k) &
1476  - pgeml(i,k))*dtcld,0.)
1477  qs(i,k) = max(qs(i,k)+(psevp(i,k)-pgacs(i,k) &
1478  + pseml(i,k))*dtcld,0.)
1479  qg(i,k) = max(qg(i,k)+(pgacs(i,k)+pgevp(i,k) &
1480  + pgeml(i,k))*dtcld,0.)
1481  xlf = xls-xl(i,k)
1482  xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k)+pgevp(i,k)) &
1483  -xlf*(pseml(i,k)+pgeml(i,k))
1484  t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
1485  endif
1486  enddo
1487  enddo
1488 !
1489 ! WSM6-F90 TAG: QSAT2
1490 ! legacy_group: G15
1491 ! process: Second saturation mixing ratio computation
1492 ! compare_vars: qs, qvs, den, denfac, t, p
1493 ! Inline expansion for fpvs
1494 ! qsat(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1495 ! qsat(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1496  hsub = xls
1497  hvap = xlv0
1498  cvap = cpv
1499  ttp=t0c+0.01
1500  dldt=cvap-cliq
1501  xa=-dldt/rv
1502  xb=xa+hvap/(rv*ttp)
1503  dldti=cvap-cice
1504  xai=-dldti/rv
1505  xbi=xai+hsub/(rv*ttp)
1506  do k = kts, kte
1507  do i = its, ite
1508  tr=ttp/t(i,k)
1509  qsat(i,k,1)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
1510  qsat(i,k,1) = min(qsat(i,k,1),0.99*p(i,k))
1511  qsat(i,k,1) = ep2 * qsat(i,k,1) / (p(i,k) - qsat(i,k,1))
1512  qsat(i,k,1) = max(qsat(i,k,1),qmin)
1513  tr=ttp/t(i,k)
1514  if(t(i,k).lt.ttp) then
1515  qsat(i,k,2)=psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
1516  else
1517  qsat(i,k,2)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
1518  endif
1519  qsat(i,k,2) = min(qsat(i,k,2),0.99*p(i,k))
1520  qsat(i,k,2) = ep2 * qsat(i,k,2) / (p(i,k) - qsat(i,k,2))
1521  qsat(i,k,2) = max(qsat(i,k,2),qmin)
1522  enddo
1523  enddo
1524 !
1525 !----------------------------------------------------------------
1526 ! pcond: condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
1527 ! if there exists additional water vapor condensated/if
1528 ! evaporation of cloud water is not enough to remove subsaturation
1529 !
1530  ! WSM6-F90 TAG: PCOND
1531  ! legacy_group: G16
1532  ! process: Condensation/evaporation update
1533  ! compare_vars: pcond, t, qv, qc, qsat
1534  do k = kts, kte
1535  do i = its, ite
1536  work1(i,k,1) = conden(t(i,k),q(i,k),qsat(i,k,1),xl(i,k),cpm(i,k))
1537  work2(i,k) = qc(i,k)+work1(i,k,1)
1538  pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k),0.)/dtcld)
1539  if(qc(i,k).gt.0..and.work1(i,k,1).lt.0.) &
1540  pcond(i,k) = max(work1(i,k,1),-qc(i,k))/dtcld
1541  q(i,k) = q(i,k)-pcond(i,k)*dtcld
1542  qc(i,k) = max(qc(i,k)+pcond(i,k)*dtcld,0.)
1543  t(i,k) = t(i,k)+pcond(i,k)*xl(i,k)/cpm(i,k)*dtcld
1544  enddo
1545  enddo
1546 !
1547 !
1548 !----------------------------------------------------------------
1549 ! padding for small values
1550 !
1551  ! WSM6-F90 TAG: CLIP
1552  ! legacy_group: G17
1553  ! process: Padding/clipping for small values
1554  ! compare_vars: qv, qc, qr, qi, qs, qg
1555  do k = kts, kte
1556  do i = its, ite
1557  if(qc(i,k).le.qmin) qc(i,k) = 0.0
1558  if(qi(i,k).le.qmin) qi(i,k) = 0.0
1559  enddo
1560  enddo
1561  enddo ! big loops
1562 
1563  if(present(rainprod2d) .and. present(evapprod2d)) then
1564  do k = kts, kte
1565  do i = its,ite
1566  rainprod2d(i,k) = praut(i,k)+pracw(i,k)+praci(i,k)+psaci(i,k)+pgaci(i,k) &
1567  + psacw(i,k)+pgacw(i,k)+paacw(i,k)+psaut(i,k)
1568  evapprod2d(i,k) = -(prevp(i,k)+psevp(i,k)+pgevp(i,k)+psdep(i,k)+pgdep(i,k))
1569  enddo
1570  enddo
1571  endif
1572 !
1573 !----------------------------------------------------------------
1574 ! CCPP checks:
1575 !
1576 
1577  errmsg = 'mp_wsm6_run OK'
1578  errflg = 0
1579 

Referenced by mp_wsm6_isohelper::mp_wsm6_run_c().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ nislfv_rain_plm()

subroutine mp_wsm6::nislfv_rain_plm ( integer, intent(in)  im,
integer, intent(in)  km,
real(kind=kind_phys), dimension(im,km), intent(in)  denl,
real(kind=kind_phys), dimension(im,km), intent(in)  denfacl,
real(kind=kind_phys), dimension(im,km), intent(in)  tkl,
real(kind=kind_phys), dimension(im,km), intent(in)  dzl,
real(kind=kind_phys), dimension(im,km), intent(inout)  wwl,
real(kind=kind_phys), dimension(im,km), intent(inout)  rql,
real(kind=kind_phys), dimension(im), intent(inout)  precip,
real(kind=kind_phys), intent(in)  dt,
integer, intent(in)  id,
integer, intent(in)  iter 
)
private
1864 !=================================================================================================================
1865 !
1866 ! for non-iteration semi-Lagrangain forward advection for cloud
1867 ! with mass conservation and positive definite advection
1868 ! 2nd order interpolation with monotonic piecewise linear method
1869 ! this routine is under assumption of decfl < 1 for semi_Lagrangian
1870 !
1871 ! dzl depth of model layer in meter
1872 ! wwl terminal velocity at model layer m/s
1873 ! rql cloud density*mixing ration
1874 ! precip precipitation
1875 ! dt time step
1876 ! id kind of precip: 0 test case; 1 raindrop
1877 ! iter how many time to guess mean terminal velocity: 0 pure forward.
1878 ! 0 : use departure wind for advection
1879 ! 1 : use mean wind for advection
1880 ! > 1 : use mean wind after iter-1 iterations
1881 !
1882 ! author: hann-ming henry juang <henry.juang@noaa.gov>
1883 ! implemented by song-you hong
1884 !
1885 
1886 !--- input arguments:
1887 integer,intent(in):: im,km,id,iter
1888 
1889  real(kind=kind_phys),intent(in):: dt
1890  real(kind=kind_phys),intent(in),dimension(im,km):: dzl,denl,denfacl,tkl
1891 
1892 !--- inout arguments:
1893  real(kind=kind_phys),intent(inout),dimension(im):: precip
1894  real(kind=kind_phys),intent(inout),dimension(im,km):: rql,wwl
1895 
1896 !---- local variables and arrays:
1897 integer:: i,k,n,m,kk,kb,kt
1898 real(kind=kind_phys):: tl,tl2,qql,dql,qqd
1899 real(kind=kind_phys):: th,th2,qqh,dqh
1900 real(kind=kind_phys):: zsum,qsum,dim,dip,c1,con1,fa1,fa2
1901 real(kind=kind_phys):: allold,allnew,zz,dzamin,cflmax,decfl
1902 real(kind=kind_phys),dimension(km):: dz,ww,qq,wd,wa,was
1903 real(kind=kind_phys),dimension(km):: den,denfac,tk
1904 real(kind=kind_phys),dimension(km):: qn,qr,tmp,tmp1,tmp2,tmp3
1905 real(kind=kind_phys),dimension(km+1):: wi,zi,za
1906 real(kind=kind_phys),dimension(km+1):: dza,qa,qmi,qpi
1907 
1908 !-----------------------------------------------------------------------------------------------------------------
1909 
1910  precip(:) = 0.0
1911 
1912  i_loop: do i=1,im
1913  dz(:) = dzl(i,:)
1914  qq(:) = rql(i,:)
1915  ww(:) = wwl(i,:)
1916  den(:) = denl(i,:)
1917  denfac(:) = denfacl(i,:)
1918  tk(:) = tkl(i,:)
1919 ! skip for no precipitation for all layers
1920  allold = 0.0
1921  do k=1,km
1922  allold = allold + qq(k)
1923  enddo
1924  if(allold.le.0.0) then
1925  cycle i_loop
1926  endif
1927 !
1928 ! compute interface values
1929  zi(1)=0.0
1930  do k=1,km
1931  zi(k+1) = zi(k)+dz(k)
1932  enddo
1933 !
1934 ! save departure wind
1935  wd(:) = ww(:)
1936  n=1
1937  100 continue
1938 ! plm is 2nd order, we can use 2nd order wi or 3rd order wi
1939 ! 2nd order interpolation to get wi
1940  wi(1) = ww(1)
1941  wi(km+1) = ww(km)
1942  do k=2,km
1943  wi(k) = (ww(k)*dz(k-1)+ww(k-1)*dz(k))/(dz(k-1)+dz(k))
1944  enddo
1945 ! 3rd order interpolation to get wi
1946  fa1 = 9./16.
1947  fa2 = 1./16.
1948  wi(1) = ww(1)
1949  wi(2) = 0.5*(ww(2)+ww(1))
1950  do k=3,km-1
1951  wi(k) = fa1*(ww(k)+ww(k-1))-fa2*(ww(k+1)+ww(k-2))
1952  enddo
1953  wi(km) = 0.5*(ww(km)+ww(km-1))
1954  wi(km+1) = ww(km)
1955 !
1956 ! terminate of top of raingroup
1957  do k=2,km
1958  if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
1959  enddo
1960 !
1961 ! diffusivity of wi
1962  con1 = 0.05
1963  do k=km,1,-1
1964  decfl = (wi(k+1)-wi(k))*dt/dz(k)
1965  if( decfl .gt. con1 ) then
1966  wi(k) = wi(k+1) - con1*dz(k)/dt
1967  endif
1968  enddo
1969 ! compute arrival point
1970  do k=1,km+1
1971  za(k) = zi(k) - wi(k)*dt
1972  enddo
1973 !
1974  do k=1,km
1975  dza(k) = za(k+1)-za(k)
1976  enddo
1977  dza(km+1) = zi(km+1) - za(km+1)
1978 !
1979 ! computer deformation at arrival point
1980  do k=1,km
1981  qa(k) = qq(k)*dz(k)/dza(k)
1982  qr(k) = qa(k)/den(k)
1983  enddo
1984  qa(km+1) = 0.0
1985 ! call maxmin(km,1,qa,' arrival points ')
1986 !
1987 ! compute arrival terminal velocity, and estimate mean terminal velocity
1988 ! then back to use mean terminal velocity
1989  if( n.le.iter ) then
1990  call slope_rain(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
1991  if( n.ge.2 ) wa(1:km)=0.5*(wa(1:km)+was(1:km))
1992  do k=1,km
1993 !#ifdef DEBUG
1994 ! print*,' slope_wsm3 ',qr(k)*1000.,den(k),denfac(k),tk(k),tmp(k),tmp1(k),tmp2(k),ww(k),wa(k)
1995 !#endif
1996 ! mean wind is average of departure and new arrival winds
1997  ww(k) = 0.5* ( wd(k)+wa(k) )
1998  enddo
1999  was(:) = wa(:)
2000  n=n+1
2001  go to 100
2002  endif
2003 !
2004 ! estimate values at arrival cell interface with monotone
2005  do k=2,km
2006  dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k))
2007  dim=(qa(k)-qa(k-1))/(dza(k-1)+dza(k))
2008  if( dip*dim.le.0.0 ) then
2009  qmi(k)=qa(k)
2010  qpi(k)=qa(k)
2011  else
2012  qpi(k)=qa(k)+0.5*(dip+dim)*dza(k)
2013  qmi(k)=2.0*qa(k)-qpi(k)
2014  if( qpi(k).lt.0.0 .or. qmi(k).lt.0.0 ) then
2015  qpi(k) = qa(k)
2016  qmi(k) = qa(k)
2017  endif
2018  endif
2019  enddo
2020  qpi(1)=qa(1)
2021  qmi(1)=qa(1)
2022  qmi(km+1)=qa(km+1)
2023  qpi(km+1)=qa(km+1)
2024 !
2025 ! interpolation to regular point
2026  qn = 0.0
2027  kb=1
2028  kt=1
2029  intp : do k=1,km
2030  kb=max(kb-1,1)
2031  kt=max(kt-1,1)
2032 ! find kb and kt
2033  if( zi(k).ge.za(km+1) ) then
2034  exit intp
2035  else
2036  find_kb : do kk=kb,km
2037  if( zi(k).le.za(kk+1) ) then
2038  kb = kk
2039  exit find_kb
2040  else
2041  cycle find_kb
2042  endif
2043  enddo find_kb
2044  find_kt : do kk=kt,km
2045  if( zi(k+1).le.za(kk) ) then
2046  kt = kk
2047  exit find_kt
2048  else
2049  cycle find_kt
2050  endif
2051  enddo find_kt
2052  kt = kt - 1
2053 ! compute q with piecewise constant method
2054  if( kt.eq.kb ) then
2055  tl=(zi(k)-za(kb))/dza(kb)
2056  th=(zi(k+1)-za(kb))/dza(kb)
2057  tl2=tl*tl
2058  th2=th*th
2059  qqd=0.5*(qpi(kb)-qmi(kb))
2060  qqh=qqd*th2+qmi(kb)*th
2061  qql=qqd*tl2+qmi(kb)*tl
2062  qn(k) = (qqh-qql)/(th-tl)
2063  zsum = dza(kb)
2064  qsum = qn(k)*zsum
2065  else if( kt.gt.kb ) then
2066  tl=(zi(k)-za(kb))/dza(kb)
2067  tl2=tl*tl
2068  qqd=0.5*(qpi(kb)-qmi(kb))
2069  qql=qqd*tl2+qmi(kb)*tl
2070  dql = qa(kb)-qql
2071  zsum = (1.-tl)*dza(kb)
2072  qsum = dql*dza(kb)
2073  if( kt-kb.gt.1 ) then
2074  do m=kb+1,kt-1
2075  zsum = zsum + dza(m)
2076  qsum = qsum + qa(m) * dza(m)
2077  enddo
2078  endif
2079  th=(zi(k+1)-za(kt))/dza(kt)
2080  th2=th*th
2081  qqd=0.5*(qpi(kt)-qmi(kt))
2082  dqh=qqd*th2+qmi(kt)*th
2083  zsum = zsum + th*dza(kt)
2084  qsum = qsum + dqh*dza(kt)
2085  qn(k) = qsum/zsum
2086  endif
2087  cycle intp
2088  endif
2089 !
2090  enddo intp
2091 !
2092 ! rain out
2093  sum_precip: do k=1,km
2094  if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then
2095  precip(i) = precip(i) + qa(k)*dza(k)
2096  cycle sum_precip
2097  else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then
2098  precip(i) = precip(i) + qa(k)*(0.0-za(k))
2099  exit sum_precip
2100  endif
2101  exit sum_precip
2102  enddo sum_precip
2103 !
2104 ! replace the new values
2105  rql(i,:) = qn(:)
2106  enddo i_loop
2107 

Referenced by mp_wsm6_run().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ nislfv_rain_plm6()

subroutine mp_wsm6::nislfv_rain_plm6 ( integer, intent(in)  im,
integer, intent(in)  km,
real(kind=kind_phys), dimension(im,km), intent(in)  denl,
real(kind=kind_phys), dimension(im,km), intent(in)  denfacl,
real(kind=kind_phys), dimension(im,km), intent(in)  tkl,
real(kind=kind_phys), dimension(im,km), intent(in)  dzl,
real(kind=kind_phys), dimension(im,km), intent(inout)  wwl,
real(kind=kind_phys), dimension(im,km), intent(inout)  rql,
real(kind=kind_phys), dimension(im,km), intent(inout)  rql2,
real(kind=kind_phys), dimension(im), intent(inout)  precip1,
real(kind=kind_phys), dimension(im), intent(inout)  precip2,
real(kind=kind_phys), intent(in)  dt,
integer, intent(in)  id,
integer, intent(in)  iter 
)
private
2112 !=================================================================================================================
2113 !
2114 ! for non-iteration semi-Lagrangain forward advection for cloud
2115 ! with mass conservation and positive definite advection
2116 ! 2nd order interpolation with monotonic piecewise linear method
2117 ! this routine is under assumption of decfl < 1 for semi_Lagrangian
2118 !
2119 ! dzl depth of model layer in meter
2120 ! wwl terminal velocity at model layer m/s
2121 ! rql cloud density*mixing ration
2122 ! precip precipitation
2123 ! dt time step
2124 ! id kind of precip: 0 test case; 1 raindrop
2125 ! iter how many time to guess mean terminal velocity: 0 pure forward.
2126 ! 0 : use departure wind for advection
2127 ! 1 : use mean wind for advection
2128 ! > 1 : use mean wind after iter-1 iterations
2129 !
2130 ! author: hann-ming henry juang <henry.juang@noaa.gov>
2131 ! implemented by song-you hong
2132 !
2133 
2134 !--- input arguments:
2135  integer,intent(in):: im,km,id,iter
2136 
2137  real(kind=kind_phys),intent(in):: dt
2138  real(kind=kind_phys),intent(in),dimension(im,km):: dzl,denl,denfacl,tkl
2139 
2140 !--- inout arguments:
2141  real(kind=kind_phys),intent(inout),dimension(im):: precip1,precip2
2142  real(kind=kind_phys),intent(inout),dimension(im,km):: rql,rql2,wwl
2143 
2144 !---- local variables and arrays:
2145  integer:: i,ist,k,n,m,kk,kb,kt
2146  real(kind=kind_phys):: tl,tl2,qql,dql,qqd
2147  real(kind=kind_phys):: th,th2,qqh,dqh
2148  real(kind=kind_phys):: zsum,qsum,dim,dip,c1,con1,fa1,fa2
2149  real(kind=kind_phys):: allold,allnew,zz,dzamin,cflmax,decfl
2150  real(kind=kind_phys),dimension(km):: dz,ww,qq,qq2,wd,wa,wa2,was
2151  real(kind=kind_phys),dimension(km):: den,denfac,tk
2152  real(kind=kind_phys),dimension(km):: qn,qr,qr2,tmp,tmp1,tmp2,tmp3
2153  real(kind=kind_phys),dimension(km+1):: wi,zi,za
2154  real(kind=kind_phys),dimension(km+1):: dza,qa,qa2,qmi,qpi
2155  real(kind=kind_phys),dimension(im):: precip
2156 
2157 !-----------------------------------------------------------------------------------------------------------------
2158 
2159  precip(:) = 0.0
2160  precip1(:) = 0.0
2161  precip2(:) = 0.0
2162 
2163  i_loop: do i=1,im
2164  dz(:) = dzl(i,:)
2165  qq(:) = rql(i,:)
2166  qq2(:) = rql2(i,:)
2167  ww(:) = wwl(i,:)
2168  den(:) = denl(i,:)
2169  denfac(:) = denfacl(i,:)
2170  tk(:) = tkl(i,:)
2171 ! skip for no precipitation for all layers
2172  allold = 0.0
2173  do k=1,km
2174  allold = allold + qq(k) + qq2(k)
2175  enddo
2176  if(allold.le.0.0) then
2177  cycle i_loop
2178  endif
2179 !
2180 ! compute interface values
2181  zi(1)=0.0
2182  do k=1,km
2183  zi(k+1) = zi(k)+dz(k)
2184  enddo
2185 !
2186 ! save departure wind
2187  wd(:) = ww(:)
2188  n=1
2189  100 continue
2190 ! plm is 2nd order, we can use 2nd order wi or 3rd order wi
2191 ! 2nd order interpolation to get wi
2192  wi(1) = ww(1)
2193  wi(km+1) = ww(km)
2194  do k=2,km
2195  wi(k) = (ww(k)*dz(k-1)+ww(k-1)*dz(k))/(dz(k-1)+dz(k))
2196  enddo
2197 ! 3rd order interpolation to get wi
2198  fa1 = 9./16.
2199  fa2 = 1./16.
2200  wi(1) = ww(1)
2201  wi(2) = 0.5*(ww(2)+ww(1))
2202  do k=3,km-1
2203  wi(k) = fa1*(ww(k)+ww(k-1))-fa2*(ww(k+1)+ww(k-2))
2204  enddo
2205  wi(km) = 0.5*(ww(km)+ww(km-1))
2206  wi(km+1) = ww(km)
2207 !
2208 ! terminate of top of raingroup
2209  do k=2,km
2210  if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
2211  enddo
2212 !
2213 ! diffusivity of wi
2214  con1 = 0.05
2215  do k=km,1,-1
2216  decfl = (wi(k+1)-wi(k))*dt/dz(k)
2217  if( decfl .gt. con1 ) then
2218  wi(k) = wi(k+1) - con1*dz(k)/dt
2219  endif
2220  enddo
2221 ! compute arrival point
2222  do k=1,km+1
2223  za(k) = zi(k) - wi(k)*dt
2224  enddo
2225 !
2226  do k=1,km
2227  dza(k) = za(k+1)-za(k)
2228  enddo
2229  dza(km+1) = zi(km+1) - za(km+1)
2230 !
2231 ! computer deformation at arrival point
2232  do k=1,km
2233  qa(k) = qq(k)*dz(k)/dza(k)
2234  qa2(k) = qq2(k)*dz(k)/dza(k)
2235  qr(k) = qa(k)/den(k)
2236  qr2(k) = qa2(k)/den(k)
2237  enddo
2238  qa(km+1) = 0.0
2239  qa2(km+1) = 0.0
2240 ! call maxmin(km,1,qa,' arrival points ')
2241 !
2242 ! compute arrival terminal velocity, and estimate mean terminal velocity
2243 ! then back to use mean terminal velocity
2244  if( n.le.iter ) then
2245  call slope_snow(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
2246  call slope_graup(qr2,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa2,1,1,1,km)
2247  do k = 1, km
2248  tmp(k) = max((qr(k)+qr2(k)), 1.e-15)
2249  if( tmp(k) .gt. 1.e-15 ) then
2250  wa(k) = (wa(k)*qr(k) + wa2(k)*qr2(k))/tmp(k)
2251  else
2252  wa(k) = 0.
2253  endif
2254  enddo
2255  if( n.ge.2 ) wa(1:km)=0.5*(wa(1:km)+was(1:km))
2256  do k=1,km
2257 !#ifdef DEBUG
2258 ! print*,' slope_wsm3 ',qr(k)*1000.,den(k),denfac(k),tk(k),tmp(k),tmp1(k),tmp2(k), &
2259 ! ww(k),wa(k)
2260 !#endif
2261 ! mean wind is average of departure and new arrival winds
2262  ww(k) = 0.5* ( wd(k)+wa(k) )
2263  enddo
2264  was(:) = wa(:)
2265  n=n+1
2266  go to 100
2267  endif
2268 
2269  ist_loop : do ist = 1, 2
2270  if (ist.eq.2) then
2271  qa(:) = qa2(:)
2272  endif
2273 !
2274  precip(i) = 0.
2275 !
2276 ! estimate values at arrival cell interface with monotone
2277  do k=2,km
2278  dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k))
2279  dim=(qa(k)-qa(k-1))/(dza(k-1)+dza(k))
2280  if( dip*dim.le.0.0 ) then
2281  qmi(k)=qa(k)
2282  qpi(k)=qa(k)
2283  else
2284  qpi(k)=qa(k)+0.5*(dip+dim)*dza(k)
2285  qmi(k)=2.0*qa(k)-qpi(k)
2286  if( qpi(k).lt.0.0 .or. qmi(k).lt.0.0 ) then
2287  qpi(k) = qa(k)
2288  qmi(k) = qa(k)
2289  endif
2290  endif
2291  enddo
2292  qpi(1)=qa(1)
2293  qmi(1)=qa(1)
2294  qmi(km+1)=qa(km+1)
2295  qpi(km+1)=qa(km+1)
2296 !
2297 ! interpolation to regular point
2298  qn = 0.0
2299  kb=1
2300  kt=1
2301  intp : do k=1,km
2302  kb=max(kb-1,1)
2303  kt=max(kt-1,1)
2304 ! find kb and kt
2305  if( zi(k).ge.za(km+1) ) then
2306  exit intp
2307  else
2308  find_kb : do kk=kb,km
2309  if( zi(k).le.za(kk+1) ) then
2310  kb = kk
2311  exit find_kb
2312  else
2313  cycle find_kb
2314  endif
2315  enddo find_kb
2316  find_kt : do kk=kt,km
2317  if( zi(k+1).le.za(kk) ) then
2318  kt = kk
2319  exit find_kt
2320  else
2321  cycle find_kt
2322  endif
2323  enddo find_kt
2324  kt = kt - 1
2325 ! compute q with piecewise constant method
2326  if( kt.eq.kb ) then
2327  tl=(zi(k)-za(kb))/dza(kb)
2328  th=(zi(k+1)-za(kb))/dza(kb)
2329  tl2=tl*tl
2330  th2=th*th
2331  qqd=0.5*(qpi(kb)-qmi(kb))
2332  qqh=qqd*th2+qmi(kb)*th
2333  qql=qqd*tl2+qmi(kb)*tl
2334  qn(k) = (qqh-qql)/(th-tl)
2335  else if( kt.gt.kb ) then
2336  tl=(zi(k)-za(kb))/dza(kb)
2337  tl2=tl*tl
2338  qqd=0.5*(qpi(kb)-qmi(kb))
2339  qql=qqd*tl2+qmi(kb)*tl
2340  dql = qa(kb)-qql
2341  zsum = (1.-tl)*dza(kb)
2342  qsum = dql*dza(kb)
2343  if( kt-kb.gt.1 ) then
2344  do m=kb+1,kt-1
2345  zsum = zsum + dza(m)
2346  qsum = qsum + qa(m) * dza(m)
2347  enddo
2348  endif
2349  th=(zi(k+1)-za(kt))/dza(kt)
2350  th2=th*th
2351  qqd=0.5*(qpi(kt)-qmi(kt))
2352  dqh=qqd*th2+qmi(kt)*th
2353  zsum = zsum + th*dza(kt)
2354  qsum = qsum + dqh*dza(kt)
2355  qn(k) = qsum/zsum
2356  endif
2357  cycle intp
2358  endif
2359 !
2360  enddo intp
2361 !
2362 ! rain out
2363  sum_precip: do k=1,km
2364  if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then
2365  precip(i) = precip(i) + qa(k)*dza(k)
2366  cycle sum_precip
2367  else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then
2368  precip(i) = precip(i) + qa(k)*(0.0-za(k))
2369  exit sum_precip
2370  endif
2371  exit sum_precip
2372  enddo sum_precip
2373 !
2374 ! replace the new values
2375  if(ist.eq.1) then
2376  rql(i,:) = qn(:)
2377  precip1(i) = precip(i)
2378  else
2379  rql2(i,:) = qn(:)
2380  precip2(i) = precip(i)
2381  endif
2382  enddo ist_loop
2383 
2384  enddo i_loop
2385 

Referenced by mp_wsm6_run().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ refl10cm_wsm6()

subroutine, public mp_wsm6::refl10cm_wsm6 ( real(kind=kind_phys), dimension(kts:kte), intent(in)  qv1d,
real(kind=kind_phys), dimension(kts:kte), intent(in)  qr1d,
real(kind=kind_phys), dimension(kts:kte), intent(in)  qs1d,
real(kind=kind_phys), dimension(kts:kte), intent(in)  qg1d,
real(kind=kind_phys), dimension(kts:kte), intent(in)  t1d,
real(kind=kind_phys), dimension(kts:kte), intent(in)  p1d,
real(kind=kind_phys), dimension(kts:kte), intent(inout)  dBZ,
integer, intent(in)  kts,
integer, intent(in)  kte 
)
2390  implicit none
2391 !=================================================================================================================
2392 
2393 !..Sub arguments
2394  integer,intent(in):: kts,kte
2395  real(kind=kind_phys),intent(in),dimension(kts:kte):: qv1d,qr1d,qs1d,qg1d,t1d,p1d
2396  real(kind=kind_phys),intent(inout),dimension(kts:kte):: dbz
2397 
2398 !..Local variables
2399  logical:: melti
2400  logical,dimension(kts:kte):: l_qr,l_qs,l_qg
2401 
2402  INTEGER:: i,k,k_0,kbot,n
2403 
2404  real(kind=kind_phys),parameter:: r=287.
2405  real(kind=kind_phys):: temp_c
2406  real(kind=kind_phys),dimension(kts:kte):: temp,pres,qv,rho
2407  real(kind=kind_phys),dimension(kts:kte):: rr,rs,rg
2408  real(kind=kind_phys),dimension(kts:kte):: ze_rain,ze_snow,ze_graupel
2409 
2410  double precision:: fmelt_s,fmelt_g
2411  double precision:: cback,x,eta,f_d
2412  double precision,dimension(kts:kte):: ilamr,ilams,ilamg
2413  double precision,dimension(kts:kte):: n0_r, n0_s, n0_g
2414  double precision:: lamr,lams,lamg
2415 
2416 !-----------------------------------------------------------------------------------------------------------------
2417 
2418  do k = kts, kte
2419  dbz(k) = -35.0
2420  enddo
2421 
2422 !+---+-----------------------------------------------------------------+
2423 !..Put column of data into local arrays.
2424 !+---+-----------------------------------------------------------------+
2425  do k = kts, kte
2426  temp(k) = t1d(k)
2427  temp_c = min(-0.001, temp(k)-273.15)
2428  qv(k) = max(1.e-10, qv1d(k))
2429  pres(k) = p1d(k)
2430  rho(k) = 0.622*pres(k)/(r*temp(k)*(qv(k)+0.622))
2431 
2432  if (qr1d(k) .gt. 1.e-9) then
2433  rr(k) = qr1d(k)*rho(k)
2434  n0_r(k) = n0r
2435  lamr = (xam_r*xcrg(3)*n0_r(k)/rr(k))**(1./xcre(1))
2436  ilamr(k) = 1./lamr
2437  l_qr(k) = .true.
2438  else
2439  rr(k) = 1.e-12
2440  l_qr(k) = .false.
2441  endif
2442 
2443  if (qs1d(k) .gt. 1.e-9) then
2444  rs(k) = qs1d(k)*rho(k)
2445  n0_s(k) = min(n0smax, n0s*exp(-alpha*temp_c))
2446  lams = (xam_s*xcsg(3)*n0_s(k)/rs(k))**(1./xcse(1))
2447  ilams(k) = 1./lams
2448  l_qs(k) = .true.
2449  else
2450  rs(k) = 1.e-12
2451  l_qs(k) = .false.
2452  endif
2453 
2454  if (qg1d(k) .gt. 1.e-9) then
2455  rg(k) = qg1d(k)*rho(k)
2456  n0_g(k) = n0g
2457  lamg = (xam_g*xcgg(3)*n0_g(k)/rg(k))**(1./xcge(1))
2458  ilamg(k) = 1./lamg
2459  l_qg(k) = .true.
2460  else
2461  rg(k) = 1.e-12
2462  l_qg(k) = .false.
2463  endif
2464  enddo
2465 
2466 !+---+-----------------------------------------------------------------+
2467 !..Locate K-level of start of melting (k_0 is level above).
2468 !+---+-----------------------------------------------------------------+
2469  melti = .false.
2470  k_0 = kts
2471  do k = kte-1, kts, -1
2472  if ( (temp(k).gt.273.15) .and. l_qr(k) &
2473  .and. (l_qs(k+1).or.l_qg(k+1)) ) then
2474  k_0 = max(k+1, k_0)
2475  melti=.true.
2476  goto 195
2477  endif
2478  enddo
2479  195 continue
2480 
2481 !+---+-----------------------------------------------------------------+
2482 !..Assume Rayleigh approximation at 10 cm wavelength. Rain (all temps)
2483 !.. and non-water-coated snow and graupel when below freezing are
2484 !.. simple. Integrations of m(D)*m(D)*N(D)*dD.
2485 !+---+-----------------------------------------------------------------+
2486 
2487  do k = kts, kte
2488  ze_rain(k) = 1.e-22
2489  ze_snow(k) = 1.e-22
2490  ze_graupel(k) = 1.e-22
2491  if (l_qr(k)) ze_rain(k) = n0_r(k)*xcrg(4)*ilamr(k)**xcre(4)
2492  if (l_qs(k)) ze_snow(k) = (0.176/0.93) * (6.0/pi)*(6.0/pi) &
2493  * (xam_s/900.0)*(xam_s/900.0) &
2494  * n0_s(k)*xcsg(4)*ilams(k)**xcse(4)
2495  if (l_qg(k)) ze_graupel(k) = (0.176/0.93) * (6.0/pi)*(6.0/pi) &
2496  * (xam_g/900.0)*(xam_g/900.0) &
2497  * n0_g(k)*xcgg(4)*ilamg(k)**xcge(4)
2498  enddo
2499 
2500 
2501 !+---+-----------------------------------------------------------------+
2502 !..Special case of melting ice (snow/graupel) particles. Assume the
2503 !.. ice is surrounded by the liquid water. Fraction of meltwater is
2504 !.. extremely simple based on amount found above the melting level.
2505 !.. Uses code from Uli Blahak (rayleigh_soak_wetgraupel and supporting
2506 !.. routines).
2507 !+---+-----------------------------------------------------------------+
2508 
2509  if (melti .and. k_0.ge.kts+1) then
2510  do k = k_0-1, kts, -1
2511 
2512 !..Reflectivity contributed by melting snow
2513  if (l_qs(k) .and. l_qs(k_0) ) then
2514  fmelt_s = max(0.005d0, min(1.0d0-rs(k)/rs(k_0), 0.99d0))
2515  eta = 0.d0
2516  lams = 1./ilams(k)
2517  do n = 1, nrbins
2518  x = xam_s * xxds(n)**xbm_s
2519  call rayleigh_soak_wetgraupel (x,dble(xocms),dble(xobms), &
2520  fmelt_s, melt_outside_s, m_w_0, m_i_0, lamda_radar, &
2521  cback, mixingrulestring_s, matrixstring_s, &
2522  inclusionstring_s, hoststring_s, &
2523  hostmatrixstring_s, hostinclusionstring_s)
2524  f_d = n0_s(k)*xxds(n)**xmu_s * dexp(-lams*xxds(n))
2525  eta = eta + f_d * cback * simpson(n) * xdts(n)
2526  enddo
2527  ze_snow(k) = sngl(lamda4 / (pi5 * k_w) * eta)
2528  endif
2529 
2530 
2531 !..Reflectivity contributed by melting graupel
2532 
2533  if (l_qg(k) .and. l_qg(k_0) ) then
2534  fmelt_g = max(0.005d0, min(1.0d0-rg(k)/rg(k_0), 0.99d0))
2535  eta = 0.d0
2536  lamg = 1./ilamg(k)
2537  do n = 1, nrbins
2538  x = xam_g * xxdg(n)**xbm_g
2539  call rayleigh_soak_wetgraupel (x,dble(xocmg),dble(xobmg), &
2540  fmelt_g, melt_outside_g, m_w_0, m_i_0, lamda_radar, &
2541  cback, mixingrulestring_g, matrixstring_g, &
2542  inclusionstring_g, hoststring_g, &
2543  hostmatrixstring_g, hostinclusionstring_g)
2544  f_d = n0_g(k)*xxdg(n)**xmu_g * dexp(-lamg*xxdg(n))
2545  eta = eta + f_d * cback * simpson(n) * xdtg(n)
2546  enddo
2547  ze_graupel(k) = sngl(lamda4 / (pi5 * k_w) * eta)
2548  endif
2549 
2550  enddo
2551  endif
2552 
2553  do k = kte, kts, -1
2554  dbz(k) = 10.*log10((ze_rain(k)+ze_snow(k)+ze_graupel(k))*1.d18)
2555  enddo
2556 
2557 
Here is the call graph for this function:

◆ rgmma()

real(kind=kind_phys) function mp_wsm6::rgmma ( real(kind=kind_phys), intent(in)  x)
private
1584 !=================================================================================================================
1585 !rgmma function: use infinite product form
1586 
1587  real(kind=kind_phys),intent(in):: x
1588 
1589  integer:: i
1590  real(kind=kind_phys),parameter:: euler=0.577215664901532
1591  real(kind=kind_phys):: y
1592 
1593 !-----------------------------------------------------------------------------------------------------------------
1594 
1595  if(x.eq.1.)then
1596  rgmma=0.
1597  else
1598  rgmma=x*exp(euler*x)
1599  do i = 1,10000
1600  y = float(i)
1601  rgmma=rgmma*(1.000+x/y)*exp(-x/y)
1602  enddo
1603  rgmma=1./rgmma
1604  endif
1605 

Referenced by WSM6::initialize_coeffs(), mp_wsm6_init(), wsm6_nislfv_rain_plm6_scratch(), and wsm6_nislfv_rain_plm_scratch().

Here is the caller graph for this function:

◆ slope_graup()

subroutine mp_wsm6::slope_graup ( real(kind=kind_phys), dimension(its:ite,kts:kte), intent(in)  qrs,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(in)  den,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(in)  denfac,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(in)  t,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(inout)  rslope,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(inout)  rslopeb,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(inout)  rslope2,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(inout)  rslope3,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(inout)  vt,
integer  its,
integer  ite,
integer  kts,
integer  kte 
)
private
1815 !=================================================================================================================
1816 
1817 !--- input arguments:
1818  integer:: its,ite,kts,kte
1819 
1820  real(kind=kind_phys),intent(in),dimension(its:ite,kts:kte):: qrs,den,denfac,t
1821 
1822 !--- inout arguments:
1823  real(kind=kind_phys),intent(inout),dimension(its:ite,kts:kte):: &
1824  rslope,rslopeb,rslope2,rslope3,vt
1825 
1826 !--- local variables and arrays:
1827  integer:: i,k
1828 
1829  real(kind=kind_phys),parameter:: t0c = 273.15
1830  real(kind=kind_phys):: lamdag,x,y
1831  real(kind=kind_phys),dimension(its:ite,kts:kte):: n0sfac
1832 
1833 !-----------------------------------------------------------------------------------------------------------------
1834 
1835 !size distributions: (x=mixing ratio, y=air density):
1836 !valid for mixing ratio > 1.e-9 kg/kg.
1837  lamdag(x,y)= sqrt(sqrt(pidn0g/(x*y))) ! (pidn0g/(x*y))**.25
1838 
1839  do k = kts, kte
1840  do i = its, ite
1841 !---------------------------------------------------------------
1842 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
1843 !---------------------------------------------------------------
1844  if(qrs(i,k).le.qcrmin)then
1845  rslope(i,k) = rslopegmax
1846  rslopeb(i,k) = rslopegbmax
1847  rslope2(i,k) = rslopeg2max
1848  rslope3(i,k) = rslopeg3max
1849  else
1850  rslope(i,k) = 1./lamdag(qrs(i,k),den(i,k))
1851  rslopeb(i,k) = rslope(i,k)**bvtg
1852  rslope2(i,k) = rslope(i,k)*rslope(i,k)
1853  rslope3(i,k) = rslope2(i,k)*rslope(i,k)
1854  endif
1855  vt(i,k) = pvtg*rslopeb(i,k)*denfac(i,k)
1856  if(qrs(i,k).le.0.0) vt(i,k) = 0.0
1857  enddo
1858  enddo
1859 

Referenced by nislfv_rain_plm6().

Here is the caller graph for this function:

◆ slope_rain()

subroutine mp_wsm6::slope_rain ( real(kind=kind_phys), dimension(its:ite,kts:kte), intent(in)  qrs,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(in)  den,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(in)  denfac,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(in)  t,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(inout)  rslope,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(inout)  rslopeb,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(inout)  rslope2,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(inout)  rslope3,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(inout)  vt,
integer  its,
integer  ite,
integer  kts,
integer  kte 
)
private
1718 !=================================================================================================================
1719 
1720 !--- input arguments:
1721  integer:: its,ite,kts,kte
1722 
1723  real(kind=kind_phys),intent(in),dimension(its:ite,kts:kte):: qrs,den,denfac,t
1724 
1725 !--- inout arguments:
1726  real(kind=kind_phys),intent(inout),dimension(its:ite,kts:kte):: &
1727  rslope,rslopeb,rslope2,rslope3,vt
1728 
1729 !--- local variables and arrays:
1730  integer:: i,k
1731 
1732  real(kind=kind_phys),parameter:: t0c = 273.15
1733  real(kind=kind_phys):: lamdar,x,y
1734  real(kind=kind_phys),dimension(its:ite,kts:kte):: n0sfac
1735 
1736 !-----------------------------------------------------------------------------------------------------------------
1737 
1738 !size distributions: (x=mixing ratio, y=air density):
1739 !valid for mixing ratio > 1.e-9 kg/kg.
1740  lamdar(x,y)= sqrt(sqrt(pidn0r/(x*y))) ! (pidn0r/(x*y))**.25
1741 
1742  do k = kts, kte
1743  do i = its, ite
1744  if(qrs(i,k).le.qcrmin)then
1745  rslope(i,k) = rslopermax
1746  rslopeb(i,k) = rsloperbmax
1747  rslope2(i,k) = rsloper2max
1748  rslope3(i,k) = rsloper3max
1749  else
1750  rslope(i,k) = 1./lamdar(qrs(i,k),den(i,k))
1751  rslopeb(i,k) = rslope(i,k)**bvtr
1752  rslope2(i,k) = rslope(i,k)*rslope(i,k)
1753  rslope3(i,k) = rslope2(i,k)*rslope(i,k)
1754  endif
1755  vt(i,k) = pvtr*rslopeb(i,k)*denfac(i,k)
1756  if(qrs(i,k).le.0.0) vt(i,k) = 0.0
1757  enddo
1758  enddo
1759 

Referenced by nislfv_rain_plm().

Here is the caller graph for this function:

◆ slope_snow()

subroutine mp_wsm6::slope_snow ( real(kind=kind_phys), dimension(its:ite,kts:kte), intent(in)  qrs,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(in)  den,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(in)  denfac,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(in)  t,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(inout)  rslope,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(inout)  rslopeb,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(inout)  rslope2,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(inout)  rslope3,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(inout)  vt,
integer  its,
integer  ite,
integer  kts,
integer  kte 
)
private
1764 !=================================================================================================================
1765 
1766 !--- input arguments:
1767  integer:: its,ite,kts,kte
1768 
1769  real(kind=kind_phys),intent(in),dimension(its:ite,kts:kte):: qrs,den,denfac,t
1770 
1771 !--- inout arguments:
1772  real(kind=kind_phys),intent(inout),dimension(its:ite,kts:kte):: &
1773  rslope,rslopeb,rslope2,rslope3,vt
1774 
1775 !--- local variables and arrays:
1776  integer:: i,k
1777 
1778  real(kind=kind_phys),parameter:: t0c = 273.15
1779  real(kind=kind_phys):: lamdas,x,y,z,supcol
1780  real(kind=kind_phys),dimension(its:ite,kts:kte):: n0sfac
1781 
1782 !-----------------------------------------------------------------------------------------------------------------
1783 
1784 !size distributions: (x=mixing ratio, y=air density):
1785 !valid for mixing ratio > 1.e-9 kg/kg.
1786  lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y))) ! (pidn0s*z/(x*y))**.25
1787 !
1788  do k = kts, kte
1789  do i = its, ite
1790  supcol = t0c-t(i,k)
1791 !---------------------------------------------------------------
1792 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
1793 !---------------------------------------------------------------
1794  n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
1795  if(qrs(i,k).le.qcrmin)then
1796  rslope(i,k) = rslopesmax
1797  rslopeb(i,k) = rslopesbmax
1798  rslope2(i,k) = rslopes2max
1799  rslope3(i,k) = rslopes3max
1800  else
1801  rslope(i,k) = 1./lamdas(qrs(i,k),den(i,k),n0sfac(i,k))
1802  rslopeb(i,k) = rslope(i,k)**bvts
1803  rslope2(i,k) = rslope(i,k)*rslope(i,k)
1804  rslope3(i,k) = rslope2(i,k)*rslope(i,k)
1805  endif
1806  vt(i,k) = pvts*rslopeb(i,k)*denfac(i,k)
1807  if(qrs(i,k).le.0.0) vt(i,k) = 0.0
1808  enddo
1809  enddo
1810 

Referenced by nislfv_rain_plm6().

Here is the caller graph for this function:

◆ slope_wsm6()

subroutine mp_wsm6::slope_wsm6 ( real(kind=kind_phys), dimension(its:ite,kts:kte,3), intent(in)  qrs,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(in)  den,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(in)  denfac,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(in)  t,
real(kind=kind_phys), dimension(its:ite,kts:kte,3), intent(inout)  rslope,
real(kind=kind_phys), dimension(its:ite,kts:kte,3), intent(inout)  rslopeb,
real(kind=kind_phys), dimension(its:ite,kts:kte,3), intent(inout)  rslope2,
real(kind=kind_phys), dimension(its:ite,kts:kte,3), intent(inout)  rslope3,
real(kind=kind_phys), dimension(its:ite,kts:kte,3), intent(inout)  vt,
integer  its,
integer  ite,
integer  kts,
integer  kte 
)
private
1638 !=================================================================================================================
1639 
1640 !--- input arguments:
1641  integer:: its,ite,kts,kte
1642 
1643  real(kind=kind_phys),intent(in),dimension(its:ite,kts:kte):: den,denfac,t
1644  real(kind=kind_phys),intent(in),dimension(its:ite,kts:kte,3):: qrs
1645 
1646 !--- inout arguments:
1647  real(kind=kind_phys),intent(inout),dimension(its:ite,kts:kte,3):: &
1648  rslope,rslopeb,rslope2,rslope3,vt
1649 
1650 !--- local variables and arrays:
1651  integer:: i,k
1652 
1653  real(kind=kind_phys),parameter:: t0c = 273.15
1654  real(kind=kind_phys):: lamdar,lamdas,lamdag,x,y,z,supcol
1655  real(kind=kind_phys),dimension(its:ite,kts:kte):: n0sfac
1656 
1657 !-----------------------------------------------------------------------------------------------------------------
1658 
1659 !size distributions: (x=mixing ratio, y=air density):
1660 !valid for mixing ratio > 1.e-9 kg/kg.
1661  lamdar(x,y)= sqrt(sqrt(pidn0r/(x*y))) ! (pidn0r/(x*y))**.25
1662  lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y))) ! (pidn0s*z/(x*y))**.25
1663  lamdag(x,y)= sqrt(sqrt(pidn0g/(x*y))) ! (pidn0g/(x*y))**.25
1664 
1665  do k = kts, kte
1666  do i = its, ite
1667  supcol = t0c-t(i,k)
1668 !---------------------------------------------------------------
1669 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
1670 !---------------------------------------------------------------
1671  n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
1672  if(qrs(i,k,1).le.qcrmin)then
1673  rslope(i,k,1) = rslopermax
1674  rslopeb(i,k,1) = rsloperbmax
1675  rslope2(i,k,1) = rsloper2max
1676  rslope3(i,k,1) = rsloper3max
1677  else
1678  rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k))
1679  rslopeb(i,k,1) = rslope(i,k,1)**bvtr
1680  rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
1681  rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
1682  endif
1683  if(qrs(i,k,2).le.qcrmin)then
1684  rslope(i,k,2) = rslopesmax
1685  rslopeb(i,k,2) = rslopesbmax
1686  rslope2(i,k,2) = rslopes2max
1687  rslope3(i,k,2) = rslopes3max
1688  else
1689  rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
1690  rslopeb(i,k,2) = rslope(i,k,2)**bvts
1691  rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
1692  rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
1693  endif
1694  if(qrs(i,k,3).le.qcrmin)then
1695  rslope(i,k,3) = rslopegmax
1696  rslopeb(i,k,3) = rslopegbmax
1697  rslope2(i,k,3) = rslopeg2max
1698  rslope3(i,k,3) = rslopeg3max
1699  else
1700  rslope(i,k,3) = 1./lamdag(qrs(i,k,3),den(i,k))
1701  rslopeb(i,k,3) = rslope(i,k,3)**bvtg
1702  rslope2(i,k,3) = rslope(i,k,3)*rslope(i,k,3)
1703  rslope3(i,k,3) = rslope2(i,k,3)*rslope(i,k,3)
1704  endif
1705  vt(i,k,1) = pvtr*rslopeb(i,k,1)*denfac(i,k)
1706  vt(i,k,2) = pvts*rslopeb(i,k,2)*denfac(i,k)
1707  vt(i,k,3) = pvtg*rslopeb(i,k,3)*denfac(i,k)
1708  if(qrs(i,k,1).le.0.0) vt(i,k,1) = 0.0
1709  if(qrs(i,k,2).le.0.0) vt(i,k,2) = 0.0
1710  if(qrs(i,k,3).le.0.0) vt(i,k,3) = 0.0
1711  enddo
1712  enddo
1713 

Referenced by mp_wsm6_run().

Here is the caller graph for this function:

Variable Documentation

◆ alpha

real(kind=kind_phys), parameter, public mp_wsm6::alpha = .12
44  real(kind=kind_phys),parameter,public :: alpha = .12 ! .122 exponen factor for n0s

Referenced by compute_max_reflectivity_dbz(), SurfaceLayer::fill_tsurf_with_sst_and_tsk(), interpolate_1d(), mp_wsm6_run(), realbdy_compute_interior_ghost_rhs(), refl10cm_wsm6(), slope_snow(), and slope_wsm6().

◆ avtg

real(kind=kind_phys), save mp_wsm6::avtg
private

Referenced by mp_wsm6_init(), and wsm6_nislfv_rain_plm6_scratch().

◆ avtr

real(kind=kind_phys), parameter, private mp_wsm6::avtr = 841.9
private
19  real(kind=kind_phys),parameter,private:: avtr = 841.9 ! a constant for terminal velocity of rain

Referenced by WSM6::Advance(), WSM6::initialize_coeffs(), mp_wsm6_init(), and mp_wsm6_run().

◆ avts

real(kind=kind_phys), parameter, private mp_wsm6::avts = 11.72
private
25  real(kind=kind_phys),parameter,private:: avts = 11.72 ! a constant for terminal velocity of snow

Referenced by WSM6::initialize_coeffs(), and mp_wsm6_init().

◆ bvtg

real(kind=kind_phys), save mp_wsm6::bvtg
private

Referenced by mp_wsm6_init(), slope_graup(), slope_wsm6(), and wsm6_nislfv_rain_plm6_scratch().

◆ bvtg1

real(kind=kind_phys), save mp_wsm6::bvtg1
private

Referenced by mp_wsm6_init().

◆ bvtg2

real(kind=kind_phys), save mp_wsm6::bvtg2
private

Referenced by mp_wsm6_init().

◆ bvtg3

real(kind=kind_phys), save mp_wsm6::bvtg3
private

Referenced by mp_wsm6_init().

◆ bvtg4

real(kind=kind_phys), save mp_wsm6::bvtg4
private

Referenced by mp_wsm6_init().

◆ bvtr

real(kind=kind_phys), parameter, private mp_wsm6::bvtr = 0.8
private
20  real(kind=kind_phys),parameter,private:: bvtr = 0.8 ! a constant for terminal velocity of rain

Referenced by WSM6::Advance(), WSM6::initialize_coeffs(), mp_wsm6_init(), slope_rain(), and slope_wsm6().

◆ bvtr1

real(kind=kind_phys), save mp_wsm6::bvtr1
private

Referenced by mp_wsm6_init().

◆ bvtr2

real(kind=kind_phys), save mp_wsm6::bvtr2
private

Referenced by mp_wsm6_init().

◆ bvtr3

real(kind=kind_phys), save mp_wsm6::bvtr3
private

Referenced by mp_wsm6_init().

◆ bvtr4

real(kind=kind_phys), save mp_wsm6::bvtr4
private

Referenced by mp_wsm6_init().

◆ bvtr6

real(kind=kind_phys), save mp_wsm6::bvtr6
private

Referenced by mp_wsm6_init().

◆ bvts

real(kind=kind_phys), parameter, private mp_wsm6::bvts = .41
private
26  real(kind=kind_phys),parameter,private:: bvts = .41 ! a constant for terminal velocity of snow

Referenced by WSM6::Advance(), WSM6::initialize_coeffs(), mp_wsm6_init(), slope_snow(), and slope_wsm6().

◆ bvts1

real(kind=kind_phys), save mp_wsm6::bvts1
private

Referenced by mp_wsm6_init().

◆ bvts2

real(kind=kind_phys), save mp_wsm6::bvts2
private

Referenced by mp_wsm6_init().

◆ bvts3

real(kind=kind_phys), save mp_wsm6::bvts3
private

Referenced by mp_wsm6_init().

◆ bvts4

real(kind=kind_phys), save mp_wsm6::bvts4
private

Referenced by mp_wsm6_init().

◆ deng

real(kind=kind_phys), save mp_wsm6::deng
private

Referenced by mp_wsm6_init(), and wsm6_nislfv_rain_plm6_scratch().

◆ dens

real(kind=kind_phys), parameter, private mp_wsm6::dens = 100.0
private
39  real(kind=kind_phys),parameter,private:: dens = 100.0 ! Density of snow

Referenced by WSM6::Advance(), ERF::erf_enforce_hse(), ERF::MakeHorizontalAverages(), mp_wsm6_init(), mp_wsm6_run(), and wsm6_nislfv_rain_plm6_scratch().

◆ dicon

real(kind=kind_phys), parameter, private mp_wsm6::dicon = 11.9
private
33  real(kind=kind_phys),parameter,private:: dicon = 11.9 ! constant for the cloud-ice diamter

Referenced by WSM6::Advance(), and mp_wsm6_run().

◆ dimax

real(kind=kind_phys), parameter, private mp_wsm6::dimax = 500.e-6
private
34  real(kind=kind_phys),parameter,private:: dimax = 500.e-6 ! limited maximum value for the cloud-ice diamter

Referenced by WSM6::Advance(), WSM6::initialize_coeffs(), mp_wsm6_init(), and mp_wsm6_run().

◆ dtcldcr

real(kind=kind_phys), parameter, private mp_wsm6::dtcldcr = 120.
private
16  real(kind=kind_phys),parameter,private:: dtcldcr = 120. ! maximum time step for minor loops

Referenced by WSM6::Advance(), and mp_wsm6_run().

◆ eacrc

real(kind=kind_phys), parameter, private mp_wsm6::eacrc = 1.0
private
38  real(kind=kind_phys),parameter,private:: eacrc = 1.0 ! Snow/cloud-water collection efficiency

Referenced by WSM6::initialize_coeffs(), and mp_wsm6_init().

◆ eacrr

real(kind=kind_phys), save mp_wsm6::eacrr
private

Referenced by mp_wsm6_init().

◆ g1pbg

real(kind=kind_phys), save mp_wsm6::g1pbg
private

Referenced by mp_wsm6_init().

◆ g1pbr

real(kind=kind_phys), save mp_wsm6::g1pbr
private

Referenced by mp_wsm6_init().

◆ g1pbs

real(kind=kind_phys), save mp_wsm6::g1pbs
private

Referenced by mp_wsm6_init().

◆ g3pbg

real(kind=kind_phys), save mp_wsm6::g3pbg
private

Referenced by mp_wsm6_init().

◆ g3pbr

real(kind=kind_phys), save mp_wsm6::g3pbr
private

Referenced by mp_wsm6_init().

◆ g3pbs

real(kind=kind_phys), save mp_wsm6::g3pbs
private

Referenced by mp_wsm6_init().

◆ g4pbg

real(kind=kind_phys), save mp_wsm6::g4pbg
private

Referenced by mp_wsm6_init().

◆ g4pbr

real(kind=kind_phys), save mp_wsm6::g4pbr
private

Referenced by mp_wsm6_init().

◆ g4pbs

real(kind=kind_phys), save mp_wsm6::g4pbs
private

Referenced by mp_wsm6_init().

◆ g5pbgo2

real(kind=kind_phys), save mp_wsm6::g5pbgo2
private

Referenced by mp_wsm6_init().

◆ g5pbro2

real(kind=kind_phys), save mp_wsm6::g5pbro2
private

Referenced by mp_wsm6_init().

◆ g5pbso2

real(kind=kind_phys), save mp_wsm6::g5pbso2
private

Referenced by mp_wsm6_init().

◆ g6pbr

real(kind=kind_phys), save mp_wsm6::g6pbr
private

Referenced by WSM6::Advance(), mp_wsm6_init(), and mp_wsm6_run().

◆ kind_phys

integer, parameter mp_wsm6::kind_phys = c_double
9  integer, parameter :: kind_phys = c_double

◆ lamdagmax

real(kind=kind_phys), save mp_wsm6::lamdagmax
private

Referenced by mp_wsm6_init(), and wsm6_nislfv_rain_plm6_scratch().

◆ lamdarmax

real(kind=kind_phys), parameter, private mp_wsm6::lamdarmax = 8.e4
private
30  real(kind=kind_phys),parameter,private:: lamdarmax = 8.e4 ! limited maximum value for slope parameter of rain

Referenced by WSM6::initialize_coeffs(), and mp_wsm6_init().

◆ lamdasmax

real(kind=kind_phys), parameter, private mp_wsm6::lamdasmax = 1.e5
private
31  real(kind=kind_phys),parameter,private:: lamdasmax = 1.e5 ! limited maximum value for slope parameter of snow

Referenced by WSM6::initialize_coeffs(), and mp_wsm6_init().

◆ n0g

real(kind=kind_phys), save mp_wsm6::n0g
private

Referenced by mp_wsm6_init(), mp_wsm6_run(), and refl10cm_wsm6().

◆ n0r

real(kind=kind_phys), parameter, private mp_wsm6::n0r = 8.e6
private
17  real(kind=kind_phys),parameter,private:: n0r = 8.e6 ! intercept parameter rain

Referenced by mp_wsm6_init(), mp_wsm6_run(), and refl10cm_wsm6().

◆ n0s

real(kind=kind_phys), parameter, public mp_wsm6::n0s = 2.e6
43  real(kind=kind_phys),parameter,public :: n0s = 2.e6 ! temperature dependent intercept parameter snow

Referenced by mp_wsm6_init(), mp_wsm6_run(), refl10cm_wsm6(), slope_snow(), and slope_wsm6().

◆ n0smax

real(kind=kind_phys), parameter, public mp_wsm6::n0smax = 1.e11
42  real(kind=kind_phys),parameter,public :: n0smax = 1.e11 ! maximum n0s (t=-90C unlimited)

Referenced by WSM6::Advance(), mp_wsm6_run(), refl10cm_wsm6(), slope_snow(), and slope_wsm6().

◆ pacrc

real(kind=kind_phys), save mp_wsm6::pacrc
private

Referenced by WSM6::Advance(), mp_wsm6_init(), and mp_wsm6_run().

◆ pacrg

real(kind=kind_phys), save mp_wsm6::pacrg
private

Referenced by WSM6::Advance(), mp_wsm6_init(), and mp_wsm6_run().

◆ pacrr

real(kind=kind_phys), save mp_wsm6::pacrr
private

Referenced by WSM6::Advance(), mp_wsm6_init(), and mp_wsm6_run().

◆ pacrs

real(kind=kind_phys), save mp_wsm6::pacrs
private

Referenced by mp_wsm6_init().

◆ peaut

real(kind=kind_phys), parameter, private mp_wsm6::peaut = .55
private
22  real(kind=kind_phys),parameter,private:: peaut = .55 ! collection efficiency

Referenced by WSM6::initialize_coeffs(), and mp_wsm6_init().

◆ pfrz1

real(kind=kind_phys), parameter, private mp_wsm6::pfrz1 = 100.
private
35  real(kind=kind_phys),parameter,private:: pfrz1 = 100. ! constant in Biggs freezing

Referenced by WSM6::Advance(), and mp_wsm6_run().

◆ pfrz2

real(kind=kind_phys), parameter, private mp_wsm6::pfrz2 = 0.66
private
36  real(kind=kind_phys),parameter,private:: pfrz2 = 0.66 ! constant in Biggs freezing

Referenced by WSM6::Advance(), and mp_wsm6_run().

◆ pi

real(kind=kind_phys), save mp_wsm6::pi
private

Referenced by mp_wsm6_init(), mp_wsm6_run(), and refl10cm_wsm6().

◆ pidn0g

real(kind=kind_phys), save mp_wsm6::pidn0g
private

Referenced by mp_wsm6_init(), slope_graup(), slope_wsm6(), and wsm6_nislfv_rain_plm6_scratch().

◆ pidn0r

real(kind=kind_phys), save mp_wsm6::pidn0r
private

Referenced by mp_wsm6_init(), slope_rain(), slope_wsm6(), and wsm6_nislfv_rain_plm_scratch().

◆ pidn0s

real(kind=kind_phys), save, public mp_wsm6::pidn0s
64  real(kind=kind_phys),public,save:: pidn0s,pidnc

Referenced by mp_wsm6_init(), slope_snow(), slope_wsm6(), and wsm6_nislfv_rain_plm6_scratch().

◆ pidnc

real(kind=kind_phys), save, public mp_wsm6::pidnc

Referenced by mp_wsm6_init().

◆ precg1

real(kind=kind_phys), save mp_wsm6::precg1
private

Referenced by WSM6::Advance(), mp_wsm6_init(), and mp_wsm6_run().

◆ precg2

real(kind=kind_phys), save mp_wsm6::precg2
private

Referenced by WSM6::Advance(), mp_wsm6_init(), and mp_wsm6_run().

◆ precr1

real(kind=kind_phys), save mp_wsm6::precr1
private

Referenced by WSM6::Advance(), mp_wsm6_init(), and mp_wsm6_run().

◆ precr2

real(kind=kind_phys), save mp_wsm6::precr2
private

Referenced by WSM6::Advance(), mp_wsm6_init(), and mp_wsm6_run().

◆ precs1

real(kind=kind_phys), save mp_wsm6::precs1
private

Referenced by WSM6::Advance(), mp_wsm6_init(), and mp_wsm6_run().

◆ precs2

real(kind=kind_phys), save mp_wsm6::precs2
private

Referenced by WSM6::Advance(), mp_wsm6_init(), and mp_wsm6_run().

◆ pvtg

real(kind=kind_phys), save mp_wsm6::pvtg
private

Referenced by WSM6::Advance(), mp_wsm6_init(), mp_wsm6_run(), slope_graup(), slope_wsm6(), and wsm6_nislfv_rain_plm6_scratch().

◆ pvtr

real(kind=kind_phys), save mp_wsm6::pvtr
private

Referenced by WSM6::Advance(), mp_wsm6_init(), mp_wsm6_run(), slope_rain(), slope_wsm6(), and wsm6_nislfv_rain_plm_scratch().

◆ pvts

real(kind=kind_phys), save mp_wsm6::pvts
private

Referenced by WSM6::Advance(), mp_wsm6_init(), mp_wsm6_run(), slope_snow(), slope_wsm6(), and wsm6_nislfv_rain_plm6_scratch().

◆ qc0

real(kind=kind_phys), save mp_wsm6::qc0
private
46  real(kind=kind_phys),save:: &
47  qc0,qck1, &
48  bvtr1,bvtr2,bvtr3,bvtr4,g1pbr, &
49  g3pbr,g4pbr,g5pbro2,pvtr,eacrr,pacrr, &
50  bvtr6,g6pbr, &
51  precr1,precr2,roqimax,bvts1, &
52  bvts2,bvts3,bvts4,g1pbs,g3pbs,g4pbs, &
53  n0g,avtg,bvtg,deng,lamdagmax, & !RAS13.3 - set these in wsm6init
54  g5pbso2,pvts,pacrs,precs1,precs2,pidn0r, &
55  xlv1,pacrc,pi, &
56  bvtg1,bvtg2,bvtg3,bvtg4,g1pbg, &
57  g3pbg,g4pbg,g5pbgo2,pvtg,pacrg, &
58  precg1,precg2,pidn0g, &
59  rslopermax,rslopesmax,rslopegmax, &
60  rsloperbmax,rslopesbmax,rslopegbmax, &
61  rsloper2max,rslopes2max,rslopeg2max, &
62  rsloper3max,rslopes3max,rslopeg3max

Referenced by WSM6::Advance(), ERF::derive_diag_profiles_stag(), mp_wsm6_init(), and mp_wsm6_run().

◆ qck1

real(kind=kind_phys), save mp_wsm6::qck1
private

Referenced by WSM6::Advance(), mp_wsm6_init(), and mp_wsm6_run().

◆ qcrmin

real(kind=kind_phys), parameter, private mp_wsm6::qcrmin = 1.e-9
private
37  real(kind=kind_phys),parameter,private:: qcrmin = 1.e-9 ! minimun values for qr, qs, and qg

Referenced by WSM6::Advance(), mp_wsm6_run(), slope_graup(), slope_rain(), slope_snow(), and slope_wsm6().

◆ qs0

real(kind=kind_phys), parameter, private mp_wsm6::qs0 = 6.e-4
private
40  real(kind=kind_phys),parameter,private:: qs0 = 6.e-4 ! threshold amount for aggretion to occur

Referenced by WSM6::Advance(), and mp_wsm6_run().

◆ r0

real(kind=kind_phys), parameter, private mp_wsm6::r0 = .8e-5
private
21  real(kind=kind_phys),parameter,private:: r0 = .8e-5 ! 8 microm in contrast to 10 micro m

Referenced by ERF::advance_dycore(), ApplySpongeZoneBCsForCC(), ApplySpongeZoneBCsForMom(), WSM6::initialize_coeffs(), make_buoyancy(), make_fast_coeffs(), make_mom_sources(), make_sources(), mp_wsm6_init(), and ParallelFor().

◆ roqimax

real(kind=kind_phys), save mp_wsm6::roqimax
private

Referenced by WSM6::Advance(), mp_wsm6_init(), and mp_wsm6_run().

◆ rslopeg2max

real(kind=kind_phys), save mp_wsm6::rslopeg2max
private

Referenced by mp_wsm6_init(), slope_graup(), slope_wsm6(), and wsm6_nislfv_rain_plm6_scratch().

◆ rslopeg3max

real(kind=kind_phys), save mp_wsm6::rslopeg3max
private

Referenced by mp_wsm6_init(), slope_graup(), slope_wsm6(), and wsm6_nislfv_rain_plm6_scratch().

◆ rslopegbmax

real(kind=kind_phys), save mp_wsm6::rslopegbmax
private

Referenced by mp_wsm6_init(), slope_graup(), slope_wsm6(), and wsm6_nislfv_rain_plm6_scratch().

◆ rslopegmax

real(kind=kind_phys), save mp_wsm6::rslopegmax
private

Referenced by mp_wsm6_init(), slope_graup(), slope_wsm6(), and wsm6_nislfv_rain_plm6_scratch().

◆ rsloper2max

real(kind=kind_phys), save mp_wsm6::rsloper2max
private

Referenced by mp_wsm6_init(), slope_rain(), slope_wsm6(), and wsm6_nislfv_rain_plm_scratch().

◆ rsloper3max

real(kind=kind_phys), save mp_wsm6::rsloper3max
private

Referenced by mp_wsm6_init(), slope_rain(), slope_wsm6(), and wsm6_nislfv_rain_plm_scratch().

◆ rsloperbmax

real(kind=kind_phys), save mp_wsm6::rsloperbmax
private

Referenced by mp_wsm6_init(), slope_rain(), slope_wsm6(), and wsm6_nislfv_rain_plm_scratch().

◆ rslopermax

real(kind=kind_phys), save mp_wsm6::rslopermax
private

Referenced by mp_wsm6_init(), slope_rain(), slope_wsm6(), and wsm6_nislfv_rain_plm_scratch().

◆ rslopes2max

real(kind=kind_phys), save mp_wsm6::rslopes2max
private

Referenced by mp_wsm6_init(), slope_snow(), slope_wsm6(), and wsm6_nislfv_rain_plm6_scratch().

◆ rslopes3max

real(kind=kind_phys), save mp_wsm6::rslopes3max
private

Referenced by mp_wsm6_init(), slope_snow(), slope_wsm6(), and wsm6_nislfv_rain_plm6_scratch().

◆ rslopesbmax

real(kind=kind_phys), save mp_wsm6::rslopesbmax
private

Referenced by mp_wsm6_init(), slope_snow(), slope_wsm6(), and wsm6_nislfv_rain_plm6_scratch().

◆ rslopesmax

real(kind=kind_phys), save mp_wsm6::rslopesmax
private

Referenced by mp_wsm6_init(), slope_snow(), slope_wsm6(), and wsm6_nislfv_rain_plm6_scratch().

◆ xlv1

real(kind=kind_phys), save mp_wsm6::xlv1
private

Referenced by mp_wsm6_init(), and mp_wsm6_run().

◆ xmyu

real(kind=kind_phys), parameter, private mp_wsm6::xmyu = 1.718e-5
private
24  real(kind=kind_phys),parameter,private:: xmyu = 1.718e-5 ! the dynamic viscosity kgm-1s-1

Referenced by WSM6::initialize_coeffs(), and mp_wsm6_init().

◆ xncr

real(kind=kind_phys), parameter, private mp_wsm6::xncr = 3.e8
private
23  real(kind=kind_phys),parameter,private:: xncr = 3.e8 ! maritime cloud in contrast to 3.e8 in tc80

Referenced by WSM6::initialize_coeffs(), mp_wsm6_init(), and mp_wsm6_run().