!WRF:MODEL_LAYER:PHYSICS
!
MODULE module_bl_mrf

CONTAINS

!-------------------------------------------------------------------          
   SUBROUTINE MRF(U3D,V3D,TH3D,T3D,QV3D,QC3D,P3D,PI3D,             &
                  RUBLTEN,RVBLTEN,RTHBLTEN,                        &
                  RQVBLTEN,RQCBLTEN,                               &
                  CP,G,ROVCP,R,ROVG,                               &
                  dz8w,z,XLV,RV,PSFC,                              &
                  ZNT,UST,ZOL,HOL,PBL,PSIM,PSIH,                   &
                  XLAND,HFX,QFX,TSK,GZ1OZ0,WSPD,BR,                &
                  DT,DTMIN,KPBL2D,                                 &
                  SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN,EOMEG,STBOLT,&
                  flag_qi,                                         &
                  ids,ide, jds,jde, kds,kde,                       &
                  ims,ime, jms,jme, kms,kme,                       &
                  its,ite, jts,jte, kts,kte,                       &
               ! Optional
                  QI3D,RQIBLTEN,                                   & 
                  regime                                           )
!-------------------------------------------------------------------
      IMPLICIT NONE
!-------------------------------------------------------------------
!-- U3D         3D u-velocity interpolated to theta points (m/s)
!-- V3D         3D v-velocity interpolated to theta points (m/s)
!-- TH3D        3D potential temperature (K)
!-- T3D         temperature (K)
!-- QV3D        3D water vapor mixing ratio (Kg/Kg)
!-- QC3D        3D cloud mixing ratio (Kg/Kg)
!-- QI3D        3D ice mixing ratio (Kg/Kg)
!-- P3D         3D pressure (Pa)
!-- PI3D        3D exner function (dimensionless)
!-- rr3D        3D dry air density (kg/m^3)
!-- RUBLTEN     U tendency due to
!               PBL parameterization (m/s^2)
!-- RVBLTEN     V tendency due to
!               PBL parameterization (m/s^2)
!-- RTHBLTEN    Theta tendency due to
!               PBL parameterization (K/s)
!-- RQVBLTEN    Qv tendency due to
!               PBL parameterization (kg/kg/s)
!-- RQCBLTEN    Qc tendency due to
!               PBL parameterization (kg/kg/s)
!-- RQIBLTEN    Qi tendency due to
!               PBL parameterization (kg/kg/s)
!-- CP          heat capacity at constant pressure for dry air (J/kg/K)
!-- G           acceleration due to gravity (m/s^2)
!-- ROVCP       R/CP
!-- R           gas constant for dry air (J/kg/K)
!-- ROVG        R/G
!-- dz8w        dz between full levels (m)
!-- z           height above sea level (m)
!-- XLV         latent heat of vaporization (J/kg)
!-- RV          gas constant for water vapor (J/kg/K)
!-- PSFC        pressure at the surface (Pa)
!-- ZNT         roughness length (m)
!-- UST         u* in similarity theory (m/s)
!-- ZOL         z/L height over Monin-Obukhov length
!-- HOL         PBL height over Monin-Obukhov length
!-- PBL         PBL height (m)
!-- REGIME      flag indicating PBL regime (stable, unstable, etc.)
!-- PSIM        similarity stability function for momentum
!-- PSIH        similarity stability function for heat
!-- XLAND       land mask (1 for land, 2 for water)
!-- HFX         upward heat flux at the surface (W/m^2)
!-- QFX         upward moisture flux at the surface (kg/m^2/s)
!-- TSK         surface temperature (K)
!-- GZ1OZ0      log(z/z0) where z0 is roughness length
!-- WSPD        wind speed at lowest model level (m/s)
!-- BR          bulk Richardson number in surface layer
!-- DT          time step (s)
!-- DTMIN       time step (minute)
!-- rvovrd      R_v divided by R_d (dimensionless)
!-- SVP1        constant for saturation vapor pressure (kPa)
!-- SVP2        constant for saturation vapor pressure (dimensionless)
!-- SVP3        constant for saturation vapor pressure (K)
!-- SVPT0       constant for saturation vapor pressure (K)
!-- EP1         constant for virtual temperature (R_v/R_d - 1) (dimensionless)
!-- EP2         constant for specific humidity calculation
!-- KARMAN      Von Karman constant
!-- EOMEG       angular velocity of earth's rotation (rad/s)
!-- STBOLT      Stefan-Boltzmann constant (W/m^2/K^4)
!-- ids         start index for i in domain
!-- ide         end index for i in domain
!-- jds         start index for j in domain
!-- jde         end index for j in domain
!-- kds         start index for k in domain
!-- kde         end index for k in domain
!-- ims         start index for i in memory
!-- ime         end index for i in memory
!-- jms         start index for j in memory
!-- jme         end index for j in memory
!-- kms         start index for k in memory
!-- kme         end index for k in memory
!-- its         start index for i in tile
!-- ite         end index for i in tile
!-- jts         start index for j in tile
!-- jte         end index for j in tile
!-- kts         start index for k in tile
!-- kte         end index for k in tile
!-------------------------------------------------------------------

      INTEGER,  INTENT(IN   )   ::      ids,ide, jds,jde, kds,kde, &
                                        ims,ime, jms,jme, kms,kme, &
                                        its,ite, jts,jte, kts,kte

!
      REAL,     INTENT(IN   )   ::      DT,DTMIN,CP,G,ROVCP,       &
                                        ROVG,R,XLV,RV             

      REAL,     INTENT(IN )     ::     SVP1,SVP2,SVP3,SVPT0 
      REAL,     INTENT(IN )     ::     EP1,EP2,KARMAN,EOMEG,STBOLT

      REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
                INTENT(IN   )   ::                           QV3D, &
                                                             QC3D, &
                                                              P3D, &
                                                             PI3D, &
                                                             TH3D, &
                                                              T3D, &
                                                             dz8w, &
                                                                z   
!
      REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
                INTENT(INOUT)   ::                        RUBLTEN, &
                                                          RVBLTEN, &
                                                         RTHBLTEN, &
                                                         RQVBLTEN, &
                                                         RQCBLTEN

      REAL,     DIMENSION( ims:ime, jms:jme )                    , &
                INTENT(IN   )   ::                          XLAND, &
                                                              HFX, &
                                                              QFX
 
      REAL,     DIMENSION( ims:ime, jms:jme )                    , &
                INTENT(INOUT)   ::                            HOL, &
                                                              UST, &
                                                              PBL, &
                                                              ZNT

      LOGICAL,  INTENT(IN)      ::                        FLAG_QI
!
!m The following 5 variables are changed to memory size from tile size--
!
     REAL,     DIMENSION( ims:ime, jms:jme ), INTENT(IN   )  ::    &
                                                             PSIM, &
                                                             PSIH

     REAL,     DIMENSION( ims:ime, jms:jme ), INTENT(INOUT)   ::   &
                                                             WSPD

     REAL,     DIMENSION( ims:ime, jms:jme ), INTENT(IN   )   ::   &
                                                           GZ1OZ0, &
                                                               BR

      REAL,     DIMENSION( ims:ime, jms:jme )                    , &
                INTENT(IN   )               ::               PSFC

      REAL,     DIMENSION( ims:ime, jms:jme )                    , &
                INTENT(IN   )   ::                            TSK

      REAL,     DIMENSION( ims:ime, jms:jme )                    , &
                INTENT(INOUT)   ::                            ZOL
                                                                      
      INTEGER,  DIMENSION( ims:ime, jms:jme )                    , &
                INTENT(OUT  )   ::                         KPBL2D 

      REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
                INTENT(IN   )   ::                            U3D, &
                                                              V3D
!
! Optional
!
      REAL,     DIMENSION( ims:ime, jms:jme )                    , &
                OPTIONAL                                         , &
                INTENT(INOUT)   ::                         REGIME

      REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
                INTENT(INOUT)   ::                       RQIBLTEN

      REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
                OPTIONAL                                         , &
                INTENT(IN   )   ::                           QI3D

! LOCAL VARS
   REAL,       DIMENSION( its:ite, kts:kte )          ::   dz8w2d, & 
                                                              z2d
                                                           

   INTEGER ::  I,J,K,NK

