      SUBROUTINE VRTDF6O(   VTG,    UTG,    QFS,    HFS,    UFS,    VFS,
     1                       ST,     SQ,     SU,     SV, UFSROL, VFSROL,
     2                     DRAG,     EF,     GC,     GT,     TS,     QG,
     3                      QFX,    TFX,UTENDGW,VTENDGW, PRESSG,      Q,
     4                        U,      V,    TSG,   TSGB,     TF,     TH,
     5                      SGJ,   SGBJ, SHTXKJ,  SHXKJ,   SHTJ,   SHBJ,
     6                      SHJ,   DSGJ,   DSHJ,    TCV,
     7                     XROW,    XFS,    XFX,  ITRAC,  NTRAC,        
     8                  DELT,ILG,IL1,IL2,ILEV,LEV,LEVS,SAVEGG,SAVEBEG,
     9                  PI,RGAS,RGOCP,CPRES,RKMIN,RKQMIN, 
     A                      EPS,   EPSS,    CDH,    CDM,  CDVLH,
     B                    CDVLM,  VMODL, RATMOH,    CDO,    RIB,
     C                       ZO,   DVDS,  DTTDS,     CL,     RI,
     D                     DVDZ,      X,   HEAT,      A,      B,
     E                        C,   TEND,   RAUS,   VINT,   WORK,
     F                      RKM,    RKQ,    RKH,   WIND,  VMODH,
     G                   SCALPH,    BVF,   BVFL,   AMPL,     DZ,
     H                      DGW, HEIGHT,    TYP,  CDVLT,     XX,
     I                      ZER,     UN,TRACSFS,  UGWTS,  VGWTS        )

