gwdps.f

Go to the documentation of this file.

      SUBROUTINE gwdps(IM,IX,IY,KM,A,B,C,U1,V1,T1,Q1,KPBL,
     &               prsi,del,prsl,prslk,phii, phil,deltim,kdt,
     &               hprime,oc,oa4,clx4,theta,sigma,gamma,elvmax, 
     &               dusfc,dvsfc,g, cp, rd, rv, imx, 
     &               nmtvr, cdmbgwd, me, lprnt, ipr)
!
!   ********************************************************************
! ----->  I M P L E M E N T A T I O N    V E R S I O N   <----------
!
!          --- Not in this code --  History of GWDP at NCEP----
!              ----------------     -----------------------
!  VERSION 3  MODIFIED FOR GRAVITY WAVES, LOCATION: .FR30(V3GWD)  *J*
!---       3.1 INCLUDES VARIABLE SATURATION FLUX PROFILE CF ISIGST
!---       3.G INCLUDES PS COMBINED W/ PH (GLAS AND GFDL)
!-----         ALSO INCLUDED IS RI  SMOOTH OVER A THICK LOWER LAYER
!-----         ALSO INCLUDED IS DECREASE IN DE-ACC AT TOP BY 1/2
!-----     THE NMC GWD INCORPORATING BOTH GLAS(P&S) AND GFDL(MIGWD)
!-----        MOUNTAIN INDUCED GRAVITY WAVE DRAG 
!-----    CODE FROM .FR30(V3MONNX) FOR MONIN3
!-----        THIS VERSION (06 MAR 1987)
!-----        THIS VERSION (26 APR 1987)    3.G
!-----        THIS VERSION (01 MAY 1987)    3.9
!-----    CHANGE TO FORTRAN 77 (FEB 1989)     --- HANN-MING HENRY JUANG
!-----    20070601 ELVMAX bug fix (*j*)
!
!   VERSION 4
!                ----- This code -----
!
!-----   MODIFIED TO IMPLEMENT THE ENHANCED LOW TROPOSPHERIC GRAVITY
!-----   WAVE DRAG DEVELOPED BY KIM AND ARAKAWA(JAS, 1995).
!        Orographic Std Dev (hprime), Convexity (OC), Asymmetry (OA4)
!        and Lx (CLX4) are input topographic statistics needed.
!
!-----   PROGRAMMED AND DEBUGGED BY HONG, ALPERT AND KIM --- JAN 1996.
!-----   debugged again - moorthi and iredell --- may 1998.
!-----
!       Further Cleanup, optimization and modification
!                                       - S. Moorthi May 98, March 99.
!-----   modified for usgs orography data (ncep office note 424)
!        and with several bugs fixed  - moorthi and hong --- july 1999.
!
!-----   Modified & implemented into NRL NOGAPS
!                                       - Young-Joon Kim, July 2000
!-----
!   VERSION lm MB  (6): oz fix 8/2003
!                ----- This code -----
!
!------   Changed to include the Lott and Miller Mtn Blocking
!         with some modifications by (*j*)  4/02
!        From a Principal Coordinate calculation using the
!        Hi Res 8 minute orography, the Angle of the
!        mtn with that to the East (x) axis is THETA, the slope
!        parameter SIGMA. The anisotropy is in GAMMA - all  are input
!        topographic statistics needed.  These are calculated off-line
!        as a function of model resolution in the fortran code ml01rg2.f,
!        with script mlb2.sh.   (*j*)
!-----   gwdps_mb.f version (following lmi) elvmax < hncrit (*j*)
!        MB3a expt to enhance elvmax mtn hgt see sigfac & hncrit
!        gwdps_GWDFIX_v6.f FIXGWD GF6.0 20070608 sigfac=4.
!-----
!----------------------------------------------------------------------C
!    USE
!        ROUTINE IS CALLED FROM GBPHYS  (AFTER CALL TO MONNIN)
!
!    PURPOSE
!        USING THE GWD PARAMETERIZATIONS OF PS-GLAS AND PH-
!        GFDL TECHNIQUE.  THE TIME TENDENCIES OF U V
!        ARE ALTERED TO INCLUDE THE EFFECT OF MOUNTAIN INDUCED
!        GRAVITY WAVE DRAG FROM SUB-GRID SCALE OROGRAPHY INCLUDING
!        CONVECTIVE BREAKING, SHEAR BREAKING AND THE PRESENCE OF
!        CRITICAL LEVELS
!
!  INPUT
!        A(IY,KM)  NON-LIN TENDENCY FOR V WIND COMPONENT
!        B(IY,KM)  NON-LIN TENDENCY FOR U WIND COMPONENT
!        C(IY,KM)  NON-LIN TENDENCY FOR TEMPERATURE
!        U1(IX,KM) ZONAL WIND M/SEC  AT T0-DT
!        V1(IX,KM) MERIDIONAL WIND M/SEC AT T0-DT
!        T1(IX,KM) TEMPERATURE DEG K AT T0-DT
!        Q1(IX,KM) SPECIFIC HUMIDITY AT T0-DT
!
!        DELTIM  TIME STEP    SECS
!        SI(N)   P/PSFC AT BASE OF LAYER N
!        SL(N)   P/PSFC AT MIDDLE OF LAYER N
!        DEL(N)  POSITIVE INCREMENT OF P/PSFC ACROSS LAYER N
!        KPBL(IM) is the index of the top layer of the PBL
!        ipr & lprnt for diagnostics
!
!  OUTPUT
!        A, B    AS AUGMENTED BY TENDENCY DUE TO GWDPS
!                OTHER INPUT VARIABLES UNMODIFIED.
!  revision log:
!    May 2013  J. Wang change cleff back to opn setting
!    Jan 2014  J. Wang merge Henry and Fangin's dissipation heat in gfs to nems 
!     
!
!   ********************************************************************
      USE machine , ONLY : kind_phys
      implicit none
      integer im, iy, ix, km, imx, kdt, ipr, me
      integer KPBL(im)                 ! Index for the PBL top layer!
      real(kind=kind_phys) deltim, G, CP, RD, RV,      cdmbgwd(2)
      real(kind=kind_phys) A(iy,km),    B(iy,km),      C(iy,km),
     &                     u1(ix,km),   v1(ix,km),     t1(ix,km),
     &                     q1(ix,km),   prsi(ix,km+1), del(ix,km),
     &                     prsl(ix,km), prslk(ix,km),  phil(ix,km),
     &                     phii(ix,km+1)
      real(kind=kind_phys) OC(im),     OA4(iy,4), CLX4(iy,4)
     &,                    hprime(im)
