      SUBROUTINE GWDFL17(UROW,VROW,TSG,TSGB,ENV,GAM,PSIM_EAST,
     1             ALPHA,DELTA,SIGX,PHIE,KOUNT,IDAMP,
     2             ITUNE,SGJ,SGBJ,DSGJ,PRESSG,
     3             UTENDGW,VTENDGW,UTENDDW,VTENDDW,UTENDLLB,VTENDLLB,
     4             FTOX,FTOY,TDOX,TDOY,
     5             RGAS,RGOCP,ZTMST,ILEV,IL1,IL2,ILG,LONSL,
     6             LMAM,DAMPFAC,ENVELOP,IFIZGWD,UGROW,VGROW,
     7             USROW,VSROW,UFROW,VFROW,UNOTND,VNOTND)
C
C     * OCT 2018  - M.LAZARE. - WORK ARRAYS {UFROW,VFROW} ADDED TO PROPERLY CALCULATE
C     *                         EVD/SPONGE TENDENCIES.
C     * SEPT 2018 - J. COLE   - APPLY EXCESS VELOCITY DRAG AT ALL LEVELS INSTEAD OF 
C                               ABOVE SOME LEVEL.
C     * APR 2018 - M.LAZARE.  - REMOVE UNUSED DEFINITION OF "DZ" IN CODE JOHN ADDED. 
C     * FEB 2016 - J.SCINOCCA - BACKPORT ALL CHANGES FROM OFFICIAL CODE INTO
C                               CODE THAT HAS BEEN USED FOR CMAM WORK OVER THE
C                               PAST 10 YEARS - MOSTLY DIAGNOSTIC OUTPUT
C                             - CMAM CHANGES INCLUDE EVALUATION OF THE MOMENTUM
C                               FLUXES ON MOMENTUN BASE LEVELS
C                             - ADDITIONAL OUTPUT DIAGNOSTICS:
C     FTOX,FTOY - (ILG,ILEV) VERTICAL PROFILE OF NET EASTWARD AND NORTHWARD GW FLUXES
C     TAUP_X,TAUP_Y,TAUM_X,TAUM_Y - (ILG,ILEV) VERTICAL PROFILE OF NET EASTWARD AND 
C                               NORTHWARD GW FLUXES FOR THE "P" AND "M" WAVES IN THE SCHEME
C     TDOX,TDOY - (ILG,ILEV) VERTICAL PROFILE OF NET EASTWARD AND NORTHWARD GW TENDENCIES
C                            DUE TO THE FREELY PROPAGATING WAVE COMPONENT 
C     TNDP_X,TNDP_Y,TNDM_X,TNDM_Y - (ILG,ILEV) VERTICAL PROFILE OF NET EASTWARD AND 
C                            NORTHWARD GW FLUXES DUE TO THE FREELY PROPAGATING WAVE
C                            COMPONENT FOR THE "P" AND "M" WAVES IN THE SCHEME 
C                             - ORIENTATION ANGLE OR REPRESENTATIVE ELLIPSE NOW
C                               RELATIVE TO DUE EAST RATHER THAN LOCAL ZONAL
C                               DIRECTION OF NATIVE GRID
C
C     BACKPORTED CHANGES INCLUDE:
C           * FEB 10/2015 - M.LAXARE.  New version for gcm18:
C           *                          - Remove unused GCROW.
C           * APR 29/2012 - J.SCINOCCA/PREVIOUS VERSION GWDFL16 FOR GCM16/GCM17:
C           *                          - UNUSED CODE REMOVED.
C           *                          - TUNING ADJUSTMENTS TO GPHIL AND FCRIT.
C           *                          - DEPOSIT OUTGOING GWD IN TOP LAYER.
C           * OCT 02/2007 - M.LAZARE.  PREVIOUS VERSION GWDFL15 FOR GCM15H/I:
C           *                          "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 GWDFL14 FOR GCM15F/G.
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           * JUN 27/2006 - M.LAZARE.  UNUSED COSJ REMOVED.
C           * JUN 19/2006 - M.LAZARE.  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 22/2006 - M.LAZARE.  NEW VERSION FOR GCM15F, BASED ON GWDFL13:
C           *                          - WORK ARRAYS ARE NOW LOCAL TO THE
C           *                            ROUTINE, IE NOT PASSED IN WITH
C           *                            "WRK" POINTERS.