!
   DO J=jts,jte
      DO k=kts,kte
      NK=kme-k
      DO i=its,ite
         dz8w2d(I,K) = dz8w(i,NK,j)
         z2d(I,K) = z(i,NK,j)
      ENDDO
      ENDDO


      CALL MRF2D(J,U3D(ims,kms,j),V3D(ims,kms,j),T3D(ims,kms,j),    &
               QV3D(ims,kms,j),QC3D(ims,kms,j),                     &
               P3D(ims,kms,j),RUBLTEN(ims,kms,j),RVBLTEN(ims,kms,j),&
               RTHBLTEN(ims,kms,j),RQVBLTEN(ims,kms,j),             &
               RQCBLTEN(ims,kms,j),                                 &
               CP,G,ROVCP,R,ROVG,                                   &
               dz8w2d,z2d,XLV,Rv,                                   &
               PSFC(ims,j),ZNT(ims,j),                              &
               UST(ims,j),ZOL(ims,j),                               &
               HOL(ims,j),PBL(ims,j),PSIM(ims,j),                   &
               PSIH(ims,j),XLAND(ims,j),HFX(ims,j),QFX(ims,j),      &
               TSK(ims,j),GZ1OZ0(ims,j),WSPD(ims,j),BR(ims,j),      &
               DT,DTMIN,KPBL2D(ims,j),                              &
               SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN,EOMEG,STBOLT,    &
               flag_qi,                                             &
               ids,ide, jds,jde, kds,kde,                           &
               ims,ime, jms,jme, kms,kme,                           &
               its,ite, jts,jte, kts,kte,                           &
            !optional
               QI2DTEN=RQIBLTEN(ims,kms,j),                         &
               REGIME=REGIME(ims,j),QI2D=QI3D(ims,kms,j)            )


      DO k=kts,kte
      DO i=its,ite
         RTHBLTEN(I,K,J)=RTHBLTEN(I,K,J)/PI3D(I,K,J)
      ENDDO
      ENDDO

    ENDDO

   END SUBROUTINE MRF

!-------------------------------------------------------------------          
   SUBROUTINE MRF2D(J,U2D,V2D,T2D,QV2D,QC2D,     P2D,              &
                  U2DTEN,V2DTEN,T2DTEN,                            &
                  QV2DTEN,QC2DTEN,                                 & 
                  CP,G,ROVCP,R,ROVG,                               &
                  dz8w2d,z2d,XLV,RV,PSFCPA,                        &
                  ZNT,UST,ZOL,HOL,PBL,PSIM,PSIH,                   &
                  XLAND,HFX,QFX,TSK,GZ1OZ0,WSPD,BR,                &
                  DT,DTMIN,KPBL1D,                                 &
                  SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN,EOMEG,STBOLT,&
                  flag_qi,                                         &
                  ids,ide, jds,jde, kds,kde,                       &
                  ims,ime, jms,jme, kms,kme,                       &
                  its,ite, jts,jte, kts,kte,                       &
               ! optional
                  regime, qi2d, QI2DTEN                            )
!-------------------------------------------------------------------
      IMPLICIT NONE
!-------------------------------------------------------------------
!     BASED ON THE "COUNTERGRADIENT" TRANSPORT TERM OF TROEN                     
!     AND MAHRT (1986) FOR THE UNSTABLE PBL.                                     
!     THIS ROUTINE USES AN IMPLICIT APPROACH FOR VERTICAL FLUX                   
!     DIVERGENCE AND DOES NOT REQUIRE "MITER" TIMESTEPS.                         
!     IT INCLUDES VERTICAL DIFFUSION IN THE STABLE ATMOSPHERE                    
!     AND MOIST VERTICAL DIFFUSION IN CLOUDS.                                    
!     SURFACE FLUXES CALCULATED AS IN HIRPBL.                                    
!     5-LAYER SOIL MODEL OPTION REQUIRED IN SLAB DUE TO LONG TIMESTEP            
!                                                                                
!     CODED BY SONG-YOU HONG (NCEP), IMPLEMENTED BY JIMY DUDHIA (NCAR)           
!     FALL 1996                                                                  
!                                                                                
!     REFERENCES:                                                                
!                                                                                
!        HONG AND PAN (1996), MON. WEA. REV.                                     
!        TROEN AND MAHRT (1986), BOUNDARY LAYER MET.                             
!                                                                                
!     CHANGES:                                                                   
!        INCREASE RLAM FROM 30 TO 150, AND CHANGE FREE ATMOSPHERE                
!        STABILITY FUNCTION TO INCREASE VERTICAL DIFFUSION                       
!        (HONG, JUNE 1997)                                                       
!                                                                                
!        PUT LOWER LIMIT ON PSI FOR STABLE CONDITIONS. THIS WILL                 
!        PREVENT FLUXES FROM BECOMING TOO SMALL (DUDHIA, OCTOBER 1997)           
!                                                                                
!        CORRECTION TO REGIME CALCULATION. THIS WILL ALLOW POINTS IN             
!        REGIME 4 MUCH MORE FREQUENTLY GIVING LARGER SURFACE FLUXES              
!        REGIME 3 NO LONGER USES HOL < 1.5 OR THVX LAPSE-RATE CHECK              
!        IN MRF SCHEME. THIS WILL MAKE REGIME 3 MUCH LESS FREQUENT.              
!                                                                                
!        ADD SURFACE PRESSURE, PS(I), ARRAY FOR EFFICIENCY                       
!                                                                                
!        FIX FOR PROBLEM WITH THIN LAYERS AND HIGH ROUGHNESS                     
!                                                                                
!        CHARNOCK CONSTANT NOW COMES FROM NAMELIST (DEFAULT SAME)                
!                                                                                
!-------------------------------------------------------------------          

      REAL      RLAM,PRMIN,PRMAX,XKZMIN,XKZMAX,RIMIN,BRCR,         &
                CFAC,PFAC,SFCFRAC,CKZ,ZFMIN,APHI5,APHI16,GAMCRT,   &
                GAMCRQ,XKA,PRT

      PARAMETER (RLAM=150.,PRMIN=0.5,PRMAX=4.)                       
      PARAMETER (XKZMIN=0.01,XKZMAX=1000.,RIMIN=-100.)                           
      PARAMETER (BRCR=0.5,CFAC=7.8,PFAC=2.0,SFCFRAC=0.1)                         
      PARAMETER (CKZ=0.001,ZFMIN=1.E-8,APHI5=5.,APHI16=16.)                      
      PARAMETER (GAMCRT=3.,GAMCRQ=2.E-3)                                         
      PARAMETER (XKA=2.4E-5)                                                     
      PARAMETER (PRT=1.)                                                         
!
      INTEGER,  INTENT(IN   )   ::      ids,ide, jds,jde, kds,kde, &
                                        ims,ime, jms,jme, kms,kme, &
                                        its,ite, jts,jte, kts,kte, &
                                        J
!
      LOGICAL,  INTENT(IN)      ::                        FLAG_QI
!
      REAL,     INTENT(IN   )   ::      DT,DTMIN,CP,G,ROVCP,       &
                                        ROVG,R,XLV,RV

      REAL,     INTENT(IN )     ::     SVP1,SVP2,SVP3,SVPT0 
      REAL,     INTENT(IN )     ::     EP1,EP2,KARMAN,EOMEG,STBOLT

      REAL,     DIMENSION( ims:ime, kms:kme )                    , &
                INTENT(IN   )   ::                           QV2D, &
                                                             QC2D, &
                                                              P2D, &
                                                              T2D
!
      REAL,     DIMENSION( ims:ime, kms:kme )                    , &
                INTENT(INOUT)   ::                         U2DTEN, &
                                                           V2DTEN, &
                                                           T2DTEN, &
                                                          QV2DTEN, &
                                                          QC2DTEN
                                                                  
      REAL,     DIMENSION( ims:ime )                             , &
                INTENT(INOUT)   ::                            HOL, &
                                                              UST, &
                                                              PBL, &
                                                              ZNT

      REAL,     DIMENSION( ims:ime )                             , &
                INTENT(IN   )   ::                          XLAND, &
                                                              HFX, &
                                                              QFX
