      SUBROUTINE VRTDF22(     TIO,    QIO,    UIO,    VIO,    XIO,
     1                        UTG,    VTG,
     2                    UFSIROL,VFSIROL,UFSOROL,VFSOROL,     
     3                    UFSROL, VFSROL,    UFS,    VFS,
     +                    UFSINST, VFSINST, 
     4                       ALMX,   ALMC,   PBLT, 
     5                     CNDROL, DEPROL,   WSUB, RKMROL,RKHROL,RKQROL,
C
C    -------------- OUTPUTS OR UPDATED INPUTS ARE ABOVE THIS LINE,
C    -------------- INPUTS ARE BELOW.
C
     6                      GTROT,  QGROT, CDMROT, QFSROT, HFSROT,     
     7                      GTAGG,  QGAGG, CDMAGG, QFSAGG, HFSAGG,
     8                    FAREROT, PRESSG,    CRH,   ZSPD,   FLND,
     9                    CQFXROW,CHFXROW,UTENDCV,VTENDCV,
     A                    UTENDGW,VTENDGW,   TSGB,     TF,    SGJ,
     B                       SGBJ, SHTXKJ,  SHXKJ,   SHTJ,   SHBJ, 
     C                        SHJ,   DSGJ,   DSHJ,   CVSG,  ZCLFC,  
     D                     ITRPHS,   ILWC,   IIWC,  IOWAT,  IOSIC,
     E                      ZTMST,    ILG,    IL1,    IL2,  NTILE,
     F                       ILEV,    LEV,   LEVS,  NTRAC,   IPAM,
     G                     ISAVLS,SAVERAD,SAVEBEG                )