C     * JAN 02/94 - M.LAZARE. AS IN PREVIOUS VRTDFS6, XCP FOR NEW,
C     *                       OPTIMIZED MODEL VERSION GCM6O:
C     *                       - JLAT,COSL REMOVED FROM CALL (NOT USED 
C     *                         IN CONJUNCTION WITH MULTIPLE LATITUDES.
C     *                       - OLD, UNUSED GWD CODE REMOVED.
C     *                       - "N-TRACER" CODE FROM HIGHER MODEL
C     *                         VERSIONS ADDED (TRAC,TRACT REMOVED 
C     *                         AND XROW,XFS,XFX,NTRAC ADDED). 
C     * JAN 30/89 - M.LAZARE. PREVIOUS VERSION VRTDFS6.
C 
C     * SMALL MIN VALUE OF RKQ APPLIED. 
C     * RIB COMPUTED WITH FIXED THICKNESS SFC LAYER (SLTHICK).
C     * UTG AND VTG NOW MEAN D(U,V)/DT DUE TO VERT.DIFFUS.
C     * 
C     * CALCULATE THE TENDENCIES OF U,V,TH,Q, DUE TO VERTICAL DIFFUSION 
C     * AND POSSIBLY GRAVITY WAVE DRAG. 
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 REAL (A-H,O-Z),
     +INTEGER (I-N)

C     * I/O ARRAYS. 
  
      REAL UTG    (ILG,ILEV),VTG    (ILG,ILEV)
      REAL U      (ILG,ILEV),V      (ILG,ILEV)
      REAL TH     (ILG,ILEV),TF     (ILG,ILEV),Q      (ILG,ILEV)
      REAL TSG    (ILG,ILEV),TSGB   (ILG,ILEV)
      REAL UTENDGW(ILG,ILEV),VTENDGW(ILG,ILEV)

      REAL XROW   (ILG,LEV,NTRAC)

      REAL SGJ    (ILG,ILEV),SGBJ   (ILG,ILEV)
      REAL SHJ    (ILG,ILEV),SHBJ   (ILG,ILEV),SHTJ   (ILG, LEV)
      REAL SHTXKJ (ILG,ILEV),SHXKJ  (ILG,ILEV)
      REAL DSGJ   (ILG,ILEV),DSHJ   (ILG,ILEV)
  
      REAL QFS    (ILG),     HFS    (ILG),     UFS    (ILG) 
      REAL VFS    (ILG),     UFSROL (ILG),     VFSROL (ILG) 
      REAL XFS    (ILG,NTRAC), XFX  (ILG,NTRAC)                         
      REAL SQ     (ILG),     ST     (ILG),     SU     (ILG) 
      REAL SV     (ILG),     TCV    (ILG) 
  
      REAL DRAG   (ILG),     EF     (ILG),     GC     (ILG) 
      REAL GT     (ILG),     TS     (ILG),     QG     (ILG) 
      REAL QFX    (ILG),     TFX    (ILG),     PRESSG (ILG) 
  
  
C     * WORK ARRAYS DEFINED THROUGH ARRAY "WRK" WITH "IV" POINTERS. 
  
      REAL EPS    (ILG),     EPSS   (ILG),     CDH    (ILG) 
      REAL CDM    (ILG),     CDVLH  (ILG),     CDVLM  (ILG) 
      REAL VMODL  (ILG),     RATMOH (ILG),     CDO    (ILG) 
      REAL RIB    (ILG),     ZO     (ILG),     DVDS   (ILG,ILEV)
      REAL DTTDS  (ILG,ILEV),CL     (ILG),     RI     (ILG) 
      REAL DVDZ   (ILG),     X      (ILG),     HEAT   (ILG) 
      REAL A      (ILG,ILEV),B      (ILG,ILEV),C      (ILG,ILEV)
      REAL TEND   (ILG,ILEV),RAUS   (ILG),     VINT   (ILG) 
      REAL WORK   (ILG,ILEV),RKM    (ILG,ILEV),RKQ    (ILG,ILEV)
      REAL RKH    (ILG,ILEV),WIND   (ILG,ILEV),VMODH  (ILG,ILEV)
      REAL SCALPH (ILG,ILEV),BVF    (ILG,ILEV),BVFL   (ILG) 
      REAL AMPL   (ILG,ILEV),DZ     (ILG,ILEV),DGW    (ILG,ILEV)
      REAL HEIGHT (ILG)     ,TYP    (ILG),     CDVLT  (ILG) 
      REAL XX     (ILG,ILEV),ZER    (ILG),     UN     (ILG) 
      REAL TRACSFS(ILG)     ,UGWTS  (ILG),     VGWTS  (ILG) 
  
      LOGICAL IGW,NOGWDRG 
  
      DATA SLTHICK/0.060/ 
      DATA RKLARG/1.E6/ 
      DATA RAYON/6.37122E06/, GRAV/9.80616/ 
      DATA ALFAH/1.0/, BEEM/10.0/, VKC/0.4/ 
      DATA VMIN/2.0/, XINGL/100.0/, DVMINS/7.9/ 
      DATA FAC,GAMRH,GAMRM /32.5, 9.33, 9.33/ 
  
C     * FAC=3**(3/2) / VKC**2, GAMRH=SQRT(87), GAMRM=SQRT(87).
C     * NOTE  ALFAH MUST BE .GE. 1. 
C-----------------------------------------------------------------------
C     * EVALUATE THE VENTILATION VELOCITY.
  
      DO 02 I=IL1,IL2 
         VMODL(I)=MAX(VMIN,SQRT(U(I,ILEV)**2+V(I,ILEV)**2)) 
   02 CONTINUE
C===================================================================
C     * EVALUATE THE SURFACE DRAG COEFFICIENTS, (CDH,CDM).
C       ------------------------------------------------- 
C     * OVER OCEAN (HOMOGENEOUS TURB.), ACCOUNT FOR THE FACT THAT 
C     * THE MODEL LEVELS OF WIND AND TEMP. ARE ABOVE 10M. 
  
      CT       =1.15E-3 
      DO 30 I=IL1,IL2 
         CDOCEAN  =CT*VKC/(VKC+SQRT(CT)*LOG((1.-SHJ(I,ILEV))*1.E3))
         FM2      =(1.+SQRT(CT)*LOG((1.-SGJ(I,ILEV))*1.E3)/VKC)**2 
         CDO(I)   =MERGE(CDOCEAN,  MAX(DRAG(I),CT),  GC(I).EQ.0.) 
         IF(GC(I).EQ.1.) CDO(I)=CDO(I)*2. 
         EPS(I)  =MERGE(0.3,  0.,  GC(I).EQ.0.) 
         RATMOH(I)=MERGE( 
     1                0.001*MAX(0.9,0.61+0.063*VMODL(I))/(CDO(I)*FM2),
     2                1.0,
     3                GC(I).EQ.0. ) 
   30 CONTINUE
  
C     * RIB IS THE BULK RICHARDSON'S NUMBER.
C     * ZO IS A NORMALIZED MIXING HEIGHT, VKC IS THE VON KARMAN CONSTANT
C     * SLTHICK IS THE SIGMA THICKNESS OF THE SFC LAYER,
C     * DEFINE TO GIVE DZ(MOM)**2/DZ(THERM) OF OLD STAGGERED MODEL. 
  
      DO 40 I=IL1,IL2 
         RIB(I)=-RGAS*(GT(I)-TH(I,ILEV)/SHXKJ(I,ILEV))*SLTHICK /
     1          (VMODL(I)**2) 
         ZO(I) =EXP(-VKC/SQRT(CDO(I)))
   40 CONTINUE
  
C     * EVALUATE THE SURFACE DRAG COEFFICIENT FOR UNSTABLE
C     * CASE (RIB < 0.0) AND FOR THE STABLE CASE IN WHICH 
C     * 0.0 <= RIB < (2./(EPS*BEEM)) .
C     * ALSO EVALUATE SCREEN-LEVEL VARIABLES BASED ON ACTUAL SURFACE
C     * ROUGHNESS HEIGHT ZRUF (MINIMUM OF 2 METRES).
  
      DO 50 I=IL1,IL2 
         X(I)=MERGE(0., 
     1        (1.-EPS(I)*0.5*BEEM*RIB(I))**2, 
     2        0.5*RIB(I)*EPS(I)*BEEM.GE.1.) 
   50 CONTINUE
  
      DO 60 I=IL1,IL2 
         CDH(I)=MERGE(
     1          CDO(I)*(1.-ALFAH*BEEM*RIB(I)/(1.+CDO(I)*ALFAH*
     2                FAC*BEEM*SQRT(ABS(RIB(I))/ZO(I))/GAMRH)), 
     3          CDO(I)*X(I)/(1.+(1.   -EPS(I))    *BEEM*RIB(I)),
     4          RIB(I).LT.0.) 
         CDM(I)=MERGE(
     1          CDO(I)*(1.-      BEEM*RIB(I)/(1.+CDO(I)*
     2                FAC*BEEM*SQRT(ABS(RIB(I))/ZO(I))/GAMRM)), 
     3          CDO(I)*X(I)/(1.+(1.   -EPS(I))    *BEEM*RIB(I)),
     4          RIB(I).LT.0.) 
         ZLEV=-RGAS*TH(I,ILEV)*LOG(SHJ(I,ILEV))/GRAV 
         ZRUF=ZO(I)*ZLEV
         ZSCRN=MAX(ZRUF,2.0)
         RATFC=LOG(ZSCRN/ZRUF)/VKC 
         RATFC1=RATFC*SQRT(CDM(I))
         IF(RIB(I).GE.0.0)  THEN
            RATIO=RATFC1
         ELSE 
            RATIO=RATFC1*CDH(I)/CDM(I)
            RATIO=MIN(RATIO,(ZSCRN/ZLEV)**(1./3.))
         ENDIF
         ST(I)=GT(I)-(MIN(RATIO,1.))*(GT(I)-TS(I))
         SU(I)=RATFC1*U(I,ILEV) 
         SV(I)=RATFC1*V(I,ILEV) 
         SQ(I)=Q(I,ILEV)+EF(I)*(QG(I)-Q(I,ILEV))*MIN(RATIO,1.)
   60 CONTINUE
  
      DO 70 I=IL1,IL2 
         CDM(I)=CDM(I)*RATMOH(I)
   70 CONTINUE
  
      DO 80 I=IL1,IL2 
         CDVLH(I)=CDH(I)*VMODL(I) 
         CDVLM(I)=CDM(I)*VMODL(I) 
   80 CONTINUE
  
C-----------------------------------------------------------------------
C     * CALCULATE THE VERTICAL DIFFUSION COEFFICIENTS, (RKH,RKM). 
C       --------------------------------------------------------
  
C     * DVDS = MOD(DV/DSIGMA) AT THE INTERFACE OF MOMENTUM LAYERS.
  
      ILEVM=ILEV-1
      DO 100 L=1,ILEVM
         DO 100 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))
  100    CONTINUE 
      DO 110 I=IL1,IL2
         DVDS(I,ILEV)=DVDS(I,ILEVM) 
  110 CONTINUE
  