!
!m The following 5 are changed to memory size---
!
     REAL,     DIMENSION( ims:ime ), INTENT(IN   )   ::      PSIM, &
                                                             PSIH

     REAL,     DIMENSION( ims:ime ), INTENT(INOUT)   ::      WSPD

     REAL,     DIMENSION( ims:ime ), INTENT(IN   )   ::    GZ1OZ0, &
                                                               BR

      REAL,     DIMENSION( ims:ime )                             , &
                INTENT(IN   )               ::             PSFCPA

      REAL,     DIMENSION( ims:ime )                             , &
                INTENT(IN   )   ::                            TSK

      REAL,     DIMENSION( ims:ime )                             , &
                INTENT(INOUT)   ::                            ZOL

      INTEGER,  DIMENSION( ims:ime )                             , &
                INTENT(OUT  )   ::                         KPBL1D

      REAL,     DIMENSION( ims:ime, kms:kme )                    , &
                INTENT(IN   )   ::                            U2D, &
                                                              V2D
                                                                      
! MODULE-LOCAL VARIABLES (DEFINED IN SUBROUTINE MRF)
!
      REAL,     DIMENSION( its:ite, kts:kte ) ,                    &
                INTENT(IN)      ::                         dz8w2d, &
                                                              z2d
!                                                                                
! 
! Optional
!
      REAL,     DIMENSION( ims:ime )                             , &
                OPTIONAL                                         , &
                INTENT(INOUT)   ::                         REGIME

      REAL,     DIMENSION( ims:ime, kms:kme )                    , &
                OPTIONAL                                         , &
                INTENT(IN   )   ::                           QI2D

      REAL,     DIMENSION( ims:ime, kms:kme )                    , &
                OPTIONAL                                         , &
                INTENT(INOUT)   ::                        QI2DTEN
    
! LOCAL VARS

      REAL,     DIMENSION( its:ite, kts:kte+1 ) ::             ZQ

      REAL,     DIMENSION( its:ite, kts:kte )   ::                 &
                                                         UX,VX,QX, &
                                                     QCX,THX,THVX, & 
                                                          DZQ,DZA, & 
                                                        TTNP,QTNP, &
                                                         QCTNP,ZA, &
                                                          UXS,VXS, &
                                                         THXS,QXS, &
                                                         QCXS,QIX, &
                                                       QITNP,QIXS, &
                                                        UTNP,VTNP
!
      REAL,    DIMENSION( its:ite )             ::     QIXSV,RHOX, &
                                                     WSPD1,GOVRTH, & 
                                                       PBL0,THXSV, &
                                                        UXSV,VXSV, &             
                                                       QXSV,QCXSV, & 
                                                     QGH,TGDSA,PS

      INTEGER                                   ::   ILXM,JLXM,KL, &
                                                   KLM,KLP1,KLPBL
!
      INTEGER, DIMENSION( its:ite )             ::     KPBL,KPBL0
!
      REAL,    DIMENSION( its:ite, kts:kte )    ::      SCR3,SCR4
!
      REAL,    DIMENSION( its:ite )             ::           DUM1, &
                                                           XKZMKL
!
      REAL,    DIMENSION( its:ite )             ::    ZL1,THERMAL, &
                                                     WSCALE,HGAMT, &
                                                       HGAMQ,BRDN, &
                                                        BRUP,PHIM, &
                                                         PHIH,CPM, &
                                                      DUSFC,DVSFC, &
                                                      DTSFC,DQSFC
                
!
      REAL,    DIMENSION( its:ite, kts:kte )    ::      XKZM,XKZH, &
                                                            A1,A2, &  
                                                            AD,AU, &
                                                               TX
!
      REAL,    DIMENSION( its:ite, kts:kte )  ::               AL
!
      LOGICAL, DIMENSION( its:ite )             ::         PBLFLG, &
                                                           SFCFLG, &
                                                           STABLE
!                                                         
      REAL, DIMENSION( its:ite )                ::           THGB

      INTEGER ::  N,I,K,KK,L,NZOL,IMVDIF

      INTEGER ::  JBGN,JEND,IBGN,IEND,NK

      REAL    ::  ZOLN,X,Y,CONT,CONQ,CONW,PL,THCON,TVCON,E1,DTSTEP
      REAL    ::  ZL,TSKV,DTHVDZ,DTHVM,VCONV,RZOL
      REAL    ::  DTTHX,PSIX,DTG,PSIQ,USTM
      REAL    ::  DT4,RDT,SPDK2,FM,FH,HOL1,GAMFAC,VPERT,PRNUM
      REAL    ::  ZFAC,XKZO,SS,RI,QMEAN,TMEAN,ALPH,CHI,ZK,RL2,DK,SRI
      REAL    ::  BRINT,DTODSD,DSIG,RDZ,DSDZT,DSDZQ,DSDZ2,TTEND,QTEND
      REAL    ::  UTEND,VTEND,QCTEND,QITEND,TGC,DTODSU

!----------------------------------------------------------------------

      KLPBL=1
      KL=kte
      ILXM=ite-1
      JLXM=jte-1
      KLM=kte-1
      KLP1=kte+1
!                                                                                
      CONT=1000.*CP/G                                                            
      CONQ=1000.*XLV/G                                                           
      CONW=1000./G                                                               

!-- IMVDIF      imvdif=1 for moist adiabat vertical diffusion

      IMVDIF=1

!     DO i=its,ite
!!PS PSFC cmb
!        PSFC(I)=PSFCPA(I)/1000.
!     ENDDO


!                                                                                
!----CONVERT GROUND TEMPERATURE TO POTENTIAL TEMPERATURE:                        
!                                                                                
      DO 5 I=its,ite                                                       
        TGDSA(I)=TSK(I)                                                        
! PS PSFC cmb
        PS(I)=PSFCPA(I)/1000.
        THGB(I)=TSK(I)*(100./PS(I))**ROVCP   
    5 CONTINUE                                                                   
!                                                                                
!-----DECOUPLE FLUX-FORM VARIABLES TO GIVE U,V,T,THETA,THETA-VIR.,               
!     T-VIR., QV, AND QC AT CROSS POINTS AND AT KTAU-1.                          
!                                                                                
!     *** NOTE ***                                                               
!         THE BOUNDARY WINDS MAY NOT BE ADEQUATELY AFFECTED BY FRICTION,         
!         SO USE ONLY INTERIOR VALUES OF UX AND VX TO CALCULATE                  
!         TENDENCIES.                                                            
!                                                                                
      DO 24 K=kts,kte                                                      
        NK=kme-K
        DO 24 I=its,ite
          UX(I,K)=U2D(I,NK)
          VX(I,K)=V2D(I,NK)
   24 CONTINUE                                                                 
!                                                                                
!.....SCR3(I,K) STORE TEMPERATURE,                                               
!     SCR4(I,K) STORE VIRTUAL TEMPERATURE.                                       
!                                                                                
      DO 30 K=kts,kte
        NK=kme-K
        DO 30 I=its,ite
! PL cmb
          PL=P2D(I,NK)/1000.
          SCR3(I,K)=T2D(I,NK)
          THCON=(100./PL)**ROVCP                                                 
          THX(I,K)=SCR3(I,K)*THCON                                               
          TX(I,K)=SCR3(I,K)                                                      
          SCR4(I,K)=SCR3(I,K)                                                    
          THVX(I,K)=THX(I,K)                                                     
          QX(I,K)=0.                                                             
   30 CONTINUE                                                                 
!                                                                                
      DO I=its,ite
         QGH(i)=0.                                                                
         CPM(i)=CP                                                                
      ENDDO
!                                                                                
!     IF(IDRY.EQ.1)GOTO 80                                                   
      DO 50 K=kts,kte
        NK=kme-K
        DO 50 I=its,ite
          QX(I,K)=QV2D(I,NK)
          TVCON=(1.+EP1*QX(I,K))                                      
          THVX(I,K)=THX(I,K)*TVCON                                               
          SCR4(I,K)=SCR3(I,K)*TVCON                                              
   50 CONTINUE                                                                 
!                                                                                
      DO 60 I=its,ite
        E1=SVP1*EXP(SVP2*(TGDSA(I)-SVPT0)/(TGDSA(I)-SVP3))                       
        QGH(I)=EP2*E1/(PS(I)-E1)                                                 
        CPM(I)=CP*(1.+0.8*QX(I,KL))                                   
   60 CONTINUE                                                                   
