c **********************************************************************
c *   An improved Mellor-Yamada turbulence closure model               *
c *                                                                    *
c *                                   Aug/2005  M. Nakanishi (N.D.A)   *
c *                        Modified:  Dec/2005  M. Nakanishi (N.D.A)   *
c *                    Minor change:  Jan/2009  M. Nakanishi (N.D.A)   *
c *                                             naka@nda.ac.jp         *
c *                                                                    *
c *   Contents:                                                        *
c *     1. mym_initialize  (to be called once initially)               *
c *        gives the closure constants and initializes the turbulent   *
c *        quantities.                                                 *
c *    (2) mym_level2      (called in the other subroutines)           *
c *        calculates the stability functions at Level 2.              *
c *    (3) mym_length      (called in the other subroutines)           *
c *        calculates the master length scale.                         *
c *     4. mym_turbulence                                              *
c *        calculates the vertical diffusivity coefficients and the    *
c *        production terms for the turbulent quantities.              *
c *     5. mym_predict                                                 *
c *        predicts the turbulent quantities at the next step.         *
c *     6. mym_condensation                                            *
c *        determines the liquid water content and the cloud fraction  *
c *        diagnostically.                                             *
c *                                                                    *
c *             call mym_initialize                                    *
c *                  |                                                 *
c *                  |<----------------+                               *
c *                  |                 |                               *
c *             call mym_condensation  |                               *
c *             call mym_turbulence    |                               *
c *             call mym_predict       |                               *
c *                  |                 |                               *
c *                  |-----------------+                               *
c *                  |                                                 *
c *                 end                                                *
c *                                                                    *
c *   Variables worthy of special mention:                             *
c *     tref   : Reference temperature                                 *
c *     tl     : Liquid water potential temperature minus tref         *
c *     qw     : Total water (water vapor+liquid water) content        *
c *     ql     : Liquid water content                                  *
c *     vt, vq : Functions for computing the buoyancy flux             *
c *                                                                    *
c *     If the water contents are unnecessary, e.g., in the case of    *
c *     ocean models, tl is the potential temperature and qw, ql, vt   *
c *     and vq are all zero.                                           *
c *                                                                    *
c *   Grid arrangement:                                                *
c *             k+1 +---------+                                        *
c *                 |         |     i = 1 - nx                         *
c *             (k) |    *    |     j = 1 - ny                         *
c *                 |         |     k = 1 - nz                         *
c *              k  +---------+                                        *
c *                 i   (i)  i+1                                       *
c *                                                                    *
c *     All the predicted variables are defined at the center (*) of   *
c *     the grid boxes. The diffusivity coefficients are, however,     *
c *     defined on the walls of the grid boxes.                        *
c *     # Upper boundary values are given at k=nz+1.                   *
c *                                                                    *
c *   References:                                                      *
c *     1. Nakanishi, M., 2001:                                        *
c *        Boundary-Layer Meteor., 99, 349-378.                        *
c *     2. Nakanishi, M. and H. Niino, 2004:                           *
c *        Boundary-Layer Meteor., 112, 1-31.                          *
c *     3. Nakanishi, M. and H. Niino, 2006:                           *
c *        Boundary-Layer Meteor., 119, 397-407.                       *
c **********************************************************************
c
c
c     SUBROUTINE  mym_initialize:
c
c     Input variables:
c       iniflag         : <>0; turbulent quantities will be initialized
c                         = 0; turbulent quantities have been already
c                              given, i.e., they will not be initialized
c       mx, my          : Maximum numbers of grid boxes
c                         in the x and y directions, respectively
c       nx, ny, nz      : Numbers of the actual grid boxes
c                         in the x, y and z directions, respectively
c       tref            : Reference temperature                      (K)
c       dz(nz+1)        : Vertical grid spacings                     (m)
c                         # dz(nz+1)=dz(nz)
c       zw(nz+1)        : Heights of the walls of the grid boxes     (m)
c                         # zw(1)=0.0 and zw(k)=zw(k-1)+dz(k-1)
c       h(mx,ny)        : G^(1/2) in the terrain-following coordinate
c                         # h=1-zg/zt, where zg is the height of the
c                           terrain and zt the top of the model domain
c       pi0(mx,my,nz+1) : Exner function at zw*h+zg             (J/kg K)
c                         defined by c_p*( p_basic/1000hPa )^kappa
c                         This is usually computed by integrating
c                         d(pi0)/dz = -h*g/tref.
c       rmo(mx,ny)      : Inverse of the Obukhov length         (m^(-1))
c       flt, flq(mx,ny) : Turbulent fluxes of sensible and latent heat,
c                         respectively, e.g., flt=-u_*Theta_* (K m/s)
c       ust(mx,ny)      : Friction velocity                        (m/s)
c       pmz(mx,ny)      : phi_m-zeta at z1*h+z0, where z1 (=0.5*dz(1))
c                         is the first grid point above the surafce, z0
c                         the roughness length and zeta=(z1*h+z0)*rmo
c       phh(mx,ny)      : phi_h at z1*h+z0
c       u, v(mx,my,nz+1): Components of the horizontal wind        (m/s)
c       tl(mx,my,nz+1)  : Liquid water potential temperature minus tref
c                                                                    (K)
c       qw(mx,my,nz+1)  : Total water content Q_w                (kg/kg)
c
c     Output variables:
c       ql(mx,my,nz+1)  : Liquid water content                   (kg/kg)
c       v?(mx,my,nz+1)  : Functions for computing the buoyancy flux
c       qke(mx,my,nz+1) : Twice the turbulent kinetic energy q^2
c                                                              (m^2/s^2)
c       tsq(mx,my,nz+1) : Variance of Theta_l                      (K^2)
c       qsq(mx,my,nz+1) : Variance of Q_w
c       cov(mx,my,nz+1) : Covariance of Theta_l and Q_w              (K)
c       el(mx,my,nz+1)  : Master length scale L                      (m)
c                         defined on the walls of the grid boxes
c       bsh             : no longer used
c       via common      : Closure constants
c
c     Work arrays:        see subroutine mym_level2
c       pd?(mx,my,nz+1) : Half of the production terms at Level 2
c                         defined on the walls of the grid boxes
c       qkw(mx,my,nz+1) : q on the walls of the grid boxes         (m/s)
c
c     # As to dtl, ...gh, see subroutine mym_turbulence.
c
      subroutine  mym_initialize (
     i            iniflag, mx, my, nx, ny, nz,
     i            tref, dz, zw, h, pi0,
     i            rmo, flt, flq, ust, pmz, phh,
     i            u, v, tl, qw,
     o            ql, vt, vq,
     o            qke, tsq, qsq, cov,
     o            el,
     o            bsh,
     w            pdk, pdt, pdq, pdc,
     w            qkw,
     w            dtl, dqw, dtv, gm, gh, sm, sh )
