      SUBROUTINE NLCLMX8(XROW,QH,TH,P,Z,ZF,DSHJ,SHXKJ,SQFLX,TQFLX,SHFLX,
     1                   THFLX,CDM,CVSG,SIGMA,PRESSG,X,
     2                   ITRPHS,DT,GRAV,RGAS,TVFA,TFREZ,
     2                   TICE,HV,HS,IIWC,ILWC,MSG,IL1,IL2,ILG,ILEV,
     3                   LEV,NTRAC,QT,HMN,DSR,CHI,RRL,CPM,DSMIX,HINT,
     4                   QINT,TVPPBL,TAU,XINT,SXFLX,
     5                   BCR,PBLT,ZER,QTN,HMNN,
     6                   QCW,ZCLF,ZCRAUT,SSH,ALMX,ALMC,VMODL,
     7                   ALMIX,  DHLDZ,   DRWDZ,   HMNMIN,    QTMIN, 
     8                   IMXL,CNDROL,DEPROL,XLWROL,XICROL,ISAVLS       )     
C
C     * MAY 03/2018 - N.MCFARLANE. Added code to prevent statcld crashes.
C     *                            Removed unused {ESW,ESI} statement functions.
C     * FEB 10/2015 - M.LAZARE/      New version 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     * JUL 31/2013 - M.LAZARE.    Previous version NLCLMX7 for gcm17:
C     *                            Cosmetic changes (since QTN.ge.0. in
C     *                            our tests) to:
C     *                            - Print out warning if QTN<0 before
C     *                              calling statcld5, then limit
C     *                              it to be non-negative (2x).
C     *                            - Pass in LVL and ICALL to statcld5
C     *                              to aid in future debugging.
C     * JUN 26/2013 - K.VONSALZEN. NEW VERSION FOR GCM17:
C     *                            - CALLS NEW STATCLD5.
C     * APR 29/2012 - K.VONSALZEN/ PREVIOUS VERSION NLCLMX6 FOR GCM16:
C     *               M.LAZARE.    - CALLS NEW STATCLD4.
C     * JAN 17/2008 - K.VONSALZEN/ PREVIOUS VERSION NLCLMX5 FOR GCM15G/H/I:
C     *               M.LAZARE.    - CALLS NEW STATCLD3.
C     *                            - BUGFIX IN EXPRESSION FOR USTAR TO
C     *                              PROPERLY INCLUDE WIND SPEED. 
C     *                            - MODIFIED CALCULATION OF ADJUSTMENT
C     *                              TIMESCALE IN MIXED LAYER (TAU).
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 DT(=ADELT FROM PHYSICS)
C     *                              TO NEW STATCLD3 IN PLACE OF
C     *                              HARD-CODED 2.*DELT.
C     *                            - GUSTINESS EFFECT NOW INCLUDED 
C     *                              SINCE IS NOW ALREADY CONTAINED IN ZSPD
C     *                              FROM PHYSICS.
C     * JAN 13/2007 - K.VONSALZEN. PREVIOUS VERSION NLCLMX4 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.
C     * JUN 19/2006 - M.LAZARE.    - CALLS NEW VERSION STATCLD2.
C     *                            - COSMETIC: 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     * MAY 06/06 - K.VONSALZEN/ PREVIOUS VERSION NLCLMX3 FOR GCM15E:
C     *             M.LAZARE.    - USE NEW COMMON SUBROUTINE STATCLD FOR 
C     *                            UNRAVELLING.
C     *                          - SINCE SIMILAR MIXING NOW ALSO DONE IN NEW
C     *                            SUBROUTINE VRTDF14 (IE ALSO CALLS
C     *                            STATCLD), PASS IN REQUIRED INPUT FIELDS
C     *                            INSTEAD OF RE-DERIVING THEM HERE.
C     *                          - REORDERING OF TERMS AND LIMITATION OF
C     *                            FLUXES AT TOP OF MIXED LAYER TO IMPROVE
C     *                            NUMERICAL ROBUSTNESS.
C     *                          - REVISIONS TO ALMIX AND CONVECTIVE
C     *                            CONTRIBUTION TO VARIANCE, CONSISTENT WITH
C     *                            COND3.
C     * DEC 14/05 - K.VONSALZEN/ PREVIOUS VERSION NLCLMX2 FOR GCM15D.
C     *             M.LAZARE.    
C     *
C     * PERFORMS NON-LOCAL MIXING IN PLANETARY BOUNDARY LAYER BASED
C     * ON SIMPLE ADJUSTMENT SCHEME. THE MIXING IS PERFORMED FOR LIQUID 
C     * WATER STATIC ENERGY AND TOTAL WATER. TEMPERATURE AND SPECIFIC 
C     * HUMIDITY ARE OBTAINED FROM UNRAVELLING PROCEDURE. A 
C     * SEMI-IMPLICIT METHOD IS USED.
C
      IMPLICIT REAL   (A-H,O-Z),
     +INTEGER (I-N)