!                                                                                
!     IF(IMOIST.EQ.1)GOTO 80
      DO 70 K=kts,kte
        NK=kme-K
        DO 70 I=its,ite
          QCX(I,K)=QC2D(I,NK)
   70 CONTINUE

      IF (flag_QI .AND. PRESENT( QI2D ) ) THEN
         DO K=kts,kte
            NK=kme-K
            DO I=its,ite
               QIX(I,K)=QI2D(I,NK)
            ENDDO
         ENDDO
      ELSE
         DO K=kts,kte
            NK=kme-K
            DO I=its,ite
               QIX(I,K)=0.
            ENDDO
         ENDDO
      ENDIF

   80 CONTINUE

!                                                                                
!-----COMPUTE THE HEIGHT OF FULL- AND HALF-SIGMA LEVELS ABOVE GROUND             
!     LEVEL, AND THE LAYER THICKNESSES.                                          
!                                                                                
      DO 90 I=its,ite
        ZQ(I,KLP1)=0.                                                            
        RHOX(I)=PS(I)*1000./(R*SCR4(I,KL))                                       
   90 CONTINUE                                                                   
!                                                                                
      DO 110 KK=kts,kte
        K=kme-KK
        DO 100 I=its,ite                                                   
          DUM1(I)=ZQ(I,K+1)                                                      
  100   CONTINUE                                                                 
!                                                                                
        DO 110 I=its,ite                                                   
           ZQ(I,K)=dz8w2d(I,K)+DUM1(I)
  110   CONTINUE                                                                 
!                                                                                
      DO 120 K=kts,kte
        DO 120 I=its,ite
          ZA(I,K)=0.5*(ZQ(I,K)+ZQ(I,K+1))                                        
          DZQ(I,K)=ZQ(I,K)-ZQ(I,K+1)                                             
  120 CONTINUE                                                                 
!                                                                                
      DO 130 K=kts,kte-1
        DO 130 I=its,ite
          DZA(I,K)=ZA(I,K)-ZA(I,K+1)                                             
  130 CONTINUE                                                                 
               
      DTSTEP=DT                                                                  
!                                                                                
      DO 160 I=its,ite
        GOVRTH(I)=G/THX(I,KL)                                                    
  160 CONTINUE                                                                   
!                                                                                
!-----INITIALIZE VERTICAL TENDENCIES AND                                         
!                                                                                
      DO I=its,ite
      DO K=kts,kte
         UTNP(i,k)=0.                                                           
         VTNP(i,k)=0.                                                           
         TTNP(i,k)=0.                                                           
      ENDDO
      ENDDO
!                                                                                
!     IF(IDRY.EQ.1)GOTO 250                                                  
      DO 230 K=kts,kte
        DO 230 I=its,ite
          QTNP(I,K)=0.                                                           
  230 CONTINUE                                                                 
!                                                                                
!     IF(IMOIST.EQ.1)GOTO 250                                                
      DO 240 K=kts,kte
        DO 240 I=its,ite
          QCTNP(I,K)=0.                                                          
          QITNP(I,K)=0.                                                          
  240 CONTINUE                                                                 
                                                                                 
  250 CONTINUE                                                                   
!                                                                                
!-----CALCULATE BULK RICHARDSON NO. OF SURFACE LAYER, ACCORDING TO               
!     AKB(1976), EQ(12).                                                         
                                                                                 
!     DO 260 I=its,ite
!       GZ1OZ0(I)=ALOG(ZA(I,KL)/ZNT(I))                                        
!       IF((XLAND(I)-1.5).GE.0)THEN                                            
!         ZL=ZNT(I)                                                            
!       ELSE                                                                     
!         ZL=0.01                                                                
!       ENDIF                                                                    
!       WSPD(I)=SQRT(UX(I,KL)*UX(I,KL)+VX(I,KL)*VX(I,KL))                        
!       TSKV=THGB(I)*(1.+EP1*QGH(I)*MAVAIL(I))                     
!       DTHVDZ=(THVX(I,KL)-TSKV)                                                 
!       IF(-DTHVDZ.GE.0)THEN                                                     
!         DTHVM=-DTHVDZ                                                          
!       ELSE                                                                     
!         DTHVM=0.                                                               
!       ENDIF                                                                    
!       VCONV=VCONVC*SQRT(DTHVM)                                                 
!       WSPD(I)=SQRT(WSPD(I)*WSPD(I)+VCONV*VCONV)                                
!       WSPD(I)=AMAX1(WSPD(I),1.)                                                
!       BR(I)=GOVRTH(I)*ZA(I,KL)*DTHVDZ/(WSPD(I)*WSPD(I))                        
! 260 CONTINUE                                                                   
                                                                                 
!!-----DIAGNOSE BASIC PARAMETERS FOR THE APPROPRIATED STABILITY CLASS:            
!!                                                                                
!!     THE STABILITY CLASSES ARE DETERMINED BY BR (BULK RICHARDSON NO.)           
!!     AND HOL (HEIGHT OF PBL/MONIN-OBUKHOV LENGTH).                              
!!                                                                                
!!     CRITERIA FOR THE CLASSES ARE AS FOLLOWS:                                   
!!                                                                                
!!        1. BR .GE. 0.2;                                                         
!!               REPRESENTS NIGHTTIME STABLE CONDITIONS (REGIME=1),               
!!                                                                                
!!        2. BR .LT. 0.2 .AND. BR .GT. 0.0;                                       
!!               REPRESENTS DAMPED MECHANICAL TURBULENT CONDITIONS                
!!               (REGIME=2),                                                      
!!                                                                                
!!        3. BR .EQ. 0.0                                                          
!!               REPRESENTS FORCED CONVECTION CONDITIONS (REGIME=3),              
!!                                                                                
!!        4. BR .LT. 0.0                                                          
!!               REPRESENTS FREE CONVECTION CONDITIONS (REGIME=4).                
!!                                                                                
!!-----                                                                           
!
!      DO 320 I=its,ite
!!----                                                                            
!!--     REMOVE REGIME 3 DEPENDENCE ON PBL HEIGHT                                 
!!--          IF(BR(I).LT.0..AND.HOL(I).GT.1.5)GOTO 310                         
!
!        IF(BR(I).LT.0.)GOTO 310                                                  
!!                                                                                
!!-----CLASS 1; STABLE (NIGHTTIME) CONDITIONS:                                    
!!                                                                                
!        IF(BR(I).LT.0.2)GOTO 270                                                 
!        REGIME(I)=1.                                                           
!        PSIM(I)=-10.*GZ1OZ0(I)                                                   
!!    LOWER LIMIT ON PSI IN STABLE CONDITIONS                                     
!        PSIM(I)=AMAX1(PSIM(I),-10.)                                              
!        PSIH(I)=PSIM(I)                                                          
!        HOL(I)=0.0
!        PBL(I)=0.0                                                             
!        GOTO 320                                                                 
!!                                                                                
!!-----CLASS 2; DAMPED MECHANICAL TURBULENCE:                                     
!!                                                                                
!  270   IF(BR(I).EQ.0.0)GOTO 280                                                 
!        REGIME(I)=2.                                                           
!        PSIM(I)=-5.0*BR(I)*GZ1OZ0(I)/(1.1-5.0*BR(I))                             
!!    LOWER LIMIT ON PSI IN STABLE CONDITIONS                                     
!        PSIM(I)=AMAX1(PSIM(I),-10.)                                              
!!.....AKB(1976), EQ(16).                                                         
!        PSIH(I)=PSIM(I)                                                          
!        HOL(I)=0.0
!        PBL(I)=0.0                                                             
!        GOTO 320                                                                 
!!                                                                                
!!-----CLASS 3; FORCED CONVECTION:                                                
!!                                                                                
!  280   REGIME(I)=3.                                                           
!        PSIM(I)=0.0                                                              
!        PSIH(I)=PSIM(I)                                                          
!                                                                                 
!! special use kte instead of kme
!
!        DO 290 KK=kts,kte-1
!          K=kte-KK
!          IF(THVX(I,K).GT.THVX(I,KL))GOTO 300                                    
!  290   CONTINUE                                                                 
!        STOP 290                                                                 
!  300   PBL(I)=ZQ(I,K+1)                                                       
!        IF(UST(I).LT.0.01)THEN                                                 
!          ZOL(I)=BR(I)*GZ1OZ0(I)                                               
!        ELSE                                                                     
!          ZOL(I)=KARMAN*GOVRTH(I)*ZA(I,KL)*MOL(I,J)/(UST(I)*UST(I)) 
!        ENDIF                                                                    
!        HOL(I)=-ZOL(I)*PBL(I)/ZA(I,KL)
!        GOTO 320                                                                 
!                                                                                 
!!-----CLASS 4; FREE CONVECTION:                                                  
!                                                                                 
!! 310     IF(THVX(I,KLM).GT.THVX(I,KL))GOTO 280                                  
!
!  310   CONTINUE                                                                 
!        REGIME(I)=4.                                                           
!        IF(UST(I).LT.0.01)THEN                                                 
!          ZOL(I)=BR(I)*GZ1OZ0(I)                                               
!        ELSE                                                                     
!          ZOL(I)=KARMAN*GOVRTH(I)*ZA(I,KL)*MOL(I,J)/(UST(I)*UST(I))
!        ENDIF                                                                    
!        ZOL(I)=AMIN1(ZOL(I),0.)                                              
!        ZOL(I)=AMAX1(ZOL(I),-9.9999)                                         
!        NZOL=INT(-ZOL(I)*100.)                                                 
!        RZOL=-ZOL(I)*100.-NZOL                                                 
!        PSIM(I)=PSIMTB(NZOL)+RZOL*(PSIMTB(NZOL+1)-PSIMTB(NZOL))                  
!        PSIH(I)=PSIHTB(NZOL)+RZOL*(PSIHTB(NZOL+1)-PSIHTB(NZOL))                  
!!---LIMIT PSIH AND PSIM IN THE CASE OF THIN LAYERS AND HIGH ROUGHNESS            
!!---  THIS PREVENTS DENOMINATOR IN FLUXES FROM GETTING TOO SMALL                 
!        PSIH(I)=AMIN1(PSIH(I),0.9*GZ1OZ0(I))                                     
!        PSIM(I)=AMIN1(PSIM(I),0.9*GZ1OZ0(I))                                     
!  320 CONTINUE                                                                   