! for lm mtn blocking
      real(kind=kind_phys) ELVMAX(im),THETA(im),SIGMA(im),GAMMA(im)
      real(kind=kind_phys) wk(im)
      real(kind=kind_phys) bnv2lm(im,km),PE(im),EK(im),ZBK(im),UP(im)
      real(kind=kind_phys) DB(im,km),ANG(im,km),UDS(im,km)
      real(kind=kind_phys) ZLEN, DBTMP, R, PHIANG, CDmb, DBIM
      real(kind=kind_phys) ENG0, ENG1
!
!     Some constants
!
      real(kind=kind_phys) pi, dw2min, rimin, ric, bnv2min, efmin
     &,                    efmax,hpmax,hpmin, rad_to_deg, deg_to_rad
      parameter(pi=3.1415926535897931)
      parameter(rad_to_deg=180.0/pi, deg_to_rad=pi/180.0)
      parameter(dw2min=1., rimin=-100., ric=0.25, bnv2min=1.0e-5)
!     PARAMETER (EFMIN=0.0, EFMAX=10.0, hpmax=200.0)
      parameter(efmin=0.0, efmax=10.0, hpmax=2400.0, hpmin=1.0)
!
      real(kind=kind_phys) FRC,    CE,     CEOFRC, frmax, CG, GMAX
     &,                    veleps, factop, rlolev, rdi
!     &,                    CRITAC, VELEPS, FACTOP, RLOLEV, RDI
      parameter(frc=1.0, ce=0.8, ceofrc=ce/frc, frmax=100., cg=0.5)
      parameter(gmax=1.0, veleps=1.0, factop=0.5)
!      parameter (GMAX=1.0, CRITAC=5.0E-4, VELEPS=1.0, FACTOP=0.5)
      parameter(rlolev=50000.0) 
!     parameter (RLOLEV=500.0) 
!     parameter (RLOLEV=0.5)
!
       real(kind=kind_phys) dpmin,hminmt,hncrit,minwnd,sigfac
! --- for lm mtn blocking
!     parameter (cdmb = 1.0)     ! non-dim sub grid mtn drag Amp (*j*)
      parameter(hncrit=8000.)   ! Max value in meters for ELVMAX (*j*)
!  hncrit set to 8000m and sigfac added to enhance elvmax mtn hgt
      parameter(sigfac=4.0)     ! MB3a expt test for ELVMAX factor (*j*)
      parameter(hminmt=50.)     ! min mtn height (*j*)
      parameter(minwnd=0.1)     ! min wind component (*j*)

!     parameter (dpmin=00.0)     ! Minimum thickness of the reference layer
!!    parameter (dpmin=05.0)     ! Minimum thickness of the reference layer
!     parameter (dpmin=20.0)     ! Minimum thickness of the reference layer
                                 ! in centibars
      parameter(dpmin=5000.0)   ! Minimum thickness of the reference layer
                                 ! in Pa
!
      real(kind=kind_phys) FDIR
      integer mdir
      parameter(mdir=8, fdir=mdir/(pi+pi))
      integer nwdir(mdir)
      data nwdir/6,7,5,8,2,3,1,4/
      save nwdir
!
      LOGICAL ICRILV(im)
!
!----   MOUNTAIN INDUCED GRAVITY WAVE DRAG
!
      real(kind=kind_phys) TAUB(im),  XN(im),     YN(im),    UBAR(im)
     &,                    vbar(im),  ulow(im),   oa(im),    clx(im)
     &,                    roll(im),  uloi(im),   dusfc(im), dvsfc(im)
     &,                    dtfac(im), xlinv(im),  delks(im), delks1(im)
!
      real(kind=kind_phys) BNV2(im,km),  TAUP(im,km+1), ri_n(im,km) 
     &,                    taud(im,km),  ro(im,km),     vtk(im,km)
     &,                    vtj(im,km),   scor(im),      velco(im,km-1)
     &,                    bnv2bar(im)
!
!     real(kind=kind_phys) VELKO(KM-1)
      Integer   kref(im), kint(im), iwk(im), ipt(im)
! for lm mtn blocking
      Integer   kreflm(im), iwklm(im)
      Integer   idxzb(im), ktrial, klevm1, nmtvr
!
      real(kind=kind_phys) gor,    gocp,  fv,    gr2,  bnv,  fr
     &,                    brvf,   cleff, tem,   tem1,  tem2, temc, temv
     &,                    wdir,   ti,    rdz,   dw2,   shr2, bvf2
     &,                    rdelks, efact, coefm, gfobnv
     &,                    scork,  rscor, hd,    fro,   rim,  sira
     &,                    dtaux,  dtauy, pkp1log, pklog
      integer kmm1, kmm2, lcap, lcapp1, kbps, kbpsp1,kbpsm1
     &, kmps, idir, nwd, i, j, k, klcap, kp1, kmpbl, npt, npr
     &, kmll