C     * INTERPOLATE DVDS TO INTERFACE OF THERMODYNAMIC LAYERS.
C     * RMS WIND SHEAR, DVMINS=(R*T/G)*(DV/DZ)=7.9 FOR DV/DZ=1.MS-1KM-1.
  
      DO 140 L=ILEV,2,-1
         DO 120 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) ) 
  120    CONTINUE 
  140 CONTINUE
  
      DO 260 L=2,ILEV 
  
C        * EVALUATE DTTDS=D(THETA)/D(SIGMA) 
C        * AT INTERFACE OF THERMODYNAMIC LAYERS,
  
         DO 210 I=IL1,IL2 
            DTTDS(I,L)=(TH(I,L)/SHXKJ(I,L)-TH(I,L-1)/SHXKJ(I,L-1))
     1                /(SHJ(I,L)-SHJ(I,L-1))
  210    CONTINUE 
  
C        * CALCULATE THE DIFFUSION COEFFICIENTS (RKH,RKM) 
C        * AT THE INTERFACE OF THERMODYNAMIC LAYERS.
  
         DO 220 I=IL1,IL2 
            RI(I)  =-RGAS*SHTXKJ(I,L)*DTTDS(I,L)
     1              /(SHBJ(I,L-1)*DVDS(I,L)**2) 
            DVDZ(I)=DVDS(I,L)*GRAV*SHBJ(I,L-1) / (RGAS*TF(I,L)) 
  220    CONTINUE 
  
