MODULE Shall_MasFlux
  USE Constants, ONLY :  i8
  
  USE Mod_GET_PRS, ONLY: GET_PHI

  IMPLICIT NONE
  PRIVATE
  INTEGER, PARAMETER :: r8 = SELECTED_REAL_KIND(13,60) ! the '60' maps to 64-bit real
  !  --- ...  Geophysics/Astronomy constants

  REAL(kind=r8),PARAMETER:: con_g      =9.80665e+0_r8     ! gravity           (m/s2)
  REAL(kind=r8),PARAMETER:: con_pi     =3.1415926535897931 ! pi
  REAL(kind=r8),PARAMETER:: con_rerth  =6.3712e+6      ! radius of earth   (m)
  !  --- ...  Thermodynamics constants


  REAL(kind=r8),PARAMETER:: con_cp     =1.0046e+3_r8      ! spec heat air @p    (J/kg/K)
  REAL(kind=r8),PARAMETER:: con_rv     =4.6150e+2_r8      ! gas constant H2O    (J/kg/K)
  REAL(kind=r8),PARAMETER:: con_hvap   =2.5000e+6_r8      ! lat heat H2O cond   (J/kg)
  REAL(kind=r8),PARAMETER:: con_hfus   =3.3358e+5_r8      ! lat heat H2O fusion (J/kg)
  REAL(kind=r8),PARAMETER:: con_rd     =2.8705e+2_r8      ! gas constant air    (J/kg/K)
  REAL(kind=r8),PARAMETER:: con_cvap   =1.8460e+3_r8      ! spec heat H2O gas   (J/kg/K)
  REAL(kind=r8),PARAMETER:: con_cliq   =4.1855e+3_r8      ! spec heat H2O liq   (J/kg/K)
  REAL(kind=r8),PARAMETER:: con_csol   =2.1060e+3_r8      ! spec heat H2O ice   (J/kg/K)
  REAL(kind=r8),PARAMETER:: con_ttp    =2.7316e+2_r8      ! temp at H2O 3pt     (K)
  REAL(kind=r8),PARAMETER:: con_psat   =6.1078e+2_r8      ! pres at H2O 3pt     (Pa)  
  real(kind=r8),PARAMETER:: con_t0c    =2.7315e+2_r8      ! temp at 0C          (K)

  !  Secondary constants

  REAL(kind=r8),PARAMETER:: con_rocp   =con_rd/con_cp
  REAL(kind=r8),PARAMETER:: con_fvirt  =con_rv/con_rd-1.0_r8 !(J/kg/K)/(J/kg/K)
  REAL(kind=r8),PARAMETER:: con_eps    =con_rd/con_rv
  REAL(kind=r8),PARAMETER:: con_epsm1  =con_rd/con_rv-1.0_r8

  !
  !   module funcphys
  INTEGER,PARAMETER:: nxpvs=7501

  REAL(r8) c1xpvs,c2xpvs,tbpvs(nxpvs)

  !   END module funcphys
   REAL (kind=r8)      , ALLOCATABLE ::si(:)
   REAL (kind=r8)      , ALLOCATABLE ::sl(:)
   INTEGER  :: jcap


   PUBLIC :: Init_Shall_MasFlux
   PUBLIC :: Run_Shall_MasFlux