!    &, kmll,kmds,ihit,jhit
      logical lprnt
!
!     parameter (cdmb = 1.0)     ! non-dim sub grid mtn drag Amp (*j*)
! non-dim sub grid mtn drag Amp (*j*)
!     cdmb = 1.0/float(IMX/192)
!     cdmb = 192.0/float(IMX)
      cdmb = 4.0 * 192.0/float(imx)
      if (cdmbgwd(1) >= 0.0) cdmb = cdmb * cdmbgwd(1)
!
      npr = 0
      DO i = 1, im
         dusfc(i) = 0.
         dvsfc(i) = 0.
      ENDDO
!
      DO k = 1, km
        DO i = 1, im
          db(i,k)  = 0.
          ang(i,k) = 0.
          uds(i,k) = 0.
        ENDDO
      ENDDO
!
      rdi  = 1.0 / rd
      gor  = g/rd
      gr2  = g*gor
      gocp = g/cp
      fv   = rv/rd - 1
!
!     NCNT   = 0
      kmm1   = km - 1
      kmm2   = km - 2
      lcap   = km
      lcapp1 = lcap + 1
!
!
      IF ( nmtvr .eq. 14) then 
! ----  for lm and gwd calculation points
        ipt = 0
        npt = 0
        DO i = 1,im
          IF ( (elvmax(i) .GT. hminmt) 
     &       .and. (hprime(i) .GT. hpmin) )  then
             npt      = npt + 1
             ipt(npt) = i
             if (ipr .eq. i) npr = npt
          ENDIF
        ENDDO
        IF (npt .eq. 0) RETURN     ! No gwd/mb calculation done!
!
!       if (lprnt) print *,' npt=',npt,' npr=',npr,' ipr=',ipr,' im=',im
!    &,' ipt(npt)=',ipt(npt)
!
! --- iwklm is the level above the height of the of the mountain.
! --- idxzb is the level of the dividing streamline.
! INITIALIZE DIVIDING STREAMLINE (DS) CONTROL VECTOR
!
        do i=1,npt
          iwklm(i)  = 2
          idxzb(i)  = 0 
          kreflm(i) = 0
        enddo
!       if (lprnt) 
!    &  print *,' in gwdps_lm.f npt,IM,IX,IY,km,me=',npt,IM,IX,IY,km,me
!
!
!
!..............................
!..............................
!
!  (*j*)  11/03:  test upper limit on KMLL=km - 1
!      then do not need hncrit -- test with large hncrit first.
!       KMLL  = km / 2 ! maximum mtnlm height : # of vertical levels / 2
        kmll = kmm1
! --- No mtn should be as high as KMLL (so we do not have to start at 
! --- the top of the model but could do calc for all levels).
!
          DO i = 1, npt
            j = ipt(i)
            elvmax(j) = min(elvmax(j) + sigfac * hprime(j), hncrit)
          ENDDO
!
        DO k = 1,kmll
          DO i = 1, npt
            j = ipt(i)
! --- interpolate to max mtn height for index, iwklm(I) wk[gz]
! --- ELVMAX is limited to hncrit because to hi res topo30 orog.
            pkp1log =  phil(j,k+1) / g
            pklog =  phil(j,k)   / g
!!!-------     ELVMAX(J) = min (ELVMAX(J) + sigfac * hprime(j), hncrit)
            if ( ( elvmax(j) .le.  pkp1log ) .and. 
     &           ( elvmax(j) .ge.   pklog  ) ) THEN
!     print *,' in gwdps_lm.f 1  =',k,ELVMAX(j),pklog,pkp1log,me
! ---        wk for diags but can be saved and reused.  
               wk(i)  = g * elvmax(j) / ( phil(j,k+1) - phil(j,k) )
               iwklm(i)  =  max(iwklm(i), k+1 ) 
!     print *,' in gwdps_lm.f 2 npt=',npt,i,j,wk(i),iwklm(i),me
            endif
!
! ---        find at prsl levels large scale environment variables
! ---        these cover all possible mtn max heights
            vtj(i,k)  = t1(j,k)  * (1.+fv*q1(j,k))
            vtk(i,k)  = vtj(i,k) / prslk(j,k)
            ro(i,k)   = rdi * prsl(j,k) / vtj(i,k) ! DENSITY Kg/M**3
          ENDDO
        ENDDO
!
! testing for highest model level of mountain top
!
!         ihit = 2
!         jhit = 0
!        do i = 1, npt
!        j=ipt(i)
!          if ( iwklm(i) .gt. ihit ) then 
!            ihit = iwklm(i)
!            jhit = j
!          endif
!        enddo
!     print *, ' mb: kdt,max(iwklm),jhit,phil,me=',
!    &          kdt,ihit,jhit,phil(jhit,ihit),me
         
        klevm1 = kmll - 1
        DO k = 1, klevm1  
          DO i = 1, npt
           j   = ipt(i)
            rdz  = g   / ( phil(j,k+1) - phil(j,k) )
! ---                               Brunt-Vaisala Frequency
            bnv2lm(i,k) = (g+g) * rdz * ( vtk(i,k+1)-vtk(i,k) )
     &                     / ( vtk(i,k+1)+vtk(i,k) )
            bnv2lm(i,k) = max( bnv2lm(i,k), bnv2min )
          ENDDO
        ENDDO
!    print *,' in gwdps_lm.f 3 npt=',npt,j,RDZ,me
!
        DO i = 1, npt
          j   = ipt(i)
          delks(i)  = 1.0 / (prsi(j,1) - prsi(j,iwklm(i)))
          delks1(i) = 1.0 / (prsl(j,1) - prsl(j,iwklm(i)))
          ubar(i)  = 0.0
          vbar(i)  = 0.0
          roll(i)  = 0.0
          pe  (i)  = 0.0
          ek  (i)  = 0.0
          bnv2bar(i) = (prsl(j,1)-prsl(j,2)) * delks1(i) * bnv2lm(i,1)
        ENDDO