c
      dimension  dz(*), zw(*), h(mx,*), pi0(mx,my,*)
      dimension  rmo(mx,*), flt(mx,*), flq(mx,*), ust(mx,*)
      dimension  pmz(mx,*), phh(mx,*)
      dimension  u  (mx,my,*), v  (mx,my,*), tl (mx,my,*), qw (mx,my,*)
      dimension  ql (mx,my,*), vt (mx,my,*), vq (mx,my,*)
      dimension  qke(mx,my,*), tsq(mx,my,*), qsq(mx,my,*), cov(mx,my,*)
      dimension  el (mx,my,*)
      dimension  pdk(mx,my,*), pdt(mx,my,*), pdq(mx,my,*), pdc(mx,my,*)
      dimension  qkw(mx,my,*)
      dimension  dtl(mx,my,*), dqw(mx,my,*), dtv(mx,my,*)
      dimension  gm (mx,my,*), gh (mx,my,*), sm (mx,my,*), sh (mx,my,*)
c
      common  /constants/ g1, g2, a1, a2, b1, b2, c1, c2, c3, c4, c5
c
      g1  =  0.235
      b1  = 24.0
      b2  = 15.0
c  Modified: Dec/22/2005, from here, (minor impact)
c?    c2  =  0.7
c?    c3  =  0.323
      c2  =  0.75
      c3  =  0.352
c  Modified: Dec/22/2005, up to here
      c4  =  0.0
      c5  =  0.2
      pr  =  0.74
      vk  =  0.4
c
      a1  = b1*( 1.0-3.0*g1 )/6.0
      c1  = g1 -1.0/( 3.0*a1*b1**(1.0/3.0) )
      a2  = a1*( g1-c1 )/( g1*pr )
      g2  = b2/b1*( 1.0-c3 ) +2.0*a1/b1*( 3.0-2.0*c2 )
c
      if ( iniflag .eq. 0 ) return
c
c   **  At first ql, vt and vq are set to zero.  **
      do k = 1,nz+1
      do j = 1,my
      do i = 1,mx
        ql(i,j,k) = 0.0
        vt(i,j,k) = 0.0
        vq(i,j,k) = 0.0
      end do
      end do
      end do
c
      call mym_level2 (
     i            mx, my, nx, ny, nz,
     i            tref, dz, h, pi0,
     i            u, v, tl, qw,
     i            ql, vt, vq,
     o            dtl, dqw, dtv, gm, gh, sm, sh )
c
c   **  Preliminary setting  **
      do j = 1,ny
      do i = 1,mx
        el (i,j,1) = 0.0
        qke(i,j,1) = ust(i,j)**2 * ( b1*pmz(i,j) )**(2.0/3.0)
c
        phm        = phh(i,j)*b2 / ( b1*pmz(i,j) )**(1.0/3.0)
        tsq(i,j,1) = phm*( flt(i,j)/ust(i,j) )**2
        qsq(i,j,1) = phm*( flq(i,j)/ust(i,j) )**2
        cov(i,j,1) = phm*( flt(i,j)/ust(i,j) )*( flq(i,j)/ust(i,j) )
      end do
      end do
c
      do k = 2,nz+1
        vkz = vk*zw(k)
      do j = 1,my
      do i = 1,mx
        el (i,j,k) = vkz*h(i,j)/( 1.0 + vkz*h(i,j)/100.0 )
        qke(i,j,k) = 0.0
c
        tsq(i,j,k) = 0.0
        qsq(i,j,k) = 0.0
        cov(i,j,k) = 0.0
      end do
      end do
      end do
c
c   **  Initialization with an iterative manner          **
c   **  lmax is the iteration count. This is arbitrary.  **
      lmax = nz +3
c
      do l = 1,lmax
c
        call mym_length (
     i            mx, my, nx, ny, nz,
     i            tref, dz, zw, h,
     i            rmo, flt, flq,
     i            vt(1,1,1), vq(1,1,1),
     i            qke,
     i            dtv,
     o            el,
     o            qkw,
     w            pdk(1,1,nz), pdt(1,1,nz) )
c
      do k = 2,nz+1
      do j = 1,my
      do i = 1,mx
        elq = el(i,j,k)*qkw(i,j,k)
        pdk(i,j,k) = elq*( sm(i,j,k)*gm (i,j,k)+sh(i,j,k)*gh (i,j,k) )
        pdt(i,j,k) = elq*  sh(i,j,k)*dtl(i,j,k)**2
        pdq(i,j,k) = elq*  sh(i,j,k)*dqw(i,j,k)**2
        pdc(i,j,k) = elq*  sh(i,j,k)*dtl(i,j,k)*dqw(i,j,k)
      end do
      end do
      end do
c
c   **  Strictly, vkz*h(i,j) -> vk*( 0.5*dz(1)*h(i,j)+z0 )  **
      vkz = vk*0.5*dz(1)
c
      do j = 1,ny
      do i = 1,mx
        elv = 0.5*( el(i,j,2)+el(i,j,1) ) / ( vkz*h(i,j) )
        qke(i,j,1) = ust(i,j)**2 * ( b1*pmz(i,j)*elv    )**(2.0/3.0)
c
        phm        = phh(i,j)*b2 / ( b1*pmz(i,j)/elv**2 )**(1.0/3.0)
        tsq(i,j,1) = phm*( flt(i,j)/ust(i,j) )**2
        qsq(i,j,1) = phm*( flq(i,j)/ust(i,j) )**2
        cov(i,j,1) = phm*( flt(i,j)/ust(i,j) )*( flq(i,j)/ust(i,j) )
      end do
      end do
