      SUBROUTINE GWDFL10D(UROW,VROW,TH,TF,TSG,ENV,GC,SGJ,SHJ,DSGJ,
     1                    SHTXKJ,SHXKJ,UTENDGW,VTENDGW,
     2                    RGAS,RGOCP,ZTMST,ILEV,LREFP,IL1,IL2,ILG,IDAMP,
     3                    LONSL,LMAM,DAMPFAC,ENVELOP,IFIZGWD,
     4                    UGROW,VGROW,USROW,VSROW,UNOTND,VNOTND)
C
C     * OCT 02/2007 - M.LAZARE.  "COSMETIC" CHANGE TO PASS IN AND
C     *                          USE ADVECTIVE TIMESCALE ("ZTMST")
C     *                          RATHER THAN "DELT", TO WORK FOR BOTH
C     *                          SPECTRAL MODEL (2.*DELT) OR GEM (DELT).
C     *                          THUS, "DEUX" CAN BE REMOVED.
C     * JUN 30/2006 - M.LAZARE.  PREVIOUS VERSION GWDFL10C FOR GCM15F.
C     *                          INPUT WIND FIELDS ARE NOW {UROW,VROW}
C     *                          INSTEAD OF {U,V}, IN CONJUNCTION WITH
C     *                          CHANGE TO DEFINE {UROW,VROW} IN THE
C     *                          PHYSICS DRIVER SUCH THAT POSSIBLY
C     *                          CONVERT FROM NATIVE REAL TO REAL*8.
C     *                          ALSO CLEANER AND MORE USER-FRIENDLY.
C     *                          IN ADDITION, INITIALIZATION TO ZERO
C     *                          OF {UTENDGW,VTENDGW} REMOVED SINCE NOW
C     *                          DONE EXPLICITLY IN PHYSICS).
C 
      IMPLICIT REAL   (A-H,O-Z),
     +INTEGER (I-N)

C     * LREFP IS THE INDEX OF THE MODEL LEVEL BELOW THE REFERENCE LEVEL 
  
C     * I/O ARRAYS PASSED FROM MAIN.
  
      REAL   UROW   (ILG,ILEV), VROW   (ILG,ILEV),    TH(ILG,ILEV),
     1       TF     (ILG,ILEV), TSG    (ILG,ILEV), 
     2       UTENDGW(ILG,ILEV), VTENDGW(ILG,ILEV)
      REAL   GC(ILG),           ENV(ILG)
  
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
C     * ARRAYS TO CALCULATE COMPONENT TENDENCIES 
C     * (ONLY USED IF IFIZGWD=1).
C
      REAL   UGROW(ILG,ILEV),   VGROW(ILG,ILEV),
     1       USROW(ILG,ILEV),   VSROW(ILG,ILEV)

C     * INTERNAL WORK ARRAYS.

C     * INITIAL WIND PROFILES USED IN SUBSEQUENT CALCULATIONS
C     * BEFORE {UROW,VROW} GET UPDATED.

      REAL    U(ILG,ILEV),     V(ILG,ILEV)

C     * OTHERS.
  
      REAL   TAU(ILG,ILEV),     BVFREQG(ILG,ILEV), VELN(ILG,ILEV),
     1       UB(ILG),           VB(ILG),           FRB(ILG),
     2       DENOFAC(ILG),      VMODB(ILG),        DNSWTCH(ILG), 
     3       DEPFAC(ILG,ILEV),  HITESQ(ILG),       DEPFACR(ILG),
     4       AMPBSQ(ILG),       FRMAX(ILG),        DENFAC(ILG),
     5       TFG(ILG,ILEV),     UTEND(ILG,ILEV),   VTEND(ILG,ILEV), 
     6       DGW(ILG,ILEV),     DGWG(ILG,ILEV),    DRAG(ILG)

      REAL   UMEAN(ILG,ILEV), VMEAN(ILG,ILEV)
      REAL   DZAV(ILG,ILEV)
      REAL   USUM(ILEV), VSUM(ILEV), DSUM(ILEV)
      REAL   UTSPNG(ILG,ILEV), VTSPNG(ILG,ILEV)
      REAL   UNOTND(ILG,ILEV),VNOTND(ILG,ILEV)