!-----COMPUTE THE FRICTIONAL VELOCITY:        
!     ZA(1982) EQS(2.60),(2.61).     

      DO 330 I=its,ite
        DTG=THX(I,KL)-THGB(I)        
        PSIX=GZ1OZ0(I)-PSIM(I)        
        IF((XLAND(I)-1.5).GE.0)THEN        
          ZL=ZNT(I)        
        ELSE        
          ZL=0.01        
        ENDIF        
        PSIQ=ALOG(KARMAN*UST(I)*ZA(I,KL)/XKA+ZA(I,KL)/ZL)-PSIH(I)        
        UST(I)=KARMAN*WSPD(I)/PSIX        
!        
        USTM=AMAX1(UST(I),0.1)        
        IF((XLAND(I)-1.5).GE.0)THEN        
          UST(I)=UST(I)        
        ELSE                     
          UST(I)=USTM          
        ENDIF                    
!       MOL(I,J)=KARMAN*DTG/(GZ1OZ0(I)-PSIH(I))/PRT
  330 CONTINUE                   
!                                                                                
      DO 420 I=its,ite
        WSPD1(I)=SQRT(UX(I,KL)*UX(I,KL)+VX(I,KL)*VX(I,KL))+1.E-9                  
  420 CONTINUE                                                                   
!                                                                                
!---- COMPUTE VERTICAL DIFFUSION                                                 
!                                                                                
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -          
!     COMPUTE PRELIMINARY VARIABLES                                              
!                                                                                
!                                                                                
      DT4=2.*DTSTEP                                                              
      RDT=1./DT4                                                                 
!                                                                                
      DO I=its,ite
        HGAMT(I)=0.                                                              
        HGAMQ(I)=0.                                                              
        WSCALE(I)=0.                                                             
        KPBL(I)=KL                                                               
        PBL(I)=ZQ(I,KL)                                                        
        KPBL0(I)=KL                                                    
        PBL0(I)=ZQ(I,KL)                                              
        PBLFLG(I)=.TRUE.                                                         
        SFCFLG(I)=.TRUE.                                                         
        IF(BR(I).GT.0.0)SFCFLG(I)=.FALSE.                                        
        ZL1(I)=ZA(I,KL)                                                          
        THERMAL(I)=THVX(I,KL)                                                    
      ENDDO                                                                      
                                                                                 
!     COMPUTE THE FIRST GUESS OF PBL HEIGHT                                      
                                                                                 
      DO I=its,ite
        STABLE(I)=.FALSE.                                                        
        BRUP(I)=BR(I)                                                            
      ENDDO                                                                      
      DO K=KLM,KLPBL,-1                                                          
        DO I=its,ite
          IF(.NOT.STABLE(I))THEN                                                 
            BRDN(I)=BRUP(I)                                                      
            SPDK2=MAX(UX(I,K)**2+VX(I,K)**2,1.)                                  
            BRUP(I)=(THVX(I,K)-THERMAL(I))*(G*ZA(I,K)/THVX(I,KL))/SPDK2          
            KPBL(I)=K                                                            
            STABLE(I)=BRUP(I).GT.BRCR                                            
          ENDIF                                                                  
        ENDDO                                                                    
      ENDDO                                                                      
!                                                                                
      DO I=its,ite
        K=KPBL(I)                                                                
        IF(BRDN(I).GE.BRCR)THEN                                                  
          BRINT=0.                                                               
        ELSEIF(BRUP(I).LE.BRCR)THEN                                              
          BRINT=1.                                                               
        ELSE                                                                     
          BRINT=(BRCR-BRDN(I))/(BRUP(I)-BRDN(I))                                 
        ENDIF                                                                    
        PBL(I)=ZA(I,K+1)+BRINT*(ZA(I,K)-ZA(I,K+1))                             
        IF(PBL(I).LT.ZQ(I,KPBL(I)+1))KPBL(I)=KPBL(I)+1                         
      ENDDO                                                                      
!                                                                                
      DO I=its,ite
        FM=GZ1OZ0(I)-PSIM(I)                                                     
        FH=GZ1OZ0(I)-PSIH(I)                                                     
        HOL(I)=MAX(BR(I)*FM*FM/FH,RIMIN)                                       
        IF(SFCFLG(I))THEN                                                        
          HOL(I)=MIN(HOL(I),-ZFMIN)                                          
        ELSE                                                                     
          HOL(I)=MAX(HOL(I),ZFMIN)                                           
        ENDIF                                                                    
!                                                                                
        HOL1=HOL(I)*PBL(I)/ZL1(I)*SFCFRAC                                    
        HOL(I)=-HOL(I)*PBL(I)/ZL1(I)                                       
        IF(SFCFLG(I))THEN                                                        
          PHIM(I)=(1.-APHI16*HOL1)**(-1./4.)                                     
          PHIH(I)=(1.-APHI16*HOL1)**(-1./2.)                                     
        ELSE                                                                     
          PHIM(I)=(1.+APHI5*HOL1)                                                
          PHIH(I)=PHIM(I)                                                        
        ENDIF                                                                    
        WSCALE(I)=UST(I)/PHIM(I)                                               
        WSCALE(I)=MIN(WSCALE(I),UST(I)*APHI16)                                 
        WSCALE(I)=MAX(WSCALE(I),UST(I)/APHI5)                                  
      ENDDO                                                                      
                                                                                 
!     COMPUTE THE SURFACE VARIABLES FOR PBL HEIGHT ESTIMATION                    
!     UNDER UNSTABLE CONDITIONS                                                  
                                                                                 
      DO I=its,ite
        IF(SFCFLG(I))THEN                                                        
          GAMFAC=CFAC/RHOX(I)/WSCALE(I)                                          
          HGAMT(I)=MIN(GAMFAC*HFX(I)/CPM(I),GAMCRT)                            
          HGAMQ(I)=MIN(GAMFAC*QFX(I),GAMCRQ)                                   
          IF((XLAND(I)-1.5).GE.0)HGAMQ(I)=0.                                   
          VPERT=HGAMT(I)+EP1*THX(I,KL)*HGAMQ(I)                                  
          VPERT=MIN(VPERT,GAMCRT)                                                
          THERMAL(I)=THERMAL(I)+MAX(VPERT,0.)                                    
          HGAMT(I)=MAX(HGAMT(I),0.0)                                             
          HGAMQ(I)=MAX(HGAMQ(I),0.0)                                             
        ELSE                                                                     
          PBLFLG(I)=.FALSE.                                                      
        ENDIF                                                                    
      ENDDO                                                                      
