File: C:\NOAA\NEMS_11731\src\atmos\phys\gwdps.f

1     !     SUBROUTINE GWDPS(IM,IX,IY,KM,A,B,U1,V1,T1,Q1,PSTAR,KPBL,
2           SUBROUTINE GWDPS(IM,IX,IY,KM,A,B,U1,V1,T1,Q1,KPBL,
3          &               PRSI,DEL,PRSL,PRSLK,PHII, PHIL,DELTIM,KDT,
4          &               HPRIME,OC,OA4,CLX4,THETA,SIGMA,GAMMA,ELVMAX, 
5          &               DUSFC,DVSFC,G, CP, RD, RV, IMX, 
6          &               nmtvr, cdmbgwd, me, lprnt, ipr)
7     !
8     !   ********************************************************************
9     ! ----->  I M P L E M E N T A T I O N    V E R S I O N   <----------
10     !
11     !          --- Not in this code --  History of GWDP at NCEP----
12     !              ----------------     -----------------------
13     !  VERSION 3  MODIFIED FOR GRAVITY WAVES, LOCATION: .FR30(V3GWD)  *J*
14     !---       3.1 INCLUDES VARIABLE SATURATION FLUX PROFILE CF ISIGST
15     !---       3.G INCLUDES PS COMBINED W/ PH (GLAS AND GFDL)
16     !-----         ALSO INCLUDED IS RI  SMOOTH OVER A THICK LOWER LAYER
17     !-----         ALSO INCLUDED IS DECREASE IN DE-ACC AT TOP BY 1/2
18     !-----     THE NMC GWD INCORPORATING BOTH GLAS(P&S) AND GFDL(MIGWD)
19     !-----        MOUNTAIN INDUCED GRAVITY WAVE DRAG 
20     !-----    CODE FROM .FR30(V3MONNX) FOR MONIN3
21     !-----        THIS VERSION (06 MAR 1987)
22     !-----        THIS VERSION (26 APR 1987)    3.G
23     !-----        THIS VERSION (01 MAY 1987)    3.9
24     !-----    CHANGE TO FORTRAN 77 (FEB 1989)     --- HANN-MING HENRY JUANG
25     !-----    20070601 ELVMAX bug fix (*j*)
26     !
27     !   VERSION 4
28     !                ----- This code -----
29     !
30     !-----   MODIFIED TO IMPLEMENT THE ENHANCED LOW TROPOSPHERIC GRAVITY
31     !-----   WAVE DRAG DEVELOPED BY KIM AND ARAKAWA(JAS, 1995).
32     !        Orographic Std Dev (hprime), Convexity (OC), Asymmetry (OA4)
33     !        and Lx (CLX4) are input topographic statistics needed.
34     !
35     !-----   PROGRAMMED AND DEBUGGED BY HONG, ALPERT AND KIM --- JAN 1996.
36     !-----   debugged again - moorthi and iredell --- may 1998.
37     !-----
38     !       Further Cleanup, optimization and modification
39     !                                       - S. Moorthi May 98, March 99.
40     !-----   modified for usgs orography data (ncep office note 424)
41     !        and with several bugs fixed  - moorthi and hong --- july 1999.
42     !
43     !-----   Modified & implemented into NRL NOGAPS
44     !                                       - Young-Joon Kim, July 2000
45     !-----
46     !   VERSION lm MB  (6): oz fix 8/2003
47     !                ----- This code -----
48     !
49     !------   Changed to include the Lott and Miller Mtn Blocking
50     !         with some modifications by (*j*)  4/02
51     !        From a Principal Coordinate calculation using the
52     !        Hi Res 8 minute orography, the Angle of the
53     !        mtn with that to the East (x) axis is THETA, the slope
54     !        parameter SIGMA. The anisotropy is in GAMMA - all  are input
55     !        topographic statistics needed.  These are calculated off-line
56     !        as a function of model resolution in the fortran code ml01rg2.f,
57     !        with script mlb2.sh.   (*j*)
58     !-----   gwdps_mb.f version (following lmi) elvmax < hncrit (*j*)
59     !        MB3a expt to enhance elvmax mtn hgt see sigfac & hncrit
60     !        gwdps_GWDFIX_v6.f FIXGWD GF6.0 20070608 sigfac=4.
61     !-----
62     !----------------------------------------------------------------------C
63     !    USE
64     !        ROUTINE IS CALLED FROM GBPHYS  (AFTER CALL TO MONNIN)
65     !
66     !    PURPOSE
67     !        USING THE GWD PARAMETERIZATIONS OF PS-GLAS AND PH-
68     !        GFDL TECHNIQUE.  THE TIME TENDENCIES OF U V
69     !        ARE ALTERED TO INCLUDE THE EFFECT OF MOUNTAIN INDUCED
70     !        GRAVITY WAVE DRAG FROM SUB-GRID SCALE OROGRAPHY INCLUDING
71     !        CONVECTIVE BREAKING, SHEAR BREAKING AND THE PRESENCE OF
72     !        CRITICAL LEVELS
73     !
74     !  INPUT
75     !        A(IY,KM)  NON-LIN TENDENCY FOR V WIND COMPONENT
76     !        B(IY,KM)  NON-LIN TENDENCY FOR U WIND COMPONENT
77     !        U1(IX,KM) ZONAL WIND M/SEC  AT T0-DT
78     !        V1(IX,KM) MERIDIONAL WIND M/SEC AT T0-DT
79     !        T1(IX,KM) TEMPERATURE DEG K AT T0-DT
80     !        Q1(IX,KM) SPECIFIC HUMIDITY AT T0-DT
81     !
82     !        DELTIM  TIME STEP    SECS
83     !        SI(N)   P/PSFC AT BASE OF LAYER N
84     !        SL(N)   P/PSFC AT MIDDLE OF LAYER N
85     !        DEL(N)  POSITIVE INCREMENT OF P/PSFC ACROSS LAYER N
86     !        KPBL(IM) is the index of the top layer of the PBL
87     !        ipr & lprnt for diagnostics
88     !
89     !  OUTPUT
90     !        A, B    AS AUGMENTED BY TENDENCY DUE TO GWDPS
91     !                OTHER INPUT VARIABLES UNMODIFIED.
92     !   ********************************************************************
93           USE MACHINE , ONLY : kind_phys
94           implicit none
95           integer im, iy, ix, km, imx, lat, kdt, ipr, me
96           integer KPBL(IM)                 ! Index for the PBL top layer!
97           real(kind=kind_phys) deltim, G, CP, RD, RV,      cdmbgwd(2)
98     !     real(kind=kind_phys) A(IY,KM),    B(IY,KM),      PSTAR(IM)
99           real(kind=kind_phys) A(IY,KM),    B(IY,KM),
100          &                     U1(IX,KM),   V1(IX,KM),     T1(IX,KM),
101          &                     Q1(IX,KM),   PRSI(IX,KM+1), DEL(IX,KM),
102          &                     PRSL(IX,KM), PRSLK(IX,KM),  PHIL(IX,KM),
103          &                     PHII(IX,KM+1)
104           real(kind=kind_phys) OC(IM),     OA4(IY,4), CLX4(IY,4)
105          &,                    HPRIME(IM)
106     ! for lm mtn blocking
107           real(kind=kind_phys) ELVMAX(IM),THETA(IM),SIGMA(IM),GAMMA(IM)
108           real(kind=kind_phys) wk(IM)
109           real(kind=kind_phys) bnv2lm(IM,KM),PE(IM),EK(IM),ZBK(IM),UP(IM)
110           real(kind=kind_phys) DB(IM,KM),ANG(IM,KM),UDS(IM,KM)
111           real(kind=kind_phys) ZLEN, DBTMP, R, PHIANG, CDmb, DBIM
112     !
113     !     Some constants
114     !
115           real(kind=kind_phys) pi, dw2min, rimin, ric, bnv2min, efmin
116          &,                    efmax,hpmax,hpmin
117           PARAMETER (PI=3.1415926535897931)
118           PARAMETER (DW2MIN=1., RIMIN=-100., RIC=0.25, BNV2MIN=1.0E-5)
119     !     PARAMETER (EFMIN=0.0, EFMAX=10.0, hpmax=200.0)
120           PARAMETER (EFMIN=0.0, EFMAX=10.0, hpmax=2400.0, hpmin=1.0)
121     !
122           real(kind=kind_phys) FRC,    CE,     CEOFRC, frmax, CG, GMAX
123          &,                    CRITAC, VELEPS, FACTOP, RLOLEV, RDI
124           parameter (FRC=1.0, CE=0.8, CEOFRC=CE/FRC, frmax=100., CG=0.5)
125           parameter (GMAX=1.0, CRITAC=5.0E-4, VELEPS=1.0, FACTOP=0.5)
126           parameter (RLOLEV=50000.0) 
127     !     parameter (RLOLEV=500.0) 
128     !     parameter (RLOLEV=0.5)
129     !
130            real(kind=kind_phys) dpmin,hminmt,hncrit,minwnd,sigfac
131     ! --- for lm mtn blocking
132     !     parameter (cdmb = 1.0)     ! non-dim sub grid mtn drag Amp (*j*)
133           parameter (hncrit=8000.)   ! Max value in meters for ELVMAX (*j*)
134     !  hncrit set to 8000m and sigfac added to enhance elvmax mtn hgt
135           parameter (sigfac=4.0)     ! MB3a expt test for ELVMAX factor (*j*)
136           parameter (hminmt=50.)     ! min mtn height (*j*)
137           parameter (minwnd=0.1)     ! min wind component (*j*)
138     
139     !     parameter (dpmin=00.0)     ! Minimum thickness of the reference layer
140     !!    parameter (dpmin=05.0)     ! Minimum thickness of the reference layer
141     !     parameter (dpmin=20.0)     ! Minimum thickness of the reference layer
142                                      ! in centibars
143           parameter (dpmin=5000.0)   ! Minimum thickness of the reference layer
144                                      ! in Pa
145     !
146           real(kind=kind_phys) FDIR
147           integer mdir
148           parameter(mdir=8, FDIR=mdir/(PI+PI))
149           integer nwdir(mdir)
150           data nwdir/6,7,5,8,2,3,1,4/
151           save nwdir
152     !
153           LOGICAL ICRILV(IM)
154     !
155     !----   MOUNTAIN INDUCED GRAVITY WAVE DRAG
156     !
157           real(kind=kind_phys) TAUB(IM),  XN(IM),     YN(IM),    UBAR(IM)
158          &,                    VBAR(IM),  ULOW(IM),   OA(IM),    CLX(IM)
159          &,                    ROLL(IM),  ULOI(IM),   DUSFC(IM), DVSFC(IM)
160          &,                    DTFAC(IM), XLINV(IM),  DELKS(IM), DELKS1(IM)
161     !
162           real(kind=kind_phys) BNV2(IM,KM),  TAUP(IM,KM+1), ri_n(IM,KM) 
163          &,                    TAUD(IM,KM),  RO(IM,KM),     VTK(IM,KM)
164          &,                    VTJ(IM,KM),   SCOR(IM),      VELCO(IM,KM-1)
165          &,                    bnv2bar(im)
166     !
167           real(kind=kind_phys) VELKO(KM-1)
168           Integer   kref(IM), kint(im), iwk(im), iwk2(im), ipt(im)
169     ! for lm mtn blocking
170           Integer   kreflm(IM), iwklm(im), iptlm(im)
171           Integer   idxzb(im), idxm1, ktrial, klevm1, nmtvr
172     !
173           real(kind=kind_phys) gor,    gocp,  fv,    xl,    gr2,  bnv,  fr
174          &,                    brvf,   cleff, tem,   tem1,  tem2, temc, temv
175          &,                    wdir,   ti,    rdz,   dw2,   shr2, bvf2
176          &,                    rdelks, wtkbj, efact, coefm, gfobnv
177          &,                    scork,  rscor, hd,    fro,   rim,  sira
178          &,                    dtaux,  dtauy, pkp1log, pklog
179           integer ncnt, kmm1, kmm2, lcap, lcapp1, kbps, kbpsp1,kbpsm1
180          &, kmps, kmpsp1, idir, nwd, i, j, k, klcap, kp1, kmpbl, npt, npr
181          &, kmll,kmds
182     !    &, kmll,kmds,ihit,jhit
183           logical lprnt
184     !
185     !     parameter (cdmb = 1.0)     ! non-dim sub grid mtn drag Amp (*j*)
186     ! non-dim sub grid mtn drag Amp (*j*)
187     !     cdmb = 1.0/float(IMX/192)
188     !     cdmb = 192.0/float(IMX)
189           cdmb = 4.0 * 192.0/float(IMX)
190           if (cdmbgwd(1) >= 0.0) cdmb = cdmb * cdmbgwd(1)
191     !
192           npr = 0
193           DO I = 1, IM
194              DUSFC(I) = 0.
195              DVSFC(I) = 0.
196           ENDDO
197     !
198           DO K = 1, KM
199             DO I = 1, IM
200               DB(I,K)  = 0.
201               ANG(I,K) = 0.
202               UDS(I,K) = 0.
203             ENDDO
204           ENDDO
205     !
206           RDI  = 1.0 / RD
207           GOR  = G/RD
208           GR2  = G*GOR
209           GOCP = G/CP
210           FV   = RV/RD - 1
211     !
212     !     NCNT   = 0
213           KMM1   = KM - 1
214           KMM2   = KM - 2
215           LCAP   = KM
216           LCAPP1 = LCAP + 1
217     !
218     !
219           IF ( NMTVR .eq. 14) then 
220     ! ----  for lm and gwd calculation points
221             ipt = 0
222             npt = 0
223             DO I = 1,IM
224               IF ( (elvmax(i) .GT. HMINMT) 
225          &       .and. (hprime(i) .GT. hpmin) )  then
226                  npt      = npt + 1
227                  ipt(npt) = i
228                  if (ipr .eq. i) npr = npt
229               ENDIF
230             ENDDO
231             IF (npt .eq. 0) RETURN     ! No gwd/mb calculation done!
232     !
233     !       if (lprnt) print *,' npt=',npt,' npr=',npr,' ipr=',ipr,' im=',im
234     !    &,' ipt(npt)=',ipt(npt)
235     !
236     ! --- iwklm is the level above the height of the of the mountain.
237     ! --- idxzb is the level of the dividing streamline.
238     ! INITIALIZE DIVIDING STREAMLINE (DS) CONTROL VECTOR
239     !
240             do i=1,npt
241               iwklm(i)  = 2
242               IDXZB(i)  = 0 
243               kreflm(i) = 0
244             enddo
245     !       if (lprnt) 
246     !    &  print *,' in gwdps_lm.f npt,IM,IX,IY,km,me=',npt,IM,IX,IY,km,me
247     !
248     !
249     ! start lm mtn blocking (mb) section
250     !
251     !..............................
252     !..............................
253     !
254     !  (*j*)  11/03:  test upper limit on KMLL=km - 1
255     !      then do not need hncrit -- test with large hncrit first.
256     !       KMLL  = km / 2 ! maximum mtnlm height : # of vertical levels / 2
257             KMLL = kmm1
258     ! --- No mtn should be as high as KMLL (so we do not have to start at 
259     ! --- the top of the model but could do calc for all levels).
260     !
261               DO I = 1, npt
262                 j = ipt(i)
263                 ELVMAX(J) = min (ELVMAX(J) + sigfac * hprime(j), hncrit)
264               ENDDO
265     !
266             DO K = 1,KMLL
267               DO I = 1, npt
268                 j = ipt(i)
269     ! --- interpolate to max mtn height for index, iwklm(I) wk[gz]
270     ! --- ELVMAX is limited to hncrit because to hi res topo30 orog.
271                 pkp1log =  phil(j,k+1) / G
272                 pklog =  phil(j,k)   / G
273     !!!-------     ELVMAX(J) = min (ELVMAX(J) + sigfac * hprime(j), hncrit)
274                 if ( ( ELVMAX(j) .le.  pkp1log ) .and. 
275          &           ( ELVMAX(j) .ge.   pklog  ) ) THEN
276     !     print *,' in gwdps_lm.f 1  =',k,ELVMAX(j),pklog,pkp1log,me
277     ! ---        wk for diags but can be saved and reused.  
278                    wk(i)  = G * ELVMAX(j) / ( phil(j,k+1) - phil(j,k) )
279                    iwklm(I)  =  MAX(iwklm(I), k+1 ) 
280     !     print *,' in gwdps_lm.f 2 npt=',npt,i,j,wk(i),iwklm(i),me
281                 endif
282     !
283     ! ---        find at prsl levels large scale environment variables
284     ! ---        these cover all possible mtn max heights
285                 VTJ(I,K)  = T1(J,K)  * (1.+FV*Q1(J,K))
286                 VTK(I,K)  = VTJ(I,K) / PRSLK(J,K)
287                 RO(I,K)   = RDI * PRSL(J,K) / VTJ(I,K) ! DENSITY TONS/M**3
288               ENDDO
289             ENDDO
290     !
291     ! testing for highest model level of mountain top
292     !
293     !         ihit = 2
294     !         jhit = 0
295     !        do i = 1, npt
296     !        j=ipt(i)
297     !          if ( iwklm(i) .gt. ihit ) then 
298     !            ihit = iwklm(i)
299     !            jhit = j
300     !          endif
301     !        enddo
302     !     print *, ' mb: kdt,max(iwklm),jhit,phil,me=',
303     !    &          kdt,ihit,jhit,phil(jhit,ihit),me
304              
305             klevm1 = KMLL - 1
306             DO K = 1, klevm1  
307               DO I = 1, npt
308                j   = ipt(i)
309                 RDZ  = g   / ( phil(j,k+1) - phil(j,k) )
310     ! ---                               Brunt-Vaisala Frequency
311                 BNV2LM(I,K) = (G+G) * RDZ * ( VTK(I,K+1)-VTK(I,K) )
312          &                     / ( VTK(I,K+1)+VTK(I,K) )
313                 bnv2lm(i,k) = max( bnv2lm(i,k), bnv2min )
314               ENDDO
315             ENDDO
316     !    print *,' in gwdps_lm.f 3 npt=',npt,j,RDZ,me
317     !
318             DO I = 1, npt
319               J   = ipt(i)
320               DELKS(I)  = 1.0 / (PRSI(J,1) - PRSI(J,iwklm(i)))
321               DELKS1(I) = 1.0 / (PRSL(J,1) - PRSL(J,iwklm(i)))
322               UBAR (I)  = 0.0
323               VBAR (I)  = 0.0
324               ROLL (I)  = 0.0
325               PE   (I)  = 0.0
326               EK   (I)  = 0.0
327               BNV2bar(I) = (PRSL(J,1)-PRSL(J,2)) * DELKS1(I) * BNV2LM(I,1)
328             ENDDO
329     
330     ! --- find the dividing stream line height 
331     ! --- starting from the level above the max mtn downward
332     ! --- iwklm(i) is the k-index of mtn elvmax elevation
333             DO Ktrial = KMLL, 1, -1
334               DO I = 1, npt
335                  IF ( Ktrial .LT. iwklm(I) .and. kreflm(I) .eq. 0 ) then
336                     kreflm(I) = Ktrial
337                  ENDIF
338               ENDDO
339             ENDDO
340     !     print *,' in gwdps_lm.f 4 npt=',npt,kreflm(npt),me
341     !
342     ! --- in the layer kreflm(I) to 1 find PE (which needs N, ELVMAX)
343     ! ---  make averages, guess dividing stream (DS) line layer.
344     ! ---  This is not used in the first cut except for testing and
345     ! --- is the vert ave of quantities from the surface to mtn top.
346     !   
347             DO I = 1, npt
348               DO K = 1, Kreflm(I)
349                 J        = ipt(i)
350                 RDELKS     = DEL(J,K) * DELKS(I)
351                 UBAR(I)    = UBAR(I)  + RDELKS * U1(J,K) ! trial Mean U below 
352                 VBAR(I)    = VBAR(I)  + RDELKS * V1(J,K) ! trial Mean V below 
353                 ROLL(I)    = ROLL(I)  + RDELKS * RO(I,K) ! trial Mean RO below 
354                 RDELKS     = (PRSL(J,K)-PRSL(J,K+1)) * DELKS1(I)
355                 BNV2bar(I) = BNV2bar(I) + BNV2lm(I,K) * RDELKS
356     ! --- these vert ave are for diags, testing and GWD to follow (*j*).
357               ENDDO
358             ENDDO
359     !     print *,' in gwdps_lm.f 5  =',i,kreflm(npt),BNV2bar(npt),me
360     !
361     ! --- integrate to get PE in the trial layer.
362     ! --- Need the first layer where PE>EK - as soon as 
363     ! --- IDXZB is not 0 we have a hit and Zb is found.
364     !
365             DO I = 1, npt
366               J = ipt(i)
367               DO K = iwklm(I), 1, -1
368                 PHIANG   =  atan2D(V1(J,K),U1(J,K))
369                 ANG(I,K) = ( THETA(J) - PHIANG )
370                 if ( ANG(I,K) .gt.  90. ) ANG(I,K) = ANG(I,K) - 180.
371                 if ( ANG(I,K) .lt. -90. ) ANG(I,K) = ANG(I,K) + 180.
372     !
373                 UDS(I,K) = 
374          &          MAX(SQRT(U1(J,K)*U1(J,K) + V1(J,K)*V1(J,K)), minwnd)
375     ! --- Test to see if we found Zb previously
376                 IF (IDXZB(I) .eq. 0 ) then
377                   PE(I) = PE(I) + BNV2lm(I,K) * 
378          &           ( G * ELVMAX(J) - phil(J,K) ) * 
379          &           ( PHII(J,K+1) - PHII(J,K) ) / (G*G)
380     ! --- KE
381     ! --- Wind projected on the line perpendicular to mtn range, U(Zb(K)).
382     ! --- kenetic energy is at the layer Zb
383     ! --- THETA ranges from -+90deg |_ to the mtn "largest topo variations"
384                   UP(I)  =  UDS(I,K) * cosD(ANG(I,K))
385                   EK(I)  = 0.5 *  UP(I) * UP(I) 
386     
387     ! --- Dividing Stream lime  is found when PE =exceeds EK.
388                   IF ( PE(I) .ge.  EK(I) ) IDXZB(I) = K
389     ! --- Then mtn blocked flow is between Zb=k(IDXZB(I)) and surface
390     !
391                 ENDIF
392               ENDDO
393             ENDDO
394     !
395     !     print *,' in gwdps_lm.f 6  =',phiang,THETA(ipt(npt)),me
396     !     print *,' in gwdps_lm.f 7  =',IDXZB(npt),PE(npt)
397     !
398     !     if (lprnt .and. npr .gt. 0) then
399     !       print *,' BNV2bar,BNV2lm=',bnv2bar(npr),BNV2lm(npr,1:klevm1)
400     !       print *,' npr,IDXZB,UDS=',npr,IDXZB(npr),UDS(npr,:)
401     !       print *,' PE,UP,EK=',PE(npr),UP(npr),EK(npr)
402     !     endif
403     !
404             DO I = 1, npt
405               J    = ipt(i)
406     ! --- Calc if N constant in layers (Zb guess) - a diagnostic only.
407               ZBK(I) =  ELVMAX(J) - SQRT(UBAR(I)**2 + VBAR(I)**2)/BNV2bar(I)
408             ENDDO
409     !
410     !     if (lprnt .and. npr .gt. 0) then
411     !       print *,' iwklm,ZBK=',iwklm(npr),ZBK(npr),IDXZB(npr)
412     !       print *,' Zb=',PHIL(ipr),IDXZB(npr))/G
413     !     print *,' in gwdps_lm.f 8 npt =',npt,ZBK(npt),UP(npt),me
414     !     endif
415     !
416     ! --- The drag for mtn blocked flow
417     ! 
418             DO I = 1, npt
419               J = ipt(i)
420               ZLEN = 0.
421     !      print *,' in gwdps_lm.f 9  =',i,j,IDXZB(i),me
422               IF ( IDXZB(I) .gt. 0 ) then 
423                 DO K = IDXZB(I), 1, -1
424                   IF ( PHIL(J,IDXZB(I)) .gt.  PHIL(J,K) ) then
425                     ZLEN = SQRT( ( PHIL(J,IDXZB(I)) - PHIL(J,K) ) / 
426          &                       ( PHIL(J,K ) + G * hprime(J) ) )
427     ! --- lm eq 14:
428                     R = (cosD(ANG(I,K))**2 + GAMMA(J) * sinD(ANG(I,K))**2) / 
429          &              (gamma(J) * cosD(ANG(I,K))**2 + sinD(ANG(I,K))**2)
430     ! --- (negitive of DB -- see sign at tendency)
431                     DBTMP = 0.25 *  CDmb * 
432          &                  MAX( 2. - 1. / R, 0. ) * sigma(J) * 
433          &                  MAX(cosD(ANG(I,K)), gamma(J)*sinD(ANG(I,K))) *
434          &                  ZLEN / hprime(J) 
435                     DB(I,K) =  DBTMP * UDS(I,K)    
436     !
437     !               if(lprnt .and. i .eq. npr) then 
438     !                 print *,' in gwdps_lmi.f 10 npt=',npt,i,j,idxzb(i)
439     !    &,           DBTMP,R' ang=',ang(i,k),' gamma=',gamma(j),' K=',K
440     !                 print *,' in gwdps_lmi.f 11   K=',k,ZLEN,cosD(ANG(I,K))
441     !                 print *,' in gwdps_lmi.f 12  DB=',DB(i,k),sinD(ANG(I,K))
442     !               endif
443                   endif
444                 ENDDO
445     !         if(lprnt) print *,' @K=1,ZLEN,DBTMP=',K,ZLEN,DBTMP
446               endif
447             ENDDO
448     ! 
449     !.............................
450     !.............................
451     ! end  mtn blocking section
452     !
453           ELSEIF ( NMTVR .ne. 14) then 
454     ! ----  for mb not present and  gwd (nmtvr .ne .14) 
455             ipt     = 0
456             npt     = 0
457             DO I = 1,IM
458               IF ( hprime(i) .GT. hpmin )  then
459                  npt      = npt + 1
460                  ipt(npt) = i
461                  if (ipr .eq. i) npr = npt
462               ENDIF
463             ENDDO
464             IF (npt .eq. 0) RETURN     ! No gwd/mb calculation done!
465     !
466     !       if (lprnt) print *,' NPR=',npr,' npt=',npt,' IPR=',IPR
467     !      &,' ipt(npt)=',ipt(npt)
468     !
469             do i=1,npt
470               IDXZB(i) = 0
471             enddo
472           ENDIF
473     !
474     !.............................
475     !.............................
476     !
477           KMPBL  = km / 2 ! maximum pbl height : # of vertical levels / 2
478     !
479     !  Scale cleff between IM=384*2 and 192*2 for T126/T170 and T62
480     !
481           if (imx .gt. 0) then
482     !       cleff = 1.0E-5 * SQRT(FLOAT(IMX)/384.0) !  this is inverse of CLEFF!
483     !       cleff = 1.0E-5 * SQRT(FLOAT(IMX)/192.0) !  this is inverse of CLEFF!
484     !       cleff = 0.5E-5 * SQRT(FLOAT(IMX)/192.0) !  this is inverse of CLEFF!
485     !       cleff = 1.0E-5 * SQRT(FLOAT(IMX)/192)/float(IMX/192)
486     !       cleff = 1.0E-5 / SQRT(FLOAT(IMX)/192.0) !  this is inverse of CLEFF!
487             cleff = 0.5E-5 / SQRT(FLOAT(IMX)/192.0) !  this is inverse of CLEFF!
488     !       cleff = 2.0E-5 * SQRT(FLOAT(IMX)/192.0) !  this is inverse of CLEFF!
489     !       cleff = 2.5E-5 * SQRT(FLOAT(IMX)/192.0) !  this is inverse of CLEFF!
490           endif
491           if (cdmbgwd(2) >= 0.0) cleff = cleff * cdmbgwd(2)
492     !
493           DO K = 1,KM
494             DO I =1,npt
495               J         = ipt(i)
496               VTJ(I,K)  = T1(J,K)  * (1.+FV*Q1(J,K))
497               VTK(I,K)  = VTJ(I,K) / PRSLK(J,K)
498               RO(I,K)   = RDI * PRSL(J,K) / VTJ(I,K) ! DENSITY TONS/M**3
499               TAUP(I,K) = 0.0
500             ENDDO
501           ENDDO
502           DO K = 1,KMM1
503             DO I =1,npt
504               J         = ipt(i)
505               TI        = 2.0 / (T1(J,K)+T1(J,K+1))
506               TEM       = TI  / (PRSL(J,K)-PRSL(J,K+1))
507               RDZ       = g   / (phil(j,k+1) - phil(j,k))
508               TEM1      = U1(J,K) - U1(J,K+1)
509               TEM2      = V1(J,K) - V1(J,K+1)
510               DW2       = TEM1*TEM1 + TEM2*TEM2
511               SHR2      = MAX(DW2,DW2MIN) * RDZ * RDZ
512               BVF2      = G*(GOCP+RDZ*(VTJ(I,K+1)-VTJ(I,K))) * TI
513               ri_n(I,K) = MAX(BVF2/SHR2,RIMIN)   ! Richardson number
514     !                                              Brunt-Vaisala Frequency
515     !         TEM       = GR2 * (PRSL(J,K)+PRSL(J,K+1)) * TEM
516     !         BNV2(I,K) = TEM * (VTK(I,K+1)-VTK(I,K))/(VTK(I,K+1)+VTK(I,K))
517               BNV2(I,K) = (G+G) * RDZ * (VTK(I,K+1)-VTK(I,K))
518          &                            / (VTK(I,K+1)+VTK(I,K))
519               bnv2(i,k) = max( bnv2(i,k), bnv2min )
520             ENDDO
521           ENDDO
522     !      print *,' in gwdps_lm.f GWD:14  =',npt,kmm1,bnv2(npt,kmm1)
523     !
524     !     Apply 3 point smoothing on BNV2
525     !
526     !     do k=1,km
527     !       do i=1,im
528     !         vtk(i,k) = bnv2(i,k)
529     !       enddo
530     !     enddo
531     !     do k=2,kmm1
532     !       do i=1,im
533     !         bnv2(i,k) = 0.25*(vtk(i,k-1)+vtk(i,k+1)) + 0.5*vtk(i,k)
534     !       enddo
535     !     enddo
536     !
537     !     Finding the first interface index above 50 hPa level
538     !
539           do i=1,npt
540             iwk(i) = 2
541           enddo
542           DO K=3,KMPBL
543             DO I=1,npt
544               j   = ipt(i)
545               tem = (prsi(j,1) - prsi(j,k))
546               if (tem .lt. dpmin) iwk(i) = k
547             enddo
548           enddo
549     !
550           KBPS = 1
551           KMPS = KM
552           DO I=1,npt
553             J         = ipt(i)
554             kref(I)   = MAX(IWK(I), KPBL(J)+1 ) ! reference level 
555             DELKS(I)  = 1.0 / (PRSI(J,1) - PRSI(J,kref(I)))
556             DELKS1(I) = 1.0 / (PRSL(J,1) - PRSL(J,kref(I)))
557             UBAR (I)  = 0.0
558             VBAR (I)  = 0.0
559             ROLL (I)  = 0.0
560             KBPS      = MAX(KBPS,  kref(I))
561             KMPS      = MIN(KMPS,  kref(I))
562     !
563             BNV2bar(I) = (PRSL(J,1)-PRSL(J,2)) * DELKS1(I) * BNV2(I,1)
564           ENDDO
565     !      print *,' in gwdps_lm.f GWD:15  =',KBPS,KMPS
566           KBPSP1 = KBPS + 1
567           KBPSM1 = KBPS - 1
568           DO K = 1,KBPS
569             DO I = 1,npt
570               IF (K .LT. kref(I)) THEN
571                 J          = ipt(i)
572                 RDELKS     = DEL(J,K) * DELKS(I)
573                 UBAR(I)    = UBAR(I)  + RDELKS * U1(J,K)   ! Mean U below kref
574                 VBAR(I)    = VBAR(I)  + RDELKS * V1(J,K)   ! Mean V below kref
575     !
576                 ROLL(I)    = ROLL(I)  + RDELKS * RO(I,K)   ! Mean RO below kref
577                 RDELKS     = (PRSL(J,K)-PRSL(J,K+1)) * DELKS1(I)
578                 BNV2bar(I) = BNV2bar(I) + BNV2(I,K) * RDELKS
579               ENDIF
580             ENDDO
581           ENDDO
582     !      print *,' in gwdps_lm.f GWD:15B =',bnv2bar(npt)
583     !
584     !     FIGURE OUT LOW-LEVEL HORIZONTAL WIND DIRECTION AND FIND 'OA'
585     !
586     !             NWD  1   2   3   4   5   6   7   8
587     !              WD  W   S  SW  NW   E   N  NE  SE
588     !
589           DO I = 1,npt
590             J      = ipt(i)
591             wdir   = atan2(UBAR(I),VBAR(I)) + pi
592             idir   = mod(nint(fdir*wdir),mdir) + 1
593             nwd    = nwdir(idir)
594             OA(I)  = (1-2*INT( (NWD-1)/4 )) * OA4(J,MOD(NWD-1,4)+1)
595             CLX(I) = CLX4(J,MOD(NWD-1,4)+1)
596           ENDDO
597     !
598     !-----XN,YN            "LOW-LEVEL" WIND PROJECTIONS IN ZONAL
599     !                                    & MERIDIONAL DIRECTIONS
600     !-----ULOW             "LOW-LEVEL" WIND MAGNITUDE -        (= U)
601     !-----BNV2             BNV2 = N**2
602     !-----TAUB             BASE MOMENTUM FLUX
603     !-----= -(RO * U**3/(N*XL)*GF(FR) FOR N**2 > 0
604     !-----= 0.                        FOR N**2 < 0
605     !-----FR               FROUDE    =   N*HPRIME / U
606     !-----G                GMAX*FR**2/(FR**2+CG/OC)
607     !
608     !-----INITIALIZE SOME ARRAYS
609     !
610           DO I = 1,npt
611             XN(I)     = 0.0
612             YN(I)     = 0.0
613             TAUB (I)  = 0.0
614             ULOW (I)  = 0.0
615             DTFAC(I)  = 1.0
616             ICRILV(I) = .FALSE. ! INITIALIZE CRITICAL LEVEL CONTROL VECTOR
617             
618     !
619     !----COMPUTE THE "LOW LEVEL" WIND MAGNITUDE (M/S)
620     !
621             ULOW(I) = MAX(SQRT(UBAR(I)*UBAR(I) + VBAR(I)*VBAR(I)), 1.0)
622             ULOI(I) = 1.0 / ULOW(I)
623           ENDDO
624     !
625           DO  K = 1,KMM1
626             DO  I = 1,npt
627               J            = ipt(i)
628               VELCO(I,K)   = 0.5 * ((U1(J,K)+U1(J,K+1))*UBAR(I)
629          &                       +  (V1(J,K)+V1(J,K+1))*VBAR(I))
630               VELCO(I,K)   = VELCO(I,K) * ULOI(I)
631     !         IF ((VELCO(I,K).LT.VELEPS) .AND. (VELCO(I,K).GT.0.)) THEN
632     !           VELCO(I,K) = VELEPS
633     !         ENDIF
634             ENDDO
635           ENDDO
636     !      
637     !
638     !   find the interface level of the projected wind where
639     !   low levels & upper levels meet above pbl
640     !
641           do i=1,npt
642             kint(i) = km
643           enddo
644           do k = 1,kmm1
645             do i = 1,npt
646               IF (K .GT. kref(I)) THEN
647                 if(velco(i,k) .lt. veleps .and. kint(i) .eq. km) then
648                   kint(i) = k+1
649                 endif
650               endif
651             enddo
652           enddo
653     !  WARNING  KINT = KREF !!!!!!!!!
654           do i=1,npt
655             kint(i) = kref(i)
656           enddo
657     !
658     !     if(lprnt) print *,' ubar=',ubar
659     !    &,' vbar=',vbar,' ulow=',ulow,' veleps=',veleps
660     !
661           DO I = 1,npt
662             J      = ipt(i)
663             BNV    = SQRT( BNV2bar(I) )
664             FR     = BNV     * ULOI(I) * min(HPRIME(J),hpmax)
665             FR     = MIN(FR, FRMAX)
666             XN(I)  = UBAR(I) * ULOI(I)
667             YN(I)  = VBAR(I) * ULOI(I)
668     !
669     !     Compute the base level stress and store it in TAUB
670     !     CALCULATE ENHANCEMENT FACTOR, NUMBER OF MOUNTAINS & ASPECT
671     !     RATIO CONST. USE SIMPLIFIED RELATIONSHIP BETWEEN STANDARD
672     !     DEVIATION & CRITICAL HGT
673     !
674             EFACT    = (OA(I) + 2.) ** (CEOFRC*FR)
675             EFACT    = MIN( MAX(EFACT,EFMIN), EFMAX )
676     !
677             COEFM    = (1. + CLX(I)) ** (OA(I)+1.)
678     !
679             XLINV(I) = COEFM * CLEFF
680     !
681             TEM      = FR    * FR * OC(J)
682             GFOBNV   = GMAX  * TEM / ((TEM + CG)*BNV)  ! G/N0
683     !
684             TAUB(I)  = XLINV(I) * ROLL(I) * ULOW(I) * ULOW(I)
685          &           * ULOW(I)  * GFOBNV  * EFACT         ! BASE FLUX Tau0
686     !
687     !         tem      = min(HPRIME(I),hpmax)
688     !         TAUB(I)  = XLINV(I) * ROLL(I) * ULOW(I) * BNV * tem * tem
689     !
690             K        = MAX(1, kref(I)-1)
691             TEM      = MAX(VELCO(I,K)*VELCO(I,K), 0.1)
692             SCOR(I)  = BNV2(I,K) / TEM  ! Scorer parameter below ref level
693           ENDDO
694     !     if(lprnt) print *,' taub=',taub
695     !                                                                       
696     !----SET UP BOTTOM VALUES OF STRESS
697     !
698           DO K = 1, KBPS
699             DO I = 1,npt
700               IF (K .LE. kref(I)) TAUP(I,K) = TAUB(I)
701             ENDDO
702           ENDDO
703     !
704     !   Now compute vertical structure of the stress.
705     !
706           DO K = KMPS, KMM1                   ! Vertical Level K Loop!
707             KP1 = K + 1
708             DO I = 1, npt
709     !
710     !-----UNSTABLE LAYER IF RI < RIC
711     !-----UNSTABLE LAYER IF UPPER AIR VEL COMP ALONG SURF VEL <=0 (CRIT LAY)
712     !---- AT (U-C)=0. CRIT LAYER EXISTS AND BIT VECTOR SHOULD BE SET (.LE.)
713     !
714               IF (K .GE. kref(I)) THEN
715                 ICRILV(I) = ICRILV(I) .OR. ( ri_n(I,K) .LT. RIC)
716          &                            .OR. (VELCO(I,K) .LE. 0.0)
717               ENDIF
718             ENDDO
719     !
720             DO I = 1,npt
721               IF (K .GE. kref(I))   THEN
722                 IF (.NOT.ICRILV(I) .AND. TAUP(I,K) .GT. 0.0 ) THEN
723                   TEMV = 1.0 / max(VELCO(I,K), 0.01)
724     !             IF (OA(I) .GT. 0. .AND.  PRSI(ipt(i),KP1).GT.RLOLEV) THEN
725                   IF (OA(I).GT.0. .AND. kp1 .lt. kint(i)) THEN
726                     SCORK   = BNV2(I,K) * TEMV * TEMV
727                     RSCOR   = MIN(1.0, SCORK / SCOR(I))
728                     SCOR(I) = SCORK
729                   ELSE 
730                     RSCOR   = 1.
731                   ENDIF
732     !
733                   BRVF = SQRT(BNV2(I,K))        ! Brunt-Vaisala Frequency
734     !             TEM1 = XLINV(I)*(RO(I,KP1)+RO(I,K))*BRVF*VELCO(I,K)*0.5
735                   TEM1 = XLINV(I)*(RO(I,KP1)+RO(I,K))*BRVF*0.5
736          &                       * max(VELCO(I,K),0.01)
737                   HD   = SQRT(TAUP(I,K) / TEM1)
738                   FRO  = BRVF * HD * TEMV
739     !
740     !    RIM is the  MINIMUM-RICHARDSON NUMBER BY SHUTTS (1985)
741     !
742                   TEM2   = SQRT(ri_n(I,K))
743                   TEM    = 1. + TEM2 * FRO
744                   RIM    = ri_n(I,K) * (1.-FRO) / (TEM * TEM)
745     !
746     !    CHECK STABILITY TO EMPLOY THE 'SATURATION HYPOTHESIS'
747     !    OF LINDZEN (1981) EXCEPT AT TROPOSPHERIC DOWNSTREAM REGIONS
748     !
749     !                                       ----------------------
750                   IF (RIM .LE. RIC .AND.
751     !    &           (OA(I) .LE. 0. .OR.  PRSI(ipt(I),KP1).LE.RLOLEV )) THEN
752          &           (OA(I) .LE. 0. .OR.  kp1 .ge. kint(i) )) THEN
753                      TEMC = 2.0 + 1.0 / TEM2
754                      HD   = VELCO(I,K) * (2.*SQRT(TEMC)-TEMC) / BRVF
755                      TAUP(I,KP1) = TEM1 * HD * HD
756                   ELSE 
757                      TAUP(I,KP1) = TAUP(I,K) * RSCOR
758                   ENDIF
759                   taup(i,kp1) = min(taup(i,kp1), taup(i,k))
760                 ENDIF
761               ENDIF
762             ENDDO
763           ENDDO
764     !
765     !     DO I=1,IM
766     !       taup(i,km+1) = taup(i,km)
767     !     ENDDO
768     !
769           IF(LCAP .LE. KM) THEN
770              DO KLCAP = LCAPP1, KM+1
771                 DO I = 1,npt
772                   SIRA          = PRSI(ipt(I),KLCAP) / PRSI(ipt(I),LCAP)
773                   TAUP(I,KLCAP) = SIRA * TAUP(I,LCAP)
774                 ENDDO
775              ENDDO
776           ENDIF
777     !
778     !     Calculate - (g/p*)*d(tau)/d(sigma) and Decel terms DTAUX, DTAUY
779     !
780           DO K = 1,KM
781             DO I = 1,npt
782               TAUD(I,K) = G * (TAUP(I,K+1) - TAUP(I,K)) / DEL(ipt(I),K)
783             ENDDO
784           ENDDO
785     !
786     !------LIMIT DE-ACCELERATION (MOMENTUM DEPOSITION ) AT TOP TO 1/2 VALUE
787     !------THE IDEA IS SOME STUFF MUST GO OUT THE 'TOP'
788     !
789           DO KLCAP = LCAP, KM
790              DO I = 1,npt
791                 TAUD(I,KLCAP) = TAUD(I,KLCAP) * FACTOP
792              ENDDO
793           ENDDO
794     !
795     !------IF THE GRAVITY WAVE DRAG WOULD FORCE A CRITICAL LINE IN THE
796     !------LAYERS BELOW SIGMA=RLOLEV DURING THE NEXT DELTIM TIMESTEP,
797     !------THEN ONLY APPLY DRAG UNTIL THAT CRITICAL LINE IS REACHED.
798     !
799           DO K = 1,KMM1
800             DO I = 1,npt
801                IF (K .GT. kref(I) .and. PRSI(ipt(i),K) .GE. RLOLEV) THEN
802                  IF(TAUD(I,K).NE.0.) THEN
803                    TEM = DELTIM * TAUD(I,K)
804                    DTFAC(I) = MIN(DTFAC(I),ABS(VELCO(I,K)/TEM))
805                  ENDIF
806                ENDIF
807             ENDDO
808           ENDDO
809     !
810     !     if(lprnt .and. npr .gt. 0) then
811     !       print *,' before  A=',A(npr,:)
812     !       print *,' before  B=',B(npr,:)
813     !     endif
814     
815           DO K = 1,KM
816             DO I = 1,npt
817               J          = ipt(i)
818               TAUD(I,K)  = TAUD(I,K) * DTFAC(I)
819               DTAUX      = TAUD(I,K) * XN(I)
820               DTAUY      = TAUD(I,K) * YN(I)
821     ! ---  lm mb (*j*)  changes overwrite GWD
822               if ( K .lt. IDXZB(I) .AND. IDXZB(I) .ne. 0 ) then
823                 DBIM = DB(I,K) / (1.+DB(I,K)*DELTIM)
824                 A(J,K)  = - DBIM * V1(J,K) + A(J,K)
825                 B(J,K)  = - DBIM * U1(J,K) + B(J,K)
826     !          if ( ABS(DBIM * U1(J,K)) .gt. .01 ) 
827     !    & print *,' in gwdps_lmi.f KDT=',KDT,I,K,DB(I,K),
828     !    &                      dbim,idxzb(I),U1(J,K),V1(J,K),me
829                 DUSFC(J)   = DUSFC(J) - DBIM * V1(J,K) * DEL(J,K)
830                 DVSFC(J)   = DVSFC(J) - DBIM * U1(J,K) * DEL(J,K)
831               else
832     !
833                 A(J,K)     = DTAUY     + A(J,K)
834                 B(J,K)     = DTAUX     + B(J,K)
835                 DUSFC(J)   = DUSFC(J)  + DTAUX * DEL(J,K)
836                 DVSFC(J)   = DVSFC(J)  + DTAUY * DEL(J,K)
837               endif
838             ENDDO
839           ENDDO
840     !     if (lprnt) then
841     !       print *,' in gwdps_lm.f after  A=',A(ipr,:)
842     !       print *,' in gwdps_lm.f after  B=',B(ipr,:)
843     !       print *,' DB=',DB(ipr,:)
844     !     endif
845           TEM    = -1.0/G
846           DO I = 1,npt
847             J          = ipt(i)
848     !       TEM    = (-1.E3/G)
849             DUSFC(J) = TEM * DUSFC(J)
850             DVSFC(J) = TEM * DVSFC(J)
851           ENDDO
852     !                                                                       
853     !    MONITOR FOR EXCESSIVE GRAVITY WAVE DRAG TENDENCIES IF NCNT>0
854     !
855     !     IF(NCNT.GT.0) THEN
856     !        IF(LAT.GE.38.AND.LAT.LE.42) THEN
857     !CMIC$ GUARD 37
858     !           DO 92 I = 1,IM
859     !              IF(IKOUNT.GT.NCNT) GO TO 92
860     !              IF(I.LT.319.OR.I.GT.320) GO TO 92
861     !              DO 91 K = 1,KM
862     !                 IF(ABS(TAUD(I,K)) .GT. CRITAC) THEN
863     !                    IF(I.LE.IM) THEN
864     !                       IKOUNT = IKOUNT+1
865     !                       PRINT 123,I,LAT,KDT
866     !                       PRINT 124,TAUB(I),BNV(I),ULOW(I),
867     !    1                  GF(I),FR(I),ROLL(I),HPRIME(I),XN(I),YN(I)
868     !                       PRINT 124,(TAUD(I,KK),KK = 1,KM)
869     !                       PRINT 124,(TAUP(I,KK),KK = 1,KM+1)
870     !                       PRINT 124,(ri_n(I,KK),KK = 1,KM)
871     !                       DO 93 KK = 1,KMM1
872     !                          VELKO(KK) =
873     !    1                  0.5*((U1(I,KK)+U1(I,KK+1))*UBAR(I)+
874     !    2                  (V1(I,KK)+V1(I,KK+1))*VBAR(I))*ULOI(I)
875     !93                     CONTINUE
876     !                       PRINT 124,(VELKO(KK),KK = 1,KMM1)
877     !                       PRINT 124,(A    (I,KK),KK = 1,KM)
878     !                       PRINT 124,(DTAUY(I,KK),KK = 1,KM)
879     !                       PRINT 124,(B    (I,KK),KK = 1,KM)
880     !                       PRINT 124,(DTAUX(I,KK),KK = 1,KM)
881     !                       GO TO 92
882     !                    ENDIF
883     !                 ENDIF
884     !91            CONTINUE
885     !92         CONTINUE
886     !CMIC$ END GUARD 37
887     !123        FORMAT('  *** MIGWD PRINT *** I=',I3,' LAT=',I3,' KDT=',I3)
888     !124        FORMAT(2X,  10E13.6)
889     !        ENDIF
890     !     ENDIF
891     !
892     !      print *,' in gwdps_lm.f 18  =',A(ipt(1),idxzb(1))
893     !    &,                          B(ipt(1),idxzb(1)),me
894           RETURN
895           END
896