c
      do k = 2,nz
      do j = 1,my
      do i = 1,mx
        b1l = b1*0.25*( el(i,j,k+1)+el(i,j,k) )
        qke(i,j,k) = ( b1l*( pdk(i,j,k+1)+pdk(i,j,k) ) )**(2.0/3.0)
c
        if ( qke(i,j,k) .le. 0.0 ) then
          b2l = 0.0
        else
          b2l = b2*( b1l/b1 ) / sqrt( qke(i,j,k) )
        end if
c
        tsq(i,j,k) = b2l*( pdt(i,j,k+1)+pdt(i,j,k) )
        qsq(i,j,k) = b2l*( pdq(i,j,k+1)+pdq(i,j,k) )
        cov(i,j,k) = b2l*( pdc(i,j,k+1)+pdc(i,j,k) )
      end do
      end do
      end do
c
      end do
c
      return
      end
c
c     SUBROUTINE  mym_level2:
c
c     Input variables:    see subroutine mym_initialize
c
c     Output variables:
c       dtl(mx,my,nz+1) : Vertical gradient of Theta_l             (K/m)
c       dqw(mx,my,nz+1) : Vertical gradient of Q_w
c       dtv(mx,my,nz+1) : Vertical gradient of Theta_V             (K/m)
c       gm (mx,my,nz+1) : G_M divided by L^2/q^2                (s^(-2))
c       gh (mx,my,nz+1) : G_H divided by L^2/q^2                (s^(-2))
c       sm (mx,my,nz+1) : Stability function for momentum, at Level 2
c       sh (mx,my,nz+1) : Stability function for heat, at Level 2
c
c       These are defined on the walls of the grid boxes.
c
      subroutine  mym_level2 (
     i            mx, my, nx, ny, nz,
     i            tref, dz, h, pi0,
     i            u, v, tl, qw,
     i            ql, vt, vq,
     o            dtl, dqw, dtv, gm, gh, sm, sh )
c
      dimension  dz(*), h(mx,*), pi0(mx,my,*)
      dimension  u  (mx,my,*), v  (mx,my,*), tl (mx,my,*), qw (mx,my,*)
      dimension  ql (mx,my,*), vt (mx,my,*), vq (mx,my,*)
      dimension  dtl(mx,my,*), dqw(mx,my,*), dtv(mx,my,*)
      dimension  gm (mx,my,*), gh (mx,my,*), sm (mx,my,*), sh (mx,my,*)
c
      common  /constants/ g1, g2, a1, a2, b1, b2, c1, c2, c3, c4, c5
c
      ev  = 2.5e6
      tv0 = 0.61*tref
      tv1 = 1.61*tref
      gtr = 9.81/tref
c
      rfc = g1/( g1+g2 )
      f1  = b1*( g1-c1 ) +3.0*a2*( 1.0    -c2 )*( 1.0-c5 )
     :                   +2.0*a1*( 3.0-2.0*c2 )
      f2  = b1*( g1+g2 ) -3.0*a1*( 1.0    -c2 )
      rf1 = b1*( g1-c1 )/f1
      rf2 = b1*  g1     /f2
      smc = a1 /a2*  f1/f2
      shc = 3.0*a2*( g1+g2 )
c
      ri1 = 0.5/smc
      ri2 = rf1*smc
      ri3 = 4.0*rf2*smc -2.0*ri2
      ri4 = ri2**2
c
      do k = 2,nz+1
        dzk = 0.5  *( dz(k)+dz(k-1) )
        afk = dz(k)/( dz(k)+dz(k-1) )
        abk = 1.0 -afk
      do j = 1,ny
      do i = 1,mx
        duz = ( u(i,j,k)-u(i,j,k-1) )**2 +( v(i,j,k)-v(i,j,k-1) )**2
        duz =   duz                    /( dzk*h(i,j) )**2
        dtz = ( tl(i,j,k)-tl(i,j,k-1) )/( dzk*h(i,j) )
        dqz = ( qw(i,j,k)-qw(i,j,k-1) )/( dzk*h(i,j) )
c
        vtt =  1.0 +vt(i,j,k)*abk +vt(i,j,k-1)*afk
        vqq =  tv0 +vq(i,j,k)*abk +vq(i,j,k-1)*afk
        dtq =  vtt*dtz +vqq*dqz
c
        dtl(i,j,k) =  dtz
        dqw(i,j,k) =  dqz
        dtv(i,j,k) =  dtq
c?      dtv(i,j,k) =  dtz +tv0*dqz
c?   :              +( ev/pi0(i,j,k)-tv1 )
c?   :              *( ql(i,j,k)-ql(i,j,k-1) )/( dzk*h(i,j) )
c
        gm (i,j,k) =  duz
        gh (i,j,k) = -dtq*gtr
c
c   **  Gradient Richardson number  **
        ri = -gh(i,j,k)/max( duz, 1.0e-10 )
c   **  Flux Richardson number  **
        rf = min( ri1*( ri+ri2-sqrt(ri**2-ri3*ri+ri4) ), rfc )
c
        sh (i,j,k) = shc*( rfc-rf )/( 1.0-rf )
        sm (i,j,k) = smc*( rf1-rf )/( rf2-rf ) * sh(i,j,k)
      end do
      end do
      end do
c
      return
      end
c
c     SUBROUTINE  mym_length:
c
c     Input variables:    see subroutine mym_initialize
c
c     Output variables:   see subroutine mym_initialize
c
c     Work arrays:
c       elt(mx,ny)      : Length scale depending on the PBL depth    (m)
c       vsc(mx,ny)      : Velocity scale q_c                       (m/s)
c                         at first, used for computing elt
c
      subroutine  mym_length (
     i            mx, my, nx, ny, nz,
     i            tref, dz, zw, h,
     i            rmo, flt, flq,
     i            vt, vq,
     i            qke,
     i            dtv,
     o            el,
     o            qkw,
     w            elt, vsc )
c
      dimension  dz(*), zw(*), h(mx,*)
      dimension  rmo(mx,*), flt(mx,*), flq(mx,*)
      dimension  vt (mx,*), vq (mx,*)
      dimension  qke(mx,my,*)
      dimension  dtv(mx,my,*)
      dimension  el (mx,my,*)
      dimension  qkw(mx,my,*)
      dimension  elt(mx,*), vsc(mx,*)