! --- find the dividing stream line height 
! --- starting from the level above the max mtn downward
! --- iwklm(i) is the k-index of mtn elvmax elevation
        DO ktrial = kmll, 1, -1
          DO i = 1, npt
             IF ( ktrial .LT. iwklm(i) .and. kreflm(i) .eq. 0 ) then
                kreflm(i) = ktrial
             ENDIF
          ENDDO
        ENDDO
!     print *,' in gwdps_lm.f 4 npt=',npt,kreflm(npt),me
!
! --- in the layer kreflm(I) to 1 find PE (which needs N, ELVMAX)
! ---  make averages, guess dividing stream (DS) line layer.
! ---  This is not used in the first cut except for testing and
! --- is the vert ave of quantities from the surface to mtn top.
!   
        DO i = 1, npt
          DO k = 1, kreflm(i)
            j        = ipt(i)
            rdelks     = del(j,k) * delks(i)
            ubar(i)    = ubar(i)  + rdelks * u1(j,k) ! trial Mean U below 
            vbar(i)    = vbar(i)  + rdelks * v1(j,k) ! trial Mean V below 
            roll(i)    = roll(i)  + rdelks * ro(i,k) ! trial Mean RO below 
            rdelks     = (prsl(j,k)-prsl(j,k+1)) * delks1(i)
            bnv2bar(i) = bnv2bar(i) + bnv2lm(i,k) * rdelks
! --- these vert ave are for diags, testing and GWD to follow (*j*).
          ENDDO
        ENDDO
!     print *,' in gwdps_lm.f 5  =',i,kreflm(npt),BNV2bar(npt),me
!
! --- integrate to get PE in the trial layer.
! --- Need the first layer where PE>EK - as soon as 
! --- IDXZB is not 0 we have a hit and Zb is found.
!
        DO i = 1, npt
          j = ipt(i)
          DO k = iwklm(i), 1, -1
            phiang   =  atan2(v1(j,k),u1(j,k))*rad_to_deg
            ang(i,k) = ( theta(j) - phiang )
            if ( ang(i,k) .gt.  90. ) ang(i,k) = ang(i,k) - 180.
            if ( ang(i,k) .lt. -90. ) ang(i,k) = ang(i,k) + 180.
            ang(i,k) = ang(i,k) * deg_to_rad
!
            uds(i,k) = 
     &          max(sqrt(u1(j,k)*u1(j,k) + v1(j,k)*v1(j,k)), minwnd)
! --- Test to see if we found Zb previously
            IF (idxzb(i) .eq. 0 ) then
              pe(i) = pe(i) + bnv2lm(i,k) * 
     &           ( g * elvmax(j) - phil(j,k) ) * 
     &           ( phii(j,k+1) - phii(j,k) ) / (g*g)
! --- KE
! --- Wind projected on the line perpendicular to mtn range, U(Zb(K)).
! --- kenetic energy is at the layer Zb
! --- THETA ranges from -+90deg |_ to the mtn "largest topo variations"
              up(i)  =  uds(i,k) * cos(ang(i,k))
              ek(i)  = 0.5 *  up(i) * up(i) 

! --- Dividing Stream lime  is found when PE =exceeds EK.
              IF ( pe(i) .ge.  ek(i) ) idxzb(i) = k
! --- Then mtn blocked flow is between Zb=k(IDXZB(I)) and surface
!
            ENDIF
          ENDDO
        ENDDO
!
!     print *,' in gwdps_lm.f 6  =',phiang,THETA(ipt(npt)),me
!     print *,' in gwdps_lm.f 7  =',IDXZB(npt),PE(npt)
!
!     if (lprnt .and. npr .gt. 0) then
!       print *,' BNV2bar,BNV2lm=',bnv2bar(npr),BNV2lm(npr,1:klevm1)
!       print *,' npr,IDXZB,UDS=',npr,IDXZB(npr),UDS(npr,:)
!       print *,' PE,UP,EK=',PE(npr),UP(npr),EK(npr)
!     endif
!
        DO i = 1, npt
          j    = ipt(i)
! --- Calc if N constant in layers (Zb guess) - a diagnostic only.
          zbk(i) = elvmax(j)
     &           - sqrt(ubar(i)*ubar(i) + vbar(i)*vbar(i))/bnv2bar(i)
        ENDDO
!
!     if (lprnt .and. npr .gt. 0) then
!       print *,' iwklm,ZBK=',iwklm(npr),ZBK(npr),IDXZB(npr)
!       print *,' Zb=',PHIL(ipr),IDXZB(npr))/G
!     print *,' in gwdps_lm.f 8 npt =',npt,ZBK(npt),UP(npt),me
!     endif
!
! --- The drag for mtn blocked flow
! 
        DO i = 1, npt
          j = ipt(i)
          zlen = 0.
!      print *,' in gwdps_lm.f 9  =',i,j,IDXZB(i),me
          IF ( idxzb(i) .gt. 0 ) then 
            DO k = idxzb(i), 1, -1
              IF ( phil(j,idxzb(i)) .gt.  phil(j,k) ) then

                zlen = sqrt( ( phil(j,idxzb(i)) - phil(j,k) ) / 
     &                       ( phil(j,k ) + g * hprime(j) ) )
! --- lm eq 14:
                r = (cos(ang(i,k))**2 + gamma(j) * sin(ang(i,k))**2) / 
     &              (gamma(j) * cos(ang(i,k))**2 + sin(ang(i,k))**2)
