source: LMDZ4/branches/LMDZ4-dev-20091210/libf/phylmd/hines_gwd.F @ 4999

Last change on this file since 4999 was 1237, checked in by Laurent Fairhead, 15 years ago

Des modifications sur la lecture des aerosols par Michael
Correction du test sur le jour de lecture des aerosols qui ne marchait
pas avec le nouveau calendrier (a revoir?)
Menage sur quelques prints
SD/MAF

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 69.5 KB
Line 
1!
2! $Id: hines_gwd.F 1237 2009-09-03 12:03:33Z lebasn $
3!
4      SUBROUTINE HINES_GWD(NLON,NLEV,DTIME,paphm1x, papm1x,
5     I      rlat,tx,ux,vx,
6     O      zustrhi,zvstrhi,
7     O      d_t_hin, d_u_hin, d_v_hin)
8
9C ########################################################################
10C Parametrization of the momentum flux deposition due to a broad band
11C spectrum of gravity waves, following Hines (1997a,b), as coded by
12C McLANDRESS (1995). Modified by McFARLANE and MANZINI (1995-1997)
13C                 MAECHAM model stand alone version
14C ########################################################################
15C
16C
17         USE dimphy
18         implicit none
19
20cym#include "dimensions.h"
21cym#include "dimphy.h"
22#include "YOEGWD.h"
23#include "YOMCST.h"
24
25      INTEGER NAZMTH
26      PARAMETER(NAZMTH=8)
27
28C     INPUT ARGUMENTS.
29C     ----- ----------
30C
31C  - 2D
32C  PAPHM1   : HALF LEVEL PRESSURE (T-DT)
33C  PAPM1    : FULL LEVEL PRESSURE (T-DT)
34C  PTM1     : TEMPERATURE (T-DT)
35C  PUM1     : ZONAL WIND (T-DT)
36C  PVM1     : MERIDIONAL WIND (T-DT)
37C
38
39C     REFERENCE.
40C     ----------
41C         SEE MODEL DOCUMENTATION
42C
43C     AUTHOR.
44C     -------
45C
46C      N. MCFARLANE   DKRZ-HAMBURG   MAY 1995
47C      STAND ALONE E. MANZINI MPI-HAMBURG FEBRUARY 1997
48C
49C      BASED ON A COMBINATION OF THE OROGRAPHIC SCHEME BY N.MCFARLANE 1987
50C      AND THE HINES SCHEME AS CODED BY C. MCLANDRESS 1995.                       
51C
52C
53C
54cym      INTEGER KLEVM1
55C
56      REAL PAPHM1(klon,klev+1), PAPM1(klon,KLEV) 
57      REAL PTM1(klon,KLEV), PUM1(klon,KLEV), PVM1(klon,KLEV)
58      REAL PRFLUX(klon)
59C1
60C1
61C1
62      REAL RLAT(klon),COSLAT(KLON)
63C
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
69C     * VERTICAL POSITIONING ARRAYS.
70
71      REAL SGJ(klon,KLEV),     SHJ(klon,KLEV),   
72     1     SHXKJ(klon,KLEV),   DSGJ(klon,KLEV)
73
74C     * LOGICAL SWITCHES TO CONTROL ROOF DRAG, ENVELOP GW DRAG AND
75C     * HINES' DOPPLER SPREADING EXTROWAVE GW DRAG.
76C     * LOZPR IS TRUE FOR ZPR ENHANCEMENT
77
78
79C     * 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     
97C     * THES ARE THE INPUT PARAMETERS FOR HINES ROUTINE AND
98C     * ARE SPECIFIED IN ROUTINE HINES_SETUP. SINCE THIS IS CALLED
99C     * ONLY AT FIRST CALL TO THIS ROUTINE THESE VARIABLES MUST BE SAVED
100C     * FOR USE AT SUBSEQUENT CALLS. THIS CAN BE AVOIDED BY CALLING
101C     * HINES_SETUP IN MAIN PROGRAM AND PASSING THE PARAMETERS AS
102C     * SUBROUTINE ARGUEMENTS.
103C
104
105      REAL    RMSCON
106      INTEGER NMESSG, IPRINT, ILRMS
107      INTEGER IFL
108C
109      INTEGER  NAZ,ICUTOFF,NSMAX,IHEATCAL
110      REAL  SLOPE,F1,F2,F3,F5,F6,KSTAR(KLON),ALT_CUTOFF,SMCO
111C
112C    PROVIDED AS INPUT
113C
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)
119c
120c     VARIABLES FOR OUTPUT
121c
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
126C
127C     * LOGICAL SWITCHES TO CONTROL PRECIP ENHANCEMENT AND
128C     * HINES' DOPPLER SPREADING EXTROWAVE GW DRAG.
129C     * LOZPR IS TRUE FOR ZPR ENHANCEMENT
130
131      LOGICAL LOZPR, LORMS(klon)
132C
133C  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
137C
138C  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   
144C
145C     
146C     PRINT *,' IT IS STARTED HINES GOING ON...'
147C
148C
149C
150C
151C*    COMPUTATIONAL CONSTANTS.
152C     ------------- ----------
153C
154C
155      d_t_hin(:,:)=0.
156     
157      RHOH2O=1000.   
158      ZPCONS = (1000.*86400.)/RHOH2O
159cym      KLEVM1=KLEV-1
160C
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
173C
174  100 CONTINUE
175C
176C    Define constants and arrays needed for the ccc/mam gwd scheme
177C    *Constants:
178
179      RGOCP=RD/RCPD
180      LREFP=KLEV-1
181      LREF=KLEV-2
182C1
183C1    *Arrays
184C1
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     
195CC
196      DO 211 JL=KIDIA,KFDIA
197      PRESSG(JL)=PAPHM1(JL,klev+1)
198  211 CONTINUE
199C
200C
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
207C
208C
209  400 CONTINUE
210
211C
212C
213C
214*/#########################################################################
215*/
216*/
217C
218C     * AUG. 14/95 - C. MCLANDRESS.
219C     * SEP.    95   N. MCFARLANE.
220C
221C     * THIS ROUTINE CALCULATES THE HORIZONTAL WIND TENDENCIES
222C     * DUE TO MCFARLANE'S OROGRAPHIC GW DRAG SCHEME, HINES'
223C     * DOPPLER SPREAD SCHEME FOR "EXTROWAVES" AND ADDS ON
224C     * ROOF DRAG. IT IS BASED ON THE ROUTINE GWDFLX8.
225C
226C     * LREFP IS THE INDEX OF THE MODEL LEVEL BELOW THE REFERENCE LEVEL
227C     * I/O ARRAYS PASSED FROM MAIN.
228C     * (PRESSG = SURFACE PRESSURE)
229C
230C
231C
232C
233C     * CONSTANTS VALUES DEFINED IN DATA STATEMENT ARE :
234C     * VMIN     = MIMINUM WIND IN THE DIRECTION OF REFERENCE LEVEL
235C     *            WIND BEFORE WE CONSIDER BREAKING TO HAVE OCCURED.
236C     * DMPSCAL  = DAMPING TIME FOR GW DRAG IN SECONDS.
237C     * TAUFAC   = 1/(LENGTH SCALE).
238C     * HMIN     = MIMINUM ENVELOPE HEIGHT REQUIRED TO PRODUCE GW DRAG.
239C     * 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
247C     * HINES EXTROWAVE GWD CONSTANTS DEFINED IN DATA STATEMENT ARE:
248C     * RMSCON = ROOT MEAN SQUARE GRAVITY WAVE WIND AT LOWEST LEVEL (M/S).
249C     * NMESSG  = UNIT NUMBER FOR PRINTED MESSAGES.
250C     * IPRINT  = 1 TO DO PRINT OUT SOME HINES ARRAYS.
251C     * IFL     = FIRST CALL FLAG TO HINES_SETUP ("SAVE" IT)
252C     * PCRIT = CRITICAL VALUE OF ZPR (MM/D)
253C     * IPLEV = LEVEL OF APPLICATION OF PRCIT
254C     * PCONS = FACTOR OF ZPR ENHANCEMENT
255C
256
257      DATA PCRIT / 5. /, PCONS / 4.75 /
258
259      IPLEV = LREFP-1
260C
261      DATA    RMSCON  / 1.00 /
262     1        IPRINT   /  2  /, NMESSG  /   6   /
263      DATA    IFL / 0 /
264C
265      LOZPR = .FALSE.
266C
267C-----------------------------------------------------------------------
268C
269C
270C
271C     * SET ERROR FLAG
272
273      IERROR = 0
274
275C     * SPECIFY VARIOUS PARAMETERS FOR HINES ROUTINE AT VERY FIRST CALL.
276C     * (NOTE THAT ARRAY K_ALPHA IS SPECIFIED SO MAKE SURE THAT
277C     * IT IS NOT OVERWRITTEN LATER ON).
278C
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
283C
284C     * START GWD CALCULATIONS.
285
286      LREF  = LREFP-1
287
288C
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
294c
295
296
297C     * INITIALIZE NECESSARY ARRAYS.
298C
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
308C
309C     * IF USING HINES SCHEME THEN CALCULATE B V FREQUENCY AT ALL POINTS
310C     * 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(-DTTDSF*SGJ(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
331C
332C
333  300 CONTINUE
334
335C     * CALCULATE GW DRAG DUE TO HINES' EXTROWAVES
336C     * SET MOLECULAR VISCOSITY TO A VERY SMALL VALUE.
337C     * IF THE MODEL TOP IS GREATER THAN 100 KM THEN THE ACTUAL
338C     * 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
351C     * 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
359C     * 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
368C
369C     * INITIALIZE SWITCHES FOR HINES GWD CALCULATION
370C
371      ILRMS = 0
372C
373      DO 345 I=KIDIA,KFDIA
374      LORMS(I) = .FALSE.
375  345 CONTINUE
376C
377C
378C     * DEFILE BOTTOM LAUNCH LEVEL
379C
380      LEVBOT = IPLEV
381C
382C     * BACKGROUND WIND MINUS VALUE AT BOTTOM LAUNCH LEVEL.
383C
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
389C
390C     * SPECIFY ROOT MEAN SQUARE WIND AT BOTTOM LAUNCH LEVEL.
391C
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
404C
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
411C
412C     * CALCULATE GWD (NOTE THAT DIFFUSION COEFFICIENT AND
413C     * HEATING RATE ONLY CALCULATED IF IHEATCAL = 1).
414C
415      IF ( ILRMS.GT.0 )       THEN                   
416C
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
429C     * ADD ON HINES' GWD TENDENCIES TO OROGRAPHIC TENDENCIES AND
430C     * 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
437C
438
439C     * END OF HINES CALCULATIONS.
440C
441      ENDIF
442C
443  500 CONTINUE
444
445
446C-----------------------------------------------------------------------
447C
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
458c     PRINT *,'UTENDGW:',UTENDGW
459
460C     PRINT *,' HINES HAS BEEN COMPLETED (LONG ISNT IT...)'
461
462      RETURN
463 999  CONTINUE
464
465C     * 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
472C
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
492C
493C  Main routine for Hines' "extrowave" gravity wave parameterization based
494C  on Hines' Doppler spread theory. This routine calculates zonal
495C  and meridional components of gravity wave drag, heating rates
496C  and diffusion coefficient on a longitude by altitude grid.
497C  No "mythical" lower boundary region calculation is made so it
498C  is assumed that lowest level winds are weak (i.e, approximately zero).
499C
500C  Aug. 13/95 - C. McLandress
501C  SEPT. /95  - N.McFarlane
502C
503C  Modifications:
504C
505C  Output arguements:
506C
507C     * DRAG_U = zonal component of gravity wave drag (m/s^2).
508C     * DRAG_V = meridional component of gravity wave drag (m/s^2).
509C     * HEAT   = gravity wave heating (K/sec).
510C     * DIFFCO = diffusion coefficient (m^2/sec)
511C     * FLUX_U = zonal component of vertical momentum flux (Pascals)
512C     * FLUX_V = meridional component of vertical momentum flux (Pascals)
513C
514C  Input arguements:
515C
516C     * VEL_U      = background zonal wind component (m/s).
517C     * VEL_V      = background meridional wind component (m/s).
518C     * BVFREQ     = background Brunt Vassala frequency (radians/sec).
519C     * DENSITY    = background density (kg/m^3)
520C     * VISC_MOL   = molecular viscosity (m^2/s)
521C     * ALT        = altitude of momentum, density, buoyancy levels (m)
522C     *              (NOTE: levels ordered so that ALT(I,1) > ALT(I,2), etc.)
523C     * RMSWIND   = root mean square gravity wave wind at lowest level (m/s).
524C     * K_ALPHA    = horizontal wavenumber of each azimuth (1/m).
525C     * IORDER     = 1 means vertical levels are indexed from top down
526C     *              (i.e., highest level indexed 1 and lowest level NLEVS);
527C     *           .NE. 1 highest level is index NLEVS.
528C     * IHEATCAL   = 1 to calculate heating rates and diffusion coefficient.
529C     * IPRINT     = 1 to print out various arrays.
530C     * ICUTOFF    = 1 to exponentially damp GWD, heating and diffusion
531C     *              arrays above ALT_CUTOFF; otherwise arrays not modified.
532C     * ALT_CUTOFF = altitude in meters above which exponential decay applied.
533C     * SMCO       = smoothing factor used to smooth cutoff vertical
534C     *              wavenumbers and total rms winds in vertical direction
535C     *              before calculating drag or heating
536C     *              (SMCO >= 1 ==> 1:SMCO:1 stencil used).
537C     * NSMAX      = number of times smoother applied ( >= 1),
538C     *            = 0 means no smoothing performed.
539C     * KSTAR      = typical gravity wave horizontal wavenumber (1/m).
540C     * SLOPE      = slope of incident vertical wavenumber spectrum
541C     *              (SLOPE must equal 1., 1.5 or 2.).
542C     * F1 to F6   = Hines's fudge factors (F4 not needed since used for
543C     *              vertical flux of vertical momentum).
544C     * NAZ        = actual number of horizontal azimuths used.
545C     * IL1        = first longitudinal index to use (IL1 >= 1).
546C     * IL2        = last longitudinal index to use (IL1 <= IL2 <= NLONS).
547C     * LEV1       = index of first level for drag calculation.
548C     * LEV2       = index of last level for drag calculation
549C     *              (i.e., LEV1 < LEV2 <= NLEVS).
550C     * NLONS      = number of longitudes.
551C     * NLEVS      = number of vertical levels.
552C     * NAZMTH     = azimuthal array dimension (NAZMTH >= NAZ).
553C
554C  Work arrays.
555C
556C     * M_ALPHA      = cutoff vertical wavenumber (1/m).
557C     * V_ALPHA      = wind component at each azimuth (m/s) and if IHEATCAL=1
558C     *                holds vertical derivative of cutoff wavenumber.
559C     * SIGMA_ALPHA  = total rms wind in each azimuth (m/s).
560C     * SIGSQH_ALPHA = portion of wind variance from waves having wave
561C     *                normals in the alpha azimuth (m/s).
562C     * SIGMA_T      = total rms horizontal wind (m/s).
563C     * AK_ALPHA     = spectral amplitude factor at each azimuth
564C     *                (i.e.,{AjKj}) in m^4/s^2.
565C     * I_ALPHA      = Hines' integral.
566C     * MMIN_ALPHA   = minimum value of cutoff wavenumber.
567C     * DENSB        = background density at bottom level.
568C     * BVFB         = buoyancy frequency at bottom level and
569C     *                work array for ICUTOFF = 1.
570C
571C     * LORMS       = .TRUE. for drag computation
572C
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)
593C
594      LOGICAL LORMS(NLONS)
595C
596C  Internal variables.
597C
598      INTEGER  LEVBOT, LEVTOP, I, N, L, LEV1P, LEV2M
599      INTEGER  ILPRT1, ILPRT2
600C-----------------------------------------------------------------------
601C
602C     PRINT *,' IN HINES_EXTRO0'
603      LEV1P = LEV1 + 1
604      LEV2M = LEV2 - 1
605C
606C  Index of lowest altitude level (bottom of drag calculation).
607C
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
614C
615C  Buoyancy and density at bottom level.
616C
617      DO 10 I = IL1,IL2
618        BVFB(I)  = BVFREQ(I,LEVBOT)
619        DENSB(I) = DENSITY(I,LEVBOT)
620 10   CONTINUE
621C
622C  initialize some variables
623C
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
637C
638C  Compute azimuthal wind components from zonal and meridional winds.
639C
640      CALL HINES_WIND ( V_ALPHA,
641     ^                  VEL_U, VEL_V, NAZ,
642     ^                  IL1, IL2, LEV1, LEV2, NLONS, NLEVS, NAZMTH )
643C
644C  Calculate cutoff vertical wavenumber and velocity variances.
645C
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)
652C
653C  Smooth cutoff wavenumbers and total rms velocity in the vertical
654C  direction NSMAX times, using FLUX_U as temporary work array.
655C   
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
674C
675C  Calculate zonal and meridional components of the
676C  momentum flux and drag.
677C
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)
683C
684C  Cutoff drag above ALT_CUTOFF, using BVFB as temporary work array.
685C
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   
694C
695C  Print out various arrays for diagnostic purposes.
696C
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
705C
706C  If not calculating heating rate and diffusion coefficient then finished.
707C
708      IF (IHEATCAL.NE.1)  RETURN
709C
710C  Calculate vertical derivative of cutoff wavenumber (store
711C  in array V_ALPHA) using centered differences at interior gridpoints
712C  and one-sided differences at first and last levels.
713C
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
730C
731C  Heating rate and diffusion coefficient.
732C
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)
738C
739C  Finished.
740C
741      RETURN
742C-----------------------------------------------------------------------
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)
751C
752C  This routine calculates the cutoff vertical wavenumber and velocity
753C  variances on a longitude by altitude grid for the Hines' Doppler
754C  spread gravity wave drag parameterization scheme.
755C  NOTE: (1) only values of four or eight can be used for # azimuths (NAZ).
756C        (2) only values of 1.0, 1.5 or 2.0 can be used for slope (SLOPE).
757C
758C  Aug. 10/95 - C. McLandress
759C
760C  Output arguements:
761C
762C     * M_ALPHA      = cutoff wavenumber at each azimuth (1/m).
763C     * SIGMA_ALPHA  = total rms wind in each azimuth (m/s).
764C     * SIGSQH_ALPHA = portion of wind variance from waves having wave
765C     *                normals in the alpha azimuth (m/s).
766C     * SIGMA_T      = total rms horizontal wind (m/s).
767C     * AK_ALPHA     = spectral amplitude factor at each azimuth
768C     *                (i.e.,{AjKj}) in m^4/s^2.
769C
770C  Input arguements:
771C
772C     * V_ALPHA  = wind component at each azimuth (m/s).
773C     * VISC_MOL = molecular viscosity (m^2/s)
774C     * DENSITY  = background density (kg/m^3).
775C     * DENSB    = background density at model bottom (kg/m^3).
776C     * BVFREQ   = background Brunt Vassala frequency (radians/sec).
777C     * BVFB     = background Brunt Vassala frequency at model bottom.
778C     * RMS_WIND = root mean square gravity wave wind at lowest level (m/s).
779C     * KSTAR    = typical gravity wave horizontal wavenumber (1/m).
780C     * SLOPE    = slope of incident vertical wavenumber spectrum
781C     *            (SLOPE = 1., 1.5 or 2.).
782C     * F1,F2,F3 = Hines's fudge factors.
783C     * NAZ      = actual number of horizontal azimuths used (4 or 8).
784C     * LEVBOT   = index of lowest vertical level.
785C     * LEVTOP   = index of highest vertical level
786C     *            (NOTE: if LEVTOP < LEVBOT then level index
787C     *             increases from top down).
788C     * IL1      = first longitudinal index to use (IL1 >= 1).
789C     * IL2      = last longitudinal index to use (IL1 <= IL2 <= NLONS).
790C     * NLONS    = number of longitudes.
791C     * NLEVS    = number of vertical levels.
792C     * NAZMTH   = azimuthal array dimension (NAZMTH >= NAZ).
793C
794C     * LORMS       = .TRUE. for drag computation
795C
796C  Input work arrays:
797C
798C     * I_ALPHA    = Hines' integral at a single level.
799C     * MMIN_ALPHA = minimum value of cutoff wavenumber.
800C
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)
817C
818      LOGICAL LORMS(NLONS)
819C
820C Internal variables.
821C
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
826cc      REAL  N_OVER_M(1000), SIGFAC(1000)
827
828      REAL  N_OVER_M(NLONS), SIGFAC(NLONS)
829      DATA  VISC_MIN / 1.E-10 /
830C-----------------------------------------------------------------------     
831C
832
833C     PRINT *,'IN HINES_WAVNUM'
834      SP1 = SLOPE + 1.
835C
836C  Indices of levels to process.
837C
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
846C
847C  Use horizontal isotropy to calculate azimuthal variances at bottom level.
848C
849      AZFAC = 1. / FLOAT(NAZ)
850      DO 20 N = 1,NAZ
851        DO 10 I = IL1,IL2
852          SIGSQH_ALPHA(I,LEVBOT,N) = AZFAC * RMS_WIND(I)**2
853 10     CONTINUE
854 20   CONTINUE
855C
856C  Velocity variances at bottom level.
857C
858      CALL HINES_SIGMA ( SIGMA_T, SIGMA_ALPHA,
859     ^                   SIGSQH_ALPHA, NAZ, LEVBOT,
860     ^                   IL1, IL2, NLONS, NLEVS, NAZMTH)
861c
862      CALL HINES_SIGMA ( SIGMATM, SIGALPMC,
863     ^                   SIGSQMCW, NAZ, LEVBOT,
864     ^                   IL1, IL2, NLONS, NLEVS, NAZMTH)
865C
866C  Calculate cutoff wavenumber and spectral amplitude factor
867C  at bottom level where it is assumed that background winds vanish
868C  and also initialize minimum value of cutoff wavnumber.
869C
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
882C
883C  Calculate quantities from the bottom upwards,
884C  starting one level above bottom.
885C
886      DO 150 L = LSTART,LEND,LINCR
887C
888C  Level beneath present level.
889C
890        LBELOW = L - LINCR
891C
892C  Calculate N/m_M where m_M is maximum permissible value of the vertical
893C  wavenumber (i.e., m > m_M are obliterated) and N is buoyancy frequency.
894C  m_M is taken as the smaller of the instability-induced
895C  wavenumber (M_SUB_M_TURB) and that imposed by molecular viscosity
896C  (M_SUB_M_MOL). Since variance at this level is not yet known
897C  use value at level below.
898C
899        DO 50 I = IL1,IL2
900        IF (LORMS(I)) THEN
901c
902        F2MFAC=SIGMATM(I,LBELOW)**2
903        F2MOD(I,LBELOW) =1.+ 2.*F2MFAC
904     ^                      / ( F2MFAC+SIGMA_T(I,LBELOW)**2 )
905c
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
917C
918C  Calculate cutoff wavenumber at this level.
919C
920        DO 70 N = 1,NAZ
921          DO 60 I = IL1,IL2
922          IF (LORMS(I)) THEN
923C
924C  Calculate trial value (since variance at this level is not yet known
925C  use value at level below). If trial value is negative or if it exceeds
926C  minimum value (not permitted) then set it to minimum value.
927C                                                                     
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
934C
935C  Reset minimum value of cutoff wavenumber if necessary.
936C
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
940C
941          ENDIF
942 60       CONTINUE
943 70     CONTINUE
944C
945C  Calculate the Hines integral at this level.
946C
947        CALL HINES_INTGRL ( I_ALPHA,
948     ^                      V_ALPHA, M_ALPHA, BVFB, SLOPE, NAZ,
949     ^                      L, IL1, IL2, NLONS, NLEVS, NAZMTH,
950     ^                      LORMS )
951
952C
953C  Calculate the velocity variances at this level.
954C
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 )
968c
969        CALL HINES_SIGMA ( SIGMATM, SIGALPMC,
970     ^                     SIGSQMCW, NAZ, L,
971     ^                     IL1, IL2, NLONS, NLEVS, NAZMTH )
972C
973C  End of level loop.
974C
975 150  CONTINUE
976C
977      RETURN
978C-----------------------------------------------------------------------
979      END
980
981      SUBROUTINE HINES_WIND (V_ALPHA,VEL_U,VEL_V,
982     1                       NAZ,IL1,IL2,LEV1,LEV2,NLONS,NLEVS,NAZMTH)
983C
984C  This routine calculates the azimuthal horizontal background wind components
985C  on a longitude by altitude grid for the case of 4 or 8 azimuths for
986C  the Hines' Doppler spread GWD parameterization scheme.
987C
988C  Aug. 7/95 - C. McLandress
989C
990C  Output arguement:
991C
992C     * V_ALPHA   = background wind component at each azimuth (m/s).
993C     *             (note: first azimuth is in eastward direction
994C     *              and rotate in counterclockwise direction.)
995C
996C  Input arguements:
997C
998C     * VEL_U     = background zonal wind component (m/s).
999C     * VEL_V     = background meridional wind component (m/s).
1000C     * NAZ       = actual number of horizontal azimuths used (must be 4 or 8).
1001C     * IL1       = first longitudinal index to use (IL1 >= 1).
1002C     * IL2       = last longitudinal index to use (IL1 <= IL2 <= NLONS).
1003C     * LEV1      = first altitude level to use (LEV1 >=1).
1004C     * LEV2      = last altitude level to use (LEV1 < LEV2 <= NLEVS).
1005C     * NLONS     = number of longitudes.
1006C     * NLEVS     = number of vertical levels.
1007C     * NAZMTH    = azimuthal array dimension (NAZMTH >= NAZ).
1008C
1009C  Constants in DATA statements.
1010C
1011C     * COS45 = cosine of 45 degrees.           
1012C     * UMIN  = minimum allowable value for zonal or meridional
1013C     *         wind component (m/s).
1014C
1015C  Subroutine arguements.
1016C
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)
1021C
1022C  Internal variables.
1023C
1024      INTEGER  I, L
1025      REAL U, V, COS45, UMIN
1026C
1027      DATA  COS45 / 0.7071068 /
1028      DATA  UMIN / 0.001 /
1029C-----------------------------------------------------------------------     
1030C
1031C  Case with 4 azimuths.
1032C
1033
1034C      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
1049C
1050C  Case with 8 azimuths.
1051C
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
1070C
1071      RETURN
1072C-----------------------------------------------------------------------
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)
1079C
1080C  Calculate zonal and meridional components of the vertical flux
1081C  of horizontal momentum and corresponding wave drag (force per unit mass)
1082C  on a longitude by altitude grid for the Hines' Doppler spread
1083C  GWD parameterization scheme.
1084C  NOTE: only 4 or 8 azimuths can be used.
1085C
1086C  Aug. 6/95 - C. McLandress
1087C
1088C  Output arguements:
1089C
1090C     * FLUX_U = zonal component of vertical momentum flux (Pascals)
1091C     * FLUX_V = meridional component of vertical momentum flux (Pascals)
1092C     * DRAG_U = zonal component of drag (m/s^2).
1093C     * DRAG_V = meridional component of drag (m/s^2).
1094C
1095C  Input arguements:
1096C
1097C     * ALT       = altitudes (m).
1098C     * DENSITY   = background density (kg/m^3).
1099C     * DENSB     = background density at bottom level (kg/m^3).
1100C     * M_ALPHA   = cutoff vertical wavenumber (1/m).
1101C     * AK_ALPHA  = spectral amplitude factor (i.e., {AjKj} in m^4/s^2).
1102C     * K_ALPHA   = horizontal wavenumber (1/m).
1103C     * SLOPE     = slope of incident vertical wavenumber spectrum.
1104C     * NAZ       = actual number of horizontal azimuths used (must be 4 or 8).
1105C     * IL1       = first longitudinal index to use (IL1 >= 1).
1106C     * IL2       = last longitudinal index to use (IL1 <= IL2 <= NLONS).
1107C     * LEV1      = first altitude level to use (LEV1 >=1).
1108C     * LEV2      = last altitude level to use (LEV1 < LEV2 <= NLEVS).
1109C     * NLONS     = number of longitudes.
1110C     * NLEVS     = number of vertical levels.
1111C     * NAZMTH    = azimuthal array dimension (NAZMTH >= NAZ).
1112C
1113C     * LORMS       = .TRUE. for drag computation
1114C
1115C  Constant in DATA statement.
1116C
1117C     * COS45 = cosine of 45 degrees.           
1118C
1119C  Subroutine arguements.
1120C
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)
1129C
1130      LOGICAL LORMS(NLONS)
1131C
1132C  Internal variables.
1133C
1134      INTEGER  I, L, LEV1P, LEV2M
1135      REAL  COS45, PROD2, PROD4, PROD6, PROD8, DENDZ, DENDZ2
1136      DATA  COS45 / 0.7071068 /   
1137C-----------------------------------------------------------------------
1138C
1139      LEV1P = LEV1 + 1
1140      LEV2M = LEV2 - 1
1141      LEV2P = LEV2 + 1
1142C
1143C  Sum over azimuths for case where SLOPE = 1.
1144C
1145      IF (SLOPE.EQ.1.)  THEN
1146C
1147C  Case with 4 azimuths.
1148C
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
1159C
1160C  Case with 8 azimuths.
1161C
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
1180C
1181      END IF
1182C
1183C  Sum over azimuths for case where SLOPE not equal to 1.
1184C
1185      IF (SLOPE.NE.1.)  THEN
1186C
1187C  Case with 4 azimuths.
1188C
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
1201C
1202C  Case with 8 azimuths.
1203C
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
1222C
1223      END IF
1224C
1225C  Calculate flux from sum.
1226C
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
1233C
1234C  Calculate drag at intermediate levels using centered differences
1235C     
1236      DO 120 L = LEV1P,LEV2M
1237        DO 110 I = IL1,IL2
1238        IF (LORMS(I)) THEN
1239ccc       DENDZ2 = DENSITY(I,L) * ( ALT(I,L+1) - ALT(I,L-1) )
1240          DENDZ2 = DENSITY(I,L) * ( ALT(I,L-1) - ALT(I,L) )
1241ccc       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
1243ccc       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
1249C
1250C  Drag at first and last levels using one-side differences.
1251C
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
1273150   CONTINUE
1274      ENDIF
1275C
1276      RETURN
1277C-----------------------------------------------------------------------
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)
1284C
1285C  This routine calculates the gravity wave induced heating and
1286C  diffusion coefficient on a longitude by altitude grid for 
1287C  the Hines' Doppler spread gravity wave drag parameterization scheme.
1288C
1289C  Aug. 6/95 - C. McLandress
1290C
1291C  Output arguements:
1292C
1293C     * HEAT   = gravity wave heating (K/sec).
1294C     * DIFFCO = diffusion coefficient (m^2/sec)
1295C
1296C  Input arguements:
1297C
1298C     * M_ALPHA     = cutoff vertical wavenumber (1/m).
1299C     * DMDZ_ALPHA  = vertical derivative of cutoff wavenumber.
1300C     * AK_ALPHA    = spectral amplitude factor of each azimuth
1301C                     (i.e., {AjKj} in m^4/s^2).
1302C     * K_ALPHA     = horizontal wavenumber of each azimuth (1/m).
1303C     * BVFREQ      = background Brunt Vassala frequency (rad/sec).
1304C     * DENSITY     = background density (kg/m^3).
1305C     * DENSB       = background density at bottom level (kg/m^3).
1306C     * SIGMA_T     = total rms horizontal wind (m/s).
1307C     * VISC_MOL    = molecular viscosity (m^2/s).
1308C     * KSTAR       = typical gravity wave horizontal wavenumber (1/m).
1309C     * SLOPE       = slope of incident vertical wavenumber spectrum.
1310C     * F2,F3,F5,F6 = Hines's fudge factors.
1311C     * NAZ         = actual number of horizontal azimuths used.
1312C     * IL1         = first longitudinal index to use (IL1 >= 1).
1313C     * IL2         = last longitudinal index to use (IL1 <= IL2 <= NLONS).
1314C     * LEV1        = first altitude level to use (LEV1 >=1).
1315C     * LEV2        = last altitude level to use (LEV1 < LEV2 <= NLEVS).
1316C     * NLONS       = number of longitudes.
1317C     * NLEVS       = number of vertical levels.
1318C     * NAZMTH      = azimuthal array dimension (NAZMTH >= NAZ).
1319C
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)
1327C
1328C Internal variables.
1329C
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
1333C
1334C specific heat at constant pressure
1335C
1336      DATA  CPGAS / 1004. /
1337C             
1338C minimum permissible viscosity
1339C
1340      DATA  VISC_MIN / 1.E-10 /       
1341C-----------------------------------------------------------------------     
1342C
1343C  Initialize heating array.
1344C
1345      DO 20 L = 1,NLEVS
1346        DO 10 I = 1,NLONS
1347          HEAT(I,L) = 0.
1348  10    CONTINUE
1349  20  CONTINUE
1350C
1351C  Perform sum over azimuths for case where SLOPE = 1.
1352C
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
1363C
1364C  Perform sum over azimuths for case where SLOPE not 1.
1365C
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
1377C
1378C  Heating and diffusion.
1379C
1380      DO 100 L = LEV1,LEV2
1381        DO 90 I = IL1,IL2
1382C
1383C  Maximum permissible value of cutoff wavenumber is the smaller
1384C  of the instability-induced wavenumber (M_SUB_M_TURB) and
1385C  that imposed by molecular viscosity (M_SUB_M_MOL).
1386C
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 )
1391C
1392          HEATNG = - HEAT(I,L) * F5 * BVFREQ(I,L) / M_SUB_M
1393     ^               * DENSB(I) / DENSITY(I,L)
1394          DIFFCO(I,L) = F6 * HEATNG**0.33333333 / M_SUB_M**1.33333333
1395          HEAT(I,L)   = HEATNG / CPGAS
1396C
1397 90     CONTINUE
1398 100  CONTINUE
1399C
1400      RETURN
1401C-----------------------------------------------------------------------
1402      END
1403
1404      SUBROUTINE HINES_SIGMA (SIGMA_T,SIGMA_ALPHA,SIGSQH_ALPHA,
1405     1                        NAZ,LEV,IL1,IL2,NLONS,NLEVS,NAZMTH)
1406C
1407C  This routine calculates the total rms and azimuthal rms horizontal
1408C  velocities at a given level on a longitude by altitude grid for
1409C  the Hines' Doppler spread GWD parameterization scheme.
1410C  NOTE: only four or eight azimuths can be used.
1411C
1412C  Aug. 7/95 - C. McLandress
1413C
1414C  Output arguements:
1415C
1416C     * SIGMA_T      = total rms horizontal wind (m/s).
1417C     * SIGMA_ALPHA  = total rms wind in each azimuth (m/s).
1418C
1419C  Input arguements:
1420C
1421C     * SIGSQH_ALPHA = portion of wind variance from waves having wave
1422C     *                normals in the alpha azimuth (m/s).
1423C     * NAZ       = actual number of horizontal azimuths used (must be 4 or 8).
1424C     * LEV       = altitude level to process.
1425C     * IL1       = first longitudinal index to use (IL1 >= 1).
1426C     * IL2       = last longitudinal index to use (IL1 <= IL2 <= NLONS).
1427C     * NLONS     = number of longitudes.
1428C     * NLEVS     = number of vertical levels.
1429C     * NAZMTH    = azimuthal array dimension (NAZMTH >= NAZ).
1430C
1431C  Subroutine arguements.
1432C
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)
1438C
1439C  Internal variables.
1440C
1441      INTEGER  I, N
1442      REAL  SUM_EVEN, SUM_ODD
1443C-----------------------------------------------------------------------     
1444C
1445C  Calculate azimuthal rms velocity for the 4 azimuth case.
1446C
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
1457C
1458C  Calculate azimuthal rms velocity for the 8 azimuth case.
1459C
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
1484C
1485C  Calculate total rms velocity.
1486C
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
1498C
1499      RETURN
1500C-----------------------------------------------------------------------     
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)
1506C
1507C  This routine calculates the vertical wavenumber integral
1508C  for a single vertical level at each azimuth on a longitude grid
1509C  for the Hines' Doppler spread GWD parameterization scheme.
1510C  NOTE: (1) only spectral slopes of 1, 1.5 or 2 are permitted.
1511C        (2) the integral is written in terms of the product QM
1512C            which by construction is always less than 1. Series
1513C            solutions are used for small |QM| and analytical solutions
1514C            for remaining values.
1515C
1516C  Aug. 8/95 - C. McLandress
1517C
1518C  Output arguement:
1519C
1520C     * I_ALPHA = Hines' integral.
1521C
1522C  Input arguements:
1523C
1524C     * V_ALPHA = azimuthal wind component (m/s).
1525C     * M_ALPHA = azimuthal cutoff vertical wavenumber (1/m).
1526C     * BVFB    = background Brunt Vassala frequency at model bottom.
1527C     * SLOPE   = slope of initial vertical wavenumber spectrum
1528C     *           (must use SLOPE = 1., 1.5 or 2.)
1529C     * NAZ     = actual number of horizontal azimuths used.
1530C     * LEV     = altitude level to process.
1531C     * IL1     = first longitudinal index to use (IL1 >= 1).
1532C     * IL2     = last longitudinal index to use (IL1 <= IL2 <= NLONS).
1533C     * NLONS   = number of longitudes.
1534C     * NLEVS   = number of vertical levels.
1535C     * NAZMTH  = azimuthal array dimension (NAZMTH >= NAZ).
1536C
1537C     * LORMS       = .TRUE. for drag computation
1538C
1539C  Constants in DATA statements:
1540C
1541C     * QMIN = minimum value of Q_ALPHA (avoids indeterminant form of integral)
1542C     * QM_MIN = minimum value of Q_ALPHA * M_ALPHA (used to avoid numerical
1543C     *          problems).
1544C
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
1550C
1551      LOGICAL LORMS(NLONS)
1552C
1553C  Internal variables.
1554C
1555      INTEGER  I, N
1556      REAL  Q_ALPHA, QM, SQRTQM, Q_MIN, QM_MIN
1557C
1558      DATA  Q_MIN / 1.0 /, QM_MIN / 0.01 /
1559C-----------------------------------------------------------------------     
1560C
1561C  For integer value SLOPE = 1.
1562C
1563      IF (SLOPE .EQ. 1.)  THEN
1564C
1565        DO 20 N = 1,NAZ
1566          DO 10 I = IL1,IL2
1567          IF (LORMS(I)) THEN
1568C
1569            Q_ALPHA = V_ALPHA(I,LEV,N) / BVFB(I)
1570            QM = Q_ALPHA * M_ALPHA(I,LEV,N)
1571C
1572C  If |QM| is small then use first 4 terms series of Taylor series
1573C  expansion of integral in order to avoid indeterminate form of integral,
1574C  otherwise use analytical form of integral.
1575C
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) = ( QM**2/2. + QM**3/3. + QM**4/4.
1581     ^                           + QM**5/5. ) / Q_ALPHA**2
1582              END IF
1583            ELSE
1584              I_ALPHA(I,N) = - ( ALOG(1.-QM) + QM ) / Q_ALPHA**2
1585            END IF
1586C
1587          ENDIF
1588 10       CONTINUE
1589 20     CONTINUE
1590C
1591      END IF
1592C
1593C  For integer value SLOPE = 2.
1594C
1595      IF (SLOPE .EQ. 2.)  THEN
1596C
1597        DO 40 N = 1,NAZ
1598          DO 30 I = IL1,IL2
1599          IF (LORMS(I)) THEN
1600C
1601            Q_ALPHA = V_ALPHA(I,LEV,N) / BVFB(I)
1602            QM = Q_ALPHA * M_ALPHA(I,LEV,N)
1603C
1604C  If |QM| is small then use first 4 terms series of Taylor series
1605C  expansion of integral in order to avoid indeterminate form of integral,
1606C  otherwise use analytical form of integral.
1607C
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) = ( QM**3/3. + QM**4/4. + QM**5/5.
1613     ^                           + QM**6/6. ) / Q_ALPHA**3
1614              END IF
1615            ELSE
1616              I_ALPHA(I,N) = - ( ALOG(1.-QM) + QM + QM**2/2.)
1617     ^                         / Q_ALPHA**3
1618            END IF
1619C
1620          ENDIF
1621 30       CONTINUE
1622 40     CONTINUE
1623C
1624      END IF
1625C
1626C  For real value SLOPE = 1.5
1627C
1628      IF (SLOPE .EQ. 1.5)  THEN
1629C
1630        DO 60 N = 1,NAZ
1631          DO 50 I = IL1,IL2
1632          IF (LORMS(I)) THEN
1633C
1634            Q_ALPHA = V_ALPHA(I,LEV,N) / BVFB(I)
1635            QM = Q_ALPHA * M_ALPHA(I,LEV,N)       
1636C
1637C  If |QM| is small then use first 4 terms series of Taylor series
1638C  expansion of integral in order to avoid indeterminate form of integral,
1639C  otherwise use analytical form of integral.
1640C
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     ^                           + QM**3/4.5 + QM**4/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
1660C
1661          ENDIF
1662 50       CONTINUE
1663 60     CONTINUE
1664C
1665      END IF
1666C
1667C  If integral is negative (which in principal should not happen) then
1668C  print a message and some info since execution will abort when calculating
1669C  the variances.
1670C
1671c      DO 80 N = 1,NAZ
1672c        DO 70 I = IL1,IL2
1673c          IF (I_ALPHA(I,N).LT.0.)  THEN
1674c            WRITE (6,*)
1675c            WRITE (6,*) '******************************'
1676c            WRITE (6,*) 'Hines integral I_ALPHA < 0 '
1677c            WRITE (6,*) '  longitude I=',I
1678c            WRITE (6,*) '  azimuth   N=',N
1679c            WRITE (6,*) '  level   LEV=',LEV
1680c            WRITE (6,*) '  I_ALPHA =',I_ALPHA(I,N)
1681c            WRITE (6,*) '  V_ALPHA =',V_ALPHA(I,LEV,N)
1682c            WRITE (6,*) '  M_ALPHA =',M_ALPHA(I,LEV,N)
1683c            WRITE (6,*) '  Q_ALPHA =',V_ALPHA(I,LEV,N) / BVFB(I)
1684c            WRITE (6,*) '  QM      =',V_ALPHA(I,LEV,N) / BVFB(I)
1685c     ^                                * M_ALPHA(I,LEV,N)
1686c            WRITE (6,*) '******************************'
1687c          END IF
1688c 70     CONTINUE
1689c 80   CONTINUE
1690C
1691      RETURN
1692C-----------------------------------------------------------------------
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)
1698C
1699C  This routine specifies various parameters needed for the
1700C  the Hines' Doppler spread gravity wave drag parameterization scheme.
1701C
1702C  Aug. 8/95 - C. McLandress
1703C
1704C  Output arguements:
1705C
1706C     * NAZ        = actual number of horizontal azimuths used
1707C     *              (code set up presently for only NAZ = 4 or 8).
1708C     * SLOPE      = slope of incident vertical wavenumber spectrum
1709C     *              (code set up presently for SLOPE 1., 1.5 or 2.).
1710C     * F1         = "fudge factor" used in calculation of trial value of
1711C     *              azimuthal cutoff wavenumber M_ALPHA (1.2 <= F1 <= 1.9).
1712C     * F2         = "fudge factor" used in calculation of maximum
1713C     *              permissible instabiliy-induced cutoff wavenumber
1714C     *              M_SUB_M_TURB (0.1 <= F2 <= 1.4).
1715C     * F3         = "fudge factor" used in calculation of maximum
1716C     *              permissible molecular viscosity-induced cutoff wavenumber
1717C     *              M_SUB_M_MOL (0.1 <= F2 <= 1.4).
1718C     * F5         = "fudge factor" used in calculation of heating rate
1719C     *              (1 <= F5 <= 3).
1720C     * F6         = "fudge factor" used in calculation of turbulent
1721C     *              diffusivity coefficient.
1722C     * KSTAR      = typical gravity wave horizontal wavenumber (1/m)
1723C     *              used in calculation of M_SUB_M_TURB.
1724C     * ICUTOFF    = 1 to exponentially damp off GWD, heating and diffusion
1725C     *              arrays above ALT_CUTOFF; otherwise arrays not modified.
1726C     * ALT_CUTOFF = altitude in meters above which exponential decay applied.
1727C     * SMCO       = smoother used to smooth cutoff vertical wavenumbers
1728C     *              and total rms winds before calculating drag or heating.
1729C     *              (==> a 1:SMCO:1 stencil used; SMCO >= 1.).
1730C     * NSMAX      = number of times smoother applied ( >= 1),
1731C     *            = 0 means no smoothing performed.
1732C     * IHEATCAL   = 1 to calculate heating rates and diffusion coefficient.
1733C     *            = 0 means only drag and flux calculated.
1734C     * K_ALPHA    = horizontal wavenumber of each azimuth (1/m) which
1735C     *              is set here to KSTAR.
1736C     * IERROR     = error flag.
1737C     *            = 0 no errors.
1738C     *            = 10 ==> NAZ > NAZMTH
1739C     *            = 20 ==> invalid number of azimuths (NAZ must be 4 or 8).
1740C     *            = 30 ==> invalid slope (SLOPE must be 1., 1.5 or 2.).
1741C     *            = 40 ==> invalid smoother (SMCO must be >= 1.)
1742C
1743C  Input arguements:
1744C
1745C     * NMESSG  = output unit number where messages to be printed.
1746C     * NLONS   = number of longitudes.
1747C     * NAZMTH  = azimuthal array dimension (NAZMTH >= NAZ).
1748C
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
1754C
1755C Internal variables.
1756C
1757      INTEGER  I, N
1758C-----------------------------------------------------------------------     
1759C
1760C  Specify constants.
1761C
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
1777c      SMCO       = 1.0
1778      NSMAX      = 5
1779c      NSMAX      = 2
1780      IHEATCAL   = 0
1781C
1782C  Print information to output file.
1783C
1784c      WRITE (NMESSG,6000)
1785c 6000 FORMAT (/' Subroutine HINES_SETUP:')
1786c      WRITE (NMESSG,*)  '  SLOPE = ', SLOPE
1787c      WRITE (NMESSG,*)  '  NAZ = ', NAZ
1788c      WRITE (NMESSG,*)  '  F1,F2,F3  = ', F1, F2, F3
1789c      WRITE (NMESSG,*)  '  F5,F6     = ', F5, F6
1790c      WRITE (NMESSG,*)  '  KSTAR     = ', KSTAR
1791c     >           ,'  COSLAT     = ', COSLAT
1792c      IF (ICUTOFF .EQ. 1)  THEN
1793c        WRITE (NMESSG,*) '  Drag exponentially damped above ',
1794c     &                       ALT_CUTOFF/1.E3
1795c     END IF
1796c      IF (NSMAX.LT.1 )  THEN
1797c        WRITE (NMESSG,*) '  No smoothing of cutoff wavenumbers, etc'
1798c      ELSE
1799c        WRITE (NMESSG,*) '  Cutoff wavenumbers and sig_t smoothed:'
1800c        WRITE (NMESSG,*) '    SMCO  =', SMCO
1801c        WRITE (NMESSG,*) '    NSMAX =', NSMAX
1802c     END IF
1803C
1804C  Check that things are setup correctly and log error if not
1805C
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
1811C
1812C  Use single value for azimuthal-dependent horizontal wavenumber.
1813C
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
1819C
1820      RETURN
1821C-----------------------------------------------------------------------
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)
1829C
1830C  Print out altitude profiles of various quantities from
1831C  Hines' Doppler spread gravity wave drag parameterization scheme.
1832C  (NOTE: only for NAZ = 4 or 8).
1833C
1834C  Aug. 8/95 - C. McLandress
1835C
1836C  Input arguements:
1837C
1838C     * IU_PRINT = 1 to print out values in east-west direction.
1839C     * IV_PRINT = 1 to print out values in north-south direction.
1840C     * NMESSG   = unit number for printed output.
1841C     * ILPRT1   = first longitudinal index to print.
1842C     * ILPRT2   = last longitudinal index to print.
1843C     * LEVPRT1  = first altitude level to print.
1844C     * LEVPRT2  = last altitude level to print.
1845C
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)
1854C
1855C  Internal variables.
1856C
1857      INTEGER  N_EAST, N_WEST, N_NORTH, N_SOUTH
1858      INTEGER  I, L
1859C-----------------------------------------------------------------------
1860C
1861C  Azimuthal indices of cardinal directions.
1862C
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
1873C
1874C  Print out values for range of longitudes.
1875C
1876      DO 100 I = ILPRT1,ILPRT2
1877C
1878C  Print east-west wind, sigmas, cutoff wavenumbers, flux and drag.
1879C
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
1896C
1897C  Print north-south winds, sigmas, cutoff wavenumbers, flux and drag.
1898C
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
1914C
1915 100  CONTINUE
1916C
1917      RETURN
1918C-----------------------------------------------------------------------
1919      END
1920
1921      SUBROUTINE HINES_EXP (DATA,DATA_ZMAX,ALT,ALT_EXP,IORDER,
1922     1                      IL1,IL2,LEV1,LEV2,NLONS,NLEVS)
1923C
1924C  This routine exponentially damps a longitude by altitude array
1925C  of data above a specified altitude.
1926C
1927C  Aug. 13/95 - C. McLandress
1928C
1929C  Output arguements:
1930C
1931C     * DATA = modified data array.
1932C
1933C  Input arguements:
1934C
1935C     * DATA    = original data array.
1936C     * ALT     = altitudes.
1937C     * ALT_EXP = altitude above which exponential decay applied.
1938C     * IORDER  = 1 means vertical levels are indexed from top down
1939C     *           (i.e., highest level indexed 1 and lowest level NLEVS);
1940C     *           .NE. 1 highest level is index NLEVS.
1941C     * IL1     = first longitudinal index to use (IL1 >= 1).
1942C     * IL2     = last longitudinal index to use (IL1 <= IL2 <= NLONS).
1943C     * LEV1    = first altitude level to use (LEV1 >=1).
1944C     * LEV2    = last altitude level to use (LEV1 < LEV2 <= NLEVS).
1945C     * NLONS   = number of longitudes.
1946C     * NLEVS   = number of vertical
1947C
1948C  Input work arrays:
1949C
1950C     * DATA_ZMAX = data values just above altitude ALT_EXP.
1951C
1952      INTEGER  IORDER, IL1, IL2, LEV1, LEV2, NLONS, NLEVS
1953      REAL  ALT_EXP
1954      REAL  DATA(NLONS,NLEVS), DATA_ZMAX(NLONS), ALT(NLONS,NLEVS)
1955C
1956C Internal variables.
1957C
1958      INTEGER  LEVBOT, LEVTOP, LINCR, I, L
1959      REAL  HSCALE
1960      DATA  HSCALE / 5.E3 /
1961C-----------------------------------------------------------------------     
1962C
1963C  Index of lowest altitude level (bottom of drag calculation).
1964C
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
1973C
1974C  Data values at first level above ALT_EXP.
1975C
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
1983C
1984C  Exponentially damp field above ALT_EXP to model top at L=1.
1985C
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
1993C
1994      RETURN
1995C-----------------------------------------------------------------------
1996      END
1997
1998      SUBROUTINE VERT_SMOOTH (DATA,WORK,COEFF,NSMOOTH,
1999     1                        IL1,IL2,LEV1,LEV2,NLONS,NLEVS)
2000C
2001C  Smooth a longitude by altitude array in the vertical over a
2002C  specified number of levels using a three point smoother.
2003C
2004C  NOTE: input array DATA is modified on output!
2005C
2006C  Aug. 3/95 - C. McLandress
2007C
2008C  Output arguement:
2009C
2010C     * DATA    = smoothed array (on output).
2011C
2012C  Input arguements:
2013C
2014C     * DATA    = unsmoothed array of data (on input).
2015C     * WORK    = work array of same dimension as DATA.
2016C     * COEFF   = smoothing coefficient for a 1:COEFF:1 stencil.
2017C     *           (e.g., COEFF = 2 will result in a smoother which
2018C     *           weights the level L gridpoint by two and the two
2019C     *           adjecent levels (L+1 and L-1) by one).
2020C     * NSMOOTH = number of times to smooth in vertical.
2021C     *           (e.g., NSMOOTH=1 means smoothed only once,
2022C     *           NSMOOTH=2 means smoothing repeated twice, etc.)
2023C     * IL1     = first longitudinal index to use (IL1 >= 1).
2024C     * IL2     = last longitudinal index to use (IL1 <= IL2 <= NLONS).
2025C     * LEV1    = first altitude level to use (LEV1 >=1).
2026C     * LEV2    = last altitude level to use (LEV1 < LEV2 <= NLEVS).
2027C     * NLONS   = number of longitudes.
2028C     * NLEVS   = number of vertical levels.
2029C
2030C  Subroutine arguements.
2031C
2032      INTEGER  NSMOOTH, IL1, IL2, LEV1, LEV2, NLONS, NLEVS
2033      REAL  COEFF
2034      REAL  DATA(NLONS,NLEVS), WORK(NLONS,NLEVS)
2035C
2036C  Internal variables.
2037C
2038      INTEGER  I, L, NS, LEV1P, LEV2M
2039      REAL  SUM_WTS
2040C-----------------------------------------------------------------------     
2041C
2042C  Calculate sum of weights.
2043C
2044      SUM_WTS = COEFF + 2.
2045C
2046      LEV1P = LEV1 + 1
2047      LEV2M = LEV2 - 1
2048C
2049C  Smooth NSMOOTH times
2050C
2051      DO 50 NS = 1,NSMOOTH
2052C
2053C  Copy data into work array.
2054C
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
2060C
2061C  Smooth array WORK in vertical direction and put into DATA.
2062C
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
2069C
2070 50   CONTINUE
2071C
2072      RETURN
2073C-----------------------------------------------------------------------
2074      END
2075
2076
2077
2078
2079     
2080     
Note: See TracBrowser for help on using the repository browser.