<p><b>laura@ucar.edu</b> 2011-04-08 16:33:14 -0600 (Fri, 08 Apr 2011)</p><p>added sourcecode for wsm6 cloud microphysics scheme<br>
</p><hr noshade><pre><font color="gray">Added: branches/atmos_physics/src/core_physics/physics_wrf/libmassv.F
===================================================================
--- branches/atmos_physics/src/core_physics/physics_wrf/libmassv.F (rev 0)
+++ branches/atmos_physics/src/core_physics/physics_wrf/libmassv.F 2011-04-08 22:33:14 UTC (rev 789)
@@ -0,0 +1,390 @@
+! IBM libmassv compatibility library
+!
+
+#ifndef NATIVE_MASSV
+ subroutine vdiv(z,x,y,n)
+ real*8 x(*),y(*),z(*)
+ do 10 j=1,n
+ z(j)=x(j)/y(j)
+ 10 continue
+ return
+ end
+
+ subroutine vsdiv(z,x,y,n)
+ real*4 x(*),y(*),z(*)
+ do 10 j=1,n
+ z(j)=x(j)/y(j)
+ 10 continue
+ return
+ end
+
+ subroutine vexp(y,x,n)
+ real*8 x(*),y(*)
+ do 10 j=1,n
+ y(j)=exp(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vsexp(y,x,n)
+ real*4 x(*),y(*)
+ do 10 j=1,n
+ y(j)=exp(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vlog(y,x,n)
+ real*8 x(*),y(*)
+ do 10 j=1,n
+ y(j)=log(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vslog(y,x,n)
+ real*4 x(*),y(*)
+ do 10 j=1,n
+ y(j)=log(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vrec(y,x,n)
+ real*8 x(*),y(*)
+ do 10 j=1,n
+ y(j)=1.d0/x(j)
+ 10 continue
+ return
+ end
+
+ subroutine vsrec(y,x,n)
+ real*4 x(*),y(*)
+ do 10 j=1,n
+ y(j)=1.e0/x(j)
+ 10 continue
+ return
+ end
+
+ subroutine vrsqrt(y,x,n)
+ real*8 x(*),y(*)
+ do 10 j=1,n
+ y(j)=1.d0/sqrt(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vsrsqrt(y,x,n)
+ real*4 x(*),y(*)
+ do 10 j=1,n
+ y(j)=1.e0/sqrt(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vsincos(x,y,z,n)
+ real*8 x(*),y(*),z(*)
+ do 10 j=1,n
+ x(j)=sin(z(j))
+ y(j)=cos(z(j))
+ 10 continue
+ return
+ end
+
+ subroutine vssincos(x,y,z,n)
+ real*4 x(*),y(*),z(*)
+ do 10 j=1,n
+ x(j)=sin(z(j))
+ y(j)=cos(z(j))
+ 10 continue
+ return
+ end
+
+ subroutine vsqrt(y,x,n)
+ real*8 x(*),y(*)
+ do 10 j=1,n
+ y(j)=sqrt(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vssqrt(y,x,n)
+ real*4 x(*),y(*)
+ do 10 j=1,n
+ y(j)=sqrt(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vtan(y,x,n)
+ real*8 x(*),y(*)
+ do 10 j=1,n
+ y(j)=tan(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vstan(y,x,n)
+ real*4 x(*),y(*)
+ do 10 j=1,n
+ y(j)=tan(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vatan2(z,y,x,n)
+ real*8 x(*),y(*),z(*)
+ do 10 j=1,n
+ z(j)=atan2(y(j),x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vsatan2(z,y,x,n)
+ real*4 x(*),y(*),z(*)
+ do 10 j=1,n
+ z(j)=atan2(y(j),x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vasin(y,x,n)
+ real*8 x(*),y(*)
+ do 10 j=1,n
+ y(j)=asin(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vsin(y,x,n)
+ real*8 x(*),y(*)
+ do 10 j=1,n
+ y(j)=sin(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vssin(y,x,n)
+ real*4 x(*),y(*)
+ do 10 j=1,n
+ y(j)=sin(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vacos(y,x,n)
+ real*8 x(*),y(*)
+ do 10 j=1,n
+ y(j)=acos(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vcos(y,x,n)
+ real*8 x(*),y(*)
+ do 10 j=1,n
+ y(j)=cos(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vscos(y,x,n)
+ real*4 x(*),y(*)
+ do 10 j=1,n
+ y(j)=cos(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vcosisin(y,x,n)
+ complex*16 y(*)
+ real*8 x(*)
+ do 10 j=1,n
+ y(j)=dcmplx(cos(x(j)),sin(x(j)))
+ 10 continue
+ return
+ end
+
+ subroutine vscosisin(y,x,n)
+ complex*8 y(*)
+ real*4 x(*)
+ do 10 j=1,n
+ y(j)= cmplx(cos(x(j)),sin(x(j)))
+ 10 continue
+ return
+ end
+
+ subroutine vdint(y,x,n)
+ real*8 x(*),y(*)
+ do 10 j=1,n
+! y(j)=dint(x(j))
+ y(j)=int(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vdnint(y,x,n)
+ real*8 x(*),y(*)
+ do 10 j=1,n
+! y(j)=dnint(x(j))
+ y(j)=nint(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vlog10(y,x,n)
+ real*8 x(*),y(*)
+ do 10 j=1,n
+ y(j)=log10(x(j))
+ 10 continue
+ return
+ end
+
+! subroutine vlog1p(y,x,n)
+! real*8 x(*),y(*)
+! interface
+! real*8 function log1p(%val(x))
+! real*8 x
+! end function log1p
+! end interface
+! do 10 j=1,n
+! y(j)=log1p(x(j))
+! 10 continue
+! return
+! end
+
+ subroutine vcosh(y,x,n)
+ real*8 x(*),y(*)
+ do 10 j=1,n
+ y(j)=cosh(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vsinh(y,x,n)
+ real*8 x(*),y(*)
+ do 10 j=1,n
+ y(j)=sinh(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vtanh(y,x,n)
+ real*8 x(*),y(*)
+ do 10 j=1,n
+ y(j)=tanh(x(j))
+ 10 continue
+ return
+ end
+
+! subroutine vexpm1(y,x,n)
+! real*8 x(*),y(*)
+! interface
+! real*8 function expm1(%val(x))
+! real*8 x
+! end function expm1
+! end interface
+! do 10 j=1,n
+! y(j)=expm1(x(j))
+! 10 continue
+! return
+! end
+
+
+ subroutine vsasin(y,x,n)
+ real*4 x(*),y(*)
+ do 10 j=1,n
+ y(j)=asin(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vsacos(y,x,n)
+ real*4 x(*),y(*)
+ do 10 j=1,n
+#if defined (G95)
+! no reason why g95 should fail - oh well, we don't use this routine anyways
+ y(j)=asin( sqrt(1-x(j)*x(j)) )
+#else
+ y(j)=acos(x(j))
+#endif
+ 10 continue
+ return
+ end
+
+ subroutine vscosh(y,x,n)
+ real*4 x(*),y(*)
+ do 10 j=1,n
+ y(j)=cosh(x(j))
+ 10 continue
+ return
+ end
+
+! subroutine vsexpm1(y,x,n)
+! real*4 x(*),y(*)
+! interface
+! real*8 function expm1(%val(x))
+! real*8 x
+! end function expm1
+! end interface
+! do 10 j=1,n
+! y(j)=expm1(real(x(j),8))
+! 10 continue
+! return
+! end
+
+ subroutine vslog10(y,x,n)
+ real*4 x(*),y(*)
+ do 10 j=1,n
+ y(j)=log10(x(j))
+ 10 continue
+ return
+ end
+
+! subroutine vslog1p(y,x,n)
+! real*4 x(*),y(*)
+! interface
+! real*8 function log1p(%val(x))
+! real*8 x
+! end function log1p
+! end interface
+! do 10 j=1,n
+! y(j)=log1p(real(x(j),8))
+! 10 continue
+! return
+! end
+
+
+ subroutine vssinh(y,x,n)
+ real*4 x(*),y(*)
+ do 10 j=1,n
+ y(j)=sinh(x(j))
+ 10 continue
+ return
+ end
+
+ subroutine vstanh(y,x,n)
+ real*4 x(*),y(*)
+ do 10 j=1,n
+ y(j)=tanh(x(j))
+ 10 continue
+ return
+ end
+#endif
+
+ subroutine vspow(z,y,x,n)
+ real*4 x(*),y(*),z(*)
+ do 10 j=1,n
+ z(j)=y(j)**x(j)
+ 10 continue
+ return
+ end
+
+ subroutine vpow(z,y,x,n)
+ real*8 x(*),y(*),z(*)
+ do 10 j=1,n
+ z(j)=y(j)**x(j)
+ 10 continue
+ return
+ end
+
Added: branches/atmos_physics/src/core_physics/physics_wrf/module_mp_wsm6.F
===================================================================
--- branches/atmos_physics/src/core_physics/physics_wrf/module_mp_wsm6.F (rev 0)
+++ branches/atmos_physics/src/core_physics/physics_wrf/module_mp_wsm6.F 2011-04-08 22:33:14 UTC (rev 789)
@@ -0,0 +1,2218 @@
+#if ( RWORDSIZE == 4 )
+# define VREC vsrec
+# define VSQRT vssqrt
+#else
+# define VREC vrec
+# define VSQRT vsqrt
+#endif
+
+MODULE module_mp_wsm6
+!
+!
+ REAL, PARAMETER, PRIVATE :: dtcldcr = 120. ! maximum time step for minor loops
+ REAL, PARAMETER, PRIVATE :: n0r = 8.e6 ! intercept parameter rain
+ REAL, PARAMETER, PRIVATE :: n0g = 4.e6 ! intercept parameter graupel
+ REAL, PARAMETER, PRIVATE :: avtr = 841.9 ! a constant for terminal velocity of rain
+ REAL, PARAMETER, PRIVATE :: bvtr = 0.8 ! a constant for terminal velocity of rain
+ REAL, PARAMETER, PRIVATE :: r0 = .8e-5 ! 8 microm in contrast to 10 micro m
+ REAL, PARAMETER, PRIVATE :: peaut = .55 ! collection efficiency
+ REAL, PARAMETER, PRIVATE :: xncr = 3.e8 ! maritime cloud in contrast to 3.e8 in tc80
+ REAL, PARAMETER, PRIVATE :: xmyu = 1.718e-5 ! the dynamic viscosity kgm-1s-1
+ REAL, PARAMETER, PRIVATE :: avts = 11.72 ! a constant for terminal velocity of snow
+ REAL, PARAMETER, PRIVATE :: bvts = .41 ! a constant for terminal velocity of snow
+ REAL, PARAMETER, PRIVATE :: avtg = 330. ! a constant for terminal velocity of graupel
+ REAL, PARAMETER, PRIVATE :: bvtg = 0.8 ! a constant for terminal velocity of graupel
+ REAL, PARAMETER, PRIVATE :: deng = 500. ! density of graupel
+ REAL, PARAMETER, PRIVATE :: n0smax = 1.e11 ! maximum n0s (t=-90C unlimited)
+ REAL, PARAMETER, PRIVATE :: lamdarmax = 8.e4 ! limited maximum value for slope parameter of rain
+ REAL, PARAMETER, PRIVATE :: lamdasmax = 1.e5 ! limited maximum value for slope parameter of snow
+ REAL, PARAMETER, PRIVATE :: lamdagmax = 6.e4 ! limited maximum value for slope parameter of graupel
+ REAL, PARAMETER, PRIVATE :: dicon = 11.9 ! constant for the cloud-ice diamter
+ REAL, PARAMETER, PRIVATE :: dimax = 500.e-6 ! limited maximum value for the cloud-ice diamter
+ REAL, PARAMETER, PRIVATE :: n0s = 2.e6 ! temperature dependent intercept parameter snow
+ REAL, PARAMETER, PRIVATE :: alpha = .12 ! .122 exponen factor for n0s
+ REAL, PARAMETER, PRIVATE :: pfrz1 = 100. ! constant in Biggs freezing
+ REAL, PARAMETER, PRIVATE :: pfrz2 = 0.66 ! constant in Biggs freezing
+ REAL, PARAMETER, PRIVATE :: qcrmin = 1.e-9 ! minimun values for qr, qs, and qg
+ REAL, PARAMETER, PRIVATE :: eacrc = 1.0 ! Snow/cloud-water collection efficiency
+ REAL, PARAMETER, PRIVATE :: dens = 100.0 ! Density of snow
+ REAL, PARAMETER, PRIVATE :: qs0 = 6.e-4 ! threshold amount for aggretion to occur
+ REAL, SAVE :: &
+ qc0, qck1,bvtr1,bvtr2,bvtr3,bvtr4,g1pbr, &
+ g3pbr,g4pbr,g5pbro2,pvtr,eacrr,pacrr, &
+ bvtr6,g6pbr, &
+ precr1,precr2,roqimax,bvts1, &
+ bvts2,bvts3,bvts4,g1pbs,g3pbs,g4pbs, &
+ g5pbso2,pvts,pacrs,precs1,precs2,pidn0r, &
+ pidn0s,xlv1,pacrc,pi, &
+ bvtg1,bvtg2,bvtg3,bvtg4,g1pbg, &
+ g3pbg,g4pbg,g5pbgo2,pvtg,pacrg, &
+ precg1,precg2,pidn0g, &
+ rslopermax,rslopesmax,rslopegmax, &
+ rsloperbmax,rslopesbmax,rslopegbmax, &
+ rsloper2max,rslopes2max,rslopeg2max, &
+ rsloper3max,rslopes3max,rslopeg3max
+CONTAINS
+!===================================================================
+!
+ SUBROUTINE wsm6(th, q, qc, qr, qi, qs, qg &
+ ,den, pii, p, delz &
+ ,delt,g, cpd, cpv, rd, rv, t0c &
+ ,ep1, ep2, qmin &
+ ,XLS, XLV0, XLF0, den0, denr &
+ ,cliq,cice,psat &
+ ,rain, rainncv &
+ ,snow, snowncv &
+ ,sr &
+ ,graupel, graupelncv &
+ ,ids,ide, jds,jde, kds,kde &
+ ,ims,ime, jms,jme, kms,kme &
+ ,its,ite, jts,jte, kts,kte &
+ )
+!-------------------------------------------------------------------
+ IMPLICIT NONE
+!-------------------------------------------------------------------
+ INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
+ ims,ime, jms,jme, kms,kme , &
+ its,ite, jts,jte, kts,kte
+ REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
+ INTENT(INOUT) :: &
+ th, &
+ q, &
+ qc, &
+ qi, &
+ qr, &
+ qs, &
+ qg
+ REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
+ INTENT(IN ) :: &
+ den, &
+ pii, &
+ p, &
+ delz
+ REAL, INTENT(IN ) :: delt, &
+ g, &
+ rd, &
+ rv, &
+ t0c, &
+ den0, &
+ cpd, &
+ cpv, &
+ ep1, &
+ ep2, &
+ qmin, &
+ XLS, &
+ XLV0, &
+ XLF0, &
+ cliq, &
+ cice, &
+ psat, &
+ denr
+ REAL, DIMENSION( ims:ime , jms:jme ), &
+ INTENT(INOUT) :: rain, &
+ rainncv, &
+ sr
+ REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
+ INTENT(INOUT) :: snow, &
+ snowncv
+ REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
+ INTENT(INOUT) :: graupel, &
+ graupelncv
+! LOCAL VAR
+ REAL, DIMENSION( its:ite , kts:kte ) :: t
+ REAL, DIMENSION( its:ite , kts:kte, 2 ) :: qci
+ REAL, DIMENSION( its:ite , kts:kte, 3 ) :: qrs
+ INTEGER :: i,j,k
+!-------------------------------------------------------------------
+ DO j=jts,jte
+ DO k=kts,kte
+ DO i=its,ite
+ t(i,k)=th(i,k,j)*pii(i,k,j)
+ qci(i,k,1) = qc(i,k,j)
+ qci(i,k,2) = qi(i,k,j)
+ qrs(i,k,1) = qr(i,k,j)
+ qrs(i,k,2) = qs(i,k,j)
+ qrs(i,k,3) = qg(i,k,j)
+ ENDDO
+ ENDDO
+ ! Sending array starting locations of optional variables may cause
+ ! troubles, so we explicitly change the call.
+ CALL wsm62D(t, q(ims,kms,j), qci, qrs &
+ ,den(ims,kms,j) &
+ ,p(ims,kms,j), delz(ims,kms,j) &
+ ,delt,g, cpd, cpv, rd, rv, t0c &
+ ,ep1, ep2, qmin &
+ ,XLS, XLV0, XLF0, den0, denr &
+ ,cliq,cice,psat &
+ ,j &
+ ,rain(ims,j),rainncv(ims,j) &
+ ,sr(ims,j) &
+ ,ids,ide, jds,jde, kds,kde &
+ ,ims,ime, jms,jme, kms,kme &
+ ,its,ite, jts,jte, kts,kte &
+ ,snow,snowncv &
+ ,graupel,graupelncv &
+ )
+ DO K=kts,kte
+ DO I=its,ite
+ th(i,k,j)=t(i,k)/pii(i,k,j)
+ qc(i,k,j) = qci(i,k,1)
+ qi(i,k,j) = qci(i,k,2)
+ qr(i,k,j) = qrs(i,k,1)
+ qs(i,k,j) = qrs(i,k,2)
+ qg(i,k,j) = qrs(i,k,3)
+ ENDDO
+ ENDDO
+ ENDDO
+ END SUBROUTINE wsm6
+!===================================================================
+!
+ SUBROUTINE wsm62D(t, q &
+ ,qci, qrs, den, p, delz &
+ ,delt,g, cpd, cpv, rd, rv, t0c &
+ ,ep1, ep2, qmin &
+ ,XLS, XLV0, XLF0, den0, denr &
+ ,cliq,cice,psat &
+ ,lat &
+ ,rain,rainncv &
+ ,sr &
+ ,ids,ide, jds,jde, kds,kde &
+ ,ims,ime, jms,jme, kms,kme &
+ ,its,ite, jts,jte, kts,kte &
+ ,snow,snowncv &
+ ,graupel,graupelncv &
+ )
+!-------------------------------------------------------------------
+ IMPLICIT NONE
+!-------------------------------------------------------------------
+!
+! This code is a 6-class GRAUPEL phase microphyiscs scheme (WSM6) of the
+! Single-Moment MicroPhyiscs (WSMMP). The WSMMP assumes that ice nuclei
+! number concentration is a function of temperature, and seperate assumption
+! is developed, in which ice crystal number concentration is a function
+! of ice amount. A theoretical background of the ice-microphysics and related
+! processes in the WSMMPs are described in Hong et al. (2004).
+! All production terms in the WSM6 scheme are described in Hong and Lim (2006).
+! All units are in m.k.s. and source/sink terms in kgkg-1s-1.
+!
+! WSM6 cloud scheme
+!
+! Coded by Song-You Hong and Jeong-Ock Jade Lim (Yonsei Univ.)
+! Summer 2003
+!
+! Implemented by Song-You Hong (Yonsei Univ.) and Jimy Dudhia (NCAR)
+! Summer 2004
+!
+! History : semi-lagrangian scheme sedimentation(JH), and clean up
+! Hong, August 2009
+!
+! Reference) Hong, Dudhia, Chen (HDC, 2004) Mon. Wea. Rev.
+! Hong and Lim (HL, 2006) J. Korean Meteor. Soc.
+! Dudhia, Hong and Lim (DHL, 2008) J. Meteor. Soc. Japan
+! Lin, Farley, Orville (LFO, 1983) J. Appl. Meteor.
+! Rutledge, Hobbs (RH83, 1983) J. Atmos. Sci.
+! Rutledge, Hobbs (RH84, 1984) J. Atmos. Sci.
+! Juang and Hong (JH, 2010) Mon. Wea. Rev.
+!
+ INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
+ ims,ime, jms,jme, kms,kme , &
+ its,ite, jts,jte, kts,kte, &
+ lat
+ REAL, DIMENSION( its:ite , kts:kte ), &
+ INTENT(INOUT) :: &
+ t
+ REAL, DIMENSION( its:ite , kts:kte, 2 ), &
+ INTENT(INOUT) :: &
+ qci
+ REAL, DIMENSION( its:ite , kts:kte, 3 ), &
+ INTENT(INOUT) :: &
+ qrs
+ REAL, DIMENSION( ims:ime , kms:kme ), &
+ INTENT(INOUT) :: &
+ q
+ REAL, DIMENSION( ims:ime , kms:kme ), &
+ INTENT(IN ) :: &
+ den, &
+ p, &
+ delz
+ REAL, INTENT(IN ) :: delt, &
+ g, &
+ cpd, &
+ cpv, &
+ t0c, &
+ den0, &
+ rd, &
+ rv, &
+ ep1, &
+ ep2, &
+ qmin, &
+ XLS, &
+ XLV0, &
+ XLF0, &
+ cliq, &
+ cice, &
+ psat, &
+ denr
+ REAL, DIMENSION( ims:ime ), &
+ INTENT(INOUT) :: rain, &
+ rainncv, &
+ sr
+ REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, &
+ INTENT(INOUT) :: snow, &
+ snowncv
+ REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, &
+ INTENT(INOUT) :: graupel, &
+ graupelncv
+! LOCAL VAR
+ REAL, DIMENSION( its:ite , kts:kte , 3) :: &
+ rh, &
+ qs, &
+ rslope, &
+ rslope2, &
+ rslope3, &
+ rslopeb, &
+ qrs_tmp, &
+ falk, &
+ fall, &
+ work1
+ REAL, DIMENSION( its:ite , kts:kte ) :: &
+ fallc, &
+ falkc, &
+ work1c, &
+ work2c, &
+ workr, &
+ worka
+ REAL, DIMENSION( its:ite , kts:kte ) :: &
+ den_tmp, &
+ delz_tmp
+ REAL, DIMENSION( its:ite , kts:kte ) :: &
+ pigen, &
+ pidep, &
+ pcond, &
+ prevp, &
+ psevp, &
+ pgevp, &
+ psdep, &
+ pgdep, &
+ praut, &
+ psaut, &
+ pgaut, &
+ piacr, &
+ pracw, &
+ praci, &
+ pracs, &
+ psacw, &
+ psaci, &
+ psacr, &
+ pgacw, &
+ pgaci, &
+ pgacr, &
+ pgacs, &
+ paacw, &
+ psmlt, &
+ pgmlt, &
+ pseml, &
+ pgeml
+ REAL, DIMENSION( its:ite , kts:kte ) :: &
+ qsum, &
+ xl, &
+ cpm, &
+ work2, &
+ denfac, &
+ xni, &
+ denqrs1, &
+ denqrs2, &
+ denqrs3, &
+ denqci, &
+ n0sfac
+ REAL, DIMENSION( its:ite ) :: delqrs1, &
+ delqrs2, &
+ delqrs3, &
+ delqi
+ REAL, DIMENSION( its:ite ) :: tstepsnow, &
+ tstepgraup
+ INTEGER, DIMENSION( its:ite ) :: mstep, &
+ numdt
+ LOGICAL, DIMENSION( its:ite ) :: flgcld
+ REAL :: &
+ cpmcal, xlcal, diffus, &
+ viscos, xka, venfac, conden, diffac, &
+ x, y, z, a, b, c, d, e, &
+ qdt, holdrr, holdrs, holdrg, supcol, supcolt, pvt, &
+ coeres, supsat, dtcld, xmi, eacrs, satdt, &
+ qimax, diameter, xni0, roqi0, &
+ fallsum, fallsum_qsi, fallsum_qg, &
+ vt2i,vt2r,vt2s,vt2g,acrfac,egs,egi, &
+ xlwork2, factor, source, value, &
+ xlf, pfrzdtc, pfrzdtr, supice, alpha2, delta2, delta3
+ REAL :: vt2ave
+ REAL :: holdc, holdci
+ INTEGER :: i, j, k, mstepmax, &
+ iprt, latd, lond, loop, loops, ifsat, n, idim, kdim
+! Temporaries used for inlining fpvs function
+ REAL :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
+! variables for optimization
+ REAL, DIMENSION( its:ite ) :: tvec1
+ REAL :: temp
+!
+!=================================================================
+! compute internal functions
+!
+ cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
+ xlcal(x) = xlv0-xlv1*(x-t0c)
+!----------------------------------------------------------------
+! diffus: diffusion coefficient of the water vapor
+! viscos: kinematic viscosity(m2s-1)
+! Optimizatin : A**B => exp(log(A)*(B))
+!
+ diffus(x,y) = 8.794e-5 * exp(log(x)*(1.81)) / y ! 8.794e-5*x**1.81/y
+ viscos(x,y) = 1.496e-6 * (x*sqrt(x)) /(x+120.)/y ! 1.496e-6*x**1.5/(x+120.)/y
+ xka(x,y) = 1.414e3*viscos(x,y)*y
+ diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b))
+ venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333))) &
+ /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
+ conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+d*d/(rv*e)*c/(a*a))
+!
+!
+ idim = ite-its+1
+ kdim = kte-kts+1
+!
+!----------------------------------------------------------------
+! paddint 0 for negative values generated by dynamics
+!
+ do k = kts, kte
+ do i = its, ite
+ qci(i,k,1) = max(qci(i,k,1),0.0)
+ qrs(i,k,1) = max(qrs(i,k,1),0.0)
+ qci(i,k,2) = max(qci(i,k,2),0.0)
+ qrs(i,k,2) = max(qrs(i,k,2),0.0)
+ qrs(i,k,3) = max(qrs(i,k,3),0.0)
+ enddo
+ enddo
+!
+!----------------------------------------------------------------
+! latent heat for phase changes and heat capacity. neglect the
+! changes during microphysical process calculation
+! emanuel(1994)
+!
+ do k = kts, kte
+ do i = its, ite
+ cpm(i,k) = cpmcal(q(i,k))
+ xl(i,k) = xlcal(t(i,k))
+ enddo
+ enddo
+ do k = kts, kte
+ do i = its, ite
+ delz_tmp(i,k) = delz(i,k)
+ den_tmp(i,k) = den(i,k)
+ enddo
+ enddo
+!
+!----------------------------------------------------------------
+! initialize the surface rain, snow, graupel
+!
+ do i = its, ite
+ rainncv(i) = 0.
+ if(PRESENT (snowncv) .AND. PRESENT (snow)) snowncv(i,lat) = 0.
+ if(PRESENT (graupelncv) .AND. PRESENT (graupel)) graupelncv(i,lat) = 0.
+ sr(i) = 0.
+! new local array to catch step snow and graupel
+ tstepsnow(i) = 0.
+ tstepgraup(i) = 0.
+ enddo
+!
+!----------------------------------------------------------------
+! compute the minor time steps.
+!
+ loops = max(nint(delt/dtcldcr),1)
+ dtcld = delt/loops
+ if(delt.le.dtcldcr) dtcld = delt
+!
+ do loop = 1,loops
+!
+!----------------------------------------------------------------
+! initialize the large scale variables
+!
+ do i = its, ite
+ mstep(i) = 1
+ flgcld(i) = .true.
+ enddo
+!
+! do k = kts, kte
+! do i = its, ite
+! denfac(i,k) = sqrt(den0/den(i,k))
+! enddo
+! enddo
+ do k = kts, kte
+ CALL VREC( tvec1(its), den(its,k), ite-its+1)
+ do i = its, ite
+ tvec1(i) = tvec1(i)*den0
+ enddo
+ CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
+ enddo
+!
+! Inline expansion for fpvs
+! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
+! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
+ hsub = xls
+ hvap = xlv0
+ cvap = cpv
+ ttp=t0c+0.01
+ dldt=cvap-cliq
+ xa=-dldt/rv
+ xb=xa+hvap/(rv*ttp)
+ dldti=cvap-cice
+ xai=-dldti/rv
+ xbi=xai+hsub/(rv*ttp)
+ do k = kts, kte
+ do i = its, ite
+ tr=ttp/t(i,k)
+ qs(i,k,1)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
+ qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
+ qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
+ qs(i,k,1) = max(qs(i,k,1),qmin)
+ rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
+ tr=ttp/t(i,k)
+ if(t(i,k).lt.ttp) then
+ qs(i,k,2)=psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
+ else
+ qs(i,k,2)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
+ endif
+ qs(i,k,2) = min(qs(i,k,2),0.99*p(i,k))
+ qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
+ qs(i,k,2) = max(qs(i,k,2),qmin)
+ rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
+ enddo
+ enddo
+!
+!----------------------------------------------------------------
+! initialize the variables for microphysical physics
+!
+!
+ do k = kts, kte
+ do i = its, ite
+ prevp(i,k) = 0.
+ psdep(i,k) = 0.
+ pgdep(i,k) = 0.
+ praut(i,k) = 0.
+ psaut(i,k) = 0.
+ pgaut(i,k) = 0.
+ pracw(i,k) = 0.
+ praci(i,k) = 0.
+ piacr(i,k) = 0.
+ psaci(i,k) = 0.
+ psacw(i,k) = 0.
+ pracs(i,k) = 0.
+ psacr(i,k) = 0.
+ pgacw(i,k) = 0.
+ paacw(i,k) = 0.
+ pgaci(i,k) = 0.
+ pgacr(i,k) = 0.
+ pgacs(i,k) = 0.
+ pigen(i,k) = 0.
+ pidep(i,k) = 0.
+ pcond(i,k) = 0.
+ psmlt(i,k) = 0.
+ pgmlt(i,k) = 0.
+ pseml(i,k) = 0.
+ pgeml(i,k) = 0.
+ psevp(i,k) = 0.
+ pgevp(i,k) = 0.
+ falk(i,k,1) = 0.
+ falk(i,k,2) = 0.
+ falk(i,k,3) = 0.
+ fall(i,k,1) = 0.
+ fall(i,k,2) = 0.
+ fall(i,k,3) = 0.
+ fallc(i,k) = 0.
+ falkc(i,k) = 0.
+ xni(i,k) = 1.e3
+ enddo
+ enddo
+!-------------------------------------------------------------
+! Ni: ice crystal number concentraiton [HDC 5c]
+!-------------------------------------------------------------
+ do k = kts, kte
+ do i = its, ite
+ temp = (den(i,k)*max(qci(i,k,2),qmin))
+ temp = sqrt(sqrt(temp*temp*temp))
+ xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
+ enddo
+ enddo
+!
+!----------------------------------------------------------------
+! compute the fallout term:
+! first, vertical terminal velosity for minor loops
+!----------------------------------------------------------------
+ do k = kts, kte
+ do i = its, ite
+ qrs_tmp(i,k,1) = qrs(i,k,1)
+ qrs_tmp(i,k,2) = qrs(i,k,2)
+ qrs_tmp(i,k,3) = qrs(i,k,3)
+ enddo
+ enddo
+ call slope_wsm6(qrs_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2,rslope3, &
+ work1,its,ite,kts,kte)
+!
+ do k = kte, kts, -1
+ do i = its, ite
+ workr(i,k) = work1(i,k,1)
+ qsum(i,k) = max( (qrs(i,k,2)+qrs(i,k,3)), 1.E-15)
+ IF ( qsum(i,k) .gt. 1.e-15 ) THEN
+ worka(i,k) = (work1(i,k,2)*qrs(i,k,2) + work1(i,k,3)*qrs(i,k,3)) &
+ /qsum(i,k)
+ ELSE
+ worka(i,k) = 0.
+ ENDIF
+ denqrs1(i,k) = den(i,k)*qrs(i,k,1)
+ denqrs2(i,k) = den(i,k)*qrs(i,k,2)
+ denqrs3(i,k) = den(i,k)*qrs(i,k,3)
+ if(qrs(i,k,1).le.0.0) workr(i,k) = 0.0
+ enddo
+ enddo
+ call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,workr,denqrs1, &
+ delqrs1,dtcld,1,1)
+ call nislfv_rain_plm6(idim,kdim,den_tmp,denfac,t,delz_tmp,worka, &
+ denqrs2,denqrs3,delqrs2,delqrs3,dtcld,1,1)
+ do k = kts, kte
+ do i = its, ite
+ qrs(i,k,1) = max(denqrs1(i,k)/den(i,k),0.)
+ qrs(i,k,2) = max(denqrs2(i,k)/den(i,k),0.)
+ qrs(i,k,3) = max(denqrs3(i,k)/den(i,k),0.)
+ fall(i,k,1) = denqrs1(i,k)*workr(i,k)/delz(i,k)
+ fall(i,k,2) = denqrs2(i,k)*worka(i,k)/delz(i,k)
+ fall(i,k,3) = denqrs3(i,k)*worka(i,k)/delz(i,k)
+ enddo
+ enddo
+ do i = its, ite
+ fall(i,1,1) = delqrs1(i)/delz(i,1)/dtcld
+ fall(i,1,2) = delqrs2(i)/delz(i,1)/dtcld
+ fall(i,1,3) = delqrs3(i)/delz(i,1)/dtcld
+ enddo
+ do k = kts, kte
+ do i = its, ite
+ qrs_tmp(i,k,1) = qrs(i,k,1)
+ qrs_tmp(i,k,2) = qrs(i,k,2)
+ qrs_tmp(i,k,3) = qrs(i,k,3)
+ enddo
+ enddo
+ call slope_wsm6(qrs_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2,rslope3, &
+ work1,its,ite,kts,kte)
+!
+ do k = kte, kts, -1
+ do i = its, ite
+ supcol = t0c-t(i,k)
+ n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
+ if(t(i,k).gt.t0c) then
+!---------------------------------------------------------------
+! psmlt: melting of snow [HL A33] [RH83 A25]
+! (T>T0: S->R)
+!---------------------------------------------------------------
+ xlf = xlf0
+ work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
+ if(qrs(i,k,2).gt.0.) then
+ coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
+ psmlt(i,k) = xka(t(i,k),den(i,k))/xlf*(t0c-t(i,k))*pi/2. &
+ *n0sfac(i,k)*(precs1*rslope2(i,k,2) &
+ +precs2*work2(i,k)*coeres)
+ psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep(i), &
+ -qrs(i,k,2)/mstep(i)),0.)
+ qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
+ qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k)
+ t(i,k) = t(i,k) + xlf/cpm(i,k)*psmlt(i,k)
+ endif
+!---------------------------------------------------------------
+! pgmlt: melting of graupel [HL A23] [LFO 47]
+! (T>T0: G->R)
+!---------------------------------------------------------------
+ if(qrs(i,k,3).gt.0.) then
+ coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3))
+ pgmlt(i,k) = xka(t(i,k),den(i,k))/xlf &
+ *(t0c-t(i,k))*(precg1*rslope2(i,k,3) &
+ +precg2*work2(i,k)*coeres)
+ pgmlt(i,k) = min(max(pgmlt(i,k)*dtcld/mstep(i), &
+ -qrs(i,k,3)/mstep(i)),0.)
+ qrs(i,k,3) = qrs(i,k,3) + pgmlt(i,k)
+ qrs(i,k,1) = qrs(i,k,1) - pgmlt(i,k)
+ t(i,k) = t(i,k) + xlf/cpm(i,k)*pgmlt(i,k)
+ endif
+ endif
+ enddo
+ enddo
+!---------------------------------------------------------------
+! Vice [ms-1] : fallout of ice crystal [HDC 5a]
+!---------------------------------------------------------------
+ do k = kte, kts, -1
+ do i = its, ite
+ if(qci(i,k,2).le.0.) then
+ work1c(i,k) = 0.
+ else
+ xmi = den(i,k)*qci(i,k,2)/xni(i,k)
+ diameter = max(min(dicon * sqrt(xmi),dimax), 1.e-25)
+ work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31))
+ endif
+ enddo
+ enddo
+!
+! forward semi-laglangian scheme (JH), PCM (piecewise constant), (linear)
+!
+ do k = kte, kts, -1
+ do i = its, ite
+ denqci(i,k) = den(i,k)*qci(i,k,2)
+ enddo
+ enddo
+ call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,work1c,denqci, &
+ delqi,dtcld,1,0)
+ do k = kts, kte
+ do i = its, ite
+ qci(i,k,2) = max(denqci(i,k)/den(i,k),0.)
+ enddo
+ enddo
+ do i = its, ite
+ fallc(i,1) = delqi(i)/delz(i,1)/dtcld
+ enddo
+!
+!----------------------------------------------------------------
+! rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
+!
+ do i = its, ite
+ fallsum = fall(i,kts,1)+fall(i,kts,2)+fall(i,kts,3)+fallc(i,kts)
+ fallsum_qsi = fall(i,kts,2)+fallc(i,kts)
+ fallsum_qg = fall(i,kts,3)
+ if(fallsum.gt.0.) then
+ rainncv(i) = fallsum*delz(i,kts)/denr*dtcld*1000. + rainncv(i)
+ rain(i) = fallsum*delz(i,kts)/denr*dtcld*1000. + rain(i)
+ endif
+ if(fallsum_qsi.gt.0.) then
+ tstepsnow(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. &
+ +tstepsnow(i)
+ IF ( PRESENT (snowncv) .AND. PRESENT (snow)) THEN
+ snowncv(i,lat) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. &
+ +snowncv(i,lat)
+ snow(i,lat) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + snow(i,lat)
+ ENDIF
+ endif
+ if(fallsum_qg.gt.0.) then
+ tstepgraup(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. &
+ +tstepgraup(i)
+ IF ( PRESENT (graupelncv) .AND. PRESENT (graupel)) THEN
+ graupelncv(i,lat) = fallsum_qg*delz(i,kts)/denr*dtcld*1000. &
+ + graupelncv(i,lat)
+ graupel(i,lat) = fallsum_qg*delz(i,kts)/denr*dtcld*1000. + graupel(i,lat)
+ ENDIF
+ endif
+! if(fallsum.gt.0.)sr(i)=(snowncv(i,lat) + graupelncv(i,lat))/(rainncv(i)+1.e-12)
+ if(fallsum.gt.0.)sr(i)=(tstepsnow(i) + tstepgraup(i))/(rainncv(i)+1.e-12)
+ enddo
+!
+!---------------------------------------------------------------
+! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
+! (T>T0: I->C)
+!---------------------------------------------------------------
+ do k = kts, kte
+ do i = its, ite
+ supcol = t0c-t(i,k)
+ xlf = xls-xl(i,k)
+ if(supcol.lt.0.) xlf = xlf0
+ if(supcol.lt.0.and.qci(i,k,2).gt.0.) then
+ qci(i,k,1) = qci(i,k,1) + qci(i,k,2)
+ t(i,k) = t(i,k) - xlf/cpm(i,k)*qci(i,k,2)
+ qci(i,k,2) = 0.
+ endif
+!---------------------------------------------------------------
+! pihmf: homogeneous freezing of cloud water below -40c [HL A45]
+! (T<-40C: C->I)
+!---------------------------------------------------------------
+ if(supcol.gt.40..and.qci(i,k,1).gt.0.) then
+ qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
+ t(i,k) = t(i,k) + xlf/cpm(i,k)*qci(i,k,1)
+ qci(i,k,1) = 0.
+ endif
+!---------------------------------------------------------------
+! pihtf: heterogeneous freezing of cloud water [HL A44]
+! (T0>T>-40C: C->I)
+!---------------------------------------------------------------
+ if(supcol.gt.0..and.qci(i,k,1).gt.qmin) then
+! pfrzdtc = min(pfrz1*(exp(pfrz2*supcol)-1.) &
+! *den(i,k)/denr/xncr*qci(i,k,1)**2*dtcld,qci(i,k,1))
+ supcolt=min(supcol,50.)
+ pfrzdtc = min(pfrz1*(exp(pfrz2*supcolt)-1.) &
+ *den(i,k)/denr/xncr*qci(i,k,1)*qci(i,k,1)*dtcld,qci(i,k,1))
+ qci(i,k,2) = qci(i,k,2) + pfrzdtc
+ t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtc
+ qci(i,k,1) = qci(i,k,1)-pfrzdtc
+ endif
+!---------------------------------------------------------------
+! pgfrz: freezing of rain water [HL A20] [LFO 45]
+! (T<T0, R->G)
+!---------------------------------------------------------------
+ if(supcol.gt.0..and.qrs(i,k,1).gt.0.) then
+! pfrzdtr = min(20.*pi**2*pfrz1*n0r*denr/den(i,k) &
+! *(exp(pfrz2*supcol)-1.)*rslope3(i,k,1)**2 &
+! *rslope(i,k,1)*dtcld,qrs(i,k,1))
+ temp = rslope3(i,k,1)
+ temp = temp*temp*rslope(i,k,1)
+ supcolt=min(supcol,50.)
+ pfrzdtr = min(20.*(pi*pi)*pfrz1*n0r*denr/den(i,k) &
+ *(exp(pfrz2*supcolt)-1.)*temp*dtcld, &
+ qrs(i,k,1))
+ qrs(i,k,3) = qrs(i,k,3) + pfrzdtr
+ t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtr
+ qrs(i,k,1) = qrs(i,k,1)-pfrzdtr
+ endif
+ enddo
+ enddo
+!
+!
+!----------------------------------------------------------------
+! update the slope parameters for microphysics computation
+!
+ do k = kts, kte
+ do i = its, ite
+ qrs_tmp(i,k,1) = qrs(i,k,1)
+ qrs_tmp(i,k,2) = qrs(i,k,2)
+ qrs_tmp(i,k,3) = qrs(i,k,3)
+ enddo
+ enddo
+ call slope_wsm6(qrs_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2,rslope3, &
+ work1,its,ite,kts,kte)
+!------------------------------------------------------------------
+! work1: the thermodynamic term in the denominator associated with
+! heat conduction and vapor diffusion
+! (ry88, y93, h85)
+! work2: parameter associated with the ventilation effects(y93)
+!
+ do k = kts, kte
+ do i = its, ite
+ work1(i,k,1) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k,1))
+ work1(i,k,2) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k,2))
+ work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
+ enddo
+ enddo
+!
+!===============================================================
+!
+! warm rain processes
+!
+! - follows the processes in RH83 and LFO except for autoconcersion
+!
+!===============================================================
+!
+ do k = kts, kte
+ do i = its, ite
+ supsat = max(q(i,k),qmin)-qs(i,k,1)
+ satdt = supsat/dtcld
+!---------------------------------------------------------------
+! praut: auto conversion rate from cloud to rain [HDC 16]
+! (C->R)
+!---------------------------------------------------------------
+ if(qci(i,k,1).gt.qc0) then
+ praut(i,k) = qck1*qci(i,k,1)**(7./3.)
+ praut(i,k) = min(praut(i,k),qci(i,k,1)/dtcld)
+ endif
+!---------------------------------------------------------------
+! pracw: accretion of cloud water by rain [HL A40] [LFO 51]
+! (C->R)
+!---------------------------------------------------------------
+ if(qrs(i,k,1).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
+ pracw(i,k) = min(pacrr*rslope3(i,k,1)*rslopeb(i,k,1) &
+ *qci(i,k,1)*denfac(i,k),qci(i,k,1)/dtcld)
+ endif
+!---------------------------------------------------------------
+! prevp: evaporation/condensation rate of rain [HDC 14]
+! (V->R or R->V)
+!---------------------------------------------------------------
+ if(qrs(i,k,1).gt.0.) then
+ coeres = rslope2(i,k,1)*sqrt(rslope(i,k,1)*rslopeb(i,k,1))
+ prevp(i,k) = (rh(i,k,1)-1.)*(precr1*rslope2(i,k,1) &
+ +precr2*work2(i,k)*coeres)/work1(i,k,1)
+ if(prevp(i,k).lt.0.) then
+ prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
+ prevp(i,k) = max(prevp(i,k),satdt/2)
+ else
+ prevp(i,k) = min(prevp(i,k),satdt/2)
+ endif
+ endif
+ enddo
+ enddo
+!
+!===============================================================
+!
+! cold rain processes
+!
+! - follows the revised ice microphysics processes in HDC
+! - the processes same as in RH83 and RH84 and LFO behave
+! following ice crystal hapits defined in HDC, inclduing
+! intercept parameter for snow (n0s), ice crystal number
+! concentration (ni), ice nuclei number concentration
+! (n0i), ice diameter (d)
+!
+!===============================================================
+!
+ do k = kts, kte
+ do i = its, ite
+ supcol = t0c-t(i,k)
+ n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
+ supsat = max(q(i,k),qmin)-qs(i,k,2)
+ satdt = supsat/dtcld
+ ifsat = 0
+!-------------------------------------------------------------
+! Ni: ice crystal number concentraiton [HDC 5c]
+!-------------------------------------------------------------
+! xni(i,k) = min(max(5.38e7*(den(i,k) &
+! *max(qci(i,k,2),qmin))**0.75,1.e3),1.e6)
+ temp = (den(i,k)*max(qci(i,k,2),qmin))
+ temp = sqrt(sqrt(temp*temp*temp))
+ xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
+ eacrs = exp(0.07*(-supcol))
+!
+ xmi = den(i,k)*qci(i,k,2)/xni(i,k)
+ diameter = min(dicon * sqrt(xmi),dimax)
+ vt2i = 1.49e4*diameter**1.31
+ vt2r=pvtr*rslopeb(i,k,1)*denfac(i,k)
+ vt2s=pvts*rslopeb(i,k,2)*denfac(i,k)
+ vt2g=pvtg*rslopeb(i,k,3)*denfac(i,k)
+ qsum(i,k) = max( (qrs(i,k,2)+qrs(i,k,3)), 1.E-15)
+ if(qsum(i,k) .gt. 1.e-15) then
+ vt2ave=(vt2s*qrs(i,k,2)+vt2g*qrs(i,k,3))/(qsum(i,k))
+ else
+ vt2ave=0.
+ endif
+ if(supcol.gt.0.and.qci(i,k,2).gt.qmin) then
+ if(qrs(i,k,1).gt.qcrmin) then
+!-------------------------------------------------------------
+! praci: Accretion of cloud ice by rain [HL A15] [LFO 25]
+! (T<T0: I->R)
+!-------------------------------------------------------------
+ acrfac = 2.*rslope3(i,k,1)+2.*diameter*rslope2(i,k,1) &
+ +diameter**2*rslope(i,k,1)
+ praci(i,k) = pi*qci(i,k,2)*n0r*abs(vt2r-vt2i)*acrfac/4.
+ praci(i,k) = min(praci(i,k),qci(i,k,2)/dtcld)
+!-------------------------------------------------------------
+! piacr: Accretion of rain by cloud ice [HL A19] [LFO 26]
+! (T<T0: R->S or R->G)
+!-------------------------------------------------------------
+ piacr(i,k) = pi**2*avtr*n0r*denr*xni(i,k)*denfac(i,k) &
+ *g6pbr*rslope3(i,k,1)*rslope3(i,k,1) &
+ *rslopeb(i,k,1)/24./den(i,k)
+ piacr(i,k) = min(piacr(i,k),qrs(i,k,1)/dtcld)
+ endif
+!-------------------------------------------------------------
+! psaci: Accretion of cloud ice by snow [HDC 10]
+! (T<T0: I->S)
+!-------------------------------------------------------------
+ if(qrs(i,k,2).gt.qcrmin) then
+ acrfac = 2.*rslope3(i,k,2)+2.*diameter*rslope2(i,k,2) &
+ +diameter**2*rslope(i,k,2)
+ psaci(i,k) = pi*qci(i,k,2)*eacrs*n0s*n0sfac(i,k) &
+ *abs(vt2ave-vt2i)*acrfac/4.
+ psaci(i,k) = min(psaci(i,k),qci(i,k,2)/dtcld)
+ endif
+!-------------------------------------------------------------
+! pgaci: Accretion of cloud ice by graupel [HL A17] [LFO 41]
+! (T<T0: I->G)
+!-------------------------------------------------------------
+ if(qrs(i,k,3).gt.qcrmin) then
+ egi = exp(0.07*(-supcol))
+ acrfac = 2.*rslope3(i,k,3)+2.*diameter*rslope2(i,k,3) &
+ +diameter**2*rslope(i,k,3)
+ pgaci(i,k) = pi*egi*qci(i,k,2)*n0g*abs(vt2ave-vt2i)*acrfac/4.
+ pgaci(i,k) = min(pgaci(i,k),qci(i,k,2)/dtcld)
+ endif
+ endif
+!-------------------------------------------------------------
+! psacw: Accretion of cloud water by snow [HL A7] [LFO 24]
+! (T<T0: C->S, and T>=T0: C->R)
+!-------------------------------------------------------------
+ if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
+ psacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2)*rslopeb(i,k,2) &
+ *qci(i,k,1)*denfac(i,k),qci(i,k,1)/dtcld)
+ endif
+!-------------------------------------------------------------
+! pgacw: Accretion of cloud water by graupel [HL A6] [LFO 40]
+! (T<T0: C->G, and T>=T0: C->R)
+!-------------------------------------------------------------
+ if(qrs(i,k,3).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
+ pgacw(i,k) = min(pacrg*rslope3(i,k,3)*rslopeb(i,k,3) &
+ *qci(i,k,1)*denfac(i,k),qci(i,k,1)/dtcld)
+ endif
+!-------------------------------------------------------------
+! paacw: Accretion of cloud water by averaged snow/graupel
+! (T<T0: C->G or S, and T>=T0: C->R)
+!-------------------------------------------------------------
+ if(qrs(i,k,2).gt.qcrmin.and.qrs(i,k,3).gt.qcrmin) then
+ paacw(i,k) = (qrs(i,k,2)*psacw(i,k)+qrs(i,k,3)*pgacw(i,k)) &
+ /(qsum(i,k))
+ endif
+!-------------------------------------------------------------
+! pracs: Accretion of snow by rain [HL A11] [LFO 27]
+! (T<T0: S->G)
+!-------------------------------------------------------------
+ if(qrs(i,k,2).gt.qcrmin.and.qrs(i,k,1).gt.qcrmin) then
+ if(supcol.gt.0) then
+ acrfac = 5.*rslope3(i,k,2)*rslope3(i,k,2)*rslope(i,k,1) &
+ +2.*rslope3(i,k,2)*rslope2(i,k,2)*rslope2(i,k,1) &
+ +.5*rslope2(i,k,2)*rslope2(i,k,2)*rslope3(i,k,1)
+ pracs(i,k) = pi**2*n0r*n0s*n0sfac(i,k)*abs(vt2r-vt2ave) &
+ *(dens/den(i,k))*acrfac
+ pracs(i,k) = min(pracs(i,k),qrs(i,k,2)/dtcld)
+ endif
+!-------------------------------------------------------------
+! psacr: Accretion of rain by snow [HL A10] [LFO 28]
+! (T<T0:R->S or R->G) (T>=T0: enhance melting of snow)
+!-------------------------------------------------------------
+ acrfac = 5.*rslope3(i,k,1)*rslope3(i,k,1)*rslope(i,k,2) &
+ +2.*rslope3(i,k,1)*rslope2(i,k,1)*rslope2(i,k,2) &
+ +.5*rslope2(i,k,1)*rslope2(i,k,1)*rslope3(i,k,2)
+ psacr(i,k) = pi**2*n0r*n0s*n0sfac(i,k)*abs(vt2ave-vt2r) &
+ *(denr/den(i,k))*acrfac
+ psacr(i,k) = min(psacr(i,k),qrs(i,k,1)/dtcld)
+ endif
+!-------------------------------------------------------------
+! pgacr: Accretion of rain by graupel [HL A12] [LFO 42]
+! (T<T0: R->G) (T>=T0: enhance melting of graupel)
+!-------------------------------------------------------------
+ if(qrs(i,k,3).gt.qcrmin.and.qrs(i,k,1).gt.qcrmin) then
+ acrfac = 5.*rslope3(i,k,1)*rslope3(i,k,1)*rslope(i,k,3) &
+ +2.*rslope3(i,k,1)*rslope2(i,k,1)*rslope2(i,k,3) &
+ +.5*rslope2(i,k,1)*rslope2(i,k,1)*rslope3(i,k,3)
+ pgacr(i,k) = pi**2*n0r*n0g*abs(vt2ave-vt2r)*(denr/den(i,k)) &
+ *acrfac
+ pgacr(i,k) = min(pgacr(i,k),qrs(i,k,1)/dtcld)
+ endif
+!
+!-------------------------------------------------------------
+! pgacs: Accretion of snow by graupel [HL A13] [LFO 29]
+! (S->G): This process is eliminated in V3.0 with the
+! new combined snow/graupel fall speeds
+!-------------------------------------------------------------
+ if(qrs(i,k,3).gt.qcrmin.and.qrs(i,k,2).gt.qcrmin) then
+ pgacs(i,k) = 0.
+ endif
+ if(supcol.le.0) then
+ xlf = xlf0
+!-------------------------------------------------------------
+! pseml: Enhanced melting of snow by accretion of water [HL A34]
+! (T>=T0: S->R)
+!-------------------------------------------------------------
+ if(qrs(i,k,2).gt.0.) &
+ pseml(i,k) = min(max(cliq*supcol*(paacw(i,k)+psacr(i,k)) &
+ /xlf,-qrs(i,k,2)/dtcld),0.)
+!-------------------------------------------------------------
+! pgeml: Enhanced melting of graupel by accretion of water [HL A24] [RH84 A21-A22]
+! (T>=T0: G->R)
+!-------------------------------------------------------------
+ if(qrs(i,k,3).gt.0.) &
+ pgeml(i,k) = min(max(cliq*supcol*(paacw(i,k)+pgacr(i,k)) &
+ /xlf,-qrs(i,k,3)/dtcld),0.)
+ endif
+ if(supcol.gt.0) then
+!-------------------------------------------------------------
+! pidep: Deposition/Sublimation rate of ice [HDC 9]
+! (T<T0: V->I or I->V)
+!-------------------------------------------------------------
+ if(qci(i,k,2).gt.0.and.ifsat.ne.1) then
+ pidep(i,k) = 4.*diameter*xni(i,k)*(rh(i,k,2)-1.)/work1(i,k,2)
+ supice = satdt-prevp(i,k)
+ if(pidep(i,k).lt.0.) then
+ pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
+ pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
+ else
+ pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
+ endif
+ if(abs(prevp(i,k)+pidep(i,k)).ge.abs(satdt)) ifsat = 1
+ endif
+!-------------------------------------------------------------
+! psdep: deposition/sublimation rate of snow [HDC 14]
+! (T<T0: V->S or S->V)
+!-------------------------------------------------------------
+ if(qrs(i,k,2).gt.0..and.ifsat.ne.1) then
+ coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
+ psdep(i,k) = (rh(i,k,2)-1.)*n0sfac(i,k)*(precs1*rslope2(i,k,2) &
+ + precs2*work2(i,k)*coeres)/work1(i,k,2)
+ supice = satdt-prevp(i,k)-pidep(i,k)
+ if(psdep(i,k).lt.0.) then
+ psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
+ psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
+ else
+ psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
+ endif
+ if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt)) &
+ ifsat = 1
+ endif
+!-------------------------------------------------------------
+! pgdep: deposition/sublimation rate of graupel [HL A21] [LFO 46]
+! (T<T0: V->G or G->V)
+!-------------------------------------------------------------
+ if(qrs(i,k,3).gt.0..and.ifsat.ne.1) then
+ coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3))
+ pgdep(i,k) = (rh(i,k,2)-1.)*(precg1*rslope2(i,k,3) &
+ +precg2*work2(i,k)*coeres)/work1(i,k,2)
+ supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
+ if(pgdep(i,k).lt.0.) then
+ pgdep(i,k) = max(pgdep(i,k),-qrs(i,k,3)/dtcld)
+ pgdep(i,k) = max(max(pgdep(i,k),satdt/2),supice)
+ else
+ pgdep(i,k) = min(min(pgdep(i,k),satdt/2),supice)
+ endif
+ if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)+pgdep(i,k)).ge. &
+ abs(satdt)) ifsat = 1
+ endif
+!-------------------------------------------------------------
+! pigen: generation(nucleation) of ice from vapor [HL 50] [HDC 7-8]
+! (T<T0: V->I)
+!-------------------------------------------------------------
+ if(supsat.gt.0.and.ifsat.ne.1) then
+ supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)-pgdep(i,k)
+ xni0 = 1.e3*exp(0.1*supcol)
+ roqi0 = 4.92e-11*xni0**1.33
+ pigen(i,k) = max(0.,(roqi0/den(i,k)-max(qci(i,k,2),0.))/dtcld)
+ pigen(i,k) = min(min(pigen(i,k),satdt),supice)
+ endif
+!
+!-------------------------------------------------------------
+! psaut: conversion(aggregation) of ice to snow [HDC 12]
+! (T<T0: I->S)
+!-------------------------------------------------------------
+ if(qci(i,k,2).gt.0.) then
+ qimax = roqimax/den(i,k)
+ psaut(i,k) = max(0.,(qci(i,k,2)-qimax)/dtcld)
+ endif
+!
+!-------------------------------------------------------------
+! pgaut: conversion(aggregation) of snow to graupel [HL A4] [LFO 37]
+! (T<T0: S->G)
+!-------------------------------------------------------------
+ if(qrs(i,k,2).gt.0.) then
+ alpha2 = 1.e-3*exp(0.09*(-supcol))
+ pgaut(i,k) = min(max(0.,alpha2*(qrs(i,k,2)-qs0)),qrs(i,k,2)/dtcld)
+ endif
+ endif
+!
+!-------------------------------------------------------------
+! psevp: Evaporation of melting snow [HL A35] [RH83 A27]
+! (T>=T0: S->V)
+!-------------------------------------------------------------
+ if(supcol.lt.0.) then
+ if(qrs(i,k,2).gt.0..and.rh(i,k,1).lt.1.) then
+ coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
+ psevp(i,k) = (rh(i,k,1)-1.)*n0sfac(i,k)*(precs1 &
+ *rslope2(i,k,2)+precs2*work2(i,k) &
+ *coeres)/work1(i,k,1)
+ psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)/dtcld),0.)
+ endif
+!-------------------------------------------------------------
+! pgevp: Evaporation of melting graupel [HL A25] [RH84 A19]
+! (T>=T0: G->V)
+!-------------------------------------------------------------
+ if(qrs(i,k,3).gt.0..and.rh(i,k,1).lt.1.) then
+ coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3))
+ pgevp(i,k) = (rh(i,k,1)-1.)*(precg1*rslope2(i,k,3) &
+ +precg2*work2(i,k)*coeres)/work1(i,k,1)
+ pgevp(i,k) = min(max(pgevp(i,k),-qrs(i,k,3)/dtcld),0.)
+ endif
+ endif
+ enddo
+ enddo
+!
+!
+!----------------------------------------------------------------
+! check mass conservation of generation terms and feedback to the
+! large scale
+!
+ do k = kts, kte
+ do i = its, ite
+!
+ delta2=0.
+ delta3=0.
+ if(qrs(i,k,1).lt.1.e-4.and.qrs(i,k,2).lt.1.e-4) delta2=1.
+ if(qrs(i,k,1).lt.1.e-4) delta3=1.
+ if(t(i,k).le.t0c) then
+!
+! cloud water
+!
+ value = max(qmin,qci(i,k,1))
+ source = (praut(i,k)+pracw(i,k)+paacw(i,k)+paacw(i,k))*dtcld
+ if (source.gt.value) then
+ factor = value/source
+ praut(i,k) = praut(i,k)*factor
+ pracw(i,k) = pracw(i,k)*factor
+ paacw(i,k) = paacw(i,k)*factor
+ endif
+!
+! cloud ice
+!
+ value = max(qmin,qci(i,k,2))
+ source = (psaut(i,k)-pigen(i,k)-pidep(i,k)+praci(i,k)+psaci(i,k) &
+ +pgaci(i,k))*dtcld
+ if (source.gt.value) then
+ factor = value/source
+ psaut(i,k) = psaut(i,k)*factor
+ pigen(i,k) = pigen(i,k)*factor
+ pidep(i,k) = pidep(i,k)*factor
+ praci(i,k) = praci(i,k)*factor
+ psaci(i,k) = psaci(i,k)*factor
+ pgaci(i,k) = pgaci(i,k)*factor
+ endif
+!
+! rain
+!
+ value = max(qmin,qrs(i,k,1))
+ source = (-praut(i,k)-prevp(i,k)-pracw(i,k)+piacr(i,k)+psacr(i,k) &
+ +pgacr(i,k))*dtcld
+ if (source.gt.value) then
+ factor = value/source
+ praut(i,k) = praut(i,k)*factor
+ prevp(i,k) = prevp(i,k)*factor
+ pracw(i,k) = pracw(i,k)*factor
+ piacr(i,k) = piacr(i,k)*factor
+ psacr(i,k) = psacr(i,k)*factor
+ pgacr(i,k) = pgacr(i,k)*factor
+ endif
+!
+! snow
+!
+ value = max(qmin,qrs(i,k,2))
+ source = -(psdep(i,k)+psaut(i,k)-pgaut(i,k)+paacw(i,k)+piacr(i,k) &
+ *delta3+praci(i,k)*delta3-pracs(i,k)*(1.-delta2) &
+ +psacr(i,k)*delta2+psaci(i,k)-pgacs(i,k) )*dtcld
+ if (source.gt.value) then
+ factor = value/source
+ psdep(i,k) = psdep(i,k)*factor
+ psaut(i,k) = psaut(i,k)*factor
+ pgaut(i,k) = pgaut(i,k)*factor
+ paacw(i,k) = paacw(i,k)*factor
+ piacr(i,k) = piacr(i,k)*factor
+ praci(i,k) = praci(i,k)*factor
+ psaci(i,k) = psaci(i,k)*factor
+ pracs(i,k) = pracs(i,k)*factor
+ psacr(i,k) = psacr(i,k)*factor
+ pgacs(i,k) = pgacs(i,k)*factor
+ endif
+!
+! graupel
+!
+ value = max(qmin,qrs(i,k,3))
+ source = -(pgdep(i,k)+pgaut(i,k) &
+ +piacr(i,k)*(1.-delta3)+praci(i,k)*(1.-delta3) &
+ +psacr(i,k)*(1.-delta2)+pracs(i,k)*(1.-delta2) &
+ +pgaci(i,k)+paacw(i,k)+pgacr(i,k)+pgacs(i,k))*dtcld
+ if (source.gt.value) then
+ factor = value/source
+ pgdep(i,k) = pgdep(i,k)*factor
+ pgaut(i,k) = pgaut(i,k)*factor
+ piacr(i,k) = piacr(i,k)*factor
+ praci(i,k) = praci(i,k)*factor
+ psacr(i,k) = psacr(i,k)*factor
+ pracs(i,k) = pracs(i,k)*factor
+ paacw(i,k) = paacw(i,k)*factor
+ pgaci(i,k) = pgaci(i,k)*factor
+ pgacr(i,k) = pgacr(i,k)*factor
+ pgacs(i,k) = pgacs(i,k)*factor
+ endif
+!
+ work2(i,k)=-(prevp(i,k)+psdep(i,k)+pgdep(i,k)+pigen(i,k)+pidep(i,k))
+! update
+ q(i,k) = q(i,k)+work2(i,k)*dtcld
+ qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
+ +paacw(i,k)+paacw(i,k))*dtcld,0.)
+ qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
+ +prevp(i,k)-piacr(i,k)-pgacr(i,k) &
+ -psacr(i,k))*dtcld,0.)
+ qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+praci(i,k) &
+ +psaci(i,k)+pgaci(i,k)-pigen(i,k)-pidep(i,k)) &
+ *dtcld,0.)
+ qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k)+paacw(i,k) &
+ -pgaut(i,k)+piacr(i,k)*delta3 &
+ +praci(i,k)*delta3+psaci(i,k)-pgacs(i,k) &
+ -pracs(i,k)*(1.-delta2)+psacr(i,k)*delta2) &
+ *dtcld,0.)
+ qrs(i,k,3) = max(qrs(i,k,3)+(pgdep(i,k)+pgaut(i,k) &
+ +piacr(i,k)*(1.-delta3) &
+ +praci(i,k)*(1.-delta3)+psacr(i,k)*(1.-delta2) &
+ +pracs(i,k)*(1.-delta2)+pgaci(i,k)+paacw(i,k) &
+ +pgacr(i,k)+pgacs(i,k))*dtcld,0.)
+ xlf = xls-xl(i,k)
+ xlwork2 = -xls*(psdep(i,k)+pgdep(i,k)+pidep(i,k)+pigen(i,k)) &
+ -xl(i,k)*prevp(i,k)-xlf*(piacr(i,k)+paacw(i,k) &
+ +paacw(i,k)+pgacr(i,k)+psacr(i,k))
+ t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
+ else
+!
+! cloud water
+!
+ value = max(qmin,qci(i,k,1))
+ source=(praut(i,k)+pracw(i,k)+paacw(i,k)+paacw(i,k))*dtcld
+ if (source.gt.value) then
+ factor = value/source
+ praut(i,k) = praut(i,k)*factor
+ pracw(i,k) = pracw(i,k)*factor
+ paacw(i,k) = paacw(i,k)*factor
+ endif
+!
+! rain
+!
+ value = max(qmin,qrs(i,k,1))
+ source = (-paacw(i,k)-praut(i,k)+pseml(i,k)+pgeml(i,k)-pracw(i,k) &
+ -paacw(i,k)-prevp(i,k))*dtcld
+ if (source.gt.value) then
+ factor = value/source
+ praut(i,k) = praut(i,k)*factor
+ prevp(i,k) = prevp(i,k)*factor
+ pracw(i,k) = pracw(i,k)*factor
+ paacw(i,k) = paacw(i,k)*factor
+ pseml(i,k) = pseml(i,k)*factor
+ pgeml(i,k) = pgeml(i,k)*factor
+ endif
+!
+! snow
+!
+ value = max(qcrmin,qrs(i,k,2))
+ source=(pgacs(i,k)-pseml(i,k)-psevp(i,k))*dtcld
+ if (source.gt.value) then
+ factor = value/source
+ pgacs(i,k) = pgacs(i,k)*factor
+ psevp(i,k) = psevp(i,k)*factor
+ pseml(i,k) = pseml(i,k)*factor
+ endif
+!
+! graupel
+!
+ value = max(qcrmin,qrs(i,k,3))
+ source=-(pgacs(i,k)+pgevp(i,k)+pgeml(i,k))*dtcld
+ if (source.gt.value) then
+ factor = value/source
+ pgacs(i,k) = pgacs(i,k)*factor
+ pgevp(i,k) = pgevp(i,k)*factor
+ pgeml(i,k) = pgeml(i,k)*factor
+ endif
+ work2(i,k)=-(prevp(i,k)+psevp(i,k)+pgevp(i,k))
+! update
+ q(i,k) = q(i,k)+work2(i,k)*dtcld
+ qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
+ +paacw(i,k)+paacw(i,k))*dtcld,0.)
+ qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
+ +prevp(i,k)+paacw(i,k)+paacw(i,k)-pseml(i,k) &
+ -pgeml(i,k))*dtcld,0.)
+ qrs(i,k,2) = max(qrs(i,k,2)+(psevp(i,k)-pgacs(i,k) &
+ +pseml(i,k))*dtcld,0.)
+ qrs(i,k,3) = max(qrs(i,k,3)+(pgacs(i,k)+pgevp(i,k) &
+ +pgeml(i,k))*dtcld,0.)
+ xlf = xls-xl(i,k)
+ xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k)+pgevp(i,k)) &
+ -xlf*(pseml(i,k)+pgeml(i,k))
+ t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
+ endif
+ enddo
+ enddo
+!
+! Inline expansion for fpvs
+! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
+! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
+ hsub = xls
+ hvap = xlv0
+ cvap = cpv
+ ttp=t0c+0.01
+ dldt=cvap-cliq
+ xa=-dldt/rv
+ xb=xa+hvap/(rv*ttp)
+ dldti=cvap-cice
+ xai=-dldti/rv
+ xbi=xai+hsub/(rv*ttp)
+ do k = kts, kte
+ do i = its, ite
+ tr=ttp/t(i,k)
+ qs(i,k,1)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
+ qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
+ qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
+ qs(i,k,1) = max(qs(i,k,1),qmin)
+ tr=ttp/t(i,k)
+ if(t(i,k).lt.ttp) then
+ qs(i,k,2)=psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
+ else
+ qs(i,k,2)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
+ endif
+ qs(i,k,2) = min(qs(i,k,2),0.99*p(i,k))
+ qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
+ qs(i,k,2) = max(qs(i,k,2),qmin)
+ enddo
+ enddo
+!
+!----------------------------------------------------------------
+! pcond: condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
+! if there exists additional water vapor condensated/if
+! evaporation of cloud water is not enough to remove subsaturation
+!
+ do k = kts, kte
+ do i = its, ite
+ work1(i,k,1) = conden(t(i,k),q(i,k),qs(i,k,1),xl(i,k),cpm(i,k))
+ work2(i,k) = qci(i,k,1)+work1(i,k,1)
+ pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k),0.)/dtcld)
+ if(qci(i,k,1).gt.0..and.work1(i,k,1).lt.0.) &
+ pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
+ q(i,k) = q(i,k)-pcond(i,k)*dtcld
+ qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
+ t(i,k) = t(i,k)+pcond(i,k)*xl(i,k)/cpm(i,k)*dtcld
+ enddo
+ enddo
+!
+!
+!----------------------------------------------------------------
+! padding for small values
+!
+ do k = kts, kte
+ do i = its, ite
+ if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
+ if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
+ enddo
+ enddo
+ enddo ! big loops
+ END SUBROUTINE wsm62d
+! ...................................................................
+ REAL FUNCTION rgmma(x)
+!-------------------------------------------------------------------
+ IMPLICIT NONE
+!-------------------------------------------------------------------
+! rgmma function: use infinite product form
+ REAL :: euler
+ PARAMETER (euler=0.577215664901532)
+ REAL :: x, y
+ INTEGER :: i
+ if(x.eq.1.)then
+ rgmma=0.
+ else
+ rgmma=x*exp(euler*x)
+ do i=1,10000
+ y=float(i)
+ rgmma=rgmma*(1.000+x/y)*exp(-x/y)
+ enddo
+ rgmma=1./rgmma
+ endif
+ END FUNCTION rgmma
+!
+!--------------------------------------------------------------------------
+ REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
+!--------------------------------------------------------------------------
+ IMPLICIT NONE
+!--------------------------------------------------------------------------
+ REAL t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti, &
+ xai,xbi,ttp,tr
+ INTEGER ice
+! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ ttp=t0c+0.01
+ dldt=cvap-cliq
+ xa=-dldt/rv
+ xb=xa+hvap/(rv*ttp)
+ dldti=cvap-cice
+ xai=-dldti/rv
+ xbi=xai+hsub/(rv*ttp)
+ tr=ttp/t
+ if(t.lt.ttp.and.ice.eq.1) then
+ fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
+ else
+ fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
+ endif
+! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ END FUNCTION fpvs
+!-------------------------------------------------------------------
+ SUBROUTINE wsm6init(den0,denr,dens,cl,cpv,allowed_to_read)
+!-------------------------------------------------------------------
+ IMPLICIT NONE
+!-------------------------------------------------------------------
+!.... constants which may not be tunable
+ REAL, INTENT(IN) :: den0,denr,dens,cl,cpv
+ LOGICAL, INTENT(IN) :: allowed_to_read
+!
+ pi = 4.*atan(1.)
+ xlv1 = cl-cpv
+!
+ qc0 = 4./3.*pi*denr*r0**3*xncr/den0 ! 0.419e-3 -- .61e-3
+ qck1 = .104*9.8*peaut/(xncr*denr)**(1./3.)/xmyu*den0**(4./3.) ! 7.03
+!
+ bvtr1 = 1.+bvtr
+ bvtr2 = 2.5+.5*bvtr
+ bvtr3 = 3.+bvtr
+ bvtr4 = 4.+bvtr
+ bvtr6 = 6.+bvtr
+ g1pbr = rgmma(bvtr1)
+ g3pbr = rgmma(bvtr3)
+ g4pbr = rgmma(bvtr4) ! 17.837825
+ g6pbr = rgmma(bvtr6)
+ g5pbro2 = rgmma(bvtr2) ! 1.8273
+ pvtr = avtr*g4pbr/6.
+ eacrr = 1.0
+ pacrr = pi*n0r*avtr*g3pbr*.25*eacrr
+ precr1 = 2.*pi*n0r*.78
+ precr2 = 2.*pi*n0r*.31*avtr**.5*g5pbro2
+ roqimax = 2.08e22*dimax**8
+!
+ bvts1 = 1.+bvts
+ bvts2 = 2.5+.5*bvts
+ bvts3 = 3.+bvts
+ bvts4 = 4.+bvts
+ g1pbs = rgmma(bvts1) !.8875
+ g3pbs = rgmma(bvts3)
+ g4pbs = rgmma(bvts4) ! 12.0786
+ g5pbso2 = rgmma(bvts2)
+ pvts = avts*g4pbs/6.
+ pacrs = pi*n0s*avts*g3pbs*.25
+ precs1 = 4.*n0s*.65
+ precs2 = 4.*n0s*.44*avts**.5*g5pbso2
+ pidn0r = pi*denr*n0r
+ pidn0s = pi*dens*n0s
+!
+ pacrc = pi*n0s*avts*g3pbs*.25*eacrc
+!
+ bvtg1 = 1.+bvtg
+ bvtg2 = 2.5+.5*bvtg
+ bvtg3 = 3.+bvtg
+ bvtg4 = 4.+bvtg
+ g1pbg = rgmma(bvtg1)
+ g3pbg = rgmma(bvtg3)
+ g4pbg = rgmma(bvtg4)
+ pacrg = pi*n0g*avtg*g3pbg*.25
+ g5pbgo2 = rgmma(bvtg2)
+ pvtg = avtg*g4pbg/6.
+ precg1 = 2.*pi*n0g*.78
+ precg2 = 2.*pi*n0g*.31*avtg**.5*g5pbgo2
+ pidn0g = pi*deng*n0g
+!
+ rslopermax = 1./lamdarmax
+ rslopesmax = 1./lamdasmax
+ rslopegmax = 1./lamdagmax
+ rsloperbmax = rslopermax ** bvtr
+ rslopesbmax = rslopesmax ** bvts
+ rslopegbmax = rslopegmax ** bvtg
+ rsloper2max = rslopermax * rslopermax
+ rslopes2max = rslopesmax * rslopesmax
+ rslopeg2max = rslopegmax * rslopegmax
+ rsloper3max = rsloper2max * rslopermax
+ rslopes3max = rslopes2max * rslopesmax
+ rslopeg3max = rslopeg2max * rslopegmax
+!
+ END SUBROUTINE wsm6init
+!------------------------------------------------------------------------------
+ subroutine slope_wsm6(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
+ vt,its,ite,kts,kte)
+ IMPLICIT NONE
+ INTEGER :: its,ite, jts,jte, kts,kte
+ REAL, DIMENSION( its:ite , kts:kte,3) :: &
+ qrs, &
+ rslope, &
+ rslopeb, &
+ rslope2, &
+ rslope3, &
+ vt
+ REAL, DIMENSION( its:ite , kts:kte) :: &
+ den, &
+ denfac, &
+ t
+ REAL, PARAMETER :: t0c = 273.15
+ REAL, DIMENSION( its:ite , kts:kte ) :: &
+ n0sfac
+ REAL :: lamdar, lamdas, lamdag, x, y, z, supcol
+ integer :: i, j, k
+!----------------------------------------------------------------
+! size distributions: (x=mixing ratio, y=air density):
+! valid for mixing ratio > 1.e-9 kg/kg.
+ lamdar(x,y)= sqrt(sqrt(pidn0r/(x*y))) ! (pidn0r/(x*y))**.25
+ lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y))) ! (pidn0s*z/(x*y))**.25
+ lamdag(x,y)= sqrt(sqrt(pidn0g/(x*y))) ! (pidn0g/(x*y))**.25
+!
+ do k = kts, kte
+ do i = its, ite
+ supcol = t0c-t(i,k)
+!---------------------------------------------------------------
+! n0s: Intercept parameter for snow [m-4] [HDC 6]
+!---------------------------------------------------------------
+ n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
+ if(qrs(i,k,1).le.qcrmin)then
+ rslope(i,k,1) = rslopermax
+ rslopeb(i,k,1) = rsloperbmax
+ rslope2(i,k,1) = rsloper2max
+ rslope3(i,k,1) = rsloper3max
+ else
+ rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k))
+ rslopeb(i,k,1) = rslope(i,k,1)**bvtr
+ rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
+ rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
+ endif
+ if(qrs(i,k,2).le.qcrmin)then
+ rslope(i,k,2) = rslopesmax
+ rslopeb(i,k,2) = rslopesbmax
+ rslope2(i,k,2) = rslopes2max
+ rslope3(i,k,2) = rslopes3max
+ else
+ rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
+ rslopeb(i,k,2) = rslope(i,k,2)**bvts
+ rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
+ rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
+ endif
+ if(qrs(i,k,3).le.qcrmin)then
+ rslope(i,k,3) = rslopegmax
+ rslopeb(i,k,3) = rslopegbmax
+ rslope2(i,k,3) = rslopeg2max
+ rslope3(i,k,3) = rslopeg3max
+ else
+ rslope(i,k,3) = 1./lamdag(qrs(i,k,3),den(i,k))
+ rslopeb(i,k,3) = rslope(i,k,3)**bvtg
+ rslope2(i,k,3) = rslope(i,k,3)*rslope(i,k,3)
+ rslope3(i,k,3) = rslope2(i,k,3)*rslope(i,k,3)
+ endif
+ vt(i,k,1) = pvtr*rslopeb(i,k,1)*denfac(i,k)
+ vt(i,k,2) = pvts*rslopeb(i,k,2)*denfac(i,k)
+ vt(i,k,3) = pvtg*rslopeb(i,k,3)*denfac(i,k)
+ if(qrs(i,k,1).le.0.0) vt(i,k,1) = 0.0
+ if(qrs(i,k,2).le.0.0) vt(i,k,2) = 0.0
+ if(qrs(i,k,3).le.0.0) vt(i,k,3) = 0.0
+ enddo
+ enddo
+ END subroutine slope_wsm6
+!-----------------------------------------------------------------------------
+ subroutine slope_rain(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
+ vt,its,ite,kts,kte)
+ IMPLICIT NONE
+ INTEGER :: its,ite, jts,jte, kts,kte
+ REAL, DIMENSION( its:ite , kts:kte) :: &
+ qrs, &
+ rslope, &
+ rslopeb, &
+ rslope2, &
+ rslope3, &
+ vt, &
+ den, &
+ denfac, &
+ t
+ REAL, PARAMETER :: t0c = 273.15
+ REAL, DIMENSION( its:ite , kts:kte ) :: &
+ n0sfac
+ REAL :: lamdar, x, y, z, supcol
+ integer :: i, j, k
+!----------------------------------------------------------------
+! size distributions: (x=mixing ratio, y=air density):
+! valid for mixing ratio > 1.e-9 kg/kg.
+ lamdar(x,y)= sqrt(sqrt(pidn0r/(x*y))) ! (pidn0r/(x*y))**.25
+!
+ do k = kts, kte
+ do i = its, ite
+ if(qrs(i,k).le.qcrmin)then
+ rslope(i,k) = rslopermax
+ rslopeb(i,k) = rsloperbmax
+ rslope2(i,k) = rsloper2max
+ rslope3(i,k) = rsloper3max
+ else
+ rslope(i,k) = 1./lamdar(qrs(i,k),den(i,k))
+ rslopeb(i,k) = rslope(i,k)**bvtr
+ rslope2(i,k) = rslope(i,k)*rslope(i,k)
+ rslope3(i,k) = rslope2(i,k)*rslope(i,k)
+ endif
+ vt(i,k) = pvtr*rslopeb(i,k)*denfac(i,k)
+ if(qrs(i,k).le.0.0) vt(i,k) = 0.0
+ enddo
+ enddo
+ END subroutine slope_rain
+!------------------------------------------------------------------------------
+ subroutine slope_snow(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
+ vt,its,ite,kts,kte)
+ IMPLICIT NONE
+ INTEGER :: its,ite, jts,jte, kts,kte
+ REAL, DIMENSION( its:ite , kts:kte) :: &
+ qrs, &
+ rslope, &
+ rslopeb, &
+ rslope2, &
+ rslope3, &
+ vt, &
+ den, &
+ denfac, &
+ t
+ REAL, PARAMETER :: t0c = 273.15
+ REAL, DIMENSION( its:ite , kts:kte ) :: &
+ n0sfac
+ REAL :: lamdas, x, y, z, supcol
+ integer :: i, j, k
+!----------------------------------------------------------------
+! size distributions: (x=mixing ratio, y=air density):
+! valid for mixing ratio > 1.e-9 kg/kg.
+ lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y))) ! (pidn0s*z/(x*y))**.25
+!
+ do k = kts, kte
+ do i = its, ite
+ supcol = t0c-t(i,k)
+!---------------------------------------------------------------
+! n0s: Intercept parameter for snow [m-4] [HDC 6]
+!---------------------------------------------------------------
+ n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
+ if(qrs(i,k).le.qcrmin)then
+ rslope(i,k) = rslopesmax
+ rslopeb(i,k) = rslopesbmax
+ rslope2(i,k) = rslopes2max
+ rslope3(i,k) = rslopes3max
+ else
+ rslope(i,k) = 1./lamdas(qrs(i,k),den(i,k),n0sfac(i,k))
+ rslopeb(i,k) = rslope(i,k)**bvts
+ rslope2(i,k) = rslope(i,k)*rslope(i,k)
+ rslope3(i,k) = rslope2(i,k)*rslope(i,k)
+ endif
+ vt(i,k) = pvts*rslopeb(i,k)*denfac(i,k)
+ if(qrs(i,k).le.0.0) vt(i,k) = 0.0
+ enddo
+ enddo
+ END subroutine slope_snow
+!----------------------------------------------------------------------------------
+ subroutine slope_graup(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
+ vt,its,ite,kts,kte)
+ IMPLICIT NONE
+ INTEGER :: its,ite, jts,jte, kts,kte
+ REAL, DIMENSION( its:ite , kts:kte) :: &
+ qrs, &
+ rslope, &
+ rslopeb, &
+ rslope2, &
+ rslope3, &
+ vt, &
+ den, &
+ denfac, &
+ t
+ REAL, PARAMETER :: t0c = 273.15
+ REAL, DIMENSION( its:ite , kts:kte ) :: &
+ n0sfac
+ REAL :: lamdag, x, y, z, supcol
+ integer :: i, j, k
+!----------------------------------------------------------------
+! size distributions: (x=mixing ratio, y=air density):
+! valid for mixing ratio > 1.e-9 kg/kg.
+ lamdag(x,y)= sqrt(sqrt(pidn0g/(x*y))) ! (pidn0g/(x*y))**.25
+!
+ do k = kts, kte
+ do i = its, ite
+!---------------------------------------------------------------
+! n0s: Intercept parameter for snow [m-4] [HDC 6]
+!---------------------------------------------------------------
+ if(qrs(i,k).le.qcrmin)then
+ rslope(i,k) = rslopegmax
+ rslopeb(i,k) = rslopegbmax
+ rslope2(i,k) = rslopeg2max
+ rslope3(i,k) = rslopeg3max
+ else
+ rslope(i,k) = 1./lamdag(qrs(i,k),den(i,k))
+ rslopeb(i,k) = rslope(i,k)**bvtg
+ rslope2(i,k) = rslope(i,k)*rslope(i,k)
+ rslope3(i,k) = rslope2(i,k)*rslope(i,k)
+ endif
+ vt(i,k) = pvtg*rslopeb(i,k)*denfac(i,k)
+ if(qrs(i,k).le.0.0) vt(i,k) = 0.0
+ enddo
+ enddo
+ END subroutine slope_graup
+!---------------------------------------------------------------------------------
+!-------------------------------------------------------------------
+ SUBROUTINE nislfv_rain_plm(im,km,denl,denfacl,tkl,dzl,wwl,rql,precip,dt,id,iter)
+!-------------------------------------------------------------------
+!
+! for non-iteration semi-Lagrangain forward advection for cloud
+! with mass conservation and positive definite advection
+! 2nd order interpolation with monotonic piecewise linear method
+! this routine is under assumption of decfl < 1 for semi_Lagrangian
+!
+! dzl depth of model layer in meter
+! wwl terminal velocity at model layer m/s
+! rql cloud density*mixing ration
+! precip precipitation
+! dt time step
+! id kind of precip: 0 test case; 1 raindrop
+! iter how many time to guess mean terminal velocity: 0 pure forward.
+! 0 : use departure wind for advection
+! 1 : use mean wind for advection
+! > 1 : use mean wind after iter-1 iterations
+!
+! author: hann-ming henry juang <henry.juang@noaa.gov>
+! implemented by song-you hong
+!
+ implicit none
+ integer im,km,id
+ real dt
+ real dzl(im,km),wwl(im,km),rql(im,km),precip(im)
+ real denl(im,km),denfacl(im,km),tkl(im,km)
+!
+ integer i,k,n,m,kk,kb,kt,iter
+ real tl,tl2,qql,dql,qqd
+ real th,th2,qqh,dqh
+ real zsum,qsum,dim,dip,c1,con1,fa1,fa2
+ real allold, allnew, zz, dzamin, cflmax, decfl
+ real dz(km), ww(km), qq(km), wd(km), wa(km), was(km)
+ real den(km), denfac(km), tk(km)
+ real wi(km+1), zi(km+1), za(km+1)
+ real qn(km), qr(km),tmp(km),tmp1(km),tmp2(km),tmp3(km)
+ real dza(km+1), qa(km+1), qmi(km+1), qpi(km+1)
+!
+ precip(:) = 0.0
+!
+ i_loop : do i=1,im
+! -----------------------------------
+ dz(:) = dzl(i,:)
+ qq(:) = rql(i,:)
+ ww(:) = wwl(i,:)
+ den(:) = denl(i,:)
+ denfac(:) = denfacl(i,:)
+ tk(:) = tkl(i,:)
+! skip for no precipitation for all layers
+ allold = 0.0
+ do k=1,km
+ allold = allold + qq(k)
+ enddo
+ if(allold.le.0.0) then
+ cycle i_loop
+ endif
+!
+! compute interface values
+ zi(1)=0.0
+ do k=1,km
+ zi(k+1) = zi(k)+dz(k)
+ enddo
+!
+! save departure wind
+ wd(:) = ww(:)
+ n=1
+ 100 continue
+! plm is 2nd order, we can use 2nd order wi or 3rd order wi
+! 2nd order interpolation to get wi
+ wi(1) = ww(1)
+ wi(km+1) = ww(km)
+ do k=2,km
+ wi(k) = (ww(k)*dz(k-1)+ww(k-1)*dz(k))/(dz(k-1)+dz(k))
+ enddo
+! 3rd order interpolation to get wi
+ fa1 = 9./16.
+ fa2 = 1./16.
+ wi(1) = ww(1)
+ wi(2) = 0.5*(ww(2)+ww(1))
+ do k=3,km-1
+ wi(k) = fa1*(ww(k)+ww(k-1))-fa2*(ww(k+1)+ww(k-2))
+ enddo
+ wi(km) = 0.5*(ww(km)+ww(km-1))
+ wi(km+1) = ww(km)
+!
+! terminate of top of raingroup
+ do k=2,km
+ if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
+ enddo
+!
+! diffusivity of wi
+ con1 = 0.05
+ do k=km,1,-1
+ decfl = (wi(k+1)-wi(k))*dt/dz(k)
+ if( decfl .gt. con1 ) then
+ wi(k) = wi(k+1) - con1*dz(k)/dt
+ endif
+ enddo
+! compute arrival point
+ do k=1,km+1
+ za(k) = zi(k) - wi(k)*dt
+ enddo
+!
+ do k=1,km
+ dza(k) = za(k+1)-za(k)
+ enddo
+ dza(km+1) = zi(km+1) - za(km+1)
+!
+! computer deformation at arrival point
+ do k=1,km
+ qa(k) = qq(k)*dz(k)/dza(k)
+ qr(k) = qa(k)/den(k)
+ enddo
+ qa(km+1) = 0.0
+! call maxmin(km,1,qa,' arrival points ')
+!
+! compute arrival terminal velocity, and estimate mean terminal velocity
+! then back to use mean terminal velocity
+ if( n.le.iter ) then
+ call slope_rain(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
+ if( n.ge.2 ) wa(1:km)=0.5*(wa(1:km)+was(1:km))
+ do k=1,km
+!#ifdef DEBUG
+! print*,' slope_wsm3 ',qr(k)*1000.,den(k),denfac(k),tk(k),tmp(k),tmp1(k),tmp2(k),ww(k),wa(k)
+!#endif
+! mean wind is average of departure and new arrival winds
+ ww(k) = 0.5* ( wd(k)+wa(k) )
+ enddo
+ was(:) = wa(:)
+ n=n+1
+ go to 100
+ endif
+!
+! estimate values at arrival cell interface with monotone
+ do k=2,km
+ dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k))
+ dim=(qa(k)-qa(k-1))/(dza(k-1)+dza(k))
+ if( dip*dim.le.0.0 ) then
+ qmi(k)=qa(k)
+ qpi(k)=qa(k)
+ else
+ qpi(k)=qa(k)+0.5*(dip+dim)*dza(k)
+ qmi(k)=2.0*qa(k)-qpi(k)
+ if( qpi(k).lt.0.0 .or. qmi(k).lt.0.0 ) then
+ qpi(k) = qa(k)
+ qmi(k) = qa(k)
+ endif
+ endif
+ enddo
+ qpi(1)=qa(1)
+ qmi(1)=qa(1)
+ qmi(km+1)=qa(km+1)
+ qpi(km+1)=qa(km+1)
+!
+! interpolation to regular point
+ qn = 0.0
+ kb=1
+ kt=1
+ intp : do k=1,km
+ kb=max(kb-1,1)
+ kt=max(kt-1,1)
+! find kb and kt
+ if( zi(k).ge.za(km+1) ) then
+ exit intp
+ else
+ find_kb : do kk=kb,km
+ if( zi(k).le.za(kk+1) ) then
+ kb = kk
+ exit find_kb
+ else
+ cycle find_kb
+ endif
+ enddo find_kb
+ find_kt : do kk=kt,km
+ if( zi(k+1).le.za(kk) ) then
+ kt = kk
+ exit find_kt
+ else
+ cycle find_kt
+ endif
+ enddo find_kt
+ kt = kt - 1
+! compute q with piecewise constant method
+ if( kt.eq.kb ) then
+ tl=(zi(k)-za(kb))/dza(kb)
+ th=(zi(k+1)-za(kb))/dza(kb)
+ tl2=tl*tl
+ th2=th*th
+ qqd=0.5*(qpi(kb)-qmi(kb))
+ qqh=qqd*th2+qmi(kb)*th
+ qql=qqd*tl2+qmi(kb)*tl
+ qn(k) = (qqh-qql)/(th-tl)
+ else if( kt.gt.kb ) then
+ tl=(zi(k)-za(kb))/dza(kb)
+ tl2=tl*tl
+ qqd=0.5*(qpi(kb)-qmi(kb))
+ qql=qqd*tl2+qmi(kb)*tl
+ dql = qa(kb)-qql
+ zsum = (1.-tl)*dza(kb)
+ qsum = dql*dza(kb)
+ if( kt-kb.gt.1 ) then
+ do m=kb+1,kt-1
+ zsum = zsum + dza(m)
+ qsum = qsum + qa(m) * dza(m)
+ enddo
+ endif
+ th=(zi(k+1)-za(kt))/dza(kt)
+ th2=th*th
+ qqd=0.5*(qpi(kt)-qmi(kt))
+ dqh=qqd*th2+qmi(kt)*th
+ zsum = zsum + th*dza(kt)
+ qsum = qsum + dqh*dza(kt)
+ qn(k) = qsum/zsum
+ endif
+ cycle intp
+ endif
+!
+ enddo intp
+!
+! rain out
+ sum_precip: do k=1,km
+ if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then
+ precip(i) = precip(i) + qa(k)*dza(k)
+ cycle sum_precip
+ else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then
+ precip(i) = precip(i) + qa(k)*(0.0-za(k))
+ exit sum_precip
+ endif
+ exit sum_precip
+ enddo sum_precip
+!
+! replace the new values
+ rql(i,:) = qn(:)
+!
+! ----------------------------------
+ enddo i_loop
+!
+ END SUBROUTINE nislfv_rain_plm
+!-------------------------------------------------------------------
+ SUBROUTINE nislfv_rain_plm6(im,km,denl,denfacl,tkl,dzl,wwl,rql,rql2, precip1, precip2,dt,id,iter)
+!-------------------------------------------------------------------
+!
+! for non-iteration semi-Lagrangain forward advection for cloud
+! with mass conservation and positive definite advection
+! 2nd order interpolation with monotonic piecewise linear method
+! this routine is under assumption of decfl < 1 for semi_Lagrangian
+!
+! dzl depth of model layer in meter
+! wwl terminal velocity at model layer m/s
+! rql cloud density*mixing ration
+! precip precipitation
+! dt time step
+! id kind of precip: 0 test case; 1 raindrop
+! iter how many time to guess mean terminal velocity: 0 pure forward.
+! 0 : use departure wind for advection
+! 1 : use mean wind for advection
+! > 1 : use mean wind after iter-1 iterations
+!
+! author: hann-ming henry juang <henry.juang@noaa.gov>
+! implemented by song-you hong
+!
+ implicit none
+ integer im,km,id
+ real dt
+ real dzl(im,km),wwl(im,km),rql(im,km),rql2(im,km),precip(im),precip1(im),precip2(im)
+ real denl(im,km),denfacl(im,km),tkl(im,km)
+!
+ integer i,k,n,m,kk,kb,kt,iter,ist
+ real tl,tl2,qql,dql,qqd
+ real th,th2,qqh,dqh
+ real zsum,qsum,dim,dip,c1,con1,fa1,fa2
+ real allold, allnew, zz, dzamin, cflmax, decfl
+ real dz(km), ww(km), qq(km), qq2(km), wd(km), wa(km), wa2(km), was(km)
+ real den(km), denfac(km), tk(km)
+ real wi(km+1), zi(km+1), za(km+1)
+ real qn(km), qr(km),qr2(km),tmp(km),tmp1(km),tmp2(km),tmp3(km)
+ real dza(km+1), qa(km+1), qa2(km+1),qmi(km+1), qpi(km+1)
+!
+ precip(:) = 0.0
+ precip1(:) = 0.0
+ precip2(:) = 0.0
+!
+ i_loop : do i=1,im
+! -----------------------------------
+ dz(:) = dzl(i,:)
+ qq(:) = rql(i,:)
+ qq2(:) = rql2(i,:)
+ ww(:) = wwl(i,:)
+ den(:) = denl(i,:)
+ denfac(:) = denfacl(i,:)
+ tk(:) = tkl(i,:)
+! skip for no precipitation for all layers
+ allold = 0.0
+ do k=1,km
+ allold = allold + qq(k)
+ enddo
+ if(allold.le.0.0) then
+ cycle i_loop
+ endif
+!
+! compute interface values
+ zi(1)=0.0
+ do k=1,km
+ zi(k+1) = zi(k)+dz(k)
+ enddo
+!
+! save departure wind
+ wd(:) = ww(:)
+ n=1
+ 100 continue
+! plm is 2nd order, we can use 2nd order wi or 3rd order wi
+! 2nd order interpolation to get wi
+ wi(1) = ww(1)
+ wi(km+1) = ww(km)
+ do k=2,km
+ wi(k) = (ww(k)*dz(k-1)+ww(k-1)*dz(k))/(dz(k-1)+dz(k))
+ enddo
+! 3rd order interpolation to get wi
+ fa1 = 9./16.
+ fa2 = 1./16.
+ wi(1) = ww(1)
+ wi(2) = 0.5*(ww(2)+ww(1))
+ do k=3,km-1
+ wi(k) = fa1*(ww(k)+ww(k-1))-fa2*(ww(k+1)+ww(k-2))
+ enddo
+ wi(km) = 0.5*(ww(km)+ww(km-1))
+ wi(km+1) = ww(km)
+!
+! terminate of top of raingroup
+ do k=2,km
+ if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
+ enddo
+!
+! diffusivity of wi
+ con1 = 0.05
+ do k=km,1,-1
+ decfl = (wi(k+1)-wi(k))*dt/dz(k)
+ if( decfl .gt. con1 ) then
+ wi(k) = wi(k+1) - con1*dz(k)/dt
+ endif
+ enddo
+! compute arrival point
+ do k=1,km+1
+ za(k) = zi(k) - wi(k)*dt
+ enddo
+!
+ do k=1,km
+ dza(k) = za(k+1)-za(k)
+ enddo
+ dza(km+1) = zi(km+1) - za(km+1)
+!
+! computer deformation at arrival point
+ do k=1,km
+ qa(k) = qq(k)*dz(k)/dza(k)
+ qa2(k) = qq2(k)*dz(k)/dza(k)
+ qr(k) = qa(k)/den(k)
+ qr2(k) = qa2(k)/den(k)
+ enddo
+ qa(km+1) = 0.0
+ qa2(km+1) = 0.0
+! call maxmin(km,1,qa,' arrival points ')
+!
+! compute arrival terminal velocity, and estimate mean terminal velocity
+! then back to use mean terminal velocity
+ if( n.le.iter ) then
+ call slope_snow(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
+ call slope_graup(qr2,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa2,1,1,1,km)
+ do k = 1, km
+ tmp(k) = max((qr(k)+qr2(k)), 1.E-15)
+ IF ( tmp(k) .gt. 1.e-15 ) THEN
+ wa(k) = (wa(k)*qr(k) + wa2(k)*qr2(k))/tmp(k)
+ ELSE
+ wa(k) = 0.
+ ENDIF
+ enddo
+ if( n.ge.2 ) wa(1:km)=0.5*(wa(1:km)+was(1:km))
+ do k=1,km
+!#ifdef DEBUG
+! print*,' slope_wsm3 ',qr(k)*1000.,den(k),denfac(k),tk(k),tmp(k),tmp1(k),tmp2(k), &
+! ww(k),wa(k)
+!#endif
+! mean wind is average of departure and new arrival winds
+ ww(k) = 0.5* ( wd(k)+wa(k) )
+ enddo
+ was(:) = wa(:)
+ n=n+1
+ go to 100
+ endif
+ ist_loop : do ist = 1, 2
+ if (ist.eq.2) then
+ qa(:) = qa2(:)
+ endif
+!
+ precip(i) = 0.
+!
+! estimate values at arrival cell interface with monotone
+ do k=2,km
+ dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k))
+ dim=(qa(k)-qa(k-1))/(dza(k-1)+dza(k))
+ if( dip*dim.le.0.0 ) then
+ qmi(k)=qa(k)
+ qpi(k)=qa(k)
+ else
+ qpi(k)=qa(k)+0.5*(dip+dim)*dza(k)
+ qmi(k)=2.0*qa(k)-qpi(k)
+ if( qpi(k).lt.0.0 .or. qmi(k).lt.0.0 ) then
+ qpi(k) = qa(k)
+ qmi(k) = qa(k)
+ endif
+ endif
+ enddo
+ qpi(1)=qa(1)
+ qmi(1)=qa(1)
+ qmi(km+1)=qa(km+1)
+ qpi(km+1)=qa(km+1)
+!
+! interpolation to regular point
+ qn = 0.0
+ kb=1
+ kt=1
+ intp : do k=1,km
+ kb=max(kb-1,1)
+ kt=max(kt-1,1)
+! find kb and kt
+ if( zi(k).ge.za(km+1) ) then
+ exit intp
+ else
+ find_kb : do kk=kb,km
+ if( zi(k).le.za(kk+1) ) then
+ kb = kk
+ exit find_kb
+ else
+ cycle find_kb
+ endif
+ enddo find_kb
+ find_kt : do kk=kt,km
+ if( zi(k+1).le.za(kk) ) then
+ kt = kk
+ exit find_kt
+ else
+ cycle find_kt
+ endif
+ enddo find_kt
+ kt = kt - 1
+! compute q with piecewise constant method
+ if( kt.eq.kb ) then
+ tl=(zi(k)-za(kb))/dza(kb)
+ th=(zi(k+1)-za(kb))/dza(kb)
+ tl2=tl*tl
+ th2=th*th
+ qqd=0.5*(qpi(kb)-qmi(kb))
+ qqh=qqd*th2+qmi(kb)*th
+ qql=qqd*tl2+qmi(kb)*tl
+ qn(k) = (qqh-qql)/(th-tl)
+ else if( kt.gt.kb ) then
+ tl=(zi(k)-za(kb))/dza(kb)
+ tl2=tl*tl
+ qqd=0.5*(qpi(kb)-qmi(kb))
+ qql=qqd*tl2+qmi(kb)*tl
+ dql = qa(kb)-qql
+ zsum = (1.-tl)*dza(kb)
+ qsum = dql*dza(kb)
+ if( kt-kb.gt.1 ) then
+ do m=kb+1,kt-1
+ zsum = zsum + dza(m)
+ qsum = qsum + qa(m) * dza(m)
+ enddo
+ endif
+ th=(zi(k+1)-za(kt))/dza(kt)
+ th2=th*th
+ qqd=0.5*(qpi(kt)-qmi(kt))
+ dqh=qqd*th2+qmi(kt)*th
+ zsum = zsum + th*dza(kt)
+ qsum = qsum + dqh*dza(kt)
+ qn(k) = qsum/zsum
+ endif
+ cycle intp
+ endif
+!
+ enddo intp
+!
+! rain out
+ sum_precip: do k=1,km
+ if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then
+ precip(i) = precip(i) + qa(k)*dza(k)
+ cycle sum_precip
+ else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then
+ precip(i) = precip(i) + qa(k)*(0.0-za(k))
+ exit sum_precip
+ endif
+ exit sum_precip
+ enddo sum_precip
+!
+! replace the new values
+ if(ist.eq.1) then
+ rql(i,:) = qn(:)
+ precip1(i) = precip(i)
+ else
+ rql2(i,:) = qn(:)
+ precip2(i) = precip(i)
+ endif
+ enddo ist_loop
+!
+! ----------------------------------
+ enddo i_loop
+!
+ END SUBROUTINE nislfv_rain_plm6
+END MODULE module_mp_wsm6
</font>
</pre>