CONTAINS

  SUBROUTINE Init_Shall_MasFlux (kMax,trunc,si_in,sl_in)
    IMPLICIT NONE
    INTEGER      , INTENT(IN   ) ::  kMax
    INTEGER      , INTENT(IN   ) ::  trunc
    REAL(kind=r8), INTENT(IN   ) ::  si_in(kMax+1)
    REAL(kind=r8), INTENT(IN   ) ::  sl_in(kMax)
    
    jcap=trunc
    
    ALLOCATE (si(kMax+1));si=si_in
    ALLOCATE (sl(kMax))  ;sl=sl_in

    CALL gfuncphys()

  END SUBROUTINE Init_Shall_MasFlux

  SUBROUTINE Run_Shall_MasFlux(nCols,kMax,dt,ps,tgrs,qgrs,ugrs,vgrs,omgb,qliq,qice,dudt,dvdt,kbot,ktop,kuo,noshal,mask,hpbl,sens,latheat)
    INTEGER      , INTENT(IN   ) :: nCols
    INTEGER      , INTENT(IN   ) :: kMax
    REAL(KIND=r8), INTENT(IN   ) :: dt

    REAL(KIND=r8), INTENT(IN   ) :: ps   (nCols)        !cb
    REAL(KIND=r8), INTENT(INOUT) :: tgrs (nCols,kMax)!     tgrs     - real, layer mean temperature ( k )             ix,levs !
    REAL(KIND=r8), INTENT(INOUT) :: qgrs (nCols,kMax)!     qgrs     - real, layer mean tracer concentration     ix,levs,ntrac!
    REAL(KIND=r8), INTENT(IN   ) :: ugrs (nCols,kMax)!     ugrs,vgrs- real, u/v component of layer wind              ix,levs !
    REAL(KIND=r8), INTENT(IN   ) :: vgrs (nCols,kMax)
    REAL(KIND=r8), INTENT(IN   ) :: omgb (nCols,kMax) ! (Pa/s) !REAL(KIND=r8)   ,    INTENT(IN   ) :: omgb      (iMax,kMax) ! (Pa/s)
    REAL(KIND=r8), INTENT(INOUT) :: qliq (nCols,kMax)!     qgrs     - real, layer mean tracer concentration     ix,levs,ntrac!
    REAL(KIND=r8), INTENT(INOUT) :: qice (nCols,kMax)!     qgrs     - real, layer mean tracer concentration     ix,levs,ntrac!
    REAL(KIND=r8), INTENT(inout) :: dudt (nCols,kMax)
    REAL(KIND=r8), INTENT(inout) :: dvdt (nCols,kMax)
    INTEGER      , INTENT(INOUT) :: kbot (nCols)
    INTEGER      , INTENT(INOUT) :: ktop (nCols)
    INTEGER      , INTENT(IN   ) :: kuo(nCols)    !   kuo                     ! convection yes(1) or not(0) for shallow convection
    INTEGER      , INTENT(INOUT) :: noshal (nCols) 
    INTEGER(KIND=i8), INTENT(IN   ) :: mask    (1:nCols) 
    REAL(kind=r8), INTENT(IN   ) :: hpbl (nCols)  ! hpbl [m]
    REAL(KIND=r8), INTENT(in   ) :: sens    (1:nCols)
    REAL(KIND=r8), INTENT(in   ) :: latheat (1:nCols)

    REAL(kind=r8) :: heat (nCols)    !     hflx     - real, sensible heat flux  			im   !
    REAL(kind=r8) :: evap (nCols)    !     evap     - real, evaperation from latent heat flux		im   !
    REAL(kind=r8) :: rcs   (nCols)
    REAL(KIND=r8) :: prsi  (nCols,kMax+1) !     prsi     - real, pressure at layer interfaces             ix,levs+1
    REAL(KIND=r8) :: prsl  (nCols,kMax)   !     prsl     - real, mean layer presure                       ix,levs !
    REAL(KIND=r8) :: prsik (nCols,kMax+1)
    REAL(KIND=r8) :: prslk (nCols,kMax)
    REAL(KIND=r8) :: phii  (nCols,kMax+1)
    REAL(KIND=r8) :: phil  (nCols,kMax)
    REAL(KIND=r8) :: del   (nCols,kMax)
    REAL(KIND=r8) :: ql    (nCols,kMax,2)! 1-ice    ! 2-water
    REAL(kind=r8) :: dot   (nCols,kMax)!     vvel     - real, layer mean vertical velocity             ix,levs !
    REAL(kind=r8) :: rho   (nCols,kMax)
    REAL(kind=r8) :: rn    (nCols)
    REAL(kind=r8) :: q1    (nCols,kMax)
    REAL(kind=r8) :: t1    (nCols,kMax)
    REAL(kind=r8) :: u1    (nCols,kMax)
    REAL(kind=r8) :: v1    (nCols,kMax)
    REAL(kind=r8) :: slimsk(nCols)!           slmsk    - real, sea/land/ice mask (=0/1/2)                  im   !
    INTEGER       :: ncloud   !     ncld     - integer, number of cloud species                  1    !
    REAL(kind=r8) :: ud_mf (nCols,kMax)
    REAL(kind=r8) :: dt_mf (nCols,kMax)
    REAL(kind=r8) :: shalldudt (nCols,kMax)
    REAL(kind=r8) :: shalldvdt (nCols,kMax)

    REAL(KIND=r8) :: delt
    INTEGER       :: ntrac,i,k
    ntrac=0
    ncloud=1
    rcs=1.0_r8
    delt=2.0_r8*dt
    ! grell mask
    DO i=1,nCols      
    
       rcs  (i) = 1.0_r8       ! wind conversion
       heat (i) = sens       (i)   !     hflx     - real, sensible heat fluxim   !
       evap (i) = latheat    (i)   !     evap     - real, evaperation from latent heat fluxim   !
       IF(mask(i).GT.0_i8)THEN
          !mask2(i)=0 ! land
          slimsk(i) = 1.0_r8
       ELSE
          slimsk(i) = 0.0_r8
          !mask2(i)=1 ! seaice/water/ocean
       END IF
    END DO

    DO k=1,kMax
       DO i=1,nCols
         shalldudt(i,k)= 0.0_r8
         shalldvdt(i,k)= 0.0_r8
         ql(i,k,1) = qice (i,k)
         ql(i,k,2) = qliq (i,k)
         q1(i,k)   = qgrs (i,k) 
         t1(i,k)   = tgrs (i,k)
         u1(i,k)   = ugrs (i,k)
         v1(i,k)   = vgrs (i,k)
       END DO
    END DO
       CALL sig2press(nCols    ,& !
                      kMax     ,&!
                      ps       ,&! cb
                      sl       ,&!
                      si       ,&!
                      prsi  (1:nCols,1:kMax+1)   ,&!!cb
                      prsl  (1:nCols,1:kMax)     ,&!!cb
                      prsik (1:nCols,1:kMax+1)   ,&!!cb
                      prslk (1:nCols,1:kMax)       )!cb
  
       CALL GET_PHI(&
                    nCols                     , &!INTEGER         , INTENT(IN   ) :: ix
                    kMax                      , &!INTEGER         , INTENT(IN   ) :: levs
                    1                         , &!INTEGER         , INTENT(IN   ) :: ntrac
                    tgrs  (1:nCols,1:kMax)    , &!, &!REAL(kind=r8), INTENT(IN   ) :: T(ix,levs)
                    qgrs  (1:nCols,1:kMax)    , &!, &!REAL(kind=r8), INTENT(IN   ) :: q(ix,levs,ntrac)
                    prsi  (1:nCols,1:kMax+1)  , &!REAL(kind=r8), INTENT(IN   ) :: prsi(ix,levs+1)
                    prsik (1:nCols,1:kMax+1)  , &!REAL(kind=r8), INTENT(IN   ) :: prki(ix,levs+1)
                    prslk (1:nCols,1:kMax)    , &!REAL(kind=r8), INTENT(IN   ) :: prkl(ix,levs)
                    phii  (1:nCols,1:kMax+1)  , &!REAL(kind=r8), INTENT(out   ) :: phii(ix,levs+1)
                    phil  (1:nCols,1:kMax)    , &!REAL(kind=r8), INTENT(out   ) :: phil(ix,levs)
                    del   (1:nCols,1:kMax)      )
    DO k=1,kMax
       DO i=1,nCols
         !; p_at_units = "Pa" 
         !; t_at_units = "K" 
         !   RGAS = 287. ; J/(kg-K) => m2/(s2 K) 
         !; omega_at_units = "Pa/sec" 
         !   rho = p/(RGAS*t) ; density => kg/m3 
         rho(i,k) = (1000.0_r8*prsl(i,k))/(con_rd*tgrs (i,k))
         !   GRAV = 9.8 ; m/s2 

         !   w = -omega/(rho*GRAV) 

         dot(i,k) = -omgb (i,k)/(rho(i,k)*  con_g ) 
       END DO
    END DO
    !    prepare input, erase output
    !
    DO i=1,nCols
       noshal(i)=0
       kbot(i)=1
       ktop(i)=1
    END DO

    CALL shalcnv(                     &
       nCols                        , &! INTEGER      , INTENT(IN   ) :: im
       kMax                         , &! INTEGER      , INTENT(IN   ) :: km
       jcap                         , &! INTEGER      , INTENT(IN   ) :: jcap
       delt                         , &! REAL(kind=r8), INTENT(IN   ) :: delt
       del    (1:nCols,1:kMax)      , &! REAL(kind=r8), INTENT(IN   ) :: del   (im,km)
       prsl   (1:nCols,1:kMax)      , &! REAL(kind=r8), INTENT(IN   ) :: prsl  (im,km)
       ps     (1:nCols)             , &! REAL(kind=r8), INTENT(IN   ) :: ps    (im)
       phil   (1:nCols,1:kMax)      , &! REAL(kind=r8), INTENT(IN   ) :: phil  (im,km)
       ql     (1:nCols,1:kMax,1:2)  , &! REAL(kind=r8), INTENT(INOUT) :: ql    (im,km,2)
       q1     (1:nCols,1:kMax)      , &! REAL(kind=r8), INTENT(INOUT) :: q1    (im,km)
       t1     (1:nCols,1:kMax)      , &! REAL(kind=r8), INTENT(INOUT) :: t1    (im,km)
       u1     (1:nCols,1:kMax)      , &! REAL(kind=r8), INTENT(INOUT) :: u1    (im,km)
       v1     (1:nCols,1:kMax)      , &! REAL(kind=r8), INTENT(INOUT) :: v1    (im,km)
       rcs    (1:nCols)             , &! REAL(kind=r8), INTENT(IN   ) :: rcs   (im)
       rn     (1:nCols)             , &! REAL(kind=r8), INTENT(OUT  ) :: rn    (im)
       kbot   (1:nCols)             , &! INTEGER      , INTENT(INOUT) :: kbot  (im)
       ktop   (1:nCols)             , &! INTEGER      , INTENT(INOUT) :: ktop  (im)
       kuo    (1:nCols)             , &! INTEGER      , INTENT(in) :: ktop  (im)
       noshal   (1:nCols)             , &! INTEGER      , INTENT(INOUT) :: noshal  (im) 
       slimsk (1:nCols)             , &! REAL(kind=r8), INTENT(IN   ) :: slimsk(im)
       dot    (1:nCols,1:kMax)      , &! REAL(kind=r8), INTENT(IN   ) :: dot   (im,km)
       ncloud                       , &! INTEGER      , INTENT(IN   ) :: ncloud 
       hpbl   (1:nCols)             , &! REAL(kind=r8), INTENT(IN   ) :: hpbl  (im)
       heat   (1:nCols)             , &! REAL(kind=r8), INTENT(IN   ) :: heat  (im)
       evap   (1:nCols)             , &! REAL(kind=r8), INTENT(IN   ) :: evap  (im)
       ud_mf  (1:nCols,1:kMax)      , &! REAL(kind=r8), INTENT(OUT  ) :: ud_mf (im,km)
       dt_mf  (1:nCols,1:kMax)        )! REAL(kind=r8), INTENT(OUT  ) :: dt_mf (im,km)

    DO k=1,kMax
       DO i=1,nCols
          qice (i,k) = ql(i,k,1) 
          qliq (i,k) = ql(i,k,2) 
          qgrs (i,k) = q1(i,k)
          tgrs (i,k) = t1(i,k)
          shalldudt (i,k) = (u1(i,k) - ugrs (i,k) )/delt
          shalldvdt (i,k) = (v1(i,k) - vgrs (i,k) )/delt
          dudt (i,k) = dudt (i,k) + shalldudt (i,k)
          dvdt (i,k) = dvdt (i,k) + shalldvdt (i,k)
       END DO
    END DO


  END SUBROUTINE Run_Shall_MasFlux

  SUBROUTINE shalcnv( &
       im        , &! INTEGER      , INTENT(IN   ) :: im
       km        , &! INTEGER      , INTENT(IN   ) :: km
       jcap      , &! INTEGER      , INTENT(IN   ) :: jcap
       delt      , &! REAL(kind=r8), INTENT(IN   ) :: delt
       del       , &! REAL(kind=r8), INTENT(IN   ) :: del   (im,km)
       prsl      , &! REAL(kind=r8), INTENT(IN   ) :: prsl  (im,km)
       ps        , &! REAL(kind=r8), INTENT(IN   ) :: ps    (im)
       phil      , &! REAL(kind=r8), INTENT(IN   ) :: phil  (im,km)
       ql        , &! REAL(kind=r8), INTENT(INOUT) :: ql    (im,km,2)
       q1        , &! REAL(kind=r8), INTENT(INOUT) :: q1    (im,km)
       t1        , &! REAL(kind=r8), INTENT(INOUT) :: t1    (im,km)
       u1        , &! REAL(kind=r8), INTENT(INOUT) :: u1    (im,km)
       v1        , &! REAL(kind=r8), INTENT(INOUT) :: v1    (im,km)
       rcs       , &! REAL(kind=r8), INTENT(IN   ) :: rcs   (im)
       rn        , &! REAL(kind=r8), INTENT(OUT  ) :: rn    (im)
       kbot      , &! INTEGER      , INTENT(INOUT) :: kbot  (im)
       ktop      , &! INTEGER      , INTENT(INOUT) :: ktop  (im)
       kuo       , &! INTEGER      , INTENT(INOUT) :: ktop  (im)
       noshal      , &! INTEGER      , INTENT(INOUT) :: noshal  (im) 
       slimsk    , &! REAL(kind=r8), INTENT(IN   ) :: slimsk(im)
       dot       , &! REAL(kind=r8), INTENT(IN   ) :: dot   (im,km)
       ncloud    , &! INTEGER      , INTENT(IN   ) :: ncloud 
       hpbl      , &! REAL(kind=r8), INTENT(IN   ) :: hpbl  (im)
       heat      , &! REAL(kind=r8), INTENT(IN   ) :: heat  (im)
       evap      , &! REAL(kind=r8), INTENT(IN   ) :: evap  (im)
       ud_mf     , &! REAL(kind=r8), INTENT(OUT  ) :: ud_mf (im,km)
       dt_mf       )! REAL(kind=r8), INTENT(OUT  ) :: dt_mf (im,km)
    !
    IMPLICIT NONE
    !
    INTEGER      , INTENT(IN   ) :: im
    INTEGER      , INTENT(IN   ) :: km
    INTEGER      , INTENT(IN   ) :: jcap
    REAL(kind=r8), INTENT(IN   ) :: delt
    REAL(kind=r8), INTENT(IN   ) :: del   (im,km)
    REAL(kind=r8), INTENT(IN   ) :: prsl  (im,km)
    REAL(kind=r8), INTENT(IN   ) :: ps    (im)
    REAL(kind=r8), INTENT(IN   ) :: phil  (im,km)
    REAL(kind=r8), INTENT(INOUT) :: ql    (im,km,2)
    REAL(kind=r8), INTENT(INOUT) :: q1    (im,km)
    REAL(kind=r8), INTENT(INOUT) :: t1    (im,km)
    REAL(kind=r8), INTENT(INOUT) :: u1    (im,km)
    REAL(kind=r8), INTENT(INOUT) :: v1    (im,km)
    REAL(kind=r8), INTENT(IN   ) :: rcs   (im)
    REAL(kind=r8), INTENT(OUT  ) :: rn    (im)
    INTEGER      , INTENT(INOUT) :: kbot(im)
    INTEGER      , INTENT(INOUT) :: ktop(im)
    INTEGER      , INTENT(in   ) :: kuo    (im)        !   kuo                     ! convection yes(1) or not(0) for shallow convection
    INTEGER      , INTENT(INOUT) :: noshal(im) 
    REAL(kind=r8), INTENT(IN   ) :: slimsk(im)
    REAL(kind=r8), INTENT(IN   ) :: dot   (im,km)!     vvel     - real, layer mean vertical velocity             im,levs !
    INTEGER      , INTENT(IN   ) :: ncloud !     ncld [1 - 0]    - integer, indicate cloud types                     1    !
    REAL(kind=r8), INTENT(IN   ) :: hpbl  (im)  ! hpbl [m]
    REAL(kind=r8), INTENT(IN   ) :: heat  (im)    !     hflx     - real, sensible heat flux                          im   !
    REAL(kind=r8), INTENT(IN   ) :: evap  (im)    !     evap     - real, evaperation from latent heat flux           im   !

    ! hchuang code change mass flux output
    REAL(kind=r8), INTENT(OUT  ) :: ud_mf (im,km)
    REAL(kind=r8), INTENT(OUT  ) :: dt_mf (im,km)
    !
    INTEGER       :: i
    INTEGER       :: indx