! --- (negitive of DB -- see sign at tendency)

                dbtmp = 0.25 *  cdmb *
     &                  max( 2. - 1. / r, 0. ) * sigma(j) *
     &                  max(cos(ang(i,k)), gamma(j)*sin(ang(i,k))) *
     &                  zlen / hprime(j) 
                db(i,k) =  dbtmp * uds(i,k)    
!
!               if(lprnt .and. i .eq. npr) then 
!                 print *,' in gwdps_lmi.f 10 npt=',npt,i,j,idxzb(i)
!    &,           DBTMP,R' ang=',ang(i,k),' gamma=',gamma(j),' K=',K
!                 print *,' in gwdps_lmi.f 11   K=',k,ZLEN,cos(ANG(I,K))
!                 print *,' in gwdps_lmi.f 12  DB=',DB(i,k),sin(ANG(I,K))
!               endif
              endif
            ENDDO
!         if(lprnt) print *,' @K=1,ZLEN,DBTMP=',K,ZLEN,DBTMP
          endif
        ENDDO
! 
!.............................
!.............................
! end  mtn blocking section
!
      ELSEIF ( nmtvr .ne. 14) then 
! ----  for mb not present and  gwd (nmtvr .ne .14) 
        ipt     = 0
        npt     = 0
        DO i = 1,im
          IF ( hprime(i) .GT. hpmin )  then
             npt      = npt + 1
             ipt(npt) = i
             if (ipr .eq. i) npr = npt
          ENDIF
        ENDDO
        IF (npt .eq. 0) RETURN     ! No gwd/mb calculation done!
!
!       if (lprnt) print *,' NPR=',npr,' npt=',npt,' IPR=',IPR
!      &,' ipt(npt)=',ipt(npt)
!
        do i=1,npt
          idxzb(i) = 0
        enddo
      ENDIF
!
!.............................
!.............................
!
      kmpbl  = km / 2 ! maximum pbl height : # of vertical levels / 2
!
!  Scale cleff between IM=384*2 and 192*2 for T126/T170 and T62
!
      if (imx .gt. 0) then
!       cleff = 1.0E-5 * SQRT(FLOAT(IMX)/384.0) !  this is inverse of CLEFF!
!       cleff = 1.0E-5 * SQRT(FLOAT(IMX)/192.0) !  this is inverse of CLEFF!
!       cleff = 0.5E-5 * SQRT(FLOAT(IMX)/192.0) !  this is inverse of CLEFF!
!       cleff = 1.0E-5 * SQRT(FLOAT(IMX)/192)/float(IMX/192)
!       cleff = 1.0E-5 / SQRT(FLOAT(IMX)/192.0) !  this is inverse of CLEFF!
        cleff = 0.5e-5 / sqrt(float(imx)/192.0) !  this is inverse of CLEFF!
! hmhj for ndsl
! jw        cleff = 0.1E-5 / SQRT(FLOAT(IMX)/192.0) !  this is inverse of CLEFF!
!       cleff = 2.0E-5 * SQRT(FLOAT(IMX)/192.0) !  this is inverse of CLEFF!
!       cleff = 2.5E-5 * SQRT(FLOAT(IMX)/192.0) !  this is inverse of CLEFF!
      endif
      if (cdmbgwd(2) >= 0.0) cleff = cleff * cdmbgwd(2)
!
      DO k = 1,km
        DO i =1,npt
          j         = ipt(i)
          vtj(i,k)  = t1(j,k)  * (1.+fv*q1(j,k))
          vtk(i,k)  = vtj(i,k) / prslk(j,k)
          ro(i,k)   = rdi * prsl(j,k) / vtj(i,k) ! DENSITY TONS/M**3
          taup(i,k) = 0.0
        ENDDO
      ENDDO
      DO k = 1,kmm1
        DO i =1,npt
          j         = ipt(i)
          ti        = 2.0 / (t1(j,k)+t1(j,k+1))
          tem       = ti  / (prsl(j,k)-prsl(j,k+1))
          rdz       = g   / (phil(j,k+1) - phil(j,k))
          tem1      = u1(j,k) - u1(j,k+1)
          tem2      = v1(j,k) - v1(j,k+1)
          dw2       = tem1*tem1 + tem2*tem2
          shr2      = max(dw2,dw2min) * rdz * rdz
          bvf2      = g*(gocp+rdz*(vtj(i,k+1)-vtj(i,k))) * ti
          ri_n(i,k) = max(bvf2/shr2,rimin)   ! Richardson number
!                                              Brunt-Vaisala Frequency
!         TEM       = GR2 * (PRSL(J,K)+PRSL(J,K+1)) * TEM
!         BNV2(I,K) = TEM * (VTK(I,K+1)-VTK(I,K))/(VTK(I,K+1)+VTK(I,K))
          bnv2(i,k) = (g+g) * rdz * (vtk(i,k+1)-vtk(i,k))
     &                            / (vtk(i,k+1)+vtk(i,k))
          bnv2(i,k) = max( bnv2(i,k), bnv2min )
        ENDDO
      ENDDO
!      print *,' in gwdps_lm.f GWD:14  =',npt,kmm1,bnv2(npt,kmm1)
!
!     Apply 3 point smoothing on BNV2
!
!     do k=1,km
!       do i=1,im
!         vtk(i,k) = bnv2(i,k)
!       enddo
!     enddo
!     do k=2,kmm1
!       do i=1,im
!         bnv2(i,k) = 0.25*(vtk(i,k-1)+vtk(i,k+1)) + 0.5*vtk(i,k)
!       enddo
!     enddo
!
!     Finding the first interface index above 50 hPa level
!
      do i=1,npt
        iwk(i) = 2
      enddo
      DO k=3,kmpbl
        DO i=1,npt
          j   = ipt(i)
          tem = (prsi(j,1) - prsi(j,k))
          if (tem .lt. dpmin) iwk(i) = k
        enddo
      enddo