c
      data     qmin/ 0.0 /
      data     zmax, cns , alp1, alp2, alp3, alp4
     :       / 1.0 , 2.7 , 0.23, 1.0 , 5.0 , 100./
c
      vk  = 0.4
      tv0 = 0.61*tref
      gtr = 9.81/tref
c
      do k = 2,nz+1
        afk = dz(k)/( dz(k)+dz(k-1) )
        abk = 1.0 -afk
      do j = 1,ny
      do i = 1,mx
        qkw(i,j,k) = sqrt(max(qke(i,j,k)*abk+qke(i,j,k-1)*afk,1.0e-10))
      end do
      end do
      end do
c
      do j = 1,ny
      do i = 1,mx
        elt(i,j) = 1.0e-5
        vsc(i,j) = 1.0e-5
      end do
      end do
c
c   **  Strictly, zwk*h(i,j) -> ( zwk*h(i,j)+z0 )  **
      do k = 2,nz
        zwk = zw(k)
        dzk = 0.5*( dz(k)+dz(k-1) )
      do j = 1,ny
      do i = 1,mx
        qdz = max( qkw(i,j,k)-qmin, 0.0 )*dzk*h(i,j)
        elt(i,j) = elt(i,j) +qdz*zwk*h(i,j)
        vsc(i,j) = vsc(i,j) +qdz
      end do
      end do
      end do
c
      do j = 1,ny
      do i = 1,mx
        vflx = ( vt(i,j)+1.0 )*flt(i,j) +( vq(i,j)+tv0 )*flq(i,j)
        elt(i,j) =  alp1*elt(i,j)/vsc(i,j)
        vsc(i,j) = ( gtr*elt(i,j)*max( vflx, 0.0 ) )**(1.0/3.0)
c
c   **  Strictly, el(i,j,1) is not zero.  **
        el(i,j,1) = 0.0
      end do
      end do
c
      do k = 2,nz+1
        zwk = zw(k)
      do j = 1,ny
      do i = 1,mx
c   **  Length scale limited by the buoyancy effect  **
        if ( dtv(i,j,k) .gt. 0.0 ) then
          bv  = sqrt( gtr*dtv(i,j,k) )
          elb = alp2*qkw(i,j,k) / bv
     :       *( 1.0 + alp3/alp2*sqrt( vsc(i,j)/( bv*elt(i,j) ) ) )
          elf = alp2*qkw(i,j,k) / bv
        else
          elb = 1.0e10
          elf = elb
        end if
c
c   **  Length scale in the surface layer  **
        if ( rmo(i,j) .gt. 0.0 ) then
          els =  vk*zwk*h(i,j)
     :        /( 1.0 + cns*min( zwk*h(i,j)*rmo(i,j), zmax ) )
        else
          els =  vk*zwk*h(i,j)
     :        *( 1.0 - alp4*    zwk*h(i,j)*rmo(i,j)         )**0.2
        end if
c
c  Added: Jan/13/2009, min( ..., elf ) for the free atmosphere
c?      el(i,j,k) =      elb/( elb/elt(i,j)+elb/els+1.0 )
        el(i,j,k) = min( elb/( elb/elt(i,j)+elb/els+1.0 ), elf )
      end do
      end do
      end do
c
      return
      end
c
c     SUBROUTINE  mym_turbulence:
c
c     Input variables:    see subroutine mym_initialize
c       levflag         : <>3;  Level 2.5
c                         = 3;  Level 3
c
c     # ql, vt, vq, qke, tsq, qsq and cov are changed to input variables.
c
c     Output variables:   see subroutine mym_initialize
c       dfm(mx,my,nz+1) : Diffusivity coefficient for momentum,
c                         divided by dz (not dz*h(i,j))            (m/s)
c       dfh(mx,my,nz+1) : Diffusivity coefficient for heat,
c                         divided by dz (not dz*h(i,j))            (m/s)
c       dfq(mx,my,nz+1) : Diffusivity coefficient for q^2,
c                         divided by dz (not dz*h(i,j))            (m/s)
c       tcd(mx,my,nz)   : Countergradient diffusion term for Theta_l
c                                                                  (K/s)
c       qcd(mx,my,nz)   : Countergradient diffusion term for Q_w
c                                                              (kg/kg s)
c       pd?(mx,my,nz+1) : Half of the production terms
c
c       Only tcd and qcd are defined at the center of the grid boxes
c
c     # DO NOT forget that tcd and qcd are added on the right-hand side
c       of the equations for Theta_l and Q_w, respectively.
c
c     Work arrays:        see subroutine mym_initialize and level2
c
c     # dtl, dqw, dtv, gm and gh are allowed to share storage units with
c       dfm, dfh, dfq, tcd and qcd, respectively, for saving memory.
c
      subroutine  mym_turbulence (
     i            levflag, mx, my, nx, ny, nz,
     i            tref, dz, zw, h, pi0,
     i            rmo, flt, flq,
     i            u, v, tl, qw,
     i            ql, vt, vq,
     i            qke, tsq, qsq, cov,
     o            el,
     o            bsh,
     o            dfm, dfh, dfq, tcd, qcd,
     o            pdk, pdt, pdq, pdc,
     w            qkw,
     w            dtl, dqw, dtv, gm, gh, sm, sh )
c
      double precision  q2sq, t2sq, r2sq, c2sq, elsq, gmel, ghel
      double precision  q3sq, t3sq, r3sq, c3sq, dlsq, qdiv
      double precision  e1, e2, e3, e4, enum, eden, wden
c
      dimension  dz(*), zw(*), h(mx,*), pi0(mx,my,*)
      dimension  rmo(mx,*), flt(mx,*), flq(mx,*)
      dimension  u  (mx,my,*), v  (mx,my,*), tl (mx,my,*), qw (mx,my,*)
      dimension  ql (mx,my,*), vt (mx,my,*), vq (mx,my,*)
      dimension  qke(mx,my,*), tsq(mx,my,*), qsq(mx,my,*), cov(mx,my,*)
      dimension  el (mx,my,*)
      dimension  dfm(mx,my,*), dfh(mx,my,*), dfq(mx,my,*)
      dimension  tcd(mx,my,*), qcd(mx,my,*)
      dimension  pdk(mx,my,*), pdt(mx,my,*), pdq(mx,my,*), pdc(mx,my,*)
      dimension  qkw(mx,my,*)
      dimension  dtl(mx,my,*), dqw(mx,my,*), dtv(mx,my,*)
      dimension  gm (mx,my,*), gh (mx,my,*), sm (mx,my,*), sh (mx,my,*)
