Context Navigation

hines_gwd.F @ 1255

Last change on this file since 1255 was 1237, checked in by Laurent Fairhead, 15 years ago
Des modifications sur la lecture des aerosols par Michael Correction du test sur le jour de lecture des aerosols qui ne marchait pas avec le nouveau calendrier (a revoir?) Menage sur quelques prints SD/MAF
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 1237 2009-09-03 12:03:33Z emillour $
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. / FLOAT(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