!
      kbps = 1
      kmps = km
      DO i=1,npt
        j         = ipt(i)
        kref(i)   = max(iwk(i), kpbl(j)+1 ) ! reference level 
        delks(i)  = 1.0 / (prsi(j,1) - prsi(j,kref(i)))
        delks1(i) = 1.0 / (prsl(j,1) - prsl(j,kref(i)))
        ubar(i)  = 0.0
        vbar(i)  = 0.0
        roll(i)  = 0.0
        kbps      = max(kbps,  kref(i))
        kmps      = min(kmps,  kref(i))
!
        bnv2bar(i) = (prsl(j,1)-prsl(j,2)) * delks1(i) * bnv2(i,1)
      ENDDO
!      print *,' in gwdps_lm.f GWD:15  =',KBPS,KMPS
      kbpsp1 = kbps + 1
      kbpsm1 = kbps - 1
      DO k = 1,kbps
        DO i = 1,npt
          IF (k .LT. kref(i)) THEN
            j          = ipt(i)
            rdelks     = del(j,k) * delks(i)
            ubar(i)    = ubar(i)  + rdelks * u1(j,k)   ! Mean U below kref
            vbar(i)    = vbar(i)  + rdelks * v1(j,k)   ! Mean V below kref
!
            roll(i)    = roll(i)  + rdelks * ro(i,k)   ! Mean RO below kref
            rdelks     = (prsl(j,k)-prsl(j,k+1)) * delks1(i)
            bnv2bar(i) = bnv2bar(i) + bnv2(i,k) * rdelks
          ENDIF
        ENDDO
      ENDDO
!      print *,' in gwdps_lm.f GWD:15B =',bnv2bar(npt)
!
!     FIGURE OUT LOW-LEVEL HORIZONTAL WIND DIRECTION AND FIND 'OA'
!
!             NWD  1   2   3   4   5   6   7   8
!              WD  W   S  SW  NW   E   N  NE  SE
!
      DO i = 1,npt
        j      = ipt(i)
        wdir   = atan2(ubar(i),vbar(i)) + pi
        idir   = mod(nint(fdir*wdir),mdir) + 1
        nwd    = nwdir(idir)
        oa(i)  = (1-2*int( (nwd-1)/4 )) * oa4(j,mod(nwd-1,4)+1)
        clx(i) = clx4(j,mod(nwd-1,4)+1)
      ENDDO
