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, rainonly, snowonly, 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
1623 !=================================================================================================================
1624 
1625  integer,intent(in):: ice
1626  real(kind=kind_phys),intent(in):: cice,cliq,cvap,hsub,hvap,psat,rd,rv,t0c
1627  real(kind=kind_phys),intent(in):: t
1628 
1629  real(kind=kind_phys):: tr,ttp,dldt,dldti,xa,xb,xai,xbi
1630 
1631 !-----------------------------------------------------------------------------------------------------------------
1632 
1633  ttp=t0c+0.01
1634  dldt=cvap-cliq
1635  xa=-dldt/rv
1636  xb=xa+hvap/(rv*ttp)
1637  dldti=cvap-cice
1638  xai=-dldti/rv
1639  xbi=xai+hsub/(rv*ttp)
1640  tr=ttp/t
1641  if(t.lt.ttp.and.ice.eq.1) then
1642  fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
1643  else
1644  fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
1645  endif
1646 

◆ 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,
real(kind=kind_phys), dimension(its:), intent(inout), optional  rainonly,
real(kind=kind_phys), dimension(its:), intent(inout), optional  snowonly,
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

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

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

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

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

Referenced by mp_wsm6_init().

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
1828 !=================================================================================================================
1829 
1830 !--- input arguments:
1831  integer:: its,ite,kts,kte
1832 
1833  real(kind=kind_phys),intent(in),dimension(its:ite,kts:kte):: qrs,den,denfac,t
1834 
1835 !--- inout arguments:
1836  real(kind=kind_phys),intent(inout),dimension(its:ite,kts:kte):: &
1837  rslope,rslopeb,rslope2,rslope3,vt
1838 
1839 !--- local variables and arrays:
1840  integer:: i,k
1841 
1842  real(kind=kind_phys),parameter:: t0c = 273.15
1843  real(kind=kind_phys):: lamdag,x,y
1844  real(kind=kind_phys),dimension(its:ite,kts:kte):: n0sfac
1845 
1846 !-----------------------------------------------------------------------------------------------------------------
1847 
1848 !size distributions: (x=mixing ratio, y=air density):
1849 !valid for mixing ratio > 1.e-9 kg/kg.
1850  lamdag(x,y)= sqrt(sqrt(pidn0g/(x*y))) ! (pidn0g/(x*y))**.25
1851 
1852  do k = kts, kte
1853  do i = its, ite
1854 !---------------------------------------------------------------
1855 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
1856 !---------------------------------------------------------------
1857  if(qrs(i,k).le.qcrmin)then
1858  rslope(i,k) = rslopegmax
1859  rslopeb(i,k) = rslopegbmax
1860  rslope2(i,k) = rslopeg2max
1861  rslope3(i,k) = rslopeg3max
1862  else
1863  rslope(i,k) = 1./lamdag(qrs(i,k),den(i,k))
1864  rslopeb(i,k) = rslope(i,k)**bvtg
1865  rslope2(i,k) = rslope(i,k)*rslope(i,k)
1866  rslope3(i,k) = rslope2(i,k)*rslope(i,k)
1867  endif
1868  vt(i,k) = pvtg*rslopeb(i,k)*denfac(i,k)
1869  if(qrs(i,k).le.0.0) vt(i,k) = 0.0
1870  enddo
1871  enddo
1872 

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
1731 !=================================================================================================================
1732 
1733 !--- input arguments:
1734  integer:: its,ite,kts,kte
1735 
1736  real(kind=kind_phys),intent(in),dimension(its:ite,kts:kte):: qrs,den,denfac,t
1737 
1738 !--- inout arguments:
1739  real(kind=kind_phys),intent(inout),dimension(its:ite,kts:kte):: &
1740  rslope,rslopeb,rslope2,rslope3,vt
1741 
1742 !--- local variables and arrays:
1743  integer:: i,k
1744 
1745  real(kind=kind_phys),parameter:: t0c = 273.15
1746  real(kind=kind_phys):: lamdar,x,y
1747  real(kind=kind_phys),dimension(its:ite,kts:kte):: n0sfac
1748 
1749 !-----------------------------------------------------------------------------------------------------------------
1750 
1751 !size distributions: (x=mixing ratio, y=air density):
1752 !valid for mixing ratio > 1.e-9 kg/kg.
1753  lamdar(x,y)= sqrt(sqrt(pidn0r/(x*y))) ! (pidn0r/(x*y))**.25
1754 
1755  do k = kts, kte
1756  do i = its, ite
1757  if(qrs(i,k).le.qcrmin)then
1758  rslope(i,k) = rslopermax
1759  rslopeb(i,k) = rsloperbmax
1760  rslope2(i,k) = rsloper2max
1761  rslope3(i,k) = rsloper3max
1762  else
1763  rslope(i,k) = 1./lamdar(qrs(i,k),den(i,k))
1764  rslopeb(i,k) = rslope(i,k)**bvtr
1765  rslope2(i,k) = rslope(i,k)*rslope(i,k)
1766  rslope3(i,k) = rslope2(i,k)*rslope(i,k)
1767  endif
1768  vt(i,k) = pvtr*rslopeb(i,k)*denfac(i,k)
1769  if(qrs(i,k).le.0.0) vt(i,k) = 0.0
1770  enddo
1771  enddo
1772 

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
1777 !=================================================================================================================
1778 
1779 !--- input arguments:
1780  integer:: its,ite,kts,kte
1781 
1782  real(kind=kind_phys),intent(in),dimension(its:ite,kts:kte):: qrs,den,denfac,t
1783 
1784 !--- inout arguments:
1785  real(kind=kind_phys),intent(inout),dimension(its:ite,kts:kte):: &
1786  rslope,rslopeb,rslope2,rslope3,vt
1787 
1788 !--- local variables and arrays:
1789  integer:: i,k
1790 
1791  real(kind=kind_phys),parameter:: t0c = 273.15
1792  real(kind=kind_phys):: lamdas,x,y,z,supcol
1793  real(kind=kind_phys),dimension(its:ite,kts:kte):: n0sfac
1794 
1795 !-----------------------------------------------------------------------------------------------------------------
1796 
1797 !size distributions: (x=mixing ratio, y=air density):
1798 !valid for mixing ratio > 1.e-9 kg/kg.
1799  lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y))) ! (pidn0s*z/(x*y))**.25
1800 !
1801  do k = kts, kte
1802  do i = its, ite
1803  supcol = t0c-t(i,k)
1804 !---------------------------------------------------------------
1805 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
1806 !---------------------------------------------------------------
1807  n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
1808  if(qrs(i,k).le.qcrmin)then
1809  rslope(i,k) = rslopesmax
1810  rslopeb(i,k) = rslopesbmax
1811  rslope2(i,k) = rslopes2max
1812  rslope3(i,k) = rslopes3max
1813  else
1814  rslope(i,k) = 1./lamdas(qrs(i,k),den(i,k),n0sfac(i,k))
1815  rslopeb(i,k) = rslope(i,k)**bvts
1816  rslope2(i,k) = rslope(i,k)*rslope(i,k)
1817  rslope3(i,k) = rslope2(i,k)*rslope(i,k)
1818  endif
1819  vt(i,k) = pvts*rslopeb(i,k)*denfac(i,k)
1820  if(qrs(i,k).le.0.0) vt(i,k) = 0.0
1821  enddo
1822  enddo
1823 

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

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 mp_wsm6_run(), refl10cm_wsm6(), slope_snow(), and slope_wsm6().

