Context Navigation

hines_gwd.F @ 1105

Last change on this file since 1105 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

Rev	Line
[1001]	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