!    INTEGER       :: jmn
    INTEGER       :: k
    INTEGER       :: kk
!    INTEGER       :: latd
!    INTEGER       :: lond
    INTEGER       :: km1
    INTEGER       :: kpbl(im)
    !
!    REAL(kind=r8) :: alpha
    REAL(kind=r8) :: alphal
    REAL(kind=r8) :: alphas
    REAL(kind=r8) :: dellat
    REAL(kind=r8) :: desdt
!    REAL(kind=r8) :: deta
!    REAL(kind=r8) :: detad
!    REAL(kind=r8) :: dg
!    REAL(kind=r8) :: dh
!    REAL(kind=r8) :: dhh
!    REAL(kind=r8) :: dlnsig
    REAL(kind=r8) :: dp
    REAL(kind=r8) :: dq
    REAL(kind=r8) :: dqsdp
    REAL(kind=r8) :: dqsdt
    REAL(kind=r8) :: dt
    REAL(kind=r8) :: dt2
    REAL(kind=r8) :: dtmax
    REAL(kind=r8) :: dtmin
    REAL(kind=r8) :: dv1h
    REAL(kind=r8) :: dv1q
    REAL(kind=r8) :: dv2h
    REAL(kind=r8) :: dv2q
    REAL(kind=r8) :: dv1u
    REAL(kind=r8) :: dv1v
    REAL(kind=r8) :: dv2u
    REAL(kind=r8) :: dv2v
    REAL(kind=r8) :: dv3q
    REAL(kind=r8) :: dv3h
    REAL(kind=r8) :: dv3u
    REAL(kind=r8) :: dv3v
    REAL(kind=r8) :: clam
    REAL(kind=r8) :: dz
    REAL(kind=r8) :: e1
    REAL(kind=r8) :: es
    REAL(kind=r8) :: etah
    REAL(kind=r8) :: evef
    REAL(kind=r8) :: evfact
    REAL(kind=r8) :: evfactl
    REAL(kind=r8) :: factor
    REAL(kind=r8) :: fjcap
    REAL(kind=r8) :: gamma
    REAL(kind=r8) :: pprime
    REAL(kind=r8) :: betaw
    REAL(kind=r8) :: qc
    REAL(kind=r8) :: qlk
    REAL(kind=r8) :: qrch
    REAL(kind=r8) :: qs
!    REAL(kind=r8) ::  rain
!    REAL(kind=r8) :: rfact
    REAL(kind=r8) :: shear
    REAL(kind=r8) :: tem1
!    REAL(kind=r8) :: tem2
    REAL(kind=r8) :: val
    REAL(kind=r8) :: val1
    REAL(kind=r8) :: val2
    REAL(kind=r8) :: w1
    REAL(kind=r8) :: w1l
    REAL(kind=r8) :: w1s
    REAL(kind=r8) :: w2
    REAL(kind=r8) :: w2l
    REAL(kind=r8) :: w2s
    REAL(kind=r8) :: w3
    REAL(kind=r8) :: w3l
    REAL(kind=r8) :: w3s
    REAL(kind=r8) :: w4
    REAL(kind=r8) :: w4l
    REAL(kind=r8) :: w4s
    REAL(kind=r8) :: tem
    REAL(kind=r8) :: ptem
    REAL(kind=r8) :: ptem1
    REAL(kind=r8) :: pgcon
    !
    INTEGER       :: kb(im)
    INTEGER       :: kbcon(im)
    INTEGER       :: kbcon1(im)
    INTEGER       :: ktcon(im)
    INTEGER       :: kbm(im)
    INTEGER       :: kmax(im)
    !
    REAL(kind=r8) :: delhbar(im)
    REAL(kind=r8) :: delq(im)
    REAL(kind=r8) :: delq2(im)

    REAL(kind=r8) :: delqbar(im)
    REAL(kind=r8) :: delqev(im)
    REAL(kind=r8) :: deltbar(im)

    REAL(kind=r8) :: deltv(im)
    REAL(kind=r8) :: edt(im)

    REAL(kind=r8) :: wstar(im)
    REAL(kind=r8) :: sflx(im)

    REAL(kind=r8) :: pdot(im)
    REAL(kind=r8) :: po(im,km)

    REAL(kind=r8) :: qcond(im)
    REAL(kind=r8) :: qevap(im)
    REAL(kind=r8) :: hmax(im)

    REAL(kind=r8) :: rntot(im)
    REAL(kind=r8) :: vshear(im)

!    REAL(kind=r8) :: xk(im)
    REAL(kind=r8) :: xlamud(im)

    REAL(kind=r8) :: xmb(im)
    REAL(kind=r8) :: xmbmax(im)

    REAL(kind=r8) :: delubar(im)
    REAL(kind=r8) :: delvbar(im)

    REAL(kind=r8) :: cincr
    !1
    !   physical parameters

    REAL(kind=r8), PARAMETER :: c0=0.002_r8
!    REAL(kind=r8), PARAMETER :: cpoel=con_cp/con_hvap
!    REAL(kind=r8), PARAMETER :: delta=con_fvirt
    REAL(kind=r8), PARAMETER :: c1=3.e-4_r8
    REAL(kind=r8), PARAMETER :: el2orc=con_hvap*con_hvap/(con_rv*con_cp)
    REAL(kind=r8), PARAMETER :: elocp=con_hvap/con_cp
    REAL(kind=r8), PARAMETER :: h1=0.33333333_r8
    REAL(kind=r8), PARAMETER :: dthk=25.0_r8
    REAL(kind=r8), PARAMETER :: g =con_g
    REAL(kind=r8), PARAMETER :: fact1 =(con_cvap-con_cliq)/con_rv
    REAL(kind=r8), PARAMETER :: fact2 =con_hvap/con_rv-fact1*con_t0c