C 
      REAL   XROW(ILG,LEV,NTRAC)
      REAL   QH   (ILG,ILEV),  TH   (ILG,ILEV),  DSHJ (ILG,ILEV),
     1       Z    (ILG,ILEV),  P    (ILG,ILEV),  ZF   (ILG,ILEV),
     2       SHXKJ(ILG,ILEV),  CVSG (ILG,ILEV),  SIGMA(ILG,ILEV),
     3       ALMX (ILG,ILEV),  ALMC(ILG,ILEV)
      REAL   SQFLX  (ILG),     TQFLX  (ILG), 
     1       SHFLX  (ILG),     THFLX  (ILG),
     2       CDM    (ILG),     PRESSG (ILG),     X      (ILG),
     3       PBLT   (ILG),     ZER    (ILG),     QTN    (ILG),
     4       HMNN   (ILG),     QCW    (ILG),     ZCLF   (ILG), 
     5       ZCRAUT (ILG),     SSH    (ILG),     ALMIX  (ILG),  
     6       DHLDZ  (ILG),     DRWDZ  (ILG),     HMNMIN (ILG),
     7       QTMIN  (ILG),   VMODL  (ILG)
      REAL   DSR  (ILG,ILEV),  CHI  (ILG,ILEV),  CPM   (ILG,ILEV),
     1       RRL  (ILG,ILEV),  QT   (ILG,ILEV),  HMN   (ILG,ILEV)
      REAL  CNDROL(ILG,ILEV), DEPROL(ILG,ILEV),  XLWROL(ILG,ILEV),
     1      XICROL(ILG,ILEV)
      REAL   DSMIX  (ILG),     HINT   (ILG),     QINT   (ILG),
     1       TVPPBL (ILG),     TAU    (ILG),     BCR    (ILG)
      REAL   XINT(ILG,NTRAC),SXFLX(ILG,NTRAC)
      REAL   ACOR(ILG)

      INTEGER IMXL(ILG)
      INTEGER ITRPHS(NTRAC)
C
C     * INTERNAL WORK ARRAY.
C 
      REAL  QCWVAR (ILG)
C
      PARAMETER ( QCWCRIT=0. )
      PARAMETER ( TAUPRP=1.0 )
      PARAMETER ( YSEC=0.999 )
C
      COMMON /EPS  / A,B,EPS1,EPS2                                     
      COMMON /ESTWI/ RW1,RW2,RW3,RI1,RI2,RI3
C
      DATA ZERO,ONE /0., 1./
      DATA PZ3  /0.03/
C=======================================================================
C     * CONSTANTS.
C
      EXPM1=EXP(-1.)
      CSIGMA=0.2
      ICVSG=0
      ISUBG=1
      YFR=0.50
      ZSECFRL=1.E-7
      ZEPCLC=1.E-3
C
C     * INITIAL PROFILES OF THERMODYNAMIC PROPERTIES:
C     * TOTAL WATER MIXING RATIO AND LIQUID WATER STATIC ENERGY,
C  
      DO 50 L=MSG+1,ILEV
      DO 50 IL=IL1,IL2 
        QCWX=(XROW(IL,L+1,ILWC)+XROW(IL,L+1,IIWC))*DSR(IL,L)
        QV=QH(IL,L)*DSR(IL,L)
        QT(IL,L)=QV+QCWX
        DST=CPM(IL,L)*TH(IL,L)+GRAV*Z(IL,L)
        HMN(IL,L)=DST-RRL(IL,L)*QCWX
 50   CONTINUE  