!                                                                                
      DO I=its,ite
        IF(PBLFLG(I))THEN                                                        
          KPBL(I)=KL                                                             
          PBL(I)=ZQ(I,KL)                                                      
        ENDIF                                                                    
      ENDDO                                                                      
!                                                                                
!     ENHANCE THE PBL HEIGHT BY CONSIDERING THE THERMAL                          
!                                                                                
      DO I=its,ite
        IF(PBLFLG(I))THEN                                                        
          STABLE(I)=.FALSE.                                                      
          BRUP(I)=BR(I)                                                          
        ENDIF                                                                    
      ENDDO                                                                      
      DO K=KLM,KLPBL,-1                                                          
        DO I=its,ite
          IF(.NOT.STABLE(I).AND.PBLFLG(I))THEN                                   
            BRDN(I)=BRUP(I)                                                      
            SPDK2=MAX((UX(I,K)**2+VX(I,K)**2),1.)                                
            BRUP(I)=(THVX(I,K)-THERMAL(I))*(G*ZA(I,K)/THVX(I,KL))/SPDK2          
            KPBL(I)=K                                                            
            STABLE(I)=BRUP(I).GT.BRCR                                            
          ENDIF                                                                  
        ENDDO                                                                    
      ENDDO                                                                      
!                                                                                
      DO I=its,ite
        IF(PBLFLG(I))THEN                                                        
          K=KPBL(I)                                                              
          IF(BRDN(I).GE.BRCR)THEN                                                
            BRINT=0.                                                             
          ELSEIF(BRUP(I).LE.BRCR)THEN                                            
            BRINT=1.                                                             
          ELSE                                                                   
            BRINT=(BRCR-BRDN(I))/(BRUP(I)-BRDN(I))                               
          ENDIF                                                                  
          PBL(I)=ZA(I,K+1)+BRINT*(ZA(I,K)-ZA(I,K+1))                           
          IF(PBL(I).LT.ZQ(I,KPBL(I)+1))KPBL(I)=KPBL(I)+1                       
          IF(KPBL(I).LE.1)PBLFLG(I)=.FALSE.                                      
        ENDIF                                                                    
      ENDDO                                                                      
!                                                                                
!     DIAGNOSTIC PBL HEIGHT WITH BRCR EFFECTIVELY ZERO (PBL0)                    
!                                                                                
      DO I=its,ite
        IF(PBLFLG(I))THEN                                                        
          STABLE(I)=.FALSE.                                                      
          BRUP(I)=BR(I)                                                          
        ENDIF                                                                    
      ENDDO                                                                      
      DO K=KLM,KLPBL,-1                                                          
        DO I=its,ite
          IF(.NOT.STABLE(I).AND.PBLFLG(I))THEN                                   
            BRDN(I)=BRUP(I)                                                      
            SPDK2=MAX((UX(I,K)**2+VX(I,K)**2),1.)                                
            BRUP(I)=(THVX(I,K)-THERMAL(I))*(G*ZA(I,K)/THVX(I,KL))/SPDK2          
            KPBL0(I)=K                                                           
            STABLE(I)=BRUP(I).GT.0.0                                             
          ENDIF                                                                  
                                                                                 
        ENDDO                                                                    
      ENDDO                                                                      
!                                                                                
      DO I=its,ite
        IF(PBLFLG(I))THEN                                                        
          K=KPBL0(I)                                                             
          IF(BRDN(I).GE.0.0)THEN                                                 
            BRINT=0.                                                             
          ELSEIF(BRUP(I).LE.0.0)THEN                                             
            BRINT=1.                                                             
          ELSE                                                                   
            BRINT=(0.0-BRDN(I))/(BRUP(I)-BRDN(I))                                
          ENDIF                                                                  
          PBL0(I)=ZA(I,K+1)+BRINT*(ZA(I,K)-ZA(I,K+1))                            
          IF(PBL0(I).LT.ZQ(I,KPBL0(I)+1))KPBL0(I)=KPBL0(I)+1                     
          IF(KPBL0(I).LE.1)PBLFLG(I)=.FALSE.                                     
        ENDIF                                                                    
      ENDDO                                                                      

!                                                                                
!     COMPUTE DIFFUSION COEFFICIENTS BELOW PBL                                   
!                                                                                
      DO K=kte,KLPBL,-1 
        DO I=its,ite
          IF(KPBL(I).LT.K)THEN                                                   
            PRNUM=(PHIH(I)/PHIM(I)+CFAC*KARMAN*SFCFRAC)                          
            PRNUM=MIN(PRNUM,PRMAX)                                               
            PRNUM=MAX(PRNUM,PRMIN)                                               
            ZFAC=MAX((1.-(ZQ(I,K)-ZL1(I))/(PBL(I)-ZL1(I))),ZFMIN)              
            XKZO=CKZ*DZA(I,K-1)                                                    
            XKZM(I,K)=XKZO+WSCALE(I)*KARMAN*ZQ(I,K)*ZFAC**PFAC                   
            XKZH(I,K)=XKZM(I,K)/PRNUM                                            
            XKZM(I,K)=MIN(XKZM(I,K),XKZMAX)                                      
            XKZM(I,K)=MAX(XKZM(I,K),XKZMIN)                                      
            XKZH(I,K)=MIN(XKZH(I,K),XKZMAX)                                      
            XKZH(I,K)=MAX(XKZH(I,K),XKZMIN)                                      
          ENDIF                                                                  
        ENDDO                                                                    
      ENDDO                                                                      
!                                                                                
!     COMPUTE DIFFUSION COEFFICIENTS OVER PBL (FREE ATMOSPHERE)                  
!                                                                                
      DO K=kts+1,kte
        DO I=its,ite
          XKZO=CKZ*DZA(I,K-1)                                                      
          IF(K.LE.KPBL(I))THEN                                                   
            SS=((UX(I,K-1)-UX(I,K))*(UX(I,K-1)-UX(I,K))+(VX(I,K-1)-   & 
               VX(I,K))*(VX(I,K-1)-VX(I,K)))/(DZA(I,K-1)*DZA(I,K-1))+ &          
               1.E-9                                                             
            RI=GOVRTH(I)*(THVX(I,K-1)-THVX(I,K))/(SS*DZA(I,K-1))                 
            IF(IMVDIF.EQ.1)THEN                              
              IF((QCX(I,K)+QIX(I,K)).GT.0.01E-3.AND.(QCX(I,K-1)+      & 
                QIX(I,K-1)).GT.0.01E-3)THEN                                      
!      IN CLOUD                                                                  
                QMEAN=0.5*(QX(I,K)+QX(I,K-1))                                    
                TMEAN=0.5*(SCR3(I,K)+SCR3(I,K-1))                                
                ALPH=XLV*QMEAN/R/TMEAN                                           
                CHI=XLV*XLV*QMEAN/CP/RV/TMEAN/TMEAN                              
                RI=(1.+ALPH)*(RI-G*G/SS/TMEAN/CP*((CHI-ALPH)/(1.+CHI)))          
              ENDIF                                                              
            ENDIF                                                                
            ZK=KARMAN*ZQ(I,K)                                                    
            RL2=(ZK*RLAM/(RLAM+ZK))**2                                           
            DK=RL2*SQRT(SS)                                                      
            IF(RI.LT.0.)THEN                                                     
! UNSTABLE REGIME                                                                
              SRI=SQRT(-RI)                                                      
              XKZM(I,K)=XKZO+DK*(1+8.*(-RI)/(1+1.746*SRI))                       
              XKZH(I,K)=XKZO+DK*(1+8.*(-RI)/(1+1.286*SRI))                       
            ELSE                                                                 
! STABLE REGIME                                                                  
              XKZH(I,K)=XKZO+DK/(1+5.*RI)**2                                     
              PRNUM=1.0+2.1*RI                                                   
              PRNUM=MIN(PRNUM,PRMAX)                                             
              XKZM(I,K)=(XKZH(I,K)-XKZO)*PRNUM+XKZO                              
            ENDIF                                                                
!                                                                                
            XKZM(I,K)=MIN(XKZM(I,K),XKZMAX)                                      
            XKZM(I,K)=MAX(XKZM(I,K),XKZMIN)                                      
            XKZH(I,K)=MIN(XKZH(I,K),XKZMAX)                                      
            XKZH(I,K)=MAX(XKZH(I,K),XKZMIN)                                      
          ENDIF                                                                  