C
C     * MAR 31/04 - J.SCINOCCA.  - DO ROOF DRAG FOR MAM AS WELL IF IN
C     *                            SPINUP PERIOD ("IDAMP" PASSED IN
C     *                            AND USED, WITH QVCRIT ADJUSTED LOWER).
C     *                          - LOWER GPHIL TO 0.5 FOR MAM RUNS.
C     * FEB  4/00 - J.SCINOCCA.  PREVIOUS VERSION GWDFL12.
C     * MAY 17/99 - M. LAZARE/   PREVIOUS VERSION GWDFL11 FOR GCM12.
C     *             J. SCINOCCA. 
C
C NEW VERSION OF GRAVITY WAVE DRAG PARAMETERIZATION.  THIS SCHEME REQUIRES 
C 6 QUANTITIES WHICH CHARACTERIZE THE SUB-GRID TOPOGRAPHY IN EACH
C GAUSSIAN GRID CELL.  GIVEN A REPRESENTATION OF THE SUB-GRID TOPOGRAPHIC
C VARIATIONS IN A GG CELL H(X,Y) WE CAN DEFINE:
C
C
C ENV   -  RMS DEVIATIONS  (UNITS M) 
C                
C             ENV = ( < H(X,Y)2 > )^(1/2) 
C 
C             WHERE THE ANGLE BRACKETS DENOTE AN AVERAGE OVER THE GG CELL
C             <F(X,Y)> = (1/LX)*(1/LY)INT{INT{ F(X,Y) DX}DY} 
C
C PSIM_EAST  -  DIRECTION OF PRINCIPLE AXIS (MINOR AXIS) realative to due east
C                    -PI/2 <= PSIM_EAST <= PI/2
C
C             PSIM_EAST = (1/2) ARCTAN(M/L), 
C
C             WHERE  M = < H_X H_Y >, AND, L = (1/2)[ <(H_X)2> - <(H_Y)2>]
C
C GAM   -  ANISOTROPY (OR ECCENTRICITY OF REPRESENTATIVE ELLIPSE)
C                      0 <= GAM <= 1.
C
C             GAM = ( < (H_Y')2 > / < (H_X')2 > )^(1/2)
C
C             WHERE (X',Y') ARE NEW COORDINATES ROTATED SO THAT
C             X' POINTS IN ALONG PRINCIPLE DIRECTION (MINOR AXIS)
C
C ALPHA -  FRACTION OF ENV2 WHICH IS 2D 
C                   0 <= ALPHA <= 1
C
C             ALPHA = INT{ INT{ H(X',Y') DY'}2 DX'}    INT{ H(X')2 DX'}
C                     -----------------------------  =  -----------------
C                     INT{ INT{ H(X',Y')2 DY'} DX'}          ENV2
C
C          FOR ALPHA NON-ZERO, H(X') REPRESENTS THE PROFILE OF THE 2D RIDGE
C
C DELTA -  DIRECTION ALONG PRINCIPLE AXIS OF INCREASING SLOPE OF 2D TOPOGRAPHY
C                   -1 <= DELTA <= 1.
C
C             DELTA = INT{ H_X'(X') DX'}
C                     ------------------
C                     INT{|H_X'(X')|DX'}
C
C             DELTA =  0 => GG 2D RIDGE IS LOCALIZED IN CELL 
C                   <= 0 => GG 2D RIDGE PREDOMINANTLY
C                              DOWNSLOPE ALONG PRINCIPLE AXIS
C                   >= 0 => GG 2D RIDGE PREDOMINANTLY
C                              UPSLOPE ALONG PRINCIPLE AXIS
C
C SIGX  - RMS VALUE OF SLOPE ALONG PRINCIPLE AXIS
C
C             SIGX = < (H_X')2 >^(1/2)
C
C
C THE NEW SCHEME HAS THE FOLLOWING THREE PRINCIPLE COMPONENTS:
C
C
C 1) VERTICAL FLUX OF MOMENTUM ASSOCIATED WITH FREELY PROPAGATING GRAVITY
C    WAVES CARRIED BY TWO WAVES WITH AMPLITUDES AND DIRECTIONS THAT ARE
C    DETERMINED FROM LINEAR THEORY OF FLOW OVER A REPRESENTATIVE ELLIPTICAL
C    BARRIER (PHILLIPS, 1984).  THE ELLIPTICAL BARRIER WHICH REPRESENTS THE
C    SUB-GRID TOPOGRAPHY IN EACH GAUSSIAN GRID CELL IS SPECIFIED BY:
C    ENV, GAM, PSIM_EAST, AND SIGX
C
C    THE MOMENTUM FLUX CARRIED BY EACH WAVE IS DETERMINED BY PERFORMING
C    THE AZIMUTHAL INTEGRAL OF THE PRESSURE-DRAG DENSITY OVER THE HALF-SPACE
C    TO EACH SIDE OF THE WIND DIRECTION.  THE HORIZONTAL MOMENTUM FLUX 
C    PER UNIT AREA IN EACH WAVE AT THE SURFACE IS GIVEN BY VECTORS OF THE FORM:
C
C    (TAUP_X,TAUP_Y) = - RHO*N*U*GPHIL*B*(2*ENV)**2*(VDP,VTP)/L**2
C
C    (TAUM_X,TAUM_X) = - RHO*N*U*GPHIL*B*(2*ENV)**2*(VDM,VTM)/L**2
C
C     WHERE:
C            RHO - DENSITY
C              N - BV FREQ
C              U - MAGNITUDE OF WIND
C          GPHIL - (EQ 6.6 PHILLIPS, 1984) IS A CONSTANT WHICH DEPENDS
C                   ON MOUNTAIN SHARPNESS. 
C                  IT HAS A VALUE IN THE RANGE 1<= GPHIL <=2.
C              B - MAJOR AXIS LENGTH OF REPRESENTATIVE ELLIPSE
C          2*ENV - EFFECTIVE HEIGHT OF TOPOGRAPHY
C      (VDP,VTP) 
C      (VDM,VTM) - VECTORS REPRESENTING THE AZIMUTHAL INTEGRAL
C                  OF PRESSURE FORCE.  THESE REQUIRE THE NUMERICAL
C                  INTEGRATION OF INCOMPLETE ELLIPTIC INTEGRALS.
C                  RATHER THAN CALCULATE THESE VECTORS WHEN REQUIRED,
C                  WE FILL IN A LOOK_UP TABLE FOR (VDP,VTP) AS A
C                  FUNCTION OF GAMMA AND ANGLE BETWEEN THE WIND AND THE
C                  MINOR AXIS OF THE REPRESENTATIVE ELLIPSE AT THE 
C                  BEGINNING OF EACH RUN.  THE VALUE OF (VDM,VTM)
C                  IS ALSO OBTAINED FROM THE LOOK_UP TABLE AND IS
C                  DERIVED FROM SYMMETRY.
C           L**2 - AREA OF GAUSSIAN GRID BOX
C
C    THE ABOVE FORMULAS FOR (TAUP_X,TAUP_Y) AND (TAUM_X,TAUM_X) MUST BE
C    MULTIPLIED BY L**2/4AB - THE NUMBER OF ELLIPSES IN THE GG CELL.
C    THIS IS AS IN LOTT AND MILLER (1997).  THIS ALSO REDUCES THE STRESS
C    FOR THE CASE OF B>>L (IE 2D RIDGE).  THIS REQUIRES THE EVALUATION OF
C    1/4A FOR THE NORMALIZATION OF (TAUP_X,TAUP_Y) AND (TAUM_X,TAUM_X).
C    WE TAKE A ~ ENV/SIGX
C
C    GIVEN (TAUP_X,TAUP_Y) AND (TAUM_X,TAUM_X), WE DERIVE TWO MONOCHROMATIC
C    WAVES IN THE SAME DIRECTIONS THAT CARRY THE SAME MOMENTUM FLUX.
C    GIVEN THESE, WE CAN PROPAGATE THEM UPWARDS IN THE SAME MANNER DESCRIBED
C    BY THE OLD SCHEME.  FOR THE CASE WHEN THE TWO WAVES POINT IN THE
C    SAME DIRECTION (2D RIDGE) IT IS DESIRED THAT THE SCHEME GIVE THE SAME
C    RESULT AS IF 1 WAVE CARRIED ALL THE MOMENTUM FLUX.  AS A CONSEQUENCE,
C    WE DO NOT ALLOW THE TWO WAVES TO INTERACT AS THEY PROPAGATE UPWARDS
C    AND WE MAKE THEIR INITIAL AMPLITUDES THE SAME AND EQUAL TO THE
C    LAUNCH AMPLITUDE U*FR/N.
C
C    THE FORM OF THE TWO MONOCHROMATIC WAVES ARE:
C
C    (TAUP_X,TAUP_Y) = -(MUP/2)*RHO*N*UP*(HP)**2*(XP,YP)
C
C    (TAUM_X,TAUM_X) = -(MUM/2)*RHO*N*UM*(HM)**2*(XM,YM)
C
C    WHERE:
C
C    (XP,YP),(XM,YM) - UNIT VECTORS IN (VDP,VTP) AND (VDM,VTM) DIRECTIONS RESP.
C        UP = INNER PRODUCT OF LAUNCH WIND (U,V) AND (XP,YP)
C        UM = INNER PRODUCT OF LAUNCH WIND (U,V) AND (XM,YM)
C
C    AND MUP,MUM,HP,HM DETERMINED FROM THE ABOVE TWO CONDITIONS PLUS
C
C        N*HP/UP = FR,       N*HM/UM=FR
C
C    RESULTING IN:
C
C    MUP/2 = (GPHIL/4A)*|(VDP,VTP)|/COS(THP)**3
C    MUM/2 = (GPHIL/4A)*|(VDM,VTM)|/COS(THM)**3
C       HP = H*COS(THP)
C       HM = H*COS(THM)
C
C    WHERE:
C
C    THP - ANGLE BETWEEN WIND AND P-WAVE
C    THM - ANGLE BETWEEN WIND AND M-WAVE
C      H - LAUNCH HEIGHT OF TOPOG
C                  
C
C 2) LOW-LEVEL BLOCKING DRAG ASSOCIATED WITH SUB-GRID TOPOGRAPHIC BARRIER.  
C    THIS PARAMETERIZATION REQUIRES ENV, GAM, PSIM_EAST, ALPHA, DELTA AND SIGX. 
C    THE SURFACE STRESS OF THE LOW-LEVEL BLOCKING IS EVALUATED AS FORM DRAG:
C
C               - 0.5*RHO*CDBLK*HBEFF*U*|U| *(B/L2)*(L2/4AB)
C
C    WHERE RHO   - DENSITY
C          CDBLK - DRAG COEFFICIENT O(1)
C          HBEFF - HEIGHT OF BLOCKING (AT MOST HBEFF = 2*ENV - FCRIT*UAV/NAV)
C          U     - MAGNITUDE OF WIND 
C          A,B   - MAJOR,MINOR AXIS OF ELLIPSE 
C          L2   - AREA OF GAUSSIAN GRID CELL
C          L2/4AB - NUMBER OF ELLIPSES IN GG BOX
C                    NOTE: A ~ ENV/SIGX  SO THAT NORMALIZATION FACTOR
C                    (B/L2)*(L2/4AB) = 1/4A ~ 0.25*SIGX/ENV
C
C    THIS DRAG WILL IS APPLIED IN THREE STEPS:
C    1) THE COMPONENT IF THE FLOW ALONG THE PRINCIPLE AXIS (IE NORMAL
C       TO THE POTENTIAL RIDGE) IS FIRST BLOCKED BY THAT PORTION OF 
C       THE SUB-GRID TOPOGRAPHY WHICH IS 2D IN NATURE USING THE BLOCKING
C       COEFFICIENT CDBLK2D.  THE 2D BLOCKED REGION IS SPECIFIED BY
C       HB2D = ALPHA*(2*ENV - FCRIT*UAV/NAV)
C
C    2) THE COMPONENT OF THE FLOW ALONG THE PRINCIPLE AXIS 
C       IS NEXT BLOCKED OVER THE DEPTH HB2D TO HB2D+HB3D
C       WHERE HB3D = (1-ALPHA)*(2*ENV - FCRIT*UAV/NAV) WITH BLOCKING 
C       COEFFICIENT CDBLK3D
C
C    3) THE COMPONENT OF THE FLOW NORMAL TO THE PRINCIPLE AXIS IS
C       THEN BLOCKED UP TO A HEIGHT HB3D WITH BLOCKING COEFFICIENT CDBLK3D
C
C 3) ENHANCED WAVE-DRAG ASSOCIATED WITH LOW-LEVEL WAVE BREAKING.
C    BETWEEN VALUES OF FCRIT <= FR <= 3*FCRIT CONSTANT DRAG IS
C    APPLIED BELOW 3/4 OF A HYDROSTATIC VERTICAL WAVELENGTH.
C    THERE ARE TWO CONDITIONS THAT MUST BE SATISFIED IN ORDER
C    FOR THIS PARAMETERIZATION TO TURN ON.
C
C    A) FCRIT <= FR <= 3*FCRIT
C    B) WIND MUST NOT BE DIRECTED UPSLOPE FOR CASE WHEN ALPHA AND DELTA 
C       ARE LARGE (>0.75)
C
C    SURFACE STRESS ASSOCIATED WITH DOWNSLOPE FLOW GIVEN BY:
C
C       TAU_D = ( 1 + EN(ALPHA,F))* TAU_FREE  - TAUC
C
C       WHERE:
C       TAU_FREE - IS LINEAR-THEORY SURFACE STRESS FOR FLOW OVER ELLIPSE
C       TAUC     - IS CRITICAL LINEAR-THEORY SURFACE STRESS FOR F=FCRIT
C       EN       - IS ENHANCEMENT FACTOR OF DRAG GIVING A MAXIMUM VALUE
C                  OF 1 FOR 3D (ALPHA=0) TOPOGRAPHY AND A MAXIMUM VALUE
C                  OF 3 FOR 2D (ALPHA=1) TOPOGRAPHY.
C       
C       EN IS A FUNCTION OF F SUCH THAT TAU_D IS ZERO AT F=FCRIT AND 
C       F=3FCRIT AND HAVING ITS MAXIMUM VALUE AT 1.5FCRIT
C
C
C THREE SWITCHES ARE USED:  IWFLG, IBFLG, IDFLG
C
C   * IWFLG - FLAG WHICH DETERMINES THE MODE OF OPERATION 
C             IWFLG < 0  => MOMENTUM FLUX DETERMINED BY JUST ENV FIELD
C                           ASSUMING ISOTROPY OF SUB-GRID TOPOGRAPHY.
C                           THIS IS OLD WAY OF DOING GWD (NOT TESTED)
C             IWFLG >= 0 => MOMENTUM FLUX DETERMINED BY FLOW OVER
C                           ELLIPSE WHICH IS REPRESENTATIVE OF SUB-GRID
C                           TOPOGRAPHY.  ELLIPSE SPECIFIED BY GAM,ENV
C                           PSIM_EAST
C             =-1 => 1 WAVE COUNTERFLOW DIRECTION
C             = 0 => 1 WAVE COUNTERFLOW DIRECTION
C             = 1 => 1 WAVE, DIRECTION DETERMINED BY ORIENTATION OF
C                      REPRESENTATIVE ELLIPSE AND DIRECTION OF WIND
C             = 2 => 2 WAVES, DIRECTIONS DETERMINED BY ORIENTATION OF
C                      REPRESENTATIVE ELLIPSE AND DIRECTION OF WIND.
C                      THE VECTORIAL SUM OF THE MOMENTUM FLUX FROM
C                      THE TWO WAVES IS EQUIVALENT TO THE IWFLG=1 CASE
C
C   * IBFLG - FLAG WHICH TURNS ON LOW-LEVEL BLOCKING
C             
C             IBFLG = 1 -> ON
C                   = 0 -> OFF
C
C   * IDFLG - FLAG WHICH TURNS ON LOW-LEVEL BREAKING
C             
C             IDFLG = 1 -> ON
C                   = 0 -> OFF
C
C OTHER CHANGES:
C
C - DEPOSITION OF MOMENTUM IS RETURNED TO METHOD DESCRIBED IN MCFARLANE 1987
C - REYNOLDS STRESS IS 3-POINT AVERAGED IN THE VERTICAL TO SMOOTH OUT 
C   DISCONTINUITY WHERE DRAG SWITCHES ON IN THE VERTICAL.  THIS ALSO
C   ALLOWS US NOT TO USE THE IMPLICIT SCHEME IN THE DEPOSITION.
C - THE EFFICIENCY FACTOR TAUFAC IS ONLY USED WHEN THE OPTION IWFLG=-1 IS
C   USED.  OTHERWISE, THE NEW FACTOR OUT FRONT OF THE DRAG IS:
C                   (L**2/4AB) * (B*GPHIL/L**2) = GPHIL/4A ~ GPHIL*0.25*SIGX/ENV
C
C       WHERE GPHIL (EQ 6.6 PHILLIPS, 1984) IS A CONSTANT WHICH DEPENDS
C             ON MOUNTAIN SHARPNESS IT HAS A VALUE IN THE RANGE 1<= GPHIL <=2.
C - TOPOGRAPHIC WAVES ARE NO LONGER LAUNCHED FROM A SINGLE SPECIFIED 
C   REFERENCE LEVEL (LREF=ILEV-1).  INSTEAD, THE VELOCITY AND 
C   BUOYANCY FREQUENCY ARE AVERAGED IN THE VERTICAL UP TO 2*ENV 
C   TO GIVE UAV,VAV, AND NAV. THESE VALUES ARE THEN USED TO DETERMINE 
C   THE NONDIMENSIONAL HEIGHT OF THE SUB-GRID TOPOGRAPHY:
C                     2*ENV*NAV/SQRT(UAV**2+VAV**2).
C   IF THE NONDIMENSIONAL TOPOGRAPHY IS SUB CRITICAL (< FCRIT) THEN WAVES ARE
C   LAUNCHED AT ILEV-1 AND THE REYNOLDS STRESS CAN BEGIN TO CHANGE VALUE ON ILEV-2.
C   IF THE NONDIMENSIONAL HEIGHT OF THE TOPOGRAPHY IS GREATER THAN THE FCRIT 
C   THEN WE LAUNCH WITH AMPLITUDE FCRIT*SQRT(UAV**2+VAV**2)/NAV FROM THE
C   NEXT MODEL LEVEL BELOW THE HEIGHT 2*ENV-FCRIT*SQRT(UAV**2+VAV**2)/NAV.
C   THE REYNOLDS STRESS IS ALLOWED TO THEN BEGIN CHANGING VALUE ON THE NEXT 
C   HIGHEST MODEL LEVEL (ALTHOUGH THE 3-POINT VERTICAL SMOOTHER LOWERS THE
C   TURN-ON OF THE DRAG BY ONE LEVEL).

      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)

      PARAMETER(IWFLG=2,IBFLG=1,IDFLG=1)
  