C
C=======================================================================
C     * DETERMINE DEPTH OF MIXED LAYER FROM ADJUSTMENT SCHEME FOR
C     * MIXED LAYER WITH FLUXES ACROSS TOP AND BOTTOM LAYERS.
C
      L=ILEV
      DO 100 IL=IL1,IL2
        TVPOT=TH(IL,L)*(1.+TVFA*QH(IL,L))/SHXKJ(IL,L)
        TVPPBL(IL)=TVPOT*ZF(IL,L)
        HINT(IL)=HMN(IL,L)*DSHJ(IL,L)
        QINT(IL)=QT (IL,L)*DSHJ(IL,L)
        DSMIX(IL)=DSHJ(IL,L)
        IMXL(IL)=L
        BCR(IL)=9.E+20
        HMNMIN(IL)=HMN(IL,L)
        QTMIN(IL)=QT(IL,L)
 100  CONTINUE
C
C     * PERFORM SEMI-IMPLICIT TIME INTEGRATION TO DETERMINE
C     * THERMODYNAMIC PROPERTIES AT THE TOP OF THE VIRTUAL MIXED 
C     * LAYER. USE BUOYANCY IN THAT LAYER (RELATIVE TO LAYER ABOVE)
C     * TO DETERMINE DEPTH OF MIXED LAYER.
C
      DO 200 L=ILEV-1,MSG+2,-1
       DO 125 IL=IL1,IL2
        IF ( IMXL(IL).EQ.(L+1) ) THEN
C
C         * ADJUSTMENT TIME SCALE IN MIXED LAYER.
C
          RHO=100.*P(IL,ILEV)/(RGAS*TH(IL,ILEV))
          DRDQV=1.
          RFX=SQFLX(IL)*DRDQV/RHO
          TPOT=TH(IL,ILEV)/SHXKJ(IL,ILEV)
          VPTFLX=(SHFLX(IL)/(RHO*CPM(IL,ILEV)))*(1.+TVFA*QH(IL,ILEV))
     1           +TVFA*RFX*TPOT
          DPBL=ZF(IL,L+1)
          WSTAR=MAX(GRAV*VPTFLX*DPBL/(TVPPBL(IL)/DPBL),ZERO)**(1./3.)
          USTAR=SQRT(CDM(IL))*VMODL(IL)
          TAU(IL)=TAUPRP*DPBL/MAX(SQRT(0.6*WSTAR**2+USTAR**2),PZ3)
C
C         * LIQ. WATER STATIC ENERGY AND TOTAL WATER AT TOP OF MIXED LAYER.
C
          FAC=GRAV/(PRESSG(IL)*DSMIX(IL))
          WEIGH=MIN(MAX(TAU(IL)/(TAU(IL)+DT),ZERO),ONE)
          HMNMIN(IL)=MIN(HMNMIN(IL),HMN(IL,L+1))
          THFLXT=MIN(THFLX(IL),SHFLX(IL)+YSEC*HMNMIN(IL)/(DT*FAC))
          QTMIN(IL)=MIN(QTMIN(IL),QT(IL,L+1))
          TQFLXT=MIN(TQFLX(IL),SQFLX(IL)+YSEC*QTMIN(IL)/(DT*FAC))
          HMNN(IL)=WEIGH*HMN(IL,L+1)+(1.-WEIGH)*HINT(IL)/DSMIX(IL)
     1            +DT*FAC*(SHFLX(IL)-THFLXT)
          QTN (IL)=WEIGH*QT (IL,L+1)+(1.-WEIGH)*QINT(IL)/DSMIX(IL)
     1            +DT*FAC*(SQFLX(IL)-TQFLXT)
          IF ( QTN(IL).LE.0. ) CALL WRN('NLCLMX8',-1)
c       QT(IL,L)=MAX(QT(IL,L),0.) 
C
          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
C
          ALMIX(IL)=MAX(ALMX(IL,L),ALMC(IL,L),10.)
          X(IL)=1.
        ELSE
          X(IL)=0. 
        ENDIF
  125  CONTINUE 
C
       IDUM=0
       ICALL=3
       CALL STATCLD5(QCW,ZCLF,SIGMA(1,L),ZCRAUT,QCWVAR,SSH,
     1               CVSG(1,L),QTN,HMNN,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,DT,ILEV,ILG,IL1,IL2,ICVSG,ISUBG,
     5               L,ICALL                                 )
C
       DO 150 IL=IL1,IL2
         IF ( IMXL(IL).EQ.(L+1) ) THEN
C
C         * CLOUD LIQUID AND ICE WATER CONTENTS FROM DIAGNOSTIC
C         * RELATIONSHIP.
C
          QLWC=(1.-CHI(IL,L))*QCW(IL)
          QIWC=CHI(IL,L)*QCW(IL)
