gwdps.f

 
 
      module gwdps
 
      contains
 
      subroutine gwdps_run(                                             &
     &           IM,KM,A,B,C,ten_q,U1,V1,T1,Q1,KPBL,                    &
     &           PRSI,DEL,PRSL,PRSLK,PHII, PHIL,DELTIM,KDT,             &
     &           HPRIME,OC,OA4,CLX4,THETA,SIGMA,GAMMA,ELVMAX,           &
     &           DUSFC,DVSFC,dtaux2d_ms,dtauy2d_ms,dtaux2d_bl,          &
     &           dtauy2d_bl,dusfc_ms,dvsfc_ms,dusfc_bl,dvsfc_bl,        &
     &           G, CP, RD, RV, IMX,                                    &
     &           nmtvr, cdmbgwd, me, lprnt, ipr, rdxzb, ldiag_ugwp,     &
     &           errmsg, errflg)
 
!
!   ********************************************************************
! ----->  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(IM,KM)  NON-LIN TENDENCY FOR V WIND COMPONENT
!        B(IM,KM)  NON-LIN TENDENCY FOR U WIND COMPONENT
!        C(IM,KM)  NON-LIN TENDENCY FOR TEMPERATURE
!        U1(IM,KM) ZONAL WIND M/SEC  AT T0-DT
!        V1(IM,KM) MERIDIONAL WIND M/SEC AT T0-DT
!        T1(IM,KM) TEMPERATURE DEG K AT T0-DT
!        Q1(IM,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
!
      ! Interface variables
      integer, intent(in) :: im, km, imx, kdt, ipr, me
      integer, intent(in) :: KPBL(:) ! Index for the PBL top layer!
      real(kind=kind_phys), intent(in) ::                               &
     &                     deltim, g, cp, rd, rv, cdmbgwd(:)
      real(kind=kind_phys), intent(out) ::                              &
     &                     a(:,:), b(:,:), c(:,:)
      real(kind=kind_phys), intent(out) :: ten_q(:,:,:)
      real(kind=kind_phys), intent(in) ::                               &
     &                     u1(:,:),   v1(:,:),     t1(:,:),             &
     &                     q1(:,:),   prsi(:,:),   del(:,:),            &
     &                     prsl(:,:), prslk(:,:),  phil(:,:),           &
     &                     phii(:,:)
      real(kind=kind_phys), intent(in) ::                               &
     &                     oc(:), oa4(:,:), clx4(:,:), hprime(:)
      real(kind=kind_phys), intent(inout) :: elvmax(:)
      real(kind=kind_phys), intent(in) ::                               &
     &                     theta(:), sigma(:), gamma(:)
      real(kind=kind_phys), intent(inout) :: dusfc(:), dvsfc(:),        &
     &                     rdxzb(:)
      real(kind=kind_phys), intent(inout), optional :: dtaux2d_ms(:,:), &
     &                     dtauy2d_ms(:,:), dtaux2d_bl(:,:),            &
     &                     dtauy2d_bl(:,:)
      real(kind=kind_phys), intent(inout), optional :: dusfc_ms(:),     &
     &     dvsfc_ms(:), dusfc_bl(:), dvsfc_bl(:)
      integer, intent(in) :: nmtvr
      logical, intent(in) :: lprnt
      logical, intent(in) :: ldiag_ugwp
      character(len=*), intent(out) :: errmsg
      integer,          intent(out) :: errflg
!
      ! Local variables
! for lm mtn blocking
      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, rtrm, phiang, cdmb, dbim, zr, cdmbo4
      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.)   
!  hncrit set to 8000m and sigfac added to enhance elvmax mtn hgt
 
      parameter(sigfac=4.0)     
      parameter(hminmt=50.)     
      parameter(minwnd=0.1)     
 
!     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
                                 
      parameter(dpmin=5000.0)   
!
      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)                          &
     &,                    dtfac(im), xlinv(im),  delks(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),  cdsigohp(im)
!
!     real(kind=kind_phys) VELKO(KM-1)
      integer   kref(IM), kint(im), iwk(im), ipt(im)
! for lm mtn blocking
      integer   iwklm(im)
!     integer   kreflm(IM), iwklm(im)
      integer   idxzb(im), ktrial, klevm1
!
      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, onebg          &
     &,                    scork,  rscor, hd,    fro,   rim,  sira      &
     &,                    dtaux,  dtauy, pkp1log, pklog                &
     &,                    cosang, sinang, cos2a, sin2a, oneocpdt
!
      integer kmm1, kmm2, lcap, lcapp1, kbps, kbpsp1,kbpsm1             &
     &, kmps, idir, nwd, i, j, k, klcap, kp1, kmpbl, npt, npr, kmll