C     * I/O ARRAYS PASSED FROM MAIN.
  
      REAL   UROW(ILG,ILEV),        VROW(ILG,ILEV),      TSG(ILG,ILEV),
     1       TSGB(ILG,ILEV),        SGJ(ILG,ILEV),   
     2       SGBJ(ILG,ILEV),      DSGJ(ILG,ILEV),   PRESSG(ILG)


      REAL FTOX(ILG,ILEV),FTOY(ILG,ILEV)
      REAL TDOX(ILG,ILEV),TDOY(ILG,ILEV)

c     keep these as work arrays.  If we want to output them later 
c     they can be added to the call with pack and row arrays for output
      REAL TNDP_X(ILG,ILEV),TNDP_Y(ILG,ILEV)
      REAL TNDM_X(ILG,ILEV),TNDM_Y(ILG,ILEV)
      REAL TAUP_X(ILG,ILEV),TAUP_Y(ILG,ILEV)
      REAL TAUM_X(ILG,ILEV),TAUM_Y(ILG,ILEV)

         

C     * ARRAYS WHICH SPECIFY SUB-GRID SCALE TOPOGRAPHY

C     ENV  - SUB-GRID STANDARD DEVIATION FIELD
C     GAM  - ASPECT RATIO OF REPRESENTATIVE ELLIPSE 0<= GAM <=1
C     PSIM_EAST - ANGLE THAT MINOR AXIS IS ROTATED FROM X AXIS
C            (IE DUE EAST) -PI/2<= PSIM_EAST <= PI/2
C     ALPHA - FRACTION OF ENV**2 THAT IS ASSOCIATED WITH 2D RIDGE
C             WHICH IS NORMAL TO PSIM_EAST DIRECTION 0<= ALPHA <=1
C     DELTA - NORMALIZED NET SLOPE OF 2D RIDGE  -1 <= DELTA <= 1
C     SIGX - RMS SLOPE ALONG PRINCIPLE AXIS

      REAL ENV(ILG),GAM(ILG),PSIM(ILG),ALPHA(ILG),DELTA(ILG),
     1     SIGX(ILG),PSIM_EAST(ILG)
      REAL PHIE(ILG)

      REAL     UTENDGW(ILG,ILEV),  VTENDGW(ILG,ILEV),
     1         UTENDDW(ILG,ILEV),  VTENDDW(ILG,ILEV),
     2        UTENDLLB(ILG,ILEV), VTENDLLB(ILG,ILEV),
     3            UM(ILG,ILEV),     VM(ILG,ILEV)

C     * WORK ARRAYS TO CALCULATE COMPONENT TENDENCIES (ONLY USED IF
C     * IFIZGWD=1).

      REAL UGROW(ILG,ILEV),   VGROW(ILG,ILEV),
     1     USROW(ILG,ILEV),   VSROW(ILG,ILEV),
     2     UFROW(ILG,ILEV),   VFROW(ILG,ILEV)

C     * WORK ARRAYS PASSED IN FROM MAIN

C     * VERTICAL ELEVATION ARRAYS

      REAL PB(ILG,ILEV),    ZB(ILG,ILEV)
      REAL ZM(ILG,ILEV) 
      REAL  P(ILG,ILEV),   GOPS(ILG) 

C     * ARRAYS USED FOR GENERAL TWO-WAVE SPECTRUM

      REAL UVDP(ILG),  UVDM(ILG),  UVTP(ILG),  UVTM(ILG)
      REAL   KP(ILG),    KM(ILG),  HPSQ(ILG),  HMSQ(ILG)
      REAL PHI,PSI,PIH,PWLO,PWHI,GWLO,GWHI
      REAL TAUP(ILG,ILEV),TAUM(ILG,ILEV),WRK1(ILG)

C     * UPSI - COMPONENT OF (UAV,VAV) WIND ALONG PRINCIPLE AXIS
C     * TAUW - STRESS DUE TO LOW-LEVEL WAVEBREAKING
C     * TAUB - STRESS DUE TO LOW-LEVEL BLOCKING
C     * THIS STRESS IS MADE TO GO TO ZERO AT 3/4 OF A HYDROSTATIC
C     *  VERTICAL WAVELENGTH, WHICH OCCURS AT LEVEL IZW IN THE MODEL

      REAL UPSI(ILG),VPSI(ILG),TAUW(ILG,ILEV),TAUB(ILG,ILEV)
      
C     * ARRAYS TO DETERMINE NONDIMENSIONAL MOUNTAIN HEIGHT
C     * EQUAL TO NOUAV*HWAV
C     * SD2 - 2 TIMES STANDARD DEVIATION ENV
C     * IZT - INDEX OF NEXT MOMENTUM BASE LEVEL ABOVE 2.*ENV
C     * IZB - INDEX OF NEXT MOMENTUM BASE LEVEL ABOVE (SD2-HWAV)
C     * NAV   - VERTICAL AVERAGE OF   N   FROM (0.,ZT)
C     * UAV   - VERTICAL AVERAGE OF   U   FROM (0.,ZT)
C     * VAV   - VERTICAL AVERAGE OF   V   FROM (0.,ZT)
C     * RAV   - VERTICAL AVERAGE OF  RHO  FROM (0.,ZT)
C     * NOUAV - RATIO OF (NAV//SQRT(UAV**2+VAV**2))
C     * HWAV  - THAT PORTION OF SD2 THAT LAUNCHES WAVES
C     * HBLK  - THAT PORTION OF SD2 (SD2-HWAV) THAT LIES BELOW LAUNCHING
C     *         LEVEL.  USEFUL FOR BLOCKING. 

      REAL NOUAV(ILG),NAV(ILG),UAV(ILG),VAV(ILG),RAV(ILG),SD2(ILG),
     1     HWAV(ILG),HBLK(ILG)

      INTEGER IZT(ILG),IZB(ILG)

C     * FLAG DRAG(IL)=1. FOR GWD SCHEME ACTIVE AT LONGITUDE IL AND
C      *     DRAG(IL)=-1.  FOR GWD SCHEME NOT ACTIVE AT LONGITUDE IL

      REAL DRAG(ILG)
  
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),  
     1     UB(ILG),           VB(ILG),
     2    AMPBSQ(ILG),   HITESQ(ILG),DENFACP(ILG),UMAV(ILG)

      REAL PSIM_SUM(ILG)

C     * LOOK-UP TABLES DP AND TP

C       DP,TP - ALONG-FLOW AND NORMAL-FLOW COMPONENTS OF FIRST GRAVITY WAVE 
C             LAUNCHED.  THE VALUES DEPEND ON THE ASPECT RATIO OF THE
C             REPRESENTATIVE ELLIPSE (GAM) AND THE ANGLE BETWEEN THE 
C             MINOR AXIS OF THE ELLIPSE AND THE DIRECTION OF THE FLOW 
C             (UAV,VAV).
C     * THE SECOND WAVE MAY BE OBTAINED FROM THE SAME LOOK-UP TABLES DP,TP
C       WITH DP(PSI,GAM) -> DP(-PSI,GAM) AND
C          TP(PSI,GAM) -> -TP(-PSI,GAM)
C       NPSI,NGAM - NUMBER OF POINTS USED IN DP,TP
C       DPSI,DGAM - INCREMENTS IN EACH DIRECTION 
C       PSIST,GAMST - START VALUES OF PSI AND GAM

      COMMON /GWDTABL/PI,NPSI,NGAM,DPSI,DGAM,PSIST,GAMST,
     1                DP(100,100),TP(100,100)
      REAL*8 PI,DPSI,DGAM,PSIST,GAMST,DP,TP

C     
C     * EXTRA WORK ARRAYS USED FOR OPTIMIZATION
C       IZW,UDU,UDV,TV1,TV2,TV3,TV4,TV5,TV6,TF1

      INTEGER IZW(ILG)
      REAL TV1(ILG),TV2(ILG),TV3(ILG),TV4(ILG),TV5(ILG),TV6(ILG)
      REAL TF1(ILG,ILEV)
      REAL UDU(ILG),UDV(ILG)