!                                                                                
        ENDDO                                                                    
      ENDDO                                                                      
                                                                                 
!     COMPUTE TRIDIAGONAL MATRIX ELEMENTS FOR HEAT AND MOISTURE                  
                                                                                 
      DO I=its,ite
      DO K=kts,kte
         AU(i,k)=0.
         AL(i,k)=0.
         AD(i,k)=0.
         A1(i,k)=0.
         A2(i,k)=0.
      ENDDO
      ENDDO
 
      DO I=its,ite
        AD(I,1)=1.                                                               
        A1(I,1)=SCR3(I,KL)+HFX(I)/(RHOX(I)*CPM(I))/ZQ(I,KL)*DT4                
        A2(I,1)=QX(I,KL)+QFX(I)/(RHOX(I))/ZQ(I,KL)*DT4                         
      ENDDO                                                                      
!                                                                                
      DO K=kte,kts+1,-1
        KK=kme-K                                                                
        DO I=its,ite
          DTODSD=DT4/dz8w2d(I,K)                                                   
          DTODSU=DT4/dz8w2d(I,K-1)                                                 
          DSIG=z2d(I,K)-z2d(I,K-1)                                                 
          DSIG=-DSIG
          RDZ=1./DZA(I,K-1)                                                      
          IF(PBLFLG(I).AND.KPBL(I).LT.K)THEN                                     
            DSDZT=DSIG*XKZH(I,K)*RDZ*(G/CP-HGAMT(I)/PBL(I))                    
            DSDZQ=DSIG*XKZH(I,K)*RDZ*(-HGAMQ(I)/PBL(I))                        
            A2(I,KK)=A2(I,KK)+DTODSD*DSDZQ                                       
            A2(I,KK+1)=QX(I,K-1)-DTODSU*DSDZQ                                    
          ELSE                                                                   
            DSDZT=DSIG*XKZH(I,K)*RDZ*(G/CP)                                      
            A2(I,KK+1)=QX(I,K-1)                                                 
          ENDIF                                                                  
          DSDZ2=DSIG*XKZH(I,K)*RDZ*RDZ                                           
          AU(I,KK)=-DTODSD*DSDZ2                                                 
          AL(I,KK)=-DTODSU*DSDZ2                                                 
          AD(I,KK)=AD(I,KK)-AU(I,KK)                                             
          AD(I,KK+1)=1.-AL(I,KK)                                                 
          A1(I,KK)=A1(I,KK)+DTODSD*DSDZT                                         
          A1(I,KK+1)=SCR3(I,K-1)-DTODSU*DSDZT                                    
        ENDDO                                                                    
      ENDDO                                                                      
                                                                                 
!     SOLVE TRIDIAGONAL PROBLEM FOR HEAT AND MOISTURE                            
      
      CALL TRIDI2(AL,AD,AU,A1,A2,AU,A1,A2,                 &
                  its,ite,kts,kte                          )

!     RECOVER TENDENCIES OF HEAT AND MOISTURE                                    
                                                                                 
      DO K=kte,kts,-1
        KK=kme-K
        DO I=its,ite
          TTEND=(A1(I,KK)-SCR3(I,K))*RDT                                         
          QTEND=(A2(I,KK)-QX(I,K))*RDT                                           
          TTNP(I,K)=TTNP(I,K)+TTEND                                              
          QTNP(I,K)=QTNP(I,K)+QTEND                                              
        ENDDO                                                                    
      ENDDO                                                                      
                                                                                 
!     COMPUTE TRIDIAGONAL MATRIX ELEMENTS FOR MOMENTUM                           
                                                                                 
      DO I=its,ite
      DO K=kts,kte
         AU(i,k)=0.
         AL(i,k)=0.
         AD(i,k)=0.
         A1(i,k)=0.
         A2(i,k)=0.
      ENDDO
      ENDDO
 
      DO I=its,ite
        AD(I,1)=1.                                                               
        A1(I,1)=UX(I,KL)-UX(I,KL)/WSPD1(I)*UST(I)*UST(I)/ZQ(I,KL)  &
                *DT4*(WSPD1(I)/WSPD(I))**2                            
        A2(I,1)=VX(I,KL)-VX(I,KL)/WSPD1(I)*UST(I)*UST(I)/ZQ(I,KL)  &
                *DT4*(WSPD1(I)/WSPD(I))**2                          
      ENDDO                                                                      
!                                                                                
      DO K=kte,kts+1,-1
        KK=kme-K
        DO I=its,ite
          DTODSD=DT4/dz8w2d(I,K)                                                   
          DTODSU=DT4/dz8w2d(I,K-1)                                                 
          DSIG=z2d(I,K)-z2d(I,K-1)                                                 
          DSIG=-DSIG
          RDZ=1./DZA(I,K-1)                                                      
          DSDZ2=DSIG*XKZM(I,K)*RDZ*RDZ                                           
          AU(I,KK)=-DTODSD*DSDZ2                                                 
          AL(I,KK)=-DTODSU*DSDZ2                                                 
          AD(I,KK)=AD(I,KK)-AU(I,KK)                                             
          AD(I,KK+1)=1.-AL(I,KK)                                                 
          A1(I,KK+1)=UX(I,K-1)                                                   
          A2(I,KK+1)=VX(I,K-1)                                                   
        ENDDO                                                                    
      ENDDO                                                                      
                                                                                 
!     SOLVE TRIDIAGONAL PROBLEM FOR MOMENTUM                                     
                                                                                 
      CALL TRIDI2(AL,AD,AU,A1,A2,AU,A1,A2,                 &
                  its,ite,kts,kte                          )
                                                                                 
!     RECOVER TENDENCIES OF MOMENTUM                                             
                                                                                 
      DO K=kte,kts,-1 
        KK=kme-K
        DO I=its,ite
          UTEND=(A1(I,KK)-UX(I,K))*RDT                                           
          VTEND=(A2(I,KK)-VX(I,K))*RDT                                           
          UTNP(I,K)=UTNP(I,K)+UTEND                                              
          VTNP(I,K)=VTNP(I,K)+VTEND                                              
        ENDDO                                                                    
      ENDDO                                                                      
                                                                                 
!     COMPUTE TRIDIAGONAL MATRIX ELEMENTS FOR CLOUD                              
                                                                                 
      DO I=its,ite
      DO K=kts,kte
         AU(i,k)=0.
         AL(i,k)=0.
         AD(i,k)=0.
         A1(i,k)=0.
         A2(i,k)=0.
      ENDDO
      ENDDO
 
!     IF(IMOIST.EQ.1)GOTO 690                                                
      DO I=its,ite
        AD(I,1)=1.                                                               
        A1(I,1)=QCX(I,KL)                                                        
        A2(I,1)=QIX(I,KL)                                                        
      ENDDO                                                                      
!                                                                                
      DO K=kte,kts+1,-1
        KK=kme-K
        DO I=its,ite
          DTODSD=DT4/dz8w2d(I,K)                                                   
          DTODSU=DT4/dz8w2d(I,K-1)                                                 
          DSIG=z2d(I,K)-z2d(I,K-1)                                                 
          DSIG=-DSIG
          RDZ=1./DZA(I,K-1)                                                      
          A1(I,KK+1)=QCX(I,K-1)                                                  
          A2(I,KK+1)=QIX(I,K-1)                                                  
          DSDZ2=DSIG*XKZH(I,K)*RDZ*RDZ                                           
          AU(I,KK)=-DTODSD*DSDZ2                                                 
          AL(I,KK)=-DTODSU*DSDZ2                                                 
          AD(I,KK)=AD(I,KK)-AU(I,KK)                                             
          AD(I,KK+1)=1.-AL(I,KK)                                                 
        ENDDO                                                                    
      ENDDO                                                                      
                                                                                 
!     SOLVE TRIDIAGONAL PROBLEM FOR CLOUD                                        
                                                                                 
      CALL TRIDI2(AL,AD,AU,A1,A2,AU,A1,A2,                 &
                  its,ite,kts,kte                          )