c
      common  /constants/ g1, g2, a1, a2, b1, b2, c1, c2, c3, c4, c5
c
      tv0 = 0.61*tref
      gtr = 9.81/tref
c
      cc2 =  1.0-c2
      cc3 =  1.0-c3
      e1c =  3.0*a2*b2*cc3
      e2c =  9.0*a1*a2*cc2
      e3c =  9.0*a2*a2*cc2*( 1.0-c5 )
      e4c = 12.0*a1*a2*cc2
      e5c =  6.0*a1*a1
c
      call mym_level2 (
     i            mx, my, nx, ny, nz,
     i            tref, dz, h, pi0,
     i            u, v, tl, qw,
     i            ql, vt, vq,
     o            dtl, dqw, dtv, gm, gh, sm, sh )
c
      call mym_length (
     i            mx, my, nx, ny, nz,
     i            tref, dz, zw, h,
     i            rmo, flt, flq,
     i            vt(1,1,1), vq(1,1,1),
     i            qke,
     i            dtv,
     o            el,
     o            qkw,
     w            pdk(1,1,nz), pdt(1,1,nz) )
c
      do k = 2,nz+1
        dzk = 0.5  *( dz(k)+dz(k-1) )
        afk = dz(k)/( dz(k)+dz(k-1) )
        abk = 1.0 -afk
      do j = 1,ny
      do i = 1,mx
        elsq = el (i,j,k)**2
        q2sq = b1*elsq*( sm(i,j,k)*gm(i,j,k)+sh(i,j,k)*gh(i,j,k) )
        q3sq = qkw(i,j,k)**2
c
c  Modified: Dec/22/2005, from here, (dlsq -> elsq)
        gmel = gm (i,j,k)*elsq
        ghel = gh (i,j,k)*elsq
c  Modified: Dec/22/2005, up to here
c
c     **  Since qkw is set to more than 0.0, q3sq > 0.0.  **
        if ( q3sq .lt. q2sq ) then
          qdiv = sqrt( q3sq/q2sq )
          sm(i,j,k) = sm(i,j,k) * qdiv
          sh(i,j,k) = sh(i,j,k) * qdiv
c
          e1   = q3sq - e1c*ghel * qdiv**2
          e2   = q3sq - e2c*ghel * qdiv**2
          e3   = e1   + e3c*ghel * qdiv**2
          e4   = e1   - e4c*ghel * qdiv**2
          eden = e2*e4 + e3*e5c*gmel * qdiv**2
          eden = max( eden, 1.0d-20 )
        else
          e1   = q3sq - e1c*ghel
          e2   = q3sq - e2c*ghel
          e3   = e1   + e3c*ghel
          e4   = e1   - e4c*ghel
          eden = e2*e4 + e3*e5c*gmel
          eden = max( eden, 1.0d-20 )
c
          qdiv = 1.0
          sm(i,j,k) = q3sq*a1*( e3-3.0*c1*e4       )/eden
          sh(i,j,k) = q3sq*a2*( e2+3.0*c1*e5c*gmel )/eden
        end if
c
c   **  Level 3 : start  **
        if ( levflag .eq. 3 ) then
          t2sq = qdiv*b2*elsq*sh(i,j,k)*dtl(i,j,k)**2
          r2sq = qdiv*b2*elsq*sh(i,j,k)*dqw(i,j,k)**2
          c2sq = qdiv*b2*elsq*sh(i,j,k)*dtl(i,j,k)*dqw(i,j,k)
          t3sq = max( tsq(i,j,k)*abk+tsq(i,j,k-1)*afk, 0.0 )
          r3sq = max( qsq(i,j,k)*abk+qsq(i,j,k-1)*afk, 0.0 )
          c3sq =      cov(i,j,k)*abk+cov(i,j,k-1)*afk
c
c  Modified: Dec/22/2005, from here
          c3sq = sign( min( abs(c3sq), sqrt(t3sq*r3sq) ), c3sq )
c
          vtt  = 1.0 +vt(i,j,k)*abk +vt(i,j,k-1)*afk
          vqq  = tv0 +vq(i,j,k)*abk +vq(i,j,k-1)*afk
          t2sq = vtt*t2sq +vqq*c2sq
          r2sq = vtt*c2sq +vqq*r2sq
          c2sq = max( vtt*t2sq+vqq*r2sq, 0.0d0 )
          t3sq = vtt*t3sq +vqq*c3sq
          r3sq = vtt*c3sq +vqq*r3sq
          c3sq = max( vtt*t3sq+vqq*r3sq, 0.0d0 )
c
          cw25 = e1*( e2 + 3.0*c1*e5c*gmel*qdiv**2 )/( 3.0*eden )
c
c     **  Limitation on q, instead of L/q  **
                                           dlsq =  elsq
          if ( q3sq/dlsq .lt. -gh(i,j,k) ) q3sq = -dlsq*gh(i,j,k)
c
c     **  Limitation on c3sq (0.12 =< cw =< 0.76) **
          e2   = q3sq - e2c*ghel * qdiv**2
          e3   = q3sq + e3c*ghel * qdiv**2
          e4   = q3sq - e4c*ghel * qdiv**2
          eden = e2*e4  + e3 *e5c*gmel * qdiv**2
c
          wden = cc3*gtr**2 * dlsq**2/elsq * qdiv**2
     :        *( e2*e4c - e3c*e5c*gmel * qdiv**2 )
c
          if ( wden .ne. 0.0 ) then
            clow = q3sq*( 0.12-cw25 )*eden/wden
            cupp = q3sq*( 0.76-cw25 )*eden/wden
c
            if ( wden .gt. 0.0 ) then
              c3sq  = min( max( c3sq, c2sq+clow ), c2sq+cupp )
            else
              c3sq  = max( min( c3sq, c2sq+clow ), c2sq+cupp )
            end if
          end if
c
          e1   = e2 + e5c*gmel * qdiv**2
          eden = max( eden, 1.0d-20 )