C
C     * MAY 03/2018 - N.MCFARLANE.   Added code to prevent STATCLD crashes.
C     *                              Removed unused {ESTEFF,HTVOCP} statement functions.
C     * FEB 24/2018 - M.LAZARE.      Bugfix to ensure convective contribution
C     *                              is seen in stresses passed to coupler:
C     *                              - Arrays {UTENDCV,VTENDCV} passed in from
C     *                                physics and used.
C     *                              - New internal arrays {UCVTS,VCVTS} defined
C     *                                and used, to hold vertically-integrated tendencies
C     *                                due to convection, analagous to {UGWTS,VGWTS}.
C     * JUN 22/2016 - M.LAZARE.      New tiled version for gcm19:
C     *                              - {GT,QG,QFS,HFS,CDM} are now tiled
C     *                                input and {UFSROL,VFSROL} are
C     *                                tiled output. Number of tiles
C     *                                (NTILE) and tile fractional area
C     *                                (FAREROT) are passed in.
C     *                              - {TH,Q,U,V,XROW} are now internal
C     *                                work arrays, set to incoming
C     *                                values {TIO,QIO,UIO,VIO,XIO} at
C     *                                the start of the tiling loop. For
C     *                                each pass through the tiling loop,
C     *                                tendencies are calculated for
C     *                                each tile and after the tiling loop,
C     *                                new values of {TIO,QIO,UIO,VIO,XIO}
C     *                                are determined by aggregating over
C     *                                the tiled tendencies weighted by
C     *                                the tile fractional area.
C     *                              - {UTG,VTG} are new output fields passed
C     *                                out after taking into account 
C     *                                contributions from GWD and convection
C     *                                (in incoming {UTENDGW,VTENDGW} and
C     *                                this routine. The setting of
C     *                                {UTENDGW,VTENDGW} to {UTG,VTG} in the
C     *                                physics is therefore removed.
C     *                              - the diagnostic {UFS,VFS} is determined
C     *                                by also aggregating over the tile
C     *                                stresses and then the {UGWTS,VGWTS}
C     *                                contribution is added in.
C     *                              - all non-tile calculations are done
C     *                                first, before the tiled loop starts.
C     *                              - we have left in the option to revert
C     *                                to the old way with no tiling. Under
C     *                                that option, defined by ITILVRT=0,
C     *                                we use the non-tiled input fields
C     *                                and set the tile loop index to one.
C     *                              - RKXMIN now local and not passed in.
C     *                              - "PARAMS","PARAM1" AND "PARAM3" common
C     *                                blocks added, so {TFREZ,HS,HV,PI,CPRES}
C     *                                removed from call.
C     *                              - ITRAC removed from call after two
C     *                                IF conditional sections changed to always
C     *                                be true.
C     *                              - SAVEGG changed to real name SAVERAD
C     *                                (variable name in call).
C     *                              - {HFS,QFS} used directly instead of
C     *                                copied {SHFLX,SQFLX} in nlclmx8 call,
C     *                                so these two work arrays removed.
C     *                              - {PBLT,ALMX,ALMC} aggregated over tiles.
C     *                              - {XLWROL,XICROL} now internal work arrays.
C     * FEB 10/2015 - M.LAZARE/      Previous version vrtdf21 for gcm18:
C     *               K.VONSALZEN:   - Comment-out lower bound on total
C     *                                water (don't think we need it
C     *                                any more due to implemented
C     *                                improvements and bugfixes elsewhere).
C     *                              - Revised call to new NLCLMX8.
C     * JUL 31/2013 - M.LAZARE.      Previous version VRTDF20 for gcm17:
C     *                              - STATCLD5 not done above 10 Hpa.
C     *                              - Print out warning if QTN<0 before
C     *                                calling statcld5, then limit
C     *                                it to be non-negative.
C     *                              - Pass in LVL and ICALL to statcld5
C     *                                to aid in future debugging.
C     * JUN 26/2013 - K.VONSALZEN/   Previous version for gcm17:
C     *               M.LAZARE.      - Revise calls for new STATCLD5 and
C     *                                NLCLMX7.
C     *                              - Add the calculation of the subgrid-
C     *                                scale component of the vertical
C     *                                velocity (to determine aerosol
C     *                                activation). This is passed out
C     *                                as new i/o array WSUB.
C     *                              - Pass in IPAM to determine if bulk
C     *                                or pla aerosols are being used and
C     *                                choose between different calculations
C     *                                of cloudy asymptotic mixing length
C     *                                (ALMC) depending on choice of IPAM.
C     *                              - Don't call statcld above 10 mbs to
C     *                                avoid spurious crashes.
C     * APR 28/2012 - K.VONSALZEN/   PREVIOUS VERSION VRTDF19 FOR GCM16:
C     *               M.LAZARE:      - RKMIN,RKQMIN INCREASED FROM 
C     *                                0.01 TO 0.1.
C     *                              - STATISTICAL CLOUD SCHEME ONLY
C     *                                CALLED FOR L>LEV1 (DEFINED IN RAD
C     *                                AND PASSED IN ITOPLW COMMON BLOCK)
C     *                                AND EST<1.01*P.
C     *                              - REVISION OF SURFACE WIND STRESS
C     *                                CALCULATION TO REMOVE DOUBLE
C     *                                COUNTING OF VERTICAL DIFFUSION
C     *                                GOING INTO {UFS,VFS}.
C     *                              - CALLS NEW NLCLMX6 and STATCLD4.
C     * APR 21/2010 - M.LAZARE.      PREVIOUS VERSION VRTDF18 FOR GCM15I:
C     *                              - REMOVE UNUSED CALCULATION OF "OMET"
C     *                                AND ACCORDINGLY FROM PHYSICS SO
C     *                                IT CAN BE USED IN THE NEW CORE15PI
C     *                                TO SAVE THE NEW VERTICAL VELOCITY
C     *                                DIAGNOSTIC FIELD.  
C     * FEB 17/2009 - K.VONSALZEN/   PREVIOUS VERSION VRTDF17 FOR GCM15H:
C     *               M.LAZARE.      - MINIMUM RK'S REDUCED BY ORDER
C     *                                OF MAGNITUDE.
C     *                              - UNIFIED CLEAR/CLOUD MIXING LENGTHS.
C     * JAN 17/2008 - K.VONSALZEN/   PREVIOUS VERSION VRTDF16 FOR GCM15G:
C     *               M.LAZARE.      - CALLS NEW NLCLMX5 AND STATCLD3.
C     *                              - BUGFIX FOR CALCULATION OF ZF.
C     *                              - USES ROECKEL FRACTIONAL PROBABILITY
C     *                                OF WATER PHASE FUNCTION.
C     *                              - USES CLEAR-SKY AND ALL-SKY MIXING
C     *                                LENGTHS (ALMC AND ALMX, RESPECTIVELY)
C     *                                TO DEFINE A CONSISTENT TOTAL
C     *                                MIXING LENGTH (ALMIX) PASSED TO
C     *                                STATCLD3.
C     *                              - PASSES IN ADELT=2.*DELT FROM PHYSICS
C     *                                AS ZTMST WHICH IS USED THROUGHOUT
C     *                                ROUTINE IN PLACE OF HARD-CODED 2.*DELT.
C     *                              - CONSISTENT PBL USEAGE (LPBL=NINT(PBLT)). 
C     *                              - GUSTINESS EFFECT NOW INCLUDED 
C     *                                SINCE IS NOW ALREADY CONTAINED IN ZSPD
C     *                                FROM PHYSICS.
C     * JAN 13/2007 - K.VONSALZEN.   PREVIOUS VERSTION VRTDF15 FOR GCM15F.
C     *                              - MODIFIED CALL TO STATCLD2, IN
C     *                                CONJUNCTION WITH CHANGES TO ADD
C     *                                "QCWVAR".  
C     * NOV 28/2006 - M.LAZARE/      - ADD CALCULATION FOR CNDROL,DEPROL
C     *               K.VONSALZEN.     UNDER CONTROL OF "ISAVLS".
C     *                              - MOVE MIXING OF TRACERS TO BEFORE 
C     *                                CALCULATION OF MLSE AND TOTAL
C     *                                WATER, SO PROFILES ARE WELL-MIXED
C     *                                BEFORE BEING PROCESSED. 
C     *                                NOTE THAT NOW CLOUD WATER AND ICE
C     *                                ARE MIXED AS WELL, WHICH REQUIRED
C     *                                THAT RKX HAD TO BE DEFINED
C     *                                REGARDLESS OF WHETHER TRACER IS
C     *                                ADVECTED OR NOT. 
C     * NOV 24/2006 - M.LAZARE.      - REMOVE XFS FROM ROUTINE, AND
C     *                                DON'T PASS OUT RKMIN,RKQMIN
C     *                                (ALL NOT NEEDED). 
C     * JUL 30/2006 - M.LAZARE.      - BUGFIX: "IL"->"I" IN CALCULATION 
C     *                                OF ALWC.
C     * JUN 30/2006 - M.LAZARE.      - UPDATE AND PASS OUT {UTENDGW,VTENDGW}
C     *                                INSTEAD OF {UTG,VTG} SINCE THESE
C     *                                LATTER (TOTAL) TENDENCIES ARE
C     *                                CALCULATED SUBSEQUENTLY NOW IN THE
C     *                                PHYSICS DRIVER. 
C     *                              - USE VARIABLE INSTEAD OF 
C     *                                CONSTANT IN INTRINSICS SUCH AS "MAX",
C     *                                SO THAT CAN COMPILE IN 32-BIT MODE
C     *                                WITH REAL*8.  
C     *                              - CALLS NEW STATCLD2 AND NLCLMX4.
C     *                              - WORK ARRAYS NOW LOCAL.
C     *                              - DEFINE ALMIX,DHLDZ,DRWDZ
C     *                                OUTSIDE OF CONDITIONAL TEST ON
C     *                                ES<P, TO AVOID NAN'S ENTERING
C     *                                STATCLD2.
C     * MAY 06/2006 - K.VONSALZEN/   PREVIOUS VERSION VRTDF14 FOR GCM15E:
C     *               M.LAZARE.      - MIXES TOTAL WATER AND LIQUID WATER
C     *                                STATIC ENERGY RATHER THAN 
C     *                                TEMPERATURE AND MOISTURE, THEN
C     *                                UNRAVELS CONSISTENTLY VIA THE
C     *                                NEW COMMON SUBROUTINE STATCLD TO
C     *                                GET ALL THERMODYNAMIC FIELDS.
C     *                              - CALLS NEW NLCLMX3 INSTEAD OF NLCLMX2.
C     *                              - MIXING IS NOW DONE FOR COMPLETE
C     *                                VERTICAL DOMAIN, SO "MSGPHYS" IS
C     *                                NO LONGER PASSED IN; RATHER "MSG=0"
C     *                                IS DEFINED AND PASSED TO NEW
C     *                                NLCLMX3.
C     * FEB 03/2006 - K.VONSALZEN/   PREVIOUS VERSION VRTDF13 FOR GCM15D.
C     *               M.LAZARE. 
C                                                                              
C     * CALCULATE THE TENDENCIES OF U,V,LWSE,QT,X DUE TO VERTICAL DIFFUSION.
C     *                                                                        
C     *            MOMENTUM                         THERMODYNAMICS             
C     *            ========                         ==============             
C     *                                   D2SG(1)=0.                           
C     * SIGMA=0. /////////////////////////////////////////////////// SIGMA=0.  
C     *                A                B                A                     
C     * SG (1) . . . . A . . . . . . . .B                A                     
C     *                A=DSG(1)         B=D1SG(1)        A=DSH(1)              
C     *                A                B. . . . . . . . A . . . . . SH (1)    
C     * SGB(1) --------------------------                A                     
C     *                B                A=D2SG(2)        A                     
C     *                B                ---------------------------- SHB(1)    
C     * SG (2) . . . . B . . . . . . . .B                B                     
C     *                B=DSG(2)         B=D1SG(2)        B=DSH(2)              
C     *                B                B. . . . . . . . B . . . . . SH (2)    
C     * SGB(2) --------------------------                B                     
C     *                C                A=D2SG(3)        B                     
C     *                C                ---------------------------- SHB(2)    
C     * SG (3) . . . . C . . . . . . . .B                C                     
C     *                C=DSG(3)         B=D1SG(3)        C=DSH(3)              
C     *                C                B. . . . . . . . C . . . . . SH (3)    
C     * SGB(3) --------------------------                C                     
C     *                D                A=D2SG(4)        C                     
C     *                D                ---------------------------- SHB(3)    
C     * SG (4) . . . . D . . . . . . . .B                D                     
C     *                D=DSG(4)         B=D1SG(4)        D=DSH(4)              
C     *                D                B. . . . . . . . D . . . . . SH (4)    
C     * SGB(4) --------------------------                D                     
C     *                E                A=D2SG(5)        D                     
C     *                E                ---------------------------- SHB(4)    
C     * SG (5) . . . . E . . . . . . . .B                E                     
C     *                E=DSG(5)         B=D1SG(5)        E=DSH(5)              
C     *                E                B. . . . . . . . E . . . . . SH (5)    
C     *                E                B                E                     
C     * SGB(5)=1. ////////////////////////////////////////////////// SHB(5)=1. 
C     *                                                                       
      IMPLICIT NONE
C                       
C     * I/O ARRAYS.                                                           
C
      REAL   TIO    (ILG,ILEV),  QIO    (ILG,ILEV)
      REAL   UIO    (ILG,ILEV),  VIO    (ILG,ILEV)
      REAL   XIO    (ILG,LEV,NTRAC)

      REAL   UTG    (ILG,ILEV),  VTG    (ILG,ILEV)

      REAL   UFSIROL(ILG),       VFSIROL(ILG)
      REAL   UFSOROL(ILG),       VFSOROL(ILG)

      REAL   UFSROL (ILG),       VFSROL (ILG)                  
      REAL   UFSINST (ILG),      VFSINST (ILG)
      REAL   UFS    (ILG),       VFS    (ILG)
      REAL   PBLT   (ILG)

      REAL   ALMX   (ILG,ILEV),  ALMC   (ILG,ILEV)
      REAL   CNDROL (ILG,ILEV),  DEPROL (ILG,ILEV),  WSUB   (ILG,ILEV)
      REAL   RKMROL(ILG,ILEV), RKHROL(ILG,ILEV), RKQROL(ILG,ILEV)
C
C     * INPUT FIELDS.
C
      REAL   GTROT  (ILG,NTILE), QGROT  (ILG,NTILE)
      REAL   CDMROT (ILG,NTILE), QFSROT (ILG,NTILE)
      REAL   HFSROT (ILG,NTILE), FAREROT(ILG,NTILE)

      REAL   GTAGG  (ILG),       QGAGG  (ILG)
      REAL   CDMAGG (ILG),       QFSAGG (ILG),     HFSAGG  (ILG)
      REAL   PRESSG (ILG),       CRH    (ILG),     ZSPD    (ILG)
      REAL   CQFXROW(ILG),       CHFXROW(ILG),     FLND    (ILG)

      REAL   UTENDGW(ILG,ILEV),  VTENDGW(ILG,ILEV)
      REAL   UTENDCV(ILG,ILEV),  VTENDCV(ILG,ILEV)
      REAL   TSGB   (ILG,ILEV),  TF     (ILG,ILEV)
      REAL   SGJ    (ILG,ILEV),  SGBJ   (ILG,ILEV)
      REAL   SHTXKJ (ILG,ILEV),  SHXKJ  (ILG,ILEV)
      REAL   SHTJ   (ILG, LEV),  SHBJ   (ILG,ILEV)
      REAL   SHJ    (ILG,ILEV),  DSGJ   (ILG,ILEV)
      REAL   DSHJ   (ILG,ILEV),  CVSG   (ILG,ILEV)
      REAL   ZCLFC  (ILG,ILEV)
                                                                              
      INTEGER ITRPHS(NTRAC)
C                                                                           
C     * INTERNAL WORK ARRAYS.
C
      REAL  ,  DIMENSION(ILG, LEV,NTRAC) :: XROW,CHGX
      REAL  ,  DIMENSION(ILG,ILEV,NTRAC) :: RKX
      REAL  ,  DIMENSION(ILG,ILEV)  :: TH,Q,U,V,CHGT,CHGQ,CHGU,CHGV
      REAL  ,  DIMENSION(ILG,ILEV)  :: DVDS,TEND,DTTDS,RI,A,B,C
      REAL  ,  DIMENSION(ILG,ILEV)  :: WORK,RKM,RKQ,RKH,ZF,P
      REAL  ,  DIMENSION(ILG,ILEV)  :: RDZ,CHI,Z,QT,RRL,HMN
      REAL  ,  DIMENSION(ILG,ILEV)  :: QL,DST,SIGMA,DSR,CPM,PF
      REAL  ,  DIMENSION(ILG,ILEV)  :: ALMXT,ALMXIN,ALMCIN,ALMXX,ALMCX
      REAL  ,  DIMENSION(ILG,ILEV)  :: CNDROLX,DEPROLX,XLWROL,XICROL
      REAL  ,  DIMENSION(ILG,ILEV)  :: WSUBX
      REAL  ,  DIMENSION(ILG,NTRAC) :: XINT,SXFLX,XREF

      REAL  ,  DIMENSION(ILG)   :: GT,QG,CDM,QFS,HFS,FARE
      REAL  ,  DIMENSION(ILG)   :: DSMIX,HINT,CDVLH,CDVLM,VMODL,QINT
      REAL  ,  DIMENSION(ILG)   :: QTG,BCR,HMNG,CL,DVDZ,X,PBLTX
      REAL  ,  DIMENSION(ILG)   :: HEAT,RAUS,VINT,XINGLF,TVREF
      REAL  ,  DIMENSION(ILG)   :: CDVLT,ZER,UN,WRKRL,WRKR
      REAL  ,  DIMENSION(ILG)   :: UCVTS,VCVTS,UGWTS,VGWTS
      REAL  ,  DIMENSION(ILG)   :: QCW,ZCLF,ZCRAUT,SSH
      REAL  ,  DIMENSION(ILG)   :: ALMIX,DHLDZ,DRWDZ,ATUNE
      REAL  ,  DIMENSION(NTRAC) :: RKXMIN

      INTEGER, DIMENSION(ILG) :: IPBL,IPBLC
C
C     * SCALARS PASSED IN CALL:
C
      REAL      ZTMST,SAVEBEG,SAVERAD
      INTEGER   ILWC,IIWC,ILG,IL1,IL2,NTILE,ILEV,LEV,LEVS,NTRAC
      INTEGER   IOWAT,IOSIC,IPAM,ISAVLS
C
C     * LOCAL SCALARS.
C
      REAL      ZSECFRL,ZTHOMI,TAUADJ,ROG,ALWC,THETVP,THETVM,PBLX,
     1          FACT,RINEUT,RIINF,FACT0,SARI,ALL,XIMIN,XIMINT,
     2          ALU,ALD,FAC,XINGLH,XINGLM,ELH,ELM,EPSSH,EPSSM,
     3          FACTH,FACTM,XH,XM,ATMP,PRANDTL,SCALF,FACTS,RKMN,
     4          RKHMN,RKQMN,TODT,EST,QCWX,DZ,FRACW,FACMOM,UTMP,VTMP,
     5          QCWI,QCWL,FARETOT,ACOR,ACORT,ACOR1
      INTEGER   NTILEND,L,I,IL,ILEVM,LPBL,N,NT,ICVSG,ISUBG,IDUM,ICALL
C
C     * COMMON BLOCK SCALARS.
C     
      REAL WW,TWX,AX,ASQ,GRAV,RGAS,RGOCP,RGOASQ,CPRES,RGASV,CPRESV
      COMMON /PARAMS/ WW,     TWX,   AX,     ASQ,  GRAV, RGAS,  RGOCP,
     1                RGOASQ, CPRES, RGASV, CPRESV

      REAL PI,RVORD,TFREZ,HS,HV,DAYLNT
      COMMON /PARAM1/ PI,     RVORD, TFREZ, HS,   HV,   DAYLNT

      REAL RGASX,RGASVX,GRAVX,SBC,VKC,CT,VMIN
      COMMON /CLASS2/ RGASX,RGASVX,GRAVX,SBC,VKC,CT,VMIN

      REAL AP,BP,EPS1,EPS2
      COMMON /EPS   / AP,BP,EPS1,EPS2

      REAL T1S,T2S,AI,BI,AW,BW,SLP
      COMMON /HTCP   / T1S, T2S, AI, BI, AW, BW, SLP      

      REAL RW1,RW2,RW3,RI1,RI2,RI3
      COMMON /ESTWI/ RW1,RW2,RW3,RI1,RI2,RI3
C
C     * THERE MUST BE NO CLOUD WITHIN FIRST "LEV1" LAYERS FOR RADIATION
C     * TO WORK PROPERLY. THIS IS DEFINED IN THE "TOPLW" SUBROUTINE
C     * CALLED AT THE BEGINNING OF THE MODEL.
C
      INTEGER LEV1
      COMMON /ITOPLW/ LEV1
C
      REAL    GAMRH,GAMRM,BEEM,ALFAH,RAYON,DVMINS,TVFA,TICE
      INTEGER MSG
      DATA      GAMRH,      GAMRM,      BEEM,       ALFAH
     1       /  6.00,       6.00,       10.0,       1.0    /                  
      DATA RAYON/6.37122E06/,  DVMINS/7.9/                              
      DATA TVFA / 0.608 /, TICE/253./, MSG/0/
C
C     * STATEMENT FUNCTION TO CALCULATE SATURATION VAPOUR PRESSURE
C     * OVER WATER OR ICE.
C
      REAL    TTT,UUU,ESW,ESI
      ESW(TTT)    = EXP(RW1+RW2/TTT)*TTT**RW3
      ESI(TTT)    = EXP(RI1+RI2/TTT)*TTT**RI3
C
C     * COMPUTES THE RATIO OF LATENT HEAT OF VAPORIZATION OF     
C     * WATER OR ICE TO THE SPECIFIC HEAT OF AIR AT CONSTANT
C     * PRESSURE CP.     
C
      REAL    TW,TI
      TW(TTT)     = AW-BW*TTT  
      TI(TTT)     = AI-BI*TTT  
C
      REAL    ZERO,ONE,PRANDTL_MIN,PRANDTL_MAX,XLMIN
      DATA ZERO,ONE /0., 1./
      DATA PRANDTL_MIN,PRANDTL_MAX /1., 3./
      DATA XLMIN /10./
C
C     * SWITCH TO CONTROL TILING.
C
      INTEGER ITILVRT
      DATA ITILVRT /1/
C--------------------------------------------------------------------          
C     * SET TILING LOOP INDEX BASED ON CHOICE OF "ITILVRT".
C
      IF(ITILVRT.EQ.1) THEN
         NTILEND=NTILE
      ELSE
         NTILEND=1
      ENDIF
C
      ROG=RGAS/GRAV 
      ILEVM=ILEV-1                                                    
C
      DO IL=IL1,IL2
        UN   (IL)=1.                                                          
        ZER  (IL)=0.    
        VMODL(IL)=MAX(VMIN,ZSPD(IL))    
      ENDDO
C
C     * CALCULATE LOCAL PRESSURE (MBS) AND HEIGHTS (M).
C                             
      DO 10 IL=IL1,IL2                                                  
        P (IL,1) = SHJ (IL,1)*PRESSG(IL)*0.01                       
        PF(IL,1) = SHTJ(IL,1)*PRESSG(IL)*0.01
   10 CONTINUE

      DO 20 L=2,ILEV                                                    
      DO 20 IL=IL1,IL2                                                  
        P (IL,L) = SHJ (IL,L  )*PRESSG(IL)*0.01                   
        PF(IL,L) = SHBJ(IL,L-1)*PRESSG(IL)*0.01  
   20 CONTINUE                                                          
C
      DO 30 IL=IL1,IL2                                                  
        Z (IL,ILEV) = ROG*GTAGG(IL)   *LOG(0.01*PRESSG(IL)/P (IL,ILEV))  
        ZF(IL,ILEV) = ROG*TIO(IL,ILEV)*LOG(0.01*PRESSG(IL)/PF(IL,ILEV))  
   30 CONTINUE                                                          

      DO 40 L=ILEVM,1,-1
      DO 40 IL=IL1,IL2                                                  
        Z (IL,L)=Z (IL,L+1)+ROG*TF(IL,L+1)*LOG(P (IL,L+1)/P (IL,L))   
        ZF(IL,L)=ZF(IL,L+1)+ROG*TIO(IL,L) *LOG(PF(IL,L+1)/PF(IL,L))   
   40 CONTINUE
C
C     * PARAMETERS FOR MOIST PHYSICS.
C
      ZSECFRL=1.E-7
      ZTHOMI=238.16
      DO L=1,ILEV
      DO IL=IL1,IL2 
        DSR(IL,L)=1.
        IF( XIO(IL,L+1,IIWC).GT.ZSECFRL ) THEN
          CHI(IL,L)=1./(1.+XIO(IL,L+1,ILWC)/XIO(IL,L+1,IIWC))
        ELSE
          IF ( XIO(IL,L+1,ILWC).GT.ZSECFRL ) THEN
            CHI(IL,L)=0.
          ELSE
C 
C           * COMPUTE THE FRACTIONAL PROBABILITY OF WATER PHASE      
C           * EXISTING AS A FUNCTION OF TEMPERATURE (FROM ROCKEL,     
C           * RASCHKE AND WEYRES, BEITR. PHYS. ATMOSPH., 1991.)       
C
            FRACW = MERGE( 1.,    
     1                  0.0059+0.9941*EXP(-0.003102*(T1S-TIO(IL,L))**2),
     2                  TIO(IL,L).GE.T1S )  
            CHI(IL,L)=1.-FRACW
          ENDIF
        ENDIF 
        RRL(IL,L)=(1.-CHI(IL,L))*HV+CHI(IL,L)*HS
        CPM(IL,L)=CPRES
      ENDDO
      ENDDO
C
C     * CALCULATE VERTICALLY-INTEGRATED GRAVITY-WAVE DRAG EFFECTS.
C     * INITIALIZE OTHER FIELDS.
C
      DO I=IL1,IL2
         UGWTS(I)=0.
         VGWTS(I)=0.
         UCVTS(I)=0.
         VCVTS(I)=0.
         RAUS (I)=PRESSG(I)/GRAV
         CDVLH(I)=0.
         CDVLT(I)=0.
         PBLT (I)=0.
      ENDDO
C
      DO L=1,ILEV
      DO I=IL1,IL2
         UGWTS(I) = UGWTS(I)+DSGJ(I,L)*UTENDGW(I,L)*RAUS(I)
         VGWTS(I) = VGWTS(I)+DSGJ(I,L)*VTENDGW(I,L)*RAUS(I)
         UCVTS(I) = UCVTS(I)+DSGJ(I,L)*UTENDCV(I,L)*RAUS(I)
         VCVTS(I) = VCVTS(I)+DSGJ(I,L)*VTENDCV(I,L)*RAUS(I)
         UTG(I,L) = UTENDGW(I,L)+UTENDCV(I,L)
         VTG(I,L) = VTENDGW(I,L)+VTENDCV(I,L)
      ENDDO
      ENDDO
C
C     * INITIALIZE DIAGNOSED STRESSES TO THAT FROM GWD AND CONVECTION ABOVE.
C
      DO I=IL1,IL2
         UFS   (I)  = UFS   (I) + (UGWTS(I)+UCVTS(I))*SAVERAD
         VFS   (I)  = VFS   (I) + (VGWTS(I)+VCVTS(I))*SAVERAD
         UFSROL(I)  = UFSROL(I) + UCVTS(I)*SAVEBEG
         VFSROL(I)  = VFSROL(I) + VCVTS(I)*SAVEBEG
         UFSINST(I) = UCVTS(I)
         VFSINST(I) = VCVTS(I)
      ENDDO
C
C     * INITIALIZE "CHG" ARRAYS.
C
      DO N=1,NTRAC
      DO L=1,LEV
      DO IL=IL1,IL2
        CHGX(IL,L,N)= 0.
      ENDDO
      ENDDO
      ENDDO

      DO L=1,ILEV
      DO IL=IL1,IL2
        CHGT(IL,L)  = 0.
        CHGQ(IL,L)  = 0.
        CHGU(IL,L)  = 0.
        CHGV(IL,L)  = 0.
      ENDDO
      ENDDO
C
C     * STORE INCOMING {ALMX,ALMC} INTO {ALMXIN,ALMCIN} WHICH ARE
C     * USED AS INPUT (NON-TILED, AS PER MICROPHYSICS) TO NLCLMX8.
C     * THEN ZERO OUT {ALMX,ALMC} SINCE THESE WILL BE AGGREGATED OVER TILES.
C
      DO L=1,ILEV
      DO IL=IL1,IL2
        ALMXIN(IL,L) = ALMX(IL,L)
        ALMCIN(IL,L) = ALMC(IL,L)
        ALMX  (IL,L) = 0.
        ALMC  (IL,L) = 0.
        WSUB  (IL,L) = 0.
        RKMROL(IL,L) = 0.
        RKHROL(IL,L) = 0.
        RKQROL(IL,L) = 0.        
      ENDDO
      ENDDO
C
C     * INITIALIZE CND/DEP ARRAYS.
C
      IF ( ISAVLS.NE.0 ) THEN
        DO L=1,ILEV
        DO IL=IL1,IL2 
          CNDROL(IL,L) = 0.
          DEPROL(IL,L) = 0.
        ENDDO
        ENDDO
      ENDIF
C
C     * NOTE THAT THE FOLLOWING IS ALL CONTAINED IN THE TILING LOOP,
C     * DUE TO DEPENDANCE ON TILED INPUT FIELDS!
C
      DO 800 NT=1,NTILEND
C
C     * SET INPUT FIELDS ACCORDING TO VALUE OF "ITILVRT".
C
      FARETOT=0.
      DO IL=IL1,IL2
        IF(ITILVRT.EQ.0)                          THEN
          HFS (IL) = HFSAGG (IL)
          QFS (IL) = QFSAGG (IL)
          CDM (IL) = CDMAGG (IL)
          GT  (IL) = GTAGG  (IL)
          QG  (IL) = QGAGG  (IL)
          FARE(IL) = 1.
          FARETOT  = FARETOT + FARE(IL)
        ELSE
          FARE(IL) = FAREROT(IL,NT)
          IF(FARE(IL).EQ.0.) THEN
C
C           * USE AGGREGATE VALUES AS DUMMY INPUT TO AVOID REQUIRING
C           * "IF" CONDITIONS ON "FARE" FOR EACH SUBSEQUENT DO-LOOP.
C           * ANY ZERO-TILE VALUES SUBSEQUENTLY COMPUTED WILL BE IGNORED.
C
            HFS (IL) = HFSAGG (IL)
            QFS (IL) = QFSAGG (IL)
            CDM (IL) = CDMAGG (IL)
            GT  (IL) = GTAGG  (IL)
            QG  (IL) = QGAGG  (IL)
          ELSE
            HFS (IL) = HFSROT (IL,NT)
            QFS (IL) = QFSROT (IL,NT)
            CDM (IL) = CDMROT (IL,NT)
            GT  (IL) = GTROT  (IL,NT)
            QG  (IL) = QGROT  (IL,NT)
          ENDIF
          FARETOT  = FARETOT + FARE(IL)
        ENDIF
      ENDDO
C
C     * SKIP CALCULATIONS IF NO TILED POINTS IN THIS PASS.
C
      IF(FARETOT.EQ.0.) GO TO 800     
C
C     * RESET WORK FIELDS TO INPUT FIELDS.
C
      DO N=1,NTRAC
      DO L=1,LEV
      DO IL=IL1,IL2
        XROW(IL,L,N)= XIO(IL,L,N)
      ENDDO
      ENDDO
      ENDDO

      DO L=1,ILEV
      DO IL=IL1,IL2
        TH  (IL,L)  = TIO(IL,L)
        Q   (IL,L)  = QIO(IL,L)
        U   (IL,L)  = UIO(IL,L)
        V   (IL,L)  = VIO(IL,L)
      ENDDO
      ENDDO
C
C     * PERFORM NON-LOCAL MIXING.
C
      CALL NLCLMX8(XROW,Q,TH,P,Z,ZF,DSHJ,SHXKJ,
     1             QFS,CQFXROW,HFS,CHFXROW,CDM,CVSG,SIGMA,
     2             PRESSG,CRH,ITRPHS,ZTMST,GRAV,RGAS,TVFA,TFREZ,
     3             TICE,HV,HS,IIWC,ILWC,MSG,IL1,IL2,ILG,ILEV,
     4             LEV,NTRAC,QT,HMN,DSR,CHI,RRL,CPM,DSMIX,HINT,
     5             QINT,WRKR,WRKRL,XINT,SXFLX,
     6             BCR,PBLTX,ZER,QTG,HMNG,
     7             QCW,ZCLF,ZCRAUT,SSH,ALMXIN,ALMCIN,VMODL,
     8             ALMIX,DHLDZ,DRWDZ,WORK(1,1),WORK(1,2),
     9             IPBL,CNDROLX,DEPROLX,XLWROL,XICROL,ISAVLS)
C
C-----------------------------------------------------------------------      
C     * CALCULATE THE VERTICAL DIFFUSION COEFFICIENTS, (RKH,RKM).             
C       --------------------------------------------------------              
C
      TAUADJ=ZTMST
      DO 110 I=IL1,IL2
         CDVLM(I)=CDM(I)*VMODL(I)                                              
         HEAT(I)=TAUADJ*GRAV*SHXKJ(I,ILEV)/(RGAS*TH(I,ILEV))                   
         FACMOM=1./(1.+HEAT(I)*CDVLM(I)/DSHJ(I,ILEV))
         UTMP=U(I,ILEV)*FACMOM
         VTMP=V(I,ILEV)*FACMOM
         DVDS(I,ILEV-1) = SQRT((UTMP-U(I,ILEV-1))**2+           
     1                (VTMP-V(I,ILEV-1))**2)/(SGJ(I,ILEV)-SGJ(I,ILEV-1)) 
  110 CONTINUE                                 
C                                                                             
C     * EVALUATE DTTDS=D(THETA)/D(SIGMA) AT INTERFACE OF THERMODYNAMIC        
C     * LAYERS.                                                               
C                                                                             
      DO 120 L=2,ILEV                                                         
      DO 120 I=IL1,IL2                                                        
         ALWC=XROW(I,L+1,ILWC)+XROW(I,L+1,IIWC)
         THETVP=TH(I,L)*(1.+TVFA*Q(I,L)/(1.-Q(I,L))-(1.+TVFA)*ALWC)
     1         /SHXKJ(I,L)                           
         ALWC=XROW(I,L,ILWC)+XROW(I,L,IIWC)
         THETVM=TH(I,L-1)*(1.+TVFA*Q(I,L-1)/(1.-Q(I,L))-(1.+TVFA)*ALWC)
     1         /SHXKJ(I,L-1)                     
         DTTDS(I,L)=(THETVP-THETVM)                                            
     1             /(SHJ(I,L)-SHJ(I,L-1))   
  120 CONTINUE                                                                
C                                                                             
C     * DVDS = MOD(DV/DSIGMA) AT THE **BOTTOM** INTERFACE OF MOMENTUM LAYERS. 
C                                                                             
      DO 130 L=1,ILEV-2
      DO 130 I=IL1,IL2                                                        
         DVDS(I,L)=SQRT((U(I,L+1)-U(I,L))**2+                                 
     1                  (V(I,L+1)-V(I,L))**2)/(SGJ(I,L+1)-SGJ(I,L))           
  130 CONTINUE                                                                
C                                                                             
      DO 140 I=IL1,IL2                                                        
         PBLTX(I) = REAL(IPBL(I))
         IPBLC(I)=IPBL(I)
         IPBL(I)=1                                                            
         DVDS(I,ILEV)=DVDS(I,ILEVM)                                           
  140 CONTINUE                                                                
C                                                                             
C     * INTERPOLATE DVDS TO ** TOP ** INTERFACE OF THERMODYNAMIC LAYERS.      
C     * RMS WIND SHEAR, DVMINS=(R*T/G)*(DV/DZ)=7.9 FOR DV/DZ=1.MS-1KM-1.      
C     * DETERMINE PLANETARY BOUNDARY LAYER TOP AS LEVEL INDEX ABOVE           
C     * WHICH THE RICHARDSON NUMBER EXCEEDS THE CRITICAL VALUE OF 1.00.       
C                                                                             
      DO 150 L=ILEV,2,-1                                                      
      DO 150 I=IL1,IL2                                                        
         DVDS(I,L)=MAX( DVMINS/SHBJ(I,L-1),                                    
     1                   (DVDS(I,L-1)*(SGBJ(I,  L)-SHBJ(I,L-1))                
     2                   +DVDS(I,L  )*(SHBJ(I,L-1)-SGBJ(I,L-1)))               
     3                               / DSGJ(I,  L) )                           
         RI  (I,L)=-RGAS*SHTXKJ(I,L)*DTTDS(I,L)                               
     1              /(SHBJ(I,L-1)*DVDS(I,L)**2)                               
         IF(RI(I,L).LT.1.00 .AND. IPBL(I).EQ.1)                 THEN          
            PBLX=REAL(L-1)
            PBLTX(I) = MIN(PBLTX(I), PBLX)
         ELSE                                                                 
            IPBL(I)=0                                                         
         ENDIF                                                                
  150 CONTINUE                                                                
C                                    
      FACT=VKC*RGAS*273.                        
      RINEUT=1.
      RIINF=0.25
      DO 260 L=2,ILEV                                         
C                                                                              
C        * CALCULATE THE DIFFUSION COEFFICIENTS (RKH,RKM)                    
C        * AT THE ** TOP ** INTERFACE OF THERMODYNAMIC LAYERS.             
C                                                                             
C        * HEAT: FINITE STABILITY CUTOFF.                                     
C        * LOW XINGL USED SINCE PREVIOUSLY MIXED "INSTANTANEOUSLY".         
C
C        * MOMENTUM: FINITE STABILITY CUTOFF.
C        * HIGH XINGL USED SINCE NOT PREVIOUSLY MIXED "INSTANTANEOUSLY".     
C
C        * UNSTABLE CASE (RI <= 0.) - STABLE CASE (RI > 0.) .                  
C
         DO 220 I=IL1,IL2                                                    
            DVDZ(I)=DVDS(I,L)*GRAV*SHBJ(I,L-1) / (RGAS*TF(I,L))             
  220    CONTINUE
C
         DO 250 I=IL1,IL2                       
            FACT0=FACT*LOG(SHBJ(I,L-1))                             
            SARI=SQRT(ABS(RI(I,L)))
            LPBL=NINT(PBLTX(I))
C
C           * LOWER CUTOFFS.
C
            ALL=0.5*VKC*ZF(I,LPBL)
            XIMINT =MAX(75.*ALL/(75.+ALL),XLMIN)
            ALL=0.5*VKC*Z (I,L)
            XIMIN  =MAX(75.*ALL/(75.+ALL),XLMIN)
C
C           * EFFECTIVE MIXING LENGTHS AS FUNCTION OF MIXING
C           * LENGTHS FOR UP- AND DOWNWARD MIXING.
C
            ALU=VKC*Z(I,L)
            ALD=VKC*(ZF(I,LPBL)-Z(I,L))
            IF ( ALD.GT.0. ) THEN
              XINGLH=MAX(ALU*ALD/(ALU+ALD),XIMIN)
              XINGLM=MAX(ALU*ALD/(ALU+ALD),XIMIN)
              ALMCX(I,L)=XINGLH
            ELSE
              FAC=(P(I,L)/PF(I,LPBL))**2.
              XINGLH=10.+(XIMINT-10.)*FAC
              XINGLM=10.+(XIMINT-10.)*FAC
              IF ( IPAM.EQ.1 ) THEN
                ALMCX(I,L)=FLND(I)*600.+(1.-FLND(I))*100.
                ALMCX(I,L)=MIN(ALMCX(I,L)*FAC,VKC*(Z(I,L)-ZF(I,LPBL)))
              ELSE
                ALMCX(I,L)=MIN(600.*FAC,VKC*(Z(I,L)-ZF(I,LPBL)))
              ENDIF
            END IF
            ALMXX(I,L)=XINGLH
            ALMXT(I,L)=MAX(ALMXX(I,L),ALMCX(I,L),10.)
            XINGLH=(1.-ZCLFC(I,L))*XINGLH
     1              +ZCLFC(I,L)*ALMXT(I,L)
            XINGLM=XINGLH
C
C           DIFFUSION COEFFICIENTS.
C
            ELH   = XINGLH**2         
            EPSSH = 0.
            FACTH = 0.5*RI(I,L)*EPSSH*BEEM
            IF(FACTH.GT.1.)    THEN
              XH  = 0.
            ELSE
              XH  = (1.-FACTH)**2                   
            ENDIF
            ELM   = XINGLM**2         
            EPSSM = 0.
            FACTM = 0.5*RI(I,L)*EPSSM*BEEM
            IF(FACTM.GT.1.)    THEN
              XM  = 0.
            ELSE
              XM  = (1.-FACTM)**2                   
            ENDIF
            IF(RI(I,L).LT.0.)                                       THEN
              RKH(I,L) = ELH*DVDZ(I)*(1.-ALFAH*GAMRH*BEEM*RI(I,L) /       
     1                   (GAMRH+ALFAH*BEEM*SARI))
              RKM(I,L) = ELM*DVDZ(I)*(1.-      GAMRM*BEEM*RI(I,L) /     
     1                   (GAMRM+      BEEM*SARI))
            ELSE
              RKH(I,L) = ELH*DVDZ(I)*XH/(1.+(1.-EPSSH)*BEEM*RI(I,L))      
              RKM(I,L) = ELM*DVDZ(I)*XM/(1.+(1.-EPSSM)*BEEM*RI(I,L))   
            ENDIF 
C
C           * PRANDTL NUMBER SCALING FOR EDDY MOMENTUM DIFFUSIVITY,
C           * BASED ON APPROACH SUGGESTED BY SCHUMANN AND GERZ (1995).
C 
            ATMP=0.
            IF ( RI(I,L).GT.0. ) THEN
              ATMP=RINEUT*EXP(-RI(I,L)/(RINEUT*RIINF))
            ENDIF
            PRANDTL=ATMP+RI(I,L)/RIINF
            PRANDTL=MIN(MAX(PRANDTL,PRANDTL_MIN),PRANDTL_MAX)
            RKM(I,L)=RKM(I,L)*PRANDTL 
  250    CONTINUE
  260 CONTINUE                             
C                                                                             
      DO 270 I=IL1,IL2                                                        
         RKM(I,1)=RKM(I,2)                
         RKH(I,1)=RKH(I,2)                           
         ALMCX(I,1)=XLMIN
         ALMXX(I,1)=XLMIN
  270 CONTINUE          
C                                                       
C     * DEFINE RKQ ALSO AT ** TOP ** INTERFACE OF THERMODYNAMIC LAYERS.       
C                                                                             
      DO 315 L=1,ILEV                                                         
      DO 315 I=IL1,IL2                                             
         RKQ(I,L)=RKH(I,L)                                                
  315 CONTINUE                  
C
C     * DIAGNOSE SUBGRID-SCALE COMPONENT OF THE VERTICAL VELOCITY
C     * (STANDARD DEVIATION, GHAN ET AL., 1997). ACCORDING TO PENG
C     * ET AL. (2005), THE REPRESENTATIVE VERTICAL VELOCITY FOR
C     * AEROSOL ACTIVATION IS OBTAINED BY APPLYING A SCALING FACTOR 
C     * TO THE VELOCITY STANDARD DEVIATION (SCALF).
C
      ATUNE=FLND*1.+(1.-FLND)*0.05
      SCALF=0.7
      FACTS=SCALF*SQRT(2.*PI)
      DO L=1,ILEVM
        WSUBX(IL1:IL2,L)=FACTS*ATUNE(IL1:IL2)*RKH(IL1:IL2,L)
     1                  /(ZF(IL1:IL2,L)-ZF(IL1:IL2,L+1))
      ENDDO
      WSUBX(IL1:IL2,ILEV)=FACTS*ATUNE(IL1:IL2)*RKH(IL1:IL2,ILEV)
     1                   /ZF(IL1:IL2,ILEV)
      WSUBX(IL1:IL2,:)=MIN(MAX(0.01,WSUBX(IL1:IL2,:)),.5)
C
C     * INTERPOLATE RKM TO ** BOTTOM ** INTERFACE OF MOMENTUM LAYERS.
C     * THIS IS REQUIRED FOR USE IN THE ROUTINE "ABCVDM6" WHICH WAS
C     * ORIGINALLY WRITTEN FOR GCMI STAGGERED LAYER SCHEME AND WHICH
C     * REQUIRES RKM TO BE DEFINED AT THE BASE OF THE MOMENTUM LAYER.
C
      DO 321 L=1,ILEVM                                                    
         DO 320 I=IL1,IL2                                                 
            RKM(I,L)=(RKM(I,L  )*(SHTJ  (I,L+1)-SGBJ  (I,L))               
     1               +RKM(I,L+1)*(SGBJ  (I,  L)-SHTJ  (I,L)))    
     2                          / DSHJ  (I,  L)                 
  320    CONTINUE                                                         
  321 CONTINUE                                                            
C                                                                          
C     * SET A LOWER BOUND TO RK'S.                                         
C
      RKMN=0.1
      RKHMN=0.1
      RKQMN=0.1
      RKXMIN(:)=0.1
C                                                                         
      DO 330 L=1,ILEV                                                     
      DO 330 I=IL1,IL2                                                     
         RKM(I,L)=MAX(RKMN ,RKM(I,L))
         RKH(I,L)=MAX(RKHMN,RKH(I,L))
         RKQ(I,L)=MAX(RKQMN,RKQ(I,L))
! Save the eddy coefficients for output diagnostics
         RKMROL(I,L) = RKMROL(I,L) + RKM(I,L)*FARE(I)
         RKHROL(I,L) = RKHROL(I,L) + RKH(I,L)*FARE(I)
         RKQROL(I,L) = RKQROL(I,L) + RKQ(I,L)*FARE(I)
  330 CONTINUE
C
C     * DEFINE THE MINIUMUM ("BACKGROUND") DIFFUSIVITIES FOR ALL THE
C     * TRACERS.
C
      DO 350 N=1,NTRAC
        DO L=1,ILEV
        DO I=IL1,IL2
           RKX(I,L,N)=MAX(RKXMIN(N),RKQ(I,L))
        ENDDO
        ENDDO
 350  CONTINUE
C-----------------------------------------------------------------------
C     * CALCULATE THE TENDENCIES AND REFINE EVALUATION OF SFC FLUXES.
C       -------------------------------------------------------------

C     * CALCULATE THE COEFFICIENT MATRIX FOR MOMENTUM DIFFUSION.
C     * GET TENDENCIES FROM MOMENTUM VERTICAL DIFFUSION.
C     * EVALUATE UFS AND VFS; SINCE THESE ARE TO BE SCALED BY THE SAVING
C     * INTERVAL, THEIR ACCUMULATION DONE INSIDE IMPLVDH MUST BE RE-
C     * DONE OUTSIDE.

      TODT=ZTMST
      CALL ABCVDM6 (A,B,C,CL,CDVLM,GRAV,IL1,IL2,ILG,ILEV,
     1              RGAS,RKM,SGJ,SGBJ,SHJ,TSGB,TODT,.TRUE.)

      CALL IMPLVD7 (A,B,C,CL,U,CL,IL1,IL2,ILG,ILEV,TODT,
     1              TEND,DSGJ,RAUS,WORK,VINT)
C
C     * COMBINE THE TENDENCIES DUE TO MOMENTUM DIFFUSION.
C
      DO 498 I=IL1,IL2
        IF(FARE(I).GT.0.)                                      THEN
          FACT       = VINT(I)*RAUS(I)
          UFS    (I) = UFS    (I)  + FARE(I)*FACT*SAVERAD
          UFSROL (I) = UFSROL (I)  + FARE(I)*FACT*SAVEBEG
          UFSINST(I) = UFSINST (I) + FARE(I)*FACT
          IF(ITILVRT.EQ.1)                                   THEN
            IF(NT.EQ.IOWAT)                                THEN
              UFSOROL(I) = UFSOROL(I) + (UCVTS(I)+FACT)*SAVEBEG
            ELSE IF(NT.EQ.IOSIC)                           THEN
              UFSIROL(I) = UFSIROL(I) + (UCVTS(I)+FACT)*SAVEBEG      
            ENDIF
          ENDIF
        ENDIF
  498 CONTINUE

      DO 500 L=1,ILEV
      DO 500 I=IL1,IL2
         UTG(I,L) = UTG(I,L) + TEND(I,L)*FARE(I)
         U  (I,L) = U  (I,L) + TODT*TEND(I,L)
  500 CONTINUE

      CALL IMPLVD7 (A,B,C,CL,V,CL,IL1,IL2,ILG,ILEV,TODT,
     1              TEND,DSGJ,RAUS,WORK,VINT)

      DO 505 I=IL1,IL2
        IF(FARE(I).GT.0.)                                      THEN
          FACT       = VINT(I)*RAUS(I)
          VFS    (I) = VFS    (I)  + FARE(I)*FACT*SAVERAD
          VFSROL (I) = VFSROL (I)  + FARE(I)*FACT*SAVEBEG
          VFSINST(I) = VFSINST (I) + FARE(I)*FACT
          IF(ITILVRT.EQ.1)                                   THEN
            IF(NT.EQ.IOWAT)                                THEN
              VFSOROL(I) = VFSOROL(I) + (VCVTS(I)+FACT)*SAVEBEG
            ELSE IF(NT.EQ.IOSIC)                           THEN
              VFSIROL(I) = VFSIROL(I) + (VCVTS(I)+FACT)*SAVEBEG      
            ENDIF
          ENDIF
        ENDIF
  505 CONTINUE

      DO 510 L=1,ILEV
      DO 510 I=IL1,IL2
         VTG(I,L) = VTG(I,L) + TEND(I,L)*FARE(I)
         V  (I,L) = V  (I,L) + TODT*TEND(I,L)
  510 CONTINUE 
C--------------------------------------------------------------
C     * LIQUID WATER STATIC ENERGY AND TOTAL WATER PROFILES. 
C
      DO L=1,ILEV
      DO IL=IL1,IL2 
        EST=(1.-CHI(IL,L))*ESW(TH(IL,L))+CHI(IL,L)*ESI(TH(IL,L))
        IF ( EST < P(IL,L) ) THEN 
          QCWX=(XROW(IL,L+1,ILWC)+XROW(IL,L+1,IIWC))*DSR(IL,L)
          QT(IL,L)=Q(IL,L)*DSR(IL,L)+QCWX
          HMN(IL,L)=CPM(IL,L)*TH(IL,L)+GRAV*Z(IL,L)-RRL(IL,L)*QCWX
        ELSE
          QT(IL,L)=Q(IL,L)*DSR(IL,L)
          HMN(IL,L)=CPM(IL,L)*TH(IL,L)+GRAV*Z(IL,L)
        ENDIF
        IF (QT(IL,L).LE.0.) CALL WRN('VRTDF22',-1)
        QT(IL,L)=MAX(QT(IL,L),0.)
      ENDDO
      ENDDO
C
C     * LIQUID WATER STATIC ENERGY AND TOTAL WATER AT SURFACE.
C
      L=ILEV
      DO IL=IL1,IL2 
        QTG(IL)=QG(IL)*DSR(IL,L)
        HMNG(IL)=CPM(IL,L)*GT(IL)
      ENDDO
C
C     * CALCULATE THE COEFFICIENT MATRIX FOR TRACER DIFFUSION.            
C     * GET TENDENCIES FROM VERTICAL TRACER DIFFUSION.                    
C     * EVALUATE TRACSFS.                                                 
C
      DO 532 N=1,NTRAC  
       IF (ITRPHS(N).NE.0 .OR. N.EQ.ILWC .OR. N.EQ.IIWC)    THEN
                                                                              
         CALL ABCVDQ6 (                                                    
     1                 A,B,C,CL, UN ,CDVLT,GRAV,                           
     2                 IL1,IL2,ILG,ILEV,ILEV+1,ILEV,                       
     3                 RGAS,RKX(1,1,N),SHTJ,SHJ,DSHJ,                      
     4                 TH(1,ILEV),TF,TODT)                                 
                                                                      
         CALL IMPLVD7(A,B,C,CL,XROW(1,2,N),ZER,IL1,IL2,ILG,ILEV,           
     1                TODT,TEND,DSHJ,RAUS,WORK,VINT)                       
C                                                                              
C        * NOW APPLY THE TENDENCY COMPUTED IN IMPLVD7:
C
         DO L=1,ILEV
         DO I=IL1,IL2
            XROW(I,L+1,N)=XROW(I,L+1,N)+TODT*TEND(I,L)
         ENDDO
         ENDDO
       ENDIF
 532  CONTINUE
C
C     * SAVE CLOUD WATER/ICE PROFILES AFTER MIXING AS PASSIVE TRACERS.
C
      IF ( ISAVLS.NE.0 ) THEN
        DO L=1,ILEV
        DO IL=IL1,IL2 
          XLWROL(IL,L)=XROW(IL,L+1,ILWC)
          XICROL(IL,L)=XROW(IL,L+1,IIWC)
        ENDDO
        ENDDO
      ENDIF
C                                                                            
C     * CALCULATE THE COEFFICIENT MATRIX FOR DIFFUSION OF LIQUID
C     * WATER STATIC ENERGY AND TOTAL WATER. GET TENDENCIES FROM 
C     * VERTICAL DIFFUSION.     
C
      CALL ABCVDQ6 (                                                    
     1              A,B,C,CL,ZER ,CDVLH,GRAV,                           
     2              IL1,IL2,ILG,ILEV,ILEV+1,ILEV,                       
     3              RGAS,RKH,SHTJ,SHJ,DSHJ,                      
     4              TH(1,ILEV),TF,TODT)                                 
                                                                      
      CALL IMPLVD7(A,B,C,CL,HMN,HMNG,IL1,IL2,ILG,ILEV,           
     1             TODT,TEND,DSHJ,RAUS,WORK,VINT)                       
C                                                                              
C     * NOW APPLY THE TENDENCY COMPUTED IN IMPLVD7:
C
      DO L=1,ILEV                                                   
      DO IL=IL1,IL2                                                  
        HMN(IL,L)=HMN(IL,L)+TODT*TEND(IL,L)                  
      ENDDO
      ENDDO
      CALL ABCVDQ6 (                                                    
     1              A,B,C,CL,ZER ,CDVLH,GRAV,                           
     2              IL1,IL2,ILG,ILEV,ILEV+1,ILEV,                       
     3              RGAS,RKH,SHTJ,SHJ,DSHJ,                      
     4              TH(1,ILEV),TF,TODT)                                 
                                                                      
      CALL IMPLVD7(A,B,C,CL,QT,QTG,IL1,IL2,ILG,ILEV,           
     1             TODT,TEND,DSHJ,RAUS,WORK,VINT)                       
C                                                                              
C     * NOW APPLY THE TENDENCY COMPUTED IN IMPLVD7:
C
      DO L=1,ILEV
      DO IL=IL1,IL2
        ACOR=1.
        IF(QT(IL,L)+TODT*TEND(IL,L).LT.0.) THEN
        ACORT=-QT(IL,L)/(TEND(IL,L)*TODT)
        ACOR1=MAX(MIN(ACOR,ACORT),0.)
        ACOR=(3.-2.*ACOR1)*ACOR1**3
        ENDIF
        QT(IL,L)=QT(IL,L)+TODT*ACOR*TEND(IL,L)
      ENDDO
      ENDDO
C
C     * UNRAVEL LIQUID WATER STATIC ENERGY AND TOTAL WATER TO OBTAIN
C     * SPECIFIC HUMIDITY, TEMPERATURE, AND CLOUD WATER.
C
      ICVSG=0
      ISUBG=1
      DO L=1,ILEV
        DO IL=IL1,IL2
          IF(L.GT.1)                            THEN
            DZ       =Z   (IL,L-1)-Z  (IL,L)
            DHLDZ(IL)=(HMN(IL,L-1)-HMN(IL,L))/DZ
            DRWDZ(IL)=(QT (IL,L-1)-QT (IL,L))/DZ
            ALMIX(IL)=ALMXT(IL,L)
          ELSE
            DHLDZ(IL)=0.
            DRWDZ(IL)=0.
            ALMIX(IL)=10.
          ENDIF 
C
          EST=(1.-CHI(IL,L))*ESW(TH(IL,L))+CHI(IL,L)*ESI(TH(IL,L))
          IF(L.GT.LEV1 .AND.
     1       EST.LT.P(IL,L) .AND. P(IL,L).GE.10.)    THEN
            X(IL)=1.
          ELSE
            X(IL)=0.
          ENDIF 
        ENDDO
C
        IDUM=0
        ICALL=2       
        CALL STATCLD5(QCW,ZCLF,SIGMA(1,L),ZCRAUT,WORK(1,1),SSH,
     1                CVSG(1,L),QT(1,L),HMN(1,L),CHI(1,L),CPM(1,L),
     2                P(1,L),Z(1,L),RRL(1,L),ZER,X,
     3                ALMIX,DHLDZ,DRWDZ,IDUM,
     4                GRAV,ZTMST,ILEV,ILG,IL1,IL2,ICVSG,ISUBG,
     5                L,ICALL                                )
C
        DO IL=IL1,IL2
         IF(X(IL).GT.0.)                                  THEN 
          QCWI=CHI(IL,L)*QCW(IL)/DSR(IL,L)
          QCWL=(1.-CHI(IL,L))*QCW(IL)/DSR(IL,L)
C
C         * UPDATE TEMPERATURE, SPECIFIC HUMIDITY, AND CLOUD WATER.
C
          TH(IL,L)=(HMN(IL,L)-GRAV*Z(IL,L)+RRL(IL,L)*QCW(IL))/CPM(IL,L)
          Q (IL,L)=(QT(IL,L)-QCW(IL))/DSR(IL,L)
          XROW(IL,L+1,ILWC)=QCWL
          XROW(IL,L+1,IIWC)=QCWI
         ELSE
C
C         * UPDATE TEMPERATURE, SPECIFIC HUMIDITY, AND CLOUD WATER IN DRY
C         * LAYER NEAR TOP OF MODEL DOMAIN ASSUMING THAT THERE IS NO
C         * CONDENSATE.
C
          TH(IL,L)=(HMN(IL,L)-GRAV*Z(IL,L))/CPM(IL,L)
          Q (IL,L)=QT(IL,L)
          XROW(IL,L+1,ILWC)=0.
          XROW(IL,L+1,IIWC)=0.
         ENDIF
        ENDDO
      ENDDO
C
C     * SAVE CONDENSATION/DEPOSITION TENDENCIES.
C
      IF ( ISAVLS.NE.0 ) THEN
        DO 575 L=1,ILEV
        DO 575 IL=IL1,IL2 
          CNDROLX(IL,L)=CNDROLX(IL,L)
     1                +(XROW(IL,L+1,ILWC)-XLWROL(IL,L))/TODT
          DEPROLX(IL,L)=DEPROLX(IL,L)
     1                +(XROW(IL,L+1,IIWC)-XICROL(IL,L))/TODT
  575   CONTINUE
      ENDIF
C
C     * AGGREGATE CHANGES OVER ALL TILES.
C
      DO N=1,NTRAC
      DO L=1,LEV
      DO IL=IL1,IL2
        CHGX(IL,L,N)= CHGX(IL,L,N) + (XROW(IL,L,N)-XIO(IL,L,N))*FARE(IL)
      ENDDO
      ENDDO
      ENDDO

      DO L=1,ILEV
      DO IL=IL1,IL2
        CHGT(IL,L)  = CHGT(IL,L) + (TH(IL,L)-TIO(IL,L))*FARE(IL)
        CHGQ(IL,L)  = CHGQ(IL,L) + (Q (IL,L)-QIO(IL,L))*FARE(IL)
        CHGU(IL,L)  = CHGU(IL,L) + (U (IL,L)-UIO(IL,L))*FARE(IL)
        CHGV(IL,L)  = CHGV(IL,L) + (V (IL,L)-VIO(IL,L))*FARE(IL)
C
        ALMX(IL,L)  = ALMX(IL,L) + ALMXX(IL,L)*FARE(IL)
        ALMC(IL,L)  = ALMC(IL,L) + ALMCX(IL,L)*FARE(IL)
        WSUB(IL,L)  = WSUB(IL,L) + WSUBX(IL,L)*FARE(IL)
        IF ( ISAVLS.NE.0 ) THEN
          CNDROL(IL,L) = CNDROL(IL,L) + CNDROLX(IL,L)*FARE(IL)
          DEPROL(IL,L) = DEPROL(IL,L) + DEPROLX(IL,L)*FARE(IL)
        ENDIF
      ENDDO
      ENDDO
C
C     * CALCULATE THE AGGREGATE PBL VERTICAL INDEX.
C
      DO IL=IL1,IL2
        PBLT(IL) = PBLT(IL) + PBLTX(IL)*FARE(IL)
      ENDDO

 800  CONTINUE   ! end of tiling loop!
C
C     * CALCULATE THE NEW PROGNOSTIC THERMODYNAMIC VALUES.
C
      DO N=1,NTRAC
      DO L=1,LEV
      DO IL=IL1,IL2
        XIO(IL,L,N) = XIO(IL,L,N) + CHGX(IL,L,N)
      ENDDO
      ENDDO
      ENDDO

      DO L=1,ILEV
      DO IL=IL1,IL2
        TIO(IL,L)  = TIO(IL,L) + CHGT(IL,L)
        QIO(IL,L)  = QIO(IL,L) + CHGQ(IL,L)
        UIO(IL,L)  = UIO(IL,L) + CHGU(IL,L)
        VIO(IL,L)  = VIO(IL,L) + CHGV(IL,L)
      ENDDO
      ENDDO
                                                                                  
      RETURN                                                                      
C-----------------------------------------------------------------------          
      END
