Context Navigation

hines_gwd.F @ 2261

Last change on this file since 2261 was 1403, checked in by Laurent Fairhead, 14 years ago

Merged LMDZ4V5.0-dev branch changes r1292:r1399 to trunk.

Validation:
Validation consisted in compiling the HEAD revision of the trunk,
LMDZ4V5.0-dev branch and the merged sources and running different
configurations on local and SX8 machines comparing results.

Local machine: bench configuration, 32x24x11, gfortran

IPSLCM5A configuration (comparison between trunk and merged sources):
- numerical convergence on dynamical fields over 3 days
- start files are equivalent (except for RN and PB fields)
- daily history files equivalent

MH07 configuration, new physics package (comparison between LMDZ4V5.0-dev branch and merged sources):
- numerical convergence on dynamical fields over 3 days
- start files are equivalent (except for RN and PB fields)
- daily history files equivalent

SX8 machine (brodie), 96x95x39 on 4 processors:

IPSLCM5A configuration:
- start files are equivalent (except for RN and PB fields)
- monthly history files equivalent

MH07 configuration:
- start files are equivalent (except for RN and PB fields)
- monthly history files equivalent

Changes to the makegcm and create_make_gcm scripts to take into account
main programs in F90 files

Fusion de la branche LMDZ4V5.0-dev (r1292:r1399) au tronc principal

Validation:
La validation a consisté à compiler la HEAD de le trunk et de la banche
LMDZ4V5.0-dev et les sources fusionnées et de faire tourner le modéle selon
différentes configurations en local et sur SX8 et de comparer les résultats

En local: 32x24x11, config bench/gfortran

pour une config IPSLCM5A (comparaison tronc/fusion):
- convergence numérique sur les champs dynamiques après 3 jours
- restart et restartphy égaux (à part sur RN et Pb)
- fichiers histoire égaux

pour une config nlle physique (MH07) (comparaison LMDZ4v5.0-dev/fusion):
- convergence numérique sur les champs dynamiques après 3 jours
- restart et restartphy égaux
- fichiers histoire équivalents

Sur brodie, 96x95x39 sur 4 proc:

pour une config IPSLCM5A:
- restart et restartphy égaux (à part sur RN et PB)
- pas de différence dans les fichiers histmth.nc

pour une config MH07
- restart et restartphy égaux (à part sur RN et PB)
- pas de différence dans les fichiers histmth.nc

Changement sur makegcm et create_make-gcm pour pouvoir prendre en compte des
programmes principaux en *F90

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

Rev	Line
[1279]	1	!
	2	! $Id: hines_gwd.F 1403 2010-07-01 09:02:53Z fhourdin $
	3	!
[1001]	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
[1403]	849	AZFAC = 1. / REAL(NAZ)
[1001]	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
[1279]	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
[1001]	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