      SUBROUTINE GWDFL6O (U,V,TH,TF,TSG,ENV,GC,SGJ,SHJ,DSGJ,
     1                   SHTXKJ,SHXKJ,UTENDGW,VTENDGW,UROW,VROW,
     2                   RGAS,RGOCP,DELT,ILEV,LREFP,IL1,IL2,ILG,
     3                   DAMPFAC,ENVELOP, 
     4                   BVFREQ,BVFREQG,VELN,UB,VB,UBG,VBG, 
     5                   VMODB,VMODBG,DEPFAC,HEIGHT,HITESQ, 
     6                   AMPBSQ,AMPDIF,DENFAC,TFG,UTEND,VTEND,
     7                   DGW,DGWG,DRAG) 

C     * JANUARY 13/95 - M. LAZARE, ELLISON CHAN -
C     *                 OPTIMIZATIONS DONE IN GWDFLX8 MERGED WITH GWDFLX6
C     *                 SO LATTER'S RESULTS REPRODUCED EXACTLY.
C     *
C     * OCTOBER 20/88 - N.MCFARLANE. PREVIOUS VERSION GWDFLX6.
C 
C     * LREFP IS THE INDEX OF THE MODEL LEVEL BELOW THE REFERENCE LEVEL 
  
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)

C     * I/O ARRAYS PASSED FROM MAIN.
  
      REAL U(ILG,ILEV),       V(ILG,ILEV),       TH(ILG,ILEV),
     1     TF(ILG,ILEV),      TSG(ILG,ILEV),     GC(ILG), 
     2     UTENDGW(ILG,ILEV), VTENDGW(ILG,ILEV), ENV(ILG),
     3     UROW(ILG,ILEV),    VROW(ILG,ILEV)
  
C     * VERTICAL POSITIONNING ARRAYS. 
  
      REAL SGJ(ILG,ILEV),     SHJ(ILG,ILEV),     SHTXKJ(ILG,ILEV),
     1     SHXKJ(ILG,ILEV),   DSGJ(ILG,ILEV)
  
C     * LOGICAL SWITCHES TO CONTROL ROOF DRAG AND ENVELOP GW DRAG.
  
      LOGICAL DAMPFAC, ENVELOP
  
C     * WORK ARRAYS (ALSO PASSED FROM MAIN).
  
      REAL BVFREQ(ILG,ILEV),  BVFREQG(ILG,ILEV), VELN(ILG,ILEV),
     1     UB(ILG),           VB(ILG),           UBG(ILG), 
     2     VBG(ILG),          VMODB(ILG),        VMODBG(ILG),  
     3     DEPFAC(ILG,ILEV),  HEIGHT(ILG),       HITESQ(ILG),
     4     AMPBSQ(ILG),       AMPDIF(ILG),       DENFAC(ILG),
     5     TFG(ILG,ILEV),     UTEND(ILG,ILEV),   VTEND(ILG,ILEV), 
     6     DGW(ILG,ILEV),     DGWG(ILG,ILEV),    DRAG(ILG)
  
C     * CONSTANTS VALUES DEFINED IN DATA STATEMENT ARE :  
  
C     * VMIN     = MIMINUM WIND IN THE DIRECTION OF REFERENCE LEVEL 
C     *            WIND BEFORE WE CONSIDER BREAKING TO HAVE OCCURED.
C     * DMPSCAL  = DAMPING TIME FOR GW DRAG IN SECONDS. 
C     * TAUFAC   = 1/(LENGTH SCALE).
C     * HMIN     = MIMINUM ENVELOPE HEIGHT REQUIRED TO PRODUCE GW DRAG. 
C     * V0       = VALUE OF WIND THAT APPROXIMATES ZERO.
  
      DATA    VMIN  /    2.0 /, V0       / 1.E-10 /,
     1        TAUFAC/ 1.5E-5 /, HMIN     /   10.0 /,
     2        ZERO  /    0.0 /, UN       /    1.0 /, DEUX  /    2.0 /,
     3        GRAV  /9.80616 /, DMPSCAL  / 6.5E+6 / 
  
C-----------------------------------------------------------------------
      ILEVM = ILEV - 1
      LEN   = IL2 - IL1 + 1 
      LREF  = LREFP-1 
  