C
      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     WORK ARRAYS FOR WINDS ON BASE LEVELS
      REAL U(ILG,ILEV),V(ILG,ILEV),BVM(ILG,ILEV)

      REAL CPSIM(ILG),SPSIM(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 OCCURRED.
C     * DMPSCAL  = DAMPING TIME FOR GW DRAG IN SECONDS. 
C     * TAUFAC   = 1/(LENGTH SCALE).
C     * HMIN     = MIMINUM ENVELOPE HEIGHT REQUIRED TO PRODUCE GW RAG. 
C     * V0       = VALUE OF WIND THAT APPROXIMATES ZERO.
C     * RADEA    = RADIUS OF EARTH (M)
C     * GPHIL    = MOUNTAIN SHARPNESS NUMBER
C     * CDBLK    = DRAG COEFFICIENT FOR LOW-LEVEL BLOCKING
  
      DATA    VMIN  /    2.0 /, V0       / 1.E-10 /, 
     1        TAUFAC/ 1.0E-5 /, HMIN     /   3.0 /,
     2        ZERO  /    0.0 /, UN       /    1.0 /,
     3        GRAV  /9.80616 /, DMPSCAL  / 6.5E+6 /, RADEA /6.37E+6/,
     4        GPHIL /1.0/, CDBLK2D /1.0/, CDBLK3D /0.5/, 
     5        ENHANCE /3.0/
C
      DATA    BVFMIN /1.E-10/
C-----------------------------------------------------------------------


C     MC FOR OGW -> FOR CONSERVATION MAKE IMCONS=1 OTHERWISE IMCONS=0
      IMCONS=1
      OMEGA2=2.*7.2722E-5

      PIH=PI/2

C     ORIENTATION OF ELLIPSE RELATIVE TO DUE EAST, PSIM_EAST, AND ANGLE OF 
C     DUE EAST RELATIVE TO THE LOCAL ZONAL DIRECTION, PHIE, USED TO SET THE
C     THE ORIENTATION ANGLE OF THE ELLIPSE RELATIVE TO THE LOCAL ZONAL DIRECTION
      PSIM_SUM=PSIM_EAST+PHIE
      DO IL=IL1,IL2
         IF(PSIM_SUM(IL).GT.PIH) THEN
            PSIM(IL)=PSIM_SUM(IL)-PI
         ELSEIF(PSIM_SUM(IL).LT.-PIH) THEN
            PSIM(IL)=PSIM_SUM(IL)+PI
         ELSE
            PSIM(IL)=PSIM_SUM(IL)
         ENDIF
      ENDDO

      

C     * CHOOSE TUNING CONSTANTS BASED ON RESOLUTION BEING RUN:
C     * ITUNE=1 => T47  ;   ITUNE=2 => T63
C     * QVCRIT SET TO LOW VALUE IN SPINUP PERIOD FOR MAM RUNS (IDAMP=1).
C
      IF(IDAMP.EQ.1)    THEN
        QVCRIT=50.
      ELSE
        IF(ITUNE.EQ.1)    THEN
          QVCRIT=100.
        ELSE
          QVCRIT=90.
        ENDIF
      ENDIF
C


C
C     * OVERRIDE GCM13 SETTING OF GPHIL FOR MAM RUNS
C
      IF(LMAM) GPHIL=1.0
C
      ILEVM = ILEV - 1
      DPG=DPSI*DGAM
      GRSQ=GRAV**2
C
      FCRIT=0.22
C     FCRIT=0.375
      FCSQ=FCRIT**2
      BRKLEV=3*PIH


C     GET U,V ON MOMENTUM BASE LEVELS

      DO 55 L=1,ILEV-1  
        DO 50 I=IL1,IL2                                        
           FCT1=LOG(SGJ (I,L)/SGBJ(I,L))
           FCT2=LOG(SGBJ(I,L)/SGJ(I,L+1))
           FCT3=1./LOG(SGJ (I,L)/SGJ(I,L+1))
           U(I,L) = (UROW(I,L+1)*FCT1+UROW(I,L)*FCT2)*FCT3
           V(I,L) = (VROW(I,L+1)*FCT1+VROW(I,L)*FCT2)*FCT3
 50     CONTINUE
 55   CONTINUE                                              
      DO 60 I=IL1,IL2
         U(I,ILEV)=U(I,ILEV-1)
         V(I,ILEV)=V(I,ILEV-1)
 60   CONTINUE



C
C     * TEST FOR APPLICATION OF GWD SCHEME 
C     * DRAG REFERENCES THE POINTS WHERE GW CALCULATIONS WILL BE DONE,
C     * THAT IS (A- REMOVE TEST FOR OVER LAND IE GC=-1), B- IF WE
C     * ASK FOR IT AND C- IF ENVELOPE HEIGHT >= HMIN )
C     * DRAG=1. FOR POINTS THAT SATISFY THE ABOVE, ELSE DRAG=-1.
C

      LENGW = 0

      IF(ENVELOP)   THEN
        DO 110 I=IL1,IL2
          IF(ENV(I).GE.HMIN)          THEN
            LENGW   = LENGW + 1
            DRAG(I) = 1.
          ELSE
            DRAG(I) = -1.
          ENDIF
          SD2(I)=2.*ENV(I)
          HBLK(I)=0.
          TV1(I)=-1.
          IZT(I)=ILEV-1
          IZB(I)=ILEV-2
  110   CONTINUE
      ENDIF


C
C     * INITIALIZE NECESSARY ARRAYS.
C

      UTENDGW = ZERO
      VTENDGW = ZERO
      UTENDDW = ZERO
      VTENDDW = ZERO
      UTENDLLB = ZERO
      VTENDLLB = ZERO
      UM = UROW
      VM = VROW


      IF (LENGW.EQ.0)                                        GO TO 500

C      
C     * GET VERTICAL ELEVATIONS Z IN METERS
C
      DO 130 IL=IL1,IL2
        PB(IL,1) = SGBJ(IL,1) * PRESSG(IL)
        P(IL,1)  = SGJ(IL,1)  * PRESSG(IL)
        GOPS(IL) = GRAV/PRESSG(IL)
  130 CONTINUE

      DO 140 L=2,ILEV
      DO 140 IL=IL1,IL2
        PB(IL,L)   = SGBJ(IL,L)   * PRESSG(IL)
        P(IL,L)    = SGJ(IL,L)    * PRESSG(IL)
  140 CONTINUE

      ROG=RGAS/GRAV
      DO 150 IL=IL1,IL2
       ZB(IL,ILEV) = 0.
  150 CONTINUE

      DO 160 L=ILEVM,1,-1
      DO 160 IL=IL1,IL2
       ZB(IL,L) = ZB(IL,L+1)+ ROG*TSG(IL,L+1)*LOG(PB(IL,L+1)/PB(IL,L))
  160 CONTINUE

C     GET MID-LEVEL HEIGHTS
      DO 162 IL=IL1,IL2                                               
        ZM(IL,ILEV) = ROG*TSGB(IL,ILEV)*LOG(PRESSG(IL)/P(IL,ILEV))    
 162  CONTINUE                                                        
                                                                      
      DO 165 L=ILEVM,1,-1                                             
      DO 165 IL=IL1,IL2                                               
        ZM(IL,L) = ZM(IL,L+1) + ROG*TSGB(IL,L)*LOG(P(IL,L+1)/P(IL,L))
 165  CONTINUE                                                        
      

C
C     * GET LEVEL INDEX OF BASE (IZT) JUST BELOW SD2
C       IZT WILL BE AT LEAST ILEV-1 IN ORDER THAT WE AVERAGE U,V,BVFREQ ETC.
C       OVER AT LEAST THE FIRST TWO LEVELS
C

      DO 170 L=ILEV-2,2,-1
      DO 170 IL=IL1,IL2
        IF(DRAG(IL).GT.0..AND.ZM(IL,L+1).LT.SD2(IL)) THEN
           IZT(IL)=L
        ENDIF
  170 CONTINUE
         


C      
C     * CALCULATE  B V FREQUENCY 
C
 
      DO 180 L=2,ILEV
      DO 180 I=IL1,IL2
        IF(DRAG(I).GT.0.) THEN
          DTDS=-SGJ(I,L)/TSG(I,L)*(TSGB(I,L)-TSGB(I,L-1))/DSGJ(I,L)
     1          +RGOCP
          BVFREQ(I,L)=MAX(GRSQ*(DTDS/(RGAS*TSG(I,L))),BVFMIN)
          BVM(I,L)=SQRT(BVFREQ(I,L))
        ENDIF
  180 CONTINUE
      DO 190 I=IL1,IL2
        IF(DRAG(I).GT.0.) THEN
          BVM(I,1)=BVM(I,2)
        ENDIF
  190 CONTINUE



C     GET BVFREQ ON MOMENTUM BASE LEVELS

      DO 193 L=1,ILEV-1  
        DO 192 I=IL1,IL2                                        
           IF(DRAG(I).GT.0.) THEN
              FCT1=LOG(SGJ (I,L)/SGBJ(I,L))
              FCT2=LOG(SGBJ(I,L)/SGJ(I,L+1))
              FCT3=1./LOG(SGJ (I,L)/SGJ(I,L+1))
              BVFREQ(I,L) = (BVM(I,L+1)*FCT1+BVM(I,L)*FCT2)*FCT3
           ENDIF
 192    CONTINUE
 193  CONTINUE                                              
      DO 195 I=IL1,IL2
         IF(DRAG(I).GT.0.) THEN
            BVFREQ(I,ILEV)=BVFREQ(I,ILEV-1)
         ENDIF
 195  CONTINUE

      DO 200 L=2,ILEV
      DO 200 I=IL1,IL2
        IF(DRAG(I).GT.0.) THEN
         RATIO=5.*LOG(SGBJ(I,L)/SGBJ(I,L-1))
         BVFREQ(I,L) = (BVFREQ(I,L-1) + RATIO*BVFREQ(I,L))
     1                   /(1.+RATIO)
        ENDIF
  200 CONTINUE



      DO 202 L=2,ILEV                                         
      DO 202 I=IL1,IL2                                        
        IF(DRAG(I).GT.0.) THEN                                
          RATIO=5.*LOG(SGJ(I,L)/SGJ(I,L-1))                   
          BVM(I,L) = (BVM(I,L-1) + RATIO*BVM(I,L))   
     1                   /(1.+RATIO)                          
        ENDIF                                                 
 202  CONTINUE                                                

C
C     * EVALUATE AVERAGES UP TO HEIGHT SD2
C
      DO 205 IL=IL1,IL2
        IF(DRAG(IL).GT.0.) THEN
           DZ=ZM(IL,ILEV)                  
           TV1(IL)=DZ                           
           UAV(IL)=DZ*U(IL,ILEV)                   
           VAV(IL)=DZ*V(IL,ILEV)                   
           NAV(IL)=DZ*BVFREQ(IL,ILEV)              
           RAV(IL)=DZ*PB(IL,ILEV)/(RGAS*TSGB(IL,ILEV))
        ENDIF
  205 CONTINUE

      DO 210 L=ILEV-1,2,-1
      DO 210 IL=IL1,IL2
        IF(DRAG(IL).GT.0..AND.L.GE.IZT(IL)) THEN
          DZ=ZM(IL,L)-ZM(IL,L+1)
          TV1(IL)=TV1(IL)+DZ
          UAV(IL)=UAV(IL)+DZ*U(IL,L)
          VAV(IL)=VAV(IL)+DZ*V(IL,L)
          NAV(IL)=NAV(IL)+DZ*BVFREQ(IL,L)
          RAV(IL)=RAV(IL)+DZ*PB(IL,L)/(RGAS*TSGB(IL,L))
        ENDIF
  210 CONTINUE

      DO 215 IL=IL1,IL2
        IF(DRAG(IL).GT.0.) THEN
          ZAV=TV1(IL)
          UAV(IL)=UAV(IL)/ZAV
          VAV(IL)=VAV(IL)/ZAV
          NAV(IL)=NAV(IL)/ZAV
          RAV(IL)=RAV(IL)/ZAV
          NOUAV(IL)=NAV(IL)/SQRT(UAV(IL)**2+VAV(IL)**2)
        ENDIF
  215 CONTINUE