!
!-----XN,YN            "LOW-LEVEL" WIND PROJECTIONS IN ZONAL
!                                    & MERIDIONAL DIRECTIONS
!-----ULOW             "LOW-LEVEL" WIND MAGNITUDE -        (= U)
!-----BNV2             BNV2 = N**2
!-----TAUB             BASE MOMENTUM FLUX
!-----= -(RO * U**3/(N*XL)*GF(FR) FOR N**2 > 0
!-----= 0.                        FOR N**2 < 0
!-----FR               FROUDE    =   N*HPRIME / U
!-----G                GMAX*FR**2/(FR**2+CG/OC)
!
!-----INITIALIZE SOME ARRAYS
!
      DO i = 1,npt
        xn(i)     = 0.0
        yn(i)     = 0.0
        taub(i)  = 0.0
        ulow(i)  = 0.0
        dtfac(i)  = 1.0
        icrilv(i) = .false. ! INITIALIZE CRITICAL LEVEL CONTROL VECTOR
        
!
!----COMPUTE THE "LOW LEVEL" WIND MAGNITUDE (M/S)
!
        ulow(i) = max(sqrt(ubar(i)*ubar(i) + vbar(i)*vbar(i)), 1.0)
        uloi(i) = 1.0 / ulow(i)
      ENDDO
!
      DO  k = 1,kmm1
        DO  i = 1,npt
          j            = ipt(i)
          velco(i,k)   = 0.5 * ((u1(j,k)+u1(j,k+1))*ubar(i)
     &                       +  (v1(j,k)+v1(j,k+1))*vbar(i))
          velco(i,k)   = velco(i,k) * uloi(i)
!         IF ((VELCO(I,K).LT.VELEPS) .AND. (VELCO(I,K).GT.0.)) THEN
!           VELCO(I,K) = VELEPS
!         ENDIF
        ENDDO
      ENDDO
!      
!
!   find the interface level of the projected wind where
!   low levels & upper levels meet above pbl
!
!     do i=1,npt
!       kint(i) = km
!     enddo
!     do k = 1,kmm1
!       do i = 1,npt
!         IF (K .GT. kref(I)) THEN
!           if(velco(i,k) .lt. veleps .and. kint(i) .eq. km) then
!             kint(i) = k+1
!           endif
!         endif
!       enddo
!     enddo
!  WARNING  KINT = KREF !!!!!!!!!
      do i=1,npt
        kint(i) = kref(i)
      enddo
!
!     if(lprnt) print *,' ubar=',ubar
!    &,' vbar=',vbar,' ulow=',ulow,' veleps=',veleps
!
      DO i = 1,npt
        j      = ipt(i)
        bnv    = sqrt( bnv2bar(i) )
        fr     = bnv     * uloi(i) * min(hprime(j),hpmax)
        fr     = min(fr, frmax)
        xn(i)  = ubar(i) * uloi(i)
        yn(i)  = vbar(i) * uloi(i)
!
!     Compute the base level stress and store it in TAUB
!     CALCULATE ENHANCEMENT FACTOR, NUMBER OF MOUNTAINS & ASPECT
!     RATIO CONST. USE SIMPLIFIED RELATIONSHIP BETWEEN STANDARD
!     DEVIATION & CRITICAL HGT
!


        efact    = (oa(i) + 2.) ** (ceofrc*fr)
        efact    = min( max(efact,efmin), efmax )
!
        coefm    = (1. + clx(i)) ** (oa(i)+1.)
!
        xlinv(i) = coefm * cleff
!
        tem      = fr    * fr * oc(j)
        gfobnv   = gmax  * tem / ((tem + cg)*bnv)  ! G/N0
!
        taub(i)  = xlinv(i) * roll(i) * ulow(i) * ulow(i)
     &           * ulow(i)  * gfobnv  * efact         ! BASE FLUX Tau0
!
!         tem      = min(HPRIME(I),hpmax)
!         TAUB(I)  = XLINV(I) * ROLL(I) * ULOW(I) * BNV * tem * tem
!
        k        = max(1, kref(i)-1)
        tem      = max(velco(i,k)*velco(i,k), 0.1)
        scor(i)  = bnv2(i,k) / tem  ! Scorer parameter below ref level
      ENDDO
!     if(lprnt) print *,' taub=',taub
!                                                                       
!----SET UP BOTTOM VALUES OF STRESS
!
      DO k = 1, kbps
        DO i = 1,npt
          IF (k .LE. kref(i)) taup(i,k) = taub(i)
        ENDDO
      ENDDO
!
!   Now compute vertical structure of the stress.
!
      DO k = kmps, kmm1                   ! Vertical Level K Loop!
        kp1 = k + 1
        DO i = 1, npt
!
!-----UNSTABLE LAYER IF RI < RIC
!-----UNSTABLE LAYER IF UPPER AIR VEL COMP ALONG SURF VEL <=0 (CRIT LAY)
!---- AT (U-C)=0. CRIT LAYER EXISTS AND BIT VECTOR SHOULD BE SET (.LE.)
!
          IF (k .GE. kref(i)) THEN
            icrilv(i) = icrilv(i) .OR. ( ri_n(i,k) .LT. ric)
     &                            .OR. (velco(i,k) .LE. 0.0)
          ENDIF
        ENDDO
!
        DO i = 1,npt
          IF (k .GE. kref(i))   THEN
            IF (.NOT.icrilv(i) .AND. taup(i,k) .GT. 0.0 ) THEN
              temv = 1.0 / max(velco(i,k), 0.01)
!             IF (OA(I) .GT. 0. .AND.  PRSI(ipt(i),KP1).GT.RLOLEV) THEN
              IF (oa(i).GT.0. .AND. kp1 .lt. kint(i)) THEN
                scork   = bnv2(i,k) * temv * temv
                rscor   = min(1.0, scork / scor(i))
                scor(i) = scork
              ELSE 
                rscor   = 1.
              ENDIF
!

              brvf = sqrt(bnv2(i,k))        ! Brunt-Vaisala Frequency
!             TEM1 = XLINV(I)*(RO(I,KP1)+RO(I,K))*BRVF*VELCO(I,K)*0.5
              tem1 = xlinv(i)*(ro(i,kp1)+ro(i,k))*brvf*0.5
     &                       * max(velco(i,k),0.01)
              hd   = sqrt(taup(i,k) / tem1)
              fro  = brvf * hd * temv
!
!    RIM is the  MINIMUM-RICHARDSON NUMBER BY SHUTTS (1985)
!

              tem2   = sqrt(ri_n(i,k))
              tem    = 1. + tem2 * fro
              rim    = ri_n(i,k) * (1.-fro) / (tem * tem)
!
!    CHECK STABILITY TO EMPLOY THE 'SATURATION HYPOTHESIS'
!    OF LINDZEN (1981) EXCEPT AT TROPOSPHERIC DOWNSTREAM REGIONS
!
!                                       ----------------------
              IF (rim .LE. ric .AND.
!    &           (OA(I) .LE. 0. .OR.  PRSI(ipt(I),KP1).LE.RLOLEV )) THEN
     &           (oa(i) .LE. 0. .OR.  kp1 .ge. kint(i) )) THEN
                 temc = 2.0 + 1.0 / tem2
                 hd   = velco(i,k) * (2.*sqrt(temc)-temc) / brvf
                 taup(i,kp1) = tem1 * hd * hd
              ELSE 
                 taup(i,kp1) = taup(i,k) * rscor
              ENDIF
              taup(i,kp1) = min(taup(i,kp1), taup(i,k))
            ENDIF
          ENDIF
        ENDDO
      ENDDO
!
!     DO I=1,IM
!       taup(i,km+1) = taup(i,km)
!     ENDDO
!
      IF(lcap .LE. km) THEN
         DO klcap = lcapp1, km+1
            DO i = 1,npt
              sira          = prsi(ipt(i),klcap) / prsi(ipt(i),lcap)
              taup(i,klcap) = sira * taup(i,lcap)
            ENDDO
         ENDDO
      ENDIF
!
!     Calculate - (g/p*)*d(tau)/d(sigma) and Decel terms DTAUX, DTAUY
!
      DO k = 1,km
        DO i = 1,npt
          taud(i,k) = g * (taup(i,k+1) - taup(i,k)) / del(ipt(i),k)
        ENDDO
      ENDDO
!
!------LIMIT DE-ACCELERATION (MOMENTUM DEPOSITION ) AT TOP TO 1/2 VALUE
!------THE IDEA IS SOME STUFF MUST GO OUT THE 'TOP'
!
      DO klcap = lcap, km
         DO i = 1,npt
            taud(i,klcap) = taud(i,klcap) * factop
         ENDDO
      ENDDO