C        * HEAT: FINITE STABILITY CUTOFF DEPENDING ON LOCAL SIGMA.
C        * UNSTABLE CASE (RI <= 0.) - STABLE CASE (RI > 0.) . 
  
         DO 250 I=IL1,IL2 
            EL  =(VKC*RGAS*273.*LOG(SHBJ(I,L-1)) / 
     1           (GRAV-VKC*RGAS*273.*LOG(SHBJ(I,L-1))/XINGL))**2 
            EPSS(I) = 0.3 
            X(I)=MERGE(0.,
     1           (1.-EPSS(I)*0.5*BEEM*RI(I))**2,
     2           0.5*RI(I)*EPSS(I)*BEEM.GT.1.)
            RKH(I,L)=MERGE( 
     1               EL*DVDZ(I)*(1.-ALFAH*GAMRH*BEEM*RI(I) /
     2                     (GAMRH+ALFAH*BEEM*SQRT(ABS(RI(I))))),
     3               EL*DVDZ(I)*X(I)/(1.+(1.   -EPSS(I))* 
     4                     BEEM*RI(I)), 
     5               RI(I).LT.0.) 
  250    CONTINUE 
  
C        * MOMENTUM: NO STABILITY CUTOFF. 
C        * UNSTABLE CASE (RI <= 0.) - STABLE CASE (RI > 0.) . 
  
         DO 255 I=IL1,IL2 
            EL  =(VKC*RGAS*273.*LOG(SHBJ(I,L-1)) / 
     1           (GRAV-VKC*RGAS*273.*LOG(SHBJ(I,L-1))/XINGL))**2 
            EPSS(I)=0.
            X(I)=MERGE(0.,
     1           (1.-EPSS(I)*0.5*BEEM*RI(I))**2,
     2           0.5*RI(I)*EPSS(I)*BEEM.GT.1.)
            RKM(I,L)=MERGE( 
     1               EL*DVDZ(I)*(1.-      GAMRM*BEEM*RI(I) /
     2                    (GAMRM+      BEEM*SQRT(ABS(RI(I))))), 
     3               EL*DVDZ(I)*X(I)/(1.+(1.   -EPSS(I))* 
     4                     BEEM*RI(I)), 
     5               RI(I).LT.0.) 
  255    CONTINUE 
  260 CONTINUE
  
      DO 270 I=IL1,IL2
         RKM(I,1)=RKM(I,2)
  270 CONTINUE
  