C      
C     * PERFORM LAST CHECK ON MAGNITUDE OF WIND TO SET VALUE OF DRAG.
C

      DO 220 IL=IL1,IL2
        IF(DRAG(IL).GT.0.) THEN
          UMAV(IL)=SQRT(UAV(IL)**2+VAV(IL)**2)
          IF(UMAV(IL).LT.VMIN) THEN
            LENGW=LENGW-1
            DRAG(IL)=-1.
          ELSE
            HWAV(IL)=MIN(SD2(IL),FCRIT/NOUAV(IL))
            HBLK(IL)=SD2(IL)-HWAV(IL)
          ENDIF
        ENDIF
  220 CONTINUE


      IF (LENGW.EQ.0)                                        GO TO 500

C
C     * DETERMINE NEXT MODEL LEVEL ABOVE HEIGHT OF POTENTIAL 
C     * BLOCKED LAYER. NOTE THAT FOR CASES IN WHICH THE LAUNCH INVERSE
C     * FROUDE NUMBER SD2*NOUAV IS LESS THAT FCRIT (IE NO BLOCKING) IZB=ILEV-2. 
C     * THAT IS, THE VERTICAL FLUX OF HORIZONTAL MOMENTUM IS CONSTANT IN AT 
C     * LEAST THE FIRST TWO LEVELS.
C

      DO 230 L=ILEV-2,2,-1 
      DO 230 IL=IL1,IL2
        IF(DRAG(IL).GT.0..AND.
     1     ZM(IL,L+1).LT.HBLK(IL)) THEN
           IZB(IL)=L-1
        ENDIF
 230  CONTINUE