C     * VMOD IS THE REFERENCE LEVEL WIND AND ( UB, VB)
C     * ARE IT'S UNIT VECTOR COMPONENTS.
  
      DO 100 I=IL1,IL2
        UB(I)    = U(I,LREF)
        VB(I)    = V(I,LREF)
        VMODB(I) = SQRT (UB(I)**2 + VB(I)**2) 
        VMODB(I) = MAX (VMODB(I), UN) 
        UB(I)    = UB(I)/VMODB(I) 
        VB(I)    = VB(I)/VMODB(I) 
  100 CONTINUE
  
C     * DRAG REFERENCES THE POINTS WHERE GW CALCULATIONS WILL BE DONE,
C     * THAT IS (A- IF OVER LAND, B- IF BOTTOM WIND > VMIN, C- IF WE
C     * ASK FOR IT AND C- IF ENVELOPPE HEIGHT >= HMIN ) 
C     * DRAG=1. FOR POINTS THAT SATISFY THE ABOVE, ELSE DRAG=0.
  
      LENGW = 0

      IF(ENVELOP)                                                  THEN
        DO 110 I=IL1,IL2
          IF(GC(I).EQ.-1. .AND.
     1       VMODB(I).GT.VMIN .AND. ENV(I).GE.HMIN)          THEN 
            LENGW   = LENGW + 1
            DRAG(I) = 1.
          ELSE
            DRAG(I) = 0.
          ENDIF 
  110   CONTINUE
      ENDIF

C     * INITIALIZE NECESSARY ARRAYS.
  
      DO 120 L=1,ILEV 
      DO 120 I=IL1,IL2
        UTENDGW(I,L) = ZERO 
        VTENDGW(I,L) = ZERO 
        UROW(I,L)    = U(I,L) 
        VROW(I,L)    = V(I,L) 
  120 CONTINUE
  
      IF (LENGW.EQ.0)                                 GO TO 300

C     * CALCULATE  B V FREQUENCY IN DRAG REGIONS
  
      DO 140 L=2,ILEV 
      DO 140 I=IL1,IL2
        IF(DRAG(I).EQ.1.)                                           THEN
          DTTDSF=(TH(I,L)/SHXKJ(I,L)-TH(I,L-1)/
     1           SHXKJ(I,L-1))/(SHJ(I,L)-SHJ(I,L-1)) 
          DTTDSF=MIN(DTTDSF, -5./SGJ(I,L)) 
          BVFREQG(I,L)=SQRT(-DTTDSF*SGJ(I,L)*(SGJ(I,L)**RGOCP)/RGAS)
     1                     *GRAV/TSG(I,L) 
        ENDIF 
  140 CONTINUE

C     * DEFINE BVFREQ AT TOP AND SMOOTH BVFREQ.
C     * DEFINE SQUARE OF WIND MAGNITUDE RELATIVE TO REFERENCE LEVEL.

      DO 160 L=1,ILEV 
      DO 160 I=IL1,IL2
        IF(DRAG(I).EQ.1.)                                           THEN
          IF(L.EQ.1)                        THEN
            BVFREQG(I,L) = BVFREQG(I,L+1) 
          ENDIF
          IF(L.GT.1)                        THEN 
            RATIO=5.*LOG(SGJ(I,L)/SGJ(I,L-1))
            BVFREQG(I,L) = (BVFREQG(I,L-1) + RATIO*BVFREQG(I,L))
     1                     /(1.+RATIO)
          ENDIF
          VELTAN=U(I,L)*UB(I)+V(I,L)*VB(I)
          VELN(I,L)=MAX(VELTAN,V0)
        ENDIF 
  160 CONTINUE
  
C     * CALCULATE EFFECTIVE SQUARE LAUNCHING
C     * HEIGHT, REFERENCE B V FREQUENCY, ETC. 
C     * ENV(I) IS THE SUB-GRID SCALE VARIANCE FIELD. 
  
      DO 200 I=IL1,IL2
        IF(DRAG(I).EQ.1.)                                           THEN
          SIGB=SGJ(I,LREF) 
          BVFB = BVFREQG(I,LREF) 
          HSQMAX=    (VMODB(I)/BVFB)**2 
          HSQ=MIN(4.*ENV(I),HSQMAX/2.)
          HSCAL=RGAS*TF(I,LREF)/GRAV
          DFAC=BVFB*SIGB*VMODB(I)/HSCAL 
          AMPBSQ(I)=DFAC
          HITESQ(I)=HSQ
        ENDIF
  200 CONTINUE
  
