Context Navigation

hines_gwd.F @ 1176

Last change on this file since 1176 was 1001, checked in by Laurent Fairhead, 16 years ago

Modifs sur le parallelisme: masquage dans la physique
Inclusion strato
mise en coherence etat0
le mode offline fonctionne maintenant en parallele,
les fichiers de la dynamiques sont correctement sortis et peuvent etre reconstruit avec rebuild
la version parallele de la dynamique peut s'executer sans MPI (sur 1 proc)
L'OPENMP fonctionne maintenant sans la parallelisation MPI.

YM
LF

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

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

Download in other formats:

Original Format