C     * DEFINE RKQ ALSO AT INTERFACE OF THERMODYNAMIC LAYERS. 
  
      DO 315 L=2,ILEV 
         DO 315 I=IL1,IL2 
            RKQ(I,L)=RKH(I,L) 
  315 CONTINUE
  
C     * INTERPOLATE RKM TO INTERFACE OF MOMENTUM LAYERS.
  
      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     * SET A LOWER BOUND TO RK.
  
      DO 330 L=1,ILEVM
         DO 330 I=IL1,IL2 
            RKM(I,L  )=MAX(RKMIN ,RKM(I,L  )) 
            RKH(I,L+1)=MAX(RKMIN ,RKH(I,L+1)) 
            RKQ(I,L+1)=MAX(RKQMIN,RKQ(I,L+1)) 
  330 CONTINUE
C-----------------------------------------------------------------------
C     * FIND TOP OF BOUNDARY LAYER (AFTER ADDING SURFACE HEAT FLUX).
C       --------------------------
  
C     * GOING UPWARD TO NEXT INTERFACE OF THERMODYNAMIC LAYER.
  
      DO 340 I=IL1,IL2
         CONST=2.*GRAV*DELT*SHXKJ(I,ILEV)/RGAS
         HEAT(I)=CONST*CDVLH(I)/TH(I,ILEV)
  340 CONTINUE
  
      DO 400 I=IL1,IL2
         HEATI=HEAT(I)
         IF(TH(I,ILEV)/SHXKJ(I,ILEV).GE.GT(I)) GO TO 400
  
         DO 370 L=ILEVM,1,-1
            TTADJ=TH(I,L)/SHXKJ(I,L)
            IF(TTADJ.GE.GT(I)) GO TO 380
            HEATIL=HEATI*(GT(I)-TTADJ)
  
            QUE=0.
            DO 350 K=L+1,ILEV 
               QUE = QUE+(TTADJ*SHXKJ(I,K)-TH(I,K))*DSHJ(I,K) 
  350       CONTINUE
  
            IF(QUE.GT.HEATIL) GO TO 380 
  370    CONTINUE 
  
C     * SO TOP P.B.L. LIES BETWEEN L+1/2 AND L+3/2. 
C     * RESET K TO LARGE VALUE FROM L+2 DOWNWARD. 
  
  380    LPBL=L+2 
         IF(LPBL.EQ.ILEV+1) GO TO 400 
         DO 390 L=LPBL,ILEV 
            RKH(I,L)=RKLARG 
  390    CONTINUE 
  
  400 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 IMPLVD6 MUST BE RE-