!
!------IF THE GRAVITY WAVE DRAG WOULD FORCE A CRITICAL LINE IN THE
!------LAYERS BELOW SIGMA=RLOLEV DURING THE NEXT DELTIM TIMESTEP,
!------THEN ONLY APPLY DRAG UNTIL THAT CRITICAL LINE IS REACHED.
!
      DO k = 1,kmm1
        DO i = 1,npt
           IF (k .GT. kref(i) .and. prsi(ipt(i),k) .GE. rlolev) THEN
             IF(taud(i,k).NE.0.) THEN
               tem = deltim * taud(i,k)
               dtfac(i) = min(dtfac(i),abs(velco(i,k)/tem))
             ENDIF
           ENDIF
        ENDDO
      ENDDO
!
!     if(lprnt .and. npr .gt. 0) then
!       print *,' before  A=',A(npr,:)
!       print *,' before  B=',B(npr,:)
!     endif

      DO k = 1,km
        DO i = 1,npt
          j          = ipt(i)
          taud(i,k)  = taud(i,k) * dtfac(i)
          dtaux      = taud(i,k) * xn(i)
          dtauy      = taud(i,k) * yn(i)
          eng0       = 0.5*(u1(j,k)*u1(j,k)+v1(j,k)*v1(j,k))
! ---  lm mb (*j*)  changes overwrite GWD
          if ( k .lt. idxzb(i) .AND. idxzb(i) .ne. 0 ) then
            dbim = db(i,k) / (1.+db(i,k)*deltim)
            a(j,k)  = - dbim * v1(j,k) + a(j,k)
            b(j,k)  = - dbim * u1(j,k) + b(j,k)
            eng1    = eng0*(1.0-dbim*deltim)*(1.0-dbim*deltim)
!          if ( ABS(DBIM * U1(J,K)) .gt. .01 ) 
!    & print *,' in gwdps_lmi.f KDT=',KDT,I,K,DB(I,K),
!    &                      dbim,idxzb(I),U1(J,K),V1(J,K),me
            dusfc(j)   = dusfc(j) - dbim * v1(j,k) * del(j,k)
            dvsfc(j)   = dvsfc(j) - dbim * u1(j,k) * del(j,k)
          else
!
            a(j,k)     = dtauy     + a(j,k)
            b(j,k)     = dtaux     + b(j,k)
            eng1       = 0.5*(
     &                   (u1(j,k)+dtaux*deltim)*(u1(j,k)+dtaux*deltim)
     &                 + (v1(j,k)+dtauy*deltim)*(v1(j,k)+dtauy*deltim))
            dusfc(j)   = dusfc(j)  + dtaux * del(j,k)
            dvsfc(j)   = dvsfc(j)  + dtauy * del(j,k)
          endif
          c(j,k) = c(j,k) + max(eng0-eng1,0.)/cp/deltim
        ENDDO
      ENDDO
!     if (lprnt) then
!       print *,' in gwdps_lm.f after  A=',A(ipr,:)
!       print *,' in gwdps_lm.f after  B=',B(ipr,:)
!       print *,' DB=',DB(ipr,:)
!     endif
      tem    = -1.0/g
      DO i = 1,npt
        j          = ipt(i)
!       TEM    = (-1.E3/G)
        dusfc(j) = tem * dusfc(j)
        dvsfc(j) = tem * dvsfc(j)
      ENDDO
!                                                                       
!    MONITOR FOR EXCESSIVE GRAVITY WAVE DRAG TENDENCIES IF NCNT>0
!
!     IF(NCNT.GT.0) THEN
!        IF(LAT.GE.38.AND.LAT.LE.42) THEN
!CMIC$ GUARD 37
!           DO 92 I = 1,IM
!              IF(IKOUNT.GT.NCNT) GO TO 92
!              IF(I.LT.319.OR.I.GT.320) GO TO 92
!              DO 91 K = 1,KM
!                 IF(ABS(TAUD(I,K)) .GT. CRITAC) THEN
!                    IF(I.LE.IM) THEN
!                       IKOUNT = IKOUNT+1
!                       PRINT 123,I,LAT,KDT
!                       PRINT 124,TAUB(I),BNV(I),ULOW(I),
!    1                  GF(I),FR(I),ROLL(I),HPRIME(I),XN(I),YN(I)
!                       PRINT 124,(TAUD(I,KK),KK = 1,KM)
!                       PRINT 124,(TAUP(I,KK),KK = 1,KM+1)
!                       PRINT 124,(ri_n(I,KK),KK = 1,KM)
!                       DO 93 KK = 1,KMM1
!                          VELKO(KK) =
!    1                  0.5*((U1(I,KK)+U1(I,KK+1))*UBAR(I)+
!    2                  (V1(I,KK)+V1(I,KK+1))*VBAR(I))*ULOI(I)
!93                     CONTINUE
!                       PRINT 124,(VELKO(KK),KK = 1,KMM1)
!                       PRINT 124,(A    (I,KK),KK = 1,KM)
!                       PRINT 124,(DTAUY(I,KK),KK = 1,KM)
!                       PRINT 124,(B    (I,KK),KK = 1,KM)
!                       PRINT 124,(DTAUX(I,KK),KK = 1,KM)
!                       GO TO 92
!                    ENDIF
!                 ENDIF
!91            CONTINUE
!92         CONTINUE
!CMIC$ END GUARD 37
!123        FORMAT('  *** MIGWD PRINT *** I=',I3,' LAT=',I3,' KDT=',I3)
!124        FORMAT(2X,  10E13.6)
!        ENDIF
!     ENDIF
!
!      print *,' in gwdps_lm.f 18  =',A(ipt(1),idxzb(1))
!    &,                          B(ipt(1),idxzb(1)),me
      RETURN
      END