C     * CALCULATE TERMS FROM THE BOTTOM-UP. 
  
      DO 240 L=LREF,1,-1                    
      DO 240 I=IL1,IL2
        IF(DRAG(I).EQ.1.)                                           THEN
          WIND=VELN(I,L) 
          BVF=BVFREQG(I,L) 
          HSCAL=RGAS*TF(I,L)/GRAV 
          HSQMAX=0.5*(WIND/BVF)**2
          IF (VELN(I,L) .LT. UN) HSQMAX=ZERO
          DFAC=BVF*SGJ(I,L)*WIND/HSCAL
          RATIO=AMPBSQ(I)/DFAC
          HSQ=MIN(RATIO*HITESQ(I),HSQMAX)
          HITESQ(I)=HSQ
          AMPBSQ(I)=DFAC
          DEPFAC(I,L)=TAUFAC*DFAC*HSQ 
        ENDIF 
  240 CONTINUE
  
C     * CALCULATE GW TENDENCIES (TOP AND BOTTOM LAYERS FIRST).
C     * BOTTOM LAYER KEEPS INITIALIZED VALUE OF ZERO.
C     * APPLY GW DRAG ON (UROW,VROW) WORK ARRAYS TO BE PASSED TO VRTDFS.
  
      DO 260 I=IL1,IL2
        IF(DRAG(I).EQ.1.)                                           THEN
          DEPFAC(I,LREFP)=DEPFAC(I,LREF) 
          WIND=VELN(I,1)
          ETA=ZERO 
          DFLXP=DEPFAC(I,2)-DEPFAC(I,1)
          IF(DFLXP.GT.ZERO) ETA=UN 
          DFAC=6.*DELT*DEPFAC(I,1)*ETA/WIND
          DENOM=2.*DSGJ(I,1)+DFAC
          TENDFAC=DFLXP/DENOM
          DENFAC(I)=DFAC*TENDFAC 
          UTENDGW(I,1)=-TENDFAC*UB(I) 
          VTENDGW(I,1)=-TENDFAC*VB(I) 
          UROW(I,1) = UROW(I,1) + DEUX*DELT*UTENDGW(I,1) 
          VROW(I,1) = VROW(I,1) + DEUX*DELT*VTENDGW(I,1) 
        ENDIF 
  260 CONTINUE
  
      DO 270 L=2,LREF 
      DO 270 I=IL1,IL2
        IF(DRAG(I).EQ.1.)                                           THEN
          WIND=VELN(I,L)
          ETA=ZERO 
          DFLXP=DEPFAC(I,L+1)-DEPFAC(I,L)
          DFLXM=DEPFAC(I,L)-DEPFAC(I,L-1)
          IF(DFLXP.GT.ZERO) ETA=UN 
          DFAC=6.*DELT*DEPFAC(I,L)*ETA/WIND
          DENOM=2.*DSGJ(I,L)+DFAC
          TENDFAC=(DFLXP+DFLXM+DENFAC(I))/DENOM
          DENFAC(I)=DFAC*TENDFAC 
          UTENDGW(I,L) = -TENDFAC*UB(I) 
          VTENDGW(I,L) = -TENDFAC*VB(I) 
          UROW(I,L) = UROW(I,L) + DEUX*DELT*UTENDGW(I,L) 
          VROW(I,L) = VROW(I,L) + DEUX*DELT*VTENDGW(I,L) 
        ENDIF 
  270 CONTINUE
  
C     * APPLY ROOF DRAG.
  
  300 CONTINUE

      IF (DAMPFAC)                                                THEN
          DAMPRF = UN
      ELSE
          DAMPRF = ZERO
      ENDIF
 
      DO 310 I=IL1,IL2
          VELMOD       = MAX (SQRT (UROW(I,1)**2 + VROW(I,1)**2), UN)
          DFAC         = DAMPRF / (DMPSCAL*MIN((75./VELMOD)**4,UN))
          DENOM        = UN + DEUX*DELT*DFAC
          UROW(I,1)    = UROW(I,1)/DENOM
          VROW(I,1)    = VROW(I,1)/DENOM
          UTENDGW(I,1) = UTENDGW(I,1) - DFAC*UROW(I,1)
          VTENDGW(I,1) = VTENDGW(I,1) - DFAC*VROW(I,1)
  310 CONTINUE
  
C-----------------------------------------------------------------------------
      RETURN
      END 