C     * DONE OUTSIDE. NO SUCH PROCESS IS REQUIRED FOR THE UFSROL, 
C     * VFSROL FIELDS SAVED EVERY ISBEG STEPS.
  
      DO 490 I=IL1,IL2
         RAUS(I)=PRESSG(I)/GRAV 
  490 CONTINUE
      TODT=2.*DELT
      NOGWDRG=.TRUE.
      IGW=.FALSE.
  
      CALL ABCVDM6 (A,B,C,CL,CDVLM,GRAV,IL1,IL2,ILG,ILEV, 
     1              RGAS,RKM,SGJ,SGBJ,SHJ,TSGB,TODT,NOGWDRG)
  
      CALL IMPLVD6 (A,B,C,CL,U,CL,IL1,IL2,ILG,ILEV,TODT,
     1              TEND,UFS,DSGJ,RAUS,WORK,VINT) 
  
      DO 498 I=IL1,IL2
         UFS(I)=UFS(I)-VINT(I)*RAUS(I)*(1.-SAVEGG)
         UFSROL(I)=UFSROL(I)+VINT(I)*RAUS(I)*SAVEBEG
  498 CONTINUE
  
      DO 500 L=1,ILEV 
         DO 500 I=IL1,IL2 
            UTG(I,L)=UTG(I,L)+(TEND(I,L)+UTENDGW(I,L))
  500 CONTINUE
  
      CALL IMPLVD6 (A,B,C,CL,V,CL,IL1,IL2,ILG,ILEV,TODT,
     1              TEND,VFS,DSGJ,RAUS,WORK,VINT) 
  
      DO 505 I=IL1,IL2
         VFS(I)=VFS(I)-VINT(I)*RAUS(I)*(1.-SAVEGG)
         VFSROL(I)=VFSROL(I)+VINT(I)*RAUS(I)*SAVEBEG
  505 CONTINUE
  
      DO 510 L=1,ILEV 
         DO 510 I=IL1,IL2 
            VTG(I,L)=VTG(I,L)+(TEND(I,L)+VTENDGW(I,L))
  510 CONTINUE
  
      IF(.NOT.IGW) THEN 
C 
C        * IF USING ENVELOPE LAUNCHING HEIGHTS (ENVELOP=.TRUE. SO
C        * IGW=.FALSE.) 
C        * UPDATE (UFS,VFS) TO INCLUDE GRAVITY-WAVE DRAG EFFECTS. 
C        * THIS IS NOT DONE FOR (UFSROL,VFSROL) SINCE THESE ARE USED FOR
C        * THE OCEAN MODEL. THE DIFFERENCE BETWEEN THE TWO, IN A TIME-
C        * AVERAGED SENSE, GIVES THE CONTRIBUTION OF GRAVITY-WAVE DRAG
C        * ALONE, WHILE THE LATTER GIVES THE CONTRIBUTION OF TURBULENCE.
C 
         DO 506 I=IL1,IL2
            UGWTS(I)=0.
            VGWTS(I)=0.
  506    CONTINUE
C 
         DO 508 L=1,ILEV 
            DO 507 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) 
  507       CONTINUE 
  508    CONTINUE
C 
         DO 509 I=IL1,IL2 
            UFS(I)=UFS(I)+UGWTS(I)*SAVEGG 
            VFS(I)=VFS(I)+VGWTS(I)*SAVEGG 
  509    CONTINUE 
      ENDIF 
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      IF (ITRAC.GT.0)   THEN                                            
                                                                        
C         * CALCULATE THE COEFFICIENT MATRIX FOR TRACER DIFFUSION.                
C         * GET TENDENCIES FROM VERTICAL TRACER DIFFUSION.                        
                                                                        
          DO 515 I=IL1,IL2                                              
             UN(I)   =1.                                                
             ZER(I)  =0.                                                
             CDVLT(I)=0.                                                
  515     CONTINUE                                                      
                                                                        
          CALL ABCVDQ6 (                                                
     1                  A,B,C,CL, UN ,CDVLT,GRAV,                       
     2                  IL1,IL2,ILG,ILEV,ILEV+1,ILEV,                   
     3                  RGAS,RKH,SHTJ,SHJ,DSHJ,                         
     4                  TH(1,ILEV),TF,TODT)                             
                                                                        
          DO 532 N=1,NTRAC                                              
                                                                        