C
      LOGICAL LMAM  

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     * QVCRIT   = THRESHOLD WIND SPEED ABOVE WHICH ROOF DRAG IS
C     *            INITIATED.
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 /, QVC / 110.0 /,
     1        TAUFAC/ 1.0E-5 /, HMIN     /   10.0 /,
     2        ZERO  /    0.0 /, UN       /    1.0 /, 
     3        GRAV  /9.80616 /, DMPSCAL  / 6.5E+6 / 
C-----------------------------------------------------------------------
C     * RESET QVCRIT TO LOWER VALUE IF IN SPINUP PERIOD.
C
      IF(IDAMP.EQ.1) THEN
         QVCRIT=70.
      ELSE
         QVCRIT=QVC
      ENDIF
C
      ILEVM = ILEV - 1
      LEN   = IL2 - IL1 + 1 
      LREF  = LREFP-1 
      LREFM = LREFP-2                                                   
  
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)**2.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
        U      (I,L) = UROW(I,L)
        V      (I,L) = VROW(I,L)
  120 CONTINUE
  
      IF (LENGW.EQ.0)                                        GO TO 275

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. 
C     
      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)**2,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
               TAU(I,L)  =TAUFAC*DFAC*HSQ
          ENDIF
  240 CONTINUE
C     
      DO 250 IL=IL1,IL2
        IF(DRAG(IL).EQ.1.) THEN
          TAU(IL,LREFP)=TAU(IL,LREF)
          DEPFAC(IL,1)=TAU(IL,1)
        ENDIF
  250 CONTINUE
C
C     * DO 3-POINT SMOOTHING OF STRESS VECTOR.
C
      DO 255 L=2,LREF
      DO 255 IL=IL1,IL2
        IF(DRAG(IL).EQ.1.) THEN
          DEPFAC(IL,L)=(.25*TAU(IL,L-1)+
     1                 .5*TAU(IL,L)+.25*TAU(IL,L+1))
        ENDIF
  255 CONTINUE
C  
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.
C  
      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)
C         IF(DFLXP.GT.ZERO) ETA=UN 
          DFAC=3.*ZTMST*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) + ZTMST*UTENDGW(I,1) 
          VROW(I,1) = VROW(I,1) + ZTMST*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)
C         IF(DFLXP.GT.ZERO) ETA=UN 
          DFAC=3.*ZTMST*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) + ZTMST*UTENDGW(I,L) 
          VROW(I,L) = VROW(I,L) + ZTMST*VTENDGW(I,L) 
        ENDIF 
  270 CONTINUE

  275 CONTINUE

      IF(IFIZGWD.NE.0)                                      THEN

C       * STORE OROGRAPHIC GWD ADJUSTED WIND TO LATER DETERMINE
C       * ITS TENDENCY. 

        DO 280 L=1,ILEV
        DO 280 I=IL1,IL2
          UGROW(I,L)=UROW(I,L)
          VGROW(I,L)=VROW(I,L) 
  280   CONTINUE
      ENDIF
C
C     * USROW,UVROW CONTAIN NON-ORO GWD TENDENCY 
C     * APPLY NON-ORO TENDENCY HERE 
C
      DO 285 L=1,ILEV  
      DO 285 I=IL1,IL2  
          UROW(I,L) = UROW(I,L) + ZTMST*UNOTND(I,L)
          VROW(I,L) = VROW(I,L) + ZTMST*VNOTND(I,L)
          UTENDGW(I,L)=UTENDGW(I,L)+UNOTND(I,L)
          VTENDGW(I,L)=VTENDGW(I,L)+VNOTND(I,L) 
  285 CONTINUE
C
      IF(IFIZGWD.NE.0)                                      THEN

C       * STORE NON-OROGRAPHIC GWD ADJUSTED WIND TO LATER DETERMINE
C       * ITS TENDENCY. 

        DO 290 L=1,ILEV
        DO 290 I=IL1,IL2
          USROW(I,L)=UROW(I,L)
          VSROW(I,L)=VROW(I,L) 
  290   CONTINUE
      ENDIF
  
