source: lmdz_wrf/trunk/WRFV3/phys/module_mp_wsm5_accel.F @ 1562

Last change on this file since 1562 was 1, checked in by lfita, 11 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 130.2 KB
Line 
1#if ( RWORDSIZE == 4 )
2#  define VREC vsrec
3#  define VSQRT vssqrt
4#else
5#  define VREC vrec
6#  define VSQRT vsqrt
7#endif
8
9!Including inline expansion statistical function
10MODULE module_mp_wsm5
11!
12!
13   REAL, PARAMETER, PRIVATE :: dtcldcr     = 120. ! maximum time step for minor loops
14   REAL, PARAMETER, PRIVATE :: n0r = 8.e6         ! intercept parameter rain
15   REAL, PARAMETER, PRIVATE :: avtr = 841.9       ! a constant for terminal velocity of rain
16   REAL, PARAMETER, PRIVATE :: bvtr = 0.8         ! a constant for terminal velocity of rain
17   REAL, PARAMETER, PRIVATE :: r0 = .8e-5         ! 8 microm  in contrast to 10 micro m
18   REAL, PARAMETER, PRIVATE :: peaut = .55        ! collection efficiency
19   REAL, PARAMETER, PRIVATE :: xncr = 3.e8        ! maritime cloud in contrast to 3.e8 in tc80
20   REAL, PARAMETER, PRIVATE :: xmyu = 1.718e-5    ! the dynamic viscosity kgm-1s-1
21   REAL, PARAMETER, PRIVATE :: avts = 11.72       ! a constant for terminal velocity of snow
22   REAL, PARAMETER, PRIVATE :: bvts = .41         ! a constant for terminal velocity of snow
23   REAL, PARAMETER, PRIVATE :: n0smax =  1.e11    ! maximum n0s (t=-90C unlimited)
24   REAL, PARAMETER, PRIVATE :: lamdarmax = 8.e4   ! limited maximum value for slope parameter of rain
25   REAL, PARAMETER, PRIVATE :: lamdasmax = 1.e5   ! limited maximum value for slope parameter of snow
26   REAL, PARAMETER, PRIVATE :: lamdagmax = 6.e4   ! limited maximum value for slope parameter of graupel
27   REAL, PARAMETER, PRIVATE :: dicon = 11.9       ! constant for the cloud-ice diamter
28   REAL, PARAMETER, PRIVATE :: dimax = 500.e-6    ! limited maximum value for the cloud-ice diamter
29   REAL, PARAMETER, PRIVATE :: n0s = 2.e6         ! temperature dependent intercept parameter snow
30   REAL, PARAMETER, PRIVATE :: alpha = .12        ! .122 exponen factor for n0s
31   REAL, PARAMETER, PRIVATE :: pfrz1 = 100.       ! constant in Biggs freezing
32   REAL, PARAMETER, PRIVATE :: pfrz2 = 0.66       ! constant in Biggs freezing
33   REAL, PARAMETER, PRIVATE :: qcrmin = 1.e-9     ! minimun values for qr, qs, and qg
34   REAL, PARAMETER, PRIVATE :: eacrc = 1.0        ! Snow/cloud-water collection efficiency
35   REAL, SAVE ::                                      &
36             qc0, qck1,bvtr1,bvtr2,bvtr3,bvtr4,g1pbr, &
37             g3pbr,g4pbr,g5pbro2,pvtr,eacrr,pacrr,    &
38             precr1,precr2,xmmax,roqimax,bvts1,       &
39             bvts2,bvts3,bvts4,g1pbs,g3pbs,g4pbs,     &
40             g5pbso2,pvts,pacrs,precs1,precs2,pidn0r, &
41             pidn0s,xlv1,pacrc,                       &
42             rslopermax,rslopesmax,rslopegmax,        &
43             rsloperbmax,rslopesbmax,rslopegbmax,     &
44             rsloper2max,rslopes2max,rslopeg2max,     &
45             rsloper3max,rslopes3max,rslopeg3max
46!
47! Specifies code-inlining of fpvs function in WSM52D below. JM 20040507
48!
49CONTAINS
50!===================================================================
51!
52  SUBROUTINE wsm5(th, q, qc, qr, qi, qs                            &
53                 ,den, pii, p, delz                                &
54                 ,delt,g, cpd, cpv, rd, rv, t0c                    &
55                 ,ep1, ep2, qmin                                   &
56                 ,XLS, XLV0, XLF0, den0, denr                      &
57                 ,cliq,cice,psat                                   &
58                 ,rain, rainncv                                    &
59                 ,snow, snowncv                                    &
60                 ,sr                                               &
61                 ,ids,ide, jds,jde, kds,kde                        &
62                 ,ims,ime, jms,jme, kms,kme                        &
63                 ,its,ite, jts,jte, kts,kte                        &
64                                                                   )
65#ifdef _OPENMP
66  use omp_lib
67#endif
68!-------------------------------------------------------------------
69  IMPLICIT NONE
70!-------------------------------------------------------------------
71!
72!  This code is a 5-class mixed ice microphyiscs scheme (WSM5) of the WRF
73!  Single-Moment MicroPhyiscs (WSMMP). The WSMMP assumes that ice nuclei
74!  number concentration is a function of temperature, and seperate assumption
75!  is developed, in which ice crystal number concentration is a function
76!  of ice amount. A theoretical background of the ice-microphysics and related
77!  processes in the WSMMPs are described in Hong et al. (2004).
78!  Production terms in the WSM6 scheme are described in Hong and Lim (2006).
79!  All units are in m.k.s. and source/sink terms in kgkg-1s-1.
80!
81!  WSM5 cloud scheme
82!
83!  Coded by Song-You Hong (Yonsei Univ.)
84!             Jimy Dudhia (NCAR) and Shu-Hua Chen (UC Davis)
85!             Summer 2002
86!
87!  Implemented by Song-You Hong (Yonsei Univ.) and Jimy Dudhia (NCAR)
88!             Summer 2003
89!
90!  Reference) Hong, Dudhia, Chen (HDC, 2004) Mon. Wea. Rev.
91!             Rutledge, Hobbs (RH83, 1983) J. Atmos. Sci.
92!             Hong and Lim (HL, 2006) J. Korean Meteor. Soc.
93!
94  INTEGER,      INTENT(IN   )    ::   ids,ide, jds,jde, kds,kde , &
95                                      ims,ime, jms,jme, kms,kme , &
96                                      its,ite, jts,jte, kts,kte
97  REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
98        INTENT(INOUT) ::                                          &
99                                                             th,  &
100                                                              q,  &
101                                                              qc, &
102                                                              qi, &
103                                                              qr, &
104                                                              qs
105  REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
106        INTENT(IN   ) ::                                          &
107                                                             den, &
108                                                             pii, &
109                                                               p, &
110                                                            delz
111  REAL, INTENT(IN   ) ::                                    delt, &
112                                                               g, &
113                                                              rd, &
114                                                              rv, &
115                                                             t0c, &
116                                                            den0, &
117                                                             cpd, &
118                                                             cpv, &
119                                                             ep1, &
120                                                             ep2, &
121                                                            qmin, &
122                                                             XLS, &
123                                                            XLV0, &
124                                                            XLF0, &
125                                                            cliq, &
126                                                            cice, &
127                                                            psat, &
128                                                            denr
129  REAL, DIMENSION( ims:ime , jms:jme ),                           &
130        INTENT(INOUT) ::                                    rain, &
131                                                         rainncv, &
132                                                              sr
133
134  REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL,                 &
135        INTENT(INOUT) ::                                    snow, &
136                                                         snowncv
137
138! LOCAL VAR
139  REAL, DIMENSION( its:ite , kts:kte ) ::   t
140  REAL, DIMENSION( its:ite , kts:kte, 2 ) ::   qci, qrs
141  CHARACTER*256 :: emess
142  INTEGER :: mkx_test
143  INTEGER ::               i,j,k
144
145#ifdef _ACCEL_PROF
146  INTEGER :: l
147  real*8 wsm3_t(8,256), wsm5_t(8,256), t1, t2
148  common /wsm_times/ wsm3_t(8,256), wsm5_t(8,256)
149#endif
150
151
152!-------------------------------------------------------------------
153
154#ifdef _ACCEL_PROF
155call cpu_time(t1)
156#endif
157
158#ifndef RUN_ON_GPU
159
160#ifdef _ACCEL
161
162      !  Need to send th, pii, qc, qi, qr, qs
163      !  Don't send t
164
165      CALL wsm52D(th, pii, q, qc, qr, qi, qs                       &
166                  ,den                                             &
167                  ,p, delz                                         &
168                  ,delt,g, cpd, cpv, rd, rv, t0c                   &
169                  ,ep1, ep2, qmin                                  &
170                  ,XLS, XLV0, XLF0, den0, denr                     &
171                  ,cliq,cice,psat                                  &
172                  ,rain,rainncv                                    &
173                  ,sr                                              &
174                  ,ids,ide, jds,jde, kds,kde                       &
175                  ,ims,ime, jms,jme, kms,kme                       &
176                  ,its,ite, jts,jte, kts,kte                       &
177                  ,snow,snowncv                                    &
178                                                                   )
179
180#else
181
182      DO j=jts,jte
183         DO k=kts,kte
184         DO i=its,ite
185            t(i,k)=th(i,k,j)*pii(i,k,j)
186            qci(i,k,1) = qc(i,k,j)
187            qci(i,k,2) = qi(i,k,j)
188            qrs(i,k,1) = qr(i,k,j)
189            qrs(i,k,2) = qs(i,k,j)
190         ENDDO
191         ENDDO
192
193         !  Sending array starting locations of optional variables may cause
194         !  troubles, so we explicitly change the call.
195
196         CALL wsm52D(t, q(ims,kms,j), qci, qrs                     &
197                    ,den(ims,kms,j)                                &
198                    ,p(ims,kms,j), delz(ims,kms,j)                 &
199                    ,delt,g, cpd, cpv, rd, rv, t0c                 &
200                    ,ep1, ep2, qmin                                &
201                    ,XLS, XLV0, XLF0, den0, denr                   &
202                    ,cliq,cice,psat                                &
203                    ,j                                             &
204                    ,rain(ims,j),rainncv(ims,j)                    &
205                    ,sr(ims,j)                                     &
206                    ,ids,ide, jds,jde, kds,kde                     &
207                    ,ims,ime, jms,jme, kms,kme                     &
208                    ,its,ite, jts,jte, kts,kte                     &
209                    ,snow,snowncv                                  &
210                                                                   )
211
212         DO K=kts,kte
213         DO I=its,ite
214            th(i,k,j)=t(i,k)/pii(i,k,j)
215            qc(i,k,j) = qci(i,k,1)
216            qi(i,k,j) = qci(i,k,2)
217            qr(i,k,j) = qrs(i,k,1)
218            qs(i,k,j) = qrs(i,k,2)
219         ENDDO
220         ENDDO
221      ENDDO
222#endif
223#else
224      CALL get_wsm5_gpu_levels ( mkx_test )
225      IF ( mkx_test .LT. kte ) THEN
226        WRITE(emess,*)'Number of levels compiled for GPU WSM5 too small. ',      &
227                      mkx_test,' < ',kte
228        CALL wrf_error_fatal(emess)
229      ENDIF
230      CALL wsm5_host (                                                           &
231                    th(its:ite,kts:kte,jts:jte), pii(its:ite,kts:kte,jts:jte)    &
232                   ,q(its:ite,kts:kte,jts:jte), qc(its:ite,kts:kte,jts:jte)      &
233                   ,qi(its:ite,kts:kte,jts:jte), qr(its:ite,kts:kte,jts:jte)     &
234                   ,qs(its:ite,kts:kte,jts:jte), den(its:ite,kts:kte,jts:jte)    &
235                   ,p(its:ite,kts:kte,jts:jte), delz(its:ite,kts:kte,jts:jte)    &
236                   ,delt                                                         &
237                   ,rain(its:ite,jts:jte),rainncv(its:ite,jts:jte)               &
238                   ,snow(its:ite,jts:jte),snowncv(its:ite,jts:jte)               &
239                   ,sr(its:ite,jts:jte)                                          &
240                   ,its, ite,  jts, jte,  kts, kte                               &
241                   ,its, ite,  jts, jte,  kts, kte                               &
242                   ,its, ite,  jts, jte,  kts, kte                               &
243          )
244#endif
245
246
247#ifdef _ACCEL_PROF
248  call cpu_time(t2)
249#ifdef _OPENMP
250  l = omp_get_thread_num() + 1
251#else
252  l = 1
253#endif
254  wsm5_t(1,l) = wsm5_t(1,l) + (t2 - t1)
255#endif
256
257  END SUBROUTINE wsm5
258
259
260#ifdef _ACCEL
261
262!===================================================================
263!
264  SUBROUTINE wsm52D(th, pii, q, qc, qr, qi, qqs, den, p, delz     &
265                   ,delt,g, cpd, cpv, rd, rv, t0c                 &
266                   ,ep1, ep2, qmin                                &
267                   ,XLS, XLV0, XLF0, den0, denr                   &
268                   ,cliq,cice,psat                                &
269                   ,rain,rainncv                                  &
270                   ,sr                                            &
271                   ,ids,ide, jds,jde, kds,kde                     &
272                   ,ims,ime, jms,jme, kms,kme                     &
273                   ,its,ite, jts,jte, kts,kte                     &
274                   ,snow,snowncv                                  &
275                                                                  )
276!-------------------------------------------------------------------
277  IMPLICIT NONE
278!-------------------------------------------------------------------
279  INTEGER,      INTENT(IN   )    ::   ids,ide, jds,jde, kds,kde , &
280                                      ims,ime, jms,jme, kms,kme , &
281                                      its,ite, jts,jte, kts,kte
282  REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
283        INTENT(INOUT) ::                                          &
284                                                              th
285  REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
286        INTENT(IN) ::                                             &
287                                                             pii
288  REAL, DIMENSION( ims:ime , kms:kme, jms:jme ),                  &
289        INTENT(INOUT) ::                                          &
290                                                              qc, &
291                                                              qr, &
292                                                              qi, &
293                                                             qqs
294  REAL, DIMENSION( ims:ime , kms:kme, jms:jme ),                  &
295        INTENT(INOUT) ::                                          &
296                                                               q
297  REAL, DIMENSION( ims:ime , kms:kme, jms:jme ),                  &
298        INTENT(IN   ) ::                                          &
299                                                             den, &
300                                                               p, &
301                                                            delz
302  REAL, INTENT(IN   ) ::                                    delt, &
303                                                               g, &
304                                                             cpd, &
305                                                             cpv, &
306                                                             t0c, &
307                                                            den0, &
308                                                              rd, &
309                                                              rv, &
310                                                             ep1, &
311                                                             ep2, &
312                                                            qmin, &
313                                                             XLS, &
314                                                            XLV0, &
315                                                            XLF0, &
316                                                            cliq, &
317                                                            cice, &
318                                                            psat, &
319                                                            denr
320  REAL, DIMENSION( ims:ime, jms:jme ),                            &
321        INTENT(INOUT) ::                                    rain, &
322                                                         rainncv, &
323                                                              sr
324
325  REAL, DIMENSION( ims:ime, jms:jme ),     OPTIONAL,              &
326        INTENT(INOUT) ::                                    snow, &
327                                                         snowncv
328
329! LOCAL VAR
330  REAL, DIMENSION( its:ite , kts:kte , 2) ::                      &
331                                                              rh, &
332                                                              qs, &
333                                                          rslope, &
334                                                         rslope2, &
335                                                         rslope3, &
336                                                         rslopeb, &
337                                                            falk, &
338                                                            fall, &
339                                                           work1
340  REAL, DIMENSION( its:ite , kts:kte, jts:jte ) ::                &
341                                                               t
342  REAL, DIMENSION( its:ite , kts:kte , 2 ) ::                     &
343                                                             qci, &
344                                                             qrs
345  REAL, DIMENSION( its:ite , kts:kte ) ::                         &
346                                                           falkc, &
347                                                           fallc, &
348                                                              xl, &
349                                                             cpm, &
350                                                          denfac, &
351                                                             xni, &
352                                                          n0sfac, &
353                                                           work2, &
354                                                          work1c, &
355                                                          work2c
356
357  REAL, DIMENSION( its:ite , kts:kte ) ::                         &
358                                                           pigen, &
359                                                           pidep, &
360                                                           psdep, &
361                                                           praut, &
362                                                           psaut, &
363                                                           prevp, &
364                                                           psevp, &
365                                                           pracw, &
366                                                           psacw, &
367                                                           psaci, &
368                                                           pcond, &
369                                                           psmlt
370  INTEGER ::                                                      &
371                                                           mstep, &
372                                                           numdt
373  REAL ::                                                 rmstep
374  REAL dtcldden, rdelz, rdtcld
375  LOGICAL ::                                              flgcld
376
377#define WSM_NO_CONDITIONAL_IN_VECTOR
378#ifdef WSM_NO_CONDITIONAL_IN_VECTOR
379  REAL ::                     xal, xbl
380#endif
381
382  REAL  ::  pi,                                                   &
383            cpmcal, xlcal, lamdar, lamdas, diffus,                &
384            viscos, xka, venfac, conden, diffac,                  &
385            x, y, z, a, b, c, d, e,                               &
386            qdt, holdrr, holdrs, supcol, supcolt, pvt,            &
387            coeres, supsat, dtcld, xmi, eacrs, satdt,             &
388            vt2i,vt2s,acrfac,                                     &
389            qimax, diameter, xni0, roqi0,                         &
390            fallsum, fallsum_qsi, xlwork2, factor, source,        &
391            value, xlf, pfrzdtc, pfrzdtr, supice,  holdc, holdci
392! variables for optimization
393  REAL, DIMENSION( its:ite )           ::                    tvec1
394  REAL ::                                                    temp
395  INTEGER :: i, j, k,                                             &
396            iprt, latd, lond, loop, loops, ifsat, n
397! Temporaries used for inlining fpvs function
398  REAL  :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
399  REAL  :: logtr
400!
401!=================================================================
402!   compute internal functions
403!
404      cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
405      xlcal(x) = xlv0-xlv1*(x-t0c)
406!----------------------------------------------------------------
407!     size distributions: (x=mixing ratio, y=air density):
408!     valid for mixing ratio > 1.e-9 kg/kg.
409!
410! Optimizatin : A**B => exp(log(A)*(B))
411      lamdar(x,y)=   sqrt(sqrt(pidn0r/(x*y)))      ! (pidn0r/(x*y))**.25
412      lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y)))    ! (pidn0s*z/(x*y))**.25
413!
414!----------------------------------------------------------------
415!     diffus: diffusion coefficient of the water vapor
416!     viscos: kinematic viscosity(m2s-1)
417!     diffus(x,y) = 8.794e-5 * exp(log(x)*(1.81)) / y        ! 8.794e-5*x**1.81/y
418!     viscos(x,y) = 1.496e-6 * (x*sqrt(x)) /(x+120.)/y  ! 1.496e-6*x**1.5/(x+120.)/y
419!     xka(x,y) = 1.414e3*viscos(x,y)*y
420!     diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b))
421!     venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333)))         &
422!                    /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
423!     conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+d*d/(rv*e)*c/(a*a))
424!
425!
426      pi = 4. * atan(1.)
427!
428!----------------------------------------------------------------
429!     paddint 0 for negative values generated by dynamics
430!
431
432!
433! Moved outside of accelerator region
434!
435      loops = max(nint(delt/dtcldcr),1)
436      dtcld = delt/loops
437      if(delt.le.dtcldcr) dtcld = delt
438!
439
440!....!$acc        local(t) &
441
442      IF ( PRESENT (snowncv) .AND. PRESENT (snow)) THEN
443
444!$acc region &
445!$acc        local(t) &
446!$acc        copyin(delz(:,:,:),p(:,:,:),den(:,:,:),pii(:,:,:)) &
447!$acc        copyout(snowncv(:,:),rainncv(:,:),sr(:,:)) &
448!$acc        copy(qqs(:,:,:),qr(:,:,:),qi(:,:,:),qc(:,:,:)) &
449!$acc        copy(th(:,:,:),q(:,:,:),snow(:,:),rain(:,:))
450!$acc do &
451!$acc        private(rh,qs,rslope,rslope2,rslope3,rslopeb,falk,fall) &
452!$acc        private(work1,qci,qrs,falkc,fallc,xl,cpm,denfac,xni) &
453!$acc        private(n0sfac,work2,work1c,work2c,pigen,pidep,psdep) &
454!$acc        private(praut,psaut,prevp,psevp) &
455!$acc        private(pracw,psacw,psaci,pcond,psmlt) &
456!$acc        parallel
457      do j = jts, jte
458!$acc do &
459!$acc        private(numdt,mstep) &
460!$acc        kernel vector
461        do i = its, ite
462        do k = kts, kte
463          t(i,k,j)=th(i,k,j)*pii(i,k,j)
464          qci(i,k,1) = max(qc(i,k,j),0.0)
465          qci(i,k,2) = max(qi(i,k,j),0.0)
466          qrs(i,k,1) = max(qr(i,k,j),0.0)
467          qrs(i,k,2) = max(qqs(i,k,j),0.0)
468        enddo
469!
470!----------------------------------------------------------------
471!     latent heat for phase changes and heat capacity. neglect the
472!     changes during microphysical process calculation
473!     emanuel(1994)
474!
475      do k = kts, kte
476          cpm(i,k) = cpmcal(q(i,k,j))
477          xl(i,k) = xlcal(t(i,k,j))
478      enddo
479!
480!----------------------------------------------------------------
481!     compute the minor time steps.
482!
483!     loops = max(nint(delt/dtcldcr),1)
484!     dtcld = delt/loops
485!     if(delt.le.dtcldcr) dtcld = delt
486!
487      do loop = 1,loops
488!
489!----------------------------------------------------------------
490!     initialize the large scale variables
491!
492        mstep = 1
493        flgcld = .true.
494!
495      do k = kts, kte
496          denfac(i,k) = sqrt(den0/den(i,k,j))
497      enddo
498!     do k = kts, kte
499!       CALL VREC( tvec1(its), den(its,k,j), ite-its+1)
500!       do i = its, ite
501!         tvec1(i) = tvec1(i)*den0
502!       enddo
503!       CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
504!     enddo
505!
506! Inline expansion for fpvs
507!         qs(i,k,1) = fpvs(t(i,k,j),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
508!         qs(i,k,2) = fpvs(t(i,k,j),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
509      hsub = xls
510      hvap = xlv0
511      cvap = cpv
512      ttp=t0c+0.01
513      dldt=cvap-cliq
514      xa=-dldt/rv
515      xb=xa+hvap/(rv*ttp)
516      dldti=cvap-cice
517      xai=-dldti/rv
518      xbi=xai+hsub/(rv*ttp)
519
520! this is for compilers where the conditional inhibits vectorization
521#ifdef WSM_NO_CONDITIONAL_IN_VECTOR
522      do k = kts, kte
523          if(t(i,k,j).lt.ttp) then
524            xal = xai
525            xbl = xbi
526          else
527            xal = xa
528            xbl = xb
529          endif
530          tr=ttp/t(i,k,j)
531          logtr=log(tr)
532          qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
533          qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k,j) - qs(i,k,1))
534          qs(i,k,1) = max(qs(i,k,1),qmin)
535          rh(i,k,1) = max(q(i,k,j) / qs(i,k,1),qmin)
536          qs(i,k,2)=psat*exp(logtr*(xal)+xbl*(1.-tr))
537          qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k,j) - qs(i,k,2))
538          qs(i,k,2) = max(qs(i,k,2),qmin)
539          rh(i,k,2) = max(q(i,k,j) / qs(i,k,2),qmin)
540      enddo
541#else
542      do k = kts, kte
543          tr=ttp/t(i,k,j)
544          logtr=log(tr)
545          qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
546          qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k,j) - qs(i,k,1))
547          qs(i,k,1) = max(qs(i,k,1),qmin)
548          rh(i,k,1) = max(q(i,k,j) / qs(i,k,1),qmin)
549          if(t(i,k,j).lt.ttp) then
550            qs(i,k,2)=psat*exp(logtr*(xai)+xbi*(1.-tr))
551          else
552            qs(i,k,2)=psat*exp(logtr*(xa)+xb*(1.-tr))
553          endif
554          qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k,j) - qs(i,k,2))
555          qs(i,k,2) = max(qs(i,k,2),qmin)
556          rh(i,k,2) = max(q(i,k,j) / qs(i,k,2),qmin)
557      enddo
558#endif
559
560!
561!----------------------------------------------------------------
562!     initialize the variables for microphysical physics
563!
564!
565      do k = kts, kte
566          prevp(i,k) = 0.
567          psdep(i,k) = 0.
568          praut(i,k) = 0.
569          psaut(i,k) = 0.
570          pracw(i,k) = 0.
571          psaci(i,k) = 0.
572          psacw(i,k) = 0.
573          pigen(i,k) = 0.
574          pidep(i,k) = 0.
575          pcond(i,k) = 0.
576          psmlt(i,k) = 0.
577          psevp(i,k) = 0.
578          falk(i,k,1) = 0.
579          falk(i,k,2) = 0.
580          fall(i,k,1) = 0.
581          fall(i,k,2) = 0.
582          fallc(i,k) = 0.
583          falkc(i,k) = 0.
584          xni(i,k) = 1.e3
585      enddo
586!
587!----------------------------------------------------------------
588!     compute the fallout term:
589!     first, vertical terminal velosity for minor loops
590!
591      do k = kts, kte
592          supcol = t0c-t(i,k,j)
593!---------------------------------------------------------------
594! n0s: Intercept parameter for snow [m-4] [HDC 6]
595!---------------------------------------------------------------
596          n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
597          if(qrs(i,k,1).le.qcrmin)then
598            rslope(i,k,1) = rslopermax
599            rslopeb(i,k,1) = rsloperbmax
600            rslope2(i,k,1) = rsloper2max
601            rslope3(i,k,1) = rsloper3max
602          else
603            rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k,j))
604            rslopeb(i,k,1) = exp(log(rslope(i,k,1))*(bvtr))
605            rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
606            rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
607          endif
608          if(qrs(i,k,2).le.qcrmin)then
609            rslope(i,k,2) = rslopesmax
610            rslopeb(i,k,2) = rslopesbmax
611            rslope2(i,k,2) = rslopes2max
612            rslope3(i,k,2) = rslopes3max
613          else
614            rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k,j),n0sfac(i,k))
615            rslopeb(i,k,2) = exp(log(rslope(i,k,2))*(bvts))
616            rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
617            rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
618          endif
619!-------------------------------------------------------------
620! Ni: ice crystal number concentraiton   [HDC 5c]
621!-------------------------------------------------------------
622!         xni(i,k) = min(max(5.38e7*(den(i,k,j)                                  &
623!                   *max(qci(i,k,2),qmin))**0.75,1.e3),1.e6)
624          temp = (den(i,k,j)*max(qci(i,k,2),qmin))
625          temp = sqrt(sqrt(temp*temp*temp))
626          xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
627      enddo
628!
629      numdt = 1
630      do k = kte, kts, -1
631          work1(i,k,1) = pvtr*rslopeb(i,k,1)*denfac(i,k)/delz(i,k,j)
632          work1(i,k,2) = pvts*rslopeb(i,k,2)*denfac(i,k)/delz(i,k,j)
633          numdt = max(nint(max(work1(i,k,1),work1(i,k,2))*dtcld+.5),1)
634          if(numdt.ge.mstep) mstep = numdt
635      enddo
636        rmstep = 1./mstep
637!
638      do n = 1, mstep
639        k = kte
640!             falk(i,k,1) = den(i,k,j)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
641!             falk(i,k,2) = den(i,k,j)*qrs(i,k,2)*work1(i,k,2)/mstep(i)
642              falk(i,k,1) = den(i,k,j)*qrs(i,k,1)*work1(i,k,1)*rmstep
643              falk(i,k,2) = den(i,k,j)*qrs(i,k,2)*work1(i,k,2)*rmstep
644              fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
645              fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
646!             qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcld/den(i,k,j),0.)
647!             qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcld/den(i,k,j),0.)
648              dtcldden = dtcld/den(i,k,j)
649              qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcldden,0.)
650              qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcldden,0.)
651!           endif
652        do k = kte-1, kts, -1
653              falk(i,k,1) = den(i,k,j)*qrs(i,k,1)*work1(i,k,1)*rmstep
654              falk(i,k,2) = den(i,k,j)*qrs(i,k,2)*work1(i,k,2)*rmstep
655              fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
656              fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
657              dtcldden = dtcld/den(i,k,j)
658              rdelz = 1./delz(i,k,j)
659              qrs(i,k,1) = max(qrs(i,k,1)-(falk(i,k,1)-falk(i,k+1,1)           &
660                          *delz(i,k+1,j)*rdelz)*dtcldden,0.)
661              qrs(i,k,2) = max(qrs(i,k,2)-(falk(i,k,2)-falk(i,k+1,2)           &
662                          *delz(i,k+1,j)*rdelz)*dtcldden,0.)
663        enddo
664        do k = kte, kts, -1
665              if(t(i,k,j).gt.t0c.and.qrs(i,k,2).gt.0.) then
666!----------------------------------------------------------------
667! psmlt: melting of snow [HL A33] [RH83 A25]
668!       (T>T0: S->R)
669!----------------------------------------------------------------
670                xlf = xlf0
671!               work2(i,k)= venfac(p(i,k),t(i,k,j),den(i,k,j))
672                work2(i,k)= (exp(log(((1.496e-6*((t(i,k,j))*sqrt(t(i,k,j)))        &
673                            /((t(i,k,j))+120.)/(den(i,k,j)))/(8.794e-5             &
674                            *exp(log(t(i,k,j))*(1.81))/p(i,k,j))))                 &
675                            *((.3333333)))/sqrt((1.496e-6*((t(i,k,j))            &
676                            *sqrt(t(i,k,j)))/((t(i,k,j))+120.)/(den(i,k,j))))        &
677                            *sqrt(sqrt(den0/(den(i,k,j)))))
678                coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
679!               psmlt(i,k) = xka(t(i,k,j),den(i,k,j))/xlf*(t0c-t(i,k,j))*pi/2.       &
680!                           *n0sfac(i,k)*(precs1*rslope2(i,k,2)+precs2         &
681!                           *work2(i,k)*coeres)
682                psmlt(i,k) = (1.414e3*(1.496e-6*((t(i,k,j))*sqrt(t(i,k,j)))        &         
683                            /((t(i,k,j))+120.)/(den(i,k,j)) )*(den(i,k,j)))          &
684                            /xlf*(t0c-t(i,k,j))*pi/2.                            &
685                            *n0sfac(i,k)*(precs1*rslope2(i,k,2)+precs2         &
686                            *work2(i,k)*coeres)
687                psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep,                &
688                            -qrs(i,k,2)/mstep),0.)
689                qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
690                qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k)
691                t(i,k,j) = t(i,k,j) + xlf/cpm(i,k)*psmlt(i,k)
692              endif
693        enddo
694      enddo
695
696
697!---------------------------------------------------------------
698! Vice [ms-1] : fallout of ice crystal [HDC 5a]
699!---------------------------------------------------------------
700      mstep = 1
701      numdt = 1
702      do k = kte, kts, -1
703          if(qci(i,k,2).le.0.) then
704            work2c(i,k) = 0.
705          else
706            xmi = den(i,k,j)*qci(i,k,2)/xni(i,k)
707!           diameter  = min(dicon * sqrt(xmi),dimax)
708            diameter  = max(min(dicon * sqrt(xmi),dimax), 1.e-25)
709            work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31))
710            work2c(i,k) = work1c(i,k)/delz(i,k,j)
711          endif
712          numdt = max(nint(work2c(i,k)*dtcld+.5),1)
713          if(numdt.ge.mstep) mstep = numdt
714      enddo
715!
716      do n = 1, mstep
717        k = kte
718            falkc(i,k) = den(i,k,j)*qci(i,k,2)*work2c(i,k)/mstep
719            holdc = falkc(i,k)
720            fallc(i,k) = fallc(i,k)+falkc(i,k)
721            holdci = qci(i,k,2)
722            qci(i,k,2) = max(qci(i,k,2)-falkc(i,k)*dtcld/den(i,k,j),0.)
723!         endif
724        do k = kte-1, kts, -1
725              falkc(i,k) = den(i,k,j)*qci(i,k,2)*work2c(i,k)/mstep
726              holdc = falkc(i,k)
727              fallc(i,k) = fallc(i,k)+falkc(i,k)
728              holdci = qci(i,k,2)
729              qci(i,k,2) = max(qci(i,k,2)-(falkc(i,k)-falkc(i,k+1)             &
730                          *delz(i,k+1,j)/delz(i,k,j))*dtcld/den(i,k,j),0.)
731!           endif
732        enddo
733      enddo
734!
735!
736!----------------------------------------------------------------
737!      rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
738!
739        fallsum = fall(i,1,1)+fall(i,1,2)+fallc(i,1)
740        fallsum_qsi = fall(i,1,2)+fallc(i,1)
741        rainncv(i,j) = 0.
742        if(fallsum.gt.0.) then
743          rainncv(i,j) = fallsum*delz(i,1,j)/denr*dtcld*1000.
744          rain(i,j) = fallsum*delz(i,1,j)/denr*dtcld*1000. + rain(i,j)
745        endif
746        snowncv(i,j) = 0.
747        if(fallsum_qsi.gt.0.) then
748          snowncv(i,j) = fallsum_qsi*delz(i,kts,j)/denr*dtcld*1000.
749          snow(i,j) = fallsum_qsi*delz(i,kts,j)/denr*dtcld*1000. + snow(i,j)
750        endif
751        sr(i,j) = 0.
752        if(fallsum.gt.0.)sr(i,j)=fallsum_qsi*delz(i,kts,j)/denr*dtcld*1000.        &   
753        /(rainncv(i,j)+1.e-12)
754!
755!---------------------------------------------------------------
756! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
757!       (T>T0: I->C)
758!---------------------------------------------------------------
759      do k = kts, kte
760          supcol = t0c-t(i,k,j)
761          xlf = xls-xl(i,k)
762          if(supcol.lt.0.) xlf = xlf0
763          if(supcol.lt.0.and.qci(i,k,2).gt.0.) then
764            qci(i,k,1) = qci(i,k,1) + qci(i,k,2)
765            t(i,k,j) = t(i,k,j) - xlf/cpm(i,k)*qci(i,k,2)
766            qci(i,k,2) = 0.
767          endif
768!---------------------------------------------------------------
769! pihmf: homogeneous freezing of cloud water below -40c [HL A45]
770!        (T<-40C: C->I)
771!---------------------------------------------------------------
772          if(supcol.gt.40..and.qci(i,k,1).gt.0.) then
773            qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
774            t(i,k,j) = t(i,k,j) + xlf/cpm(i,k)*qci(i,k,1)
775            qci(i,k,1) = 0.
776          endif
777!---------------------------------------------------------------
778! pihtf: heterogeneous freezing of cloud water [HL A44]
779!        (T0>T>-40C: C->I)
780!---------------------------------------------------------------
781          if(supcol.gt.0..and.qci(i,k,1).gt.0.) then
782            supcolt=min(supcol,50.)
783!           pfrzdtc = min(pfrz1*(exp(pfrz2*supcol)-1.)                         &
784!              *den(i,k,j)/denr/xncr*qci(i,k,1)**2*dtcld,qci(i,k,1))
785            pfrzdtc = min(pfrz1*(exp(pfrz2*supcolt)-1.)                        &
786            *den(i,k,j)/denr/xncr*qci(i,k,1)*qci(i,k,1)*dtcld,qci(i,k,1))
787            qci(i,k,2) = qci(i,k,2) + pfrzdtc
788            t(i,k,j) = t(i,k,j) + xlf/cpm(i,k)*pfrzdtc
789            qci(i,k,1) = qci(i,k,1)-pfrzdtc
790          endif
791!---------------------------------------------------------------
792! psfrz: freezing of rain water [HL A20] [LFO 45]
793!        (T<T0, R->S)
794!---------------------------------------------------------------
795          if(supcol.gt.0..and.qrs(i,k,1).gt.0.) then
796            supcolt=min(supcol,50.)
797!           pfrzdtr = min(20.*pi**2*pfrz1*n0r*denr/den(i,k,j)                    &
798!                 *(exp(pfrz2*supcol)-1.)*rslope(i,k,1)**7*dtcld,              &
799!                 qrs(i,k,1))
800            temp = rslope(i,k,1)
801            temp = temp*temp*temp*temp*temp*temp*temp
802            pfrzdtr = min(20.*(pi*pi)*pfrz1*n0r*denr/den(i,k,j)                  &
803                  *(exp(pfrz2*supcolt)-1.)*temp*dtcld,                         &
804                  qrs(i,k,1))
805            qrs(i,k,2) = qrs(i,k,2) + pfrzdtr
806            t(i,k,j) = t(i,k,j) + xlf/cpm(i,k)*pfrzdtr
807            qrs(i,k,1) = qrs(i,k,1)-pfrzdtr
808          endif
809      enddo
810!
811!----------------------------------------------------------------
812!     rsloper: reverse of the slope parameter of the rain(m)
813!     xka:    thermal conductivity of air(jm-1s-1k-1)
814!     work1:  the thermodynamic term in the denominator associated with
815!             heat conduction and vapor diffusion
816!             (ry88, y93, h85)
817!     work2: parameter associated with the ventilation effects(y93)
818!
819      do k = kts, kte
820          if(qrs(i,k,1).le.qcrmin)then
821            rslope(i,k,1) = rslopermax
822            rslopeb(i,k,1) = rsloperbmax
823            rslope2(i,k,1) = rsloper2max
824            rslope3(i,k,1) = rsloper3max
825          else
826!           rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k,j))
827            rslope(i,k,1) = 1./(sqrt(sqrt(pidn0r/((qrs(i,k,1))*(den(i,k,j))))))
828            rslopeb(i,k,1) = exp(log(rslope(i,k,1))*(bvtr))
829            rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
830            rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
831          endif
832          if(qrs(i,k,2).le.qcrmin)then
833            rslope(i,k,2) = rslopesmax
834            rslopeb(i,k,2) = rslopesbmax
835            rslope2(i,k,2) = rslopes2max
836            rslope3(i,k,2) = rslopes3max
837          else
838!            rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k,j),n0sfac(i,k))
839            rslope(i,k,2) = 1./(sqrt(sqrt(pidn0s*(n0sfac(i,k))/((qrs(i,k,2))   & 
840                          *(den(i,k,j))))))
841            rslopeb(i,k,2) = exp(log(rslope(i,k,2))*(bvts))
842            rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
843            rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
844          endif
845      enddo
846!
847      do k = kts, kte
848!         work1(i,k,1) = diffac(xl(i,k),p(i,k,j),t(i,k,j),den(i,k,j),qs(i,k,1))
849          work1(i,k,1) = ((((den(i,k,j))*(xl(i,k))*(xl(i,k)))*((t(i,k,j))+120.)    &
850                       *(den(i,k,j)))/(1.414e3*(1.496e-6*((t(i,k,j))*sqrt(t(i,k,j))))&
851                       *(den(i,k,j))*(rv*(t(i,k,j))*(t(i,k,j)))))                    &
852                      +  p(i,k,j)/((qs(i,k,1))*(8.794e-5*exp(log(t(i,k,j))*(1.81))))
853!         work1(i,k,2) = diffac(xls,p(i,k,j),t(i,k,j),den(i,k,j),qs(i,k,2))
854          work1(i,k,2) = ((((den(i,k,j))*(xls)*(xls))*((t(i,k,j))+120.)*(den(i,k,j)))&
855                       /(1.414e3*(1.496e-6*((t(i,k,j))*sqrt(t(i,k,j))))*(den(i,k,j)) &
856                       *(rv*(t(i,k,j))*(t(i,k,j))))                                &
857                      + p(i,k,j)/(qs(i,k,2)*(8.794e-5*exp(log(t(i,k,j))*(1.81)))))
858!         work2(i,k) = venfac(p(i,k,j),t(i,k,j),den(i,k,j))
859          work2(i,k) = (exp(.3333333*log(((1.496e-6 * ((t(i,k,j))*sqrt(t(i,k,j)))) &
860                     *p(i,k,j))/(((t(i,k,j))+120.)*den(i,k,j)*(8.794e-5              &
861                     *exp(log(t(i,k,j))*(1.81))))))*sqrt(sqrt(den0/(den(i,k,j))))) &
862                      /sqrt((1.496e-6*((t(i,k,j))*sqrt(t(i,k,j))))                 &
863                      /(((t(i,k,j))+120.)*den(i,k,j)))
864      enddo
865!
866!===============================================================
867!
868! warm rain processes
869!
870! - follows the processes in RH83 and LFO except for autoconcersion
871!
872!===============================================================
873!
874      do k = kts, kte
875          supsat = max(q(i,k,j),qmin)-qs(i,k,1)
876          satdt = supsat/dtcld
877!---------------------------------------------------------------
878! praut: auto conversion rate from cloud to rain [HDC 16]
879!        (C->R)
880!---------------------------------------------------------------
881          if(qci(i,k,1).gt.qc0) then
882            praut(i,k) = qck1*exp(log(qci(i,k,1))*((7./3.)))
883            praut(i,k) = min(praut(i,k),qci(i,k,1)/dtcld)
884          endif
885!---------------------------------------------------------------
886! pracw: accretion of cloud water by rain [HL A40] [LFO 51]
887!        (C->R)
888!---------------------------------------------------------------
889          if(qrs(i,k,1).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
890            pracw(i,k) = min(pacrr*rslope3(i,k,1)*rslopeb(i,k,1)               &
891                         *qci(i,k,1)*denfac(i,k),qci(i,k,1)/dtcld)
892          endif
893!---------------------------------------------------------------
894! prevp: evaporation/condensation rate of rain [HDC 14]
895!        (V->R or R->V)
896!---------------------------------------------------------------
897          if(qrs(i,k,1).gt.0.) then
898            coeres = rslope2(i,k,1)*sqrt(rslope(i,k,1)*rslopeb(i,k,1))
899            prevp(i,k) = (rh(i,k,1)-1.)*(precr1*rslope2(i,k,1)                 &
900                         +precr2*work2(i,k)*coeres)/work1(i,k,1)
901            if(prevp(i,k).lt.0.) then
902              prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
903              prevp(i,k) = max(prevp(i,k),satdt/2)
904            else
905              prevp(i,k) = min(prevp(i,k),satdt/2)
906            endif
907          endif
908      enddo
909!
910!===============================================================
911!
912! cold rain processes
913!
914! - follows the revised ice microphysics processes in HDC
915! - the processes same as in RH83 and RH84  and LFO behave
916!   following ice crystal hapits defined in HDC, inclduing
917!   intercept parameter for snow (n0s), ice crystal number
918!   concentration (ni), ice nuclei number concentration
919!   (n0i), ice diameter (d)
920!
921!===============================================================
922!
923      rdtcld = 1./dtcld
924      do k = kts, kte
925          supcol = t0c-t(i,k,j)
926          supsat = max(q(i,k,j),qmin)-qs(i,k,2)
927          satdt = supsat/dtcld
928          ifsat = 0
929!-------------------------------------------------------------
930! Ni: ice crystal number concentraiton   [HDC 5c]
931!-------------------------------------------------------------
932!         xni(i,k) = min(max(5.38e7*(den(i,k,j)                                  &
933!                      *max(qci(i,k,2),qmin))**0.75,1.e3),1.e6)
934          temp = (den(i,k,j)*max(qci(i,k,2),qmin))
935          temp = sqrt(sqrt(temp*temp*temp))
936          xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
937          eacrs = exp(0.07*(-supcol))
938!
939          if(supcol.gt.0) then
940            if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,2).gt.qmin) then
941              xmi = den(i,k,j)*qci(i,k,2)/xni(i,k)
942              diameter  = min(dicon * sqrt(xmi),dimax)
943              vt2i = 1.49e4*diameter**1.31
944              vt2s = pvts*rslopeb(i,k,2)*denfac(i,k)
945!-------------------------------------------------------------
946! psaci: Accretion of cloud ice by rain [HDC 10]
947!        (T<T0: I->S)
948!-------------------------------------------------------------
949              acrfac = 2.*rslope3(i,k,2)+2.*diameter*rslope2(i,k,2)            &
950                      +diameter**2*rslope(i,k,2)
951              psaci(i,k) = pi*qci(i,k,2)*eacrs*n0s*n0sfac(i,k)                 &
952                           *abs(vt2s-vt2i)*acrfac/4.
953            endif
954          endif
955!-------------------------------------------------------------
956! psacw: Accretion of cloud water by snow  [HL A7] [LFO 24]
957!        (T<T0: C->S, and T>=T0: C->R)
958!-------------------------------------------------------------
959          if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
960            psacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2)                  &
961                           *rslopeb(i,k,2)*qci(i,k,1)*denfac(i,k)              &
962!                          ,qci(i,k,1)/dtcld)
963                           ,qci(i,k,1)*rdtcld)
964          endif
965          if(supcol .gt. 0) then
966!-------------------------------------------------------------
967! pidep: Deposition/Sublimation rate of ice [HDC 9]
968!       (T<T0: V->I or I->V)
969!-------------------------------------------------------------
970            if(qci(i,k,2).gt.0.and.ifsat.ne.1) then
971              xmi = den(i,k,j)*qci(i,k,2)/xni(i,k)
972              diameter = dicon * sqrt(xmi)
973              pidep(i,k) = 4.*diameter*xni(i,k)*(rh(i,k,2)-1.)/work1(i,k,2)
974              supice = satdt-prevp(i,k)
975              if(pidep(i,k).lt.0.) then
976!               pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
977!               pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
978                pidep(i,k) = max(max(pidep(i,k),satdt*.5),supice)
979                pidep(i,k) = max(pidep(i,k),-qci(i,k,2)*rdtcld)
980              else
981!               pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
982                pidep(i,k) = min(min(pidep(i,k),satdt*.5),supice)
983              endif
984              if(abs(prevp(i,k)+pidep(i,k)).ge.abs(satdt)) ifsat = 1
985            endif
986!-------------------------------------------------------------
987! psdep: deposition/sublimation rate of snow [HDC 14]
988!        (V->S or S->V)
989!-------------------------------------------------------------
990            if(qrs(i,k,2).gt.0..and.ifsat.ne.1) then
991              coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
992              psdep(i,k) = (rh(i,k,2)-1.)*n0sfac(i,k)                          &
993                           *(precs1*rslope2(i,k,2)+precs2                      &
994                           *work2(i,k)*coeres)/work1(i,k,2)
995              supice = satdt-prevp(i,k)-pidep(i,k)
996              if(psdep(i,k).lt.0.) then
997!               psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
998!               psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
999                psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)*rdtcld)
1000                psdep(i,k) = max(max(psdep(i,k),satdt*.5),supice)
1001              else
1002!             psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
1003                psdep(i,k) = min(min(psdep(i,k),satdt*.5),supice)
1004              endif
1005              if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt))          &
1006                ifsat = 1
1007            endif
1008!-------------------------------------------------------------
1009! pigen: generation(nucleation) of ice from vapor [HL A50] [HDC 7-8]
1010!       (T<T0: V->I)
1011!-------------------------------------------------------------
1012            if(supsat.gt.0.and.ifsat.ne.1) then
1013              supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
1014              xni0 = 1.e3*exp(0.1*supcol)
1015              roqi0 = 4.92e-11*exp(log(xni0)*(1.33))
1016              pigen(i,k) = max(0.,(roqi0/den(i,k,j)-max(qci(i,k,2),0.))          &
1017!                        /dtcld)
1018                         *rdtcld)
1019              pigen(i,k) = min(min(pigen(i,k),satdt),supice)
1020            endif
1021!
1022!-------------------------------------------------------------
1023! psaut: conversion(aggregation) of ice to snow [HDC 12]
1024!       (T<T0: I->S)
1025!-------------------------------------------------------------
1026            if(qci(i,k,2).gt.0.) then
1027              qimax = roqimax/den(i,k,j)
1028!             psaut(i,k) = max(0.,(qci(i,k,2)-qimax)/dtcld)
1029              psaut(i,k) = max(0.,(qci(i,k,2)-qimax)*rdtcld)
1030            endif
1031          endif
1032!-------------------------------------------------------------
1033! psevp: Evaporation of melting snow [HL A35] [RH83 A27]
1034!       (T>T0: S->V)
1035!-------------------------------------------------------------
1036          if(supcol.lt.0.) then
1037            if(qrs(i,k,2).gt.0..and.rh(i,k,1).lt.1.)                           &
1038              psevp(i,k) = psdep(i,k)*work1(i,k,2)/work1(i,k,1)
1039!              psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)/dtcld),0.)
1040              psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)*rdtcld),0.)
1041          endif
1042      enddo
1043!
1044!
1045!----------------------------------------------------------------
1046!     check mass conservation of generation terms and feedback to the
1047!     large scale
1048!
1049      do k = kts, kte
1050          if(t(i,k,j).le.t0c) then
1051!
1052!     cloud water
1053!
1054            value = max(qmin,qci(i,k,1))
1055            source = (praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
1056            if (source.gt.value) then
1057              factor = value/source
1058              praut(i,k) = praut(i,k)*factor
1059              pracw(i,k) = pracw(i,k)*factor
1060              psacw(i,k) = psacw(i,k)*factor
1061            endif
1062!
1063!     cloud ice
1064!
1065            value = max(qmin,qci(i,k,2))
1066            source = (psaut(i,k)+psaci(i,k)-pigen(i,k)-pidep(i,k))*dtcld
1067            if (source.gt.value) then
1068              factor = value/source
1069              psaut(i,k) = psaut(i,k)*factor
1070              psaci(i,k) = psaci(i,k)*factor
1071              pigen(i,k) = pigen(i,k)*factor
1072              pidep(i,k) = pidep(i,k)*factor
1073            endif
1074!
1075!     rain
1076!
1077!
1078            value = max(qmin,qrs(i,k,1))
1079            source = (-praut(i,k)-pracw(i,k)-prevp(i,k))*dtcld
1080            if (source.gt.value) then
1081              factor = value/source
1082              praut(i,k) = praut(i,k)*factor
1083              pracw(i,k) = pracw(i,k)*factor
1084              prevp(i,k) = prevp(i,k)*factor
1085            endif
1086!
1087!    snow
1088!
1089            value = max(qmin,qrs(i,k,2))
1090            source = (-psdep(i,k)-psaut(i,k)-psaci(i,k)-psacw(i,k))*dtcld 
1091            if (source.gt.value) then
1092              factor = value/source
1093              psdep(i,k) = psdep(i,k)*factor
1094              psaut(i,k) = psaut(i,k)*factor
1095              psaci(i,k) = psaci(i,k)*factor
1096              psacw(i,k) = psacw(i,k)*factor
1097            endif
1098!
1099            work2(i,k)=-(prevp(i,k)+psdep(i,k)+pigen(i,k)+pidep(i,k))
1100!     update
1101            q(i,k,j) = q(i,k,j)+work2(i,k)*dtcld
1102            qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k)                 &
1103                        +psacw(i,k))*dtcld,0.)
1104            qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k)                 &
1105                        +prevp(i,k))*dtcld,0.)
1106            qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+psaci(i,k)                 &
1107                        -pigen(i,k)-pidep(i,k))*dtcld,0.)
1108            qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k)                 &
1109                        +psaci(i,k)+psacw(i,k))*dtcld,0.)
1110            xlf = xls-xl(i,k)
1111            xlwork2 = -xls*(psdep(i,k)+pidep(i,k)+pigen(i,k))                  &
1112                      -xl(i,k)*prevp(i,k)-xlf*psacw(i,k)
1113            t(i,k,j) = t(i,k,j)-xlwork2/cpm(i,k)*dtcld
1114          else
1115!
1116!     cloud water
1117!
1118            value = max(qmin,qci(i,k,1))
1119            source=(praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
1120            if (source.gt.value) then
1121              factor = value/source
1122              praut(i,k) = praut(i,k)*factor
1123              pracw(i,k) = pracw(i,k)*factor
1124              psacw(i,k) = psacw(i,k)*factor
1125            endif
1126!
1127!     rain
1128!
1129            value = max(qmin,qrs(i,k,1))
1130            source = (-praut(i,k)-pracw(i,k)-prevp(i,k)-psacw(i,k))*dtcld
1131            if (source.gt.value) then
1132              factor = value/source
1133              praut(i,k) = praut(i,k)*factor
1134              pracw(i,k) = pracw(i,k)*factor
1135              prevp(i,k) = prevp(i,k)*factor
1136              psacw(i,k) = psacw(i,k)*factor
1137            endif 
1138!
1139!     snow
1140!
1141            value = max(qcrmin,qrs(i,k,2))
1142            source=(-psevp(i,k))*dtcld
1143            if (source.gt.value) then
1144              factor = value/source
1145              psevp(i,k) = psevp(i,k)*factor
1146            endif
1147            work2(i,k)=-(prevp(i,k)+psevp(i,k))
1148!     update
1149            q(i,k,j) = q(i,k,j)+work2(i,k)*dtcld
1150            qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k)                 &
1151                        +psacw(i,k))*dtcld,0.)
1152            qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k)                 &
1153                        +prevp(i,k) +psacw(i,k))*dtcld,0.)
1154            qrs(i,k,2) = max(qrs(i,k,2)+psevp(i,k)*dtcld,0.)
1155            xlf = xls-xl(i,k)
1156            xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k))
1157            t(i,k,j) = t(i,k,j)-xlwork2/cpm(i,k)*dtcld
1158          endif
1159      enddo
1160!
1161! Inline expansion for fpvs
1162!         qs(i,k,1) = fpvs(t(i,k,j),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1163!         qs(i,k,2) = fpvs(t(i,k,j),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1164      hsub = xls
1165      hvap = xlv0
1166      cvap = cpv
1167      ttp=t0c+0.01
1168      dldt=cvap-cliq
1169      xa=-dldt/rv
1170      xb=xa+hvap/(rv*ttp)
1171      dldti=cvap-cice
1172      xai=-dldti/rv
1173      xbi=xai+hsub/(rv*ttp)
1174      do k = kts, kte
1175          tr=ttp/t(i,k,j)
1176          logtr = log(tr)
1177          qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
1178          qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k,j) - qs(i,k,1))
1179          qs(i,k,1) = max(qs(i,k,1),qmin)
1180      enddo
1181!
1182!----------------------------------------------------------------
1183!  pcond: condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
1184!     if there exists additional water vapor condensated/if
1185!     evaporation of cloud water is not enough to remove subsaturation
1186!
1187      do k = kts, kte
1188!         work1(i,k,1) = conden(t(i,k,j),q(i,k,j),qs(i,k,1),xl(i,k),cpm(i,k))
1189          work1(i,k,1) = ((max(q(i,k,j),qmin)-(qs(i,k,1)))/(1.+(xl(i,k))         &   
1190                        *(xl(i,k))/(rv*(cpm(i,k)))*(qs(i,k,1))                 &
1191                        /((t(i,k,j))*(t(i,k,j)))))
1192          work2(i,k) = qci(i,k,1)+work1(i,k,1)
1193          pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k,j),0.)/dtcld)
1194          if(qci(i,k,1).gt.0..and.work1(i,k,1).lt.0.)                          &
1195            pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
1196          q(i,k,j) = q(i,k,j)-pcond(i,k)*dtcld
1197          qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
1198          t(i,k,j) = t(i,k,j)+pcond(i,k)*xl(i,k)/cpm(i,k)*dtcld
1199      enddo
1200!
1201!
1202!----------------------------------------------------------------
1203!     padding for small values
1204!
1205      do k = kts, kte
1206          if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
1207          if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
1208      enddo
1209      enddo                  ! big loops
1210
1211         DO K=kts,kte
1212            th(i,k,j)=t(i,k,j)/pii(i,k,j)
1213            qc(i,k,j) = qci(i,k,1)
1214            qi(i,k,j) = qci(i,k,2)
1215            qr(i,k,j) = qrs(i,k,1)
1216            qqs(i,k,j) = qrs(i,k,2)
1217         ENDDO
1218
1219
1220
1221      ENDDO     !  i loop
1222      enddo     !  j loop
1223!$acc end region
1224
1225    ELSE
1226
1227!
1228! Moved outside of accelerator region
1229!
1230      loops = max(nint(delt/dtcldcr),1)
1231      dtcld = delt/loops
1232      if(delt.le.dtcldcr) dtcld = delt
1233
1234!$acc region &
1235!$acc        local(t) &
1236!$acc        copyin(delz(:,:,:),p(:,:,:),den(:,:,:),pii(:,:,:)) &
1237!$acc        copyout(rainncv(:,:),sr(:,:)) &
1238!$acc        copy(qqs(:,:,:),qr(:,:,:),qi(:,:,:),qc(:,:,:)) &
1239!$acc        copy(th(:,:,:),q(:,:,:),rain(:,:))
1240!$acc do &
1241!$acc        private(rh,qs,rslope,rslope2,rslope3,rslopeb,falk,fall) &
1242!$acc        private(work1,qci,qrs,falkc,fallc,xl,cpm,denfac,xni) &
1243!$acc        private(n0sfac,work2,work1c,work2c,pigen,pidep,psdep) &
1244!$acc        private(praut,psaut,prevp,psevp) &
1245!$acc        private(pracw,psacw,psaci,pcond,psmlt) &
1246!$acc        parallel
1247      do j = jts, jte
1248!$acc do &
1249!$acc        private(numdt,mstep) &
1250!$acc        kernel vector
1251        do i = its, ite
1252        do k = kts, kte
1253          t(i,k,j)=th(i,k,j)*pii(i,k,j)
1254          qci(i,k,1) = max(qc(i,k,j),0.0)
1255          qci(i,k,2) = max(qi(i,k,j),0.0)
1256          qrs(i,k,1) = max(qr(i,k,j),0.0)
1257          qrs(i,k,2) = max(qqs(i,k,j),0.0)
1258        enddo
1259!
1260!----------------------------------------------------------------
1261!     latent heat for phase changes and heat capacity. neglect the
1262!     changes during microphysical process calculation
1263!     emanuel(1994)
1264!
1265      do k = kts, kte
1266          cpm(i,k) = cpmcal(q(i,k,j))
1267          xl(i,k) = xlcal(t(i,k,j))
1268      enddo
1269!
1270!----------------------------------------------------------------
1271!     compute the minor time steps.
1272!
1273!     loops = max(nint(delt/dtcldcr),1)
1274!     dtcld = delt/loops
1275!     if(delt.le.dtcldcr) dtcld = delt
1276!
1277      do loop = 1,loops
1278!
1279!----------------------------------------------------------------
1280!     initialize the large scale variables
1281!
1282        mstep = 1
1283        flgcld = .true.
1284!
1285      do k = kts, kte
1286          denfac(i,k) = sqrt(den0/den(i,k,j))
1287      enddo
1288!     do k = kts, kte
1289!       CALL VREC( tvec1(its), den(its,k,j), ite-its+1)
1290!       do i = its, ite
1291!         tvec1(i) = tvec1(i)*den0
1292!       enddo
1293!       CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
1294!     enddo
1295!
1296! Inline expansion for fpvs
1297!         qs(i,k,1) = fpvs(t(i,k,j),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1298!         qs(i,k,2) = fpvs(t(i,k,j),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1299      hsub = xls
1300      hvap = xlv0
1301      cvap = cpv
1302      ttp=t0c+0.01
1303      dldt=cvap-cliq
1304      xa=-dldt/rv
1305      xb=xa+hvap/(rv*ttp)
1306      dldti=cvap-cice
1307      xai=-dldti/rv
1308      xbi=xai+hsub/(rv*ttp)
1309
1310! this is for compilers where the conditional inhibits vectorization
1311#ifdef WSM_NO_CONDITIONAL_IN_VECTOR
1312      do k = kts, kte
1313          if(t(i,k,j).lt.ttp) then
1314            xal = xai
1315            xbl = xbi
1316          else
1317            xal = xa
1318            xbl = xb
1319          endif
1320          tr=ttp/t(i,k,j)
1321          logtr=log(tr)
1322          qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
1323          qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k,j) - qs(i,k,1))
1324          qs(i,k,1) = max(qs(i,k,1),qmin)
1325          rh(i,k,1) = max(q(i,k,j) / qs(i,k,1),qmin)
1326          qs(i,k,2)=psat*exp(logtr*(xal)+xbl*(1.-tr))
1327          qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k,j) - qs(i,k,2))
1328          qs(i,k,2) = max(qs(i,k,2),qmin)
1329          rh(i,k,2) = max(q(i,k,j) / qs(i,k,2),qmin)
1330      enddo
1331#else
1332      do k = kts, kte
1333          tr=ttp/t(i,k,j)
1334          logtr=log(tr)
1335          qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
1336          qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k,j) - qs(i,k,1))
1337          qs(i,k,1) = max(qs(i,k,1),qmin)
1338          rh(i,k,1) = max(q(i,k,j) / qs(i,k,1),qmin)
1339          if(t(i,k,j).lt.ttp) then
1340            qs(i,k,2)=psat*exp(logtr*(xai)+xbi*(1.-tr))
1341          else
1342            qs(i,k,2)=psat*exp(logtr*(xa)+xb*(1.-tr))
1343          endif
1344          qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k,j) - qs(i,k,2))
1345          qs(i,k,2) = max(qs(i,k,2),qmin)
1346          rh(i,k,2) = max(q(i,k,j) / qs(i,k,2),qmin)
1347      enddo
1348#endif
1349!
1350!----------------------------------------------------------------
1351!     initialize the variables for microphysical physics
1352!
1353!
1354      do k = kts, kte
1355          prevp(i,k) = 0.
1356          psdep(i,k) = 0.
1357          praut(i,k) = 0.
1358          psaut(i,k) = 0.
1359          pracw(i,k) = 0.
1360          psaci(i,k) = 0.
1361          psacw(i,k) = 0.
1362          pigen(i,k) = 0.
1363          pidep(i,k) = 0.
1364          pcond(i,k) = 0.
1365          psmlt(i,k) = 0.
1366          psevp(i,k) = 0.
1367          falk(i,k,1) = 0.
1368          falk(i,k,2) = 0.
1369          fall(i,k,1) = 0.
1370          fall(i,k,2) = 0.
1371          fallc(i,k) = 0.
1372          falkc(i,k) = 0.
1373          xni(i,k) = 1.e3
1374      enddo
1375!
1376!----------------------------------------------------------------
1377!     compute the fallout term:
1378!     first, vertical terminal velosity for minor loops
1379!
1380      do k = kts, kte
1381          supcol = t0c-t(i,k,j)
1382!---------------------------------------------------------------
1383! n0s: Intercept parameter for snow [m-4] [HDC 6]
1384!---------------------------------------------------------------
1385          n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
1386          if(qrs(i,k,1).le.qcrmin)then
1387            rslope(i,k,1) = rslopermax
1388            rslopeb(i,k,1) = rsloperbmax
1389            rslope2(i,k,1) = rsloper2max
1390            rslope3(i,k,1) = rsloper3max
1391          else
1392            rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k,j))
1393            rslopeb(i,k,1) = exp(log(rslope(i,k,1))*(bvtr))
1394            rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
1395            rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
1396          endif
1397          if(qrs(i,k,2).le.qcrmin)then
1398            rslope(i,k,2) = rslopesmax
1399            rslopeb(i,k,2) = rslopesbmax
1400            rslope2(i,k,2) = rslopes2max
1401            rslope3(i,k,2) = rslopes3max
1402          else
1403            rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k,j),n0sfac(i,k))
1404            rslopeb(i,k,2) = exp(log(rslope(i,k,2))*(bvts))
1405            rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
1406            rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
1407          endif
1408!-------------------------------------------------------------
1409! Ni: ice crystal number concentraiton   [HDC 5c]
1410!-------------------------------------------------------------
1411!         xni(i,k) = min(max(5.38e7*(den(i,k,j)                                  &
1412!                   *max(qci(i,k,2),qmin))**0.75,1.e3),1.e6)
1413          temp = (den(i,k,j)*max(qci(i,k,2),qmin))
1414          temp = sqrt(sqrt(temp*temp*temp))
1415          xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
1416      enddo
1417!
1418      numdt = 1
1419      do k = kte, kts, -1
1420          work1(i,k,1) = pvtr*rslopeb(i,k,1)*denfac(i,k)/delz(i,k,j)
1421          work1(i,k,2) = pvts*rslopeb(i,k,2)*denfac(i,k)/delz(i,k,j)
1422          numdt = max(nint(max(work1(i,k,1),work1(i,k,2))*dtcld+.5),1)
1423          if(numdt.ge.mstep) mstep = numdt
1424      enddo
1425        rmstep = 1./mstep
1426!
1427      do n = 1, mstep
1428        k = kte
1429!             falk(i,k,1) = den(i,k,j)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
1430!             falk(i,k,2) = den(i,k,j)*qrs(i,k,2)*work1(i,k,2)/mstep(i)
1431              falk(i,k,1) = den(i,k,j)*qrs(i,k,1)*work1(i,k,1)*rmstep
1432              falk(i,k,2) = den(i,k,j)*qrs(i,k,2)*work1(i,k,2)*rmstep
1433              fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
1434              fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
1435!             qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcld/den(i,k,j),0.)
1436!             qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcld/den(i,k,j),0.)
1437              dtcldden = dtcld/den(i,k,j)
1438              qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcldden,0.)
1439              qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcldden,0.)
1440!           endif
1441        do k = kte-1, kts, -1
1442              falk(i,k,1) = den(i,k,j)*qrs(i,k,1)*work1(i,k,1)*rmstep
1443              falk(i,k,2) = den(i,k,j)*qrs(i,k,2)*work1(i,k,2)*rmstep
1444              fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
1445              fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
1446              dtcldden = dtcld/den(i,k,j)
1447              rdelz = 1./delz(i,k,j)
1448              qrs(i,k,1) = max(qrs(i,k,1)-(falk(i,k,1)-falk(i,k+1,1)           &
1449                          *delz(i,k+1,j)*rdelz)*dtcldden,0.)
1450              qrs(i,k,2) = max(qrs(i,k,2)-(falk(i,k,2)-falk(i,k+1,2)           &
1451                          *delz(i,k+1,j)*rdelz)*dtcldden,0.)
1452        enddo
1453        do k = kte, kts, -1
1454              if(t(i,k,j).gt.t0c.and.qrs(i,k,2).gt.0.) then
1455!----------------------------------------------------------------
1456! psmlt: melting of snow [HL A33] [RH83 A25]
1457!       (T>T0: S->R)
1458!----------------------------------------------------------------
1459                xlf = xlf0
1460!               work2(i,k)= venfac(p(i,k),t(i,k,j),den(i,k,j))
1461                work2(i,k)= (exp(log(((1.496e-6*((t(i,k,j))*sqrt(t(i,k,j)))        &
1462                            /((t(i,k,j))+120.)/(den(i,k,j)))/(8.794e-5             &
1463                            *exp(log(t(i,k,j))*(1.81))/p(i,k,j))))                 &
1464                            *((.3333333)))/sqrt((1.496e-6*((t(i,k,j))            &
1465                            *sqrt(t(i,k,j)))/((t(i,k,j))+120.)/(den(i,k,j))))        &
1466                            *sqrt(sqrt(den0/(den(i,k,j)))))
1467                coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
1468!               psmlt(i,k) = xka(t(i,k,j),den(i,k,j))/xlf*(t0c-t(i,k,j))*pi/2.       &
1469!                           *n0sfac(i,k)*(precs1*rslope2(i,k,2)+precs2         &
1470!                           *work2(i,k)*coeres)
1471                psmlt(i,k) = (1.414e3*(1.496e-6*((t(i,k,j))*sqrt(t(i,k,j)))        &         
1472                            /((t(i,k,j))+120.)/(den(i,k,j)) )*(den(i,k,j)))          &
1473                            /xlf*(t0c-t(i,k,j))*pi/2.                            &
1474                            *n0sfac(i,k)*(precs1*rslope2(i,k,2)+precs2         &
1475                            *work2(i,k)*coeres)
1476                psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep,                &
1477                            -qrs(i,k,2)/mstep),0.)
1478                qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
1479                qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k)
1480                t(i,k,j) = t(i,k,j) + xlf/cpm(i,k)*psmlt(i,k)
1481              endif
1482        enddo
1483      enddo
1484!---------------------------------------------------------------
1485! Vice [ms-1] : fallout of ice crystal [HDC 5a]
1486!---------------------------------------------------------------
1487      mstep = 1
1488      numdt = 1
1489      do k = kte, kts, -1
1490          if(qci(i,k,2).le.0.) then
1491            work2c(i,k) = 0.
1492          else
1493            xmi = den(i,k,j)*qci(i,k,2)/xni(i,k)
1494!           diameter  = min(dicon * sqrt(xmi),dimax)
1495            diameter  = max(min(dicon * sqrt(xmi),dimax), 1.e-25)
1496            work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31))
1497            work2c(i,k) = work1c(i,k)/delz(i,k,j)
1498          endif
1499          numdt = max(nint(work2c(i,k)*dtcld+.5),1)
1500          if(numdt.ge.mstep) mstep = numdt
1501      enddo
1502!
1503      do n = 1, mstep
1504        k = kte
1505            falkc(i,k) = den(i,k,j)*qci(i,k,2)*work2c(i,k)/mstep
1506            holdc = falkc(i,k)
1507            fallc(i,k) = fallc(i,k)+falkc(i,k)
1508            holdci = qci(i,k,2)
1509            qci(i,k,2) = max(qci(i,k,2)-falkc(i,k)*dtcld/den(i,k,j),0.)
1510        do k = kte-1, kts, -1
1511              falkc(i,k) = den(i,k,j)*qci(i,k,2)*work2c(i,k)/mstep
1512              holdc = falkc(i,k)
1513              fallc(i,k) = fallc(i,k)+falkc(i,k)
1514              holdci = qci(i,k,2)
1515              qci(i,k,2) = max(qci(i,k,2)-(falkc(i,k)-falkc(i,k+1)             &
1516                          *delz(i,k+1,j)/delz(i,k,j))*dtcld/den(i,k,j),0.)
1517        enddo
1518      enddo
1519!
1520!
1521!----------------------------------------------------------------
1522!      rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
1523!
1524        fallsum = fall(i,1,1)+fall(i,1,2)+fallc(i,1)
1525        fallsum_qsi = fall(i,1,2)+fallc(i,1)
1526        rainncv(i,j) = 0.
1527        if(fallsum.gt.0.) then
1528          rainncv(i,j) = fallsum*delz(i,1,j)/denr*dtcld*1000.
1529          rain(i,j) = fallsum*delz(i,1,j)/denr*dtcld*1000. + rain(i,j)
1530        endif
1531        sr(i,j) = 0.
1532        if(fallsum.gt.0.)sr(i,j)=fallsum_qsi*delz(i,kts,j)/denr*dtcld*1000.        &   
1533        /(rainncv(i,j)+1.e-12)
1534!
1535!---------------------------------------------------------------
1536! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
1537!       (T>T0: I->C)
1538!---------------------------------------------------------------
1539      do k = kts, kte
1540          supcol = t0c-t(i,k,j)
1541          xlf = xls-xl(i,k)
1542          if(supcol.lt.0.) xlf = xlf0
1543          if(supcol.lt.0.and.qci(i,k,2).gt.0.) then
1544            qci(i,k,1) = qci(i,k,1) + qci(i,k,2)
1545            t(i,k,j) = t(i,k,j) - xlf/cpm(i,k)*qci(i,k,2)
1546            qci(i,k,2) = 0.
1547          endif
1548!---------------------------------------------------------------
1549! pihmf: homogeneous freezing of cloud water below -40c [HL A45]
1550!        (T<-40C: C->I)
1551!---------------------------------------------------------------
1552          if(supcol.gt.40..and.qci(i,k,1).gt.0.) then
1553            qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
1554            t(i,k,j) = t(i,k,j) + xlf/cpm(i,k)*qci(i,k,1)
1555            qci(i,k,1) = 0.
1556          endif
1557!---------------------------------------------------------------
1558! pihtf: heterogeneous freezing of cloud water [HL A44]
1559!        (T0>T>-40C: C->I)
1560!---------------------------------------------------------------
1561          if(supcol.gt.0..and.qci(i,k,1).gt.0.) then
1562            supcolt=min(supcol,50.)
1563!           pfrzdtc = min(pfrz1*(exp(pfrz2*supcol)-1.)                         &
1564!              *den(i,k,j)/denr/xncr*qci(i,k,1)**2*dtcld,qci(i,k,1))
1565            pfrzdtc = min(pfrz1*(exp(pfrz2*supcolt)-1.)                        &
1566            *den(i,k,j)/denr/xncr*qci(i,k,1)*qci(i,k,1)*dtcld,qci(i,k,1))
1567            qci(i,k,2) = qci(i,k,2) + pfrzdtc
1568            t(i,k,j) = t(i,k,j) + xlf/cpm(i,k)*pfrzdtc
1569            qci(i,k,1) = qci(i,k,1)-pfrzdtc
1570          endif
1571!---------------------------------------------------------------
1572! psfrz: freezing of rain water [HL A20] [LFO 45]
1573!        (T<T0, R->S)
1574!---------------------------------------------------------------
1575          if(supcol.gt.0..and.qrs(i,k,1).gt.0.) then
1576            supcolt=min(supcol,50.)
1577!           pfrzdtr = min(20.*pi**2*pfrz1*n0r*denr/den(i,k,j)                    &
1578!                 *(exp(pfrz2*supcol)-1.)*rslope(i,k,1)**7*dtcld,              &
1579!                 qrs(i,k,1))
1580            temp = rslope(i,k,1)
1581            temp = temp*temp*temp*temp*temp*temp*temp
1582            pfrzdtr = min(20.*(pi*pi)*pfrz1*n0r*denr/den(i,k,j)                  &
1583                  *(exp(pfrz2*supcolt)-1.)*temp*dtcld,                         &
1584                  qrs(i,k,1))
1585            qrs(i,k,2) = qrs(i,k,2) + pfrzdtr
1586            t(i,k,j) = t(i,k,j) + xlf/cpm(i,k)*pfrzdtr
1587            qrs(i,k,1) = qrs(i,k,1)-pfrzdtr
1588          endif
1589      enddo
1590!
1591!----------------------------------------------------------------
1592!     rsloper: reverse of the slope parameter of the rain(m)
1593!     xka:    thermal conductivity of air(jm-1s-1k-1)
1594!     work1:  the thermodynamic term in the denominator associated with
1595!             heat conduction and vapor diffusion
1596!             (ry88, y93, h85)
1597!     work2: parameter associated with the ventilation effects(y93)
1598!
1599      do k = kts, kte
1600          if(qrs(i,k,1).le.qcrmin)then
1601            rslope(i,k,1) = rslopermax
1602            rslopeb(i,k,1) = rsloperbmax
1603            rslope2(i,k,1) = rsloper2max
1604            rslope3(i,k,1) = rsloper3max
1605          else
1606!           rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k,j))
1607            rslope(i,k,1) = 1./(sqrt(sqrt(pidn0r/((qrs(i,k,1))*(den(i,k,j))))))
1608            rslopeb(i,k,1) = exp(log(rslope(i,k,1))*(bvtr))
1609            rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
1610            rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
1611          endif
1612          if(qrs(i,k,2).le.qcrmin)then
1613            rslope(i,k,2) = rslopesmax
1614            rslopeb(i,k,2) = rslopesbmax
1615            rslope2(i,k,2) = rslopes2max
1616            rslope3(i,k,2) = rslopes3max
1617          else
1618!            rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k,j),n0sfac(i,k))
1619            rslope(i,k,2) = 1./(sqrt(sqrt(pidn0s*(n0sfac(i,k))/((qrs(i,k,2))   & 
1620                          *(den(i,k,j))))))
1621            rslopeb(i,k,2) = exp(log(rslope(i,k,2))*(bvts))
1622            rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
1623            rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
1624          endif
1625      enddo
1626!
1627      do k = kts, kte
1628!         work1(i,k,1) = diffac(xl(i,k),p(i,k,j),t(i,k,j),den(i,k,j),qs(i,k,1))
1629          work1(i,k,1) = ((((den(i,k,j))*(xl(i,k))*(xl(i,k)))*((t(i,k,j))+120.)    &
1630                       *(den(i,k,j)))/(1.414e3*(1.496e-6*((t(i,k,j))*sqrt(t(i,k,j))))&
1631                       *(den(i,k,j))*(rv*(t(i,k,j))*(t(i,k,j)))))                    &
1632                      +  p(i,k,j)/((qs(i,k,1))*(8.794e-5*exp(log(t(i,k,j))*(1.81))))
1633!         work1(i,k,2) = diffac(xls,p(i,k,j),t(i,k,j),den(i,k,j),qs(i,k,2))
1634          work1(i,k,2) = ((((den(i,k,j))*(xls)*(xls))*((t(i,k,j))+120.)*(den(i,k,j)))&
1635                       /(1.414e3*(1.496e-6*((t(i,k,j))*sqrt(t(i,k,j))))*(den(i,k,j)) &
1636                       *(rv*(t(i,k,j))*(t(i,k,j))))                                &
1637                      + p(i,k,j)/(qs(i,k,2)*(8.794e-5*exp(log(t(i,k,j))*(1.81)))))
1638!         work2(i,k) = venfac(p(i,k,j),t(i,k,j),den(i,k,j))
1639          work2(i,k) = (exp(.3333333*log(((1.496e-6 * ((t(i,k,j))*sqrt(t(i,k,j)))) &
1640                     *p(i,k,j))/(((t(i,k,j))+120.)*den(i,k,j)*(8.794e-5              &
1641                     *exp(log(t(i,k,j))*(1.81))))))*sqrt(sqrt(den0/(den(i,k,j))))) &
1642                      /sqrt((1.496e-6*((t(i,k,j))*sqrt(t(i,k,j))))                 &
1643                      /(((t(i,k,j))+120.)*den(i,k,j)))
1644      enddo
1645!
1646!===============================================================
1647!
1648! warm rain processes
1649!
1650! - follows the processes in RH83 and LFO except for autoconcersion
1651!
1652!===============================================================
1653!
1654      do k = kts, kte
1655          supsat = max(q(i,k,j),qmin)-qs(i,k,1)
1656          satdt = supsat/dtcld
1657!---------------------------------------------------------------
1658! praut: auto conversion rate from cloud to rain [HDC 16]
1659!        (C->R)
1660!---------------------------------------------------------------
1661          if(qci(i,k,1).gt.qc0) then
1662            praut(i,k) = qck1*exp(log(qci(i,k,1))*((7./3.)))
1663            praut(i,k) = min(praut(i,k),qci(i,k,1)/dtcld)
1664          endif
1665!---------------------------------------------------------------
1666! pracw: accretion of cloud water by rain [HL A40] [LFO 51]
1667!        (C->R)
1668!---------------------------------------------------------------
1669          if(qrs(i,k,1).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
1670            pracw(i,k) = min(pacrr*rslope3(i,k,1)*rslopeb(i,k,1)               &
1671                         *qci(i,k,1)*denfac(i,k),qci(i,k,1)/dtcld)
1672          endif
1673!---------------------------------------------------------------
1674! prevp: evaporation/condensation rate of rain [HDC 14]
1675!        (V->R or R->V)
1676!---------------------------------------------------------------
1677          if(qrs(i,k,1).gt.0.) then
1678            coeres = rslope2(i,k,1)*sqrt(rslope(i,k,1)*rslopeb(i,k,1))
1679            prevp(i,k) = (rh(i,k,1)-1.)*(precr1*rslope2(i,k,1)                 &
1680                         +precr2*work2(i,k)*coeres)/work1(i,k,1)
1681            if(prevp(i,k).lt.0.) then
1682              prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
1683              prevp(i,k) = max(prevp(i,k),satdt/2)
1684            else
1685              prevp(i,k) = min(prevp(i,k),satdt/2)
1686            endif
1687          endif
1688      enddo
1689!
1690!===============================================================
1691!
1692! cold rain processes
1693!
1694! - follows the revised ice microphysics processes in HDC
1695! - the processes same as in RH83 and RH84  and LFO behave
1696!   following ice crystal hapits defined in HDC, inclduing
1697!   intercept parameter for snow (n0s), ice crystal number
1698!   concentration (ni), ice nuclei number concentration
1699!   (n0i), ice diameter (d)
1700!
1701!===============================================================
1702!
1703      rdtcld = 1./dtcld
1704      do k = kts, kte
1705          supcol = t0c-t(i,k,j)
1706          supsat = max(q(i,k,j),qmin)-qs(i,k,2)
1707          satdt = supsat/dtcld
1708          ifsat = 0
1709!-------------------------------------------------------------
1710! Ni: ice crystal number concentraiton   [HDC 5c]
1711!-------------------------------------------------------------
1712!         xni(i,k) = min(max(5.38e7*(den(i,k,j)                                  &
1713!                      *max(qci(i,k,2),qmin))**0.75,1.e3),1.e6)
1714          temp = (den(i,k,j)*max(qci(i,k,2),qmin))
1715          temp = sqrt(sqrt(temp*temp*temp))
1716          xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
1717          eacrs = exp(0.07*(-supcol))
1718!
1719          if(supcol.gt.0) then
1720            if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,2).gt.qmin) then
1721              xmi = den(i,k,j)*qci(i,k,2)/xni(i,k)
1722              diameter  = min(dicon * sqrt(xmi),dimax)
1723              vt2i = 1.49e4*diameter**1.31
1724              vt2s = pvts*rslopeb(i,k,2)*denfac(i,k)
1725!-------------------------------------------------------------
1726! psaci: Accretion of cloud ice by rain [HDC 10]
1727!        (T<T0: I->S)
1728!-------------------------------------------------------------
1729              acrfac = 2.*rslope3(i,k,2)+2.*diameter*rslope2(i,k,2)            &
1730                      +diameter**2*rslope(i,k,2)
1731              psaci(i,k) = pi*qci(i,k,2)*eacrs*n0s*n0sfac(i,k)                 &
1732                           *abs(vt2s-vt2i)*acrfac/4.
1733            endif
1734          endif
1735!-------------------------------------------------------------
1736! psacw: Accretion of cloud water by snow  [HL A7] [LFO 24]
1737!        (T<T0: C->S, and T>=T0: C->R)
1738!-------------------------------------------------------------
1739          if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
1740            psacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2)                  &
1741                           *rslopeb(i,k,2)*qci(i,k,1)*denfac(i,k)              &
1742!                          ,qci(i,k,1)/dtcld)
1743                           ,qci(i,k,1)*rdtcld)
1744          endif
1745          if(supcol .gt. 0) then
1746!-------------------------------------------------------------
1747! pidep: Deposition/Sublimation rate of ice [HDC 9]
1748!       (T<T0: V->I or I->V)
1749!-------------------------------------------------------------
1750            if(qci(i,k,2).gt.0.and.ifsat.ne.1) then
1751              xmi = den(i,k,j)*qci(i,k,2)/xni(i,k)
1752              diameter = dicon * sqrt(xmi)
1753              pidep(i,k) = 4.*diameter*xni(i,k)*(rh(i,k,2)-1.)/work1(i,k,2)
1754              supice = satdt-prevp(i,k)
1755              if(pidep(i,k).lt.0.) then
1756!               pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
1757!               pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
1758                pidep(i,k) = max(max(pidep(i,k),satdt*.5),supice)
1759                pidep(i,k) = max(pidep(i,k),-qci(i,k,2)*rdtcld)
1760              else
1761!               pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
1762                pidep(i,k) = min(min(pidep(i,k),satdt*.5),supice)
1763              endif
1764              if(abs(prevp(i,k)+pidep(i,k)).ge.abs(satdt)) ifsat = 1
1765            endif
1766!-------------------------------------------------------------
1767! psdep: deposition/sublimation rate of snow [HDC 14]
1768!        (V->S or S->V)
1769!-------------------------------------------------------------
1770            if(qrs(i,k,2).gt.0..and.ifsat.ne.1) then
1771              coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
1772              psdep(i,k) = (rh(i,k,2)-1.)*n0sfac(i,k)                          &
1773                           *(precs1*rslope2(i,k,2)+precs2                      &
1774                           *work2(i,k)*coeres)/work1(i,k,2)
1775              supice = satdt-prevp(i,k)-pidep(i,k)
1776              if(psdep(i,k).lt.0.) then
1777!               psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
1778!               psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
1779                psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)*rdtcld)
1780                psdep(i,k) = max(max(psdep(i,k),satdt*.5),supice)
1781              else
1782!             psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
1783                psdep(i,k) = min(min(psdep(i,k),satdt*.5),supice)
1784              endif
1785              if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt))          &
1786                ifsat = 1
1787            endif
1788!-------------------------------------------------------------
1789! pigen: generation(nucleation) of ice from vapor [HL A50] [HDC 7-8]
1790!       (T<T0: V->I)
1791!-------------------------------------------------------------
1792            if(supsat.gt.0.and.ifsat.ne.1) then
1793              supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
1794              xni0 = 1.e3*exp(0.1*supcol)
1795              roqi0 = 4.92e-11*exp(log(xni0)*(1.33))
1796              pigen(i,k) = max(0.,(roqi0/den(i,k,j)-max(qci(i,k,2),0.))          &
1797!                        /dtcld)
1798                         *rdtcld)
1799              pigen(i,k) = min(min(pigen(i,k),satdt),supice)
1800            endif
1801!
1802!-------------------------------------------------------------
1803! psaut: conversion(aggregation) of ice to snow [HDC 12]
1804!       (T<T0: I->S)
1805!-------------------------------------------------------------
1806            if(qci(i,k,2).gt.0.) then
1807              qimax = roqimax/den(i,k,j)
1808!             psaut(i,k) = max(0.,(qci(i,k,2)-qimax)/dtcld)
1809              psaut(i,k) = max(0.,(qci(i,k,2)-qimax)*rdtcld)
1810            endif
1811          endif
1812!-------------------------------------------------------------
1813! psevp: Evaporation of melting snow [HL A35] [RH83 A27]
1814!       (T>T0: S->V)
1815!-------------------------------------------------------------
1816          if(supcol.lt.0.) then
1817            if(qrs(i,k,2).gt.0..and.rh(i,k,1).lt.1.)                           &
1818              psevp(i,k) = psdep(i,k)*work1(i,k,2)/work1(i,k,1)
1819!              psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)/dtcld),0.)
1820              psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)*rdtcld),0.)
1821          endif
1822      enddo
1823!
1824!
1825!----------------------------------------------------------------
1826!     check mass conservation of generation terms and feedback to the
1827!     large scale
1828!
1829      do k = kts, kte
1830          if(t(i,k,j).le.t0c) then
1831!
1832!     cloud water
1833!
1834            value = max(qmin,qci(i,k,1))
1835            source = (praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
1836            if (source.gt.value) then
1837              factor = value/source
1838              praut(i,k) = praut(i,k)*factor
1839              pracw(i,k) = pracw(i,k)*factor
1840              psacw(i,k) = psacw(i,k)*factor
1841            endif
1842!
1843!     cloud ice
1844!
1845            value = max(qmin,qci(i,k,2))
1846            source = (psaut(i,k)+psaci(i,k)-pigen(i,k)-pidep(i,k))*dtcld
1847            if (source.gt.value) then
1848              factor = value/source
1849              psaut(i,k) = psaut(i,k)*factor
1850              psaci(i,k) = psaci(i,k)*factor
1851              pigen(i,k) = pigen(i,k)*factor
1852              pidep(i,k) = pidep(i,k)*factor
1853            endif
1854!
1855!     rain
1856!
1857!
1858            value = max(qmin,qrs(i,k,1))
1859            source = (-praut(i,k)-pracw(i,k)-prevp(i,k))*dtcld
1860            if (source.gt.value) then
1861              factor = value/source
1862              praut(i,k) = praut(i,k)*factor
1863              pracw(i,k) = pracw(i,k)*factor
1864              prevp(i,k) = prevp(i,k)*factor
1865            endif
1866!
1867!    snow
1868!
1869            value = max(qmin,qrs(i,k,2))
1870            source = (-psdep(i,k)-psaut(i,k)-psaci(i,k)-psacw(i,k))*dtcld 
1871            if (source.gt.value) then
1872              factor = value/source
1873              psdep(i,k) = psdep(i,k)*factor
1874              psaut(i,k) = psaut(i,k)*factor
1875              psaci(i,k) = psaci(i,k)*factor
1876              psacw(i,k) = psacw(i,k)*factor
1877            endif
1878!
1879            work2(i,k)=-(prevp(i,k)+psdep(i,k)+pigen(i,k)+pidep(i,k))
1880!     update
1881            q(i,k,j) = q(i,k,j)+work2(i,k)*dtcld
1882            qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k)                 &
1883                        +psacw(i,k))*dtcld,0.)
1884            qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k)                 &
1885                        +prevp(i,k))*dtcld,0.)
1886            qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+psaci(i,k)                 &
1887                        -pigen(i,k)-pidep(i,k))*dtcld,0.)
1888            qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k)                 &
1889                        +psaci(i,k)+psacw(i,k))*dtcld,0.)
1890            xlf = xls-xl(i,k)
1891            xlwork2 = -xls*(psdep(i,k)+pidep(i,k)+pigen(i,k))                  &
1892                      -xl(i,k)*prevp(i,k)-xlf*psacw(i,k)
1893            t(i,k,j) = t(i,k,j)-xlwork2/cpm(i,k)*dtcld
1894          else
1895!
1896!     cloud water
1897!
1898            value = max(qmin,qci(i,k,1))
1899            source=(praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
1900            if (source.gt.value) then
1901              factor = value/source
1902              praut(i,k) = praut(i,k)*factor
1903              pracw(i,k) = pracw(i,k)*factor
1904              psacw(i,k) = psacw(i,k)*factor
1905            endif
1906!
1907!     rain
1908!
1909            value = max(qmin,qrs(i,k,1))
1910            source = (-praut(i,k)-pracw(i,k)-prevp(i,k)-psacw(i,k))*dtcld
1911            if (source.gt.value) then
1912              factor = value/source
1913              praut(i,k) = praut(i,k)*factor
1914              pracw(i,k) = pracw(i,k)*factor
1915              prevp(i,k) = prevp(i,k)*factor
1916              psacw(i,k) = psacw(i,k)*factor
1917            endif 
1918!
1919!     snow
1920!
1921            value = max(qcrmin,qrs(i,k,2))
1922            source=(-psevp(i,k))*dtcld
1923            if (source.gt.value) then
1924              factor = value/source
1925              psevp(i,k) = psevp(i,k)*factor
1926            endif
1927            work2(i,k)=-(prevp(i,k)+psevp(i,k))
1928!     update
1929            q(i,k,j) = q(i,k,j)+work2(i,k)*dtcld
1930            qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k)                 &
1931                        +psacw(i,k))*dtcld,0.)
1932            qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k)                 &
1933                        +prevp(i,k) +psacw(i,k))*dtcld,0.)
1934            qrs(i,k,2) = max(qrs(i,k,2)+psevp(i,k)*dtcld,0.)
1935            xlf = xls-xl(i,k)
1936            xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k))
1937            t(i,k,j) = t(i,k,j)-xlwork2/cpm(i,k)*dtcld
1938          endif
1939      enddo
1940!
1941! Inline expansion for fpvs
1942!         qs(i,k,1) = fpvs(t(i,k,j),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1943!         qs(i,k,2) = fpvs(t(i,k,j),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1944      hsub = xls
1945      hvap = xlv0
1946      cvap = cpv
1947      ttp=t0c+0.01
1948      dldt=cvap-cliq
1949      xa=-dldt/rv
1950      xb=xa+hvap/(rv*ttp)
1951      dldti=cvap-cice
1952      xai=-dldti/rv
1953      xbi=xai+hsub/(rv*ttp)
1954      do k = kts, kte
1955          tr=ttp/t(i,k,j)
1956          logtr = log(tr)
1957          qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
1958          qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k,j) - qs(i,k,1))
1959          qs(i,k,1) = max(qs(i,k,1),qmin)
1960      enddo
1961!
1962!----------------------------------------------------------------
1963!  pcond: condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
1964!     if there exists additional water vapor condensated/if
1965!     evaporation of cloud water is not enough to remove subsaturation
1966!
1967      do k = kts, kte
1968!         work1(i,k,1) = conden(t(i,k,j),q(i,k,j),qs(i,k,1),xl(i,k),cpm(i,k))
1969          work1(i,k,1) = ((max(q(i,k,j),qmin)-(qs(i,k,1)))/(1.+(xl(i,k))         &   
1970                        *(xl(i,k))/(rv*(cpm(i,k)))*(qs(i,k,1))                 &
1971                        /((t(i,k,j))*(t(i,k,j)))))
1972          work2(i,k) = qci(i,k,1)+work1(i,k,1)
1973          pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k,j),0.)/dtcld)
1974          if(qci(i,k,1).gt.0..and.work1(i,k,1).lt.0.)                          &
1975            pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
1976          q(i,k,j) = q(i,k,j)-pcond(i,k)*dtcld
1977          qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
1978          t(i,k,j) = t(i,k,j)+pcond(i,k)*xl(i,k)/cpm(i,k)*dtcld
1979      enddo
1980!
1981!
1982!----------------------------------------------------------------
1983!     padding for small values
1984!
1985      do k = kts, kte
1986          if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
1987          if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
1988      enddo
1989      enddo                  ! big loops
1990
1991         DO K=kts,kte
1992            th(i,k,j)=t(i,k,j)/pii(i,k,j)
1993            qc(i,k,j) = qci(i,k,1)
1994            qi(i,k,j) = qci(i,k,2)
1995            qr(i,k,j) = qrs(i,k,1)
1996            qqs(i,k,j) = qrs(i,k,2)
1997         ENDDO
1998
1999      ENDDO     !  i loop
2000      enddo     !  j loop
2001!$acc end region
2002
2003    ENDIF
2004
2005  END SUBROUTINE wsm52d
2006
2007
2008#else
2009
2010!===================================================================
2011!
2012  SUBROUTINE wsm52D(t, q, qci, qrs, den, p, delz                  &
2013                   ,delt,g, cpd, cpv, rd, rv, t0c                 &
2014                   ,ep1, ep2, qmin                                &
2015                   ,XLS, XLV0, XLF0, den0, denr                   &
2016                   ,cliq,cice,psat                                &
2017                   ,lat                                           &
2018                   ,rain,rainncv                                  &
2019                   ,sr                                            &
2020                   ,ids,ide, jds,jde, kds,kde                     &
2021                   ,ims,ime, jms,jme, kms,kme                     &
2022                   ,its,ite, jts,jte, kts,kte                     &
2023                   ,snow,snowncv                                  &
2024                                                                  )
2025!-------------------------------------------------------------------
2026  IMPLICIT NONE
2027!-------------------------------------------------------------------
2028  INTEGER,      INTENT(IN   )    ::   ids,ide, jds,jde, kds,kde , &
2029                                      ims,ime, jms,jme, kms,kme , &
2030                                      its,ite, jts,jte, kts,kte,  &
2031                                      lat
2032  REAL, DIMENSION( its:ite , kts:kte ),                           &
2033        INTENT(INOUT) ::                                          &
2034                                                               t
2035  REAL, DIMENSION( its:ite , kts:kte, 2 ),                        &
2036        INTENT(INOUT) ::                                          &
2037                                                             qci, &
2038                                                             qrs
2039  REAL, DIMENSION( ims:ime , kms:kme ),                           &
2040        INTENT(INOUT) ::                                          &
2041                                                               q
2042  REAL, DIMENSION( ims:ime , kms:kme ),                           &
2043        INTENT(IN   ) ::                                          &
2044                                                             den, &
2045                                                               p, &
2046                                                            delz
2047  REAL, INTENT(IN   ) ::                                    delt, &
2048                                                               g, &
2049                                                             cpd, &
2050                                                             cpv, &
2051                                                             t0c, &
2052                                                            den0, &
2053                                                              rd, &
2054                                                              rv, &
2055                                                             ep1, &
2056                                                             ep2, &
2057                                                            qmin, &
2058                                                             XLS, &
2059                                                            XLV0, &
2060                                                            XLF0, &
2061                                                            cliq, &
2062                                                            cice, &
2063                                                            psat, &
2064                                                            denr
2065  REAL, DIMENSION( ims:ime ),                                     &
2066        INTENT(INOUT) ::                                    rain, &
2067                                                         rainncv, &
2068                                                              sr
2069
2070  REAL, DIMENSION( ims:ime, jms:jme ),     OPTIONAL,              &
2071        INTENT(INOUT) ::                                    snow, &
2072                                                         snowncv
2073
2074! LOCAL VAR
2075  REAL, DIMENSION( its:ite , kts:kte , 2) ::                      &
2076                                                              rh, &
2077                                                              qs, &
2078                                                          rslope, &
2079                                                         rslope2, &
2080                                                         rslope3, &
2081                                                         rslopeb, &
2082                                                            falk, &
2083                                                            fall, &
2084                                                           work1
2085
2086  REAL, DIMENSION( its:ite , kts:kte ) ::                         &
2087                                                           falkc, &
2088                                                           fallc, &
2089                                                              xl, &
2090                                                             cpm, &
2091                                                          denfac, &
2092                                                             xni, &
2093                                                          n0sfac, &
2094                                                           work2, &
2095                                                          work1c, &
2096                                                          work2c
2097
2098  REAL, DIMENSION( its:ite , kts:kte ) ::                         &
2099                                                           pigen, &
2100                                                           pidep, &
2101                                                           psdep, &
2102                                                           praut, &
2103                                                           psaut, &
2104                                                           prevp, &
2105                                                           psevp, &
2106                                                           pracw, &
2107                                                           psacw, &
2108                                                           psaci, &
2109                                                           pcond, &
2110                                                           psmlt
2111  INTEGER, DIMENSION( its:ite ) ::                                &
2112                                                           mstep, &
2113                                                           numdt
2114  REAL, DIMENSION(its:ite) ::                             rmstep
2115  REAL dtcldden, rdelz, rdtcld
2116  LOGICAL, DIMENSION( its:ite ) ::                        flgcld
2117
2118#define WSM_NO_CONDITIONAL_IN_VECTOR
2119#ifdef WSM_NO_CONDITIONAL_IN_VECTOR
2120  REAL, DIMENSION(its:ite) :: xal, xbl
2121#endif
2122
2123  REAL  ::  pi,                                                   &
2124            cpmcal, xlcal, lamdar, lamdas, diffus,                &
2125            viscos, xka, venfac, conden, diffac,                  &
2126            x, y, z, a, b, c, d, e,                               &
2127            qdt, holdrr, holdrs, supcol, supcolt, pvt,            &
2128            coeres, supsat, dtcld, xmi, eacrs, satdt,             &
2129            vt2i,vt2s,acrfac,                                     &
2130            qimax, diameter, xni0, roqi0,                         &
2131            fallsum, fallsum_qsi, xlwork2, factor, source,        &
2132            value, xlf, pfrzdtc, pfrzdtr, supice,  holdc, holdci
2133! variables for optimization
2134  REAL, DIMENSION( its:ite )           ::                    tvec1
2135  REAL ::                                                    temp
2136  INTEGER :: i, j, k, mstepmax,                                   &
2137            iprt, latd, lond, loop, loops, ifsat, n
2138! Temporaries used for inlining fpvs function
2139  REAL  :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
2140  REAL  :: logtr
2141!
2142!=================================================================
2143!   compute internal functions
2144!
2145      cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
2146      xlcal(x) = xlv0-xlv1*(x-t0c)
2147!----------------------------------------------------------------
2148!     size distributions: (x=mixing ratio, y=air density):
2149!     valid for mixing ratio > 1.e-9 kg/kg.
2150!
2151! Optimizatin : A**B => exp(log(A)*(B))
2152      lamdar(x,y)=   sqrt(sqrt(pidn0r/(x*y)))      ! (pidn0r/(x*y))**.25
2153      lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y)))    ! (pidn0s*z/(x*y))**.25
2154!
2155!----------------------------------------------------------------
2156!     diffus: diffusion coefficient of the water vapor
2157!     viscos: kinematic viscosity(m2s-1)
2158!     diffus(x,y) = 8.794e-5 * exp(log(x)*(1.81)) / y        ! 8.794e-5*x**1.81/y
2159!     viscos(x,y) = 1.496e-6 * (x*sqrt(x)) /(x+120.)/y  ! 1.496e-6*x**1.5/(x+120.)/y
2160!     xka(x,y) = 1.414e3*viscos(x,y)*y
2161!     diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b))
2162!     venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333)))         &
2163!                    /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
2164!     conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+d*d/(rv*e)*c/(a*a))
2165!
2166!
2167      pi = 4. * atan(1.)
2168!
2169!----------------------------------------------------------------
2170!     paddint 0 for negative values generated by dynamics
2171!
2172      do k = kts, kte
2173        do i = its, ite
2174          qci(i,k,1) = max(qci(i,k,1),0.0)
2175          qrs(i,k,1) = max(qrs(i,k,1),0.0)
2176          qci(i,k,2) = max(qci(i,k,2),0.0)
2177          qrs(i,k,2) = max(qrs(i,k,2),0.0)
2178        enddo
2179      enddo
2180!
2181!----------------------------------------------------------------
2182!     latent heat for phase changes and heat capacity. neglect the
2183!     changes during microphysical process calculation
2184!     emanuel(1994)
2185!
2186      do k = kts, kte
2187        do i = its, ite
2188          cpm(i,k) = cpmcal(q(i,k))
2189          xl(i,k) = xlcal(t(i,k))
2190        enddo
2191      enddo
2192!
2193!----------------------------------------------------------------
2194!     compute the minor time steps.
2195!
2196      loops = max(nint(delt/dtcldcr),1)
2197      dtcld = delt/loops
2198      if(delt.le.dtcldcr) dtcld = delt
2199!
2200      do loop = 1,loops
2201!
2202!----------------------------------------------------------------
2203!     initialize the large scale variables
2204!
2205      do i = its, ite
2206        mstep(i) = 1
2207        flgcld(i) = .true.
2208      enddo
2209!
2210!     do k = kts, kte
2211!       do i = its, ite
2212!         denfac(i,k) = sqrt(den0/den(i,k))
2213!       enddo
2214!     enddo
2215      do k = kts, kte
2216        CALL VREC( tvec1(its), den(its,k), ite-its+1)
2217        do i = its, ite
2218          tvec1(i) = tvec1(i)*den0
2219        enddo
2220        CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
2221      enddo
2222!
2223! Inline expansion for fpvs
2224!         qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
2225!         qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
2226      hsub = xls
2227      hvap = xlv0
2228      cvap = cpv
2229      ttp=t0c+0.01
2230      dldt=cvap-cliq
2231      xa=-dldt/rv
2232      xb=xa+hvap/(rv*ttp)
2233      dldti=cvap-cice
2234      xai=-dldti/rv
2235      xbi=xai+hsub/(rv*ttp)
2236
2237! this is for compilers where the conditional inhibits vectorization
2238#ifdef WSM_NO_CONDITIONAL_IN_VECTOR
2239      do k = kts, kte
2240        do i = its, ite
2241          if(t(i,k).lt.ttp) then
2242            xal(i) = xai
2243            xbl(i) = xbi
2244          else
2245            xal(i) = xa
2246            xbl(i) = xb
2247          endif
2248        enddo
2249        do i = its, ite
2250          tr=ttp/t(i,k)
2251          logtr=log(tr)
2252          qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
2253          qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
2254          qs(i,k,1) = max(qs(i,k,1),qmin)
2255          rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
2256          qs(i,k,2)=psat*exp(logtr*(xal(i))+xbl(i)*(1.-tr))
2257          qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
2258          qs(i,k,2) = max(qs(i,k,2),qmin)
2259          rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
2260        enddo
2261      enddo
2262#else
2263      do k = kts, kte
2264        do i = its, ite
2265          tr=ttp/t(i,k)
2266          logtr=log(tr)
2267          qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
2268          qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
2269          qs(i,k,1) = max(qs(i,k,1),qmin)
2270          rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
2271          if(t(i,k).lt.ttp) then
2272            qs(i,k,2)=psat*exp(logtr*(xai)+xbi*(1.-tr))
2273          else
2274            qs(i,k,2)=psat*exp(logtr*(xa)+xb*(1.-tr))
2275          endif
2276          qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
2277          qs(i,k,2) = max(qs(i,k,2),qmin)
2278          rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
2279        enddo
2280      enddo
2281#endif
2282!
2283!----------------------------------------------------------------
2284!     initialize the variables for microphysical physics
2285!
2286!
2287      do k = kts, kte
2288        do i = its, ite
2289          prevp(i,k) = 0.
2290          psdep(i,k) = 0.
2291          praut(i,k) = 0.
2292          psaut(i,k) = 0.
2293          pracw(i,k) = 0.
2294          psaci(i,k) = 0.
2295          psacw(i,k) = 0.
2296          pigen(i,k) = 0.
2297          pidep(i,k) = 0.
2298          pcond(i,k) = 0.
2299          psmlt(i,k) = 0.
2300          psevp(i,k) = 0.
2301          falk(i,k,1) = 0.
2302          falk(i,k,2) = 0.
2303          fall(i,k,1) = 0.
2304          fall(i,k,2) = 0.
2305          fallc(i,k) = 0.
2306          falkc(i,k) = 0.
2307          xni(i,k) = 1.e3
2308        enddo
2309      enddo
2310!
2311!----------------------------------------------------------------
2312!     compute the fallout term:
2313!     first, vertical terminal velosity for minor loops
2314!
2315      do k = kts, kte
2316        do i = its, ite
2317          supcol = t0c-t(i,k)
2318!---------------------------------------------------------------
2319! n0s: Intercept parameter for snow [m-4] [HDC 6]
2320!---------------------------------------------------------------
2321          n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
2322          if(qrs(i,k,1).le.qcrmin)then
2323            rslope(i,k,1) = rslopermax
2324            rslopeb(i,k,1) = rsloperbmax
2325            rslope2(i,k,1) = rsloper2max
2326            rslope3(i,k,1) = rsloper3max
2327          else
2328            rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k))
2329            rslopeb(i,k,1) = exp(log(rslope(i,k,1))*(bvtr))
2330            rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
2331            rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
2332          endif
2333          if(qrs(i,k,2).le.qcrmin)then
2334            rslope(i,k,2) = rslopesmax
2335            rslopeb(i,k,2) = rslopesbmax
2336            rslope2(i,k,2) = rslopes2max
2337            rslope3(i,k,2) = rslopes3max
2338          else
2339            rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
2340            rslopeb(i,k,2) = exp(log(rslope(i,k,2))*(bvts))
2341            rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
2342            rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
2343          endif
2344!-------------------------------------------------------------
2345! Ni: ice crystal number concentraiton   [HDC 5c]
2346!-------------------------------------------------------------
2347!         xni(i,k) = min(max(5.38e7*(den(i,k)                                  &
2348!                   *max(qci(i,k,2),qmin))**0.75,1.e3),1.e6)
2349          temp = (den(i,k)*max(qci(i,k,2),qmin))
2350          temp = sqrt(sqrt(temp*temp*temp))
2351          xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
2352        enddo
2353      enddo
2354!
2355      mstepmax = 1
2356      numdt = 1
2357      do k = kte, kts, -1
2358        do i = its, ite
2359          work1(i,k,1) = pvtr*rslopeb(i,k,1)*denfac(i,k)/delz(i,k)
2360          work1(i,k,2) = pvts*rslopeb(i,k,2)*denfac(i,k)/delz(i,k)
2361          numdt(i) = max(nint(max(work1(i,k,1),work1(i,k,2))*dtcld+.5),1)
2362          if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
2363        enddo
2364      enddo
2365      do i = its, ite
2366        if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
2367        rmstep(i) = 1./mstep(i)
2368      enddo
2369!
2370      do n = 1, mstepmax
2371        k = kte
2372        do i = its, ite
2373          if(n.le.mstep(i)) then
2374!             falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
2375!             falk(i,k,2) = den(i,k)*qrs(i,k,2)*work1(i,k,2)/mstep(i)
2376              falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)*rmstep(i)
2377              falk(i,k,2) = den(i,k)*qrs(i,k,2)*work1(i,k,2)*rmstep(i)
2378              fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
2379              fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
2380!             qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcld/den(i,k),0.)
2381!             qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcld/den(i,k),0.)
2382              dtcldden = dtcld/den(i,k)
2383              qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcldden,0.)
2384              qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcldden,0.)
2385            endif
2386          enddo
2387        do k = kte-1, kts, -1
2388          do i = its, ite
2389            if(n.le.mstep(i)) then
2390              falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)*rmstep(i)
2391              falk(i,k,2) = den(i,k)*qrs(i,k,2)*work1(i,k,2)*rmstep(i)
2392              fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
2393              fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
2394              dtcldden = dtcld/den(i,k)
2395              rdelz = 1./delz(i,k)
2396              qrs(i,k,1) = max(qrs(i,k,1)-(falk(i,k,1)-falk(i,k+1,1)           &
2397                          *delz(i,k+1)*rdelz)*dtcldden,0.)
2398              qrs(i,k,2) = max(qrs(i,k,2)-(falk(i,k,2)-falk(i,k+1,2)           &
2399                          *delz(i,k+1)*rdelz)*dtcldden,0.)
2400            endif
2401          enddo
2402        enddo
2403        do k = kte, kts, -1
2404          do i = its, ite
2405            if(n.le.mstep(i)) then
2406              if(t(i,k).gt.t0c.and.qrs(i,k,2).gt.0.) then
2407!----------------------------------------------------------------
2408! psmlt: melting of snow [HL A33] [RH83 A25]
2409!       (T>T0: S->R)
2410!----------------------------------------------------------------
2411                xlf = xlf0
2412!               work2(i,k)= venfac(p(i,k),t(i,k),den(i,k))
2413                work2(i,k)= (exp(log(((1.496e-6*((t(i,k))*sqrt(t(i,k)))        &
2414                            /((t(i,k))+120.)/(den(i,k)))/(8.794e-5             &
2415                            *exp(log(t(i,k))*(1.81))/p(i,k))))                 &
2416                            *((.3333333)))/sqrt((1.496e-6*((t(i,k))            &
2417                            *sqrt(t(i,k)))/((t(i,k))+120.)/(den(i,k))))        &
2418                            *sqrt(sqrt(den0/(den(i,k)))))
2419                coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
2420!               psmlt(i,k) = xka(t(i,k),den(i,k))/xlf*(t0c-t(i,k))*pi/2.       &
2421!                           *n0sfac(i,k)*(precs1*rslope2(i,k,2)+precs2         &
2422!                           *work2(i,k)*coeres)
2423                psmlt(i,k) = (1.414e3*(1.496e-6*((t(i,k))*sqrt(t(i,k)))        &         
2424                            /((t(i,k))+120.)/(den(i,k)) )*(den(i,k)))          &
2425                            /xlf*(t0c-t(i,k))*pi/2.                            &
2426                            *n0sfac(i,k)*(precs1*rslope2(i,k,2)+precs2         &
2427                            *work2(i,k)*coeres)
2428                psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep(i),                &
2429                            -qrs(i,k,2)/mstep(i)),0.)
2430                qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
2431                qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k)
2432                t(i,k) = t(i,k) + xlf/cpm(i,k)*psmlt(i,k)
2433              endif
2434            endif
2435          enddo
2436        enddo
2437      enddo
2438!---------------------------------------------------------------
2439! Vice [ms-1] : fallout of ice crystal [HDC 5a]
2440!---------------------------------------------------------------
2441      mstepmax = 1
2442      mstep = 1
2443      numdt = 1
2444      do k = kte, kts, -1
2445        do i = its, ite
2446          if(qci(i,k,2).le.0.) then
2447            work2c(i,k) = 0.
2448          else
2449            xmi = den(i,k)*qci(i,k,2)/xni(i,k)
2450!           diameter  = min(dicon * sqrt(xmi),dimax)
2451            diameter  = max(min(dicon * sqrt(xmi),dimax), 1.e-25)
2452            work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31))
2453            work2c(i,k) = work1c(i,k)/delz(i,k)
2454          endif
2455          numdt(i) = max(nint(work2c(i,k)*dtcld+.5),1)
2456          if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
2457        enddo
2458      enddo
2459      do i = its, ite
2460        if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
2461      enddo
2462!
2463      do n = 1, mstepmax
2464        k = kte
2465        do i = its, ite
2466          if(n.le.mstep(i)) then
2467            falkc(i,k) = den(i,k)*qci(i,k,2)*work2c(i,k)/mstep(i)
2468            holdc = falkc(i,k)
2469            fallc(i,k) = fallc(i,k)+falkc(i,k)
2470            holdci = qci(i,k,2)
2471            qci(i,k,2) = max(qci(i,k,2)-falkc(i,k)*dtcld/den(i,k),0.)
2472          endif
2473        enddo
2474        do k = kte-1, kts, -1
2475          do i = its, ite
2476            if(n.le.mstep(i)) then
2477              falkc(i,k) = den(i,k)*qci(i,k,2)*work2c(i,k)/mstep(i)
2478              holdc = falkc(i,k)
2479              fallc(i,k) = fallc(i,k)+falkc(i,k)
2480              holdci = qci(i,k,2)
2481              qci(i,k,2) = max(qci(i,k,2)-(falkc(i,k)-falkc(i,k+1)             &
2482                          *delz(i,k+1)/delz(i,k))*dtcld/den(i,k),0.)
2483            endif
2484          enddo
2485        enddo
2486      enddo
2487!
2488!
2489!----------------------------------------------------------------
2490!      rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
2491!
2492      do i = its, ite
2493        fallsum = fall(i,1,1)+fall(i,1,2)+fallc(i,1)
2494        fallsum_qsi = fall(i,1,2)+fallc(i,1)
2495        rainncv(i) = 0.
2496        if(fallsum.gt.0.) then
2497          rainncv(i) = fallsum*delz(i,1)/denr*dtcld*1000.
2498          rain(i) = fallsum*delz(i,1)/denr*dtcld*1000. + rain(i)
2499        endif
2500        IF ( PRESENT (snowncv) .AND. PRESENT (snow)) THEN
2501        snowncv(i,lat) = 0.
2502        if(fallsum_qsi.gt.0.) then
2503          snowncv(i,lat) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000.
2504          snow(i,lat) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + snow(i,lat)
2505        endif
2506        ENDIF
2507        sr(i) = 0.
2508        if(fallsum.gt.0.)sr(i)=fallsum_qsi*delz(i,kts)/denr*dtcld*1000.        &   
2509        /(rainncv(i)+1.e-12)
2510      enddo
2511!
2512!---------------------------------------------------------------
2513! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
2514!       (T>T0: I->C)
2515!---------------------------------------------------------------
2516      do k = kts, kte
2517        do i = its, ite
2518          supcol = t0c-t(i,k)
2519          xlf = xls-xl(i,k)
2520          if(supcol.lt.0.) xlf = xlf0
2521          if(supcol.lt.0.and.qci(i,k,2).gt.0.) then
2522            qci(i,k,1) = qci(i,k,1) + qci(i,k,2)
2523            t(i,k) = t(i,k) - xlf/cpm(i,k)*qci(i,k,2)
2524            qci(i,k,2) = 0.
2525          endif
2526!---------------------------------------------------------------
2527! pihmf: homogeneous freezing of cloud water below -40c [HL A45]
2528!        (T<-40C: C->I)
2529!---------------------------------------------------------------
2530          if(supcol.gt.40..and.qci(i,k,1).gt.0.) then
2531            qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
2532            t(i,k) = t(i,k) + xlf/cpm(i,k)*qci(i,k,1)
2533            qci(i,k,1) = 0.
2534          endif
2535!---------------------------------------------------------------
2536! pihtf: heterogeneous freezing of cloud water [HL A44]
2537!        (T0>T>-40C: C->I)
2538!---------------------------------------------------------------
2539          if(supcol.gt.0..and.qci(i,k,1).gt.0.) then
2540            supcolt=min(supcol,50.)
2541!           pfrzdtc = min(pfrz1*(exp(pfrz2*supcol)-1.)                         &
2542!              *den(i,k)/denr/xncr*qci(i,k,1)**2*dtcld,qci(i,k,1))
2543            pfrzdtc = min(pfrz1*(exp(pfrz2*supcolt)-1.)                        &
2544            *den(i,k)/denr/xncr*qci(i,k,1)*qci(i,k,1)*dtcld,qci(i,k,1))
2545            qci(i,k,2) = qci(i,k,2) + pfrzdtc
2546            t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtc
2547            qci(i,k,1) = qci(i,k,1)-pfrzdtc
2548          endif
2549!---------------------------------------------------------------
2550! psfrz: freezing of rain water [HL A20] [LFO 45]
2551!        (T<T0, R->S)
2552!---------------------------------------------------------------
2553          if(supcol.gt.0..and.qrs(i,k,1).gt.0.) then
2554            supcolt=min(supcol,50.)
2555!           pfrzdtr = min(20.*pi**2*pfrz1*n0r*denr/den(i,k)                    &
2556!                 *(exp(pfrz2*supcol)-1.)*rslope(i,k,1)**7*dtcld,              &
2557!                 qrs(i,k,1))
2558            temp = rslope(i,k,1)
2559            temp = temp*temp*temp*temp*temp*temp*temp
2560            pfrzdtr = min(20.*(pi*pi)*pfrz1*n0r*denr/den(i,k)                  &
2561                  *(exp(pfrz2*supcolt)-1.)*temp*dtcld,                         &
2562                  qrs(i,k,1))
2563            qrs(i,k,2) = qrs(i,k,2) + pfrzdtr
2564            t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtr
2565            qrs(i,k,1) = qrs(i,k,1)-pfrzdtr
2566          endif
2567        enddo
2568      enddo
2569!
2570!----------------------------------------------------------------
2571!     rsloper: reverse of the slope parameter of the rain(m)
2572!     xka:    thermal conductivity of air(jm-1s-1k-1)
2573!     work1:  the thermodynamic term in the denominator associated with
2574!             heat conduction and vapor diffusion
2575!             (ry88, y93, h85)
2576!     work2: parameter associated with the ventilation effects(y93)
2577!
2578      do k = kts, kte
2579        do i = its, ite
2580          if(qrs(i,k,1).le.qcrmin)then
2581            rslope(i,k,1) = rslopermax
2582            rslopeb(i,k,1) = rsloperbmax
2583            rslope2(i,k,1) = rsloper2max
2584            rslope3(i,k,1) = rsloper3max
2585          else
2586!           rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k))
2587            rslope(i,k,1) = 1./(sqrt(sqrt(pidn0r/((qrs(i,k,1))*(den(i,k))))))
2588            rslopeb(i,k,1) = exp(log(rslope(i,k,1))*(bvtr))
2589            rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
2590            rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
2591          endif
2592          if(qrs(i,k,2).le.qcrmin)then
2593            rslope(i,k,2) = rslopesmax
2594            rslopeb(i,k,2) = rslopesbmax
2595            rslope2(i,k,2) = rslopes2max
2596            rslope3(i,k,2) = rslopes3max
2597          else
2598!            rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
2599            rslope(i,k,2) = 1./(sqrt(sqrt(pidn0s*(n0sfac(i,k))/((qrs(i,k,2))   & 
2600                          *(den(i,k))))))
2601            rslopeb(i,k,2) = exp(log(rslope(i,k,2))*(bvts))
2602            rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
2603            rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
2604          endif
2605        enddo
2606      enddo
2607!
2608      do k = kts, kte
2609        do i = its, ite
2610!         work1(i,k,1) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k,1))
2611          work1(i,k,1) = ((((den(i,k))*(xl(i,k))*(xl(i,k)))*((t(i,k))+120.)    &
2612                       *(den(i,k)))/(1.414e3*(1.496e-6*((t(i,k))*sqrt(t(i,k))))&
2613                       *(den(i,k))*(rv*(t(i,k))*(t(i,k)))))                    &
2614                      +  p(i,k)/((qs(i,k,1))*(8.794e-5*exp(log(t(i,k))*(1.81))))
2615!         work1(i,k,2) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k,2))
2616          work1(i,k,2) = ((((den(i,k))*(xls)*(xls))*((t(i,k))+120.)*(den(i,k)))&
2617                       /(1.414e3*(1.496e-6*((t(i,k))*sqrt(t(i,k))))*(den(i,k)) &
2618                       *(rv*(t(i,k))*(t(i,k))))                                &
2619                      + p(i,k)/(qs(i,k,2)*(8.794e-5*exp(log(t(i,k))*(1.81)))))
2620!         work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
2621          work2(i,k) = (exp(.3333333*log(((1.496e-6 * ((t(i,k))*sqrt(t(i,k)))) &
2622                     *p(i,k))/(((t(i,k))+120.)*den(i,k)*(8.794e-5              &
2623                     *exp(log(t(i,k))*(1.81))))))*sqrt(sqrt(den0/(den(i,k))))) &
2624                      /sqrt((1.496e-6*((t(i,k))*sqrt(t(i,k))))                 &
2625                      /(((t(i,k))+120.)*den(i,k)))
2626        enddo
2627      enddo
2628!
2629!===============================================================
2630!
2631! warm rain processes
2632!
2633! - follows the processes in RH83 and LFO except for autoconcersion
2634!
2635!===============================================================
2636!
2637      do k = kts, kte
2638        do i = its, ite
2639          supsat = max(q(i,k),qmin)-qs(i,k,1)
2640          satdt = supsat/dtcld
2641!---------------------------------------------------------------
2642! praut: auto conversion rate from cloud to rain [HDC 16]
2643!        (C->R)
2644!---------------------------------------------------------------
2645          if(qci(i,k,1).gt.qc0) then
2646            praut(i,k) = qck1*exp(log(qci(i,k,1))*((7./3.)))
2647            praut(i,k) = min(praut(i,k),qci(i,k,1)/dtcld)
2648          endif
2649!---------------------------------------------------------------
2650! pracw: accretion of cloud water by rain [HL A40] [LFO 51]
2651!        (C->R)
2652!---------------------------------------------------------------
2653          if(qrs(i,k,1).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
2654            pracw(i,k) = min(pacrr*rslope3(i,k,1)*rslopeb(i,k,1)               &
2655                         *qci(i,k,1)*denfac(i,k),qci(i,k,1)/dtcld)
2656          endif
2657!---------------------------------------------------------------
2658! prevp: evaporation/condensation rate of rain [HDC 14]
2659!        (V->R or R->V)
2660!---------------------------------------------------------------
2661          if(qrs(i,k,1).gt.0.) then
2662            coeres = rslope2(i,k,1)*sqrt(rslope(i,k,1)*rslopeb(i,k,1))
2663            prevp(i,k) = (rh(i,k,1)-1.)*(precr1*rslope2(i,k,1)                 &
2664                         +precr2*work2(i,k)*coeres)/work1(i,k,1)
2665            if(prevp(i,k).lt.0.) then
2666              prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
2667              prevp(i,k) = max(prevp(i,k),satdt/2)
2668            else
2669              prevp(i,k) = min(prevp(i,k),satdt/2)
2670            endif
2671          endif
2672        enddo
2673      enddo
2674!
2675!===============================================================
2676!
2677! cold rain processes
2678!
2679! - follows the revised ice microphysics processes in HDC
2680! - the processes same as in RH83 and RH84  and LFO behave
2681!   following ice crystal hapits defined in HDC, inclduing
2682!   intercept parameter for snow (n0s), ice crystal number
2683!   concentration (ni), ice nuclei number concentration
2684!   (n0i), ice diameter (d)
2685!
2686!===============================================================
2687!
2688      rdtcld = 1./dtcld
2689      do k = kts, kte
2690        do i = its, ite
2691          supcol = t0c-t(i,k)
2692          supsat = max(q(i,k),qmin)-qs(i,k,2)
2693          satdt = supsat/dtcld
2694          ifsat = 0
2695!-------------------------------------------------------------
2696! Ni: ice crystal number concentraiton   [HDC 5c]
2697!-------------------------------------------------------------
2698!         xni(i,k) = min(max(5.38e7*(den(i,k)                                  &
2699!                      *max(qci(i,k,2),qmin))**0.75,1.e3),1.e6)
2700          temp = (den(i,k)*max(qci(i,k,2),qmin))
2701          temp = sqrt(sqrt(temp*temp*temp))
2702          xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
2703          eacrs = exp(0.07*(-supcol))
2704!
2705          if(supcol.gt.0) then
2706            if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,2).gt.qmin) then
2707              xmi = den(i,k)*qci(i,k,2)/xni(i,k)
2708              diameter  = min(dicon * sqrt(xmi),dimax)
2709              vt2i = 1.49e4*diameter**1.31
2710              vt2s = pvts*rslopeb(i,k,2)*denfac(i,k)
2711!-------------------------------------------------------------
2712! psaci: Accretion of cloud ice by rain [HDC 10]
2713!        (T<T0: I->S)
2714!-------------------------------------------------------------
2715              acrfac = 2.*rslope3(i,k,2)+2.*diameter*rslope2(i,k,2)            &
2716                      +diameter**2*rslope(i,k,2)
2717              psaci(i,k) = pi*qci(i,k,2)*eacrs*n0s*n0sfac(i,k)                 &
2718                           *abs(vt2s-vt2i)*acrfac/4.
2719            endif
2720          endif
2721!-------------------------------------------------------------
2722! psacw: Accretion of cloud water by snow  [HL A7] [LFO 24]
2723!        (T<T0: C->S, and T>=T0: C->R)
2724!-------------------------------------------------------------
2725          if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
2726            psacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2)                  &
2727                           *rslopeb(i,k,2)*qci(i,k,1)*denfac(i,k)              &
2728!                          ,qci(i,k,1)/dtcld)
2729                           ,qci(i,k,1)*rdtcld)
2730          endif
2731          if(supcol .gt. 0) then
2732!-------------------------------------------------------------
2733! pidep: Deposition/Sublimation rate of ice [HDC 9]
2734!       (T<T0: V->I or I->V)
2735!-------------------------------------------------------------
2736            if(qci(i,k,2).gt.0.and.ifsat.ne.1) then
2737              xmi = den(i,k)*qci(i,k,2)/xni(i,k)
2738              diameter = dicon * sqrt(xmi)
2739              pidep(i,k) = 4.*diameter*xni(i,k)*(rh(i,k,2)-1.)/work1(i,k,2)
2740              supice = satdt-prevp(i,k)
2741              if(pidep(i,k).lt.0.) then
2742!               pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
2743!               pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
2744                pidep(i,k) = max(max(pidep(i,k),satdt*.5),supice)
2745                pidep(i,k) = max(pidep(i,k),-qci(i,k,2)*rdtcld)
2746              else
2747!               pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
2748                pidep(i,k) = min(min(pidep(i,k),satdt*.5),supice)
2749              endif
2750              if(abs(prevp(i,k)+pidep(i,k)).ge.abs(satdt)) ifsat = 1
2751            endif
2752!-------------------------------------------------------------
2753! psdep: deposition/sublimation rate of snow [HDC 14]
2754!        (V->S or S->V)
2755!-------------------------------------------------------------
2756            if(qrs(i,k,2).gt.0..and.ifsat.ne.1) then
2757              coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
2758              psdep(i,k) = (rh(i,k,2)-1.)*n0sfac(i,k)                          &
2759                           *(precs1*rslope2(i,k,2)+precs2                      &
2760                           *work2(i,k)*coeres)/work1(i,k,2)
2761              supice = satdt-prevp(i,k)-pidep(i,k)
2762              if(psdep(i,k).lt.0.) then
2763!               psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
2764!               psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
2765                psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)*rdtcld)
2766                psdep(i,k) = max(max(psdep(i,k),satdt*.5),supice)
2767              else
2768!             psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
2769                psdep(i,k) = min(min(psdep(i,k),satdt*.5),supice)
2770              endif
2771              if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt))          &
2772                ifsat = 1
2773            endif
2774!-------------------------------------------------------------
2775! pigen: generation(nucleation) of ice from vapor [HL A50] [HDC 7-8]
2776!       (T<T0: V->I)
2777!-------------------------------------------------------------
2778            if(supsat.gt.0.and.ifsat.ne.1) then
2779              supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
2780              xni0 = 1.e3*exp(0.1*supcol)
2781              roqi0 = 4.92e-11*exp(log(xni0)*(1.33))
2782              pigen(i,k) = max(0.,(roqi0/den(i,k)-max(qci(i,k,2),0.))          &
2783!                        /dtcld)
2784                         *rdtcld)
2785              pigen(i,k) = min(min(pigen(i,k),satdt),supice)
2786            endif
2787!
2788!-------------------------------------------------------------
2789! psaut: conversion(aggregation) of ice to snow [HDC 12]
2790!       (T<T0: I->S)
2791!-------------------------------------------------------------
2792            if(qci(i,k,2).gt.0.) then
2793              qimax = roqimax/den(i,k)
2794!             psaut(i,k) = max(0.,(qci(i,k,2)-qimax)/dtcld)
2795              psaut(i,k) = max(0.,(qci(i,k,2)-qimax)*rdtcld)
2796            endif
2797          endif
2798!-------------------------------------------------------------
2799! psevp: Evaporation of melting snow [HL A35] [RH83 A27]
2800!       (T>T0: S->V)
2801!-------------------------------------------------------------
2802          if(supcol.lt.0.) then
2803            if(qrs(i,k,2).gt.0..and.rh(i,k,1).lt.1.)                           &
2804              psevp(i,k) = psdep(i,k)*work1(i,k,2)/work1(i,k,1)
2805!              psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)/dtcld),0.)
2806              psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)*rdtcld),0.)
2807          endif
2808        enddo
2809      enddo
2810!
2811!
2812!----------------------------------------------------------------
2813!     check mass conservation of generation terms and feedback to the
2814!     large scale
2815!
2816      do k = kts, kte
2817        do i = its, ite
2818          if(t(i,k).le.t0c) then
2819!
2820!     cloud water
2821!
2822            value = max(qmin,qci(i,k,1))
2823            source = (praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
2824            if (source.gt.value) then
2825              factor = value/source
2826              praut(i,k) = praut(i,k)*factor
2827              pracw(i,k) = pracw(i,k)*factor
2828              psacw(i,k) = psacw(i,k)*factor
2829            endif
2830!
2831!     cloud ice
2832!
2833            value = max(qmin,qci(i,k,2))
2834            source = (psaut(i,k)+psaci(i,k)-pigen(i,k)-pidep(i,k))*dtcld
2835            if (source.gt.value) then
2836              factor = value/source
2837              psaut(i,k) = psaut(i,k)*factor
2838              psaci(i,k) = psaci(i,k)*factor
2839              pigen(i,k) = pigen(i,k)*factor
2840              pidep(i,k) = pidep(i,k)*factor
2841            endif
2842!
2843!     rain
2844!
2845!
2846            value = max(qmin,qrs(i,k,1))
2847            source = (-praut(i,k)-pracw(i,k)-prevp(i,k))*dtcld
2848            if (source.gt.value) then
2849              factor = value/source
2850              praut(i,k) = praut(i,k)*factor
2851              pracw(i,k) = pracw(i,k)*factor
2852              prevp(i,k) = prevp(i,k)*factor
2853            endif
2854!
2855!    snow
2856!
2857            value = max(qmin,qrs(i,k,2))
2858            source = (-psdep(i,k)-psaut(i,k)-psaci(i,k)-psacw(i,k))*dtcld 
2859            if (source.gt.value) then
2860              factor = value/source
2861              psdep(i,k) = psdep(i,k)*factor
2862              psaut(i,k) = psaut(i,k)*factor
2863              psaci(i,k) = psaci(i,k)*factor
2864              psacw(i,k) = psacw(i,k)*factor
2865            endif
2866!
2867            work2(i,k)=-(prevp(i,k)+psdep(i,k)+pigen(i,k)+pidep(i,k))
2868!     update
2869            q(i,k) = q(i,k)+work2(i,k)*dtcld
2870            qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k)                 &
2871                        +psacw(i,k))*dtcld,0.)
2872            qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k)                 &
2873                        +prevp(i,k))*dtcld,0.)
2874            qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+psaci(i,k)                 &
2875                        -pigen(i,k)-pidep(i,k))*dtcld,0.)
2876            qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k)                 &
2877                        +psaci(i,k)+psacw(i,k))*dtcld,0.)
2878            xlf = xls-xl(i,k)
2879            xlwork2 = -xls*(psdep(i,k)+pidep(i,k)+pigen(i,k))                  &
2880                      -xl(i,k)*prevp(i,k)-xlf*psacw(i,k)
2881            t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
2882          else
2883!
2884!     cloud water
2885!
2886            value = max(qmin,qci(i,k,1))
2887            source=(praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
2888            if (source.gt.value) then
2889              factor = value/source
2890              praut(i,k) = praut(i,k)*factor
2891              pracw(i,k) = pracw(i,k)*factor
2892              psacw(i,k) = psacw(i,k)*factor
2893            endif
2894!
2895!     rain
2896!
2897            value = max(qmin,qrs(i,k,1))
2898            source = (-praut(i,k)-pracw(i,k)-prevp(i,k)-psacw(i,k))*dtcld
2899            if (source.gt.value) then
2900              factor = value/source
2901              praut(i,k) = praut(i,k)*factor
2902              pracw(i,k) = pracw(i,k)*factor
2903              prevp(i,k) = prevp(i,k)*factor
2904              psacw(i,k) = psacw(i,k)*factor
2905            endif 
2906!
2907!     snow
2908!
2909            value = max(qcrmin,qrs(i,k,2))
2910            source=(-psevp(i,k))*dtcld
2911            if (source.gt.value) then
2912              factor = value/source
2913              psevp(i,k) = psevp(i,k)*factor
2914            endif
2915            work2(i,k)=-(prevp(i,k)+psevp(i,k))
2916!     update
2917            q(i,k) = q(i,k)+work2(i,k)*dtcld
2918            qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k)                 &
2919                        +psacw(i,k))*dtcld,0.)
2920            qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k)                 &
2921                        +prevp(i,k) +psacw(i,k))*dtcld,0.)
2922            qrs(i,k,2) = max(qrs(i,k,2)+psevp(i,k)*dtcld,0.)
2923            xlf = xls-xl(i,k)
2924            xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k))
2925            t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
2926          endif
2927        enddo
2928      enddo
2929!
2930! Inline expansion for fpvs
2931!         qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
2932!         qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
2933      hsub = xls
2934      hvap = xlv0
2935      cvap = cpv
2936      ttp=t0c+0.01
2937      dldt=cvap-cliq
2938      xa=-dldt/rv
2939      xb=xa+hvap/(rv*ttp)
2940      dldti=cvap-cice
2941      xai=-dldti/rv
2942      xbi=xai+hsub/(rv*ttp)
2943      do k = kts, kte
2944      do i = its, ite
2945          tr=ttp/t(i,k)
2946          logtr = log(tr)
2947          qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
2948          qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
2949          qs(i,k,1) = max(qs(i,k,1),qmin)
2950        enddo
2951      enddo
2952!
2953!----------------------------------------------------------------
2954!  pcond: condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
2955!     if there exists additional water vapor condensated/if
2956!     evaporation of cloud water is not enough to remove subsaturation
2957!
2958      do k = kts, kte
2959        do i = its, ite
2960!         work1(i,k,1) = conden(t(i,k),q(i,k),qs(i,k,1),xl(i,k),cpm(i,k))
2961          work1(i,k,1) = ((max(q(i,k),qmin)-(qs(i,k,1)))/(1.+(xl(i,k))         &   
2962                        *(xl(i,k))/(rv*(cpm(i,k)))*(qs(i,k,1))                 &
2963                        /((t(i,k))*(t(i,k)))))
2964          work2(i,k) = qci(i,k,1)+work1(i,k,1)
2965          pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k),0.)/dtcld)
2966          if(qci(i,k,1).gt.0..and.work1(i,k,1).lt.0.)                          &
2967            pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
2968          q(i,k) = q(i,k)-pcond(i,k)*dtcld
2969          qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
2970          t(i,k) = t(i,k)+pcond(i,k)*xl(i,k)/cpm(i,k)*dtcld
2971        enddo
2972      enddo
2973!
2974!
2975!----------------------------------------------------------------
2976!     padding for small values
2977!
2978      do k = kts, kte
2979        do i = its, ite
2980          if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
2981          if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
2982        enddo
2983      enddo
2984      enddo                  ! big loops
2985  END SUBROUTINE wsm52d
2986
2987#endif
2988
2989! ...................................................................
2990      REAL FUNCTION rgmma(x)
2991!-------------------------------------------------------------------
2992  IMPLICIT NONE
2993!-------------------------------------------------------------------
2994!     rgmma function:  use infinite product form
2995      REAL :: euler
2996      PARAMETER (euler=0.577215664901532)
2997      REAL :: x, y
2998      INTEGER :: i
2999      if(x.eq.1.)then
3000        rgmma=0.
3001          else
3002        rgmma=x*exp(euler*x)
3003        do i=1,10000
3004          y=float(i)
3005          rgmma=rgmma*(1.000+x/y)*exp(-x/y)
3006        enddo
3007        rgmma=1./rgmma
3008      endif
3009      END FUNCTION rgmma
3010!
3011!--------------------------------------------------------------------------
3012      REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
3013!--------------------------------------------------------------------------
3014      IMPLICIT NONE
3015!--------------------------------------------------------------------------
3016      REAL t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti,         &
3017           xai,xbi,ttp,tr
3018      INTEGER ice
3019! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3020      ttp=t0c+0.01
3021      dldt=cvap-cliq
3022      xa=-dldt/rv
3023      xb=xa+hvap/(rv*ttp)
3024      dldti=cvap-cice
3025      xai=-dldti/rv
3026      xbi=xai+hsub/(rv*ttp)
3027      tr=ttp/t
3028      if(t.lt.ttp.and.ice.eq.1) then
3029        fpvs=psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
3030      else
3031        fpvs=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
3032      endif
3033! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3034      END FUNCTION fpvs
3035!-------------------------------------------------------------------
3036  SUBROUTINE wsm5init(den0,denr,dens,cl,cpv,allowed_to_read)
3037!-------------------------------------------------------------------
3038  IMPLICIT NONE
3039!-------------------------------------------------------------------
3040!.... constants which may not be tunable
3041   REAL, INTENT(IN) :: den0,denr,dens,cl,cpv
3042   LOGICAL, INTENT(IN) :: allowed_to_read
3043   REAL :: pi
3044!
3045   pi = 4.*atan(1.)
3046   xlv1 = cl-cpv
3047!
3048   qc0  = 4./3.*pi*denr*r0**3*xncr/den0  ! 0.419e-3 -- .61e-3
3049   qck1 = .104*9.8*peaut/(xncr*denr)**(1./3.)/xmyu*den0**(4./3.) ! 7.03
3050!
3051   bvtr1 = 1.+bvtr
3052   bvtr2 = 2.5+.5*bvtr
3053   bvtr3 = 3.+bvtr
3054   bvtr4 = 4.+bvtr
3055   g1pbr = rgmma(bvtr1)
3056   g3pbr = rgmma(bvtr3)
3057   g4pbr = rgmma(bvtr4)            ! 17.837825
3058   g5pbro2 = rgmma(bvtr2)          ! 1.8273
3059   pvtr = avtr*g4pbr/6.
3060   eacrr = 1.0
3061   pacrr = pi*n0r*avtr*g3pbr*.25*eacrr
3062   precr1 = 2.*pi*n0r*.78
3063   precr2 = 2.*pi*n0r*.31*avtr**.5*g5pbro2
3064   xmmax = (dimax/dicon)**2
3065   roqimax = 2.08e22*dimax**8
3066!
3067   bvts1 = 1.+bvts
3068   bvts2 = 2.5+.5*bvts
3069   bvts3 = 3.+bvts
3070   bvts4 = 4.+bvts
3071   g1pbs = rgmma(bvts1)    !.8875
3072   g3pbs = rgmma(bvts3)
3073   g4pbs = rgmma(bvts4)    ! 12.0786
3074   g5pbso2 = rgmma(bvts2)
3075   pvts = avts*g4pbs/6.
3076   pacrs = pi*n0s*avts*g3pbs*.25
3077   precs1 = 4.*n0s*.65
3078   precs2 = 4.*n0s*.44*avts**.5*g5pbso2
3079   pidn0r =  pi*denr*n0r
3080   pidn0s =  pi*dens*n0s
3081   pacrc = pi*n0s*avts*g3pbs*.25*eacrc
3082!
3083   rslopermax = 1./lamdarmax
3084   rslopesmax = 1./lamdasmax
3085   rsloperbmax = rslopermax ** bvtr
3086   rslopesbmax = rslopesmax ** bvts
3087   rsloper2max = rslopermax * rslopermax
3088   rslopes2max = rslopesmax * rslopesmax
3089   rsloper3max = rsloper2max * rslopermax
3090   rslopes3max = rslopes2max * rslopesmax
3091!
3092  END SUBROUTINE wsm5init
3093END MODULE module_mp_wsm5
Note: See TracBrowser for help on using the repository browser.