!    &, kmll,kmds,ihit,jhit
!
      ! Initialize CCPP error handling variables
      errmsg = ''
      errflg = 0
      
      a = 0.0
      b = 0.0
      c = 0.0
      ten_q = 0.0
!
!     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)
      cdmbo4 = 0.25 * cdmb
!
      npr = 0
!
      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
      onebg    = 1.0 / g
      gor      = g/rd
      gr2      = g*gor
      gocp     = g/cp
      fv       = rv/rd - 1
      oneocpdt = 1.0 / (cp*deltim)
!
!     NCNT   = 0
      kmm1   = km - 1
      kmm2   = km - 2
      lcap   = km
      lcapp1 = lcap + 1
!
      rdxzb(:)  = 0 
!
      IF ( nmtvr == 14) then 
! ----  for lm and gwd calculation points
        ipt = 0
        npt = 0
        DO i = 1,im
          IF (elvmax(i) > hminmt .and. hprime(i) > hpmin)  then
            npt      = npt + 1
            ipt(npt) = i
!           if (lprnt .and. ipr == i) npr = npt
          ENDIF
        ENDDO
        IF (npt == 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,IY,km,me=',npt,IM,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)
            cdsigohp(i) = cdmbo4 * sigma(j) / hprime(j)
          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) * onebg
            pklog   =  phil(j,k)   * onebg
!!!-------     ELVMAX(J) = min (ELVMAX(J) + sigfac * hprime(j), hncrit)
 
            if (elvmax(j) <=  pkp1log .and. elvmax(j) >=  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.0+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
          kp1 = k + 1
          DO i = 1, npt
           j   = ipt(i)
           rdz = g / (phil(j,kp1) - phil(j,k))
! ---                               Brunt-Vaisala Frequency
            bnv2lm(i,k) = (g+g) * rdz * (vtk(i,kp1) - vtk(i,k))
     &                                / (vtk(i,kp1) + 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 / (PRSI(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) = (prsi(j,1)-prsl(j,1)) * delks(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 < iwklm(I) .and. kreflm(I) == 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, iwklm(i)-1
            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
            if (k < iwklm(i)-1) then
              rdelks   = (prsl(j,k)-prsl(j,k+1)) * delks(i)
            else
              rdelks   = (prsl(j,k)-prsi(j,k+1)) * delks(i)
            endif
            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) >  90. ) ang(i,k) = ang(i,k) - 180.
            if ( ang(i,k) < -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) == 0 ) then
              pe(i) = pe(i) + bnv2lm(i,k) * (g*elvmax(j) - phil(j,k))
     &                                    * (phii(j,k+1) - phii(j,k))
     &                                    * (onebg*onebg)
! --- 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) >= ek(i)) THEN
                 idxzb(i) = k
                 rdxzb(j) = real(k,kind=kind_phys)
              ENDIF
! --- 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) > 0) then
            DO k = idxzb(i), 1, -1
              IF (phil(j,idxzb(i)) >  phil(j,k)) then
 
                zlen = sqrt( (phil(j,idxzb(i)) - phil(j,k))
     &                     / (phil(j,k )       + g*hprime(j)) )
! --- lm eq 14:
 
                cosang = cos(ang(i,k))
                sinang = sin(ang(i,k))
                cos2a  = cosang * cosang
                sin2a  = sinang * sinang
                tem    = cos2a + gamma(j)*sin2a
                                                       ! Here Rtrm is 1.0/R
                                                       ! --------------------
                if (abs(tem) > 1.e-06) then
                  rtrm   = (gamma(j)*cos2a + sin2a) / tem
                elseif (tem > 0.0) then
                  rtrm   = (gamma(j)*cos2a + sin2a) * 1.0e6
                else
                  rtrm   = - (gamma(j)*cos2a + sin2a) * 1.0e6
                endif
                zr =  max( 2.0 - rtrm, 0. )
 
! --- (negitive of DB -- see sign at tendency)
 
                db(i,k) = cdsigohp(i) * zr * ro(i,k) * zlen
     &                  * max(cosang, gamma(j)*sinang) * 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 /= 14) then
