Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_GWD.F @ 2367

Last change on this file since 2367 was 1, checked in by lfita, 11 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 63.7 KB

Line
1	!
2	!-- Module for Gravity Wave Drag (GWD) and Mountain Blocking (MB)
3	!
4	!-- Initially incorporated into the WRF NMM from the GFS by B. Ferrier
5	! in April/May 2007.
6	!
7	! Search for "ORIGINAL DOCUMENTATION BLOCK" for further description.
8	!
9	!#######################################################################
10	!
11	MODULE module_gwd
12	!
13	! USE MODULE_DM ! to get processor element
14	USE MODULE_EXT_INTERNAL ! to assign fortan unit number
15	!
16	!-- Contains subroutines GWD_init, GWD_driver, and GWD_col
17	!
18	!#######################################################################
19	!
20	INTEGER, PARAMETER :: KIND_PHYS=SELECTED_REAL_KIND(13,60) ! the '60' maps to 64-bit real
21	INTEGER,PRIVATE,SAVE :: IMX, NMTVR, IDBG, JDBG
22	REAL,PRIVATE,SAVE :: RAD_TO_DEG !-- Convert radians to degrees
23	REAL,PRIVATE,SAVE :: DEG_TO_RAD !-- Convert degrees to radians
24	REAL (KIND=KIND_PHYS),PRIVATE,SAVE :: DELTIM,RDELTIM
25	REAL(kind=kind_phys),PRIVATE,PARAMETER :: SIGFAC=0.0 !-- Key tunable parameter
26	!dbg real,private,save :: dumin,dumax,dvmin,dvmax !dbg
27	!
28	CONTAINS
29	!
30	!-- Initialize variables used in GWD + MB
31	!
32	SUBROUTINE GWD_init (DTPHS,DELX,DELY,CEN_LAT,CEN_LON,RESTRT &
33	& ,GLAT,GLON,CROT,SROT,HANGL &
34	& ,IDS,IDE,JDS,JDE,KDS,KDE &
35	& ,IMS,IME,JMS,JME,KMS,KME &
36	& ,ITS,ITE,JTS,JTE,KTS,KTE )
37
38	!
39	IMPLICIT NONE
40	!
41	!== INPUT:
42	!-- DELX, DELY - DX, DY grid resolutions in zonal, meridional directions (m)
43	!-- CEN_LAT, CEN_LON - central latitude, longitude (degrees)
44	!-- RESTRT - logical flag for restart file (true) or WRF input file (false)
45	!-- GLAT, GLON - central latitude, longitude at mass points (radians)
46	!-- CROT, SROT - cosine and sine of the angle between Earth and model coordinates
47	!-- HANGL - angle of the mountain range w/r/t east (convert to degrees)
48	!
49	!-- Saved variables within module:
50	!-- IMX - in the GFS it is an equivalent number of points along a latitude
51	! circle (e.g., IMX=3600 for a model resolution of 0.1 deg)
52	! => Calculated at start of model integration in GWD_init
53	!-- NMTVR - number of input 2D orographic fields
54	!-- GRAV = gravitational acceleration
55	!-- DELTIM - physics time step (s)
56	!-- RDELTIM - reciprocal of physics time step (s)
57	!
58	!== INPUT indices:
59	!-- ids start index for i in domain
60	!-- ide end index for i in domain
61	!-- jds start index for j in domain
62	!-- jde end index for j in domain
63	!-- kds start index for k in domain
64	!-- kde end index for k in domain
65	!-- ims start index for i in memory
66	!-- ime end index for i in memory
67	!-- jms start index for j in memory
68	!-- jme end index for j in memory
69	!-- kms start index for k in memory
70	!-- kme end index for k in memory
71	!-- its start index for i in tile
72	!-- ite end index for i in tile
73	!-- jts start index for j in tile
74	!-- jte end index for j in tile
75	!-- kts start index for k in tile
76	!-- kte end index for k in tile
77	!
78	REAL, INTENT(IN) :: DTPHS,DELX,DELY,CEN_LAT,CEN_LON
79	LOGICAL, INTENT(IN) :: RESTRT
80	REAL, INTENT(IN), DIMENSION (ims:ime,jms:jme) :: GLON,GLAT
81	REAL, INTENT(OUT), DIMENSION (ims:ime,jms:jme) :: CROT,SROT
82	REAL, INTENT(INOUT), DIMENSION (ims:ime,jms:jme) :: HANGL
83	INTEGER, INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
84	& ,IMS,IME,JMS,JME,KMS,KME &
85	& ,ITS,ITE,JTS,JTE,KTS,KTE
86	!
87	!-- Local variables:
88	!
89	REAL, PARAMETER :: POS1=1.,NEG1=-1.
90	REAL :: DX,DTR,LAT0,LoN0,CLAT0,SLAT0,CLAT,DLON,X,Y,TLON,ROT
91	INTEGER :: I,J
92
93	!dbg
94	!dbg real :: xdbg,ydbg,d_x,d_y,dist,dist_min
95	!dbg data xdbg,ydbg / -118.3,36 / ! 118.3W 36 N
96
97	!
98	!-----------------------------------------------------------------------
99	!
100	DX=SQRT((DELX)2+(DELY)2) !-- Model resolution in degrees
101	!-- IMX is the number of grid points along a latitude circle in the GFS
102	IMX=INT(360./DX+.5)
103
104	!dbg IMX=1152 !dbg -- Match the grid point printed from GFS run
105
106	NMTVR=14 !-- 14 input fields for orography
107	DELTIM=DTPHS
108	RDELTIM=1./DTPHS
109	!
110	!-- Calculate angle of rotation (ROT) between Earth and model coordinates,
111	! but pass back out cosine (CROT) and sine (SROT) of this angle
112	!
113	DTR=ACOS(-1.)/180. !-- convert from degrees to radians
114	DEG_TO_RAD=DTR !-- save conversion from degrees to radians
115	!
116	LAT0=DTR*CEN_LON !-- central latitude of grid in radians
117	LoN0=DTR*CEN_LAT !-- central longitude of grid in radians
118	!
119	DTR=1./DTR !-- convert from radians to degrees
120	RAD_TO_DEG=DTR !-- save conversion from radians to degrees
121	!
122	CLAT0=COS(LAT0)
123	SLAT0=SIN(LAT0)
124	DO J=JTS,JTE
125	DO I=ITS,ITE
126	CLAT=COS(GLAT(I,J))
127	DLON=GLON(I,J)-LoN0
128	X=CLAT0CLATCOS(DLON)+SLAT0*SIN(GLAT(I,J))
129	Y=-CLAT*SIN(DLON)
130	TLON=ATAN(Y/X) !-- model longitude
131	X=SLAT0*SIN(TLON)/CLAT
132	Y=MIN(POS1, MAX(NEG1, X) )
133	ROT=ASIN(Y) !-- angle between geodetic & model coordinates
134	CROT(I,J)=COS(ROT)
135	SROT(I,J)=SIN(ROT)
136	ENDDO !-- I
137	ENDDO !-- J
138	IF (.NOT.RESTRT) THEN
139	!-- Convert from radians to degrees for WRF input files only
140	DO J=JTS,JTE
141	DO I=ITS,ITE
142	HANGL(I,J)=DTR*HANGL(I,J) !-- convert to degrees (+/-90 deg)
143	ENDDO !-- I
144	ENDDO !-- J
145	ENDIF
146	!dbg
147	!dbg dumin=-1.
148	!dbg dumax=1.
149	!dbg dvmin=-1.
150	!dbg dvmax=1.
151	!dbg print *,'delx=',delx,' dely=',dely,' dx=',dx,' imx=',imx
152	!dbg dtr=1./dtr !-- convert from degrees back to radians
153	!dbg dist_min=dtr*DX !-- grid length in radians
154	!dbg xdbg=dtr*xdbg !-- convert xdbg to radians
155	!dbg ydbg=dtr*ydbg !-- convert ydbg to radians
156	!dbg idbg=-100
157	!dbg jdbg=-100
158	!dbg print *,'dtr,dx,dist_min,xdbg,ydbg=',dtr,dx,dist_min,xdbg,ydbg
159	!dbg do j=jts,jte
160	!dbg do i=its,ite
161	!dbg !-- Find i,j for xdbg, ydbg
162	!dbg d_x=cos(glat(i,j))*(glon(i,j)-xdbg)
163	!dbg d_y=(glat(i,j)-ydbg)
164	!dbg dist=sqrt(d_xd_x+d_yd_y)
165	!dbg !! print *,'i,j,glon,glat,d_x,d_y,dist=',i,j,glon(i,j),glat(i,j),d_x,d_y,dist
166	!dbg if (dist < dist_min) then
167	!dbg dist_min=dist
168	!dbg idbg=i
169	!dbg jdbg=j
170	!dbg print *,'dist_min,idbg,jdbg=',dist_min,idbg,jdbg
171	!dbg endif
172	!dbg enddo !-- I
173	!dbg enddo !-- J
174	!dbg if (idbg>0 .and. jdbg>0) print *,'idbg=',idbg,' jdbg=',jdbg
175
176	!
177	END SUBROUTINE GWD_init
178	!
179	!-----------------------------------------------------------------------
180	!
181	SUBROUTINE GWD_driver(U,V,T,Q,Z,DP,PINT,PMID,EXNR, KPBL, ITIME &
182	& ,HSTDV,HCNVX,HASYW,HASYS,HASYSW,HASYNW &
183	& ,HLENW,HLENS,HLENSW,HLENNW &
184	& ,HANGL,HANIS,HSLOP,HZMAX,CROT,SROT &
185	& ,DUDT,DVDT,UGWDsfc,VGWDsfc &
186	& ,IDS,IDE,JDS,JDE,KDS,KDE &
187	& ,IMS,IME,JMS,JME,KMS,KME &
188	& ,ITS,ITE,JTS,JTE,KTS,KTE )
189	!
190	!== INPUT:
191	!-- U, V - zonal (U), meridional (V) winds at mass points (m/s)
192	!-- T, Q - temperature (C), specific humidity (kg/kg)
193	!-- DP - pressure thickness (Pa)
194	!-- Z - geopotential height (m)
195	!-- PINT, PMID - interface and midlayer pressures, respectively (Pa)
196	!-- EXNR - (p/p0)**(Rd/Cp)
197	!-- KPBL - vertical index at PBL top
198	!-- ITIME - model time step (=NTSD)
199	!-- HSTDV - orographic standard deviation
200	!-- HCNVX - normalized 4th moment of the orographic convexity
201	!-- Template for the next two sets of 4 arrays:
202	! NWD 1 2 3 4 5 6 7 8
203	! WD W S SW NW E N NE SE
204	!-- Orographic asymmetry (HASYx, x=1-4) for upstream & downstream flow (4 planes)
205	!-- * HASYW - orographic asymmetry for upstream & downstream flow in W-E plane
206	!-- * HASYS - orographic asymmetry for upstream & downstream flow in S-N plane
207	!-- * HASYSW - orographic asymmetry for upstream & downstream flow in SW-NE plane
208	!-- * HASYNW - orographic asymmetry for upstream & downstream flow in NW-SE plane
209	!-- Orographic length scale or mountain width (4 planes)
210	!-- * HLENW - orographic length scale for upstream & downstream flow in W-E plane
211	!-- * HLENS - orographic length scale for upstream & downstream flow in S-N plane
212	!-- * HLENSW - orographic length scale for upstream & downstream flow in SW-NE plane
213	!-- * HLENNW - orographic length scale for upstream & downstream flow in NW-SE plane
214	!-- HANGL - angle (degrees) of the mountain range w/r/t east
215	!-- HANIS - anisotropy/aspect ratio of orography
216	!-- HSLOP - slope of orography
217	!-- HZMAX - max height above mean orography
218	!-- CROT, SROT - cosine & sine of the angle between Earth & model coordinates
219	!
220	!== OUTPUT:
221	!-- DUDT, DVDT - zonal, meridional wind tendencies
222	!-- UGWDsfc, VGWDsfc - zonal, meridional surface wind stresses (N/m**2)
223	!
224	!== INPUT indices:
225	!-- ids start index for i in domain
226	!-- ide end index for i in domain
227	!-- jds start index for j in domain
228	!-- jde end index for j in domain
229	!-- kds start index for k in domain
230	!-- kde end index for k in domain
231	!-- ims start index for i in memory
232	!-- ime end index for i in memory
233	!-- jms start index for j in memory
234	!-- jme end index for j in memory
235	!-- kms start index for k in memory
236	!-- kme end index for k in memory
237	!-- its start i index for tile
238	!-- ite end i index for tile
239	!-- jts start j index for tile
240	!-- jte end j index for tile
241	!-- kts start index for k in tile
242	!-- kte end index for k in tile
243	!
244	!-- INPUT variables:
245	!
246	REAL, INTENT(IN), DIMENSION (ims:ime, kms:kme, jms:jme) :: &
247	& U,V,T,Q,Z,DP,PINT,PMID,EXNR
248	REAL, INTENT(IN), DIMENSION (ims:ime, jms:jme) :: HSTDV,HCNVX &
249	& ,HASYW,HASYS,HASYSW,HASYNW,HLENW,HLENS,HLENSW,HLENNW,HANGL &
250	& ,HANIS,HSLOP,HZMAX,CROT,SROT
251	INTEGER, INTENT(IN), DIMENSION (ims:ime, jms:jme) :: KPBL
252	INTEGER, INTENT(IN) :: ids,ide,jds,jde,kds,kde &
253	&, ims,ime,jms,jme,kms,kme &
254	&, its,ite,jts,jte,kts,kte,ITIME
255
256	!
257	!-- OUTPUT variables:
258	!
259	REAL, INTENT(OUT), DIMENSION (ims:ime, kms:kme, jms:jme) :: &
260	& DUDT,DVDT
261	REAL, INTENT(OUT), DIMENSION (ims:ime, jms:jme) :: UGWDsfc,VGWDsfc
262	!
263	!-- Local variables
264	!-- DUsfc, DVsfc - zonal, meridional wind stresses (diagnostics)
265	!
266	INTEGER, PARAMETER :: IM=1 !-- Reduces changes in subroutine GWPDS
267	REAL(KIND=KIND_PHYS), PARAMETER :: G=9.806, GHALF=.5*G &
268	&, THRESH=1.E-6, dtlarge=1. !dbg
269	INTEGER, DIMENSION (IM) :: LPBL
270	REAL(KIND=KIND_PHYS), DIMENSION (IM,4) :: OA4,CLX4
271	REAL(KIND=KIND_PHYS), DIMENSION (IM) :: DUsfc,DVsfc &
272	&, HPRIME,OC,THETA,SIGMA,GAMMA,ELVMAX
273	REAL(KIND=KIND_PHYS), DIMENSION (IM,KTS:KTE) :: DUDTcol,DVDTcol &
274	&, Ucol,Vcol,Tcol,Qcol,DPcol,Pcol,EXNcol,PHIcol
275	REAL(KIND=KIND_PHYS), DIMENSION (IM,KTS:KTE+1) :: PINTcol,PHILIcol
276	INTEGER :: I,J,IJ,K,Imid,Jmid
277	REAL :: Ugeo,Vgeo,Umod,Vmod, TERRtest,TERRmin
278	REAL(KIND=KIND_PHYS) :: TEST
279	CHARACTER(LEN=255) :: message
280
281	!dbg
282	logical :: lprnt !dbg
283	character (len=26) :: label
284	integer :: kpblmin,kpblmax, mype, iunit
285	real :: hzmaxmin,hzmaxmax,hanglmin,hanglmax,hslopmin,hslopmax,hanismin,hanismax &
286	,hstdvmin,hstdvmax,hcnvxmin,hcnvxmax,hasywmin,hasywmax,hasysmin,hasysmax &
287	,hasyswmin,hasyswmax,hasynwmin,hasynwmax,hlenwmin,hlenwmax,hlensmin,hlensmax &
288	,hlenswmin,hlenswmax,hlennwmin,hlennwmax,zsmin,zsmax,delu,delv
289	! Added this declaration
290	real :: helvmin,helvmax
291	!dbg
292
293	!
294	!-------------------------- Executable below -------------------------
295	!
296
297	lprnt=.false.
298	!dbg
299	if (itime <= 1) then
300	CALL WRF_GET_MYPROC(MYPE) !-- Get processor number
301	kpblmin=100
302	kpblmax=-100
303	hzmaxmin=1.e6
304	hzmaxmax=-1.e6
305	hanglmin=1.e6
306	hanglmax=-1.e6
307	hslopmin=1.e6
308	hslopmax=-1.e6
309	hanismin=1.e6
310	hanismax=-1.e6
311	hstdvmin=1.e6
312	hstdvmax=-1.e6
313	hcnvxmin=1.e6
314	hcnvxmax=-1.e6
315	hasywmin=1.e6
316	hasywmax=-1.e6
317	hasysmin=1.e6
318	hasysmax=-1.e6
319	hasyswmin=1.e6
320	hasyswmax=-1.e6
321	hasynwmin=1.e6
322	hasynwmax=-1.e6
323	hlenwmin=1.e6
324	hlenwmax=-1.e6
325	hlensmin=1.e6
326	hlensmax=-1.e6
327	hlenswmin=1.e6
328	hlenswmax=-1.e6
329	hlennwmin=1.e6
330	hlennwmax=-1.e6
331	zsmin=1.e6
332	zsmax=-1.e6
333	! Added initialization of helvmin and helvmax
334	helvmin=1.e6
335	helvmax=-1.e6
336	do j=jts,jte
337	do i=its,ite
338	kpblmin=min(kpblmin,kpbl(i,j))
339	kpblmax=max(kpblmax,kpbl(i,j))
340	helvmin=min(helvmin,hzmax(i,j))
341	helvmax=max(helvmax,hzmax(i,j))
342	hanglmin=min(hanglmin,hangl(i,j))
343	hanglmax=max(hanglmax,hangl(i,j))
344	hslopmin=min(hslopmin,hslop(i,j))
345	hslopmax=max(hslopmax,hslop(i,j))
346	hanismin=min(hanismin,hanis(i,j))
347	hanismax=max(hanismax,hanis(i,j))
348	hstdvmin=min(hstdvmin,hstdv(i,j))
349	hstdvmax=max(hstdvmax,hstdv(i,j))
350	hcnvxmin=min(hcnvxmin,hcnvx(i,j))
351	hcnvxmax=max(hcnvxmax,hcnvx(i,j))
352	hasywmin=min(hasywmin,hasyw(i,j))
353	hasywmax=max(hasywmax,hasyw(i,j))
354	hasysmin=min(hasysmin,hasys(i,j))
355	hasysmax=max(hasysmax,hasys(i,j))
356	hasyswmin=min(hasyswmin,hasysw(i,j))
357	hasyswmax=max(hasyswmax,hasysw(i,j))
358	hasynwmin=min(hasynwmin,hasynw(i,j))
359	hasynwmax=max(hasynwmax,hasynw(i,j))
360	hlenwmin=min(hlenwmin,hlenw(i,j))
361	hlenwmax=max(hlenwmax,hlenw(i,j))
362	hlensmin=min(hlensmin,hlens(i,j))
363	hlensmax=max(hlensmax,hlens(i,j))
364	hlenswmin=min(hlenswmin,hlensw(i,j))
365	hlenswmax=max(hlenswmax,hlensw(i,j))
366	hlennwmin=min(hlennwmin,hlennw(i,j))
367	hlennwmax=max(hlennwmax,hlennw(i,j))
368	zsmin=min(zsmin,z(i,1,j))
369	zsmax=max(zsmax,z(i,1,j))
370	enddo
371	enddo
372	write(message,*) 'Maximum and minimum values within GWD-driver for MYPE=',MYPE
373	write(message,"(i3,2(a,e12.5))") mype,': deltim=',deltim,' rdeltim=',rdeltim
374	CALL wrf_message(trim(message))
375	write(message,"(i3,2(a,i3))") mype,': kpblmin=',kpblmin,' kpblmax=',kpblmax
376	CALL wrf_message(trim(message))
377	write(message,"(i3,2(a,e12.5))") mype,': helvmin=',helvmin,' helvmax=',helvmax
378	CALL wrf_message(trim(message))
379	write(message,"(i3,2(a,e12.5))") mype,': hanglmin=',hanglmin,' hanglmax=',hanglmax
380	CALL wrf_message(trim(message))
381	write(message,"(i3,2(a,e12.5))") mype,': hslopmin=',hslopmin,' hslopmax=',hslopmax
382	CALL wrf_message(trim(message))
383	write(message,"(i3,2(a,e12.5))") mype,': hanismin=',hanismin,' hanismax=',hanismax
384	CALL wrf_message(trim(message))
385	write(message,"(i3,2(a,e12.5))") mype,': hstdvmin=',hstdvmin,' hstdvmax=',hstdvmax
386	CALL wrf_message(trim(message))
387	write(message,"(i3,2(a,e12.5))") mype,': hcnvxmin=',hcnvxmin,' hcnvxmax=',hcnvxmax
388	CALL wrf_message(trim(message))
389	write(message,"(i3,2(a,e12.5))") mype,': hasywmin=',hasywmin,' hasywmax=',hasywmax
390	CALL wrf_message(trim(message))
391	write(message,"(i3,2(a,e12.5))") mype,': hasysmin=',hasysmin,' hasysmax=',hasysmax
392	CALL wrf_message(trim(message))
393	write(message,"(i3,2(a,e12.5))") mype,': hasyswmin=',hasyswmin,' hasyswmax=',hasyswmax
394	CALL wrf_message(trim(message))
395	write(message,"(i3,2(a,e12.5))") mype,': hasynwmin=',hasynwmin,' hasynwmax=',hasynwmax
396	CALL wrf_message(trim(message))
397	write(message,"(i3,2(a,e12.5))") mype,': hlenwmin=',hlenwmin,' hlenwmax=',hlenwmax
398	CALL wrf_message(trim(message))
399	write(message,"(i3,2(a,e12.5))") mype,': hlensmin=',hlensmin,' hlensmax=',hlensmax
400	CALL wrf_message(trim(message))
401	write(message,"(i3,2(a,e12.5))") mype,': hlenswmin=',hlenswmin,' hlenswmax=',hlenswmax
402	CALL wrf_message(trim(message))
403	write(message,"(i3,2(a,e12.5))") mype,': hlennwmin=',hlennwmin,' hlennwmax=',hlennwmax
404	CALL wrf_message(trim(message))
405	write(message,"(i3,2(a,e12.5))") mype,': zsmin=',zsmin,' zsmax=',zsmax
406	CALL wrf_message(trim(message))
407
408	endif ! if (itime <= 1) then
409	!dbg
410
411	!
412	!-- Initialize variables
413	!
414	DO J=JMS,JME
415	DO K=KMS,KME
416	DO I=IMS,IME
417	DUDT(I,K,J)=0.
418	DVDT(I,K,J)=0.
419	ENDDO
420	ENDDO
421	ENDDO
422	!
423	DO J=JMS,JME
424	DO I=IMS,IME
425	UGWDsfc(I,J)=0.
426	VGWDsfc(I,J)=0.
427	ENDDO
428	ENDDO
429	!
430	!-- For debugging, find approximate center point within each tile
431	!
432	!dbg Imid=.5*(ITS+ITE)
433	!dbg Jmid=.5*(JTS+JTE)
434	!
435	DO J=JTS,JTE
436	DO I=ITS,ITE
437	if (kpbl(i,j)<kts .or. kpbl(i,j)>kte) go to 100
438	!
439	!-- Initial test to see if GWD calculations should be made, otherwise skip
440	!
441	TERRtest=HZMAX(I,J)+SIGFAC*HSTDV(I,J)
442	TERRmin=Z(I,2,J)-Z(I,1,J)
443	IF (TERRtest < TERRmin) GO TO 100
444	!
445	!-- For debugging:
446	!
447	!dbg lprnt=.false.
448	!dbg if (i==idbg .and. j==jdbg .and. itime<=1) lprnt=.true.
449	!dbg ! 200 CONTINUE
450	!
451	DO K=KTS,KTE
452	DUDTcol(IM,K)=0.
453	DVDTcol(IM,K)=0.
454	!
455	!-- Transform/rotate winds from model to geodetic (Earth) coordinates
456	!
457	Ucol(IM,K)=U(I,K,J)CROT(I,J)+V(I,K,J)SROT(I,J)
458	Vcol(IM,K)=V(I,K,J)CROT(I,J)-U(I,K,J)SROT(I,J)
459	!
460	Tcol(IM,K)=T(I,K,J)
461	Qcol(IM,K)=Q(I,K,J)
462	!
463	!-- Convert from Pa to centibars, which is what's used in subroutine GWD_col
464	!
465	DPcol(IM,K)=.001*DP(I,K,J)
466	PINTcol(IM,K)=.001*PINT(I,K,J)
467	Pcol(IM,K)=.001*PMID(I,K,J)
468	EXNcol(IM,K)=EXNR(I,K,J)
469	!
470	!-- Next 2 fields are geopotential above the surface at the lower interface
471	! and at midlayer
472	!
473	PHILIcol(IM,K)=G*(Z(I,K,J)-Z(I,1,J))
474	PHIcol(IM,K)=GHALF(Z(I,K,J)+Z(I,K+1,J))-GZ(I,1,J)
475	ENDDO !- K
476	!
477	PINTcol(IM,KTE+1)=.001*PINT(I,KTE+1,J)
478	PHILIcol(IM,KTE+1)=G*(Z(I,KTE+1,J)-Z(I,1,J))
479	!
480	!-- Terrain-specific inputs:
481	!
482	HPRIME(IM)=HSTDV(I,J) !-- standard deviation of orography
483	OC(IM)=HCNVX(I,J) !-- Normalized convexity
484	OA4(IM,1)=HASYW(I,J) !-- orographic asymmetry in W-E plane
485	OA4(IM,2)=HASYS(I,J) !-- orographic asymmetry in S-N plane
486	OA4(IM,3)=HASYSW(I,J) !-- orographic asymmetry in SW-NE plane
487	OA4(IM,4)=HASYNW(I,J) !-- orographic asymmetry in NW-SE plane
488	CLX4(IM,1)=HLENW(I,J) !-- orographic length scale in W-E plane
489	CLX4(IM,2)=HLENS(I,J) !-- orographic length scale in S-N plane
490	CLX4(IM,3)=HLENSW(I,J) !-- orographic length scale in SW-NE plane
491	CLX4(IM,4)=HLENNW(I,J) !-- orographic length scale in NW-SE plane
492	THETA(IM)=HANGL(I,J) !
493	SIGMA(IM)=HSLOP(I,J) !
494	GAMMA(IM)=HANIS(I,J) !
495	ELVMAX(IM)=HZMAX(I,J) !
496	LPBL(IM)=KPBL(I,J) !
497	!dbg IF (LPBL(IM)<KTS .OR. LPBL(IM)>KTE) &
498	!dbg & print *,'Wacky values for KPBL: I,J,N,LPBL=',I,J,ITIME,LPBL(IM)
499	!
500	!-- Output (diagnostics)
501	!
502	DUsfc(IM)=0. !-- U wind stress
503	DVsfc(IM)=0. !-- V wind stress
504	!
505	!dbg
506	!dbg if (LPRNT) then
507	!dbg !
508	!dbg !-- Following code is for ingesting GFS-derived inputs for final testing
509	!dbg !
510	!dbg CALL INT_GET_FRESH_HANDLE(iunit)
511	!dbg close(iunit)
512	!dbg open(unit=iunit,file='gfs_gwd.input',form='formatted',iostat=ier)
513	!dbg read(iunit,*) ! skip line
514	!dbg read(iunit,*) (Ucol(im,k), k=kts,kte)
515	!dbg read(iunit,*) ! skip line
516	!dbg read(iunit,*) (Vcol(im,k), k=kts,kte)
517	!dbg read(iunit,*) ! skip line
518	!dbg read(iunit,*) (Tcol(im,k), k=kts,kte)
519	!dbg read(iunit,*) ! skip line
520	!dbg read(iunit,*) (Qcol(im,k), k=kts,kte)
521	!dbg read(iunit,*) ! skip line
522	!dbg read(iunit,*) (PINTcol(im,k), k=kts,kte+1)
523	!dbg read(iunit,*) ! skip line
524	!dbg read(iunit,*) (DPcol(im,k), k=kts,kte)
525	!dbg read(iunit,*) ! skip line
526	!dbg read(iunit,*) (Pcol(im,k), k=kts,kte)
527	!dbg read(iunit,*) ! skip line
528	!dbg read(iunit,*) (EXNcol(im,k), k=kts,kte)
529	!dbg read(iunit,*) ! skip line
530	!dbg read(iunit,*) (PHILIcol(im,k), k=kts,kte+1)
531	!dbg read(iunit,*) ! skip line
532	!dbg read(iunit,*) (PHIcol(im,k), k=kts,kte)
533	!dbg read(iunit,*) ! skip line
534	!dbg read(iunit,*) hprime(im)
535	!dbg read(iunit,*) ! skip line
536	!dbg read(iunit,*) oc(im)
537	!dbg read(iunit,*) ! skip line
538	!dbg read(iunit,*) (oa4(im,k), k=1,4)
539	!dbg read(iunit,*) ! skip line
540	!dbg read(iunit,*) (clx4(im,k), k=1,4)
541	!dbg read(iunit,*) ! skip line
542	!dbg read(iunit,*) theta(im)
543	!dbg read(iunit,*) ! skip line
544	!dbg read(iunit,*) sigma(im)
545	!dbg read(iunit,*) ! skip line
546	!dbg read(iunit,*) gamma(im)
547	!dbg read(iunit,*) ! skip line
548	!dbg read(iunit,*) elvmax(im)
549	!dbg read(iunit,*) ! skip line
550	!dbg read(iunit,*) lpbl(im)
551	!dbg close(iunit)
552	write(label,"('GWD:i,j,n=',2i5,i6)") I,J,ITIME
553	!dbg write(6,"(2a)") LABEL,' in GWD_driver: K U V T Q Pi DP P EX PHIi PHI'
554	!dbg do k=kts,kte
555	!dbg write(6,"(i3,10e12.4)") k,Ucol(im,k),Vcol(im,k),Tcol(im,k) &
556	!dbg ,Qcol(im,k),PINTcol(im,k),DPcol(im,k),Pcol(im,k),EXNcol(im,k) &
557	!dbg ,PHILIcol(im,k),PHIcol(im,k)
558	!dbg enddo
559	!dbg write(6,"(2(a,e12.4),2(a,4e12.4/),4(a,e12.4),a,i3) )") &
560	!dbg 'GWD_driver: hprime(im)=',hprime(im),' oc(im)=',oc(im) &
561	!dbg ,' oa4(im,1-4)=',(oa4(im,k),k=1,4) &
562	!dbg ,' clx4(im,1-4)=',(clx4(im,k),k=1,4) &
563	!dbg ,' theta=',theta(im),' sigma(im)=',sigma(im) &
564	!dbg ,' gamma(im)=',gamma(im),' elvmax(im)=',elvmax(im) &
565	!dbg ,' lpbl(im)=',lpbl(im)
566	!dbg endif
567	!dbg
568	!=======================================================================
569	!
570	CALL GWD_col(DVDTcol,DUDTcol, DUsfc,DVsfc & ! Output
571	&, Ucol,Vcol,Tcol,Qcol,PINTcol,DPcol,Pcol,EXNcol & ! Met input
572	&, PHILIcol,PHIcol & ! Met input
573	&, HPRIME,OC,OA4,CLX4,THETA,SIGMA,GAMMA,ELVMAX & ! Topo input
574	&, LPBL,IM,KTS,KTE,LABEL,LPRNT ) ! Indices + debugging
575	!
576	!=======================================================================
577	!
578	!dbg
579	!dbg !
580	!dbg ! IF (.NOT.LPRNT) THEN
581	!dbg ! TEST=0.
582	!dbg ! DO K=KTS,KTE
583	!dbg ! TEST=MAX( TEST, ABS(DUDTcol(IM,K)), ABS(DVDTcol(IM,K)) )
584	!dbg ! if ( ABS(DUDTcol(IM,K)) > RDELTIM) print *,'k,DUDTcol=',k,DUDTcol(IM,K)
585	!dbg ! if ( ABS(DVDTcol(IM,K)) > RDELTIM) print *,'k,DVDTcol=',k,DVDTcol(IM,K)
586	!dbg ! ENDDO
587	!dbg ! IF (TEST>RDELTIM) THEN
588	!dbg ! LPRNT=.TRUE.
589	!dbg ! Imid=I
590	!dbg ! Jmid=J
591	!dbg ! GO TO 200
592	!dbg ! ENDIF
593	!dbg ! ENDIF
594	!dbg ! TEST=ABS(DUsfc(IM))+ABS(DVsfc(IM))
595	!dbg ! if (.not.lprnt) then
596	!dbg ! do k=kts,kte
597	!dbg ! du=DUDTcol(IM,K)*DELTIM
598	!dbg ! dv=DVDTcol(IM,K)*DELTIM
599	!dbg ! if (du < dumin) then
600	!dbg ! lprnt=.true.
601	!dbg ! dumin=1.5*du
602	!dbg ! endif
603	!dbg ! if (du > dumax) then
604	!dbg ! lprnt=.true.
605	!dbg ! dumax=1.5*du
606	!dbg ! endif
607	!dbg ! if (dv < dvmin) then
608	!dbg ! lprnt=.true.
609	!dbg ! dvmin=1.5*dv
610	!dbg ! endif
611	!dbg ! if (dv > dvmax) then
612	!dbg ! lprnt=.true.
613	!dbg ! dvmax=1.5*dv
614	!dbg ! endif
615	!dbg ! enddo
616	!dbg ! if (lprnt) go to 200
617	!dbg ! else
618	!dbg if (lprnt) then
619	!dbg print *,'DUsfc,DVsfc,CROT,SROT,DELTIM=',DUsfc(IM),DVsfc(IM) &
620	!dbg &,CROT(I,J),SROT(I,J),DELTIM
621	!dbg print ,' K P \| Ucol Ugeo DUDTcolDT U Umod DU=DUDT*DT ' &
622	!dbg ,'\| Vcol Vgeo DVDTcolDT V Vmod DV=DVDTDT'
623	!dbg ENDIF
624
625	DO K=KTS,KTE
626	TEST=ABS(DUDTcol(IM,K))+ABS(DVDTcol(IM,K))
627	IF (TEST > THRESH) THEN
628
629	!dbg
630	!dbg if (lprnt) then
631	!dbg !dbg DUDTcol(IM,K)=0. !-- Test rotation
632	!dbg !dbg DVDTcol(IM,K)=0. !-- Test rotation
633	!dbg !-- Now replace with the original winds before they were written over by
634	!dbg ! the values from the GFS
635	!dbg Ucol(IM,K)=U(I,K,J)CROT(I,J)+V(I,K,J)SROT(I,J)
636	!dbg Vcol(IM,K)=V(I,K,J)CROT(I,J)-U(I,K,J)SROT(I,J)
637	!dbg endif
638
639	!
640	!-- First update winds in geodetic coordinates
641	!
642	Ugeo=Ucol(IM,K)+DUDTcol(IM,K)*DELTIM
643	Vgeo=Vcol(IM,K)+DVDTcol(IM,K)*DELTIM
644	!
645	!-- Transform/rotate winds from geodetic back to model coordinates
646	!
647	Umod=UgeoCROT(I,J)-VgeoSROT(I,J)
648	Vmod=UgeoSROT(I,J)+VgeoCROT(I,J)
649	!
650	!-- Calculate wind tendencies from the updated model winds
651	!
652	DUDT(I,K,J)=RDELTIM*(Umod-U(I,K,J))
653	DVDT(I,K,J)=RDELTIM*(Vmod-V(I,K,J))
654	!
655	!dbg
656
657	test=abs(dudt(i,k,j))+abs(dvdt(i,k,j))
658	if (test > dtlarge) write(6,"(2a,i2,2(a,e12.4))") &
659	label,' => k=',k,' dudt=',dudt(i,k,j),' dvdt=',dvdt(i,k,j)
660
661	!dbg if (lprnt) write(6,"(i2,f8.2,2(' \| ',6f10.4)") K,10.*Pcol(IM,K) &
662	!dbg ,Ucol(IM,K),Ugeo,DUDTcol(IM,K)DELTIM,U(I,K,J),Umod,DUDT(I,K,J)DELTIM &
663	!dbg ,Vcol(IM,K),Vgeo,DVDTcol(IM,K)DELTIM,V(I,K,J),Vmod,DVDT(I,K,J)DELTIM
664	!dbg
665	ENDIF !- IF (TEST > THRESH) THEN
666	!
667	ENDDO !- K
668	!
669	!-- Transform/rotate surface wind stresses from geodetic to model coordinates
670	!
671	UGWDsfc(I,J)=DUsfc(IM)CROT(I,J)-DVsfc(IM)SROT(I,J)
672	VGWDsfc(I,J)=DUsfc(IM)SROT(I,J)+DVsfc(IM)CROT(I,J)
673	!
674	100 CONTINUE
675	ENDDO !- I
676	ENDDO !- J
677	!
678	END SUBROUTINE GWD_driver
679	!
680	!-----------------------------------------------------------------------
681	!
682	!-- "A", "B" (from GFS) in GWD_col are DVDTcol, DUDTcol, respectively in GWD_driver
683	!
684	SUBROUTINE GWD_col (A,B, DUsfc,DVsfc & !-- Output
685	&, U1,V1,T1,Q1, PRSI,DEL,PRSL,PRSLK, PHII,PHIL & !-- Met inputs
686	&, HPRIME,OC,OA4,CLX4,THETA,SIGMA,GAMMA,ELVMAX & !-- Topo inputs
687	&, KPBL,IM,KTS,KTE, LABEL,LPRNT ) !-- Input indices, debugging
688	!
689	!=== Output fields
690	!
691	!-- A (DUDT), B (DVDT) - output zonal & meridional wind tendencies in Earth coordinates (m s^-2)
692	!-- DUsfc, DVsfc - surface zonal meridional wind stresses in Earth coordinates (m s^-1?)
693	!
694	!=== Input fields
695	!
696	!-- U1, V1 - zonal, meridional wind (m/s)
697	!-- T1 - temperature (deg K)
698	!-- Q1 - specific humidity (kg/kg)
699	!-- PRSI - lower interface pressure in centibars (1000 Pa)
700	!-- DEL - pressure thickness of layer in centibars (1000 Pa)
701	!-- PRSL - midlayer pressure in centibars (1000 Pa)
702	!-- PRSLK - Exner function, (P/P0)**(Rd/CP)
703	!-- PHII - lower interface geopotential in mks units
704	!-- PHIL - midlayer geopotential in mks units
705	!-- KDT - number of time steps into integration for diagnostics
706	!-- HPRIME - orographic standard deviation
707	!-- OC - normalized 4th moment of the orographic convexity
708	!-- OA4 - orographic asymmetry for upstream & downstream flow measured
709	! along 4 vertical planes (W-E, S-N, SW-NE, NW-SE)
710	!-- CLX4 - orographic length scale or mountain width measured along
711	! 4 vertical planes (W-E, S-N, SW-NE, NW-SE)
712	!-- THETA - angle of the mountain range w/r/t east
713	!-- SIGMA - slope of orography
714	!-- GAMMA - anisotropy/aspect ratio
715	!-- ELVMAX - max height above mean orography
716	!-- KPBL(IM) - vertical index at the top of the PBL
717	!-- KM - number of vertical levels
718	!
719	!== For diagnostics
720	!-- LABEL - character string for diagnostic prints
721	!-- LPRNT - logical flag for prints
722	!
723	!#######################################################################
724	!################## ORIGINAL DOCUMENTATION BLOCK #####################
725	!###### The following comments are from the original GFS code ########
726	!#######################################################################
727	! ********************************************************************
728	! -----> I M P L E M E N T A T I O N V E R S I O N <----------
729	!
730	! --- Not in this code -- History of GWDP at NCEP----
731	! ---------------- -----------------------
732	! VERSION 3 MODIFIED FOR GRAVITY WAVES, LOCATION: .FR30(V3GWD) J
733	!--- 3.1 INCLUDES VARIABLE SATURATION FLUX PROFILE CF ISIGST
734	!--- 3.G INCLUDES PS COMBINED W/ PH (GLAS AND GFDL)
735	!----- ALSO INCLUDED IS RI SMOOTH OVER A THICK LOWER LAYER
736	!----- ALSO INCLUDED IS DECREASE IN DE-ACC AT TOP BY 1/2
737	!----- THE NMC GWD INCORPORATING BOTH GLAS(P&S) AND GFDL(MIGWD)
738	!----- MOUNTAIN INDUCED GRAVITY WAVE DRAG
739	!----- CODE FROM .FR30(V3MONNX) FOR MONIN3
740	!----- THIS VERSION (06 MAR 1987)
741	!----- THIS VERSION (26 APR 1987) 3.G
742	!----- THIS VERSION (01 MAY 1987) 3.9
743	!----- CHANGE TO FORTRAN 77 (FEB 1989) --- HANN-MING HENRY JUANG
744	!-----
745	!
746	! VERSION 4
747	! ----- This code -----
748	!
749	!----- MODIFIED TO IMPLEMENT THE ENHANCED LOW TROPOSPHERIC GRAVITY
750	!----- WAVE DRAG DEVELOPED BY KIM AND ARAKAWA(JAS, 1995).
751	! Orographic Std Dev (hprime), Convexity (OC), Asymmetry (OA4)
752	! and Lx (CLX4) are input topographic statistics needed.
753	!
754	!----- PROGRAMMED AND DEBUGGED BY HONG, ALPERT AND KIM --- JAN 1996.
755	!----- debugged again - moorthi and iredell --- may 1998.
756	!-----
757	! Further Cleanup, optimization and modification
758	! - S. Moorthi May 98, March 99.
759	!----- modified for usgs orography data (ncep office note 424)
760	! and with several bugs fixed - moorthi and hong --- july 1999.
761	!
762	!----- Modified & implemented into NRL NOGAPS
763	! - Young-Joon Kim, July 2000
764	!-----
765	! VERSION lm MB (6): oz fix 8/2003
766	! ----- This code -----
767	!
768	!------ Changed to include the Lott and Miller Mtn Blocking
769	! with some modifications by (j) 4/02
770	! From a Principal Coordinate calculation using the
771	! Hi Res 8 minute orography, the Angle of the
772	! mtn with that to the East (x) axis is THETA, the slope
773	! parameter SIGMA. The anisotropy is in GAMMA - all are input
774	! topographic statistics needed. These are calculated off-line
775	! as a function of model resolution in the fortran code ml01rg2.f,
776	! with script mlb2.sh. (j)
777	!----- gwdps_mb.f version (following lmi) elvmax < hncrit (j)
778	! MB3a expt to enhance elvmax mtn hgt see sigfac & hncrit
779	!-----
780	!----------------------------------------------------------------------C
781	!==== Below in "!GFS " are the original subroutine call and comments from
782	! /nwprod/sorc/global_fcst.fd/gwdps_v.f as of April 2007
783	!GFS SUBROUTINE GWDPS(IM,IX,IY,KM,A,B,U1,V1,T1,Q1,KPBL,
784	!GFS & PRSI,DEL,PRSL,PRSLK,PHII, PHIL,RCL,DELTIM,KDT,
785	!GFS & HPRIME,OC,OA4,CLX4,THETA,SIGMA,GAMMA,ELVMAX,
786	!GFS & DUsfc,DVsfc,G, CP, RD, RV, IMX,
787	!GFS & nmtvr, me, lprnt, ipr)
788	!GFS !------------------------------------------------------------------
789	!GFS ! USE
790	!GFS ! ROUTINE IS CALLED FROM GBPHYS (AFTER CALL TO MONNIN)
791	!GFS !
792	!GFS ! PURPOSE
793	!GFS ! USING THE GWD PARAMETERIZATIONS OF PS-GLAS AND PH-
794	!GFS ! GFDL TECHNIQUE. THE TIME TENDENCIES OF U V
795	!GFS ! ARE ALTERED TO INCLUDE THE EFFECT OF MOUNTAIN INDUCED
796	!GFS ! GRAVITY WAVE DRAG FROM SUB-GRID SCALE OROGRAPHY INCLUDING
797	!GFS ! CONVECTIVE BREAKING, SHEAR BREAKING AND THE PRESENCE OF
798	!GFS ! CRITICAL LEVELS
799	!GFS !
800	!GFS ! INPUT
801	!GFS ! A(IY,KM) NON-LIN TENDENCY FOR V WIND COMPONENT
802	!GFS ! B(IY,KM) NON-LIN TENDENCY FOR U WIND COMPONENT
803	!GFS ! U1(IX,KM) ZONAL WIND / SQRT(RCL) M/SEC AT T0-DT
804	!GFS ! V1(IX,KM) MERIDIONAL WIND / SQRT(RCL) M/SEC AT T0-DT
805	!GFS ! T1(IX,KM) TEMPERATURE DEG K AT T0-DT
806	!GFS ! Q1(IX,KM) SPECIFIC HUMIDITY AT T0-DT
807	!GFS !
808	!GFS ! RCL A scaling factor = RECIPROCAL OF SQUARE OF COS(LAT)
809	!GFS ! FOR MRF GFS.
810	!GFS ! DELTIM TIME STEP SECS
811	!GFS ! SI(N) P/PSFC AT BASE OF LAYER N
812	!GFS ! SL(N) P/PSFC AT MIDDLE OF LAYER N
813	!GFS ! DEL(N) POSITIVE INCREMENT OF P/PSFC ACROSS LAYER N
814	!GFS ! KPBL(IM) is the index of the top layer of the PBL
815	!GFS ! ipr & lprnt for diagnostics
816	!GFS !
817	!GFS ! OUTPUT
818	!GFS ! A, B AS AUGMENTED BY TENDENCY DUE TO GWDPS
819	!GFS ! OTHER INPUT VARIABLES UNMODIFIED.
820	!GFS ! ********************************************************************
821	!
822	IMPLICIT NONE
823	!
824	!-- INPUT:
825	!
826	INTEGER, INTENT(IN) :: IM,KTS,KTE
827	REAL(kind=kind_phys), INTENT(IN), DIMENSION(IM,KTS:KTE) :: &
828	& U1,V1,T1,Q1,DEL,PRSL,PRSLK,PHIL
829	REAL(kind=kind_phys), INTENT(IN), DIMENSION(IM,KTS:KTE+1) :: &
830	& PRSI,PHII
831	REAL(kind=kind_phys), INTENT(IN), DIMENSION(IM,4) :: OA4,CLX4
832	REAL(kind=kind_phys), INTENT(IN), DIMENSION(IM) :: &
833	& HPRIME,OC,THETA,SIGMA,GAMMA,ELVMAX
834	INTEGER, INTENT(IN), DIMENSION(IM) :: KPBL
835	CHARACTER (LEN=26), INTENT(IN) :: LABEL
836	LOGICAL, INTENT(IN) :: LPRNT
837	!
838	!-- OUTPUT:
839	!
840	REAL(kind=kind_phys), INTENT(INOUT), DIMENSION(IM,KTS:KTE) :: A,B
841	REAL(kind=kind_phys), INTENT(INOUT), DIMENSION(IM) :: DUsfc,DVsfc
842	!
843	!-----------------------------------------------------------------------
844	!-- LOCAL variables:
845	!-----------------------------------------------------------------------
846	!
847	! Some constants
848	!
849	!
850	REAL(kind=kind_phys), PARAMETER :: PI=3.1415926535897931 &
851	&, G=9.806, CP=1004.6, RD=287.04, RV=461.6 &
852	&, FV=RV/RD-1., RDI=1./RD, GOR=G/RD, GR2=G*GOR, GOCP=G/CP &
853	&, ROG=1./G, ROG2=ROG*ROG &
854	&, DW2MIN=1., RIMIN=-100., RIC=0.25, BNV2MIN=1.0E-5 &
855	&, EFMIN=0.0, EFMAX=10.0, hpmax=200.0 & ! or hpmax=2500.0
856	&, FRC=1.0, CE=0.8, CEOFRC=CE/FRC, frmax=100. &
857	&, CG=0.5, GMAX=1.0, CRITAC=5.0E-4, VELEPS=1.0 &
858	&, FACTOP=0.5, RLOLEV=500.0, HZERO=0., HONE=1. & ! or RLOLEV=0.5
859	&, HE_4=.0001, HE_2=.01 &
860	!
861	!-- Lott & Miller mountain blocking => aka "lm mtn blocking"
862	!
863	&, cdmb = 1.0 & ! non-dim sub grid mtn drag Amp (j)
864	! hncrit set to 8000m and sigfac added to enhance elvmax mtn hgt
865	&, hncrit=8000. & ! Max value in meters for ELVMAX (j)
866	!module top! &, sigfac=3.0 & ! MB3a expt test for ELVMAX factor (j) => control value is 0.1
867	!module top &, sigfac=0. & ! MB3a expt test for ELVMAX factor (j) => control value is 0.1
868	&, hminmt=50. & ! min mtn height (j)
869	&, hstdmin=25. & ! min orographic std dev in height
870	&, minwnd=0.1 & ! min wind component (j)
871	&, dpmin=5.0 ! Minimum thickness of the reference layer (centibars)
872	! values of dpmin=0, 20 have also been used
873	!
874	integer, parameter :: mdir=8
875	real(kind=kind_phys), parameter :: FDIR=mdir/(PI+PI)
876	!
877	!-- Template:
878	! NWD 1 2 3 4 5 6 7 8
879	! WD W S SW NW E N NE SE
880	!
881	integer,save :: nwdir(mdir)
882	data nwdir /6,7,5,8,2,3,1,4/
883	!
884	LOGICAL ICRILV(IM)
885	!
886	!---- MOUNTAIN INDUCED GRAVITY WAVE DRAG
887	!
888	!
889	! for lm mtn blocking
890	real(kind=kind_phys), DIMENSION(IM) :: WK,PE,EK,ZBK,UP,TAUB,XN &
891	& ,YN,UBAR,VBAR,ULOW,OA,CLX,ROLL,ULOI,DTFAC,XLINV,DELKS,DELKS1 &
892	& ,SCOR,BNV2bar, ELEVMX ! ,PSTAR
893	!
894	real(kind=kind_phys), DIMENSION(IM,KTS:KTE) :: &
895	& BNV2LM,DB,ANG,UDS,BNV2,RI_N,TAUD,RO,VTK,VTJ
896	real(kind=kind_phys), DIMENSION(IM,KTS:KTE-1) :: VELCO
897	real(kind=kind_phys), DIMENSION(IM,KTS:KTE+1) :: TAUP
898	real(kind=kind_phys), DIMENSION(KTE-1) :: VELKO
899	!
900	integer, DIMENSION(IM) :: &
901	& kref,kint,iwk,iwk2,ipt,kreflm,iwklm,iptlm,idxzb
902	!
903	! for lm mtn blocking
904	!
905	real(kind=kind_phys) :: ZLEN, DBTMP, R, PHIANG, DBIM &
906	&, xl, rcsks, bnv, fr &
907	&, brvf, cleff, tem, tem1, tem2, temc, temv &
908	&, wdir, ti, rdz, dw2, shr2, bvf2 &
909	&, rdelks, wtkbj, efact, coefm, gfobnv &
910	&, scork, rscor, hd, fro, rim, sira &
911	&, dtaux, dtauy, pkp1log, pklog
912	!
913	integer :: ncnt, kmm1, kmm2, lcap, lcapp1, kbps, kbpsp1,kbpsm1 &
914	&, kmps, kmpsp1, idir, nwd, i, j, k, klcap, kp1, kmpbl, npt, npr &
915	&, idxm1, ktrial, klevm1, kmll,kmds, KM &
916	! &, ihit,jhit &
917	&, ME !-- processor element for debugging
918
919	real :: rcl,rcs !dbg
920
921	!
922	!-----------------------------------------------------------------------
923	!
924	KM = KTE
925	npr = 0
926	DO I = 1, IM
927	DUsfc(I) = 0.
928	DVsfc(I) = 0.
929	!
930	!-- ELEVMX is a local array that could be changed below
931	!
932	ELEVMX(I) = ELVMAX(I)
933	ENDDO
934	!
935	!-- Note that A, B already set to zero as DUDTcol, DVDTcol in subroutine GWD_driver
936	!
937	ipt = 0
938	npt = 0
939	IF (NMTVR >= 14) then
940	DO I = 1,IM
941	IF (elvmax(i) > HMINMT .AND. hprime(i) > HE_4) then
942	npt = npt + 1
943	ipt(npt) = i
944	ENDIF
945	ENDDO
946	ELSE
947	DO I = 1,IM
948	IF (hprime(i) > HE_4) then
949	npt = npt + 1
950	ipt(npt) = i
951	ENDIF
952	ENDDO
953	ENDIF !-- IF (NMTVR >= 14) then
954	!
955
956	!dbg
957	rcl=1.
958	rcs=1.
959	!dbg if (lprnt) then
960	!dbg !-- Match what's in the GFS:
961	!dbg !dbg rcl=1.53028780126139008 ! match GFS point at 36N
962	!dbg !dbg rcs=sqrt(rcl)
963	!dbg i=im
964	!dbg write(6,"(a,a)") LABEL,' in GWD_col: K U V T Q Pi DP P EX PHIi PHI'
965	!dbg do k=kts,kte
966	!dbg write(6,"(i3,10e12.4)") k,U1(i,k),V1(i,k),T1(i,k),Q1(i,k),PRSI(i,k) &
967	!dbg ,DEL(i,k),PRSL(i,k),PRSLK(i,k),PHII(i,k),PHIL(i,k)
968	!dbg enddo
969	!dbg write(6,"(2(a,e12.4),2(a,4e12.4/),4(a,e12.4),a,i3) )") &
970	!dbg 'GWD_col: hprime(i)=',hprime(i),' oc(i)=',oc(i) &
971	!dbg ,' oa4(i,1-4)=',(oa4(i,k),k=1,4) &
972	!dbg ,' clx4(i,1-4)=',(clx4(i,k),k=1,4) &
973	!dbg ,' theta(i)=',theta(i),' sigma(i)=',sigma(i) &
974	!dbg ,' gamma(i)=',gamma(i),' elvmax(i)=',elvmax(i) &
975	!dbg ,' lpbl(i)=',kpbl(i)
976	!dbg endif
977	!dbg if (lprnt) CALL WRF_GET_MYPROC(ME)
978
979	!
980	!-- Note important criterion for immediately exiting routine!
981	!
982	IF (npt <= 0) RETURN ! No gwd/mb calculation done!
983	!
984	do i=1,npt
985	IDXZB(i) = 0
986	enddo
987	!
988	DO K = 1, KM
989	DO I = 1, IM
990	DB(I,K) = 0.
991	ANG(I,K) = 0.
992	UDS(I,K) = 0.
993	ENDDO
994	ENDDO
995	!
996	!
997	! NCNT = 0
998	KMM1 = KM - 1
999	KMM2 = KM - 2
1000	LCAP = KM
1001	LCAPP1 = LCAP + 1
1002	!
1003	!
1004	IF (NMTVR .eq. 14) then
1005	! ---- for lm and gwd calculation points
1006	!
1007	! --- iwklm is the level above the height of the mountain.
1008	! --- idxzb is the level of the dividing streamline.
1009	! INITIALIZE DIVIDING STREAMLINE (DS) CONTROL VECTOR
1010	!
1011	do i=1,npt
1012	iwklm(i) = 2
1013	kreflm(i) = 0
1014	enddo
1015	!
1016	!
1017	! start lm mtn blocking (mb) section
1018	!
1019	!..............................
1020	!..............................
1021	!
1022	! (j) 11/03: test upper limit on KMLL=km - 1
1023	! then do not need hncrit -- test with large hncrit first.
1024	! KMLL = km / 2 ! maximum mtnlm height : # of vertical levels / 2
1025	KMLL = kmm1
1026	! --- No mtn should be as high as KMLL (so we do not have to start at
1027	! --- the top of the model but could do calc for all levels).
1028	!
1029
1030	!dbg
1031	!dbg if (lprnt) print *,'k pkp1log pklog vtj(i,k) vtk(i,k) ro(i,k)'
1032
1033	DO I = 1, npt
1034	j = ipt(i)
1035	ELEVMX(J) = min (ELEVMX(J) + sigfac * hprime(j), hncrit)
1036
1037	!dbg
1038	!dbg if (lprnt) print *,'k=',k,' elevmx(j)=',elevmx(j),' elvmax(j)=',elvmax(j) &
1039	!dbg ,' sigfachprime(j)=',sigfachprime(j)
1040
1041	ENDDO
1042
1043	DO K = 1,KMLL
1044	DO I = 1, npt
1045	j = ipt(i)
1046	! --- interpolate to max mtn height for index, iwklm(I) wk[gz]
1047	! --- ELEVMX is limited to hncrit because to hi res topo30 orog.
1048	pkp1log = phil(j,k+1) * ROG
1049	pklog = phil(j,k) * ROG
1050	if ( ( ELEVMX(j) .le. pkp1log ) .and. &
1051	& ( ELEVMX(j) .ge. pklog ) ) THEN
1052	! --- wk for diags but can be saved and reused.
1053	wk(i) = G * ELEVMX(j) / ( phil(j,k+1) - phil(j,k) )
1054	iwklm(I) = MAX(iwklm(I), k+1 )
1055
1056	!dbg if (lprnt) print *,'k,wk(i),iwklm(i)=',k,wk(i),iwklm(i) !dbg
1057
1058	endif
1059	!
1060	! --- find at prsl levels large scale environment variables
1061	! --- these cover all possible mtn max heights
1062	VTJ(I,K) = T1(J,K) * (1.+FV*Q1(J,K)) ! virtual temperature
1063	VTK(I,K) = VTJ(I,K) / PRSLK(J,K) ! potential temperature
1064	RO(I,K) = RDI * PRSL(J,K) / VTJ(I,K) ! DENSITY (1.e-3 kg m^-3)
1065
1066	!dbg if (lprnt) write(6,"(i2,5e12.4)") k,pkp1log,pklog,vtj(i,k),vtk(i,k),ro(i,k) !dbg
1067
1068	ENDDO !-- DO I = 1, npt
1069	!
1070	ENDDO !-- DO K = 1,KMLL
1071	!
1072	! testing for highest model level of mountain top
1073	!
1074	! ihit = 2
1075	! jhit = 0
1076	! do i = 1, npt
1077	! j=ipt(i)
1078	! if ( iwklm(i) .gt. ihit ) then
1079	! ihit = iwklm(i)
1080	! jhit = j
1081	! endif
1082	! enddo
1083	! if (lprnt) print *, ' mb: kdt,max(iwklm),jhit,phil,me=', &
1084	! & kdt,ihit,jhit,phil(jhit,ihit),me
1085	!
1086	!dbg if (lprnt) print *,'k rdz bnv2lm(i,k)' !dbg
1087	klevm1 = KMLL - 1
1088	DO K = 1, klevm1
1089	DO I = 1, npt
1090	j = ipt(i)
1091	RDZ = g / ( phil(j,k+1) - phil(j,k) )
1092	! --- Brunt-Vaisala Frequency
1093	BNV2LM(I,K) = (G+G) * RDZ * ( VTK(I,K+1)-VTK(I,K) ) &
1094	& / ( VTK(I,K+1)+VTK(I,K) )
1095	bnv2lm(i,k) = max( bnv2lm(i,k), bnv2min )
1096
1097	!dbg if (lprnt) write(6,"(i2,2e12.4)") k,rdz,bnv2lm(i,k) !dbg
1098
1099	ENDDO
1100	ENDDO
1101	!
1102	DO I = 1, npt
1103	J = ipt(i)
1104	DELKS(I) = 1.0 / (PRSI(J,1) - PRSI(J,iwklm(i)))
1105	DELKS1(I) = 1.0 / (PRSL(J,1) - PRSL(J,iwklm(i)))
1106	UBAR (I) = 0.0
1107	VBAR (I) = 0.0
1108	ROLL (I) = 0.0
1109	PE (I) = 0.0
1110	EK (I) = 0.0
1111	BNV2bar(I) = (PRSL(J,1)-PRSL(J,2)) * DELKS1(I) * BNV2LM(I,1)
1112	ENDDO
1113	!
1114	! --- find the dividing stream line height
1115	! --- starting from the level above the max mtn downward
1116	! --- iwklm(i) is the k-index of mtn elevmx elevation
1117	!
1118	DO Ktrial = KMLL, 1, -1
1119	DO I = 1, npt
1120	IF ( Ktrial .LT. iwklm(I) .and. kreflm(I) .eq. 0 ) then
1121	kreflm(I) = Ktrial
1122
1123	if (lprnt) print *,'Ktrial,iwklm(i),kreflm(i)=',Ktrial,iwklm(i),kreflm(I)
1124
1125	ENDIF
1126	ENDDO
1127	ENDDO
1128	!
1129	! --- in the layer kreflm(I) to 1 find PE (which needs N, ELEVMX)
1130	! --- make averages, guess dividing stream (DS) line layer.
1131	! --- This is not used in the first cut except for testing and
1132	! --- is the vert ave of quantities from the surface to mtn top.
1133	!
1134
1135	!dbg
1136	!dbg if (lprnt) print *,' k rdelks ubar vbar roll ' &
1137	!dbg ,'bnv2bar rcsks rcs'
1138
1139	DO I = 1, npt
1140	DO K = 1, Kreflm(I)
1141	J = ipt(i)
1142	RDELKS = DEL(J,K) * DELKS(I)
1143
1144	!dbg
1145	RCSKS = RCS * RDELKS
1146	UBAR(I) = UBAR(I) + RCSKS * U1(J,K) ! trial Mean U below
1147	VBAR(I) = VBAR(I) + RCSKS * V1(J,K) ! trial Mean V below
1148
1149	ROLL(I) = ROLL(I) + RDELKS * RO(I,K) ! trial Mean RO below
1150	RDELKS = (PRSL(J,K)-PRSL(J,K+1)) * DELKS1(I)
1151	BNV2bar(I) = BNV2bar(I) + BNV2lm(I,K) * RDELKS
1152	! --- these vert ave are for diags, testing and GWD to follow (j).
1153
1154	!dbg
1155	!dbg if (lprnt) write(6,"(i2,7e12.4)") k,rdelks,ubar(i),vbar(i),roll(i) &
1156	!dbg ,bnv2bar(i),rcsks,rcs
1157
1158	ENDDO
1159	ENDDO
1160
1161	!dbg
1162	!dbg if (lprnt) print *, 'kreflm(npt)=',kreflm(npt) &
1163	!dbg ,' bnv2bar(npt)=',bnv2bar(npt) &
1164	!dbg ,' ubar(npt)=',ubar(npt) &
1165	!dbg ,' vbar(npt)=',vbar(npt) &
1166	!dbg ,' roll(npt)=',roll(npt) &
1167	!dbg ,' delks(npt)=',delks(npt) &
1168	!dbg ,' delks1(npt)=',delks1(npt)
1169
1170	!
1171	! --- integrate to get PE in the trial layer.
1172	! --- Need the first layer where PE>EK - as soon as
1173	! --- IDXZB is not 0 we have a hit and Zb is found.
1174	!
1175	DO I = 1, npt
1176	J = ipt(i)
1177
1178	!dbg
1179	!dbg if (lprnt) print *,'k phiang u1(j,k) v1(j,k) theta(j)' &
1180	!dbg ,' ang(i,k) uds(i,k) pe(i) up(i) ek(i)'
1181
1182	DO K = iwklm(I), 1, -1
1183	PHIANG = RAD_TO_DEG*atan2(V1(J,K),U1(J,K))
1184	ANG(I,K) = ( THETA(J) - PHIANG )
1185	if ( ANG(I,K) .gt. 90. ) ANG(I,K) = ANG(I,K) - 180.
1186	if ( ANG(I,K) .lt. -90. ) ANG(I,K) = ANG(I,K) + 180.
1187	!
1188	!dbg UDS(I,K) = &
1189	UDS(I,K) = rcs* &
1190	& MAX(SQRT(U1(J,K)U1(J,K) + V1(J,K)V1(J,K)), minwnd)
1191	! --- Test to see if we found Zb previously
1192	IF (IDXZB(I) .eq. 0 ) then
1193	PE(I) = PE(I) + BNV2lm(I,K) * &
1194	& ( G * ELEVMX(J) - phil(J,K) ) * &
1195	& ( PHII(J,K+1) - PHII(J,K) ) * ROG2
1196
1197	!dbg
1198	!dbg & ( PHII(J,K+1) - PHII(J,K) ) / (G*G)
1199	!dbg if (lprnt) print *, &
1200	!dbg 'k=',k,' pe(i)=',pe(i),' bnv2lm(i,k)=',bnv2lm(i,k) &
1201	!dbg ,' g=',g,' elevmx(j)=',elevmx(j) &
1202	!dbg ,' gelevmx(j)-phil(j,k)=',gelevmx(j)-phil(j,k) &
1203	!dbg ,' (phii(j,k+1)-phi(j,k))/(gg)=',(PHII(J,K+1)-PHII(J,K) )ROG2
1204
1205	! --- KE
1206	! --- Wind projected on the line perpendicular to mtn range, U(Zb(K)).
1207	! --- kinetic energy is at the layer Zb
1208	! --- THETA ranges from -+90deg \|_ to the mtn "largest topo variations"
1209	UP(I) = UDS(I,K) * cos(DEG_TO_RAD*ANG(I,K))
1210	EK(I) = 0.5 * UP(I) * UP(I)
1211
1212	! --- Dividing Stream lime is found when PE =exceeds EK.
1213	IF ( PE(I) .ge. EK(I) ) IDXZB(I) = K
1214	! --- Then mtn blocked flow is between Zb=k(IDXZB(I)) and surface
1215	!
1216	ENDIF !-- IF (IDXZB(I) .eq. 0 ) then
1217
1218	!dbg
1219	!dbg if (lprnt) write(6,"(i2,9e12.4)") &
1220	!dbg k,phiang,u1(j,k),v1(j,k),theta(j),ang(i,k),uds(i,k),pe(i),up(i),ek(i)
1221
1222	ENDDO !-- DO K = iwklm(I), 1, -1
1223	ENDDO !-- DO I = 1, npt
1224	!
1225	DO I = 1, npt
1226	J = ipt(i)
1227	! --- Calc if N constant in layers (Zb guess) - a diagnostic only.
1228	ZBK(I) = ELEVMX(J) - SQRT(UBAR(I)2 + VBAR(I)2)/BNV2bar(I)
1229	ENDDO
1230	!
1231
1232	!dbg
1233	!dbg if (lprnt) print *,'iwklm=',iwklm(npt),' ZBK=',ZBK(npt)
1234
1235	!
1236	! --- The drag for mtn blocked flow
1237	!
1238
1239	!dbg
1240	!dbg if (lprnt) print ,'k phil(j,k) ghprime(j) ' &
1241	!dbg ,'phil(j,idxzb(i))-phil(j,k) zlen r dbtmp db(i,k)'
1242
1243	!
1244	DO I = 1, npt
1245	J = ipt(i)
1246	ZLEN = 0.
1247	IF ( IDXZB(I) .gt. 0 ) then
1248	DO K = IDXZB(I), 1, -1
1249	IF (PHIL(J,IDXZB(I)) > PHIL(J,K)) THEN
1250	ZLEN = SQRT( ( PHIL(J,IDXZB(I))-PHIL(J,K) ) / &
1251	& ( PHIL(J,K ) + G * hprime(J) ) )
1252	! --- lm eq 14:
1253	R = (cos(DEG_TO_RADANG(I,K))2 + GAMMA(J) sin(DEG_TO_RADANG(I,K))*2) / &
1254	& (gamma(J) * cos(DEG_TO_RADANG(I,K))2 + sin(DEG_TO_RADANG(I,K))**2)
1255	! --- (negative of DB -- see sign at tendency)
1256	DBTMP = 0.25 * CDmb * &
1257	& MAX( 2. - 1. / R, HZERO ) * sigma(J) * &
1258	& MAX(cos(DEG_TO_RADANG(I,K)), gamma(J)sin(DEG_TO_RADANG(I,K))) &
1259	& ZLEN / hprime(J)
1260	DB(I,K) = DBTMP * UDS(I,K)
1261	!
1262
1263	!dbg
1264	!dbg if (lprnt) write(6,"(i3,7e12.4)") &
1265	!dbg k,phil(j,k),g*hprime(j),phil(j,idxzb(i))-phil(j,k),zlen,r,dbtmp,db(i,k)
1266
1267	!
1268	ENDIF !-- IF (PHIL(J,IDXZB(I)) > PHIL(J,K) .AND. DEN > 0.) THEN
1269	ENDDO !-- DO K = IDXZB(I), 1, -1
1270	endif
1271	ENDDO !-- DO I = 1, npt
1272	!
1273	!.............................
1274	!.............................
1275	! end mtn blocking section
1276	!
1277	ENDIF !-- IF ( NMTVR .eq. 14) then
1278	!
1279	!.............................
1280	!.............................
1281	!
1282	KMPBL = km / 2 ! maximum pbl height : # of vertical levels / 2
1283	!
1284	! Scale cleff between IM=3842 and 1922 for T126/T170 and T62
1285	!
1286	if (imx .gt. 0) then
1287	! cleff = 1.0E-5 * SQRT(FLOAT(IMX)/384.0) ! this is inverse of CLEFF!
1288	! cleff = 1.0E-5 * SQRT(FLOAT(IMX)/192.0) ! this is inverse of CLEFF!
1289	cleff = 0.5E-5 * SQRT(FLOAT(IMX)/192.0) ! this is inverse of CLEFF!
1290	! cleff = 2.0E-5 * SQRT(FLOAT(IMX)/192.0) ! this is inverse of CLEFF!
1291	! cleff = 2.5E-5 * SQRT(FLOAT(IMX)/192.0) ! this is inverse of CLEFF!
1292	endif
1293
1294	!dbg
1295	!dbg if (lprnt) print *,'imx, cleff, rcl, rcs=',imx,cleff,rcl,rcs
1296	!dbg if (lprnt) print *,' k vtj vtk ro'
1297
1298	DO K = 1,KM
1299	DO I =1,npt
1300	J = ipt(i)
1301	VTJ(I,K) = T1(J,K) * (1.+FV*Q1(J,K))
1302	VTK(I,K) = VTJ(I,K) / PRSLK(J,K)
1303	RO(I,K) = RDI * PRSL(J,K) / VTJ(I,K) ! DENSITY
1304	TAUP(I,K) = 0.0
1305
1306	!dbg
1307	!dbg if (lprnt) write(6,"(i2,3e12.4)") k,vtj(i,k),vtk(i,k),ro(i,k) !dbg
1308
1309	ENDDO
1310	ENDDO
1311
1312	!dbg
1313	!dbg if (lprnt) print *,' k ti tem rdz tem1' &
1314	!dbg ,' tem2 dw2 shr2 bvf2 ri_n(i,k) bnv2(i,k)'
1315
1316	DO K = 1,KMM1
1317	DO I =1,npt
1318	J = ipt(i)
1319	TI = 2.0 / (T1(J,K)+T1(J,K+1))
1320	TEM = TI / (PRSL(J,K)-PRSL(J,K+1))
1321	! RDZ = GOR * PRSI(J,K+1) * TEM
1322	RDZ = g / (phil(j,k+1) - phil(j,k))
1323	TEM1 = U1(J,K) - U1(J,K+1)
1324	TEM2 = V1(J,K) - V1(J,K+1)
1325
1326	!dbg
1327	! DW2 = TEM1TEM1 + TEM2TEM2
1328	DW2 = rcl(TEM1TEM1 + TEM2*TEM2)
1329
1330	SHR2 = MAX(DW2,DW2MIN) * RDZ * RDZ
1331	BVF2 = G(GOCP+RDZ(VTJ(I,K+1)-VTJ(I,K))) * TI
1332	ri_n(I,K) = MAX(BVF2/SHR2,RIMIN) ! Richardson number
1333	! Brunt-Vaisala Frequency
1334	! TEM = GR2 * (PRSL(J,K)+PRSL(J,K+1)) * TEM
1335	! BNV2(I,K) = TEM * (VTK(I,K+1)-VTK(I,K))/(VTK(I,K+1)+VTK(I,K))
1336	BNV2(I,K) = (G+G) * RDZ * (VTK(I,K+1)-VTK(I,K)) &
1337	& / (VTK(I,K+1)+VTK(I,K))
1338	bnv2(i,k) = max( bnv2(i,k), bnv2min )
1339
1340	!dbg
1341	!dbg if (lprnt) write(6,"(i2,10e12.4)") &
1342	!dbg k,ti,tem,rdz,tem1,tem2,dw2,shr2,bvf2,ri_n(i,k),bnv2(i,k)
1343
1344	!
1345	ENDDO !-- DO K = 1,KMM1
1346	ENDDO !-- DO I =1,npt
1347	!
1348
1349	!dbg
1350	!dbg if (lprnt) print *,'kmm1,bnv2=',kmm1,bnv2(npt,kmm1)
1351
1352	!
1353	! Apply 3 point smoothing on BNV2
1354	!
1355	! do k=1,km
1356	! do i=1,im
1357	! vtk(i,k) = bnv2(i,k)
1358	! enddo
1359	! enddo
1360	! do k=2,kmm1
1361	! do i=1,im
1362	! bnv2(i,k) = 0.25(vtk(i,k-1)+vtk(i,k+1)) + 0.5vtk(i,k)
1363	! enddo
1364	! enddo
1365	!
1366	! Finding the first interface index above 50 hPa level
1367	!
1368	do i=1,npt
1369	iwk(i) = 2
1370	enddo
1371
1372	!dbg if (lprnt) print *,'kmpbl=',kmpbl !dbg
1373
1374	DO K=3,KMPBL
1375	DO I=1,npt
1376	j = ipt(i)
1377	tem = (prsi(j,1) - prsi(j,k))
1378	if (tem .lt. dpmin) iwk(i) = k
1379	enddo
1380	enddo
1381	!
1382	KBPS = 1
1383	KMPS = KM
1384	DO I=1,npt
1385	J = ipt(i)
1386	kref(I) = MAX(IWK(I), KPBL(J)+1 ) ! reference level
1387	DELKS(I) = 1.0 / (PRSI(J,1) - PRSI(J,kref(I)))
1388	DELKS1(I) = 1.0 / (PRSL(J,1) - PRSL(J,kref(I)))
1389	UBAR (I) = 0.0
1390	VBAR (I) = 0.0
1391	ROLL (I) = 0.0
1392	KBPS = MAX(KBPS, kref(I))
1393	KMPS = MIN(KMPS, kref(I))
1394	!
1395	BNV2bar(I) = (PRSL(J,1)-PRSL(J,2)) * DELKS1(I) * BNV2(I,1)
1396	ENDDO
1397	!
1398	!
1399	KBPSP1 = KBPS + 1
1400	KBPSM1 = KBPS - 1
1401	DO K = 1,KBPS
1402	DO I = 1,npt
1403	IF (K .LT. kref(I)) THEN
1404	J = ipt(i)
1405	RDELKS = DEL(J,K) * DELKS(I)
1406
1407	!dbg
1408	! UBAR(I) = UBAR(I) + RDELKS * U1(J,K) ! Mean U below kref
1409	! VBAR(I) = VBAR(I) + RDELKS * V1(J,K) ! Mean V below kref
1410	RCSKS = RCS * RDELKS
1411	UBAR(I) = UBAR(I) + RCSKS * U1(J,K) ! Mean U below kref
1412	VBAR(I) = VBAR(I) + RCSKS * V1(J,K) ! Mean V below kref
1413
1414	!
1415	ROLL(I) = ROLL(I) + RDELKS * RO(I,K) ! Mean RO below kref
1416	RDELKS = (PRSL(J,K)-PRSL(J,K+1)) * DELKS1(I)
1417	BNV2bar(I) = BNV2bar(I) + BNV2(I,K) * RDELKS
1418	ENDIF
1419	ENDDO
1420	ENDDO
1421	!
1422
1423	!dbg
1424	!dbg if(lprnt) print *,'ubar(npt)=',ubar(npt) &
1425	!dbg ,' vbar(npt)=',vbar(npt),' roll(npt)=',roll(npt) &
1426	!dbg ,' bnv2bar(npt)=',bnv2bar(npt)
1427
1428	!
1429	! FIGURE OUT LOW-LEVEL HORIZONTAL WIND DIRECTION AND FIND 'OA'
1430	!
1431	! NWD 1 2 3 4 5 6 7 8
1432	! WD W S SW NW E N NE SE
1433	!
1434	DO I = 1,npt
1435	J = ipt(i)
1436	wdir = atan2(UBAR(I),VBAR(I)) + pi
1437	idir = mod(nint(fdir*wdir),mdir) + 1
1438	nwd = nwdir(idir)
1439	OA(I) = (1-2INT( (NWD-1)/4 )) OA4(J,MOD(NWD-1,4)+1)
1440	CLX(I) = CLX4(J,MOD(NWD-1,4)+1)
1441	ENDDO
1442	!
1443	!-----XN,YN "LOW-LEVEL" WIND PROJECTIONS IN ZONAL
1444	! & MERIDIONAL DIRECTIONS
1445	!-----ULOW "LOW-LEVEL" WIND MAGNITUDE - (= U)
1446	!-----BNV2 BNV2 = N**2
1447	!-----TAUB BASE MOMENTUM FLUX
1448	!-----= -(RO * U*3/(NXL)GF(FR) FOR N*2 > 0
1449	!-----= 0. FOR N**2 < 0
1450	!-----FR FROUDE = N*HPRIME / U
1451	!-----G GMAXFR2/(FR*2+CG/OC)
1452	!
1453	!-----INITIALIZE SOME ARRAYS
1454	!
1455	DO I = 1,npt
1456	XN(I) = 0.0
1457	YN(I) = 0.0
1458	TAUB (I) = 0.0
1459	ULOW (I) = 0.0
1460	DTFAC(I) = 1.0
1461	ICRILV(I) = .FALSE. ! INITIALIZE CRITICAL LEVEL CONTROL VECTOR
1462	!
1463	!----COMPUTE THE "LOW LEVEL" WIND MAGNITUDE (M/S)
1464	!
1465	ULOW(I) = MAX(SQRT(UBAR(I)UBAR(I) + VBAR(I)VBAR(I)), HONE)
1466	ULOI(I) = 1.0 / ULOW(I)
1467	ENDDO
1468
1469	!dbg
1470	!dbg if (lprnt) print *,'idir,nwd,wdir,oa,clx,ulow,uloi=' &
1471	!dbg ,idir,nwd,wdir,oa(npt),clx(npt),ulow(npt),uloi(npt)
1472
1473	!
1474	DO K = 1,KMM1
1475	DO I = 1,npt
1476	J = ipt(i)
1477
1478	!dbg
1479	! VELCO(I,K) = 0.5* ((U1(J,K)+U1(J,K+1))*UBAR(I) &
1480	VELCO(I,K) = 0.5rcs((U1(J,K)+U1(J,K+1))*UBAR(I) &
1481	& + (V1(J,K)+V1(J,K+1))*VBAR(I))
1482
1483	VELCO(I,K) = VELCO(I,K) * ULOI(I)
1484
1485	!dbg if (lprnt) write(6,"(a,i2,a,e12.4)") 'k=',k,' velco(i,k)=',velco(i,k) !dbg
1486
1487	! IF ((VELCO(I,K).LT.VELEPS) .AND. (VELCO(I,K).GT.0.)) THEN
1488	! VELCO(I,K) = VELEPS
1489	! ENDIF
1490	ENDDO
1491	ENDDO
1492	!
1493	! find the interface level of the projected wind where
1494	! low levels & upper levels meet above pbl
1495	!
1496	do i=1,npt
1497	kint(i) = km
1498	enddo
1499	do k = 1,kmm1
1500	do i = 1,npt
1501	IF (K .GT. kref(I)) THEN
1502	if(velco(i,k) .lt. veleps .and. kint(i) .eq. km) then
1503	kint(i) = k+1
1504	endif
1505	endif
1506	enddo
1507	enddo
1508	! WARNING KINT = KREF !!!!!!!!!
1509	do i=1,npt
1510	kint(i) = kref(i)
1511	enddo
1512	!
1513	!
1514	DO I = 1,npt
1515	J = ipt(i)
1516	BNV = SQRT( BNV2bar(I) )
1517	FR = BNV * ULOI(I) * min(HPRIME(J),hpmax)
1518	FR = MIN(FR, FRMAX)
1519	XN(I) = UBAR(I) * ULOI(I)
1520	YN(I) = VBAR(I) * ULOI(I)
1521	!
1522	! Compute the base level stress and store it in TAUB
1523	! CALCULATE ENHANCEMENT FACTOR, NUMBER OF MOUNTAINS & ASPECT
1524	! RATIO CONST. USE SIMPLIFIED RELATIONSHIP BETWEEN STANDARD
1525	! DEVIATION & CRITICAL HGT
1526	!
1527	EFACT = (OA(I) + 2.) ** (CEOFRC*FR)
1528	EFACT = MIN( MAX(EFACT,EFMIN), EFMAX )
1529	!
1530	COEFM = (1. + CLX(I)) ** (OA(I)+1.)
1531	!
1532	XLINV(I) = COEFM * CLEFF
1533	!
1534	TEM = FR * FR * OC(J)
1535	GFOBNV = GMAX * TEM / ((TEM + CG)*BNV) ! G/N0
1536	!
1537	TAUB(I) = XLINV(I) * ROLL(I) * ULOW(I) * ULOW(I) &
1538	& * ULOW(I) * GFOBNV * EFACT ! BASE FLUX Tau0
1539	!
1540	! tem = min(HPRIME(I),hpmax)
1541	! TAUB(I) = XLINV(I) * ROLL(I) * ULOW(I) * BNV * tem * tem
1542	!
1543	K = MAX(1, kref(I)-1)
1544	TEM = MAX(VELCO(I,K)*VELCO(I,K), HE_4)
1545	SCOR(I) = BNV2(I,K) / TEM ! Scorer parameter below ref level
1546	ENDDO !-- DO I = 1,npt
1547	!
1548
1549	!dbg
1550	!dbg if (lprnt) write(6,"(a,i2,10(a,e12.4))") &
1551	!dbg 'kint=',kint(npt),' bnv=',bnv,' fr=',fr,' xn=',xn(npt) &
1552	!dbg ,' yn=',yn(npt),' efact=',efact,' coefm=',coefm,' xlinv(npt)=',xlinv(npt) &
1553	!dbg ,' gfobnv=',gfobnv,' taub(npt)=',taub(npt),' scor(npt)=',scor(npt)
1554
1555	!
1556	!----SET UP BOTTOM VALUES OF STRESS
1557	!
1558	DO K = 1, KBPS
1559	DO I = 1,npt
1560	IF (K .LE. kref(I)) TAUP(I,K) = TAUB(I)
1561	ENDDO
1562	ENDDO
1563
1564	!
1565	! Now compute vertical structure of the stress.
1566	!
1567	DO K = KMPS, KMM1 ! Vertical Level K Loop!
1568	KP1 = K + 1
1569	DO I = 1, npt
1570	!
1571	!-----UNSTABLE LAYER IF RI < RIC
1572	!-----UNSTABLE LAYER IF UPPER AIR VEL COMP ALONG SURF VEL <=0 (CRIT LAY)
1573	!---- AT (U-C)=0. CRIT LAYER EXISTS AND BIT VECTOR SHOULD BE SET (.LE.)
1574	!
1575	IF (K .GE. kref(I)) THEN
1576	ICRILV(I) = ICRILV(I) .OR. ( ri_n(I,K) .LT. RIC) &
1577	& .OR. (VELCO(I,K) .LE. 0.0)
1578	ENDIF
1579	ENDDO
1580	!
1581	DO I = 1,npt
1582	IF (K .GE. kref(I)) THEN
1583
1584	!dbg
1585	!dbg if (lprnt) write(6,"(2(a,i2),a,L1,3(a,e12.4))") 'k=',k,' kref(i)=',kref(i) &
1586	!dbg ,' icrilv(i)=',icrilv(i),' taup(i,k)=',taup(i,k) &
1587	!dbg ,' ri_n(i,k)=',ri_n(i,k),' velco(i,k)=',velco(i,k)
1588
1589	!
1590	IF (.NOT.ICRILV(I) .AND. TAUP(I,K) .GT. 0.0 ) THEN
1591	TEMV = 1.0 / max(VELCO(I,K), HE_2)
1592	! IF (OA(I) .GT. 0. .AND. PRSI(ipt(i),KP1).GT.RLOLEV) THEN
1593	IF (OA(I).GT.0. .AND. kp1 .lt. kint(i)) THEN
1594	SCORK = BNV2(I,K) * TEMV * TEMV
1595	RSCOR = MIN(HONE, SCORK / SCOR(I))
1596	SCOR(I) = SCORK
1597	ELSE
1598	RSCOR = 1.
1599	ENDIF
1600	!
1601	BRVF = SQRT(BNV2(I,K)) ! Brunt-Vaisala Frequency
1602	! TEM1 = XLINV(I)(RO(I,KP1)+RO(I,K))BRVFVELCO(I,K)0.5
1603	TEM1 = XLINV(I)(RO(I,KP1)+RO(I,K))BRVF*0.5 &
1604	& * max(VELCO(I,K),HE_2)
1605	HD = SQRT(TAUP(I,K) / TEM1)
1606	FRO = BRVF * HD * TEMV
1607	!
1608	! RIM is the MINIMUM-RICHARDSON NUMBER BY SHUTTS (1985)
1609	!
1610	TEM2 = SQRT(ri_n(I,K))
1611	TEM = 1. + TEM2 * FRO
1612	RIM = ri_n(I,K) * (1.-FRO) / (TEM * TEM)
1613	!
1614	! CHECK STABILITY TO EMPLOY THE SATURATION HYPOTHESIS
1615	! OF LINDZEN (1981) EXCEPT AT TROPOSPHERIC DOWNSTREAM REGIONS
1616	!
1617	! ----------------------
1618
1619	!dbg
1620	!dbg if (lprnt) write(6,"(a,i2,10(a,e12.4))") &
1621	!dbg 'k=',k,' brvf=',brvf,' tem1=',tem1,' xlinv(i)=',xlinv(i) &
1622	!dbg ,'ROavg=',.5*(RO(I,KP1)+RO(I,K)),' hd=',hd,' temv=',temv,' fro=',fro &
1623	!dbg ,' tem2=',tem2,' tem=',tem,' rim=',rim
1624
1625	!
1626	IF (RIM .LE. RIC .AND. &
1627	! & (OA(I) .LE. 0. .OR. PRSI(ipt(I),KP1).LE.RLOLEV )) THEN
1628	& (OA(I) .LE. 0. .OR. kp1 .ge. kint(i) )) THEN
1629	TEMC = 2.0 + 1.0 / TEM2
1630	HD = VELCO(I,K) * (2.*SQRT(TEMC)-TEMC) / BRVF
1631	TAUP(I,KP1) = TEM1 * HD * HD
1632	ELSE
1633	TAUP(I,KP1) = TAUP(I,K) * RSCOR
1634	ENDIF
1635	taup(i,kp1) = min(taup(i,kp1), taup(i,k))
1636
1637	!dbg
1638	!dbg if (lprnt) write(6,"(a,i2,2(a,e12.4))") 'kp1=',kp1 &
1639	!dbg ,' taup(i,kp1)=',taup(i,kp1),' taup(i,k)=',taup(i,k)
1640
1641	ENDIF !-- IF (.NOT.ICRILV(I) .AND. TAUP(I,K) .GT. 0.0 ) THEN
1642	ENDIF !-- IF (K .GE. kref(I)) THEN
1643	ENDDO !-- DO I = 1,npt
1644	ENDDO !-- DO K = KMPS, KMM1
1645	!
1646	! DO I=1,IM
1647	! taup(i,km+1) = taup(i,km)
1648	! ENDDO
1649	!
1650	IF(LCAP .LE. KM) THEN
1651
1652	!dbg if (lprnt) print *,'lcap,lcapp1,km=',lcap,lcapp1,km !dbg
1653
1654	DO KLCAP = LCAPP1, KM+1
1655	DO I = 1,npt
1656	SIRA = PRSI(ipt(I),KLCAP) / PRSI(ipt(I),LCAP)
1657	TAUP(I,KLCAP) = SIRA * TAUP(I,LCAP)
1658
1659	!dbg
1660	!dbg if (lprnt) write(6,"(a,i2,5(a,e12.4))") &
1661	!dbg 'klcap=',klcap,' sira=',sira,' prsi(ipt(i),klcap)=', prsi(ipt(i),klcap) &
1662	!dbg ,' prsi(ipt(i),lcap)=',prsi(ipt(i),lcap),' taup(i,lcap)=',taup(i,lcap) &
1663	!dbg ,' taup(i,klcap)=',taup(i,klcap)
1664
1665	!
1666	ENDDO
1667	ENDDO
1668	ENDIF
1669	!
1670	! Calculate - (g/p)d(tau)/d(sigma) and Decel terms DTAUX, DTAUY
1671	!
1672	DO I=1,npt
1673	! SCOR(I) = 1.0 / PSTAR(I)
1674
1675	!dbg
1676	! SCOR(I) = 1.0
1677	SCOR(I) = 1.0/RCS
1678
1679	ENDDO
1680	DO K = 1,KM
1681	DO I = 1,npt
1682	TAUD(I,K) = G * (TAUP(I,K+1) - TAUP(I,K)) * SCOR(I) &
1683	& / DEL(ipt(I),K)
1684
1685	!dbg
1686	!dbg if (lprnt) write(6,"(a,i2,4(a,e12.4))") 'k=',k,' taud(i,k)=',taud(i,k) &
1687	!dbg ,' D(taup)=',TAUP(I,K+1)-TAUP(I,K),' del(ipt(i),k)=',del(ipt(i),k) &
1688	!dbg ,' scor(i)=',scor(i)
1689
1690	ENDDO
1691	ENDDO
1692
1693	!
1694	!------LIMIT DE-ACCELERATION (MOMENTUM DEPOSITION ) AT TOP TO 1/2 VALUE
1695	!------THE IDEA IS SOME STUFF MUST GO OUT THE TOP
1696	!
1697	DO KLCAP = LCAP, KM
1698	DO I = 1,npt
1699	TAUD(I,KLCAP) = TAUD(I,KLCAP) * FACTOP
1700
1701	!dbg
1702	!dbg if (lprnt) write(6,"(a,i2,a,e12.4)") 'klcap=',klcap,' taud(i,klcap)=',taud(i,klcap)
1703
1704	ENDDO
1705	ENDDO
1706	!
1707	!------IF THE GRAVITY WAVE DRAG WOULD FORCE A CRITICAL LINE IN THE
1708	!------LAYERS BELOW SIGMA=RLOLEV DURING THE NEXT DELTIM TIMESTEP,
1709	!------THEN ONLY APPLY DRAG UNTIL THAT CRITICAL LINE IS REACHED.
1710	!
1711	DO K = 1,KMM1
1712	DO I = 1,npt
1713	IF (K .GT. kref(I) .and. PRSI(ipt(i),K) .GE. RLOLEV) THEN
1714	IF(TAUD(I,K).NE.0.) THEN
1715
1716	!dbg
1717	! TEM = DELTIM * TAUD(I,K)
1718	TEM = rcsDELTIM TAUD(I,K)
1719
1720	DTFAC(I) = MIN(DTFAC(I),ABS(VELCO(I,K)/TEM))
1721
1722	!dbg
1723	!dbg if (lprnt) write(6,"(a,i2,2(a,e12.4))") 'k=',k &
1724	!dbg ,' tem=',tem,' dtfac(i)=',dtfac(i)
1725
1726	ENDIF
1727	ENDIF
1728	ENDDO
1729	ENDDO
1730	! if(lprnt .and. npr .gt. 0) then
1731	! print *, before A=,A(npr,:)
1732	! print *, before B=,B(npr,:)
1733	! endif
1734	!
1735
1736	!dbg
1737	!dbg if (lprnt) write(6,"(a,i2,3(a,e12.4))") 'idxzb(npt)=',idxzb(npt) &
1738	!dbg ,' dtfac=',dtfac(npt),' xn=',xn(npt),' yn=',yn(npt)
1739
1740	!
1741	DO K = 1,KM
1742	DO I = 1,npt
1743	J = ipt(i)
1744	TAUD(I,K) = TAUD(I,K) * DTFAC(I)
1745	DTAUX = TAUD(I,K) * XN(I)
1746	DTAUY = TAUD(I,K) * YN(I)
1747	! --- lm mb (j) changes overwrite GWD
1748	if ( K .lt. IDXZB(I) .AND. IDXZB(I) .ne. 0 ) then
1749	DBIM = DB(I,K) / (1.+DB(I,K)*DELTIM)
1750	A(J,K) = - DBIM * V1(J,K) + A(J,K)
1751	B(J,K) = - DBIM * U1(J,K) + B(J,K)
1752	DUsfc(J) = DUsfc(J) - DBIM * V1(J,K) * DEL(J,K)
1753	DVsfc(J) = DVsfc(J) - DBIM * U1(J,K) * DEL(J,K)
1754
1755	!dbg
1756	!dbg if (lprnt .and. dbim > 0.) write(6,"(a,i2,4(a,e12.4))") &
1757	!dbg 'k=',k,' db(i,k)=',db(i,k),' dbim=',dbim &
1758	!dbg ,' dudt=',-dbimu1(j,k),' dvdt=',-dbimv1(j,k)
1759
1760	!
1761	else
1762	!
1763	A(J,K) = DTAUY + A(J,K)
1764	B(J,K) = DTAUX + B(J,K)
1765	DUsfc(J) = DUsfc(J) + DTAUX * DEL(J,K)
1766	DVsfc(J) = DVsfc(J) + DTAUY * DEL(J,K)
1767
1768	!dbg
1769	!dbg if (lprnt .and. dtaux+dtauy/=0.) write(6,"(a,i2,2(a,e12.4))") &
1770	!dbg ',k=',k,' dudt=dtaux=',dtaux,' dvdt=dtauy=',dtauy
1771
1772	endif
1773
1774	!dbg
1775	!dbg if (lprnt) write(6,"(a,i2,2(a,e12.4))") 'k=',k &
1776	!dbg ,' dusfc(j)=',dusfc(j),' dvsfc(j)=',dvsfc(j)
1777
1778	ENDDO !-- DO I = 1,npt
1779	ENDDO !-- DO K = 1,KM
1780	! if (lprnt) then
1781	! print *, in gwdps_lm.f after A=,A(ipr,:)
1782	! print *, in gwdps_lm.f after B=,B(ipr,:)
1783	! print *, DB=,DB(ipr,:)
1784	! endif
1785	DO I = 1,npt
1786	J = ipt(i)
1787	! TEM = (-1.E3/G) * PSTAR(I)
1788
1789	!dbg
1790	! TEM = -1.E3/G
1791	TEM = -1.E3ROGrcs
1792	DUsfc(J) = TEM * DUsfc(J)
1793	DVsfc(J) = TEM * DVsfc(J)
1794
1795	!dbg
1796	!dbg if (lprnt .and. i.eq.npr) write(6,"(3(a,e12.4))") 'tem=',tem &
1797	!dbg ,' dusfc(j)=',dusfc(j),' dvsfc(j)=',dvsfc(j)
1798
1799	ENDDO
1800	!
1801	END SUBROUTINE GWD_col
1802	!
1803	!#######################################################################
1804	!
1805	END MODULE module_gwd

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format