c  Modified: Dec/22/2005, up to here
c
          e6c  = 3.0*a2*cc3*gtr * dlsq/elsq
c
c     **  for Gamma_theta  **
          enum = qdiv*e6c*( t3sq-t2sq )
          gamt =-e1  *enum    /eden
c
c     **  for Gamma_q  **
          enum = qdiv*e6c*( r3sq-r2sq )
          gamq =-e1  *enum    /eden
c
c     **  for Sm' and Sh'd(Theta_V)/dz  **
          enum = qdiv*e6c*( c3sq-c2sq )
          smd  = dlsq*enum*gtr/eden * qdiv**2 * (e3c+e4c)*a1/a2
          gamv = e1  *enum*gtr/eden
c
          sm(i,j,k) = sm(i,j,k) +smd
c
c     **  For elh (see below), qdiv at Level 3 is reset to 1.0.  **
          qdiv = 1.0
c   **  Level 3 : end  **
c
        else
c     **  At Level 2.5, qdiv is not reset.  **
          gamt = 0.0
          gamq = 0.0
          gamv = 0.0
        end if
c
        elq = el(i,j,k)*qkw(i,j,k)
        elh = elq*qdiv
c
        pdk(i,j,k) = elq*( sm(i,j,k)*gm (i,j,k)
     :                    +sh(i,j,k)*gh (i,j,k)+gamv )
        pdt(i,j,k) = elh*( sh(i,j,k)*dtl(i,j,k)+gamt )*dtl(i,j,k)
        pdq(i,j,k) = elh*( sh(i,j,k)*dqw(i,j,k)+gamq )*dqw(i,j,k)
        pdc(i,j,k) = elh*( sh(i,j,k)*dtl(i,j,k)+gamt )*dqw(i,j,k)*0.5
     :              +elh*( sh(i,j,k)*dqw(i,j,k)+gamq )*dtl(i,j,k)*0.5
c
        tcd(i,j,k) = elq*gamt
        qcd(i,j,k) = elq*gamq
c
        dfm(i,j,k) = elq*sm (i,j,k) / dzk
        dfh(i,j,k) = elq*sh (i,j,k) / dzk
c  Modified: Dec/22/2005, from here
c   **  In sub.mym_predict, dfq for the TKE and scalar variance **
c   **  are set to 3.0*dfm and 1.0*dfm, respectively.           **
        dfq(i,j,k) =     dfm(i,j,k)
c  Modified: Dec/22/2005, up to here
      end do
      end do
      end do
c
      do j = 1,ny
      do i = 1,mx
        dfm(i,j,1) = 0.0
        dfh(i,j,1) = 0.0
        dfq(i,j,1) = 0.0
        tcd(i,j,1) = 0.0
        qcd(i,j,1) = 0.0
c
        tcd(i,j,nz+1) = 0.0
        qcd(i,j,nz+1) = 0.0
      end do
      end do
c
      do k = 1,nz
        dzk = dz(k)
      do j = 1,ny
      do i = 1,mx
        tcd(i,j,k) = ( tcd(i,j,k+1)-tcd(i,j,k) )/( dzk*h(i,j) )
        qcd(i,j,k) = ( qcd(i,j,k+1)-qcd(i,j,k) )/( dzk*h(i,j) )
      end do
      end do
      end do
c
      return
      end
c
c     SUBROUTINE  mym_predict:
c
c     Input variables:    see subroutine mym_initialize and turbulence
c       delt            : Time step                                (sec)
c       qke(mx,my,nz+1) : qke at (n)th time level
c       tsq, ...cov     : ditto
c
c     Output variables:
c       qke(mx,my,nz+1) : qke at (n+1)th time level
c       tsq, ...cov     : ditto
c
c     Work arrays:
c       qkw(mx,my,nz)   : q at the center of the grid boxes        (m/s)
c       bp (mx,my,nz)   : = 1/2*F,     see below
c       rp (mx,my,nz)   : = P-1/2*F*Q, see below
c
c     # The equation for a turbulent quantity Q can be expressed as
c          dQ/dt + Ah + Av = Dh + Dv + P - F*Q,                      (1)
c       where A is the advection, D the diffusion, P the production,
c       F*Q the dissipation and h and v denote horizontal and vertical,
c       respectively. If Q is q^2, F is 2q/B_1L.
c       Using the Crank-Nicholson scheme for Av, Dv and F*Q, a finite
c       difference equation is written as
c          Q{n+1} - Q{n} = dt  *( Dh{n}   - Ah{n}   + P{n} )
c                        + dt/2*( Dv{n}   - Av{n}   - F*Q{n}   )
c                        + dt/2*( Dv{n+1} - Av{n+1} - F*Q{n+1} ),    (2)
c       where n denotes the time level.
c       When the advection and diffusion terms are discretized as
c          dt/2*( Dv - Av ) = a(k)Q(k+1) - b(k)Q(k) + c(k)Q(k-1),    (3)
c       Eq.(2) can be rewritten as
c          - a(k)Q(k+1) + [ 1 + b(k) + dt/2*F ]Q(k) - c(k)Q(k-1)
c                 = Q{n} + dt  *( Dh{n}   - Ah{n}   + P{n} )
c                        + dt/2*( Dv{n}   - Av{n}   - F*Q{n}   ),    (4)
c       where Q on the left-hand side is at (n+1)th time level.
c
c       Time integration, including computation of a(k), b(k) and c(k),
c       is performed in subprogram time_integration. F and P are stored
c       in bp and rp, respectively. For the Crank-Nicholson scheme,
c       dt/2*bp and dt*rp-dt/2*bp*Q need to be computed.
c
c       Modify this subroutine according to your numerical integration
c       scheme (program).
c
      subroutine  mym_predict (
     i            levflag, mx, my, nx, ny, nz, delt,
     i            dz, h,
     i            flt, flq, ust, pmz, phh,
     m            qke, tsq, qsq, cov,
     i            el,
     i            dfq,
     m            pdk, pdt, pdq, pdc,
     w            qkw,
     w            bp, rp )
