Context Navigation

hines_gwd.F @ 1961

Last change on this file since 1961 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:executable set to *
Property svn:keywords set to Author Date Id Revision

File size: 69.5 KB

Line
1	!
2	! $Id: hines_gwd.F 1907 2013-11-26 13:10:46Z fhourdin $
3	!
4	SUBROUTINE HINES_GWD(NLON,NLEV,DTIME,paphm1x, papm1x,
5	I rlat,tx,ux,vx,
6	O zustrhi,zvstrhi,
7	O d_t_hin, d_u_hin, d_v_hin)
8
9	C ########################################################################
10	C Parametrization of the momentum flux deposition due to a broad band
11	C spectrum of gravity waves, following Hines (1997a,b), as coded by
12	C McLANDRESS (1995). Modified by McFARLANE and MANZINI (1995-1997)
13	C MAECHAM model stand alone version
14	C ########################################################################
15	C
16	C
17	USE dimphy
18	implicit none
19
20	cym#include "dimensions.h"
21	cym#include "dimphy.h"
22	#include "YOEGWD.h"
23	#include "YOMCST.h"
24
25	INTEGER NAZMTH
26	PARAMETER(NAZMTH=8)
27
28	C INPUT ARGUMENTS.
29	C ----- ----------
30	C
31	C - 2D
32	C PAPHM1 : HALF LEVEL PRESSURE (T-DT)
33	C PAPM1 : FULL LEVEL PRESSURE (T-DT)
34	C PTM1 : TEMPERATURE (T-DT)
35	C PUM1 : ZONAL WIND (T-DT)
36	C PVM1 : MERIDIONAL WIND (T-DT)
37	C
38
39	C REFERENCE.
40	C ----------
41	C SEE MODEL DOCUMENTATION
42	C
43	C AUTHOR.
44	C -------
45	C
46	C N. MCFARLANE DKRZ-HAMBURG MAY 1995
47	C STAND ALONE E. MANZINI MPI-HAMBURG FEBRUARY 1997
48	C
49	C BASED ON A COMBINATION OF THE OROGRAPHIC SCHEME BY N.MCFARLANE 1987
50	C AND THE HINES SCHEME AS CODED BY C. MCLANDRESS 1995.
51	C
52	C
53	C
54	cym INTEGER KLEVM1
55	C
56	REAL PAPHM1(klon,klev+1), PAPM1(klon,KLEV)
57	REAL PTM1(klon,KLEV), PUM1(klon,KLEV), PVM1(klon,KLEV)
58	REAL PRFLUX(klon)
59	C1
60	C1
61	C1
62	REAL RLAT(klon),COSLAT(KLON)
63	C
64	REAL TH(klon,KLEV),
65	2 UTENDGW(klon,KLEV), VTENDGW(klon,KLEV),
66	3 PRESSG(klon),
67	4 UHS(klon,KLEV), VHS(klon,KLEV), ZPR(klon)
68
69	C * VERTICAL POSITIONING ARRAYS.
70
71	REAL SGJ(klon,KLEV), SHJ(klon,KLEV),
72	1 SHXKJ(klon,KLEV), DSGJ(klon,KLEV)
73
74	C * LOGICAL SWITCHES TO CONTROL ROOF DRAG, ENVELOP GW DRAG AND
75	C * HINES' DOPPLER SPREADING EXTROWAVE GW DRAG.
76	C * LOZPR IS TRUE FOR ZPR ENHANCEMENT
77
78
79	C * WORK ARRAYS.
80
81	REAL M_ALPHA(klon,KLEV,NAZMTH), V_ALPHA(klon,KLEV,NAZMTH),
82	1 SIGMA_ALPHA(klon,KLEV,NAZMTH),
83	1 SIGSQH_ALPHA(klon,KLEV,NAZMTH),
84	2 DRAG_U(klon,KLEV), DRAG_V(klon,KLEV), FLUX_U(klon,KLEV),
85	3 FLUX_V(klon,KLEV), HEAT(klon,KLEV), DIFFCO(klon,KLEV),
86	4 BVFREQ(klon,KLEV), DENSITY(klon,KLEV), SIGMA_T(klon,KLEV),
87	5 VISC_MOL(klon,KLEV), ALT(klon,KLEV),
88	6 SIGSQMCW(klon,KLEV,NAZMTH),
89	6 SIGMATM(klon,KLEV),
90	7 AK_ALPHA(klon,NAZMTH), K_ALPHA(klon,NAZMTH),
91	8 MMIN_ALPHA(klon,NAZMTH), I_ALPHA(klon,NAZMTH),
92	9 RMSWIND(klon), BVFBOT(klon), DENSBOT(klon)
93	REAL SMOOTHR1(klon,KLEV), SMOOTHR2(klon,KLEV)
94	REAL SIGALPMC(klon,KLEV,NAZMTH)
95	REAL F2MOD(klon,KLEV)
96
97	C * THES ARE THE INPUT PARAMETERS FOR HINES ROUTINE AND
98	C * ARE SPECIFIED IN ROUTINE HINES_SETUP. SINCE THIS IS CALLED
99	C * ONLY AT FIRST CALL TO THIS ROUTINE THESE VARIABLES MUST BE SAVED
100	C * FOR USE AT SUBSEQUENT CALLS. THIS CAN BE AVOIDED BY CALLING
101	C * HINES_SETUP IN MAIN PROGRAM AND PASSING THE PARAMETERS AS
102	C * SUBROUTINE ARGUEMENTS.
103	C
104
105	REAL RMSCON
106	INTEGER NMESSG, IPRINT, ILRMS
107	INTEGER IFL
108	C
109	INTEGER NAZ,ICUTOFF,NSMAX,IHEATCAL
110	REAL SLOPE,F1,F2,F3,F5,F6,KSTAR(KLON),ALT_CUTOFF,SMCO
111	C
112	C PROVIDED AS INPUT
113	C
114	integer nlon,nlev
115
116	real dtime
117	real paphm1x(nlon,nlev+1), papm1x(nlon,nlev)
118	real ux(nlon,nlev), vx(nlon,nlev), tx(nlon,nlev)
119	c
120	c VARIABLES FOR OUTPUT
121	c
122
123	real d_t_hin(nlon,nlev),d_u_hin(nlon,nlev),d_v_hin(nlon,nlev)
124	real zustrhi(nlon),zvstrhi(nlon)
125
126	C
127	C * LOGICAL SWITCHES TO CONTROL PRECIP ENHANCEMENT AND
128	C * HINES' DOPPLER SPREADING EXTROWAVE GW DRAG.
129	C * LOZPR IS TRUE FOR ZPR ENHANCEMENT
130	C
131	LOGICAL LOZPR, LORMS(klon)
132	C
133	C LOCAL PARAMETERS TO MAKE THINGS WORK (TEMPORARY VARIABLE)
134
135	REAL RHOH2O,ZPCONS,RGOCP,ZLAT,DTTDSF,RATIO,HSCAL
136	INTEGER I,J,L,JL,JK,LE,LREF,LREFP,LEVBOT
137	C
138	C DATA PARAMETERS NEEDED, EXPLAINED LATER
139
140	REAL V0,VMIN,DMPSCAL,TAUFAC,HMIN,APIBT,CPART,FCRIT
141	REAL PCRIT,PCONS
142	INTEGER IPLEV,IERROR
143
144	C
145	C
146	C PRINT *,' IT IS STARTED HINES GOING ON...'
147	C
148	C
149	C
150	C
151	C* COMPUTATIONAL CONSTANTS.
152	C ------------- ----------
153	C
154	C
155	d_t_hin(:,:)=0.
156
157	RHOH2O=1000.
158	ZPCONS = (1000.*86400.)/RHOH2O
159	cym KLEVM1=KLEV-1
160	C
161
162	do jl=kidia,kfdia
163	PAPHM1(JL,1) = paphm1x(JL,klev+1)
164	do jk=1,klev
165	le=klev+1-jk
166	PAPHM1(JL,JK+1) = paphm1x(JL,le)
167	PAPM1(JL,JK) = papm1x(JL,le)
168	PTM1(JL,JK) = tx(JL,le)
169	PUM1(JL,JK) = ux(JL,le)
170	PVM1(JL,JK) = vx(JL,le)
171	enddo
172	enddo
173	C
174	100 CONTINUE
175	C
176	C Define constants and arrays needed for the ccc/mam gwd scheme
177	C *Constants:
178
179	RGOCP=RD/RCPD
180	LREFP=KLEV-1
181	LREF=KLEV-2
182	C1
183	C1 *Arrays
184	C1
185	DO 2101 JK=1,KLEV
186	DO 2102 JL=KIDIA,KFDIA
187	SHJ(JL,JK)=PAPM1(JL,JK)/PAPHM1(JL,klev+1)
188	SGJ(JL,JK)=PAPM1(JL,JK)/PAPHM1(JL,klev+1)
189	DSGJ(JL,JK)=(PAPHM1(JL,JK+1)-PAPHM1(JL,JK))/PAPHM1(JL,klev+1)
190	SHXKJ(JL,JK)=(PAPM1(JL,JK)/PAPHM1(JL,klev+1))**RGOCP
191	TH(JL,JK)= PTM1(JL,JK)
192	2102 CONTINUE
193	2101 CONTINUE
194
195	CC
196	DO 211 JL=KIDIA,KFDIA
197	PRESSG(JL)=PAPHM1(JL,klev+1)
198	211 CONTINUE
199	C
200	C
201	DO 301 JL=KIDIA,KFDIA
202	PRFLUX(JL) = 0.0
203	ZPR(JL)=ZPCONS*PRFLUX(JL)
204	ZLAT=(RLAT(JL)/180.)*RPI
205	COSLAT(Jl)=COS(ZLAT)
206	301 CONTINUE
207	C
208	C
209	400 CONTINUE
210	C
211	C
212	C
213	C
214	*/#########################################################################
215	*/
216	*/
217	C
218	C * AUG. 14/95 - C. MCLANDRESS.
219	C * SEP. 95 N. MCFARLANE.
220	C
221	C * THIS ROUTINE CALCULATES THE HORIZONTAL WIND TENDENCIES
222	C * DUE TO MCFARLANE'S OROGRAPHIC GW DRAG SCHEME, HINES'
223	C * DOPPLER SPREAD SCHEME FOR "EXTROWAVES" AND ADDS ON
224	C * ROOF DRAG. IT IS BASED ON THE ROUTINE GWDFLX8.
225	C
226	C * LREFP IS THE INDEX OF THE MODEL LEVEL BELOW THE REFERENCE LEVEL
227	C * I/O ARRAYS PASSED FROM MAIN.
228	C * (PRESSG = SURFACE PRESSURE)
229	C
230	C
231	C
232	C
233	C * CONSTANTS VALUES DEFINED IN DATA STATEMENT ARE :
234	C * VMIN = MIMINUM WIND IN THE DIRECTION OF REFERENCE LEVEL
235	C * WIND BEFORE WE CONSIDER BREAKING TO HAVE OCCURED.
236	C * DMPSCAL = DAMPING TIME FOR GW DRAG IN SECONDS.
237	C * TAUFAC = 1/(LENGTH SCALE).
238	C * HMIN = MIMINUM ENVELOPE HEIGHT REQUIRED TO PRODUCE GW DRAG.
239	C * V0 = VALUE OF WIND THAT APPROXIMATES ZERO.
240
241
242	DATA VMIN / 5.0 /, V0 / 1.E-10 /,
243	1 TAUFAC/ 5.E-6 /, HMIN / 40000. /,
244	3 DMPSCAL / 6.5E+6 /, APIBT / 1.5708 /,
245	4 CPART / 0.7 /, FCRIT / 1. /
246
247	C * HINES EXTROWAVE GWD CONSTANTS DEFINED IN DATA STATEMENT ARE:
248	C * RMSCON = ROOT MEAN SQUARE GRAVITY WAVE WIND AT LOWEST LEVEL (M/S).
249	C * NMESSG = UNIT NUMBER FOR PRINTED MESSAGES.
250	C * IPRINT = 1 TO DO PRINT OUT SOME HINES ARRAYS.
251	C * IFL = FIRST CALL FLAG TO HINES_SETUP ("SAVE" IT)
252	C * PCRIT = CRITICAL VALUE OF ZPR (MM/D)
253	C * IPLEV = LEVEL OF APPLICATION OF PRCIT
254	C * PCONS = FACTOR OF ZPR ENHANCEMENT
255	C
256
257	DATA PCRIT / 5. /, PCONS / 4.75 /
258
259	IPLEV = LREFP-1
260	C
261	DATA RMSCON / 1.00 /
262	1 IPRINT / 2 /, NMESSG / 6 /
263	DATA IFL / 0 /
264	C
265	LOZPR = .FALSE.
266	C
267	C-----------------------------------------------------------------------
268	C
269	C
270	C
271	C * SET ERROR FLAG
272
273	IERROR = 0
274
275	C * SPECIFY VARIOUS PARAMETERS FOR HINES ROUTINE AT VERY FIRST CALL.
276	C * (NOTE THAT ARRAY K_ALPHA IS SPECIFIED SO MAKE SURE THAT
277	C * IT IS NOT OVERWRITTEN LATER ON).
278	C
279	CALL HINES_SETUP (NAZ,SLOPE,F1,F2,F3,F5,F6,KSTAR,
280	1 ICUTOFF,ALT_CUTOFF,SMCO,NSMAX,IHEATCAL,
281	2 K_ALPHA,IERROR,NMESSG,klon,NAZMTH,COSLAT)
282	IF (IERROR.NE.0) GO TO 999
283	C
284	C * START GWD CALCULATIONS.
285
286	LREF = LREFP-1
287
288	C
289	DO 105 J=1,NAZMTH
290	DO 105 L=1,KLEV
291	DO 105 I=kidia,klon
292	SIGSQMCW(I,L,J) = 0.
293	105 CONTINUE
294	c
295
296
297	C * INITIALIZE NECESSARY ARRAYS.
298	C
299	DO 120 L=1,KLEV
300	DO 120 I=KIDIA,KFDIA
301	UTENDGW(I,L) = 0.
302	VTENDGW(I,L) = 0.
303
304	UHS(I,L) = 0.
305	VHS(I,L) = 0.
306
307	120 CONTINUE
308	C
309	C * IF USING HINES SCHEME THEN CALCULATE B V FREQUENCY AT ALL POINTS
310	C * AND SMOOTH BVFREQ.
311
312	DO 130 L=2,KLEV
313	DO 130 I=KIDIA,KFDIA
314	DTTDSF=(TH(I,L)/SHXKJ(I,L)-TH(I,L-1)/
315	1 SHXKJ(I,L-1))/(SHJ(I,L)-SHJ(I,L-1))
316	DTTDSF=MIN(DTTDSF, -5./SGJ(I,L))
317	BVFREQ(I,L)=SQRT(-DTTDSFSGJ(I,L)(SGJ(I,L)**RGOCP)/RD)
318	1 *RG/PTM1(I,L)
319	130 CONTINUE
320	DO 135 L=1,KLEV
321	DO 135 I=KIDIA,KFDIA
322	IF(L.EQ.1) THEN
323	BVFREQ(I,L) = BVFREQ(I,L+1)
324	ENDIF
325	IF(L.GT.1) THEN
326	RATIO=5.*LOG(SGJ(I,L)/SGJ(I,L-1))
327	BVFREQ(I,L) = (BVFREQ(I,L-1) + RATIO*BVFREQ(I,L))
328	1 /(1.+RATIO)
329	ENDIF
330	135 CONTINUE
331	C
332	C
333	300 CONTINUE
334
335	C * CALCULATE GW DRAG DUE TO HINES' EXTROWAVES
336	C * SET MOLECULAR VISCOSITY TO A VERY SMALL VALUE.
337	C * IF THE MODEL TOP IS GREATER THAN 100 KM THEN THE ACTUAL
338	C * VISCOSITY COEFFICIENT COULD BE SPECIFIED HERE.
339
340	DO 310 L=1,KLEV
341	DO 310 I=KIDIA,KFDIA
342	VISC_MOL(I,L) = 1.5E-5
343	DRAG_U(I,L) = 0.
344	DRAG_V(I,L) = 0.
345	FLUX_U(I,L) = 0.
346	FLUX_V(I,L) = 0.
347	HEAT(I,L) = 0.
348	DIFFCO(I,L) = 0.
349	310 CONTINUE
350
351	C * ALTITUDE AND DENSITY AT BOTTOM.
352
353	DO 330 I=KIDIA,KFDIA
354	HSCAL = RD * PTM1(I,KLEV) / RG
355	DENSITY(I,KLEV) = SGJ(I,KLEV) * PRESSG(I) / (RG*HSCAL)
356	ALT(I,KLEV) = 0.
357	330 CONTINUE
358
359	C * ALTITUDE AND DENSITY AT REMAINING LEVELS.
360
361	DO 340 L=KLEV-1,1,-1
362	DO 340 I=KIDIA,KFDIA
363	HSCAL = RD * PTM1(I,L) / RG
364	ALT(I,L) = ALT(I,L+1) + HSCAL * DSGJ(I,L) / SGJ(I,L)
365	DENSITY(I,L) = SGJ(I,L) * PRESSG(I) / (RG*HSCAL)
366	340 CONTINUE
367
368	C
369	C * INITIALIZE SWITCHES FOR HINES GWD CALCULATION
370	C
371	ILRMS = 0
372	C
373	DO 345 I=KIDIA,KFDIA
374	LORMS(I) = .FALSE.
375	345 CONTINUE
376	C
377	C
378	C * DEFILE BOTTOM LAUNCH LEVEL
379	C
380	LEVBOT = IPLEV
381	C
382	C * BACKGROUND WIND MINUS VALUE AT BOTTOM LAUNCH LEVEL.
383	C
384	DO 351 L=1,LEVBOT
385	DO 351 I=KIDIA,KFDIA
386	UHS(I,L) = PUM1(I,L) - PUM1(I,LEVBOT)
387	VHS(I,L) = PVM1(I,L) - PVM1(I,LEVBOT)
388	351 CONTINUE
389	C
390	C * SPECIFY ROOT MEAN SQUARE WIND AT BOTTOM LAUNCH LEVEL.
391	C
392	DO 355 I=KIDIA,KFDIA
393	RMSWIND(I) = RMSCON
394	355 CONTINUE
395
396	IF (LOZPR) THEN
397	DO 350 I=KIDIA,KFDIA
398	IF (ZPR(I) .GT. PCRIT) THEN
399	RMSWIND(I) = RMSCON
400	> +( (ZPR(I)-PCRIT)/ZPR(I) )*PCONS
401	ENDIF
402	350 CONTINUE
403	ENDIF
404	C
405	DO 356 I=KIDIA,KFDIA
406	IF (RMSWIND(I) .GT. 0.0) THEN
407	ILRMS = ILRMS+1
408	LORMS(I) = .TRUE.
409	ENDIF
410	356 CONTINUE
411	C
412	C * CALCULATE GWD (NOTE THAT DIFFUSION COEFFICIENT AND
413	C * HEATING RATE ONLY CALCULATED IF IHEATCAL = 1).
414	C
415	IF ( ILRMS.GT.0 ) THEN
416	C
417	CALL HINES_EXTRO0 (DRAG_U,DRAG_V,HEAT,DIFFCO,FLUX_U,FLUX_V,
418	1 UHS,VHS,BVFREQ,DENSITY,VISC_MOL,ALT,
419	2 RMSWIND,K_ALPHA,M_ALPHA,V_ALPHA,
420	3 SIGMA_ALPHA,SIGSQH_ALPHA,AK_ALPHA,
421	4 MMIN_ALPHA,I_ALPHA,SIGMA_T,DENSBOT,BVFBOT,
422	5 1,IHEATCAL,ICUTOFF,IPRINT,NSMAX,
423	6 SMCO,ALT_CUTOFF,KSTAR,SLOPE,
424	7 F1,F2,F3,F5,F6,NAZ,SIGSQMCW,SIGMATM,
425	8 KIDIA,klon,1,LEVBOT,KLON,KLEV,NAZMTH,
426	9 LORMS,SMOOTHR1,SMOOTHR2,
427	9 SIGALPMC,F2MOD)
428
429	C * ADD ON HINES' GWD TENDENCIES TO OROGRAPHIC TENDENCIES AND
430	C * APPLY HINES' GW DRAG ON (UROW,VROW) WORK ARRAYS.
431
432	DO 360 L=1,KLEV
433	DO 360 I=KIDIA,KFDIA
434	UTENDGW(I,L) = UTENDGW(I,L) + DRAG_U(I,L)
435	VTENDGW(I,L) = VTENDGW(I,L) + DRAG_V(I,L)
436	360 CONTINUE
437	C
438
439	C * END OF HINES CALCULATIONS.
440	C
441	ENDIF
442	C
443	500 CONTINUE
444
445
446	C-----------------------------------------------------------------------
447	C
448	do jl=kidia,kfdia
449	zustrhi(jl)=FLUX_U(jl,1)
450	zvstrhi(jl)=FLUX_v(jl,1)
451	do jk=1,klev
452	le=klev-jk+1
453	d_u_hin(jl,JK) = UTENDGW(jl,le) * dtime
454	d_v_hin(jl,JK) = VTENDGW(jl,le) * dtime
455	enddo
456	enddo
457
458	c PRINT *,'UTENDGW:',UTENDGW
459
460	C PRINT *,' HINES HAS BEEN COMPLETED (LONG ISNT IT...)'
461
462	RETURN
463	999 CONTINUE
464
465	C * IF ERROR DETECTED THEN ABORT.
466
467	WRITE (NMESSG,6000)
468	WRITE (NMESSG,6010) IERROR
469	6000 FORMAT (/' EXECUTION ABORTED IN GWDOREXV')
470	6010 FORMAT (' ERROR FLAG =',I4)
471
472	C
473	RETURN
474	END
475	*/
476	*/
477
478
479	SUBROUTINE HINES_EXTRO0 (DRAG_U,DRAG_V,HEAT,DIFFCO,FLUX_U,FLUX_V,
480	1 VEL_U,VEL_V,BVFREQ,DENSITY,VISC_MOL,ALT,
481	2 RMSWIND,K_ALPHA,M_ALPHA,V_ALPHA,
482	3 SIGMA_ALPHA,SIGSQH_ALPHA,AK_ALPHA,
483	4 MMIN_ALPHA,I_ALPHA,SIGMA_T,DENSB,BVFB,
484	5 IORDER,IHEATCAL,ICUTOFF,IPRINT,NSMAX,
485	6 SMCO,ALT_CUTOFF,KSTAR,SLOPE,
486	7 F1,F2,F3,F5,F6,NAZ,SIGSQMCW,SIGMATM,
487	8 IL1,IL2,LEV1,LEV2,NLONS,NLEVS,NAZMTH,
488	9 LORMS,SMOOTHR1,SMOOTHR2,
489	9 SIGALPMC,F2MOD)
490
491	implicit none
492	C
493	C Main routine for Hines' "extrowave" gravity wave parameterization based
494	C on Hines' Doppler spread theory. This routine calculates zonal
495	C and meridional components of gravity wave drag, heating rates
496	C and diffusion coefficient on a longitude by altitude grid.
497	C No "mythical" lower boundary region calculation is made so it
498	C is assumed that lowest level winds are weak (i.e, approximately zero).
499	C
500	C Aug. 13/95 - C. McLandress
501	C SEPT. /95 - N.McFarlane
502	C
503	C Modifications:
504	C
505	C Output arguements:
506	C
507	C * DRAG_U = zonal component of gravity wave drag (m/s^2).
508	C * DRAG_V = meridional component of gravity wave drag (m/s^2).
509	C * HEAT = gravity wave heating (K/sec).
510	C * DIFFCO = diffusion coefficient (m^2/sec)
511	C * FLUX_U = zonal component of vertical momentum flux (Pascals)
512	C * FLUX_V = meridional component of vertical momentum flux (Pascals)
513	C
514	C Input arguements:
515	C
516	C * VEL_U = background zonal wind component (m/s).
517	C * VEL_V = background meridional wind component (m/s).
518	C * BVFREQ = background Brunt Vassala frequency (radians/sec).
519	C * DENSITY = background density (kg/m^3)
520	C * VISC_MOL = molecular viscosity (m^2/s)
521	C * ALT = altitude of momentum, density, buoyancy levels (m)
522	C * (NOTE: levels ordered so that ALT(I,1) > ALT(I,2), etc.)
523	C * RMSWIND = root mean square gravity wave wind at lowest level (m/s).
524	C * K_ALPHA = horizontal wavenumber of each azimuth (1/m).
525	C * IORDER = 1 means vertical levels are indexed from top down
526	C * (i.e., highest level indexed 1 and lowest level NLEVS);
527	C * .NE. 1 highest level is index NLEVS.
528	C * IHEATCAL = 1 to calculate heating rates and diffusion coefficient.
529	C * IPRINT = 1 to print out various arrays.
530	C * ICUTOFF = 1 to exponentially damp GWD, heating and diffusion
531	C * arrays above ALT_CUTOFF; otherwise arrays not modified.
532	C * ALT_CUTOFF = altitude in meters above which exponential decay applied.
533	C * SMCO = smoothing factor used to smooth cutoff vertical
534	C * wavenumbers and total rms winds in vertical direction
535	C * before calculating drag or heating
536	C * (SMCO >= 1 ==> 1:SMCO:1 stencil used).
537	C * NSMAX = number of times smoother applied ( >= 1),
538	C * = 0 means no smoothing performed.
539	C * KSTAR = typical gravity wave horizontal wavenumber (1/m).
540	C * SLOPE = slope of incident vertical wavenumber spectrum
541	C * (SLOPE must equal 1., 1.5 or 2.).
542	C * F1 to F6 = Hines's fudge factors (F4 not needed since used for
543	C * vertical flux of vertical momentum).
544	C * NAZ = actual number of horizontal azimuths used.
545	C * IL1 = first longitudinal index to use (IL1 >= 1).
546	C * IL2 = last longitudinal index to use (IL1 <= IL2 <= NLONS).
547	C * LEV1 = index of first level for drag calculation.
548	C * LEV2 = index of last level for drag calculation
549	C * (i.e., LEV1 < LEV2 <= NLEVS).
550	C * NLONS = number of longitudes.
551	C * NLEVS = number of vertical levels.
552	C * NAZMTH = azimuthal array dimension (NAZMTH >= NAZ).
553	C
554	C Work arrays.
555	C
556	C * M_ALPHA = cutoff vertical wavenumber (1/m).
557	C * V_ALPHA = wind component at each azimuth (m/s) and if IHEATCAL=1
558	C * holds vertical derivative of cutoff wavenumber.
559	C * SIGMA_ALPHA = total rms wind in each azimuth (m/s).
560	C * SIGSQH_ALPHA = portion of wind variance from waves having wave
561	C * normals in the alpha azimuth (m/s).
562	C * SIGMA_T = total rms horizontal wind (m/s).
563	C * AK_ALPHA = spectral amplitude factor at each azimuth
564	C * (i.e.,{AjKj}) in m^4/s^2.
565	C * I_ALPHA = Hines' integral.
566	C * MMIN_ALPHA = minimum value of cutoff wavenumber.
567	C * DENSB = background density at bottom level.
568	C * BVFB = buoyancy frequency at bottom level and
569	C * work array for ICUTOFF = 1.
570	C
571	C * LORMS = .TRUE. for drag computation
572	C
573	INTEGER NAZ, NLONS, NLEVS, NAZMTH, IL1, IL2, LEV1, LEV2
574	INTEGER ICUTOFF, NSMAX, IORDER, IHEATCAL, IPRINT
575	REAL KSTAR(NLONS), F1, F2, F3, F5, F6, SLOPE
576	REAL ALT_CUTOFF, SMCO
577	REAL DRAG_U(NLONS,NLEVS), DRAG_V(NLONS,NLEVS)
578	REAL HEAT(NLONS,NLEVS), DIFFCO(NLONS,NLEVS)
579	REAL FLUX_U(NLONS,NLEVS), FLUX_V(NLONS,NLEVS)
580	REAL VEL_U(NLONS,NLEVS), VEL_V(NLONS,NLEVS)
581	REAL BVFREQ(NLONS,NLEVS), DENSITY(NLONS,NLEVS)
582	REAL VISC_MOL(NLONS,NLEVS), ALT(NLONS,NLEVS)
583	REAL RMSWIND(NLONS), BVFB(NLONS), DENSB(NLONS)
584	REAL SIGMA_T(NLONS,NLEVS), SIGSQMCW(NLONS,NLEVS,NAZMTH)
585	REAL SIGMA_ALPHA(NLONS,NLEVS,NAZMTH), SIGMATM(NLONS,NLEVS)
586	REAL SIGSQH_ALPHA(NLONS,NLEVS,NAZMTH)
587	REAL M_ALPHA(NLONS,NLEVS,NAZMTH), V_ALPHA(NLONS,NLEVS,NAZMTH)
588	REAL AK_ALPHA(NLONS,NAZMTH), K_ALPHA(NLONS,NAZMTH)
589	REAL MMIN_ALPHA(NLONS,NAZMTH) , I_ALPHA(NLONS,NAZMTH)
590	REAL SMOOTHR1(NLONS,NLEVS), SMOOTHR2(NLONS,NLEVS)
591	REAL SIGALPMC(NLONS,NLEVS,NAZMTH)
592	REAL F2MOD(NLONS,NLEVS)
593	C
594	LOGICAL LORMS(NLONS)
595	C
596	C Internal variables.
597	C
598	INTEGER LEVBOT, LEVTOP, I, N, L, LEV1P, LEV2M
599	INTEGER ILPRT1, ILPRT2
600	C-----------------------------------------------------------------------
601	C
602	C PRINT *,' IN HINES_EXTRO0'
603	LEV1P = LEV1 + 1
604	LEV2M = LEV2 - 1
605	C
606	C Index of lowest altitude level (bottom of drag calculation).
607	C
608	LEVBOT = LEV2
609	LEVTOP = LEV1
610	IF (IORDER.NE.1) THEN
611	write(6,1)
612	1 format(2x,' error: IORDER NOT ONE! ')
613	END IF
614	C
615	C Buoyancy and density at bottom level.
616	C
617	DO 10 I = IL1,IL2
618	BVFB(I) = BVFREQ(I,LEVBOT)
619	DENSB(I) = DENSITY(I,LEVBOT)
620	10 CONTINUE
621	C
622	C initialize some variables
623	C
624	DO 20 N = 1,NAZ
625	DO 20 L=LEV1,LEV2
626	DO 20 I=IL1,IL2
627	M_ALPHA(I,L,N) = 0.0
628	20 CONTINUE
629	DO 21 L=LEV1,LEV2
630	DO 21 I=IL1,IL2
631	SIGMA_T(I,L) = 0.0
632	21 CONTINUE
633	DO 22 N = 1,NAZ
634	DO 22 I=IL1,IL2
635	I_ALPHA(I,N) = 0.0
636	22 CONTINUE
637	C
638	C Compute azimuthal wind components from zonal and meridional winds.
639	C
640	CALL HINES_WIND ( V_ALPHA,
641	^ VEL_U, VEL_V, NAZ,
642	^ IL1, IL2, LEV1, LEV2, NLONS, NLEVS, NAZMTH )
643	C
644	C Calculate cutoff vertical wavenumber and velocity variances.
645	C
646	CALL HINES_WAVNUM ( M_ALPHA, SIGMA_ALPHA, SIGSQH_ALPHA, SIGMA_T,
647	^ AK_ALPHA, V_ALPHA, VISC_MOL, DENSITY, DENSB,
648	^ BVFREQ, BVFB, RMSWIND, I_ALPHA, MMIN_ALPHA,
649	^ KSTAR, SLOPE, F1, F2, F3, NAZ, LEVBOT,
650	^ LEVTOP,IL1,IL2,NLONS,NLEVS,NAZMTH, SIGSQMCW,
651	^ SIGMATM,LORMS,SIGALPMC,F2MOD)
652	C
653	C Smooth cutoff wavenumbers and total rms velocity in the vertical
654	C direction NSMAX times, using FLUX_U as temporary work array.
655	C
656	IF (NSMAX.GT.0) THEN
657	DO 80 N = 1,NAZ
658	DO 81 L=LEV1,LEV2
659	DO 81 I=IL1,IL2
660	SMOOTHR1(I,L) = M_ALPHA(I,L,N)
661	81 CONTINUE
662	CALL VERT_SMOOTH (SMOOTHR1,
663	^ SMOOTHR2, SMCO, NSMAX,
664	^ IL1, IL2, LEV1, LEV2, NLONS, NLEVS )
665	DO 83 L=LEV1,LEV2
666	DO 83 I=IL1,IL2
667	M_ALPHA(I,L,N) = SMOOTHR1(I,L)
668	83 CONTINUE
669	80 CONTINUE
670	CALL VERT_SMOOTH ( SIGMA_T,
671	^ SMOOTHR2, SMCO, NSMAX,
672	^ IL1, IL2, LEV1, LEV2, NLONS, NLEVS )
673	END IF
674	C
675	C Calculate zonal and meridional components of the
676	C momentum flux and drag.
677	C
678	CALL HINES_FLUX ( FLUX_U, FLUX_V, DRAG_U, DRAG_V,
679	^ ALT, DENSITY, DENSB, M_ALPHA,
680	^ AK_ALPHA, K_ALPHA, SLOPE, NAZ,
681	^ IL1, IL2, LEV1, LEV2, NLONS, NLEVS, NAZMTH,
682	^ LORMS)
683	C
684	C Cutoff drag above ALT_CUTOFF, using BVFB as temporary work array.
685	C
686	IF (ICUTOFF.EQ.1) THEN
687	CALL HINES_EXP ( DRAG_U,
688	^ BVFB, ALT, ALT_CUTOFF, IORDER,
689	^ IL1, IL2, LEV1, LEV2, NLONS, NLEVS )
690	CALL HINES_EXP ( DRAG_V,
691	^ BVFB, ALT, ALT_CUTOFF, IORDER,
692	^ IL1, IL2, LEV1, LEV2, NLONS, NLEVS )
693	END IF
694	C
695	C Print out various arrays for diagnostic purposes.
696	C
697	IF (IPRINT.EQ.1) THEN
698	ILPRT1 = 15
699	ILPRT2 = 16
700	CALL HINES_PRINT ( FLUX_U, FLUX_V, DRAG_U, DRAG_V, ALT,
701	^ SIGMA_T, SIGMA_ALPHA, V_ALPHA, M_ALPHA,
702	^ 1, 1, 6, ILPRT1, ILPRT2, LEV1, LEV2,
703	^ NAZ, NLONS, NLEVS, NAZMTH)
704	END IF
705	C
706	C If not calculating heating rate and diffusion coefficient then finished.
707	C
708	IF (IHEATCAL.NE.1) RETURN
709	C
710	C Calculate vertical derivative of cutoff wavenumber (store
711	C in array V_ALPHA) using centered differences at interior gridpoints
712	C and one-sided differences at first and last levels.
713	C
714	DO 130 N = 1,NAZ
715	DO 100 L = LEV1P,LEV2M
716	DO 90 I = IL1,IL2
717	V_ALPHA(I,L,N) = ( M_ALPHA(I,L+1,N) - M_ALPHA(I,L-1,N) )
718	^ / ( ALT(I,L+1) - ALT(I,L-1) )
719	90 CONTINUE
720	100 CONTINUE
721	DO 110 I = IL1,IL2
722	V_ALPHA(I,LEV1,N) = ( M_ALPHA(I,LEV1P,N) - M_ALPHA(I,LEV1,N) )
723	^ / ( ALT(I,LEV1P) - ALT(I,LEV1) )
724	110 CONTINUE
725	DO 120 I = IL1,IL2
726	V_ALPHA(I,LEV2,N) = ( M_ALPHA(I,LEV2,N) - M_ALPHA(I,LEV2M,N) )
727	^ / ( ALT(I,LEV2) - ALT(I,LEV2M) )
728	120 CONTINUE
729	130 CONTINUE
730	C
731	C Heating rate and diffusion coefficient.
732	C
733	CALL HINES_HEAT ( HEAT, DIFFCO,
734	^ M_ALPHA, V_ALPHA, AK_ALPHA, K_ALPHA,
735	^ BVFREQ, DENSITY, DENSB, SIGMA_T, VISC_MOL,
736	^ KSTAR, SLOPE, F2, F3, F5, F6, NAZ,
737	^ IL1, IL2, LEV1, LEV2, NLONS, NLEVS, NAZMTH)
738	C
739	C Finished.
740	C
741	RETURN
742	C-----------------------------------------------------------------------
743	END
744
745	SUBROUTINE HINES_WAVNUM (M_ALPHA,SIGMA_ALPHA,SIGSQH_ALPHA,SIGMA_T,
746	1 AK_ALPHA,V_ALPHA,VISC_MOL,DENSITY,DENSB,
747	2 BVFREQ,BVFB,RMS_WIND,I_ALPHA,MMIN_ALPHA,
748	3 KSTAR,SLOPE,F1,F2,F3,NAZ,LEVBOT,LEVTOP,
749	4 IL1,IL2,NLONS,NLEVS,NAZMTH,SIGSQMCW,
750	5 SIGMATM,LORMS,SIGALPMC,F2MOD)
751	C
752	C This routine calculates the cutoff vertical wavenumber and velocity
753	C variances on a longitude by altitude grid for the Hines' Doppler
754	C spread gravity wave drag parameterization scheme.
755	C NOTE: (1) only values of four or eight can be used for # azimuths (NAZ).
756	C (2) only values of 1.0, 1.5 or 2.0 can be used for slope (SLOPE).
757	C
758	C Aug. 10/95 - C. McLandress
759	C
760	C Output arguements:
761	C
762	C * M_ALPHA = cutoff wavenumber at each azimuth (1/m).
763	C * SIGMA_ALPHA = total rms wind in each azimuth (m/s).
764	C * SIGSQH_ALPHA = portion of wind variance from waves having wave
765	C * normals in the alpha azimuth (m/s).
766	C * SIGMA_T = total rms horizontal wind (m/s).
767	C * AK_ALPHA = spectral amplitude factor at each azimuth
768	C * (i.e.,{AjKj}) in m^4/s^2.
769	C
770	C Input arguements:
771	C
772	C * V_ALPHA = wind component at each azimuth (m/s).
773	C * VISC_MOL = molecular viscosity (m^2/s)
774	C * DENSITY = background density (kg/m^3).
775	C * DENSB = background density at model bottom (kg/m^3).
776	C * BVFREQ = background Brunt Vassala frequency (radians/sec).
777	C * BVFB = background Brunt Vassala frequency at model bottom.
778	C * RMS_WIND = root mean square gravity wave wind at lowest level (m/s).
779	C * KSTAR = typical gravity wave horizontal wavenumber (1/m).
780	C * SLOPE = slope of incident vertical wavenumber spectrum
781	C * (SLOPE = 1., 1.5 or 2.).
782	C * F1,F2,F3 = Hines's fudge factors.
783	C * NAZ = actual number of horizontal azimuths used (4 or 8).
784	C * LEVBOT = index of lowest vertical level.
785	C * LEVTOP = index of highest vertical level
786	C * (NOTE: if LEVTOP < LEVBOT then level index
787	C * increases from top down).
788	C * IL1 = first longitudinal index to use (IL1 >= 1).
789	C * IL2 = last longitudinal index to use (IL1 <= IL2 <= NLONS).
790	C * NLONS = number of longitudes.
791	C * NLEVS = number of vertical levels.
792	C * NAZMTH = azimuthal array dimension (NAZMTH >= NAZ).
793	C
794	C * LORMS = .TRUE. for drag computation
795	C
796	C Input work arrays:
797	C
798	C * I_ALPHA = Hines' integral at a single level.
799	C * MMIN_ALPHA = minimum value of cutoff wavenumber.
800	C
801	INTEGER NAZ, LEVBOT, LEVTOP, IL1, IL2, NLONS, NLEVS, NAZMTH
802	REAL SLOPE, KSTAR(NLONS), F1, F2, F3
803	REAL M_ALPHA(NLONS,NLEVS,NAZMTH)
804	REAL SIGMA_ALPHA(NLONS,NLEVS,NAZMTH)
805	REAL SIGALPMC(NLONS,NLEVS,NAZMTH)
806	REAL SIGSQH_ALPHA(NLONS,NLEVS,NAZMTH)
807	REAL SIGSQMCW(NLONS,NLEVS,NAZMTH)
808	REAL SIGMA_T(NLONS,NLEVS)
809	REAL SIGMATM(NLONS,NLEVS)
810	REAL AK_ALPHA(NLONS,NAZMTH)
811	REAL V_ALPHA(NLONS,NLEVS,NAZMTH)
812	REAL VISC_MOL(NLONS,NLEVS)
813	REAL F2MOD(NLONS,NLEVS)
814	REAL DENSITY(NLONS,NLEVS), DENSB(NLONS)
815	REAL BVFREQ(NLONS,NLEVS), BVFB(NLONS), RMS_WIND(NLONS)
816	REAL I_ALPHA(NLONS,NAZMTH), MMIN_ALPHA(NLONS,NAZMTH)
817	C
818	LOGICAL LORMS(NLONS)
819	C
820	C Internal variables.
821	C
822	INTEGER I, L, N, LSTART, LEND, LINCR, LBELOW
823	REAL M_SUB_M_TURB, M_SUB_M_MOL, M_TRIAL
824	REAL VISC, VISC_MIN, AZFAC, SP1
825
826	cc REAL N_OVER_M(1000), SIGFAC(1000)
827
828	REAL N_OVER_M(NLONS), SIGFAC(NLONS)
829	DATA VISC_MIN / 1.E-10 /
830	C-----------------------------------------------------------------------
831	C
832
833	C PRINT *,'IN HINES_WAVNUM'
834	SP1 = SLOPE + 1.
835	C
836	C Indices of levels to process.
837	C
838	IF (LEVBOT.GT.LEVTOP) THEN
839	LSTART = LEVBOT - 1
840	LEND = LEVTOP
841	LINCR = -1
842	ELSE
843	write(6,1)
844	1 format(2x,' error: IORDER NOT ONE! ')
845	END IF
846	C
847	C Use horizontal isotropy to calculate azimuthal variances at bottom level.
848	C
849	AZFAC = 1. / REAL(NAZ)
850	DO 20 N = 1,NAZ
851	DO 10 I = IL1,IL2
852	SIGSQH_ALPHA(I,LEVBOT,N) = AZFAC * RMS_WIND(I)**2
853	10 CONTINUE
854	20 CONTINUE
855	C
856	C Velocity variances at bottom level.
857	C
858	CALL HINES_SIGMA ( SIGMA_T, SIGMA_ALPHA,
859	^ SIGSQH_ALPHA, NAZ, LEVBOT,
860	^ IL1, IL2, NLONS, NLEVS, NAZMTH)
861	c
862	CALL HINES_SIGMA ( SIGMATM, SIGALPMC,
863	^ SIGSQMCW, NAZ, LEVBOT,
864	^ IL1, IL2, NLONS, NLEVS, NAZMTH)
865	C
866	C Calculate cutoff wavenumber and spectral amplitude factor
867	C at bottom level where it is assumed that background winds vanish
868	C and also initialize minimum value of cutoff wavnumber.
869	C
870	DO 40 N = 1,NAZ
871	DO 30 I = IL1,IL2
872	IF (LORMS(I)) THEN
873	M_ALPHA(I,LEVBOT,N) = BVFB(I) /
874	^ ( F1 * SIGMA_ALPHA(I,LEVBOT,N)
875	^ + F2 * SIGMA_T(I,LEVBOT) )
876	AK_ALPHA(I,N) = SIGSQH_ALPHA(I,LEVBOT,N)
877	^ / ( M_ALPHA(I,LEVBOT,N)**SP1 / SP1 )
878	MMIN_ALPHA(I,N) = M_ALPHA(I,LEVBOT,N)
879	ENDIF
880	30 CONTINUE
881	40 CONTINUE
882	C
883	C Calculate quantities from the bottom upwards,
884	C starting one level above bottom.
885	C
886	DO 150 L = LSTART,LEND,LINCR
887	C
888	C Level beneath present level.
889	C
890	LBELOW = L - LINCR
891	C
892	C Calculate N/m_M where m_M is maximum permissible value of the vertical
893	C wavenumber (i.e., m > m_M are obliterated) and N is buoyancy frequency.
894	C m_M is taken as the smaller of the instability-induced
895	C wavenumber (M_SUB_M_TURB) and that imposed by molecular viscosity
896	C (M_SUB_M_MOL). Since variance at this level is not yet known
897	C use value at level below.
898	C
899	DO 50 I = IL1,IL2
900	IF (LORMS(I)) THEN
901	c
902	F2MFAC=SIGMATM(I,LBELOW)**2
903	F2MOD(I,LBELOW) =1.+ 2.*F2MFAC
904	^ / ( F2MFAC+SIGMA_T(I,LBELOW)**2 )
905	c
906	VISC = AMAX1 ( VISC_MOL(I,L), VISC_MIN )
907	M_SUB_M_TURB = BVFREQ(I,L)
908	^ / ( F2 F2MOD(I,LBELOW)SIGMA_T(I,LBELOW))
909	M_SUB_M_MOL = (BVFREQ(I,L)KSTAR(I)/VISC)*0.33333333/F3
910	IF (M_SUB_M_TURB .LT. M_SUB_M_MOL) THEN
911	N_OVER_M(I) = F2 F2MOD(I,LBELOW)SIGMA_T(I,LBELOW)
912	ELSE
913	N_OVER_M(I) = BVFREQ(I,L) / M_SUB_M_MOL
914	END IF
915	ENDIF
916	50 CONTINUE
917	C
918	C Calculate cutoff wavenumber at this level.
919	C
920	DO 70 N = 1,NAZ
921	DO 60 I = IL1,IL2
922	IF (LORMS(I)) THEN
923	C
924	C Calculate trial value (since variance at this level is not yet known
925	C use value at level below). If trial value is negative or if it exceeds
926	C minimum value (not permitted) then set it to minimum value.
927	C
928	M_TRIAL = BVFB(I) / ( F1 * ( SIGMA_ALPHA(I,LBELOW,N)+
929	^ SIGALPMC(I,LBELOW,N)) + N_OVER_M(I) + V_ALPHA(I,L,N) )
930	IF (M_TRIAL.LE.0. .OR. M_TRIAL.GT.MMIN_ALPHA(I,N)) THEN
931	M_TRIAL = MMIN_ALPHA(I,N)
932	END IF
933	M_ALPHA(I,L,N) = M_TRIAL
934	C
935	C Reset minimum value of cutoff wavenumber if necessary.
936	C
937	IF (M_ALPHA(I,L,N) .LT. MMIN_ALPHA(I,N)) THEN
938	MMIN_ALPHA(I,N) = M_ALPHA(I,L,N)
939	END IF
940	C
941	ENDIF
942	60 CONTINUE
943	70 CONTINUE
944	C
945	C Calculate the Hines integral at this level.
946	C
947	CALL HINES_INTGRL ( I_ALPHA,
948	^ V_ALPHA, M_ALPHA, BVFB, SLOPE, NAZ,
949	^ L, IL1, IL2, NLONS, NLEVS, NAZMTH,
950	^ LORMS )
951
952	C
953	C Calculate the velocity variances at this level.
954	C
955	DO 80 I = IL1,IL2
956	SIGFAC(I) = DENSB(I) / DENSITY(I,L)
957	^ * BVFREQ(I,L) / BVFB(I)
958	80 CONTINUE
959	DO 100 N = 1,NAZ
960	DO 90 I = IL1,IL2
961	SIGSQH_ALPHA(I,L,N) = SIGFAC(I) * AK_ALPHA(I,N)
962	^ * I_ALPHA(I,N)
963	90 CONTINUE
964	100 CONTINUE
965	CALL HINES_SIGMA ( SIGMA_T, SIGMA_ALPHA,
966	^ SIGSQH_ALPHA, NAZ, L,
967	^ IL1, IL2, NLONS, NLEVS, NAZMTH )
968	c
969	CALL HINES_SIGMA ( SIGMATM, SIGALPMC,
970	^ SIGSQMCW, NAZ, L,
971	^ IL1, IL2, NLONS, NLEVS, NAZMTH )
972	C
973	C End of level loop.
974	C
975	150 CONTINUE
976	C
977	RETURN
978	C-----------------------------------------------------------------------
979	END
980
981	SUBROUTINE HINES_WIND (V_ALPHA,VEL_U,VEL_V,
982	1 NAZ,IL1,IL2,LEV1,LEV2,NLONS,NLEVS,NAZMTH)
983	C
984	C This routine calculates the azimuthal horizontal background wind components
985	C on a longitude by altitude grid for the case of 4 or 8 azimuths for
986	C the Hines' Doppler spread GWD parameterization scheme.
987	C
988	C Aug. 7/95 - C. McLandress
989	C
990	C Output arguement:
991	C
992	C * V_ALPHA = background wind component at each azimuth (m/s).
993	C * (note: first azimuth is in eastward direction
994	C * and rotate in counterclockwise direction.)
995	C
996	C Input arguements:
997	C
998	C * VEL_U = background zonal wind component (m/s).
999	C * VEL_V = background meridional wind component (m/s).
1000	C * NAZ = actual number of horizontal azimuths used (must be 4 or 8).
1001	C * IL1 = first longitudinal index to use (IL1 >= 1).
1002	C * IL2 = last longitudinal index to use (IL1 <= IL2 <= NLONS).
1003	C * LEV1 = first altitude level to use (LEV1 >=1).
1004	C * LEV2 = last altitude level to use (LEV1 < LEV2 <= NLEVS).
1005	C * NLONS = number of longitudes.
1006	C * NLEVS = number of vertical levels.
1007	C * NAZMTH = azimuthal array dimension (NAZMTH >= NAZ).
1008	C
1009	C Constants in DATA statements.
1010	C
1011	C * COS45 = cosine of 45 degrees.
1012	C * UMIN = minimum allowable value for zonal or meridional
1013	C * wind component (m/s).
1014	C
1015	C Subroutine arguements.
1016	C
1017	INTEGER NAZ, IL1, IL2, LEV1, LEV2
1018	INTEGER NLONS, NLEVS, NAZMTH
1019	REAL V_ALPHA(NLONS,NLEVS,NAZMTH)
1020	REAL VEL_U(NLONS,NLEVS), VEL_V(NLONS,NLEVS)
1021	C
1022	C Internal variables.
1023	C
1024	INTEGER I, L
1025	REAL U, V, COS45, UMIN
1026	C
1027	DATA COS45 / 0.7071068 /
1028	DATA UMIN / 0.001 /
1029	C-----------------------------------------------------------------------
1030	C
1031	C Case with 4 azimuths.
1032	C
1033
1034	C PRINT *,'IN HINES_WIND'
1035	IF (NAZ.EQ.4) THEN
1036	DO 20 L = LEV1,LEV2
1037	DO 10 I = IL1,IL2
1038	U = VEL_U(I,L)
1039	V = VEL_V(I,L)
1040	IF (ABS(U) .LT. UMIN) U = UMIN
1041	IF (ABS(V) .LT. UMIN) V = UMIN
1042	V_ALPHA(I,L,1) = U
1043	V_ALPHA(I,L,2) = V
1044	V_ALPHA(I,L,3) = - U
1045	V_ALPHA(I,L,4) = - V
1046	10 CONTINUE
1047	20 CONTINUE
1048	END IF
1049	C
1050	C Case with 8 azimuths.
1051	C
1052	IF (NAZ.EQ.8) THEN
1053	DO 40 L = LEV1,LEV2
1054	DO 30 I = IL1,IL2
1055	U = VEL_U(I,L)
1056	V = VEL_V(I,L)
1057	IF (ABS(U) .LT. UMIN) U = UMIN
1058	IF (ABS(V) .LT. UMIN) V = UMIN
1059	V_ALPHA(I,L,1) = U
1060	V_ALPHA(I,L,2) = COS45 * ( V + U )
1061	V_ALPHA(I,L,3) = V
1062	V_ALPHA(I,L,4) = COS45 * ( V - U )
1063	V_ALPHA(I,L,5) = - U
1064	V_ALPHA(I,L,6) = - V_ALPHA(I,L,2)
1065	V_ALPHA(I,L,7) = - V
1066	V_ALPHA(I,L,8) = - V_ALPHA(I,L,4)
1067	30 CONTINUE
1068	40 CONTINUE
1069	END IF
1070	C
1071	RETURN
1072	C-----------------------------------------------------------------------
1073	END
1074
1075	SUBROUTINE HINES_FLUX (FLUX_U,FLUX_V,DRAG_U,DRAG_V,ALT,DENSITY,
1076	1 DENSB,M_ALPHA,AK_ALPHA,K_ALPHA,SLOPE,
1077	2 NAZ,IL1,IL2,LEV1,LEV2,NLONS,NLEVS,NAZMTH,
1078	3 LORMS)
1079	C
1080	C Calculate zonal and meridional components of the vertical flux
1081	C of horizontal momentum and corresponding wave drag (force per unit mass)
1082	C on a longitude by altitude grid for the Hines' Doppler spread
1083	C GWD parameterization scheme.
1084	C NOTE: only 4 or 8 azimuths can be used.
1085	C
1086	C Aug. 6/95 - C. McLandress
1087	C
1088	C Output arguements:
1089	C
1090	C * FLUX_U = zonal component of vertical momentum flux (Pascals)
1091	C * FLUX_V = meridional component of vertical momentum flux (Pascals)
1092	C * DRAG_U = zonal component of drag (m/s^2).
1093	C * DRAG_V = meridional component of drag (m/s^2).
1094	C
1095	C Input arguements:
1096	C
1097	C * ALT = altitudes (m).
1098	C * DENSITY = background density (kg/m^3).
1099	C * DENSB = background density at bottom level (kg/m^3).
1100	C * M_ALPHA = cutoff vertical wavenumber (1/m).
1101	C * AK_ALPHA = spectral amplitude factor (i.e., {AjKj} in m^4/s^2).
1102	C * K_ALPHA = horizontal wavenumber (1/m).
1103	C * SLOPE = slope of incident vertical wavenumber spectrum.
1104	C * NAZ = actual number of horizontal azimuths used (must be 4 or 8).
1105	C * IL1 = first longitudinal index to use (IL1 >= 1).
1106	C * IL2 = last longitudinal index to use (IL1 <= IL2 <= NLONS).
1107	C * LEV1 = first altitude level to use (LEV1 >=1).
1108	C * LEV2 = last altitude level to use (LEV1 < LEV2 <= NLEVS).
1109	C * NLONS = number of longitudes.
1110	C * NLEVS = number of vertical levels.
1111	C * NAZMTH = azimuthal array dimension (NAZMTH >= NAZ).
1112	C
1113	C * LORMS = .TRUE. for drag computation
1114	C
1115	C Constant in DATA statement.
1116	C
1117	C * COS45 = cosine of 45 degrees.
1118	C
1119	C Subroutine arguements.
1120	C
1121	INTEGER NAZ, IL1, IL2, LEV1, LEV2
1122	INTEGER NLONS, NLEVS, NAZMTH
1123	REAL SLOPE
1124	REAL FLUX_U(NLONS,NLEVS), FLUX_V(NLONS,NLEVS)
1125	REAL DRAG_U(NLONS,NLEVS), DRAG_V(NLONS,NLEVS)
1126	REAL ALT(NLONS,NLEVS), DENSITY(NLONS,NLEVS), DENSB(NLONS)
1127	REAL M_ALPHA(NLONS,NLEVS,NAZMTH)
1128	REAL AK_ALPHA(NLONS,NAZMTH), K_ALPHA(NLONS,NAZMTH)
1129	C
1130	LOGICAL LORMS(NLONS)
1131	C
1132	C Internal variables.
1133	C
1134	INTEGER I, L, LEV1P, LEV2M
1135	REAL COS45, PROD2, PROD4, PROD6, PROD8, DENDZ, DENDZ2
1136	DATA COS45 / 0.7071068 /
1137	C-----------------------------------------------------------------------
1138	C
1139	LEV1P = LEV1 + 1
1140	LEV2M = LEV2 - 1
1141	LEV2P = LEV2 + 1
1142	C
1143	C Sum over azimuths for case where SLOPE = 1.
1144	C
1145	IF (SLOPE.EQ.1.) THEN
1146	C
1147	C Case with 4 azimuths.
1148	C
1149	IF (NAZ.EQ.4) THEN
1150	DO 20 L = LEV1,LEV2
1151	DO 10 I = IL1,IL2
1152	FLUX_U(I,L) = AK_ALPHA(I,1)K_ALPHA(I,1)M_ALPHA(I,L,1)
1153	^ - AK_ALPHA(I,3)K_ALPHA(I,3)M_ALPHA(I,L,3)
1154	FLUX_V(I,L) = AK_ALPHA(I,2)K_ALPHA(I,2)M_ALPHA(I,L,2)
1155	^ - AK_ALPHA(I,4)K_ALPHA(I,4)M_ALPHA(I,L,4)
1156	10 CONTINUE
1157	20 CONTINUE
1158	END IF
1159	C
1160	C Case with 8 azimuths.
1161	C
1162	IF (NAZ.EQ.8) THEN
1163	DO 40 L = LEV1,LEV2
1164	DO 30 I = IL1,IL2
1165	PROD2 = AK_ALPHA(I,2)K_ALPHA(I,2)M_ALPHA(I,L,2)
1166	PROD4 = AK_ALPHA(I,4)K_ALPHA(I,4)M_ALPHA(I,L,4)
1167	PROD6 = AK_ALPHA(I,6)K_ALPHA(I,6)M_ALPHA(I,L,6)
1168	PROD8 = AK_ALPHA(I,8)K_ALPHA(I,8)M_ALPHA(I,L,8)
1169	FLUX_U(I,L) =
1170	^ AK_ALPHA(I,1)K_ALPHA(I,1)M_ALPHA(I,L,1)
1171	^ - AK_ALPHA(I,5)K_ALPHA(I,5)M_ALPHA(I,L,5)
1172	^ + COS45 * ( PROD2 - PROD4 - PROD6 + PROD8 )
1173	FLUX_V(I,L) =
1174	^ AK_ALPHA(I,3)K_ALPHA(I,3)M_ALPHA(I,L,3)
1175	^ - AK_ALPHA(I,7)K_ALPHA(I,7)M_ALPHA(I,L,7)
1176	^ + COS45 * ( PROD2 + PROD4 - PROD6 - PROD8 )
1177	30 CONTINUE
1178	40 CONTINUE
1179	END IF
1180	C
1181	END IF
1182	C
1183	C Sum over azimuths for case where SLOPE not equal to 1.
1184	C
1185	IF (SLOPE.NE.1.) THEN
1186	C
1187	C Case with 4 azimuths.
1188	C
1189	IF (NAZ.EQ.4) THEN
1190	DO 60 L = LEV1,LEV2
1191	DO 50 I = IL1,IL2
1192	FLUX_U(I,L) =
1193	^ AK_ALPHA(I,1)K_ALPHA(I,1)M_ALPHA(I,L,1)**SLOPE
1194	^ - AK_ALPHA(I,3)K_ALPHA(I,3)M_ALPHA(I,L,3)**SLOPE
1195	FLUX_V(I,L) =
1196	^ AK_ALPHA(I,2)K_ALPHA(I,2)M_ALPHA(I,L,2)**SLOPE
1197	^ - AK_ALPHA(I,4)K_ALPHA(I,4)M_ALPHA(I,L,4)**SLOPE
1198	50 CONTINUE
1199	60 CONTINUE
1200	END IF
1201	C
1202	C Case with 8 azimuths.
1203	C
1204	IF (NAZ.EQ.8) THEN
1205	DO 80 L = LEV1,LEV2
1206	DO 70 I = IL1,IL2
1207	PROD2 = AK_ALPHA(I,2)K_ALPHA(I,2)M_ALPHA(I,L,2)**SLOPE
1208	PROD4 = AK_ALPHA(I,4)K_ALPHA(I,4)M_ALPHA(I,L,4)**SLOPE
1209	PROD6 = AK_ALPHA(I,6)K_ALPHA(I,6)M_ALPHA(I,L,6)**SLOPE
1210	PROD8 = AK_ALPHA(I,8)K_ALPHA(I,8)M_ALPHA(I,L,8)**SLOPE
1211	FLUX_U(I,L) =
1212	^ AK_ALPHA(I,1)K_ALPHA(I,1)M_ALPHA(I,L,1)**SLOPE
1213	^ - AK_ALPHA(I,5)K_ALPHA(I,5)M_ALPHA(I,L,5)**SLOPE
1214	^ + COS45 * ( PROD2 - PROD4 - PROD6 + PROD8 )
1215	FLUX_V(I,L) =
1216	^ AK_ALPHA(I,3)K_ALPHA(I,3)M_ALPHA(I,L,3)**SLOPE
1217	^ - AK_ALPHA(I,7)K_ALPHA(I,7)M_ALPHA(I,L,7)**SLOPE
1218	^ + COS45 * ( PROD2 + PROD4 - PROD6 - PROD8 )
1219	70 CONTINUE
1220	80 CONTINUE
1221	END IF
1222	C
1223	END IF
1224	C
1225	C Calculate flux from sum.
1226	C
1227	DO 100 L = LEV1,LEV2
1228	DO 90 I = IL1,IL2
1229	FLUX_U(I,L) = FLUX_U(I,L) * DENSB(I) / SLOPE
1230	FLUX_V(I,L) = FLUX_V(I,L) * DENSB(I) / SLOPE
1231	90 CONTINUE
1232	100 CONTINUE
1233	C
1234	C Calculate drag at intermediate levels using centered differences
1235	C
1236	DO 120 L = LEV1P,LEV2M
1237	DO 110 I = IL1,IL2
1238	IF (LORMS(I)) THEN
1239	ccc DENDZ2 = DENSITY(I,L) * ( ALT(I,L+1) - ALT(I,L-1) )
1240	DENDZ2 = DENSITY(I,L) * ( ALT(I,L-1) - ALT(I,L) )
1241	ccc DRAG_U(I,L) = - ( FLUX_U(I,L+1) - FLUX_U(I,L-1) ) / DENDZ2
1242	DRAG_U(I,L) = - ( FLUX_U(I,L-1) - FLUX_U(I,L) ) / DENDZ2
1243	ccc DRAG_V(I,L) = - ( FLUX_V(I,L+1) - FLUX_V(I,L-1) ) / DENDZ2
1244	DRAG_V(I,L) = - ( FLUX_V(I,L-1) - FLUX_V(I,L) ) / DENDZ2
1245
1246	ENDIF
1247	110 CONTINUE
1248	120 CONTINUE
1249	C
1250	C Drag at first and last levels using one-side differences.
1251	C
1252	DO 130 I = IL1,IL2
1253	IF (LORMS(I)) THEN
1254	DENDZ = DENSITY(I,LEV1) * ( ALT(I,LEV1) - ALT(I,LEV1P) )
1255	DRAG_U(I,LEV1) = FLUX_U(I,LEV1) / DENDZ
1256	DRAG_V(I,LEV1) = FLUX_V(I,LEV1) / DENDZ
1257	ENDIF
1258	130 CONTINUE
1259	DO 140 I = IL1,IL2
1260	IF (LORMS(I)) THEN
1261	DENDZ = DENSITY(I,LEV2) * ( ALT(I,LEV2M) - ALT(I,LEV2) )
1262	DRAG_U(I,LEV2) = - ( FLUX_U(I,LEV2M) - FLUX_U(I,LEV2) ) / DENDZ
1263	DRAG_V(I,LEV2) = - ( FLUX_V(I,LEV2M) - FLUX_V(I,LEV2) ) / DENDZ
1264	ENDIF
1265	140 CONTINUE
1266	IF (NLEVS .GT. LEV2) THEN
1267	DO 150 I = IL1,IL2
1268	IF (LORMS(I)) THEN
1269	DENDZ = DENSITY(I,LEV2P) * ( ALT(I,LEV2) - ALT(I,LEV2P) )
1270	DRAG_U(I,LEV2P) = - FLUX_U(I,LEV2) / DENDZ
1271	DRAG_V(I,LEV2P) = - FLUX_V(I,LEV2) / DENDZ
1272	ENDIF
1273	150 CONTINUE
1274	ENDIF
1275	C
1276	RETURN
1277	C-----------------------------------------------------------------------
1278	END
1279
1280	SUBROUTINE HINES_HEAT (HEAT,DIFFCO,M_ALPHA,DMDZ_ALPHA,
1281	1 AK_ALPHA,K_ALPHA,BVFREQ,DENSITY,DENSB,
1282	2 SIGMA_T,VISC_MOL,KSTAR,SLOPE,F2,F3,F5,F6,
1283	3 NAZ,IL1,IL2,LEV1,LEV2,NLONS,NLEVS,NAZMTH)
1284	C
1285	C This routine calculates the gravity wave induced heating and
1286	C diffusion coefficient on a longitude by altitude grid for
1287	C the Hines' Doppler spread gravity wave drag parameterization scheme.
1288	C
1289	C Aug. 6/95 - C. McLandress
1290	C
1291	C Output arguements:
1292	C
1293	C * HEAT = gravity wave heating (K/sec).
1294	C * DIFFCO = diffusion coefficient (m^2/sec)
1295	C
1296	C Input arguements:
1297	C
1298	C * M_ALPHA = cutoff vertical wavenumber (1/m).
1299	C * DMDZ_ALPHA = vertical derivative of cutoff wavenumber.
1300	C * AK_ALPHA = spectral amplitude factor of each azimuth
1301	C (i.e., {AjKj} in m^4/s^2).
1302	C * K_ALPHA = horizontal wavenumber of each azimuth (1/m).
1303	C * BVFREQ = background Brunt Vassala frequency (rad/sec).
1304	C * DENSITY = background density (kg/m^3).
1305	C * DENSB = background density at bottom level (kg/m^3).
1306	C * SIGMA_T = total rms horizontal wind (m/s).
1307	C * VISC_MOL = molecular viscosity (m^2/s).
1308	C * KSTAR = typical gravity wave horizontal wavenumber (1/m).
1309	C * SLOPE = slope of incident vertical wavenumber spectrum.
1310	C * F2,F3,F5,F6 = Hines's fudge factors.
1311	C * NAZ = actual number of horizontal azimuths used.
1312	C * IL1 = first longitudinal index to use (IL1 >= 1).
1313	C * IL2 = last longitudinal index to use (IL1 <= IL2 <= NLONS).
1314	C * LEV1 = first altitude level to use (LEV1 >=1).
1315	C * LEV2 = last altitude level to use (LEV1 < LEV2 <= NLEVS).
1316	C * NLONS = number of longitudes.
1317	C * NLEVS = number of vertical levels.
1318	C * NAZMTH = azimuthal array dimension (NAZMTH >= NAZ).
1319	C
1320	INTEGER NAZ, IL1, IL2, LEV1, LEV2, NLONS, NLEVS, NAZMTH
1321	REAL KSTAR(NLONS), SLOPE, F2, F3, F5, F6
1322	REAL HEAT(NLONS,NLEVS), DIFFCO(NLONS,NLEVS)
1323	REAL M_ALPHA(NLONS,NLEVS,NAZMTH), DMDZ_ALPHA(NLONS,NLEVS,NAZMTH)
1324	REAL AK_ALPHA(NLONS,NAZMTH), K_ALPHA(NLONS,NAZMTH)
1325	REAL BVFREQ(NLONS,NLEVS), DENSITY(NLONS,NLEVS), DENSB(NLONS)
1326	REAL SIGMA_T(NLONS,NLEVS), VISC_MOL(NLONS,NLEVS)
1327	C
1328	C Internal variables.
1329	C
1330	INTEGER I, L, N
1331	REAL M_SUB_M_TURB, M_SUB_M_MOL, M_SUB_M, HEATNG
1332	REAL VISC, VISC_MIN, CPGAS, SM1
1333	C
1334	C specific heat at constant pressure
1335	C
1336	DATA CPGAS / 1004. /
1337	C
1338	C minimum permissible viscosity
1339	C
1340	DATA VISC_MIN / 1.E-10 /
1341	C-----------------------------------------------------------------------
1342	C
1343	C Initialize heating array.
1344	C
1345	DO 20 L = 1,NLEVS
1346	DO 10 I = 1,NLONS
1347	HEAT(I,L) = 0.
1348	10 CONTINUE
1349	20 CONTINUE
1350	C
1351	C Perform sum over azimuths for case where SLOPE = 1.
1352	C
1353	IF (SLOPE.EQ.1.) THEN
1354	DO 50 N = 1,NAZ
1355	DO 40 L = LEV1,LEV2
1356	DO 30 I = IL1,IL2
1357	HEAT(I,L) = HEAT(I,L) + AK_ALPHA(I,N) * K_ALPHA(I,N)
1358	^ * DMDZ_ALPHA(I,L,N)
1359	30 CONTINUE
1360	40 CONTINUE
1361	50 CONTINUE
1362	END IF
1363	C
1364	C Perform sum over azimuths for case where SLOPE not 1.
1365	C
1366	IF (SLOPE.NE.1.) THEN
1367	SM1 = SLOPE - 1.
1368	DO 80 N = 1,NAZ
1369	DO 70 L = LEV1,LEV2
1370	DO 60 I = IL1,IL2
1371	HEAT(I,L) = HEAT(I,L) + AK_ALPHA(I,N) * K_ALPHA(I,N)
1372	^ * M_ALPHA(I,L,N)*SM1 DMDZ_ALPHA(I,L,N)
1373	60 CONTINUE
1374	70 CONTINUE
1375	80 CONTINUE
1376	END IF
1377	C
1378	C Heating and diffusion.
1379	C
1380	DO 100 L = LEV1,LEV2
1381	DO 90 I = IL1,IL2
1382	C
1383	C Maximum permissible value of cutoff wavenumber is the smaller
1384	C of the instability-induced wavenumber (M_SUB_M_TURB) and
1385	C that imposed by molecular viscosity (M_SUB_M_MOL).
1386	C
1387	VISC = AMAX1 ( VISC_MOL(I,L), VISC_MIN )
1388	M_SUB_M_TURB = BVFREQ(I,L) / ( F2 * SIGMA_T(I,L) )
1389	M_SUB_M_MOL = (BVFREQ(I,L)KSTAR(I)/VISC)*0.33333333/F3
1390	M_SUB_M = AMIN1 ( M_SUB_M_TURB, M_SUB_M_MOL )
1391	C
1392	HEATNG = - HEAT(I,L) * F5 * BVFREQ(I,L) / M_SUB_M
1393	^ * DENSB(I) / DENSITY(I,L)
1394	DIFFCO(I,L) = F6 * HEATNG0.33333333 / M_SUB_M1.33333333
1395	HEAT(I,L) = HEATNG / CPGAS
1396	C
1397	90 CONTINUE
1398	100 CONTINUE
1399	C
1400	RETURN
1401	C-----------------------------------------------------------------------
1402	END
1403
1404	SUBROUTINE HINES_SIGMA (SIGMA_T,SIGMA_ALPHA,SIGSQH_ALPHA,
1405	1 NAZ,LEV,IL1,IL2,NLONS,NLEVS,NAZMTH)
1406	C
1407	C This routine calculates the total rms and azimuthal rms horizontal
1408	C velocities at a given level on a longitude by altitude grid for
1409	C the Hines' Doppler spread GWD parameterization scheme.
1410	C NOTE: only four or eight azimuths can be used.
1411	C
1412	C Aug. 7/95 - C. McLandress
1413	C
1414	C Output arguements:
1415	C
1416	C * SIGMA_T = total rms horizontal wind (m/s).
1417	C * SIGMA_ALPHA = total rms wind in each azimuth (m/s).
1418	C
1419	C Input arguements:
1420	C
1421	C * SIGSQH_ALPHA = portion of wind variance from waves having wave
1422	C * normals in the alpha azimuth (m/s).
1423	C * NAZ = actual number of horizontal azimuths used (must be 4 or 8).
1424	C * LEV = altitude level to process.
1425	C * IL1 = first longitudinal index to use (IL1 >= 1).
1426	C * IL2 = last longitudinal index to use (IL1 <= IL2 <= NLONS).
1427	C * NLONS = number of longitudes.
1428	C * NLEVS = number of vertical levels.
1429	C * NAZMTH = azimuthal array dimension (NAZMTH >= NAZ).
1430	C
1431	C Subroutine arguements.
1432	C
1433	INTEGER LEV, NAZ, IL1, IL2
1434	INTEGER NLONS, NLEVS, NAZMTH
1435	REAL SIGMA_T(NLONS,NLEVS)
1436	REAL SIGMA_ALPHA(NLONS,NLEVS,NAZMTH)
1437	REAL SIGSQH_ALPHA(NLONS,NLEVS,NAZMTH)
1438	C
1439	C Internal variables.
1440	C
1441	INTEGER I, N
1442	REAL SUM_EVEN, SUM_ODD
1443	C-----------------------------------------------------------------------
1444	C
1445	C Calculate azimuthal rms velocity for the 4 azimuth case.
1446	C
1447	IF (NAZ.EQ.4) THEN
1448	DO 10 I = IL1,IL2
1449	SIGMA_ALPHA(I,LEV,1) = SQRT ( SIGSQH_ALPHA(I,LEV,1)
1450	^ + SIGSQH_ALPHA(I,LEV,3) )
1451	SIGMA_ALPHA(I,LEV,2) = SQRT ( SIGSQH_ALPHA(I,LEV,2)
1452	^ + SIGSQH_ALPHA(I,LEV,4) )
1453	SIGMA_ALPHA(I,LEV,3) = SIGMA_ALPHA(I,LEV,1)
1454	SIGMA_ALPHA(I,LEV,4) = SIGMA_ALPHA(I,LEV,2)
1455	10 CONTINUE
1456	END IF
1457	C
1458	C Calculate azimuthal rms velocity for the 8 azimuth case.
1459	C
1460	IF (NAZ.EQ.8) THEN
1461	DO 20 I = IL1,IL2
1462	SUM_ODD = ( SIGSQH_ALPHA(I,LEV,1)
1463	^ + SIGSQH_ALPHA(I,LEV,3)
1464	^ + SIGSQH_ALPHA(I,LEV,5)
1465	^ + SIGSQH_ALPHA(I,LEV,7) ) / 2.
1466	SUM_EVEN = ( SIGSQH_ALPHA(I,LEV,2)
1467	^ + SIGSQH_ALPHA(I,LEV,4)
1468	^ + SIGSQH_ALPHA(I,LEV,6)
1469	^ + SIGSQH_ALPHA(I,LEV,8) ) / 2.
1470	SIGMA_ALPHA(I,LEV,1) = SQRT ( SIGSQH_ALPHA(I,LEV,1)
1471	^ + SIGSQH_ALPHA(I,LEV,5) + SUM_EVEN )
1472	SIGMA_ALPHA(I,LEV,2) = SQRT ( SIGSQH_ALPHA(I,LEV,2)
1473	^ + SIGSQH_ALPHA(I,LEV,6) + SUM_ODD )
1474	SIGMA_ALPHA(I,LEV,3) = SQRT ( SIGSQH_ALPHA(I,LEV,3)
1475	^ + SIGSQH_ALPHA(I,LEV,7) + SUM_EVEN )
1476	SIGMA_ALPHA(I,LEV,4) = SQRT ( SIGSQH_ALPHA(I,LEV,4)
1477	^ + SIGSQH_ALPHA(I,LEV,8) + SUM_ODD )
1478	SIGMA_ALPHA(I,LEV,5) = SIGMA_ALPHA(I,LEV,1)
1479	SIGMA_ALPHA(I,LEV,6) = SIGMA_ALPHA(I,LEV,2)
1480	SIGMA_ALPHA(I,LEV,7) = SIGMA_ALPHA(I,LEV,3)
1481	SIGMA_ALPHA(I,LEV,8) = SIGMA_ALPHA(I,LEV,4)
1482	20 CONTINUE
1483	END IF
1484	C
1485	C Calculate total rms velocity.
1486	C
1487	DO 50 I = IL1,IL2
1488	SIGMA_T(I,LEV) = 0.
1489	50 CONTINUE
1490	DO 70 N = 1,NAZ
1491	DO 60 I = IL1,IL2
1492	SIGMA_T(I,LEV) = SIGMA_T(I,LEV) + SIGSQH_ALPHA(I,LEV,N)
1493	60 CONTINUE
1494	70 CONTINUE
1495	DO 80 I = IL1,IL2
1496	SIGMA_T(I,LEV) = SQRT ( SIGMA_T(I,LEV) )
1497	80 CONTINUE
1498	C
1499	RETURN
1500	C-----------------------------------------------------------------------
1501	END
1502
1503	SUBROUTINE HINES_INTGRL (I_ALPHA,V_ALPHA,M_ALPHA,BVFB,SLOPE,
1504	1 NAZ,LEV,IL1,IL2,NLONS,NLEVS,NAZMTH,
1505	2 LORMS)
1506	C
1507	C This routine calculates the vertical wavenumber integral
1508	C for a single vertical level at each azimuth on a longitude grid
1509	C for the Hines' Doppler spread GWD parameterization scheme.
1510	C NOTE: (1) only spectral slopes of 1, 1.5 or 2 are permitted.
1511	C (2) the integral is written in terms of the product QM
1512	C which by construction is always less than 1. Series
1513	C solutions are used for small \|QM\| and analytical solutions
1514	C for remaining values.
1515	C
1516	C Aug. 8/95 - C. McLandress
1517	C
1518	C Output arguement:
1519	C
1520	C * I_ALPHA = Hines' integral.
1521	C
1522	C Input arguements:
1523	C
1524	C * V_ALPHA = azimuthal wind component (m/s).
1525	C * M_ALPHA = azimuthal cutoff vertical wavenumber (1/m).
1526	C * BVFB = background Brunt Vassala frequency at model bottom.
1527	C * SLOPE = slope of initial vertical wavenumber spectrum
1528	C * (must use SLOPE = 1., 1.5 or 2.)
1529	C * NAZ = actual number of horizontal azimuths used.
1530	C * LEV = altitude level to process.
1531	C * IL1 = first longitudinal index to use (IL1 >= 1).
1532	C * IL2 = last longitudinal index to use (IL1 <= IL2 <= NLONS).
1533	C * NLONS = number of longitudes.
1534	C * NLEVS = number of vertical levels.
1535	C * NAZMTH = azimuthal array dimension (NAZMTH >= NAZ).
1536	C
1537	C * LORMS = .TRUE. for drag computation
1538	C
1539	C Constants in DATA statements:
1540	C
1541	C * QMIN = minimum value of Q_ALPHA (avoids indeterminant form of integral)
1542	C * QM_MIN = minimum value of Q_ALPHA * M_ALPHA (used to avoid numerical
1543	C * problems).
1544	C
1545	INTEGER LEV, NAZ, IL1, IL2, NLONS, NLEVS, NAZMTH
1546	REAL I_ALPHA(NLONS,NAZMTH)
1547	REAL V_ALPHA(NLONS,NLEVS,NAZMTH)
1548	REAL M_ALPHA(NLONS,NLEVS,NAZMTH)
1549	REAL BVFB(NLONS), SLOPE
1550	C
1551	LOGICAL LORMS(NLONS)
1552	C
1553	C Internal variables.
1554	C
1555	INTEGER I, N
1556	REAL Q_ALPHA, QM, SQRTQM, Q_MIN, QM_MIN
1557	C
1558	DATA Q_MIN / 1.0 /, QM_MIN / 0.01 /
1559	C-----------------------------------------------------------------------
1560	C
1561	C For integer value SLOPE = 1.
1562	C
1563	IF (SLOPE .EQ. 1.) THEN
1564	C
1565	DO 20 N = 1,NAZ
1566	DO 10 I = IL1,IL2
1567	IF (LORMS(I)) THEN
1568	C
1569	Q_ALPHA = V_ALPHA(I,LEV,N) / BVFB(I)
1570	QM = Q_ALPHA * M_ALPHA(I,LEV,N)
1571	C
1572	C If \|QM\| is small then use first 4 terms series of Taylor series
1573	C expansion of integral in order to avoid indeterminate form of integral,
1574	C otherwise use analytical form of integral.
1575	C
1576	IF (ABS(Q_ALPHA).LT.Q_MIN .OR. ABS(QM).LT.QM_MIN) THEN
1577	IF (Q_ALPHA .EQ. 0.) THEN
1578	I_ALPHA(I,N) = M_ALPHA(I,LEV,N)**2 / 2.
1579	ELSE
1580	I_ALPHA(I,N) = ( QM2/2. + QM3/3. + QM**4/4.
1581	^ + QM5/5. ) / Q_ALPHA2
1582	END IF
1583	ELSE
1584	I_ALPHA(I,N) = - ( ALOG(1.-QM) + QM ) / Q_ALPHA**2
1585	END IF
1586	C
1587	ENDIF
1588	10 CONTINUE
1589	20 CONTINUE
1590	C
1591	END IF
1592	C
1593	C For integer value SLOPE = 2.
1594	C
1595	IF (SLOPE .EQ. 2.) THEN
1596	C
1597	DO 40 N = 1,NAZ
1598	DO 30 I = IL1,IL2
1599	IF (LORMS(I)) THEN
1600	C
1601	Q_ALPHA = V_ALPHA(I,LEV,N) / BVFB(I)
1602	QM = Q_ALPHA * M_ALPHA(I,LEV,N)
1603	C
1604	C If \|QM\| is small then use first 4 terms series of Taylor series
1605	C expansion of integral in order to avoid indeterminate form of integral,
1606	C otherwise use analytical form of integral.
1607	C
1608	IF (ABS(Q_ALPHA).LT.Q_MIN .OR. ABS(QM).LT.QM_MIN) THEN
1609	IF (Q_ALPHA .EQ. 0.) THEN
1610	I_ALPHA(I,N) = M_ALPHA(I,LEV,N)**3 / 3.
1611	ELSE
1612	I_ALPHA(I,N) = ( QM3/3. + QM4/4. + QM**5/5.
1613	^ + QM6/6. ) / Q_ALPHA3
1614	END IF
1615	ELSE
1616	I_ALPHA(I,N) = - ( ALOG(1.-QM) + QM + QM**2/2.)
1617	^ / Q_ALPHA**3
1618	END IF
1619	C
1620	ENDIF
1621	30 CONTINUE
1622	40 CONTINUE
1623	C
1624	END IF
1625	C
1626	C For real value SLOPE = 1.5
1627	C
1628	IF (SLOPE .EQ. 1.5) THEN
1629	C
1630	DO 60 N = 1,NAZ
1631	DO 50 I = IL1,IL2
1632	IF (LORMS(I)) THEN
1633	C
1634	Q_ALPHA = V_ALPHA(I,LEV,N) / BVFB(I)
1635	QM = Q_ALPHA * M_ALPHA(I,LEV,N)
1636	C
1637	C If \|QM\| is small then use first 4 terms series of Taylor series
1638	C expansion of integral in order to avoid indeterminate form of integral,
1639	C otherwise use analytical form of integral.
1640	C
1641	IF (ABS(Q_ALPHA).LT.Q_MIN .OR. ABS(QM).LT.QM_MIN) THEN
1642	IF (Q_ALPHA .EQ. 0.) THEN
1643	I_ALPHA(I,N) = M_ALPHA(I,LEV,N)**2.5 / 2.5
1644	ELSE
1645	I_ALPHA(I,N) = ( QM/2.5 + QM**2/3.5
1646	^ + QM3/4.5 + QM4/5.5 )
1647	^ * M_ALPHA(I,LEV,N)**1.5 / Q_ALPHA
1648	END IF
1649	ELSE
1650	QM = ABS(QM)
1651	SQRTQM = SQRT(QM)
1652	IF (Q_ALPHA .GE. 0.) THEN
1653	I_ALPHA(I,N) = ( ALOG( (1.+SQRTQM)/(1.-SQRTQM) )
1654	^ -2.SQRTQM(1.+QM/3.) ) / Q_ALPHA**2.5
1655	ELSE
1656	I_ALPHA(I,N) = 2. * ( ATAN(SQRTQM) + SQRTQM*(QM/3.-1.) )
1657	^ / ABS(Q_ALPHA)**2.5
1658	END IF
1659	END IF
1660	C
1661	ENDIF
1662	50 CONTINUE
1663	60 CONTINUE
1664	C
1665	END IF
1666	C
1667	C If integral is negative (which in principal should not happen) then
1668	C print a message and some info since execution will abort when calculating
1669	C the variances.
1670	C
1671	c DO 80 N = 1,NAZ
1672	c DO 70 I = IL1,IL2
1673	c IF (I_ALPHA(I,N).LT.0.) THEN
1674	c WRITE (6,*)
1675	c WRITE (6,) '*****************************'
1676	c WRITE (6,*) 'Hines integral I_ALPHA < 0 '
1677	c WRITE (6,*) ' longitude I=',I
1678	c WRITE (6,*) ' azimuth N=',N
1679	c WRITE (6,*) ' level LEV=',LEV
1680	c WRITE (6,*) ' I_ALPHA =',I_ALPHA(I,N)
1681	c WRITE (6,*) ' V_ALPHA =',V_ALPHA(I,LEV,N)
1682	c WRITE (6,*) ' M_ALPHA =',M_ALPHA(I,LEV,N)
1683	c WRITE (6,*) ' Q_ALPHA =',V_ALPHA(I,LEV,N) / BVFB(I)
1684	c WRITE (6,*) ' QM =',V_ALPHA(I,LEV,N) / BVFB(I)
1685	c ^ * M_ALPHA(I,LEV,N)
1686	c WRITE (6,) '*****************************'
1687	c END IF
1688	c 70 CONTINUE
1689	c 80 CONTINUE
1690	C
1691	RETURN
1692	C-----------------------------------------------------------------------
1693	END
1694
1695	SUBROUTINE HINES_SETUP (NAZ,SLOPE,F1,F2,F3,F5,F6,KSTAR,
1696	1 ICUTOFF,ALT_CUTOFF,SMCO,NSMAX,IHEATCAL,
1697	2 K_ALPHA,IERROR,NMESSG,NLONS,NAZMTH,COSLAT)
1698	C
1699	C This routine specifies various parameters needed for the
1700	C the Hines' Doppler spread gravity wave drag parameterization scheme.
1701	C
1702	C Aug. 8/95 - C. McLandress
1703	C
1704	C Output arguements:
1705	C
1706	C * NAZ = actual number of horizontal azimuths used
1707	C * (code set up presently for only NAZ = 4 or 8).
1708	C * SLOPE = slope of incident vertical wavenumber spectrum
1709	C * (code set up presently for SLOPE 1., 1.5 or 2.).
1710	C * F1 = "fudge factor" used in calculation of trial value of
1711	C * azimuthal cutoff wavenumber M_ALPHA (1.2 <= F1 <= 1.9).
1712	C * F2 = "fudge factor" used in calculation of maximum
1713	C * permissible instabiliy-induced cutoff wavenumber
1714	C * M_SUB_M_TURB (0.1 <= F2 <= 1.4).
1715	C * F3 = "fudge factor" used in calculation of maximum
1716	C * permissible molecular viscosity-induced cutoff wavenumber
1717	C * M_SUB_M_MOL (0.1 <= F2 <= 1.4).
1718	C * F5 = "fudge factor" used in calculation of heating rate
1719	C * (1 <= F5 <= 3).
1720	C * F6 = "fudge factor" used in calculation of turbulent
1721	C * diffusivity coefficient.
1722	C * KSTAR = typical gravity wave horizontal wavenumber (1/m)
1723	C * used in calculation of M_SUB_M_TURB.
1724	C * ICUTOFF = 1 to exponentially damp off GWD, heating and diffusion
1725	C * arrays above ALT_CUTOFF; otherwise arrays not modified.
1726	C * ALT_CUTOFF = altitude in meters above which exponential decay applied.
1727	C * SMCO = smoother used to smooth cutoff vertical wavenumbers
1728	C * and total rms winds before calculating drag or heating.
1729	C * (==> a 1:SMCO:1 stencil used; SMCO >= 1.).
1730	C * NSMAX = number of times smoother applied ( >= 1),
1731	C * = 0 means no smoothing performed.
1732	C * IHEATCAL = 1 to calculate heating rates and diffusion coefficient.
1733	C * = 0 means only drag and flux calculated.
1734	C * K_ALPHA = horizontal wavenumber of each azimuth (1/m) which
1735	C * is set here to KSTAR.
1736	C * IERROR = error flag.
1737	C * = 0 no errors.
1738	C * = 10 ==> NAZ > NAZMTH
1739	C * = 20 ==> invalid number of azimuths (NAZ must be 4 or 8).
1740	C * = 30 ==> invalid slope (SLOPE must be 1., 1.5 or 2.).
1741	C * = 40 ==> invalid smoother (SMCO must be >= 1.)
1742	C
1743	C Input arguements:
1744	C
1745	C * NMESSG = output unit number where messages to be printed.
1746	C * NLONS = number of longitudes.
1747	C * NAZMTH = azimuthal array dimension (NAZMTH >= NAZ).
1748	C
1749	INTEGER NAZ, NLONS, NAZMTH, IHEATCAL, ICUTOFF
1750	INTEGER NMESSG, NSMAX, IERROR
1751	REAL KSTAR(NLONS), SLOPE, F1, F2, F3, F5, F6, ALT_CUTOFF, SMCO
1752	REAL K_ALPHA(NLONS,NAZMTH),COSLAT(NLONS)
1753	REAL KSMIN, KSMAX
1754	C
1755	C Internal variables.
1756	C
1757	INTEGER I, N
1758	C-----------------------------------------------------------------------
1759	C
1760	C Specify constants.
1761	C
1762	NAZ = 8
1763	SLOPE = 1.
1764	F1 = 1.5
1765	F2 = 0.3
1766	F3 = 1.0
1767	F5 = 3.0
1768	F6 = 1.0
1769	KSMIN = 1.E-5
1770	KSMAX = 1.E-4
1771	DO I=1,NLONS
1772	KSTAR(I) = KSMIN/( COSLAT(I)+(KSMIN/KSMAX) )
1773	ENDDO
1774	ICUTOFF = 1
1775	ALT_CUTOFF = 105.E3
1776	SMCO = 2.0
1777	c SMCO = 1.0
1778	NSMAX = 5
1779	c NSMAX = 2
1780	IHEATCAL = 0
1781	C
1782	C Print information to output file.
1783	C
1784	c WRITE (NMESSG,6000)
1785	c 6000 FORMAT (/' Subroutine HINES_SETUP:')
1786	c WRITE (NMESSG,*) ' SLOPE = ', SLOPE
1787	c WRITE (NMESSG,*) ' NAZ = ', NAZ
1788	c WRITE (NMESSG,*) ' F1,F2,F3 = ', F1, F2, F3
1789	c WRITE (NMESSG,*) ' F5,F6 = ', F5, F6
1790	c WRITE (NMESSG,*) ' KSTAR = ', KSTAR
1791	c > ,' COSLAT = ', COSLAT
1792	c IF (ICUTOFF .EQ. 1) THEN
1793	c WRITE (NMESSG,*) ' Drag exponentially damped above ',
1794	c & ALT_CUTOFF/1.E3
1795	c END IF
1796	c IF (NSMAX.LT.1 ) THEN
1797	c WRITE (NMESSG,*) ' No smoothing of cutoff wavenumbers, etc'
1798	c ELSE
1799	c WRITE (NMESSG,*) ' Cutoff wavenumbers and sig_t smoothed:'
1800	c WRITE (NMESSG,*) ' SMCO =', SMCO
1801	c WRITE (NMESSG,*) ' NSMAX =', NSMAX
1802	c END IF
1803	C
1804	C Check that things are setup correctly and log error if not
1805	C
1806	IERROR = 0
1807	IF (NAZ .GT. NAZMTH) IERROR = 10
1808	IF (NAZ.NE.4 .AND. NAZ.NE.8) IERROR = 20
1809	IF (SLOPE.NE.1. .AND. SLOPE.NE.1.5 .AND. SLOPE.NE.2.) IERROR = 30
1810	IF (SMCO .LT. 1.) IERROR = 40
1811	C
1812	C Use single value for azimuthal-dependent horizontal wavenumber.
1813	C
1814	DO 20 N = 1,NAZ
1815	DO 10 I = 1,NLONS
1816	K_ALPHA(I,N) = KSTAR(I)
1817	10 CONTINUE
1818	20 CONTINUE
1819	C
1820	RETURN
1821	C-----------------------------------------------------------------------
1822	END
1823
1824	SUBROUTINE HINES_PRINT (FLUX_U,FLUX_V,DRAG_U,DRAG_V,ALT,SIGMA_T,
1825	1 SIGMA_ALPHA,V_ALPHA,M_ALPHA,
1826	2 IU_PRINT,IV_PRINT,NMESSG,
1827	3 ILPRT1,ILPRT2,LEVPRT1,LEVPRT2,
1828	4 NAZ,NLONS,NLEVS,NAZMTH)
1829	C
1830	C Print out altitude profiles of various quantities from
1831	C Hines' Doppler spread gravity wave drag parameterization scheme.
1832	C (NOTE: only for NAZ = 4 or 8).
1833	C
1834	C Aug. 8/95 - C. McLandress
1835	C
1836	C Input arguements:
1837	C
1838	C * IU_PRINT = 1 to print out values in east-west direction.
1839	C * IV_PRINT = 1 to print out values in north-south direction.
1840	C * NMESSG = unit number for printed output.
1841	C * ILPRT1 = first longitudinal index to print.
1842	C * ILPRT2 = last longitudinal index to print.
1843	C * LEVPRT1 = first altitude level to print.
1844	C * LEVPRT2 = last altitude level to print.
1845	C
1846	INTEGER NAZ, ILPRT1, ILPRT2, LEVPRT1, LEVPRT2
1847	INTEGER NLONS, NLEVS, NAZMTH
1848	INTEGER IU_PRINT, IV_PRINT, NMESSG
1849	REAL FLUX_U(NLONS,NLEVS), FLUX_V(NLONS,NLEVS)
1850	REAL DRAG_U(NLONS,NLEVS), DRAG_V(NLONS,NLEVS)
1851	REAL ALT(NLONS,NLEVS), SIGMA_T(NLONS,NLEVS)
1852	REAL SIGMA_ALPHA(NLONS,NLEVS,NAZMTH)
1853	REAL V_ALPHA(NLONS,NLEVS,NAZMTH), M_ALPHA(NLONS,NLEVS,NAZMTH)
1854	C
1855	C Internal variables.
1856	C
1857	INTEGER N_EAST, N_WEST, N_NORTH, N_SOUTH
1858	INTEGER I, L
1859	C-----------------------------------------------------------------------
1860	C
1861	C Azimuthal indices of cardinal directions.
1862	C
1863	N_EAST = 1
1864	IF (NAZ.EQ.4) THEN
1865	N_WEST = 3
1866	N_NORTH = 2
1867	N_SOUTH = 4
1868	ELSE IF (NAZ.EQ.8) THEN
1869	N_WEST = 5
1870	N_NORTH = 3
1871	N_SOUTH = 7
1872	END IF
1873	C
1874	C Print out values for range of longitudes.
1875	C
1876	DO 100 I = ILPRT1,ILPRT2
1877	C
1878	C Print east-west wind, sigmas, cutoff wavenumbers, flux and drag.
1879	C
1880	IF (IU_PRINT.EQ.1) THEN
1881	WRITE (NMESSG,*)
1882	WRITE (NMESSG,6001) I
1883	WRITE (NMESSG,6005)
1884	6001 FORMAT ( 'Hines GW (east-west) at longitude I =',I3)
1885	6005 FORMAT (15x,' U ',2x,'sig_E',2x,'sig_T',3x,'m_E',
1886	& 4x,'m_W',4x,'fluxU',5x,'gwdU')
1887	DO 10 L = LEVPRT1,LEVPRT2
1888	WRITE (NMESSG,6701) ALT(I,L)/1.E3, V_ALPHA(I,L,N_EAST),
1889	& SIGMA_ALPHA(I,L,N_EAST), SIGMA_T(I,L),
1890	& M_ALPHA(I,L,N_EAST)*1.E3,
1891	& M_ALPHA(I,L,N_WEST)*1.E3,
1892	& FLUX_U(I,L)1.E5, DRAG_U(I,L)24.*3600.
1893	10 CONTINUE
1894	6701 FORMAT (' z=',f7.2,1x,3f7.1,2f7.3,f9.4,f9.3)
1895	END IF
1896	C
1897	C Print north-south winds, sigmas, cutoff wavenumbers, flux and drag.
1898	C
1899	IF (IV_PRINT.EQ.1) THEN
1900	WRITE(NMESSG,*)
1901	WRITE(NMESSG,6002) I
1902	6002 FORMAT ( 'Hines GW (north-south) at longitude I =',I3)
1903	WRITE(NMESSG,6006)
1904	6006 FORMAT (15x,' V ',2x,'sig_N',2x,'sig_T',3x,'m_N',
1905	& 4x,'m_S',4x,'fluxV',5x,'gwdV')
1906	DO 20 L = LEVPRT1,LEVPRT2
1907	WRITE (NMESSG,6701) ALT(I,L)/1.E3, V_ALPHA(I,L,N_NORTH),
1908	& SIGMA_ALPHA(I,L,N_NORTH), SIGMA_T(I,L),
1909	& M_ALPHA(I,L,N_NORTH)*1.E3,
1910	& M_ALPHA(I,L,N_SOUTH)*1.E3,
1911	& FLUX_V(I,L)1.E5, DRAG_V(I,L)24.*3600.
1912	20 CONTINUE
1913	END IF
1914	C
1915	100 CONTINUE
1916	C
1917	RETURN
1918	C-----------------------------------------------------------------------
1919	END
1920
1921	SUBROUTINE HINES_EXP (DATA,DATA_ZMAX,ALT,ALT_EXP,IORDER,
1922	1 IL1,IL2,LEV1,LEV2,NLONS,NLEVS)
1923	C
1924	C This routine exponentially damps a longitude by altitude array
1925	C of data above a specified altitude.
1926	C
1927	C Aug. 13/95 - C. McLandress
1928	C
1929	C Output arguements:
1930	C
1931	C * DATA = modified data array.
1932	C
1933	C Input arguements:
1934	C
1935	C * DATA = original data array.
1936	C * ALT = altitudes.
1937	C * ALT_EXP = altitude above which exponential decay applied.
1938	C * IORDER = 1 means vertical levels are indexed from top down
1939	C * (i.e., highest level indexed 1 and lowest level NLEVS);
1940	C * .NE. 1 highest level is index NLEVS.
1941	C * IL1 = first longitudinal index to use (IL1 >= 1).
1942	C * IL2 = last longitudinal index to use (IL1 <= IL2 <= NLONS).
1943	C * LEV1 = first altitude level to use (LEV1 >=1).
1944	C * LEV2 = last altitude level to use (LEV1 < LEV2 <= NLEVS).
1945	C * NLONS = number of longitudes.
1946	C * NLEVS = number of vertical
1947	C
1948	C Input work arrays:
1949	C
1950	C * DATA_ZMAX = data values just above altitude ALT_EXP.
1951	C
1952	INTEGER IORDER, IL1, IL2, LEV1, LEV2, NLONS, NLEVS
1953	REAL ALT_EXP
1954	REAL DATA(NLONS,NLEVS), DATA_ZMAX(NLONS), ALT(NLONS,NLEVS)
1955	C
1956	C Internal variables.
1957	C
1958	INTEGER LEVBOT, LEVTOP, LINCR, I, L
1959	REAL HSCALE
1960	DATA HSCALE / 5.E3 /
1961	C-----------------------------------------------------------------------
1962	C
1963	C Index of lowest altitude level (bottom of drag calculation).
1964	C
1965	LEVBOT = LEV2
1966	LEVTOP = LEV1
1967	LINCR = 1
1968	IF (IORDER.NE.1) THEN
1969	LEVBOT = LEV1
1970	LEVTOP = LEV2
1971	LINCR = -1
1972	END IF
1973	C
1974	C Data values at first level above ALT_EXP.
1975	C
1976	DO 20 I = IL1,IL2
1977	DO 10 L = LEVTOP,LEVBOT,LINCR
1978	IF (ALT(I,L) .GE. ALT_EXP) THEN
1979	DATA_ZMAX(I) = DATA(I,L)
1980	END IF
1981	10 CONTINUE
1982	20 CONTINUE
1983	C
1984	C Exponentially damp field above ALT_EXP to model top at L=1.
1985	C
1986	DO 40 L = 1,LEV2
1987	DO 30 I = IL1,IL2
1988	IF (ALT(I,L) .GE. ALT_EXP) THEN
1989	DATA(I,L) = DATA_ZMAX(I) * EXP( (ALT_EXP-ALT(I,L))/HSCALE )
1990	END IF
1991	30 CONTINUE
1992	40 CONTINUE
1993	C
1994	RETURN
1995	C-----------------------------------------------------------------------
1996	END
1997
1998	SUBROUTINE VERT_SMOOTH (DATA,WORK,COEFF,NSMOOTH,
1999	1 IL1,IL2,LEV1,LEV2,NLONS,NLEVS)
2000	C
2001	C Smooth a longitude by altitude array in the vertical over a
2002	C specified number of levels using a three point smoother.
2003	C
2004	C NOTE: input array DATA is modified on output!
2005	C
2006	C Aug. 3/95 - C. McLandress
2007	C
2008	C Output arguement:
2009	C
2010	C * DATA = smoothed array (on output).
2011	C
2012	C Input arguements:
2013	C
2014	C * DATA = unsmoothed array of data (on input).
2015	C * WORK = work array of same dimension as DATA.
2016	C * COEFF = smoothing coefficient for a 1:COEFF:1 stencil.
2017	C * (e.g., COEFF = 2 will result in a smoother which
2018	C * weights the level L gridpoint by two and the two
2019	C * adjecent levels (L+1 and L-1) by one).
2020	C * NSMOOTH = number of times to smooth in vertical.
2021	C * (e.g., NSMOOTH=1 means smoothed only once,
2022	C * NSMOOTH=2 means smoothing repeated twice, etc.)
2023	C * IL1 = first longitudinal index to use (IL1 >= 1).
2024	C * IL2 = last longitudinal index to use (IL1 <= IL2 <= NLONS).
2025	C * LEV1 = first altitude level to use (LEV1 >=1).
2026	C * LEV2 = last altitude level to use (LEV1 < LEV2 <= NLEVS).
2027	C * NLONS = number of longitudes.
2028	C * NLEVS = number of vertical levels.
2029	C
2030	C Subroutine arguements.
2031	C
2032	INTEGER NSMOOTH, IL1, IL2, LEV1, LEV2, NLONS, NLEVS
2033	REAL COEFF
2034	REAL DATA(NLONS,NLEVS), WORK(NLONS,NLEVS)
2035	C
2036	C Internal variables.
2037	C
2038	INTEGER I, L, NS, LEV1P, LEV2M
2039	REAL SUM_WTS
2040	C-----------------------------------------------------------------------
2041	C
2042	C Calculate sum of weights.
2043	C
2044	SUM_WTS = COEFF + 2.
2045	C
2046	LEV1P = LEV1 + 1
2047	LEV2M = LEV2 - 1
2048	C
2049	C Smooth NSMOOTH times
2050	C
2051	DO 50 NS = 1,NSMOOTH
2052	C
2053	C Copy data into work array.
2054	C
2055	DO 20 L = LEV1,LEV2
2056	DO 10 I = IL1,IL2
2057	WORK(I,L) = DATA(I,L)
2058	10 CONTINUE
2059	20 CONTINUE
2060	C
2061	C Smooth array WORK in vertical direction and put into DATA.
2062	C
2063	DO 40 L = LEV1P,LEV2M
2064	DO 30 I = IL1,IL2
2065	DATA(I,L) = ( WORK(I,L+1) + COEFF*WORK(I,L) + WORK(I,L-1) )
2066	& / SUM_WTS
2067	30 CONTINUE
2068	40 CONTINUE
2069	C
2070	50 CONTINUE
2071	C
2072	RETURN
2073	C-----------------------------------------------------------------------
2074	END
2075
2076
2077
2078
2079
2080

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

Download in other formats:

Original Format