!
      DO K=kte,kts,-1
        KK=kme-K                                                                
        DO I=its,ite
          QCTEND=(A1(I,KK)-QCX(I,K))*RDT                                         
          QITEND=(A2(I,KK)-QIX(I,K))*RDT                                         
          QCTNP(I,K)=QCTNP(I,K)+QCTEND                                           
          QITNP(I,K)=QITNP(I,K)+QITEND                                           
        ENDDO                                                                    
      ENDDO                                                                      
!                                                                                
!---- END OF VERTICAL DIFFUSION                                                  
!                                                                                
  690 CONTINUE                                                                   
!                                                                                
!-----CALCULATION OF NEW VALUES DUE TO VERTICAL EXCHANGE PROCESSES IS            
!     COMPLETED. THE FINAL STEP IS TO ADD THE TENDENCIES CALCULATED              
!     IN HIRPBL TO THOSE OF MM4.                                                 
                                                                                 
      DO 820 K=kts,kte
        NK=kme-K
        DO 820 I=its,ite
          U2DTEN(I,NK)=UTNP(I,K)
          V2DTEN(I,NK)=VTNP(I,K)
  820   CONTINUE                                                                 
!                                                                                
!     IF(J.EQ.1.AND.IN.GT.1)GOTO 860                                             
!SUE  JBGN=3                                                                     
!SUE  JEND=JLXM-1                                                                

! change when nest
!SUE  JBGN=2                                                                   
!SUE  JEND=JLXM                                                                

      JBGN=jts
      JEND=jte
      IBGN=its                                                                   
      IEND=ite

!     IF(J.LT.JBGN.OR.J.GT.JEND)GOTO 860                                         
!SUE  IBGN=3                                                                     
!SUE  IEND=ILXM-1                                                                

! change when nest
!SUE  IBGN=2                                                                   
!SUE  IEND=ILXM                                                                

      DO 830 K=kts,kte
        NK=kme-K
        DO 830 I=IBGN,IEND                                                       
          T2DTEN(I,NK)=TTNP(I,K)
  830   CONTINUE                                                                 
!                                                                                
!     IF(IDRY.EQ.1)GOTO 860                                                  
      DO 840 K=kts,kte
        NK=kme-K
        DO 840 I=IBGN,IEND                                                       
          QV2DTEN(I,NK)=QTNP(I,K) 
  840   CONTINUE                                                                 
                                                                                 
!     IF(IMOIST.EQ.1)GOTO 860
      DO 850 K=kts,kte
        NK=kme-K
        DO 850 I=IBGN,IEND
           QC2DTEN(I,NK)=QCTNP(I,K)
  850   CONTINUE

      IF(flag_QI .AND. PRESENT( QI2DTEN ) ) THEN
        DO K=kts,kte
        NK=kme-K
          DO I=IBGN,IEND
             QI2DTEN(I,NK)=QITNP(I,K)
          ENDDO
        ENDDO
      ENDIF

  860 CONTINUE                                                                   
!                                                                                
!-----APPLY ASSELIN FILTER TO TGD FOR LARGE TIME STEP:                           
!                                                                                
!     DO 885 I=its,ite
!       TSK(I)=TSK(I)*(PS(I)/100.)**ROVCP                                    
! 885 CONTINUE                                                                   
!                                                                                
  940 CONTINUE                                                                   
!                                                                                
! KPBL IS NEEDED FOR THE FDDA, AND SINCE THERE IS NO LONGER JUST ONE             
! LARGE "J LOOP" IT MUST BE STORED AS (I,J)...                                   
!                                                                                
! USE NEW DIAGNOSED PBL DEPTH (CALCULATED WITH brcr=0.0)
! PBL IS USED FOR OUTPUT AND NEXT-TIME-STEP BELJAARS CALC IN SFCLAY
      DO 950 I=its,ite
        KPBL1D(I)=KPBL0(I)                                                      
        PBL(I)=PBL0(I)
  950 CONTINUE                                                                   

   END SUBROUTINE MRF2D

!================================================================          
   SUBROUTINE TRIDI2(CL,CM,CU,R1,R2,AU,A1,A2,                   &
                     its,ite,kts,kte                            )
!----------------------------------------------------------------
   IMPLICIT NONE
!----------------------------------------------------------------

   INTEGER, INTENT(IN )      ::     its,ite, kts,kte
                                   
   REAL, DIMENSION( its:ite, kts+1:kte+1 )                    , &
         INTENT(IN   )  ::                                  CL

   REAL, DIMENSION( its:ite, kts:kte )                        , &
         INTENT(IN   )  ::                                  CM, &
                                                            R1, &
                                                            R2
   REAL, DIMENSION( its:ite, kts:kte )                        , &
         INTENT(INOUT)  ::                                  AU, &
                                                            CU, &
                                                            A1, &
                                                            A2

   REAL    :: FK
   INTEGER :: I,K,L,N

!----------------------------------------------------------------         

   L=ite 
   N=kte

   DO I=its,L
     FK=1./CM(I,1)                                                            
     AU(I,1)=FK*CU(I,1)                                                       
     A1(I,1)=FK*R1(I,1)                                                       
     A2(I,1)=FK*R2(I,1)                                                       
   ENDDO                                                                      
   DO K=2,N-1
     DO I=its,L
       FK=1./(CM(I,K)-CL(I,K)*AU(I,K-1))                                      
       AU(I,K)=FK*CU(I,K)                                                     
       A1(I,K)=FK*(R1(I,K)-CL(I,K)*A1(I,K-1))                                 
       A2(I,K)=FK*(R2(I,K)-CL(I,K)*A2(I,K-1))                                 
     ENDDO                                                                    
   ENDDO                                                                      
   DO I=its,L
     FK=1./(CM(I,N)-CL(I,N)*AU(I,N-1))                                        
     A1(I,N)=FK*(R1(I,N)-CL(I,N)*A1(I,N-1))                                   
     A2(I,N)=FK*(R2(I,N)-CL(I,N)*A2(I,N-1))                                   

   ENDDO                                                                      
   DO K=N-1,kts,-1                                                              
     DO I=its,L
       A1(I,K)=A1(I,K)-AU(I,K)*A1(I,K+1)                                      
       A2(I,K)=A2(I,K)-AU(I,K)*A2(I,K+1)                                      
     ENDDO                                                                    
   ENDDO                                                                      

   END SUBROUTINE TRIDI2
      
!===================================================================
   SUBROUTINE mrfinit(RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN,           &
                      RQCBLTEN,RQIBLTEN,P_QI,P_FIRST_SCALAR,       &
                      restart, allowed_to_read ,                   &
                      ids, ide, jds, jde, kds, kde,                &
                      ims, ime, jms, jme, kms, kme,                &
                      its, ite, jts, jte, kts, kte                 )
!-------------------------------------------------------------------          
   IMPLICIT NONE
!-------------------------------------------------------------------          
   LOGICAL , INTENT(IN)          :: restart , allowed_to_read
   INTEGER , INTENT(IN)          ::  ids, ide, jds, jde, kds, kde, &
                                     ims, ime, jms, jme, kms, kme, &
                                     its, ite, jts, jte, kts, kte
   INTEGER , INTENT(IN)          ::  P_QI,P_FIRST_SCALAR

   REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) ::         &
                                                         RUBLTEN, &
                                                         RVBLTEN, &
                                                         RTHBLTEN, &
                                                         RQVBLTEN, &
                                                         RQCBLTEN, & 
                                                         RQIBLTEN
   INTEGER :: i, j, k, itf, jtf, ktf

   jtf=min0(jte,jde-1)
   ktf=min0(kte,kde-1)
   itf=min0(ite,ide-1)

   IF(.not.restart)THEN
     DO j=jts,jtf
     DO k=kts,ktf
     DO i=its,itf
        RUBLTEN(i,k,j)=0.
        RVBLTEN(i,k,j)=0.
        RTHBLTEN(i,k,j)=0.
        RQVBLTEN(i,k,j)=0.
        RQCBLTEN(i,k,j)=0.
     ENDDO
     ENDDO
     ENDDO
   ENDIF

   IF (P_QI .ge. P_FIRST_SCALAR .and. .not.restart) THEN
      DO j=jts,jtf
      DO k=kts,ktf
      DO i=its,itf
         RQIBLTEN(i,k,j)=0.
      ENDDO
      ENDDO
      ENDDO
   ENDIF

   END SUBROUTINE mrfinit

!-------------------------------------------------------------------          

END MODULE module_bl_mrf