! ----  for mb not present and  gwd (nmtvr .ne .14)
        ipt = 0
        npt = 0
        DO i = 1,im
          IF ( hprime(i) > hpmin )  then
             npt      = npt + 1
             ipt(npt) = i
             if (ipr == i) npr = npt
          ENDIF
        ENDDO
        IF (npt == 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 > 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.0+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
        kp1 = k + 1
        DO i =1,npt
          j         = ipt(i)
          ti        = 2.0 / (t1(j,k)+t1(j,kp1))
          tem       = ti  / (prsl(j,k)-prsl(j,kp1))
          rdz       = g   / (phil(j,kp1) - phil(j,k))
          tem1      = u1(j,k) - u1(j,kp1)
          tem2      = v1(j,k) - v1(j,kp1)
          dw2       = tem1*tem1 + tem2*tem2
          shr2      = max(dw2,dw2min) * rdz * rdz
          bvf2      = g*(gocp+rdz*(vtj(i,kp1)-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,kp1)-vtk(i,k))
     &                            / (vtk(i,kp1)+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 < 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 / (PRSI(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) = (prsi(j,1)-prsl(j,1)) * delks(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 < 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
            if (k < kref(i)-1) then
              rdelks   = (prsl(j,k)-prsl(j,k+1)) * delks(i)
            else
              rdelks   = (prsl(j,k)-prsi(j,k+1)) * delks(i)
            endif
            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
        kp1 = k + 1
        DO  i = 1,npt
          j            = ipt(i)
          velco(i,k)   = 0.5 * ((u1(j,k)+u1(j,kp1))*ubar(i)
     &                       +  (v1(j,k)+v1(j,kp1))*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 <= 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 >= kref(i)) THEN
            icrilv(i) = icrilv(i) .OR. ( ri_n(i,k) <  ric)
     &                            .OR. (velco(i,k) <= 0.0)
          ENDIF
        ENDDO
!
        DO i = 1,npt
          IF (k >= kref(i))   THEN
            IF (.NOT.icrilv(i) .AND. taup(i,k) > 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) > 0. .AND. kp1 < 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 <= ric .AND.
!    &           (OA(I) .LE. 0. .OR.  PRSI(ipt(I),KP1).LE.RLOLEV )) THEN
     &           (oa(i) <= 0. .OR.  kp1 >= 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 <= 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 > kref(i) .and. prsi(ipt(i),k) >= rlolev) THEN
             IF(taud(i,k) /= 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 > 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))
 
          if (k < idxzb(i)) then                       ! ---  lm mb (*j*)  changes overwrite GWD
                                                       ! ---------------------------------------
            dbim     =   db(i,k) / (1.+db(i,k)*deltim)
            a(j,k)   = - dbim * v1(j,k)
            b(j,k)   = - dbim * u1(j,k)
            eng1     = eng0*(1.0-dbim*deltim)*(1.0-dbim*deltim)
 
!          if ( ABS(DBIM * U1(J,K)) > .01 )
!    & print *,' in gwdps_lmi.f KDT=',KDT,I,K,DB(I,K),
!    &                      dbim,idxzb(I),U1(J,K),V1(J,K),me
 
            tem1     = dbim * del(j,k)
            dusfc(j) = dusfc(j) + onebg * tem1 * u1(j,k)
            dvsfc(j) = dvsfc(j) + onebg * tem1 * v1(j,k)
            ! Output blocking tendencies if ldiag_ugwp=.true.
            if (ldiag_ugwp) then
               dtaux2d_bl(j,k) = - dbim * u1(j,k)
               dtauy2d_bl(j,k) = - dbim * v1(j,k)
               dusfc_bl(j) = dusfc_bl(j) + dtaux2d_bl(j,k) * del(j,k)
               dvsfc_bl(j) = dvsfc_bl(j) + dtauy2d_bl(j,k) * del(j,k)
            end if
          else                                         ! orographic GWD applied
                                                       ! ----------------------
            a(j,k)   = dtauy
            b(j,k)   = dtaux
            tem1     = u1(j,k) + dtaux*deltim
            tem2     = v1(j,k) + dtauy*deltim
            eng1     = 0.5 * (tem1*tem1+tem2*tem2)
            dusfc(j) = dusfc(j) - onebg * dtaux * del(j,k)
            dvsfc(j) = dvsfc(j) - onebg * dtauy * del(j,k)
            ! Output mesoscale GWD tendencies if ldiag_ugwp=.ture.
            if (ldiag_ugwp) then
               dtaux2d_ms(j,k) = dtaux
               dtauy2d_ms(j,k) = dtauy
               dusfc_ms(j) = dusfc_ms(j) + dtaux * del(j,k)
               dvsfc_ms(j) = dvsfc_ms(j) + dtauy * del(j,k)
            end if
          endif
          c(j,k) = max(eng0-eng1,0.) * oneocpdt
        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
 
      if (ldiag_ugwp) then
         !  Finalize dusfc and dvsfc diagnostics
         DO i = 1,npt
            j = ipt(i)
            dusfc_ms(j) = - onebg * dusfc_ms(j)
            dvsfc_ms(j) = - onebg * dvsfc_ms(j)
            dusfc_bl(j) = - onebg * dusfc_bl(j)
            dvsfc_bl(j) = - onebg * dvsfc_bl(j)
         ENDDO
      end if
!
!    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 subroutine gwdps_run
      end module gwdps