C
C         * ADJUSTED WATER VAPOUR MIXING RATIO AND TEMPERATURE.
C
          QVN=QTN(IL)-QCW(IL)
          TN =(HV*QLWC+HS*QIWC-GRAV*Z(IL,L)+HMNN(IL))/CPM(IL,L)
C
C         * BUOYANCY FOR AIR AT TOP OF VIRTUAL MIXED LAYER IN AMBIENT
C         * AIR AT THAT LEVEL.
C
          QCWE=(XROW(IL,L+1,ILWC)+XROW(IL,L+1,IIWC))*DSR(IL,L)
          QVE=QH(IL,L)*DSR(IL,L)
          TV =TH(IL,L)*(1.+TVFA*QVE-(1.+TVFA)*QCWE)
          TPV=TN      *(1.+TVFA*QVN-(1.+TVFA)*QCW(IL) )
          BUOY=GRAV*( TPV-TV )/TV
C
C         * TEST FOR BUOYANCY.
C
          IF ( BUOY.GE.0. .AND. 
     1        .NOT.(BUOY.GT.BCR(IL).AND.QCW(IL).GT.QCWCRIT)   ) THEN
            IMXL(IL)=L
C
C         * VERTICALLY INTEGRATED LIQ. WATER ENERGY AND TOTAL WATER,
C         * AND DEPTH OF MIXED LAYER.
C
            TVPOT=TH(IL,L)*(1.+TVFA*QH(IL,L))/SHXKJ(IL,L)
            TVPPBL(IL)=TVPPBL(IL)+TVPOT*(ZF(IL,L)-ZF(IL,L+1))
            HINT(IL)=HINT(IL)+HMN(IL,L)*DSHJ(IL,L)
            QINT(IL)=QINT(IL)+QT (IL,L)*DSHJ(IL,L)
            DSMIX(IL)=DSMIX(IL)+DSHJ(IL,L)
            BCR(IL)=BUOY
          ENDIF
        ENDIF
  150  CONTINUE
  200 CONTINUE
C
C     TRACERS.
C
      DO 400 N=1,NTRAC 
        IF ( ITRPHS(N).NE.0 .OR. N.EQ.ILWC .OR. N.EQ.IIWC ) THEN
          DO 420 IL=IL1,IL2 
            XINT (IL,N)=0.
            SXFLX(IL,N)=0.
 420      CONTINUE
          DO 440 L=ILEV,MSG+1,-1 
          DO 440 IL=IL1,IL2 
            IF ( L.GE.IMXL(IL) ) THEN
              XINT(IL,N)=XINT(IL,N)+XROW(IL,L+1,N)*DSHJ(IL,L)
            ENDIF
 440      CONTINUE
          DO 460 L=ILEV,MSG+1,-1 
          DO 460 IL=IL1,IL2 
            IF ( L.GE.IMXL(IL) ) THEN
C
C             * PERFORM SEMI-IMPLICIT TIME INTEGRATION TO DETERMINE
C             * THERMODYNAMIC PROPERTIES.
C
              FAC=GRAV/(PRESSG(IL)*DSMIX(IL))
              WEIGH=MIN(MAX(TAU(IL)/(TAU(IL)+DT),0.),1.)
              XROW(IL,L+1,N)=WEIGH*XROW(IL,L+1,N)+(1.-WEIGH)*XINT(IL,N)
     1                                    /DSMIX(IL)+DT*FAC*SXFLX(IL,N)
            END IF         
 460      CONTINUE
        ENDIF
 400  CONTINUE
C
C     * SAVE CLOUD WATER/ICE PROFILES AFTER MIXING AS PASSIVE TRACERS.
C
      IF ( ISAVLS.NE.0 ) THEN
        DO 470 L=MSG+1,ILEV
        DO 470 IL=IL1,IL2 
          XLWROL(IL,L)=XROW(IL,L+1,ILWC)
          XICROL(IL,L)=XROW(IL,L+1,IIWC)
 470    CONTINUE
      ENDIF