!    REAL(kind=r8), PARAMETER :: terr=0.0_r8
    REAL(kind=r8), PARAMETER :: cincrmax=180.0_r8
    REAL(kind=r8), PARAMETER :: cincrmin=120.0_r8




    !  local variables and arrays
    REAL(kind=r8) :: pfld(im,km)
    REAL(kind=r8) :: to(im,km)
    REAL(kind=r8) :: qo(im,km)
    REAL(kind=r8) :: uo(im,km)
    REAL(kind=r8) :: vo(im,km)
    REAL(kind=r8) :: qeso(im,km)
    !  cloud water
    !     real(kind=r8) qlko_ktcon(im), dellal(im,km), tvo(im,km),
    REAL(kind=r8) :: qlko_ktcon(im)
    REAL(kind=r8) :: dellal(im,km)
    REAL(kind=r8) :: dbyo(im,km)
    REAL(kind=r8) :: zo(im,km)
    REAL(kind=r8) :: xlamue(im,km)
    REAL(kind=r8) :: heo(im,km)
    REAL(kind=r8) :: heso(im,km)
    REAL(kind=r8) :: dellah(im,km)
    REAL(kind=r8) :: dellaq(im,km)

    REAL(kind=r8) :: dellau(im,km)
    REAL(kind=r8) :: dellav(im,km)
    REAL(kind=r8) :: hcko(im,km)
    REAL(kind=r8) :: ucko(im,km)
    REAL(kind=r8) :: vcko(im,km)
    REAL(kind=r8) :: qcko(im,km)

    REAL(kind=r8) :: eta(im,km)
    REAL(kind=r8) :: zi(im,km)
    REAL(kind=r8) :: pwo(im,km)

    REAL(kind=r8) :: tx1(im)
    !
    LOGICAL       :: totflg
    LOGICAL       :: cnvflg(im)
    LOGICAL       :: flg(im)
    !
    REAL(kind=r8), PARAMETER :: tf =233.16_r8
    REAL(kind=r8), PARAMETER :: tcr=263.16_r8
    REAL(kind=r8), PARAMETER :: tcrf=1.0_r8/(tcr-tf)
    !      parameter (tf=233.16_r8, tcr=263.16_r8, tcrf=1.0_r8/(tcr-tf))
    !
    !-----------------------------------------------------------------------
    !
    km1 = km - 1
    !
    !  compute surface buoyancy flux
    !
    DO i=1,im
       sflx(i) = heat(i)+con_fvirt*t1(i,1)*evap(i)
    ENDDO
    !
    !  initialize arrays
    !
    DO i=1,im
       noshal(i)=0
       IF(kuo(i) .EQ. 1) THEN
          noshal(i)=1             !   1 do not  shallow convection
       END IF
    END DO

    DO i=1,im
       cnvflg(i) = .TRUE.
       IF(kuo(i).EQ.1) cnvflg(i) = .FALSE.  !   deep convection yes(1) or not(0) for shallow convection
       IF(sflx(i).LE.0.0_r8) cnvflg(i) = .FALSE.
       IF(cnvflg(i)) THEN
          kbot(i)=km+1
          ktop(i)=0
       ENDIF
       rn(i)=0.0_r8
       kbcon(i)=km
       ktcon(i)=1
       kb(i)=km
       pdot(i) = 0.0_r8
       qlko_ktcon(i) = 0.0_r8
       edt(i)  = 0.0_r8
       vshear(i) = 0.0_r8
    ENDDO
    ! hchuang code change
    DO k = 1, km
       DO i = 1, im
          ud_mf(i,k) = 0.0_r8
          dt_mf(i,k) = 0.0_r8
       ENDDO
    ENDDO
    !!
    totflg = .TRUE.
    DO i=1,im
       totflg = totflg .AND. (.NOT. cnvflg(i))
    ENDDO
    IF(totflg) RETURN
    !!
    !
    dt2   = delt
    val   =         1200.0_r8
    dtmin = MAX(dt2, val )
    val   =         3600.0_r8
    dtmax = MAX(dt2, val )
    !  model tunable parameters are all here
    alphal  = 0.5_r8
    alphas  = 0.5_r8
    clam    = 0.3_r8
    betaw   = 0.03_r8
    !     evef    = 0.07_r8
    evfact  = 0.3_r8
    evfactl = 0.3_r8
    !
    !     pgcon   = 0.7_r8     ! Gregory et al. (1997, QJRMS)
    pgcon   = 0.55_r8    ! Zhang & Wu (2003,JAS)
    fjcap   = (float(jcap) / 126.0_r8) ** 2
    val     =           1.0_r8
    fjcap   = MAX(fjcap,val)
    w1l     = -8.e-3_r8 
    w2l     = -4.e-2_r8
    w3l     = -5.e-3_r8 
    w4l     = -5.e-4_r8
    w1s     = -2.e-4_r8
    w2s     = -2.e-3_r8
    w3s     = -1.e-3_r8
    w4s     = -2.e-5_r8
    !
    !  define top layer for search of the downdraft originating layer
    !  and the maximum thetae for updraft
    !
    DO i=1,im
       kbm(i)   = km
       kmax(i)  = km
       tx1(i)   = 1.0_r8 / ps(i)
    ENDDO
    !     
    DO k = 1, km
       DO i=1,im
          IF (prsl(i,k)*tx1(i) .GT. 0.70_r8) kbm(i)   = k + 1
          IF (prsl(i,k)*tx1(i) .GT. 0.60_r8) kmax(i)  = k + 1
       ENDDO
    ENDDO
    DO i=1,im
       kbm(i)   = MIN(kbm(i),kmax(i))
    ENDDO
    !
    !  hydrostatic height assume zero terr and compute
    !  updraft entrainment rate as an inverse function of height
    !
    DO k = 1, km
       DO i=1,im
          zo(i,k) = phil(i,k) / g
       ENDDO
    ENDDO
    DO k = 1, km1
       DO i=1,im
          zi(i,k) = 0.5_r8*(zo(i,k)+zo(i,k+1))
          xlamue(i,k) = clam / zi(i,k)
       ENDDO
    ENDDO
    DO i=1,im
       xlamue(i,km) = xlamue(i,km1)
    ENDDO
    !
    !  pbl height
    !
    DO i=1,im
       flg(i) = cnvflg(i)
       kpbl(i)= 1
    ENDDO
    DO k = 2, km1
       DO i=1,im
          IF (flg(i).AND.zo(i,k).LE.hpbl(i)) THEN
             kpbl(i) = k
          ELSE
             flg(i) = .FALSE.
          ENDIF
       ENDDO
    ENDDO
    DO i=1,im
       kpbl(i)= MIN(kpbl(i),kbm(i))
    ENDDO
    !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !   convert surface pressure to mb from cb
    !
    DO k = 1, km
       DO i = 1, im
          IF (cnvflg(i) .AND. k .LE. kmax(i)) THEN
             pfld(i,k) = prsl(i,k) * 10.0_r8
             eta(i,k)  = 1.0_r8
             hcko(i,k) = 0.0_r8
             qcko(i,k) = 0.0_r8
             ucko(i,k) = 0.0_r8
             vcko(i,k) = 0.0_r8
             dbyo(i,k) = 0.0_r8
             pwo(i,k)  = 0.0_r8
             dellal(i,k) = 0.0_r8
             to(i,k)   = t1(i,k)
             qo(i,k)   = q1(i,k)
             uo(i,k)   = u1(i,k) * rcs(i)
             vo(i,k)   = v1(i,k) * rcs(i)
          ENDIF
       ENDDO
    ENDDO
    !
    !  column variables
    !  p is pressure of the layer (mb)
    !  t is temperature at t-dt (k)..tn
    !  q is mixing ratio at t-dt (kg/kg)..qn
    !  to is temperature at t+dt (k)... this is after advection and turbulan
    !  qo is mixing ratio at t+dt (kg/kg)..q1
    !
    DO k = 1, km
       DO i=1,im
          IF (cnvflg(i) .AND. k .LE. kmax(i)) THEN
             qeso(i,k) = 0.01_r8 * fpvs(to(i,k))      ! fpvs is in pa
             qeso(i,k) = con_eps * qeso(i,k) / (pfld(i,k) + con_epsm1*qeso(i,k))
             val1      =             1.e-8_r8
             qeso(i,k) = MAX(qeso(i,k), val1)
             val2      =           1.e-10_r8
             qo(i,k)   = MAX(qo(i,k), val2 )
             !           qo(i,k)   = min(qo(i,k),qeso(i,k))
             !           tvo(i,k)  = to(i,k) + delta * to(i,k) * qo(i,k)
          ENDIF
       ENDDO
    ENDDO
    !
    !  compute moist static energy
    !
    DO k = 1, km
       DO i=1,im
          IF (cnvflg(i) .AND. k .LE. kmax(i)) THEN
             !           tem       = g * zo(i,k) + con_cp * to(i,k)
             tem       = phil(i,k) + con_cp * to(i,k)
             heo(i,k)  = tem  + con_hvap * qo(i,k)
             heso(i,k) = tem  + con_hvap * qeso(i,k)
             !           heo(i,k)  = min(heo(i,k),heso(i,k))
          ENDIF
       ENDDO
    ENDDO
    !
    !  determine level with largest moist static energy within pbl
    !  this is the level where updraft starts
    !
    DO i=1,im
       IF (cnvflg(i)) THEN
          hmax(i) = heo(i,1)
          kb(i) = 1
       ENDIF
    ENDDO
    DO k = 2, km
       DO i=1,im
          IF (cnvflg(i).AND.k.LE.kpbl(i)) THEN
             IF(heo(i,k).GT.hmax(i)) THEN
                kb(i)   = k
                hmax(i) = heo(i,k)
             ENDIF
          ENDIF
       ENDDO
    ENDDO
    !
    DO k = 1, km1
       DO i=1,im
          IF (cnvflg(i) .AND. k .LE. kmax(i)-1) THEN
             dz      = 0.5_r8 * (zo(i,k+1) - zo(i,k))
             dp      = 0.5_r8 * (pfld(i,k+1) - pfld(i,k))
             es      = 0.01_r8 * fpvs(to(i,k+1))      ! fpvs is in pa
             pprime  = pfld(i,k+1) + con_epsm1 * es
             qs      = con_eps * es / pprime
             dqsdp   = - qs / pprime
             desdt   = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
             dqsdt   = qs * pfld(i,k+1) * desdt / (es * pprime)
             gamma   = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
             dt      = (g * dz + con_hvap * dqsdp * dp) / (con_cp * (1.0_r8 + gamma))
             dq      = dqsdt * dt + dqsdp * dp
             to(i,k) = to(i,k+1) + dt
             qo(i,k) = qo(i,k+1) + dq
             po(i,k) = 0.5_r8 * (pfld(i,k) + pfld(i,k+1))
          ENDIF
       ENDDO
    ENDDO
    !
    DO k = 1, km1
       DO i=1,im
          IF (cnvflg(i) .AND. k .LE. kmax(i)-1) THEN
             qeso(i,k) = 0.01_r8 * fpvs(to(i,k))      ! fpvs is in pa
             qeso(i,k) = con_eps * qeso(i,k) / (po(i,k) + con_epsm1*qeso(i,k))
             val1      =             1.e-8_r8
             qeso(i,k) = MAX(qeso(i,k), val1)
             val2      =           1.e-10_r8
             qo(i,k)   = MAX(qo(i,k), val2 )
             !           qo(i,k)   = min(qo(i,k),qeso(i,k))
             heo(i,k)  = 0.5_r8 * g * (zo(i,k) + zo(i,k+1)) +   &
                  con_cp * to(i,k) + con_hvap * qo(i,k)
             heso(i,k) = 0.5_r8 * g * (zo(i,k) + zo(i,k+1)) +   &
                  con_cp * to(i,k) + con_hvap * qeso(i,k)
             uo(i,k)   = 0.5_r8 * (uo(i,k) + uo(i,k+1))
             vo(i,k)   = 0.5_r8 * (vo(i,k) + vo(i,k+1))
          ENDIF
       ENDDO
    ENDDO
    !
    !  look for the level of free convection as cloud base
    !
    DO i=1,im
       flg(i)   = cnvflg(i)
       IF(flg(i)) kbcon(i) = kmax(i)
    ENDDO
    DO k = 1, km1
       DO i=1,im
          IF (flg(i).AND.k.LT.kbm(i)) THEN
             IF(k.GT.kb(i).AND.heo(i,kb(i)).GT.heso(i,k)) THEN
                kbcon(i) = k
                flg(i)   = .FALSE.
             ENDIF
          ENDIF
       ENDDO
    ENDDO
    !
    DO i=1,im
       IF(cnvflg(i)) THEN
          IF(kbcon(i).EQ.kmax(i)) cnvflg(i) = .FALSE.
       ENDIF
    ENDDO
    !!
    totflg = .TRUE.
    DO i=1,im
       totflg = totflg .AND. (.NOT. cnvflg(i))
    ENDDO
    IF(totflg) RETURN
    !!
    !
    !  determine critical convective inhibition
    !  as a function of vertical velocity at cloud base.
    !
    DO i=1,im
       IF(cnvflg(i)) THEN
          pdot(i)  = 10.0_r8* dot(i,kbcon(i))
       ENDIF
    ENDDO
    DO i=1,im
       IF(cnvflg(i)) THEN
          IF(slimsk(i).EQ.1.0_r8) THEN
             w1 = w1l
             w2 = w2l
             w3 = w3l
             w4 = w4l
          ELSE
             w1 = w1s
             w2 = w2s
             w3 = w3s
             w4 = w4s
          ENDIF
          IF(pdot(i).LE.w4) THEN
             ptem = (pdot(i) - w4) / (w3 - w4)
          ELSEIF(pdot(i).GE.-w4) THEN
             ptem = - (pdot(i) + w4) / (w4 - w3)
          ELSE
             ptem = 0.0_r8
          ENDIF
          val1    =             -1.0_r8
          ptem = MAX(ptem,val1)
          val2    =             1.0_r8
          ptem = MIN(ptem,val2)
          ptem = 1.0_r8 - ptem
          ptem1= 0.5_r8*(cincrmax-cincrmin)
          cincr = cincrmax - ptem * ptem1
          tem1 = pfld(i,kb(i)) - pfld(i,kbcon(i))
          IF(tem1.GT.cincr) THEN
             cnvflg(i) = .FALSE.
          ENDIF
       ENDIF
    ENDDO
    !!
    totflg = .TRUE.
    DO i=1,im
       totflg = totflg .AND. (.NOT. cnvflg(i))
    ENDDO
    IF(totflg) RETURN
    !!
    !
    !  assume the detrainment rate for the updrafts to be same as 
    !  the entrainment rate at cloud base
    !
    DO i = 1, im
       IF(cnvflg(i)) THEN
          xlamud(i) = xlamue(i,kbcon(i))
       ENDIF
    ENDDO
    !
    !  determine updraft mass flux for the subcloud layers
    !
    DO k = km1, 1, -1
       DO i = 1, im
          IF (cnvflg(i)) THEN
             IF(k.LT.kbcon(i).AND.k.GE.kb(i)) THEN
                dz       = zi(i,k+1) - zi(i,k)
                ptem     = 0.5_r8*(xlamue(i,k)+xlamue(i,k+1))-xlamud(i)
                eta(i,k) = eta(i,k+1) / (1.0_r8 + ptem * dz)
             ENDIF
          ENDIF
       ENDDO
    ENDDO
    !
    !  compute updraft cloud property
    !
    DO i = 1, im
       IF(cnvflg(i)) THEN
          indx         = kb(i)
          hcko(i,indx) = heo(i,indx)
          qcko(i,indx) = qo(i,indx)
          ucko(i,indx) = uo(i,indx)
          vcko(i,indx) = vo(i,indx)
       ENDIF
    ENDDO
    !
    DO k = 2, km1
       DO i = 1, im
          IF (cnvflg(i)) THEN
             IF(k.GT.kb(i).AND.k.LT.kmax(i)) THEN
                dz   = zi(i,k) - zi(i,k-1)
                tem  = 0.5_r8 * (xlamue(i,k)+xlamue(i,k-1)) * dz
                tem1 = 0.5_r8 * xlamud(i) * dz
                factor = 1.0_r8 + tem - tem1
                ptem = 0.5_r8 * tem + pgcon
                ptem1= 0.5_r8 * tem - pgcon
                hcko(i,k) = ((1.0_r8-tem1)*hcko(i,k-1)+tem*0.5_r8*   &
                     (heo(i,k)+heo(i,k-1)))/factor
                ucko(i,k) = ((1.0_r8-tem1)*ucko(i,k-1)+ptem*uo(i,k)   &
                     +ptem1*uo(i,k-1))/factor
                vcko(i,k) = ((1.0_r8-tem1)*vcko(i,k-1)+ptem*vo(i,k)   &
                     +ptem1*vo(i,k-1))/factor
                dbyo(i,k) = hcko(i,k) - heso(i,k)
             ENDIF
          ENDIF
       ENDDO
    ENDDO
    !
    !   taking account into convection inhibition due to existence of
    !    dry layers below cloud base
    !
    DO i=1,im
       flg(i) = cnvflg(i)
       kbcon1(i) = kmax(i)
    ENDDO
    DO k = 2, km1
       DO i=1,im
          IF (flg(i).AND.k.LT.kbm(i)) THEN
             IF(k.GE.kbcon(i).AND.dbyo(i,k).GT.0.0_r8) THEN
                kbcon1(i) = k
                flg(i)    = .FALSE.
             ENDIF
          ENDIF
       ENDDO
    ENDDO
    DO i=1,im
       IF(cnvflg(i)) THEN
          IF(kbcon1(i).EQ.kmax(i)) cnvflg(i) = .FALSE.
       ENDIF
    ENDDO
    DO i=1,im
       IF(cnvflg(i)) THEN
          tem = pfld(i,kbcon(i)) - pfld(i,kbcon1(i))
          IF(tem.GT.dthk) THEN
             cnvflg(i) = .FALSE.
          ENDIF
       ENDIF
    ENDDO
    !!
    totflg = .TRUE.
    DO i = 1, im
       totflg = totflg .AND. (.NOT. cnvflg(i))
    ENDDO
    IF(totflg) RETURN
    !!
    !
    !  determine convective cloud top as the level of zero buoyancy
    !   but limited to the level of 700 mb
    !
    DO i = 1, im
       flg(i) = cnvflg(i)
       IF(flg(i)) ktcon(i) = kbm(i)
    ENDDO
    DO k = 2, km1
       DO i=1,im
          IF (flg(i).AND.k .LT. kbm(i)) THEN
             IF(k.GT.kbcon1(i).AND.dbyo(i,k).LT.0.0_r8) THEN
                ktcon(i) = k
                flg(i)   = .FALSE.
             ENDIF
          ENDIF
       ENDDO
    ENDDO
    !
    !  turn off shallow convection if cloud top is less than pbl top
    !
    DO i=1,im
       IF(cnvflg(i)) THEN
          kk = kpbl(i)+1
          IF(ktcon(i).LE.kk) cnvflg(i) = .FALSE.
       ENDIF
    ENDDO
    !!
    totflg = .TRUE.
    DO i = 1, im
       totflg = totflg .AND. (.NOT. cnvflg(i))
    ENDDO
    IF(totflg) RETURN
    !!
    !
    !  compute updraft mass flux for the cloud layers
    !
    DO k = 2, km1
       DO i = 1, im
          IF(cnvflg(i)) THEN
             IF(k.GT.kbcon(i).AND.k.LE.ktcon(i)) THEN
                dz       = zi(i,k) - zi(i,k-1)
                ptem     = 0.5_r8*(xlamue(i,k)+xlamue(i,k-1))-xlamud(i)
                eta(i,k) = eta(i,k-1) * (1 + ptem * dz)
             ENDIF
          ENDIF
       ENDDO
    ENDDO
    !
    !  specify upper limit of mass flux at cloud base
    !
    DO i = 1, im
       IF(cnvflg(i)) THEN
          !         xmbmax(i) = .1_r8
          !
          k = kbcon(i)
          dp = 1000.0_r8 * del(i,k)
          xmbmax(i) = dp / (g * dt2)
          !
          !         tem = dp / (g * dt2)
          !         xmbmax(i) = min(tem, xmbmax(i))
       ENDIF
    ENDDO
    !
    !  compute cloud moisture property and precipitation
    !
    DO k = 2, km
       DO i = 1, im
          IF (cnvflg(i)) THEN
             IF(k.GT.kb(i).AND.k.LT.ktcon(i)) THEN
                dz    = zi(i,k) - zi(i,k-1)
                gamma = el2orc * qeso(i,k) / (to(i,k)**2)
                qrch = qeso(i,k)    &
                     + gamma * dbyo(i,k) / (con_hvap * (1.0_r8 + gamma))
                !j
                tem  = 0.5_r8 * (xlamue(i,k)+xlamue(i,k-1)) * dz
                tem1 = 0.5_r8 * xlamud(i) * dz
                factor = 1.0_r8 + tem - tem1
                qcko(i,k) = ((1.0_r8-tem1)*qcko(i,k-1)+tem*0.5_r8*   &
                     (qo(i,k)+qo(i,k-1)))/factor
                !j
                dq = eta(i,k) * (qcko(i,k) - qrch)
                !
                !             rhbar(i) = rhbar(i) + qo(i,k) / qeso(i,k)
                !
                !  below lfc check if there is excess moisture to release latent heat
                !
                IF(dq.GT.0.0_r8) THEN
                   dp = 1000.0_r8 * del(i,k)
                   etah = 0.5_r8 * (eta(i,k) + eta(i,k-1))
                   qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
                   qc = qlk + qrch
                   pwo(i,k) = etah * c0 * dz * qlk
                   dellal(i,k) = etah * c1 * dz * qlk * g / dp
                   qcko(i,k)= qc
                ENDIF
             ENDIF
          ENDIF
       ENDDO
    ENDDO
    !
    !  this section is ready for cloud water
    !
    IF(ncloud.GT.0) THEN
       !
       !  compute liquid and vapor separation at cloud top
       !
       DO i = 1, im
          IF(cnvflg(i)) THEN
             k = ktcon(i)
             gamma = el2orc * qeso(i,k) / (to(i,k)**2)
             qrch = qeso(i,k) &
                  + gamma * dbyo(i,k) / (con_hvap * (1.0_r8 + gamma))
             dq = qcko(i,k-1) - qrch
             !
             !  check if there is excess moisture to release latent heat
             !
             IF(dq.GT.0.0_r8) THEN
                qlko_ktcon(i) = dq
                qcko(i,k-1) = qrch
             ENDIF
          ENDIF
       ENDDO
    ENDIF
    !
    !--- compute precipitation efficiency in terms of windshear
    !
    DO i = 1, im
       IF(cnvflg(i)) THEN
          vshear(i) = 0.0_r8
       ENDIF
    ENDDO
    DO k = 2, km
       DO i = 1, im
          IF (cnvflg(i)) THEN
             IF(k.GT.kb(i).AND.k.LE.ktcon(i)) THEN
                shear= SQRT((uo(i,k)-uo(i,k-1)) ** 2   &
                     + (vo(i,k)-vo(i,k-1)) ** 2)
                vshear(i) = vshear(i) + shear
             ENDIF
          ENDIF
       ENDDO
    ENDDO
    DO i = 1, im
       IF(cnvflg(i)) THEN
          vshear(i) = 1.e3_r8 * vshear(i) / (zi(i,ktcon(i))-zi(i,kb(i)))
          e1=1.591_r8-0.639_r8*vshear(i)   &
               +0.0953_r8*(vshear(i)**2)-0.00496_r8*(vshear(i)**3)
          edt(i)=1.0_r8-e1
          val =         0.9_r8
          edt(i) = MIN(edt(i),val)
          val =         0.0_r8
          edt(i) = MAX(edt(i),val)
       ENDIF
    ENDDO
    !
    !--- what would the change be, that a cloud with unit mass
    !--- will do to the environment?
    !
    DO k = 1, km
       DO i = 1, im
          IF(cnvflg(i) .AND. k .LE. kmax(i)) THEN
             dellah(i,k) = 0.0_r8
             dellaq(i,k) = 0.0_r8
             dellau(i,k) = 0.0_r8
             dellav(i,k) = 0.0_r8
          ENDIF
       ENDDO
    ENDDO
    !
    !--- changed due to subsidence and entrainment
    !
    DO k = 2, km1
       DO i = 1, im
          IF (cnvflg(i)) THEN
             IF(k.GT.kb(i).AND.k.LT.ktcon(i)) THEN
                dp = 1000.0_r8 * del(i,k)
                dz = zi(i,k) - zi(i,k-1)
                !
                dv1h = heo(i,k)
                dv2h = 0.5_r8 * (heo(i,k) + heo(i,k-1))
                dv3h = heo(i,k-1)
                dv1q = qo(i,k)
                dv2q = 0.5_r8 * (qo(i,k) + qo(i,k-1))
                dv3q = qo(i,k-1)
                dv1u = uo(i,k)
                dv2u = 0.5_r8 * (uo(i,k) + uo(i,k-1))
                dv3u = uo(i,k-1)
                dv1v = vo(i,k)
                dv2v = 0.5_r8 * (vo(i,k) + vo(i,k-1))
                dv3v = vo(i,k-1)
                !
                tem  = 0.5_r8 * (xlamue(i,k)+xlamue(i,k-1))
                tem1 = xlamud(i)
                !j
                dellah(i,k) = dellah(i,k) +                   &
                     &     ( eta(i,k)*dv1h - eta(i,k-1)*dv3h               &
                     &    -  tem*eta(i,k-1)*dv2h*dz                        &
                     &    +  tem1*eta(i,k-1)*0.5_r8*(hcko(i,k)+hcko(i,k-1))*dz &
                     &         ) *g/dp
                !j
                dellaq(i,k) = dellaq(i,k) +                   &
                     &     ( eta(i,k)*dv1q - eta(i,k-1)*dv3q               &
                     &    -  tem*eta(i,k-1)*dv2q*dz                        &
                     &    +  tem1*eta(i,k-1)*0.5_r8*(qcko(i,k)+qcko(i,k-1))*dz &
                     &         ) *g/dp
                !j
                dellau(i,k) = dellau(i,k) +                   &
                     &     ( eta(i,k)*dv1u - eta(i,k-1)*dv3u               &
                     &    -  tem*eta(i,k-1)*dv2u*dz                        &
                     &    +  tem1*eta(i,k-1)*0.5_r8*(ucko(i,k)+ucko(i,k-1))*dz &
                     &    -  pgcon*eta(i,k-1)*(dv1u-dv3u)                  &
                     &         ) *g/dp
                !j
                dellav(i,k) = dellav(i,k) +                   &
                     &     ( eta(i,k)*dv1v - eta(i,k-1)*dv3v               &
                     &    -  tem*eta(i,k-1)*dv2v*dz                        &
                     &    +  tem1*eta(i,k-1)*0.5_r8*(vcko(i,k)+vcko(i,k-1))*dz &
                     &    -  pgcon*eta(i,k-1)*(dv1v-dv3v)                  &
                     &         ) *g/dp
                !j
             ENDIF
          ENDIF
       ENDDO
    ENDDO
    !
    !------- cloud top
    !
    DO i = 1, im
       IF(cnvflg(i)) THEN
          indx = ktcon(i)
          dp = 1000.0_r8 * del(i,indx)
          dv1h = heo(i,indx-1)
          dellah(i,indx) = eta(i,indx-1) *                     &
               &                     (hcko(i,indx-1) - dv1h) * g / dp
          dv1q = qo(i,indx-1)
          dellaq(i,indx) = eta(i,indx-1) *                     &
               &                     (qcko(i,indx-1) - dv1q) * g / dp
          dv1u = uo(i,indx-1)
          dellau(i,indx) = eta(i,indx-1) *                     &
               &                     (ucko(i,indx-1) - dv1u) * g / dp
          dv1v = vo(i,indx-1)
          dellav(i,indx) = eta(i,indx-1) *                     &
               &                     (vcko(i,indx-1) - dv1v) * g / dp
          !
          !  cloud water
          !
          dellal(i,indx) = eta(i,indx-1) *                     &
               &                     qlko_ktcon(i) * g / dp
       ENDIF
    ENDDO
    !
    !  mass flux at cloud base for shallow convection
    !  (Grant, 2001)
    !
    DO i= 1, im
       IF(cnvflg(i)) THEN
          k = kbcon(i)
          !         ptem = g*sflx(i)*zi(i,k)/t1(i,1)
          ptem = g*sflx(i)*hpbl(i)/t1(i,1)
          wstar(i) = ptem**h1
          tem = po(i,k)*100.0_r8 / (con_rd*t1(i,k))
          xmb(i) = betaw*tem*wstar(i)
          xmb(i) = MIN(xmb(i),xmbmax(i))
       ENDIF
    ENDDO
    !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !
    DO k = 1, km
       DO i = 1, im
          IF (cnvflg(i) .AND. k .LE. kmax(i)) THEN
             qeso(i,k) = 0.01_r8 * fpvs(t1(i,k))      ! fpvs is in pa
             qeso(i,k) = con_eps * qeso(i,k) / (pfld(i,k) + con_epsm1*qeso(i,k))
             val     =             1.e-8_r8
             qeso(i,k) = MAX(qeso(i,k), val )
          ENDIF
       ENDDO
    ENDDO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !
    DO i = 1, im
       delhbar(i) = 0.0_r8
       delqbar(i) = 0.0_r8
       deltbar(i) = 0.0_r8
       delubar(i) = 0.0_r8
       delvbar(i) = 0.0_r8
       qcond(i) = 0.0_r8
    ENDDO
    DO k = 1, km
       DO i = 1, im
          IF (cnvflg(i)) THEN
             IF(k.GT.kb(i).AND.k.LE.ktcon(i)) THEN
                dellat = (dellah(i,k) - con_hvap * dellaq(i,k)) / con_cp
                t1(i,k) = t1(i,k) + dellat * xmb(i) * dt2
                q1(i,k) = q1(i,k) + dellaq(i,k) * xmb(i) * dt2
                tem = 1.0_r8/rcs(i)
                u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2 * tem
                v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2 * tem
                dp = 1000.0_r8 * del(i,k)
                delhbar(i) = delhbar(i) + dellah(i,k)*xmb(i)*dp/g
                delqbar(i) = delqbar(i) + dellaq(i,k)*xmb(i)*dp/g
                deltbar(i) = deltbar(i) + dellat*xmb(i)*dp/g
                delubar(i) = delubar(i) + dellau(i,k)*xmb(i)*dp/g
                delvbar(i) = delvbar(i) + dellav(i,k)*xmb(i)*dp/g
             ENDIF
          ENDIF
       ENDDO
    ENDDO
    DO k = 1, km
       DO i = 1, im
          IF (cnvflg(i)) THEN
             IF(k.GT.kb(i).AND.k.LE.ktcon(i)) THEN
                qeso(i,k) = 0.01_r8 * fpvs(t1(i,k))      ! fpvs is in pa
                qeso(i,k) = con_eps * qeso(i,k)/(pfld(i,k) + con_epsm1*qeso(i,k))
                val     =             1.e-8_r8
                qeso(i,k) = MAX(qeso(i,k), val )
             ENDIF
          ENDIF
       ENDDO
    ENDDO
    !
    DO i = 1, im
       rntot(i) = 0.0_r8
       delqev(i) = 0.0_r8
       delq2(i) = 0.0_r8
       flg(i) = cnvflg(i)
    ENDDO
    DO k = km, 1, -1
       DO i = 1, im
          IF (cnvflg(i)) THEN
             IF(k.LT.ktcon(i).AND.k.GT.kb(i)) THEN
                rntot(i) = rntot(i) + pwo(i,k) * xmb(i) * 0.001_r8 * dt2
             ENDIF
          ENDIF
       ENDDO
    ENDDO
    !
    ! evaporating rain
    !
    DO k = km, 1, -1
       DO i = 1, im
          IF (k .LE. kmax(i)) THEN
             deltv(i) = 0.0_r8
             delq(i) = 0.0_r8
             qevap(i) = 0.0_r8
             IF(cnvflg(i)) THEN
                IF(k.LT.ktcon(i).AND.k.GT.kb(i)) THEN
                   rn(i) = rn(i) + pwo(i,k) * xmb(i) * 0.001_r8 * dt2
                ENDIF
             ENDIF
             IF(flg(i).AND.k.LT.ktcon(i)) THEN
                evef = edt(i) * evfact
                IF(slimsk(i).EQ.1.0_r8) evef=edt(i) * evfactl
                !             if(slimsk(i).eq.1.0_r8) evef=.07_r8
                !             if(slimsk(i).ne.1.0_r8) evef = 0.0_r8
                qcond(i) = evef * (q1(i,k) - qeso(i,k))             &
                     / (1.0_r8 + el2orc * qeso(i,k) / t1(i,k)**2)
                dp = 1000.0_r8 * del(i,k)
                IF(rn(i).GT.0.0_r8.AND.qcond(i).LT.0.0_r8) THEN
                   qevap(i) = -qcond(i) * (1.0_r8-EXP(-0.32_r8*SQRT(dt2*rn(i))))
                   qevap(i) = MIN(qevap(i), rn(i)*1000.0_r8*con_g/dp)
                   delq2(i) = delqev(i) + 0.001_r8 * qevap(i) * dp / g
                ENDIF
                IF(rn(i).GT.0.0_r8.AND.qcond(i).LT.0.0_r8.AND.       &
                     delq2(i).GT.rntot(i)) THEN
                   qevap(i) = 1000.0_r8* g * (rntot(i) - delqev(i)) / dp
                   flg(i) = .FALSE.
                ENDIF
                IF(rn(i).GT.0.0_r8.AND.qevap(i).GT.0.0_r8) THEN
                   tem  = 0.001_r8 * dp / g
                   tem1 = qevap(i) * tem
                   IF(tem1.GT.rn(i)) THEN
                      qevap(i) = rn(i) / tem
                      rn(i) = 0.0_r8
                   ELSE
                      rn(i) = rn(i) - tem1
                   ENDIF
                   q1(i,k) = q1(i,k) + qevap(i)
                   t1(i,k) = t1(i,k) - elocp * qevap(i)
                   deltv(i) = - elocp*qevap(i)/dt2
                   delq(i) =  + qevap(i)/dt2
                   delqev(i) = delqev(i) + 0.001_r8*dp*qevap(i)/g
                ENDIF
                dellaq(i,k) = dellaq(i,k) + delq(i) / xmb(i)
                delqbar(i) = delqbar(i) + delq(i)*dp/g
                deltbar(i) = deltbar(i) + deltv(i)*dp/g
             ENDIF
          ENDIF
       ENDDO
    ENDDO
    !j
    !     do i = 1, im
    !     if(me.eq.31.and.cnvflg(i)) then
    !     if(cnvflg(i)) then
    !       print *, ' shallow delhbar, delqbar, deltbar = ',
    !    &             delhbar(i),con_hvap*delqbar(i),cp*deltbar(i)
    !       print *, ' shallow delubar, delvbar = ',delubar(i),delvbar(i)
    !       print *, ' precip =', con_hvap*rn(i)*1000./dt2
    !       print*,'pdif= ',pfld(i,kbcon(i))-pfld(i,ktcon(i))
    !     endif
    !     enddo
    !j
    DO i = 1, im
       IF(cnvflg(i)) THEN
          IF(rn(i).LT.0.0_r8.OR..NOT.flg(i)) rn(i) = 0.0_r8
          ktop(i) = ktcon(i)
          kbot(i) = kbcon(i)
          noshal(i) = 0
       ENDIF
    ENDDO
    !
    !  cloud water
    !
    IF (ncloud.GT.0) THEN
       !
       DO k = 1, km1
          DO i = 1, im
             IF (cnvflg(i)) THEN
                IF (k.GT.kb(i).AND.k.LE.ktcon(i)) THEN
                   tem  = dellal(i,k) * xmb(i) * dt2
                   tem1 = MAX(0.0_r8, MIN(1.0_r8, (tcr-t1(i,k))*tcrf))
                   IF (ql(i,k,2) .GT. -999.0_r8) THEN
                      ql(i,k,1) = ql(i,k,1) + tem * tem1              ! ice
                      ql(i,k,2) = ql(i,k,2) + tem *(1.0_r8-tem1)       ! water
                   ELSE
                      ql(i,k,1) = ql(i,k,1) + tem
                   ENDIF
                ENDIF
             ENDIF
          ENDDO
       ENDDO
       !
    ENDIF
    !
    ! hchuang code change
    !
    DO k = 1, km
       DO i = 1, im
          IF(cnvflg(i)) THEN
             IF(k.GE.kb(i) .AND. k.LT.ktop(i)) THEN
                ud_mf(i,k) = eta(i,k) * xmb(i) * dt2
             ENDIF
          ENDIF
       ENDDO
    ENDDO
    DO i = 1, im
       IF(cnvflg(i)) THEN
          k = ktop(i)-1
          dt_mf(i,k) = ud_mf(i,k)
       ENDIF
    ENDDO
    !!
    RETURN
  END  SUBROUTINE shalcnv
  
  
  !-------------------------------------------------------------------------------

      SUBROUTINE sig2press(nCols    ,& !
                           kMax     ,&!
                           pgr      ,&!
                           sl       ,&!
                           si       ,&!
                           prsi     ,&!
                           prsl     ,&!
                           prsik    ,&!
                           prslk      )
 
      IMPLICIT NONE
 
      INTEGER      , INTENT(IN   ) ::  nCols
      INTEGER      , INTENT(IN   ) ::  kMax
      REAL(kind=r8), INTENT(IN   ) ::  pgr(nCols)    !cb
      REAL(kind=r8), INTENT(IN   ) ::  sl(kMax)
      REAL(kind=r8), INTENT(IN   ) ::  si(kMax+1)
      REAL(kind=r8), INTENT(OUT  ) ::  prsi(nCols,kMax+1)
      REAL(kind=r8), INTENT(OUT  ) ::  prsl(nCols,kMax)
      REAL(kind=r8), INTENT(OUT  ) ::  prsik(nCols,kMax+1)
      REAL(kind=r8), INTENT(OUT  ) ::  prslk(nCols,kMax)
      
      REAL(kind=r8) :: pgrk(nCols)
      REAL(kind=r8) :: tem, rkapi, rkapp1
      REAL(kind=r8) :: slk(kMax) 
      REAL(kind=r8) :: sik(kMax+1)
      REAL(kind=r8) :: sikp1(kMax+1)
      INTEGER       :: i,k

      REAL(kind=r8), PARAMETER :: PT01=0.01
      REAL(kind=r8), PARAMETER:: con_rd     =2.8705e+2      ! gas constant air    (J/kg/K)
      REAL(kind=r8), PARAMETER:: con_cp     =1.0046e+3      ! spec heat air @p    (J/kg/K)
      REAL(kind=r8), PARAMETER:: con_rocp   =con_rd/con_cp
      REAL(kind=r8), PARAMETER:: rkap = con_rocp
      REAL(kind=r8), PARAMETER:: rk = con_rocp