C
C     * GET UNIT VECTOR IN DIRECTION OF (UAV,VAV)
C     * USE LOOK-UP TABLES (DP,TP), WHICH MAY BE USED TO DERIVE THE
C       COMPONENTS OF PRESSURE FORCE OF EACH WAVE ON TOPOGRAPHY 
C       (NEGATIVE OF REYNOLD'S STRESS) IN THE DIRECTION OF THE WIND (DP) 
C       AND NORMAL TO THE WIND (TP).
C       DP AND TP ARE A FUNCTION OF ECCENTRICITY, GAM, AND ORIENTATION OF 
C       TOPOGRAPHY WITH RESPECT TO THE DIRECTION OF THE WIND. 
C       DP AND TP ARE INTERPOLATED TO DETERMINE THE COMPONENTS (VDP,VTP) OF 
C       THE FIRST WAVE AND THE COMPONENTS (VDM,VTM) OF THE SECOND WAVE.
C      

      DO 250 IL=IL1,IL2
        IF(DRAG(IL).GT.0.) THEN
          TFNC=SIGX(IL)*GPHIL*0.25/ENV(IL)
          UMAG=UMAV(IL)
          UB(IL)=UAV(IL)/UMAG
          VB(IL)=VAV(IL)/UMAG
          PHI=ATAN2(VB(IL),UB(IL))
          CPSIM(IL)=COS(PSIM(IL))
          SPSIM(IL)=SIN(PSIM(IL))


C         
C     * ANGLE OF WIND FROM X-AXIS PHI, AND
C       ANGLE OF THE PRINCIPLE AXIS OF THE TOPOGRAPHY FROM X-AXIS PSIM
C       DETERMINE THE RELATIVE ANGLE OF SEPARATION PSI=PSIM-PHI
C

          PSI=PSIM(IL)-PHI
          PSISVE=PSI

C
C     * PSI MUST BE IN THE RANGE (-PI/2, PI/2)
C
          IF(PSI.GT.PIH)  PSI=PSI-PI
          IF(PSI.LT.-PIH) PSI=PSI+PI

C
C     * NOW FIND INDICES, IPSI AND IGAM,OF DP AND TP LOOK-UP 
C       TABLES THAT ARE JUST BELOW THE VALUES OF PSI AND GAM
C

          IPSI=(PSI-PSIST)/DPSI+1
          IF(IPSI.EQ.NPSI) IPSI=NPSI-1
          PWLO=PSI-((IPSI-1)*DPSI+PSIST)
          PWHI=DPSI-PWLO
          IGAM=(GAM(IL)-GAMST)/DGAM+1
          IF(IGAM.EQ.NGAM) IGAM=NGAM-1
          GWLO=GAM(IL)-((IGAM-1)*DGAM+GAMST)
          GWHI=DGAM-GWLO

C
C     * INTERPOLATE TO OBTAIN VDPU AND VTPU
C

          VDPU=( GWHI* (PWHI*DP(IGAM  ,IPSI)+PWLO*DP(IGAM  ,IPSI+1))
     1          +GWLO* (PWHI*DP(IGAM+1,IPSI)+PWLO*DP(IGAM+1,IPSI+1)) )
     2             /DPG
          VTPU=( GWHI* (PWHI*TP(IGAM,IPSI  )+PWLO*TP(IGAM,IPSI+1  ))
     1          +GWLO* (PWHI*TP(IGAM+1,IPSI)+PWLO*TP(IGAM+1,IPSI+1)) )
     2             /DPG

C
C     * USE SYMMETRY TO OBTAIN VALUES FOR SECOND WAVE VDMU AND VTMU
C       FROM SAME LOOK-UP TABLE
C

          PSI=-PSI
          IPSI=(PSI-PSIST)/DPSI+1
          IF(IPSI.EQ.NPSI) IPSI=NPSI-1
          PWLO=PSI-((IPSI-1)*DPSI+PSIST)
          PWHI=DPSI-PWLO
          VDMU=( GWHI* (PWHI*DP(IGAM,IPSI  )+PWLO*DP(IGAM,IPSI+1  ))
     1          +GWLO* (PWHI*DP(IGAM+1,IPSI)+PWLO*DP(IGAM+1,IPSI+1)) )
     2             /DPG
          VTMU=-( GWHI* (PWHI*TP(IGAM,IPSI  )+PWLO*TP(IGAM,IPSI+1  ))
     1           +GWLO* (PWHI*TP(IGAM+1,IPSI)+PWLO*TP(IGAM+1,IPSI+1)) )
     2              /DPG

C
C     * NOTE: (VDPU,VTPU), (VDMU,VTMU) ARE VECTORS DECOMPOSED INTO A DIRECTION
C            PARALLEL TO THE WIND, VDP,VDM      AND A DIRECTION 
C            NORMAL   TO THE WIND, VTP,VTM     (USING RIGHT-HAND RULE)
C     * (VDPU,VTPU) AND (VDMU,VTMU) MUST BE ROTATED THROUGH THE ANGLE PHI,
C       WHICH THE WIND MAKES WITH THE MODEL X-AXIS
C       THIS RESULTS IN (VDP,VTP) AND (VDM,VTM) WHICH ARE WRT THE MODEL 
C       (X,Y) (IE ZONAL AND MERIDIONAL DIRECTIONS)
C

          VDP= (VDPU*UB(IL)-VTPU*VB(IL))
          VTP= (VDPU*VB(IL)+VTPU*UB(IL))
          VDM= (VDMU*UB(IL)-VTMU*VB(IL))
          VTM= (VDMU*VB(IL)+VTMU*UB(IL))
C
C     * NOW MUST REPRESENT PRESSURE DRAG OF EACH WAVE FROM REPRESENTATIVE 
C       ELLIPSE, AT ANGLES THETP AND THETM SAY, BY SINGLE WAVE FORMULAS:
C
C       KP*RAV*NAV*UAV*COS(THETP)*HP**2  
C       KM*RAV*NAV*UAV*COS(THETM)*HM**2 
C
C       THAT IS, WE MUST DETERMINE KP AND HP, AND, KM AND HM 
C       IN ORDER THAT THE SINGLE-WAVE FORMULAS PROVIDE THE SAME DRAG
C       AS (VDP,VTP) AND (VDM,VTM)
C
C     * NOTE: THAT (UVDP,UVTP) AND (UVDM,UVTM) ARE THE COMPONENTS OF UNIT
C       VECTORS OF THE FIRST AND SECOND WAVE IN THE MODEL COORDINATES (X,Y)
C     * FOR THE CASE OF ONE WAVE, ONLY THE FIRST WAVE WILL BE USED (IE
C       (UVDP,UVTP) ETC.).
C
          IF(IWFLG.EQ.2) THEN
            RMP=SQRT(VDP**2+VTP**2)
            UVDP(IL)=VDP/RMP
            UVTP(IL)=VTP/RMP
            RMM=SQRT(VDM**2+VTM**2)
            UVDM(IL)=VDM/RMM
            UVTM(IL)=VTM/RMM
            CSP=VDPU/RMP
            CSM=VDMU/RMM
            KP(IL)=TFNC*RMP/CSP**3
            KM(IL)=TFNC*RMM/CSM**3
            HPSQ(IL)=(HWAV(IL)*CSP)**2
            HMSQ(IL)=(HWAV(IL)*CSM)**2
          ELSEIF(IWFLG.EQ.1) THEN
            VDP=VDP+VDM
            VTP=VTP+VTM
            RMP=SQRT(VDP**2+VTP**2)
            UVDP(IL)=VDP/RMP
            UVTP(IL)=VTP/RMP
            CSP=(VDPU+VDMU)/RMP
            KP(IL)=TFNC*RMP/CSP**3
            HPSQ(IL)=(HWAV(IL)*CSP)**2
          ELSEIF(IWFLG.EQ.0) THEN
            VDP=VDP+VDM
            VTP=VTP+VTM
            RMP=SQRT(VDP**2+VTP**2)
            UVDP(IL)=UB(IL)
            UVTP(IL)=VB(IL)
            KP(IL)=TFNC*RMP
            HPSQ(IL)=(HWAV(IL))**2
          ELSEIF(IWFLG.EQ.-1) THEN
            UVDP(IL)=UB(IL)
            UVTP(IL)=VB(IL)
            KP(IL)=TAUFAC
            HPSQ(IL)=(HWAV(IL))**2
          ENDIF  
        ENDIF
  250 CONTINUE





C
C     * FROM HERE ON WE TREAT ONLY SINGLE WAVE MOMENTUM FLUXES
C     * DETERMINE PROFILES OF (NEGATIVE) MOMENTUM FLUX FOR EACH WAVE 
C       TAUP AND TAUM
C

      DO 255 IL=IL1,IL2
        IF(DRAG(IL).GT.0.) THEN
          UMAG=UAV(IL)*UVDP(IL)+VAV(IL)*UVTP(IL)
          BFAC=RAV(IL)*NAV(IL)*UMAG
          AMPBSQ(IL)=BFAC
          HITESQ(IL)=HPSQ(IL)
          TV1(IL)=KP(IL)*GOPS(IL)*BFAC*HPSQ(IL)
        ENDIF
  255 CONTINUE


      DO 270 L=ILEV,1,-1
      DO 270 IL=IL1,IL2
        IF(DRAG(IL).GT.0.) THEN
          IF(L.GT.IZB(IL)) THEN
            TF1(IL,L)=TV1(IL)
          ELSE
            UMAG=MAX(U(IL,L)*UVDP(IL)+V(IL,L)*UVTP(IL),V0)
            BFAC=UMAG*BVFREQ(IL,L)*PB(IL,L)/(RGAS*TSGB(IL,L))
            HSQMAX=FCSQ*(UMAG/BVFREQ(IL,L))**2
            RATIO=AMPBSQ(IL)/BFAC
            HSQ=MIN(RATIO*HITESQ(IL),HSQMAX)
            HITESQ(IL)=HSQ
            AMPBSQ(IL)=BFAC
            TF1(IL,L)=KP(IL)*GOPS(IL)*BFAC*HSQ
         ENDIF
        ENDIF
  270 CONTINUE

      DO 300 L=1,ILEV                                           
      DO 300 IL=IL1,IL2                                            
        IF(DRAG(IL).GT.0.) THEN                                    
          TAUP(IL,L)=TF1(IL,L)
        ENDIF                                                      
  300 CONTINUE
           
 745  FORMAT(I4,3E24.17)

C
C     * IF TWO-WAVE SCHEME THEN DETERMINE STRESS FOR SECOND WAVE
C
      IF(IWFLG.EQ.2) THEN

      DO 305 IL=IL1,IL2
        IF(DRAG(IL).GT.0.) THEN
          UMAG=UAV(IL)*UVDM(IL)+VAV(IL)*UVTM(IL)
          BFAC=RAV(IL)*NAV(IL)*UMAG
          AMPBSQ(IL)=BFAC
          HITESQ(IL)=HMSQ(IL)
          TV1(IL)=KM(IL)*GOPS(IL)*BFAC*HMSQ(IL)
        ENDIF
  305 CONTINUE

        DO 320 L=ILEV,1,-1
        DO 320 IL=IL1,IL2
          IF(DRAG(IL).GT.0.) THEN
            IF(L.GT.IZB(IL)) THEN
              TF1(IL,L)=TV1(IL)
            ELSE
              UMAG=MAX(U(IL,L)*UVDM(IL)+V(IL,L)*UVTM(IL),V0)
              BFAC=UMAG*BVFREQ(IL,L)*PB(IL,L)/(RGAS*TSGB(IL,L))
              HSQMAX=FCSQ*(UMAG/BVFREQ(IL,L))**2
              RATIO=AMPBSQ(IL)/BFAC
              HSQ=MIN(RATIO*HITESQ(IL),HSQMAX)
              HITESQ(IL)=HSQ
              AMPBSQ(IL)=BFAC
              TF1(IL,L)=KM(IL)*GOPS(IL)*BFAC*HSQ
            ENDIF
          ENDIF
  320   CONTINUE

        DO 350 L=1,ILEV                                           
        DO 350 IL=IL1,IL2                                           
          IF(DRAG(IL).GT.0.) THEN                                   
            TAUM(IL,L)=TF1(IL,L)
          ENDIF                                                     
  350   CONTINUE                                                    

      ENDIF          



C
C     * CALCULATE GW TENDENCIES (TOP AND BOTTOM LAYERS FIRST).
C     * BOTTOM LAYER KEEPS INITIALIZED VALUE OF ZERO.
C     * APPLY GW DRAG ON (UM,VM) WORK ARRAYS TO BE PASSED TO VRTDFS.
C
  

C     BELOW IS THE DEFAULT FOR IMCONS=0 (I.E NO FORCING IN UPPER LAYER)
      DO 360 I=IL1,IL2
        IF(DRAG(I).GT.0.) THEN
          UTENDGW(I,1)=0.   
          VTENDGW(I,1)=0.          
        ENDIF 
  360 CONTINUE
  

      TNDP_X(:,1)=0.0
      TNDP_Y(:,1)=0.0
      TNDM_X(:,1)=0.0
      TNDM_Y(:,1)=0.0

      DO 370 L=2,ILEV
      DO 370 I=IL1,IL2
         IF(DRAG(I).GT.0.) THEN
            DFLX=TAUP(I,L)-TAUP(I,L-1)                   
            DENOM=SGBJ(I,L)-SGBJ(I,L-1)                       
            TENDFAC=DFLX/DENOM                
            UTENDGW(I,L) = -TENDFAC*UVDP(I)               
            VTENDGW(I,L) = -TENDFAC*UVTP(I)  
            TNDP_X(I,L)=UTENDGW(I,L)
            TNDP_Y(I,L)=VTENDGW(I,L)
            IF(IWFLG.EQ.2) THEN                           
               DFLX=TAUM(I,L)-TAUM(I,L-1)                   
               TENDFAC=DFLX/DENOM                
               UTENDGW(I,L) =UTENDGW(I,L) -TENDFAC*UVDM(I) 
               VTENDGW(I,L) =VTENDGW(I,L) -TENDFAC*UVTM(I) 
               TNDM_X(I,L)=-TENDFAC*UVDM(I)
               TNDM_Y(I,L)=-TENDFAC*UVTM(I)
            ENDIF                                         
            UM(I,L) = UM(I,L) + ZTMST*UTENDGW(I,L)
            VM(I,L) = VM(I,L) + ZTMST*VTENDGW(I,L)
         ENDIF 
 370  CONTINUE




C     DO UPPER LAYER FOR IMCONS=1
C     IN THIS CASE THE FLUX IS ZERO ABOVE THE HIGHEST MID LEVEL
C     AND IT IS ASSUMMED THAT THE HORIZONTAL WAVENUMBER IS CONSTANT
C     THE WKB COMPONENT IS ZERO ON THE UPPER LEVEL

      IF(IMCONS.EQ.1) THEN
      L=1
      DO 375 I=IL1,IL2
         IF(DRAG(I).GT.0.) THEN
            DFLX=TAUP(I,L)                   
            DENOM=DSGJ(I,L)                       
            TENDFAC=DFLX/DENOM                
            UTENDGW(I,L) = -TENDFAC*UVDP(I)               
            VTENDGW(I,L) = -TENDFAC*UVTP(I)  
            TNDP_X(I,L)=UTENDGW(I,L)
            TNDP_Y(I,L)=VTENDGW(I,L)
              
            IF(IWFLG.EQ.2) THEN                           
               DFLX=TAUM(I,L)                   
               TENDFAC=DFLX/DENOM                
               UTENDGW(I,L) =UTENDGW(I,L) -TENDFAC*UVDM(I) 
               VTENDGW(I,L) =VTENDGW(I,L) -TENDFAC*UVTM(I) 
               TNDM_X(I,L)=-TENDFAC*UVDM(I)
               TNDM_Y(I,L)=-TENDFAC*UVTM(I)
            ENDIF                                         
            UM(I,L) = UM(I,L) + ZTMST*UTENDGW(I,L)
            VM(I,L) = VM(I,L) + ZTMST*VTENDGW(I,L)
         ENDIF 
 375  CONTINUE
      ENDIF

C
C     * IF REQUIRED, GET WIND ALONG PRINCIPLE 
C       AXIS (UPSI) USED FOR BOTH LOW-LEVEL BLOCKING AND 
C       LOW-LEVEL BREAKING
C     * ALSO GET WIND NORMAL TO PA AND STORE IN GOPS
C       SINCE THIS ARRAY IS NOT REQUIRED ANY MORE
C

      IF(IBFLG+IDFLG.GT.0.AND.IWFLG.GE.0) THEN
         DO 380 IL=IL1,IL2
            IF(DRAG(IL).GT.0.) THEN
               UPSI(IL)=UAV(IL)*CPSIM(IL)+VAV(IL)*SPSIM(IL)
               VPSI(IL)=VAV(IL)*CPSIM(IL)-UAV(IL)*SPSIM(IL)
            ENDIF
 380     CONTINUE
         DO 382 L=ILEV-2,2,-1                                  
         DO 382 IL=IL1,IL2                                     
            IF(DRAG(IL).GT.0..AND.ZB(IL,L).LT.SD2(IL)) THEN     
               IZT(IL)=L                                        
            ENDIF                                               
 382     CONTINUE                                              
      ENDIF



C  
C     * LOW-LEVEL WAVE-BREAKING PARAMETERIZATION
C
      IF(IDFLG.EQ.0.OR.IWFLG.LT.0) GOTO 450

C
C     * CHECK FOR DRAG ENHANCEMENT DUE TO WAVEBREAKING NEAR
C       THE SURFACE
C     * CHECK FOR UPSLOPE FLOW ON 2D TOPOGRAPHY
C       AND DO NOT ALLOW LOW-LEVEL BREAKING
C


      DO 390 IL=IL1,IL2
        IF(DRAG(IL).GT.0.) THEN
          FR=NOUAV(IL)*SD2(IL)
          FS=FR/FCRIT
          IF(FS.GT.1..AND.FS.LT.3..AND.
     1      (ALPHA(IL).LT.0.75.OR.
     2      ABS(DELTA(IL)).LT.0.75.OR.
     3      UPSI(IL)*DELTA(IL).LT.0.)) THEN

            ENH=(1.+2.*ALPHA(IL))*ENHANCE/3.
            IF(FS.LE.1.5) THEN
              ZETA=2.*ENH*(FS-1.0)
            ELSE
              ZETA=((3.-FS)/1.5)**2*(1.+ENH-(1./1.5)**2)
     1              +(1./FS**2)-1.
            ENDIF
            FRAT=(1.+ZETA)*FR**2-FCRIT**2
            IF(IWFLG.EQ.2) THEN
              UPP=(UAV(IL)*UVDP(IL)+VAV(IL)*UVTP(IL))
              UMM=(UAV(IL)*UVDM(IL)+VAV(IL)*UVTM(IL))
              DRGX=(UPP**3*KP(IL)*UVDP(IL)+
     1              UMM**3*KM(IL)*UVDM(IL))
              DRGY=(UPP**3*KP(IL)*UVTP(IL)+
     1              UMM**3*KM(IL)*UVTM(IL))

              DRG=SQRT(DRGX**2+DRGY**2)
              UDU(IL)=DRGX/DRG
              UDV(IL)=DRGY/DRG
            ELSE
              UPP=(UAV(IL)*UVDP(IL)+VAV(IL)*UVTP(IL))
              DRGX=(UPP**3*KP(IL)*UVDP(IL))
              DRGY=(UPP**3*KP(IL)*UVTP(IL))
              DRG=SQRT(DRGX**2+DRGY**2)
              UDU(IL)=UVDP(IL)
              UDV(IL)=UVTP(IL)
            ENDIF

C
C       GET THE MAGNITUDE OF THE WIND IN THE DIRECTION
C       OF THE PRESSURE FORCE
C
            UMAG=UAV(IL)*UDU(IL)+VAV(IL)*UDV(IL)
            SGNM=SIGN(UN,UMAG)
            DRGD=DRG*(RAV(IL)/NAV(IL))*FRAT

C
C     * FIND 3/4 LAMBDA_Z WAVE-BREAKING LEVEL
C     * FIRST ASSUME THAT U AND N ARE CONST UP TO IZT
C     * SINCE FR/FCRIT -> [1.,3.] BREAKING LEVEL MUST
C       BE HIGHER THAN IZT. 
C 
C     * STORE PHASE UP TO IZT LEVEL IN TV2
C       STORE SIGN OF INNER PRODUCT OF WIND AND PRESSURE FORCE 
C             DIRECTION IN TV3
C       STORE MAGNITUDE OF PRESSURE FORCE IN TV4
C       STORE LAST LEVEL BELOW BREAKING LEVEL IN IZW

            TV2(IL)=(ZB(IL,IZT(IL))-ZB(IL,ILEV))*NAV(IL)/UMAG
            TV3(IL)=SGNM
            TV4(IL)=DRGD
            IZW(IL)=IZT(IL)

C
C     * TV1 IS FLAG TO INDICATE WHETHER THIS PARAMETERIZATION IS ACTIVE
C
            TV1(IL)=FS
          ELSE
            TV1(IL)=0.
          ENDIF
        ENDIF
  390 CONTINUE




C
C     * NOW PERFORM WKB INTEGRATION IN VERTICAL
C     * STORE MODEL LEVEL JUST BELOW BREAKING LEVEL IN TV2(IL)
C
      DO 400 L=ILEV-1,2,-1
      DO 400 IL=IL1,IL2
        IF(DRAG(IL).GT.0..AND.
     1    TV1(IL).GT.0.5) THEN
          IF(L.LT.IZT(IL).AND.IZW(IL).EQ.L+1) THEN
            UMAG=MAX(TV3(IL)
     1      *(UROW(IL,L)*UDU(IL)+VROW(IL,L)*UDV(IL)),V0)
            DZP=(ZB(IL,L)-ZB(IL,L+1))*BVM(IL,L)/UMAG
            TV2(IL)=TV2(IL)+DZP
            IF(TV2(IL).LT.BRKLEV) THEN
              IZW(IL)=L
            ELSE
              ZBREAK=ZB(IL,IZW(IL))
              ZDIFF=BRKLEV-(TV2(IL)-DZP)
              ZEX=ZDIFF*UMAG/BVM(IL,L)
              ZBREAK=ZBREAK+ZEX
              TV4(IL)=TV4(IL)/ZBREAK
C
C     * STORE ADDITIONAL FRACTIONAL DRAG FORCE IN TV5
C
              TV5(IL)=TV4(IL)*ZEX/
     1                (ZB(IL,IZW(IL)-1)-ZB(IL,IZW(IL)))
            ENDIF
          ENDIF
        ENDIF
  400 CONTINUE



      DO 410 IL=IL1,IL2
        IF(DRAG(IL).GT.0..AND.
     1    TV1(IL).GT.0.5) THEN

          L=IZW(IL)
          RFCT=RGAS*TSG(IL,L)/P(IL,L)
          UMAG=MAX(TV3(IL)*(UROW(IL,L)*UDU(IL)+VROW(IL,L)*UDV(IL)),V0)
          TND_MAX=UMAG/(ZTMST)
          TND=-RFCT*TV5(IL)
          IF(ABS(TND).GT.TND_MAX) TND=-TV3(IL)*TND_MAX
          UTENDDW(IL,L)=TND*UDU(IL)
          VTENDDW(IL,L)=TND*UDV(IL)
          UM(IL,L) = UM(IL,L) + ZTMST*UTENDDW(IL,L)
          VM(IL,L) = VM(IL,L) + ZTMST*VTENDDW(IL,L)
        ENDIF
  410 CONTINUE


      DO 420 L=ILEV,1,-1
      DO 420 IL=IL1,IL2
        IF(DRAG(IL).GT.0..AND.
     1    TV1(IL).GT.0.5) THEN

          IF(L.GT.IZW(IL)) THEN
            RFCT=RGAS*TSG(IL,L)/P(IL,L)
            UMAG=MAX(TV3(IL)*(UROW(IL,L)*UDU(IL)+VROW(IL,L)*UDV(IL)),V0)
            TND_MAX=UMAG/(ZTMST)
            TND=-RFCT*TV4(IL)
            IF(ABS(TND).GT.TND_MAX) TND=-TV3(IL)*TND_MAX
            UTENDDW(IL,L)=TND*UDU(IL)
            VTENDDW(IL,L)=TND*UDV(IL)
            UM(IL,L) = UM(IL,L) + ZTMST*UTENDDW(IL,L)
            VM(IL,L) = VM(IL,L) + ZTMST*VTENDDW(IL,L)
          ENDIF
        ENDIF
  420 CONTINUE


  450 CONTINUE



C
C     * LOW-LEVEL BLOCKING PARAMETERIZATION
C

      IF(IBFLG.EQ.0.OR.IWFLG.LT.0) GOTO 500
C
C     * DO LOW-LEVEL BLOCKING
C
      DO 460 IL=IL1,IL2
        IF(DRAG(IL).GT.0..AND.
     1     HBLK(IL).GT.HMIN) THEN

          HB2D=ALPHA(IL)*HBLK(IL)
          HB3D=(1.-ALPHA(IL))*HBLK(IL)

C
C     * GET EFFECTIVE 2D BLOCKING HEIGHT
C
          HBEFF=HB2D

C
C     * IF 2D PORTION OF VARIANCE HAS MOST OF SLOPE IN WRONG
C       DIRECTION (IE UPSI AND DELTA OPPOSITE SIGN) THEN REDUCE
C       2D BLOCKING HEIGHT APPROPRIATELY
C
          IF(UPSI(IL)*DELTA(IL).LT.0) THEN
            HBEFF=HBEFF*(1.-ABS(DELTA(IL)))
          ENDIF

C
C     * STORE MAGNITUDE OF BLOCKING FORCE ALONG P.A. IN TV4 (NOTE THIS 
C         MUST BE MULTIPLIED BY THE BULK DRAG COEFF CDBLK2D OR CDBLK3D)
C       STORE MAGNITUDE OF BLOCKING FORCE NORMAL TO P.A. IN TV6 
C         (NOTE THIS IS TOTAL MAGNITUDE UNLIKE TV4)
C       STORE DIRECTION VECTOR OF PRINCIPLE AXIS IN UDU AND UDV
C       STORE SIGN OF INNER PRODUCT OF WIND AND PRINCIPLE AXIS 
C             DIRECTION IN TV3
C       STORE SIGN OF INNER PRODUCT OF WIND AND MAJOR AXIS
C             DIRECTION IN TV5
C       STORE EFFECTIVE 2D BLOCKING HEIGHT IN TV1
C       STORE 3D BLOCKING HEIGHT IN TV2
C       STORE NEXT MODEL BASE LEVEL BELOW 2D BLOCKING HEIGHT IN IZW
C       STORE NEXT MODEL BASE LEVEL BELOW ELEVATED 3D BLOCKING HEIGHT
C             IN IZT
C       STORE NEXT MODEL BASE LEVEL BELOW 3D BLOCKING HEIGHT IN IZB

          TFNC=SIGX(IL)*0.25/ENV(IL)
          TV4(IL)=0.5*TFNC*RAV(IL)*UPSI(IL)*ABS(UPSI(IL))
          TV6(IL)=0.5*TFNC*CDBLK3D*RAV(IL)*VPSI(IL)*ABS(VPSI(IL))
          TV3(IL)=SIGN(UN,UPSI(IL))
          TV5(IL)=SIGN(UN,VPSI(IL))
          TV1(IL)=HBEFF
          TV2(IL)=HB3D
          IZW(IL)=ILEV
          IZT(IL)=ILEV
          IZB(IL)=ILEV
        ELSE
          TV1(IL)=-1.
        ENDIF
  460 CONTINUE

C
C     * FIND HEIGHT OF 2D, ELEVATED 3D, AND 3D BLOCKING LEVEL
C
      DO 470 L=ILEV-1,2,-1
      DO 470 IL=IL1,IL2
        IF(DRAG(IL).GT.0..AND.
     1     TV1(IL).GT.-0.5) THEN
     
          HBEFF=TV1(IL)
          HB3D=TV2(IL)
          HBTOT=HBEFF+HB3D 
          IF(ZB(IL,L).LE.HBEFF) THEN
            IZW(IL)=L
          ENDIF
          IF(ZB(IL,L).LE.HBTOT) THEN
            IZT(IL)=L
          ENDIF
          IF(ZB(IL,L).LE.HB3D) THEN
            IZB(IL)=L
          ENDIF
        ENDIF
  470 CONTINUE


C
C     * FIRST DO END-POINT TENDENCIES OF BLOCKING LAYERS
C
      DO 480 IL=IL1,IL2
        IF(DRAG(IL).GT.0..AND.
     1     TV1(IL).GT.-0.5) THEN

          LW=IZW(IL)
          LT=IZT(IL)
          LB=IZB(IL)
          HBEFF=TV1(IL)
          HB3D=TV2(IL)
          HBTOT=HBEFF+HB3D
          SPSI=TV3(IL)

C
C     * EVALUATE BLOCKING TENDENCY AT TOP OF 2D LAYER
C

          DRGB=CDBLK2D*TV4(IL)
          DRGTOP=DRGB*(HBEFF-ZB(IL,LW))/
     1          (ZB(IL,LW-1)-ZB(IL,LW))
          RFCT=RGAS*TSG(IL,LW)/P(IL,LW)

          UMAG=MAX(SPSI
     1         *(UROW(IL,LW)*CPSIM(IL)+VROW(IL,LW)*SPSIM(IL)),V0)


C
C     * SET THE MAXIMUM ALOWABLE TENDECY BASED ON UMAG
C       NOTE UMAG, AND THEREFORE TND_MAX, IS POSITIVE.
C
          TND_MAX=UMAG/(ZTMST)
C
C     * TND IS THE TENDENCY WE WOULD LIKE TO IMPOSE (WITH
C       THE CORRECT SIGN)
C
          TND=-RFCT*DRGTOP
C
C     * CHECK TO SEE IF IT EXCEEDS THE MAXIMUM.  IF IT
C       DOES THEN SET TND TO TND_MAX WITH CORRECT SIGN
C
          IF(ABS(TND).GT.TND_MAX) TND=-SPSI*TND_MAX
          UTENDLLB(IL,LW)=TND*CPSIM(IL)+UTENDLLB(IL,LW)
          VTENDLLB(IL,LW)=TND*SPSIM(IL)+VTENDLLB(IL,LW)

C
C     * EVALUATE BLOCKING TENDENCY AT TOP AND BOTTOM
C       OF ELEVATED 3D LAYER
C

C
C     * IF NOT EVEN ONE LEVEL EXTRA THEN LEAVE OUT 3D BLOCKING
C
          IF(LW.GT.LT) THEN

C
C     * APPLY 3D BLOCKING TENDECIES IN DIRECTION ALONG PA
C

            DRGB=CDBLK3D*TV4(IL)
            DRGBOT=DRGB*(ZB(IL,LW-1)-HBEFF)/
     1              (ZB(IL,LW-1)-ZB(IL,LW))
            DRGTOP=DRGB*(HBTOT-ZB(IL,LT))/
     1              (ZB(IL,LT-1)-ZB(IL,LT))

            RFCT=RGAS*TSG(IL,LW)/P(IL,LW)
            UMAG=MAX(SPSI
     1      *(UROW(IL,LW)*CPSIM(IL)+VROW(IL,LW)*SPSIM(IL)),V0)
            TND_MAX=UMAG/(ZTMST)
            TND=-RFCT*DRGBOT
            IF(ABS(TND).GT.TND_MAX) TND=-SPSI*TND_MAX
            UTENDLLB(IL,LW)=TND*CPSIM(IL)+UTENDLLB(IL,LW)
            VTENDLLB(IL,LW)=TND*SPSIM(IL)+VTENDLLB(IL,LW)

            RFCT=RGAS*TSG(IL,LT)/P(IL,LT)
            UMAG=MAX(SPSI
     1      *(UROW(IL,LT)*CPSIM(IL)+VROW(IL,LT)*SPSIM(IL)),V0)
            TND_MAX=UMAG/(ZTMST)
            TND=-RFCT*DRGTOP
            IF(ABS(TND).GT.TND_MAX) TND=-SPSI*TND_MAX
            UTENDLLB(IL,LT)=TND*CPSIM(IL)+UTENDLLB(IL,LT)
            VTENDLLB(IL,LT)=TND*SPSIM(IL)+VTENDLLB(IL,LT)

          ENDIF


C
C     * APPLY 3D BLOCKING TENDECIES IN DIRECTION NORMAL TO PA
C

          SPSI=TV5(IL)        
          DRGB=TV6(IL)
          DRGTOP=DRGB*(HBEFF-ZB(IL,LB))/
     1                (ZB(IL,LB-1)-ZB(IL,LB))

          RFCT=RGAS*TSG(IL,LB)/P(IL,LB)

C
C     * ROTATE UNIT VECTORS TO POINT NORMAL TO PRINCIPLE AXIS
C
          UMAG=MAX(SPSI
     1    *(-UROW(IL,LB)*SPSIM(IL)+VROW(IL,LB)*CPSIM(IL)),V0)
          TND_MAX=UMAG/(ZTMST)
          TND=-RFCT*DRGTOP
          IF(ABS(TND).GT.TND_MAX) TND=-SPSI*TND_MAX
          UTENDLLB(IL,LB)=-TND*SPSIM(IL)+UTENDLLB(IL,LB)
          VTENDLLB(IL,LB)=TND*CPSIM(IL)+VTENDLLB(IL,LB)

        ENDIF

  480 CONTINUE



C
C    
        DO 490 L=ILEV,1,-1
        DO 490 IL=IL1,IL2
        IF(DRAG(IL).GT.0..AND.
     1     TV1(IL).GT.-0.5) THEN

          LW=IZW(IL)
          LT=IZT(IL)
          LB=IZB(IL)

C
C     * 2D BLOCKING LAYER TENDENCIES
C
          IF(L.GT.LW) THEN
            RFCT=RGAS*TSG(IL,L)/P(IL,L)
            SPSI=TV3(IL)
            UMAG=MAX(SPSI
     1      *(UROW(IL,L)*CPSIM(IL)+VROW(IL,L)*SPSIM(IL)),V0)
            DRGB=CDBLK2D*TV4(IL)
            TND_MAX=UMAG/(ZTMST)
            TND=-RFCT*DRGB
            IF(ABS(TND).GT.TND_MAX) TND=-SPSI*TND_MAX
            UTENDLLB(IL,L)=TND*CPSIM(IL)+UTENDLLB(IL,L)
            VTENDLLB(IL,L)=TND*SPSIM(IL)+VTENDLLB(IL,L)
          ENDIF

C
C     * ELEVATED 3D BLOCKING LAYER TENDENCIES
C
          IF(L.LT.LW.AND.L.GT.LT) THEN
            RFCT=RGAS*TSG(IL,L)/P(IL,L)
            SPSI=TV3(IL)
            UMAG=MAX(SPSI
     1      *(UROW(IL,L)*CPSIM(IL)+VROW(IL,L)*SPSIM(IL)),V0)
            DRGB=CDBLK3D*TV4(IL)
            TND_MAX=UMAG/(ZTMST)
            TND=-RFCT*DRGB
            IF(ABS(TND).GT.TND_MAX) TND=-SPSI*TND_MAX
            UTENDLLB(IL,L)=TND*CPSIM(IL)+UTENDLLB(IL,L)
            VTENDLLB(IL,L)=TND*SPSIM(IL)+VTENDLLB(IL,L)
          ENDIF

C
C     * 3D BLOCKING LAYER TENDENCIES
C
          IF(L.GT.LB) THEN
            RFCT=RGAS*TSG(IL,L)/P(IL,L)
            SPSI=TV5(IL)
            DRGB=TV6(IL)
            UMAG=MAX(SPSI
     1      *(-UROW(IL,L)*SPSIM(IL)+VROW(IL,L)*CPSIM(IL)),V0)
            TND_MAX=UMAG/(ZTMST)
            TND=-RFCT*DRGB
            IF(ABS(TND).GT.TND_MAX) TND=-SPSI*TND_MAX
            UTENDLLB(IL,L)=-TND*SPSIM(IL)+UTENDLLB(IL,L)
            VTENDLLB(IL,L)= TND*CPSIM(IL)+VTENDLLB(IL,L)
          ENDIF

C
C     * UPDATE UM AND VM
C
          UM(IL,L) = UM(IL,L) + ZTMST*UTENDLLB(IL,L)
          VM(IL,L) = VM(IL,L) + ZTMST*VTENDLLB(IL,L)


        ENDIF

  490 CONTINUE


  500 CONTINUE


      DO 510 L=1,ILEV
      DO 510 I=IL1,IL2
         UTENDGW(I,L)=UTENDGW(I,L)+UTENDDW(I,L)+UTENDLLB(I,L)
         VTENDGW(I,L)=VTENDGW(I,L)+VTENDDW(I,L)+VTENDLLB(I,L)
  510 CONTINUE

      TDOX=0.0
      TDOY=0.0
      FTOX=0.0
      FTOY=0.0
      TAUP_X=0.0
      TAUP_Y=0.0
      TAUM_X=0.0
      TAUM_Y=0.0
      DO L=1,ILEV
         DO IL=IL1,IL2
            IF(DRAG(IL).GT.0.) THEN
            RFT=P(IL,L)/(RGAS*TSG(IL,L))
            UNITPX=UVDP(IL)
            UNITPY=UVTP(IL)
            UNITMX=UVDM(IL)
            UNITMY=UVTM(IL)
            TAUP_X(IL,L)=-(TAUP(IL,L))*UNITPX/GOPS(IL)
            TAUP_Y(IL,L)=-(TAUP(IL,L))*UNITPY/GOPS(IL)
            TAUM_X(IL,L)=-(TAUM(IL,L))*UNITMX/GOPS(IL)
            TAUM_Y(IL,L)=-(TAUM(IL,L))*UNITMY/GOPS(IL)
            FTOX(IL,L)=TAUM_X(IL,L) + TAUP_X(IL,L)
            FTOY(IL,L)=TAUM_Y(IL,L) + TAUP_Y(IL,L)

             TDOX(IL,L)=TNDP_X(IL,L)+TNDM_X(IL,L)
             TDOY(IL,L)=TNDP_Y(IL,L)+TNDM_Y(IL,L)
            ENDIF
         ENDDO
      ENDDO


 1550 FORMAT(I4,5E24.17)
 1551 FORMAT(I4,2E24.17)


 124  FORMAT(I4,4E24.17)





      IF(IFIZGWD.NE.0)                                      THEN

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

        DO 520 L=1,ILEV
        DO 520 I=IL1,IL2
          UGROW(I,L)=UM(I,L)
          VGROW(I,L)=VM(I,L)
  520   CONTINUE
      ENDIF

C     * USROW,UVROW CONTAIN NON-ORO GWD TENDENCY.
C     * APPLY NON-ORO TENDENCY HERE. 

      DO 525 L=1,ILEV  
      DO 525 I=IL1,IL2  
          UM(I,L) = UM(I,L) + ZTMST*UNOTND(I,L)
          VM(I,L) = VM(I,L) + ZTMST*VNOTND(I,L)
          UTENDGW(I,L)=UTENDGW(I,L)+UNOTND(I,L)
          VTENDGW(I,L)=VTENDGW(I,L)+VNOTND(I,L) 
  525 CONTINUE
C
      IF(IFIZGWD.NE.0)                                      THEN
C
C       * STORE NON-OROGRAPHIC GWD ADJUSTED WIND TO LATER DETERMINE
C       * ITS TENDENCY.
C

        DO 530 L=1,ILEV
        DO 530 I=IL1,IL2
          USROW(I,L)=UM(I,L)
          VSROW(I,L)=VM(I,L)
  530   CONTINUE
      ENDIF


C  
C     * APPLY ROOF DRAG.
C  
      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 (UM(I,L)**2 + VM(I,L)**2) 
          VELRAT       = MAX(VELMOD/QVCRIT-1.,ZERO)
          DFAC         = VELRAT/(ZTMST)
          DENOM        = UN + ZTMST*DFAC
          UM(I,L)    = UM(I,L)/DENOM
          VM(I,L)    = VM(I,L)/DENOM
          UTENDGW(I,L) = UTENDGW(I,L) - DFAC*UM(I,L)
          VTENDGW(I,L) = VTENDGW(I,L) - DFAC*VM(I,L)
  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) + UM(I,L)
              VSUM(L) = VSUM(L) + VM(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
            UM(I,L)    = ( UM(I,L)+ZTMST*DFAC*UMEAN(I,L) )
     1                      / DENOM
            VM(I,L)    = ( VM(I,L)+ZTMST*DFAC*VMEAN(I,L) )
     1                      / DENOM

C           SPONGE AND EVD DRAG APPLIED ON NON-ZONAL WIND ONLY.

            UTSPNG(I,L)  =  - DFAC * ( UM(I,L)-UMEAN(I,L) )  
            VTSPNG(I,L)  =  - DFAC * ( VM(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
      IF(IFIZGWD.NE.0)                                      THEN
C
C       * STORE NON-OROGRAPHIC GWD ADJUSTED WIND TO LATER DETERMINE
C       * ITS TENDENCY.
C
        DO L=1,ILEV
        DO I=IL1,IL2
          UFROW(I,L)=UM(I,L)
          VFROW(I,L)=VM(I,L)
        ENDDO
        ENDDO
      ENDIF

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