C
      DO 480 IL=IL1,IL2
        ACOR(IL)=1.
 480  CONTINUE

      DO 485 L=ILEV,MSG+1,-1
      DO 485 IL=IL1,IL2
        IF ( L.GE.IMXL(IL) ) THEN
          FAC=GRAV/(PRESSG(IL)*DSMIX(IL))
          WEIGH=MIN(MAX(TAU(IL)/(TAU(IL)+DT),ZERO),ONE)
          DQ1=DT*FAC*SQFLX(IL)
          DQ2=(1.-WEIGH)*(QINT(IL)/DSMIX(IL)-QT(IL,L))
          DQSUM=DQ1+DQ2
          IF(DQSUM+QT(IL,L).LT.0.) THEN
            CALL WRN('NLCLMX8',-3)
            ACORT=-MAX(QT(IL,L),0.)/DQSUM
            ACOR1=MAX(MIN(ACOR(IL),ACORT),0.)
            ACOR(IL)=(3.-2.*ACOR1)*ACOR1**3
          ENDIF
        ENDIF
 485  CONTINUE
C
C=======================================================================
C     * CALCULATE NEW THERMODYNAMIC VALUES IN MIXED LAYER FOR GIVEN
C     * MIXED LAYER DEPTH AND SURFACE FLUXES. EFFECTS OF FLUXES
C     * ACROSS TOP OF MIXED LAYER ARE ACCOUNTED FOR IN PARAMETERIZATIONS
C     * OF CONVECTION AND THEREFORE DO NOT AFFECT THESE CALCULATION.
C
      DO 550 L=ILEV,MSG+1,-1
       DO 525 IL=IL1,IL2
        IF ( L.GE.IMXL(IL) ) THEN
C
C         * LIQ. WATER STATIC ENERGY AND TOTAL WATER.
C
          FAC=GRAV/(PRESSG(IL)*DSMIX(IL))
          WEIGH=MIN(MAX(TAU(IL)/(TAU(IL)+DT),ZERO),ONE)
          HMNN(IL)=WEIGH*HMN(IL,L)+(1.-WEIGH)*HINT(IL)/DSMIX(IL)
     1            +DT*FAC*SHFLX(IL)
          QTNADJ=QT(IL,L)+ACOR(IL)*((1.-WEIGH)*(QINT(IL)/
     1           DSMIX(IL)-QT(IL,L))+DT*FAC*SQFLX(IL))
          QTN (IL)=QTNADJ
          IF ( QTN(IL).LE.0. ) CALL WRN('NLCLMX8',-2)
c         QTN(IL)=MAX(QTN(IL),0.) 
C
          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)=MAX(ALMX(IL,L),ALMC(IL,L),10.)
          X(IL)=1.
        ELSE
          X(IL)=0.
        ENDIF
  525 CONTINUE 
C
      IDUM=0
      ICALL=4
      CALL STATCLD5(QCW,ZCLF,SIGMA(1,L),ZCRAUT,QCWVAR,SSH,
     1               CVSG(1,L),QTN,HMNN,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,DT,ILEV,ILG,IL1,IL2,ICVSG,ISUBG,
     5               L,ICALL                                 )
C
       DO 540 IL=IL1,IL2
        IF ( L.GE.IMXL(IL) ) THEN
C
C         * CLOUD LIQUID AND ICE WATER CONTENTS FROM DIAGNOSTIC
C         * RELATIONSHIP.
C
          QLWC=(1.-CHI(IL,L))*QCW(IL)
          QIWC=CHI(IL,L)*QCW(IL)
          XROW(IL,L+1,ILWC)=QLWC/DSR(IL,L)
          XROW(IL,L+1,IIWC)=QIWC/DSR(IL,L)
C
C         * ADJUSTED WATER VAPOUR MIXING RATIO AND TEMPERATURE.
C
          QVN=QTN(IL)-QCW(IL)
          QH(IL,L)=QVN/DSR(IL,L)
          TH(IL,L)=(HV*QLWC+HS*QIWC-GRAV*Z(IL,L)+HMNN(IL))/CPM(IL,L)
        ENDIF
  540  CONTINUE
  550 CONTINUE
C
C     * SAVE CONDENSATION/DEPOSITION TENDENCIES.
C
      IF ( ISAVLS.NE.0 ) THEN
        DO 600 L=MSG+1,ILEV
        DO 600 IL=IL1,IL2 
          CNDROL(IL,L)=(XROW(IL,L+1,ILWC)-XLWROL(IL,L))/DT
          DEPROL(IL,L)=(XROW(IL,L+1,IIWC)-XICROL(IL,L))/DT
  600   CONTINUE
      ENDIF
C
      RETURN
C-----------------------------------------------------------------------
      END