C     * APPLY ROOF DRAG.
  
      IF (DAMPFAC)                                                THEN
      IF(.NOT.LMAM.OR.(LMAM.AND.IDAMP.EQ.1))       THEN
        DO 310 L=1,ILEV 
        DO 310 I=IL1,IL2
          VELMOD       = SQRT (UROW(I,L)**2 + VROW(I,L)**2) 
          IF((L.EQ.1 .OR. SGJ(I,L).LE.0.010) .AND. VELMOD.GT.VMIN)  THEN
            VELRAT       = MAX(VELMOD/QVCRIT-1.,ZERO)
            DFAC         = VELRAT/(ZTMST)
            DENOM        = UN + ZTMST*DFAC
            UROW(I,L)    = UROW(I,L)/DENOM
            VROW(I,L)    = VROW(I,L)/DENOM
            UTENDGW(I,L) = UTENDGW(I,L) - DFAC*UROW(I,L)
            VTENDGW(I,L) = VTENDGW(I,L) - DFAC*VROW(I,L)
          ENDIF 
  310   CONTINUE
      ELSE

C       * CALCULATE NUMBER OF CHAINED LATITUDES

        NLATJ = IL2/LONSL

C       * INITIALIZE ZONAL MEAN WIND TO ZERO. 

        DO 602 L=1,ILEV
        DO 602 I=IL1,IL2
          UMEAN(I,L) = ZERO
          VMEAN(I,L) = ZERO
          DZAV(I,L)  = ZERO
  602   CONTINUE

C       * CALCULATE ZONAL MEAN WINDS IF USING NONZONAL SPONGE.
C       * THE USE OF CHAINED LATITUDES MEANS THAT THE ZONAL AVERAGE 
C       * MUST BE CALCULATED FOR EACH LATITUDE.
C
        ZSPONG = 80.0
        FLONSL = FLOAT(LONSL)
        DO 620 J=1,NLATJ 
          DO 615 L=1,ILEV
            USUM(L) = ZERO
            VSUM(L) = ZERO
            DSUM(L) = ZERO
            DO 605 K=1,LONSL
              I = (J-1)*LONSL+K
              USUM(L) = USUM(L) + UROW(I,L)
              VSUM(L) = VSUM(L) + VROW(I,L)
              ZSCHT   = MAX(-7.*LOG(SGJ(I,L))-ZSPONG,ZERO)
              DSUM(L) = ZSCHT + DSUM(L)
  605       CONTINUE
C
            DO 610 K=1,LONSL
              I = (J-1)*LONSL+K
              UMEAN(I,L) = USUM(L) / FLONSL
              VMEAN(I,L) = VSUM(L) / FLONSL
              DZAV(I,L)  = DSUM(L) / FLONSL
  610       CONTINUE
  615     CONTINUE
  620   CONTINUE

C       * USE PARABOLIC-SHAPED PROFILE.

        DO 630 L=1,ILEV
        DO 630 I=IL1,IL2
            DFAC   = (3.6E-5)*(DZAV(I,L)/20.)**2
            DENOM  = 1.0 + ZTMST*DFAC
            UROW(I,L)    = ( UROW(I,L)+ZTMST*DFAC*UMEAN(I,L) )
     1                      / DENOM
            VROW(I,L)    = ( VROW(I,L)+ZTMST*DFAC*VMEAN(I,L) )
     1                      / DENOM

C           SPONGE AND EVD DRAG applied on NON-ZONAL wind only.

            UTSPNG(I,L)  =  - DFAC * ( UROW(I,L)-UMEAN(I,L) )  
            VTSPNG(I,L)  =  - DFAC * ( VROW(I,L)-VMEAN(I,L) )
            UTENDGW(I,L) = UTENDGW(I,L) + UTSPNG(I,L) 
            VTENDGW(I,L) = VTENDGW(I,L) + VTSPNG(I,L)
  630   CONTINUE

      ENDIF

      ENDIF 
C-----------------------------------------------------------------------------
      RETURN
      END 