◆ avtg

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

Referenced by mp_wsm6_init().

◆ 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 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 mp_wsm6_init().

◆ bvtg

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

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

◆ 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 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 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().

◆ 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 mp_wsm6_init(), and mp_wsm6_run().

◆ 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 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 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 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 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 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().

◆ 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 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 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 mp_wsm6_run(), refl10cm_wsm6(), slope_snow(), and slope_wsm6().

◆ pacrc

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

Referenced by mp_wsm6_init(), and mp_wsm6_run().

◆ pacrg

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

Referenced by mp_wsm6_init(), and mp_wsm6_run().

◆ pacrr

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

Referenced by 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 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 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 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(), and slope_wsm6().

◆ pidn0r

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

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

◆ 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(), and slope_wsm6().

◆ 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 mp_wsm6_init(), and mp_wsm6_run().

◆ precg2

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

Referenced by mp_wsm6_init(), and mp_wsm6_run().

◆ precr1

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

Referenced by mp_wsm6_init(), and mp_wsm6_run().

◆ precr2

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

Referenced by mp_wsm6_init(), and mp_wsm6_run().

◆ precs1

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

Referenced by mp_wsm6_init(), and mp_wsm6_run().

◆ precs2

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

Referenced by mp_wsm6_init(), and mp_wsm6_run().

◆ pvtg

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

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

◆ pvtr

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

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

◆ pvts

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

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

◆ 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 mp_wsm6_init(), and mp_wsm6_run().

◆ qck1

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

Referenced by 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 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 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 mp_wsm6_init().

◆ roqimax

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

Referenced by mp_wsm6_init(), and mp_wsm6_run().

◆ rslopeg2max

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

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

◆ rslopeg3max

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

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

◆ rslopegbmax

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

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

◆ rslopegmax

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

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

◆ rsloper2max

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

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

◆ rsloper3max

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

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

◆ rsloperbmax

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

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

◆ rslopermax

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

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

◆ rslopes2max

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

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

◆ rslopes3max

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

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

◆ rslopesbmax

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

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

◆ rslopesmax

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

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

◆ 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 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 mp_wsm6_init(), and mp_wsm6_run().