C            * INITIALISATION DES FLUX DE SURFACES A ZERO                         
                                                                        
             DO 516 I=IL1,IL2                                           
                XFS(I,N)=0.                                             
  516        CONTINUE                                                   
                                                                        
             CALL IMPLVD6(A,B,C,CL,XROW(1,2,N),ZER,IL1,IL2,ILG,ILEV,    
     1                    TODT,TEND,XFS(1,N),DSHJ,RAUS,WORK,VINT)       
                                                                        
             DO 525 L=1,ILEV                                            
             DO 525 I=IL1,IL2                                           
                XROW(I,L+1,N)=XROW(I,L+1,N)+TODT*TEND(I,L)              
  525        CONTINUE                                                   
                                                                        
             DO 530 I=IL1,IL2                                           
                XFS(I,N)=XFS(I,N)-VINT(I)*RAUS(I)*(1.-SAVEGG)           
                XFX(I,N)=-VINT(I)*RGAS*TH(I,ILEV)/(GRAV*SHXKJ(I,ILEV))  
  530        CONTINUE                                                   
  532     CONTINUE                                                      
                                                                        
      ENDIF                                                             
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     * CALCULATE THE COEFFICIENT MATRIX FOR MOISTURE DIFFUSION.
C     * GET TENDENCIES FROM VERTICAL MOISTURE DIFFUSION.
C     * EVALUATE QFS. 
  
      CALL ABCVDQ6 (
     1              A,B,C,CL, EF ,CDVLH,GRAV, 
     2              IL1,IL2,ILG,ILEV,ILEV+1,LEVS, 
     3              RGAS,RKQ,SHTJ,SHJ,DSHJ, 
     4              TH(1,ILEV),TF,TODT) 
  
      L=ILEV-LEVS+1 
      CALL IMPLVD6(A(1,L),B(1,L),C(1,L),CL,Q(1,L),QG,IL1,IL2,ILG,LEVS,
     1             TODT,TEND(1,L),QFS,DSHJ(1,L),RAUS,WORK,VINT) 
  
      DO 535 LL=L,ILEV
         DO 535 I=IL1,IL2 
            Q(I,LL)=Q(I,LL)+TODT*TEND(I,LL) 
  535 CONTINUE
  
      DO 540 I=IL1,IL2
         QFS(I)=QFS(I)-VINT(I)*RAUS(I)*(1.-SAVEGG)
         QFX(I)=-VINT(I)*RGAS*TH(I,ILEV)/(GRAV*SHXKJ(I,ILEV)) 
  540 CONTINUE
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  
C     * CALCULATE THE COEFFICIENT MATRIX FOR HEAT DIFFUSION.
C     * GET TENDENCIES FROM VERTICAL HEAT DIFFUSION.
C     * EVALUATE HFS. 
  
      CALL ABCVDH6 (
     1              A,B,C,CL,CDVLH,GRAV,
     2              IL1,IL2,ILG,ILEV,ILEV+1,
     3              RGAS,RKH,SHTJ,SHTXKJ,SHJ,SHXKJ,DSHJ,
     4              TH(1,ILEV),TF,TODT) 
  
      DO 545 I=IL1,IL2
         RAUS(I)=RAUS(I)*CPRES
  545 CONTINUE
  
      CALL IMPLVD6 (A,B,C,CL,TH,GT,IL1,IL2,ILG,ILEV,TODT, 
     1              TEND,HFS,DSHJ,RAUS,WORK,VINT) 
  
      DO 550 I=IL1,IL2
         HFS(I)=HFS(I)-VINT(I)*RAUS(I)*(1.-SAVEGG)
         TFX(I)=-VINT(I)*RGAS*TH(I,ILEV)/(GRAV*SHXKJ(I,ILEV)) 
  550 CONTINUE
  
      DO 555 L=1,ILEV 
         DO 555 I=IL1,IL2 
            TH(I,L)=TH(I,L)+TODT*TEND(I,L)
  555 CONTINUE
  
      RETURN
C-----------------------------------------------------------------------
      END 