c
      dimension  dz(*), h(mx,*)
      dimension  flt(mx,*), flq(mx,*), ust(mx,*)
      dimension  pmz(mx,*), phh(mx,*)
      dimension  qke(mx,my,*), tsq(mx,my,*), qsq(mx,my,*), cov(mx,my,*)
      dimension  el (mx,my,*)
      dimension  dfq(mx,my,*)
      dimension  pdk(mx,my,*), pdt(mx,my,*), pdq(mx,my,*), pdc(mx,my,*)
      dimension  qkw(mx,my,*)
      dimension  bp (mx,my,*), rp (mx,my,*)
c
      common  /constants/ g1, g2, a1, a2, b1, b2, c1, c2, c3, c4, c5
c
c   **  Strictly, vkz*h(i,j) -> vk*( 0.5*dz(1)*h(i,j)+z0 )  **
      vkz = 0.4*0.5*dz(1)
c
      do k = 1,nz
      do j = 1,ny
      do i = 1,mx
        qkw(i,j,k) = sqrt( max( qke(i,j,k), 0.0 ) )
      end do
      end do
      end do
c
      do j = 1,ny
      do i = 1,mx
        pdk1 = 2.0*ust(i,j)**3*pmz(i,j)/( vkz*h(i,j) )
        phm  = 2.0/ust(i,j)   *phh(i,j)/( vkz*h(i,j) )
        pdt1 = phm*flt(i,j)**2
        pdq1 = phm*flq(i,j)**2
        pdc1 = phm*flt(i,j)*flq(i,j)
c
c   **  pdk(i,j,1)+pdk(i,j,2) corresponds to pdk1.  **
        pdk(i,j,1) = pdk1 -pdk(i,j,2)
c
c  Modified: Jan/13/2009, (should be applied only to level 2.5?)
c   **  For some reason, pdt(i,j,1)=pdt(i,j,2) seems to be better  **
c   **  to obtain a better value of sgm in sub mym_condensation.   **
c?      pdt(i,j,1) = pdt1 -pdt(i,j,2)
c?      pdq(i,j,1) = pdq1 -pdq(i,j,2)
c?      pdc(i,j,1) = pdc1 -pdc(i,j,2)
        pdt(i,j,1) = pdt(i,j,2)
        pdq(i,j,1) = pdq(i,j,2)
        pdc(i,j,1) = pdc(i,j,2)
      end do
      end do
c
c   **  Prediction of twice the turbulent kinetic energy  **
      do k = 1,nz
      do j = 1,ny
      do i = 1,mx
        b1l = b1*0.5*( el(i,j,k+1)+el(i,j,k) )
        bp(i,j,k) = 2.0*qkw(i,j,k) / b1l
        rp(i,j,k) = pdk(i,j,k+1) +pdk(i,j,k)
      end do
      end do
      end do
c
c   **  dfq for the TKE is 3.0*dfm.  **
      call time_integration (
     i            0.0, mx, my, nx, ny, nz, delt,
     i            dz, h,
     m            qke,
     i            dfq, 3.0,
     i            bp, rp )
c
      if ( levflag .eq. 3 ) then
c
c   **  Prediction of the temperature variance  **
        do k = 1,nz
        do j = 1,ny
        do i = 1,mx
          b2l = b2*0.5*( el(i,j,k+1)+el(i,j,k) )
          bp(i,j,k) = 2.0*qkw(i,j,k) / b2l
          rp(i,j,k) = pdt(i,j,k+1) +pdt(i,j,k)
        end do
        end do
        end do
c
c   **  dfq for the scalar variance is 1.0*dfm.  **
        call time_integration (
     i            0.0, mx, my, nx, ny, nz, delt,
     i            dz, h,
     m            tsq,
     i            dfq, 1.0,
     i            bp, rp )
c
c   **  Prediction of the moisture variance  **
        do k = 1,nz
        do j = 1,ny
        do i = 1,mx
          b2l = b2*0.5*( el(i,j,k+1)+el(i,j,k) )
          bp(i,j,k) = 2.0*qkw(i,j,k) / b2l
          rp(i,j,k) = pdq(i,j,k+1) +pdq(i,j,k)
        end do
        end do
        end do
c
        call time_integration (
     i            0.0, mx, my, nx, ny, nz, delt,
     i            dz, h,
     m            qsq,
     i            dfq, 1.0,
     i            bp, rp )
c
c   **  Prediction of the temperature-moisture covariance  **
        do k = 1,nz
        do j = 1,ny
        do i = 1,mx
          b2l = b2*0.5*( el(i,j,k+1)+el(i,j,k) )
          bp(i,j,k) = 2.0*qkw(i,j,k) / b2l
          rp(i,j,k) = pdc(i,j,k+1) +pdc(i,j,k)
        end do
        end do
        end do
c
        call time_integration (
     i            0.0, mx, my, nx, ny, nz, delt,
     i            dz, h,
     m            cov,
     i            dfq, 1.0,
     i            bp, rp )
c
      else
        do k = 1,nz
        do j = 1,ny
        do i = 1,mx
          if ( qkw(i,j,k) .le. 0.0 ) then
            b2l = 0.0
          else
            b2l = b2*0.25*( el(i,j,k+1)+el(i,j,k) )/qkw(i,j,k)
          end if
c
          tsq(i,j,k) = b2l*( pdt(i,j,k+1)+pdt(i,j,k) )
          qsq(i,j,k) = b2l*( pdq(i,j,k+1)+pdq(i,j,k) )
          cov(i,j,k) = b2l*( pdc(i,j,k+1)+pdc(i,j,k) )
        end do
        end do
        end do
      end if
c
      return
      end
c
c     SUBROUTINE  mym_condensation:
c
c     Input variables:    see subroutine mym_initialize and turbulence
c       pi (mx,my,nz+1) : Perturbation of the Exner function    (J/kg K)
c                         defined on the walls of the grid boxes
c                         This is usually computed by integrating
c                         d(pi)/dz = h*g*tv/tref**2
c                         from the upper boundary, where tv is the
c                         virtual potential temperature minus tref.
c
c     Output variables:   see subroutine mym_initialize
c       cld(mx,my,nz)   : Cloud fraction
c
c     Work arrays:
c       qmq(mx,my,nz)   : Q_w-Q_{sl}, where Q_{sl} is the saturation
c                         specific humidity at T=Tl
c       alp(mx,my,nz)   : Functions in the condensation process
c       bet(mx,my,nz)   : ditto
c       sgm(mx,my,nz)   : Combined standard deviation sigma_s
c                         multiplied by 2/alp
c
c     # qmq, alp, bet and sgm are allowed to share storage units with
c       any four of other work arrays for saving memory.
c
c     # Results are sensitive particularly to values of cp and rd.
c       Set these values to those adopted by you.
c
      subroutine  mym_condensation (
     i            levflag, mx, my, nx, ny, nz,
     i            tref, dz, h, pi0,
     i            tl, qw,
     i            pi,
     i            tsq, qsq, cov,
     i            bsh,
     o            ql, vt, vq,
     o            cld,
     w            qmq, alp, bet, sgm )