!      REAL(kind=r8), PARAMETER :: cb2mb   = 10.0

      RKAPI  = 1.0 / RKAP
      RKAPP1 = 1.0 + RKAP
      DO k=1,kMax+1
        sik(k)   = si(k) ** rkap
        sikp1(k) = si(k) ** rkapp1
      END DO
      DO k=1,kMax
        tem      = rkapp1 * (si(k) - si(k+1))
        slk(k)   = (sikp1(k)-sikp1(k+1))/tem
        !sl(k)    = slk(k) ** rkapi
      END DO

      DO i=1,nCols
         prsi(i,kMax+1)  = si(kMax+1)*pgr(i)      ! prsi are now pressures
         pgrk(i)         = (pgr(i)*pt01) ** rk
         prsik(i,kMax+1) = sik(kMax+1) * pgrk(i)
      END DO
      DO k=1,kMax
        DO i=1,nCols
          prsi(i,k)  = si(k)*pgr(i)               ! prsi are now pressures
          prsl(i,k)  = sl(k)*pgr(i)
          prsik(i,k) = sik(k) * pgrk(i)
          prslk(i,k) = slk(k) * pgrk(i)
        END DO
      END DO
      RETURN
      END SUBROUTINE sig2press

  !-----------------------------------------------------------------------------------------  
  
    !-------------------------------------------------------------------------------
  SUBROUTINE gfuncphys()
    !$$$     Subprogram Documentation Block
    !
    ! Subprogram: gfuncphys    Compute all physics function tables
    !   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
    !
    ! Abstract: Compute all physics function tables.  Lookup tables are
    !   set up for computing saturation vapor pressure, dewpoint temperature,
    !   equivalent potential temperature, moist adiabatic temperature and humidity,
    !   pressure to the kappa, and lifting condensation level temperature.
    !
    ! Program History Log:
    ! 1999-03-01  Iredell             f90 module
    !
    ! Usage:  call gfuncphys
    !
    ! Subprograms called:
    !   gpvs        compute saturation vapor pressure table
    !
    ! Attributes:
    !   Language: Fortran 90.
    !
    !$$$
    IMPLICIT NONE
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    CALL gpvs()
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  END SUBROUTINE gfuncphys
  !-------------------------------------------------------------------------------
  !-------------------------------------------------------------------------------
  SUBROUTINE gpvs()
    !$$$     Subprogram Documentation Block
    !
    ! Subprogram: gpvs         Compute saturation vapor pressure table
    !   Author: N Phillips            W/NMC2X2   Date: 30 dec 82
    !
    ! Abstract: Computes saturation vapor pressure table as a function of
    !   temperature for the table lookup function fpvs.
    !   Exact saturation vapor pressures are calculated in subprogram fpvsx.
    !   The current implementation computes a table with a length
    !   of 7501 for temperatures ranging from 180. to 330. Kelvin.
    !
    ! Program History Log:
    !   91-05-07  Iredell
    !   94-12-30  Iredell             expand table
    ! 1999-03-01  Iredell             f90 module
    ! 2001-02-26  Iredell             ice phase
    !
    ! Usage:  call gpvs
    !
    ! Subprograms called:
    !   (fpvsx)    inlinable function to compute saturation vapor pressure
    !
    ! Attributes:
    !   Language: Fortran 90.
    !
    !$$$
    IMPLICIT NONE
    INTEGER jx
    REAL(r8) :: xmin,xmax,xinc,x,t
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xmin=180.0_r8
    xmax=330.0_r8
    xinc=(xmax-xmin)/(nxpvs-1)
    !   c1xpvs=1.-xmin/xinc
    c2xpvs=1.0_r8/xinc
    c1xpvs=1.0_r8-xmin*c2xpvs
    DO jx=1,nxpvs
       x=xmin+(jx-1)*xinc
       t=x
       tbpvs(jx)=fpvsx(t)
    ENDDO
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  END SUBROUTINE gpvs 
  !-------------------------------------------------------------------------------
  !-------------------------------------------------------------------------------
  ELEMENTAL FUNCTION fpvs(t)
    !$$$     Subprogram Documentation Block
    !
    ! Subprogram: fpvs         Compute saturation vapor pressure
    !   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
    !
    ! Abstract: Compute saturation vapor pressure from the temperature.
    !   A linear interpolation is done between values in a lookup table
    !   computed in gpvs. See documentation for fpvsx for details.
    !   Input values outside table range are reset to table extrema.
    !   The interpolation accuracy is almost 6 decimal places.
    !   On the Cray, fpvs is about 4 times faster than exact calculation.
    !   This function should be expanded inline in the calling routine.
    !
    ! Program History Log:
    !   91-05-07  Iredell             made into inlinable function
    !   94-12-30  Iredell             expand table
    ! 1999-03-01  Iredell             f90 module
    ! 2001-02-26  Iredell             ice phase
    !
    ! Usage:   pvs=fpvs(t)
    !
    !   Input argument list:
    !     t          Real(r8) temperature in Kelvin
    !
    !   Output argument list:
    !     fpvs       Real(r8) saturation vapor pressure in Pascals
    !
    ! Attributes:
    !   Language: Fortran 90.
    !
    !$$$
    IMPLICIT NONE
    REAL(r8) :: fpvs
    REAL(r8),INTENT(in):: t
    INTEGER jx
    REAL(r8) xj
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xj=MIN(MAX(c1xpvs+c2xpvs*t,1.0_r8),REAL(nxpvs,r8))
    jx=INT(MIN(xj,nxpvs-1.0_r8))
    fpvs=tbpvs(jx)+(xj-jx)*(tbpvs(jx+1)-tbpvs(jx))
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  END FUNCTION fpvs
  !-------------------------------------------------------------------------------
  !-------------------------------------------------------------------------------
  ELEMENTAL FUNCTION fpvsx(t)
    !$$$     Subprogram Documentation Block
    !
    ! Subprogram: fpvsx        Compute saturation vapor pressure
    !   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
    !
    ! Abstract: Exactly compute saturation vapor pressure from temperature.
    !   The saturation vapor pressure over either liquid and ice is computed
    !   over liquid for temperatures above the triple point,
    !   over ice for temperatures 20 degress below the triple point,
    !   and a linear combination of the two for temperatures in between.
    !   The water model assumes a perfect gas, constant specific heats
    !   for gas, liquid and ice, and neglects the volume of the condensate.
    !   The model does account for the variation of the latent heat
    !   of condensation and sublimation with temperature.
    !   The Clausius-Clapeyron equation is integrated from the triple point
    !   to get the formula
    !       pvsl=con_psat*(tr**xa)*exp(xb*(1.-tr))
    !   where tr is ttp/t and other values are physical constants.
    !   The reference for this computation is Emanuel(1994), pages 116-117.
    !   This function should be expanded inline in the calling routine.
    !
    ! Program History Log:
    !   91-05-07  Iredell             made into inlinable function
    !   94-12-30  Iredell             exact computation
    ! 1999-03-01  Iredell             f90 module
    ! 2001-02-26  Iredell             ice phase
    !
    ! Usage:   pvs=fpvsx(t)
    !
    !   Input argument list:
    !     t          Real(r8) temperature in Kelvin
    !
    !   Output argument list:
    !     fpvsx      Real(r8) saturation vapor pressure in Pascals
    !
    ! Attributes:
    !   Language: Fortran 90.
    !
    !$$$
    IMPLICIT NONE
    REAL(r8) :: fpvsx
    REAL(r8),INTENT(in) :: t
    REAL(r8),PARAMETER :: tliq=con_ttp
    REAL(r8),PARAMETER :: tice=con_ttp-20.0
    REAL(r8),PARAMETER :: dldtl=con_cvap-con_cliq
    REAL(r8),PARAMETER :: heatl=con_hvap
    REAL(r8),PARAMETER :: xponal=-dldtl/con_rv
    REAL(r8),PARAMETER :: xponbl=-dldtl/con_rv+heatl/(con_rv*con_ttp)
    REAL(r8),PARAMETER :: dldti=con_cvap-con_csol
    REAL(r8),PARAMETER :: heati=con_hvap+con_hfus
    REAL(r8),PARAMETER :: xponai=-dldti/con_rv
    REAL(r8),PARAMETER :: xponbi=-dldti/con_rv+heati/(con_rv*con_ttp)
    REAL(r8) tr,w,pvl,pvi
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    tr=con_ttp/t
    IF(t.GE.tliq) THEN
       fpvsx=con_psat*(tr**xponal)*EXP(xponbl*(1.0_r8-tr))
    ELSEIF(t.LT.tice) THEN
       fpvsx=con_psat*(tr**xponai)*EXP(xponbi*(1.0_r8-tr))
    ELSE
       w=(t-tice)/(tliq-tice)
       pvl=con_psat*(tr**xponal)*EXP(xponbl*(1.0_r8-tr))
       pvi=con_psat*(tr**xponai)*EXP(xponbi*(1.0_r8-tr))
       fpvsx=w*pvl+(1.0_r8-w)*pvi
    ENDIF
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  END FUNCTION  fpvsx
  
  
  

END MODULE Shall_MasFlux


!PROGRAM Main
!  USE Shall_MasFlux
!END PROGRAM MAin