c
      double precision  t3sq, r3sq, c3sq
c
      dimension  dz(*), h(mx,*), pi0(mx,my,*)
      dimension  tl (mx,my,*), qw (mx,my,*)
      dimension  pi (mx,my,*)
      dimension  tsq(mx,my,*), qsq(mx,my,*), cov(mx,my,*)
      dimension  ql (mx,my,*), vt (mx,my,*), vq (mx,my,*)
      dimension  cld(mx,my,*)
      dimension  qmq(mx,my,*), alp(mx,my,*), bet(mx,my,*), sgm(mx,my,*)
c
      data     rr2, rrp/ 0.7071068, 0.3989423/
c
      ev  = 2.5e6
      cp  = 1004.0
      rd  =  287.0
      rk  = cp/rd
c
      do k = 1,nz
      do j = 1,ny
      do i = 1,mx
        p2a = 0.5*( pi(i,j,k+1)+pi0(i,j,k+1)
     :             +pi(i,j, k )+pi0(i,j, k ) ) / cp
        ct  = ( tl(i,j,k)+tref )*p2a -273.15
cx      if ( ct .gt. 0.0 ) then
          a  =  17.27
          b  = 237.3
cx      else
cx        a  =  21.87
cx        b  = 265.5
cx      end if
c
c   **  3.8 = 0.622*6.11 (hPa)  **
        qsl  = 3.8*exp( a*ct/( b+ct ) ) / ( 1000.0*p2a**rk )
        dqsl = qsl*0.622*ev/( rd*( ct+273.15 )**2 )
c
        qmq(i,j,k) = qw(i,j,k) -qsl
        alp(i,j,k) = 1.0/( 1.0+dqsl*ev/cp )
        bet(i,j,k) = dqsl*p2a
c
        t3sq = max( tsq(i,j,k), 0.0 )
        r3sq = max( qsq(i,j,k), 0.0 )
        c3sq =      cov(i,j,k)
        c3sq = sign( min( abs(c3sq), sqrt(t3sq*r3sq) ), c3sq )
c
        r3sq = r3sq +bet(i,j,k)**2*t3sq -2.0*bet(i,j,k)*c3sq
        sgm(i,j,k) = sqrt( max( r3sq, 1.0d-10 ) )
      end do
      end do
      end do
c
      do k = 1,nz
      do j = 1,ny
      do i = 1,mx
        q1   = qmq(i,j,k) / sgm(i,j,k)
        cld0 = 0.5*( 1.0+erf( q1*rr2 ) )
        eq1  = rrp*exp( -0.5*q1*q1 )
        qll  = max( cld0*q1 + eq1, 0.0 )
c
        cld(i,j,k) = cld0
        ql (i,j,k) = alp(i,j,k)*sgm(i,j,k)*qll
c
        q2p  = 2.0*ev/( pi(i,j,k+1)+pi0(i,j,k+1)
     :                 +pi(i,j, k )+pi0(i,j, k ) )
        pt   = tl(i,j,k) +tref +q2p*ql(i,j,k)
        qt   = 1.0 +0.61*qw(i,j,k) -1.61*ql(i,j,k)
        rac  = alp(i,j,k)*( cld0-qll*eq1 )*( q2p*qt-1.61*pt )
c
        vt (i,j,k) =        qt-1.0    -rac*bet(i,j,k)
        vq (i,j,k) = 0.61*( pt-tref ) +rac
      end do
      end do
      end do
c
      do j = 1,ny
      do i = 1,mx
        ql(i,j,nz+1) = ql(i,j,nz)
        vt(i,j,nz+1) = vt(i,j,nz)
        vq(i,j,nz+1) = vq(i,j,nz)
      end do
      end do
c
      return
      end
c
c     SUBROUTINE  time_integration:
c
c     Input variables:
c       vbound          : Upper boundary condition
c                         >=0.0; equal to specified gradient
c                         < 0.0; set to prescribed value (=qqq(nz+1))
c
c     Sample program of the full-implicit time-integration scheme
c     based on the tridiagonal matrix for a one-dimensional equation
c
c     NOTE: dfq(i,j,1) must be zero.
c
      subroutine  time_integration (
     i            vbound, mx, my, nx, ny, nz, delt,
     i            dz, h,
     m            qqq,
     i            dfq, fac,
     i            bp, rp )
c
      dimension  dz(*), h(mx,*)
      dimension  qqq(mx,my,*)
      dimension  dfq(mx,my,*)
      dimension  bp (mx,my,*), rp (mx,my,*)
c
      dimension  e(500), f(500)
c
      do j = 1,ny
      do i = 1,nx
c
        do k = 1,nz
          km1 = max( k-1, 1 )
          dtz = fac * delt/( dz(k)*h(i,j)**2 )
c
          a   = dtz*  dfq(i,j,k+1)
          b   = dtz*( dfq(i,j,k+1)+dfq(i,j,k) ) +delt*bp(i,j,k) +1.0
          c   = dtz*  dfq(i,j, k )
          d   = qqq(i,j,k) +delt*rp(i,j,k)
c
          den  =   b-c*e(km1)
          e(k) =   a           /den
          f(k) = ( d+c*f(km1) )/den
        end do
c
        if ( vbound .ge. 0.0 ) then
          d  = vbound*0.5*( dz(nz+1)+dz(nz) ) * h(i,j)
          qqq(i,j,nz+1) = ( d+f(nz) )/( 1.0-e(nz) )
        end if
c
        do k = nz,1,-1
          qqq(i,j,k) = e(k)*qqq(i,j,k+1) +f(k)
        end do
c
      end do
      end do
c
      return
      end