source: LMDZ4/trunk/libf/phylmd/hines_gwd.F @ 1176

Last change on this file since 1176 was 1001, checked in by Laurent Fairhead, 16 years ago
  • Modifs sur le parallelisme: masquage dans la physique
  • Inclusion strato
  • mise en coherence etat0
  • le mode offline fonctionne maintenant en parallele,
  • les fichiers de la dynamiques sont correctement sortis et peuvent etre reconstruit avec rebuild
  • la version parallele de la dynamique peut s'executer sans MPI (sur 1 proc)
  • L'OPENMP fonctionne maintenant sans la parallelisation MPI.

YM
LF

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 69.5 KB
Line 
1      SUBROUTINE HINES_GWD(NLON,NLEV,DTIME,paphm1x, papm1x,
2     I      rlat,tx,ux,vx,
3     O      zustrhi,zvstrhi,
4     O      d_t_hin, d_u_hin, d_v_hin)
5
6C ########################################################################
7C Parametrization of the momentum flux deposition due to a broad band
8C spectrum of gravity waves, following Hines (1997a,b), as coded by
9C McLANDRESS (1995). Modified by McFARLANE and MANZINI (1995-1997)
10C                 MAECHAM model stand alone version
11C ########################################################################
12C
13C
14         USE dimphy
15         implicit none
16
17cym#include "dimensions.h"
18cym#include "dimphy.h"
19#include "YOEGWD.h"
20#include "YOMCST.h"
21
22      INTEGER NAZMTH
23      PARAMETER(NAZMTH=8)
24
25C     INPUT ARGUMENTS.
26C     ----- ----------
27C
28C  - 2D
29C  PAPHM1   : HALF LEVEL PRESSURE (T-DT)
30C  PAPM1    : FULL LEVEL PRESSURE (T-DT)
31C  PTM1     : TEMPERATURE (T-DT)
32C  PUM1     : ZONAL WIND (T-DT)
33C  PVM1     : MERIDIONAL WIND (T-DT)
34C
35
36C     REFERENCE.
37C     ----------
38C         SEE MODEL DOCUMENTATION
39C
40C     AUTHOR.
41C     -------
42C
43C      N. MCFARLANE   DKRZ-HAMBURG   MAY 1995
44C      STAND ALONE E. MANZINI MPI-HAMBURG FEBRUARY 1997
45C
46C      BASED ON A COMBINATION OF THE OROGRAPHIC SCHEME BY N.MCFARLANE 1987
47C      AND THE HINES SCHEME AS CODED BY C. MCLANDRESS 1995.                       
48C
49C
50C
51cym      INTEGER KLEVM1
52C
53      REAL PAPHM1(klon,klev+1), PAPM1(klon,KLEV) 
54      REAL PTM1(klon,KLEV), PUM1(klon,KLEV), PVM1(klon,KLEV)
55      REAL PRFLUX(klon)
56C1
57C1
58C1
59      REAL RLAT(klon),COSLAT(KLON)
60C
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
66C     * VERTICAL POSITIONING ARRAYS.
67
68      REAL SGJ(klon,KLEV),     SHJ(klon,KLEV),   
69     1     SHXKJ(klon,KLEV),   DSGJ(klon,KLEV)
70
71C     * LOGICAL SWITCHES TO CONTROL ROOF DRAG, ENVELOP GW DRAG AND
72C     * HINES' DOPPLER SPREADING EXTROWAVE GW DRAG.
73C     * LOZPR IS TRUE FOR ZPR ENHANCEMENT
74
75
76C     * 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     
94C     * THES ARE THE INPUT PARAMETERS FOR HINES ROUTINE AND
95C     * ARE SPECIFIED IN ROUTINE HINES_SETUP. SINCE THIS IS CALLED
96C     * ONLY AT FIRST CALL TO THIS ROUTINE THESE VARIABLES MUST BE SAVED
97C     * FOR USE AT SUBSEQUENT CALLS. THIS CAN BE AVOIDED BY CALLING
98C     * HINES_SETUP IN MAIN PROGRAM AND PASSING THE PARAMETERS AS
99C     * SUBROUTINE ARGUEMENTS.
100C
101
102      REAL    RMSCON
103      INTEGER NMESSG, IPRINT, ILRMS
104      INTEGER IFL
105C
106      INTEGER  NAZ,ICUTOFF,NSMAX,IHEATCAL
107      REAL  SLOPE,F1,F2,F3,F5,F6,KSTAR(KLON),ALT_CUTOFF,SMCO
108C
109C    PROVIDED AS INPUT
110C
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)
116c
117c     VARIABLES FOR OUTPUT
118c
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
123C
124C     * LOGICAL SWITCHES TO CONTROL PRECIP ENHANCEMENT AND
125C     * HINES' DOPPLER SPREADING EXTROWAVE GW DRAG.
126C     * LOZPR IS TRUE FOR ZPR ENHANCEMENT
127
128      LOGICAL LOZPR, LORMS(klon)
129C
130C  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
134C
135C  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   
141C
142C     
143C     PRINT *,' IT IS STARTED HINES GOING ON...'
144C
145C
146C
147C
148C*    COMPUTATIONAL CONSTANTS.
149C     ------------- ----------
150C
151C
152      d_t_hin(:,:)=0.
153     
154      RHOH2O=1000.   
155      ZPCONS = (1000.*86400.)/RHOH2O
156cym      KLEVM1=KLEV-1
157C
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
170C
171  100 CONTINUE
172C
173C    Define constants and arrays needed for the ccc/mam gwd scheme
174C    *Constants:
175
176      RGOCP=RD/RCPD
177      LREFP=KLEV-1
178      LREF=KLEV-2
179C1
180C1    *Arrays
181C1
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     
192CC
193      DO 211 JL=KIDIA,KFDIA
194      PRESSG(JL)=PAPHM1(JL,klev+1)
195  211 CONTINUE
196C
197C
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
204C
205C
206  400 CONTINUE
207
208C
209C
210C
211*/#########################################################################
212*/
213*/
214C
215C     * AUG. 14/95 - C. MCLANDRESS.
216C     * SEP.    95   N. MCFARLANE.
217C
218C     * THIS ROUTINE CALCULATES THE HORIZONTAL WIND TENDENCIES
219C     * DUE TO MCFARLANE'S OROGRAPHIC GW DRAG SCHEME, HINES'
220C     * DOPPLER SPREAD SCHEME FOR "EXTROWAVES" AND ADDS ON
221C     * ROOF DRAG. IT IS BASED ON THE ROUTINE GWDFLX8.
222C
223C     * LREFP IS THE INDEX OF THE MODEL LEVEL BELOW THE REFERENCE LEVEL
224C     * I/O ARRAYS PASSED FROM MAIN.
225C     * (PRESSG = SURFACE PRESSURE)
226C
227C
228C
229C
230C     * CONSTANTS VALUES DEFINED IN DATA STATEMENT ARE :
231C     * VMIN     = MIMINUM WIND IN THE DIRECTION OF REFERENCE LEVEL
232C     *            WIND BEFORE WE CONSIDER BREAKING TO HAVE OCCURED.
233C     * DMPSCAL  = DAMPING TIME FOR GW DRAG IN SECONDS.
234C     * TAUFAC   = 1/(LENGTH SCALE).
235C     * HMIN     = MIMINUM ENVELOPE HEIGHT REQUIRED TO PRODUCE GW DRAG.
236C     * 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
244C     * HINES EXTROWAVE GWD CONSTANTS DEFINED IN DATA STATEMENT ARE:
245C     * RMSCON = ROOT MEAN SQUARE GRAVITY WAVE WIND AT LOWEST LEVEL (M/S).
246C     * NMESSG  = UNIT NUMBER FOR PRINTED MESSAGES.
247C     * IPRINT  = 1 TO DO PRINT OUT SOME HINES ARRAYS.
248C     * IFL     = FIRST CALL FLAG TO HINES_SETUP ("SAVE" IT)
249C     * PCRIT = CRITICAL VALUE OF ZPR (MM/D)
250C     * IPLEV = LEVEL OF APPLICATION OF PRCIT
251C     * PCONS = FACTOR OF ZPR ENHANCEMENT
252C
253
254      DATA PCRIT / 5. /, PCONS / 4.75 /
255
256      IPLEV = LREFP-1
257C
258      DATA    RMSCON  / 1.00 /
259     1        IPRINT   /  2  /, NMESSG  /   6   /
260      DATA    IFL / 0 /
261C
262      LOZPR = .FALSE.
263C
264C-----------------------------------------------------------------------
265C
266C
267C
268C     * SET ERROR FLAG
269
270      IERROR = 0
271
272C     * SPECIFY VARIOUS PARAMETERS FOR HINES ROUTINE AT VERY FIRST CALL.
273C     * (NOTE THAT ARRAY K_ALPHA IS SPECIFIED SO MAKE SURE THAT
274C     * IT IS NOT OVERWRITTEN LATER ON).
275C
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
280C
281C     * START GWD CALCULATIONS.
282
283      LREF  = LREFP-1
284
285C
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
291c
292
293
294C     * INITIALIZE NECESSARY ARRAYS.
295C
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
305C
306C     * IF USING HINES SCHEME THEN CALCULATE B V FREQUENCY AT ALL POINTS
307C     * 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(-DTTDSF*SGJ(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
328C
329C
330  300 CONTINUE
331
332C     * CALCULATE GW DRAG DUE TO HINES' EXTROWAVES
333C     * SET MOLECULAR VISCOSITY TO A VERY SMALL VALUE.
334C     * IF THE MODEL TOP IS GREATER THAN 100 KM THEN THE ACTUAL
335C     * 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
348C     * 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
356C     * 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
365C
366C     * INITIALIZE SWITCHES FOR HINES GWD CALCULATION
367C
368      ILRMS = 0
369C
370      DO 345 I=KIDIA,KFDIA
371      LORMS(I) = .FALSE.
372  345 CONTINUE
373C
374C
375C     * DEFILE BOTTOM LAUNCH LEVEL
376C
377      LEVBOT = IPLEV
378C
379C     * BACKGROUND WIND MINUS VALUE AT BOTTOM LAUNCH LEVEL.
380C
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
386C
387C     * SPECIFY ROOT MEAN SQUARE WIND AT BOTTOM LAUNCH LEVEL.
388C
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
401C
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
408C
409C     * CALCULATE GWD (NOTE THAT DIFFUSION COEFFICIENT AND
410C     * HEATING RATE ONLY CALCULATED IF IHEATCAL = 1).
411C
412      IF ( ILRMS.GT.0 )       THEN                   
413C
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
426C     * ADD ON HINES' GWD TENDENCIES TO OROGRAPHIC TENDENCIES AND
427C     * 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
434C
435
436C     * END OF HINES CALCULATIONS.
437C
438      ENDIF
439C
440  500 CONTINUE
441
442
443C-----------------------------------------------------------------------
444C
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
455c     PRINT *,'UTENDGW:',UTENDGW
456
457C     PRINT *,' HINES HAS BEEN COMPLETED (LONG ISNT IT...)'
458
459      RETURN
460 999  CONTINUE
461
462C     * 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
469C
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
489C
490C  Main routine for Hines' "extrowave" gravity wave parameterization based
491C  on Hines' Doppler spread theory. This routine calculates zonal
492C  and meridional components of gravity wave drag, heating rates
493C  and diffusion coefficient on a longitude by altitude grid.
494C  No "mythical" lower boundary region calculation is made so it
495C  is assumed that lowest level winds are weak (i.e, approximately zero).
496C
497C  Aug. 13/95 - C. McLandress
498C  SEPT. /95  - N.McFarlane
499C
500C  Modifications:
501C
502C  Output arguements:
503C
504C     * DRAG_U = zonal component of gravity wave drag (m/s^2).
505C     * DRAG_V = meridional component of gravity wave drag (m/s^2).
506C     * HEAT   = gravity wave heating (K/sec).
507C     * DIFFCO = diffusion coefficient (m^2/sec)
508C     * FLUX_U = zonal component of vertical momentum flux (Pascals)
509C     * FLUX_V = meridional component of vertical momentum flux (Pascals)
510C
511C  Input arguements:
512C
513C     * VEL_U      = background zonal wind component (m/s).
514C     * VEL_V      = background meridional wind component (m/s).
515C     * BVFREQ     = background Brunt Vassala frequency (radians/sec).
516C     * DENSITY    = background density (kg/m^3)
517C     * VISC_MOL   = molecular viscosity (m^2/s)
518C     * ALT        = altitude of momentum, density, buoyancy levels (m)
519C     *              (NOTE: levels ordered so that ALT(I,1) > ALT(I,2), etc.)
520C     * RMSWIND   = root mean square gravity wave wind at lowest level (m/s).
521C     * K_ALPHA    = horizontal wavenumber of each azimuth (1/m).
522C     * IORDER     = 1 means vertical levels are indexed from top down
523C     *              (i.e., highest level indexed 1 and lowest level NLEVS);
524C     *           .NE. 1 highest level is index NLEVS.
525C     * IHEATCAL   = 1 to calculate heating rates and diffusion coefficient.
526C     * IPRINT     = 1 to print out various arrays.
527C     * ICUTOFF    = 1 to exponentially damp GWD, heating and diffusion
528C     *              arrays above ALT_CUTOFF; otherwise arrays not modified.
529C     * ALT_CUTOFF = altitude in meters above which exponential decay applied.
530C     * SMCO       = smoothing factor used to smooth cutoff vertical
531C     *              wavenumbers and total rms winds in vertical direction
532C     *              before calculating drag or heating
533C     *              (SMCO >= 1 ==> 1:SMCO:1 stencil used).
534C     * NSMAX      = number of times smoother applied ( >= 1),
535C     *            = 0 means no smoothing performed.
536C     * KSTAR      = typical gravity wave horizontal wavenumber (1/m).
537C     * SLOPE      = slope of incident vertical wavenumber spectrum
538C     *              (SLOPE must equal 1., 1.5 or 2.).
539C     * F1 to F6   = Hines's fudge factors (F4 not needed since used for
540C     *              vertical flux of vertical momentum).
541C     * NAZ        = actual number of horizontal azimuths used.
542C     * IL1        = first longitudinal index to use (IL1 >= 1).
543C     * IL2        = last longitudinal index to use (IL1 <= IL2 <= NLONS).
544C     * LEV1       = index of first level for drag calculation.
545C     * LEV2       = index of last level for drag calculation
546C     *              (i.e., LEV1 < LEV2 <= NLEVS).
547C     * NLONS      = number of longitudes.
548C     * NLEVS      = number of vertical levels.
549C     * NAZMTH     = azimuthal array dimension (NAZMTH >= NAZ).
550C
551C  Work arrays.
552C
553C     * M_ALPHA      = cutoff vertical wavenumber (1/m).
554C     * V_ALPHA      = wind component at each azimuth (m/s) and if IHEATCAL=1
555C     *                holds vertical derivative of cutoff wavenumber.
556C     * SIGMA_ALPHA  = total rms wind in each azimuth (m/s).
557C     * SIGSQH_ALPHA = portion of wind variance from waves having wave
558C     *                normals in the alpha azimuth (m/s).
559C     * SIGMA_T      = total rms horizontal wind (m/s).
560C     * AK_ALPHA     = spectral amplitude factor at each azimuth
561C     *                (i.e.,{AjKj}) in m^4/s^2.
562C     * I_ALPHA      = Hines' integral.
563C     * MMIN_ALPHA   = minimum value of cutoff wavenumber.
564C     * DENSB        = background density at bottom level.
565C     * BVFB         = buoyancy frequency at bottom level and
566C     *                work array for ICUTOFF = 1.
567C
568C     * LORMS       = .TRUE. for drag computation
569C
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)
590C
591      LOGICAL LORMS(NLONS)
592C
593C  Internal variables.
594C
595      INTEGER  LEVBOT, LEVTOP, I, N, L, LEV1P, LEV2M
596      INTEGER  ILPRT1, ILPRT2
597C-----------------------------------------------------------------------
598C
599C     PRINT *,' IN HINES_EXTRO0'
600      LEV1P = LEV1 + 1
601      LEV2M = LEV2 - 1
602C
603C  Index of lowest altitude level (bottom of drag calculation).
604C
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
611C
612C  Buoyancy and density at bottom level.
613C
614      DO 10 I = IL1,IL2
615        BVFB(I)  = BVFREQ(I,LEVBOT)
616        DENSB(I) = DENSITY(I,LEVBOT)
617 10   CONTINUE
618C
619C  initialize some variables
620C
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
634C
635C  Compute azimuthal wind components from zonal and meridional winds.
636C
637      CALL HINES_WIND ( V_ALPHA,
638     ^                  VEL_U, VEL_V, NAZ,
639     ^                  IL1, IL2, LEV1, LEV2, NLONS, NLEVS, NAZMTH )
640C
641C  Calculate cutoff vertical wavenumber and velocity variances.
642C
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)
649C
650C  Smooth cutoff wavenumbers and total rms velocity in the vertical
651C  direction NSMAX times, using FLUX_U as temporary work array.
652C   
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
671C
672C  Calculate zonal and meridional components of the
673C  momentum flux and drag.
674C
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)
680C
681C  Cutoff drag above ALT_CUTOFF, using BVFB as temporary work array.
682C
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   
691C
692C  Print out various arrays for diagnostic purposes.
693C
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
702C
703C  If not calculating heating rate and diffusion coefficient then finished.
704C
705      IF (IHEATCAL.NE.1)  RETURN
706C
707C  Calculate vertical derivative of cutoff wavenumber (store
708C  in array V_ALPHA) using centered differences at interior gridpoints
709C  and one-sided differences at first and last levels.
710C
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
727C
728C  Heating rate and diffusion coefficient.
729C
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)
735C
736C  Finished.
737C
738      RETURN
739C-----------------------------------------------------------------------
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)
748C
749C  This routine calculates the cutoff vertical wavenumber and velocity
750C  variances on a longitude by altitude grid for the Hines' Doppler
751C  spread gravity wave drag parameterization scheme.
752C  NOTE: (1) only values of four or eight can be used for # azimuths (NAZ).
753C        (2) only values of 1.0, 1.5 or 2.0 can be used for slope (SLOPE).
754C
755C  Aug. 10/95 - C. McLandress
756C
757C  Output arguements:
758C
759C     * M_ALPHA      = cutoff wavenumber at each azimuth (1/m).
760C     * SIGMA_ALPHA  = total rms wind in each azimuth (m/s).
761C     * SIGSQH_ALPHA = portion of wind variance from waves having wave
762C     *                normals in the alpha azimuth (m/s).
763C     * SIGMA_T      = total rms horizontal wind (m/s).
764C     * AK_ALPHA     = spectral amplitude factor at each azimuth
765C     *                (i.e.,{AjKj}) in m^4/s^2.
766C
767C  Input arguements:
768C
769C     * V_ALPHA  = wind component at each azimuth (m/s).
770C     * VISC_MOL = molecular viscosity (m^2/s)
771C     * DENSITY  = background density (kg/m^3).
772C     * DENSB    = background density at model bottom (kg/m^3).
773C     * BVFREQ   = background Brunt Vassala frequency (radians/sec).
774C     * BVFB     = background Brunt Vassala frequency at model bottom.
775C     * RMS_WIND = root mean square gravity wave wind at lowest level (m/s).
776C     * KSTAR    = typical gravity wave horizontal wavenumber (1/m).
777C     * SLOPE    = slope of incident vertical wavenumber spectrum
778C     *            (SLOPE = 1., 1.5 or 2.).
779C     * F1,F2,F3 = Hines's fudge factors.
780C     * NAZ      = actual number of horizontal azimuths used (4 or 8).
781C     * LEVBOT   = index of lowest vertical level.
782C     * LEVTOP   = index of highest vertical level
783C     *            (NOTE: if LEVTOP < LEVBOT then level index
784C     *             increases from top down).
785C     * IL1      = first longitudinal index to use (IL1 >= 1).
786C     * IL2      = last longitudinal index to use (IL1 <= IL2 <= NLONS).
787C     * NLONS    = number of longitudes.
788C     * NLEVS    = number of vertical levels.
789C     * NAZMTH   = azimuthal array dimension (NAZMTH >= NAZ).
790C
791C     * LORMS       = .TRUE. for drag computation
792C
793C  Input work arrays:
794C
795C     * I_ALPHA    = Hines' integral at a single level.
796C     * MMIN_ALPHA = minimum value of cutoff wavenumber.
797C
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)
814C
815      LOGICAL LORMS(NLONS)
816C
817C Internal variables.
818C
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
823cc      REAL  N_OVER_M(1000), SIGFAC(1000)
824
825      REAL  N_OVER_M(NLONS), SIGFAC(NLONS)
826      DATA  VISC_MIN / 1.E-10 /
827C-----------------------------------------------------------------------     
828C
829
830C     PRINT *,'IN HINES_WAVNUM'
831      SP1 = SLOPE + 1.
832C
833C  Indices of levels to process.
834C
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
843C
844C  Use horizontal isotropy to calculate azimuthal variances at bottom level.
845C
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
852C
853C  Velocity variances at bottom level.
854C
855      CALL HINES_SIGMA ( SIGMA_T, SIGMA_ALPHA,
856     ^                   SIGSQH_ALPHA, NAZ, LEVBOT,
857     ^                   IL1, IL2, NLONS, NLEVS, NAZMTH)
858c
859      CALL HINES_SIGMA ( SIGMATM, SIGALPMC,
860     ^                   SIGSQMCW, NAZ, LEVBOT,
861     ^                   IL1, IL2, NLONS, NLEVS, NAZMTH)
862C
863C  Calculate cutoff wavenumber and spectral amplitude factor
864C  at bottom level where it is assumed that background winds vanish
865C  and also initialize minimum value of cutoff wavnumber.
866C
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
879C
880C  Calculate quantities from the bottom upwards,
881C  starting one level above bottom.
882C
883      DO 150 L = LSTART,LEND,LINCR
884C
885C  Level beneath present level.
886C
887        LBELOW = L - LINCR
888C
889C  Calculate N/m_M where m_M is maximum permissible value of the vertical
890C  wavenumber (i.e., m > m_M are obliterated) and N is buoyancy frequency.
891C  m_M is taken as the smaller of the instability-induced
892C  wavenumber (M_SUB_M_TURB) and that imposed by molecular viscosity
893C  (M_SUB_M_MOL). Since variance at this level is not yet known
894C  use value at level below.
895C
896        DO 50 I = IL1,IL2
897        IF (LORMS(I)) THEN
898c
899        F2MFAC=SIGMATM(I,LBELOW)**2
900        F2MOD(I,LBELOW) =1.+ 2.*F2MFAC
901     ^                      / ( F2MFAC+SIGMA_T(I,LBELOW)**2 )
902c
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
914C
915C  Calculate cutoff wavenumber at this level.
916C
917        DO 70 N = 1,NAZ
918          DO 60 I = IL1,IL2
919          IF (LORMS(I)) THEN
920C
921C  Calculate trial value (since variance at this level is not yet known
922C  use value at level below). If trial value is negative or if it exceeds
923C  minimum value (not permitted) then set it to minimum value.
924C                                                                     
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
931C
932C  Reset minimum value of cutoff wavenumber if necessary.
933C
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
937C
938          ENDIF
939 60       CONTINUE
940 70     CONTINUE
941C
942C  Calculate the Hines integral at this level.
943C
944        CALL HINES_INTGRL ( I_ALPHA,
945     ^                      V_ALPHA, M_ALPHA, BVFB, SLOPE, NAZ,
946     ^                      L, IL1, IL2, NLONS, NLEVS, NAZMTH,
947     ^                      LORMS )
948
949C
950C  Calculate the velocity variances at this level.
951C
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 )
965c
966        CALL HINES_SIGMA ( SIGMATM, SIGALPMC,
967     ^                     SIGSQMCW, NAZ, L,
968     ^                     IL1, IL2, NLONS, NLEVS, NAZMTH )
969C
970C  End of level loop.
971C
972 150  CONTINUE
973C
974      RETURN
975C-----------------------------------------------------------------------
976      END
977
978      SUBROUTINE HINES_WIND (V_ALPHA,VEL_U,VEL_V,
979     1                       NAZ,IL1,IL2,LEV1,LEV2,NLONS,NLEVS,NAZMTH)
980C
981C  This routine calculates the azimuthal horizontal background wind components
982C  on a longitude by altitude grid for the case of 4 or 8 azimuths for
983C  the Hines' Doppler spread GWD parameterization scheme.
984C
985C  Aug. 7/95 - C. McLandress
986C
987C  Output arguement:
988C
989C     * V_ALPHA   = background wind component at each azimuth (m/s).
990C     *             (note: first azimuth is in eastward direction
991C     *              and rotate in counterclockwise direction.)
992C
993C  Input arguements:
994C
995C     * VEL_U     = background zonal wind component (m/s).
996C     * VEL_V     = background meridional wind component (m/s).
997C     * NAZ       = actual number of horizontal azimuths used (must be 4 or 8).
998C     * IL1       = first longitudinal index to use (IL1 >= 1).
999C     * IL2       = last longitudinal index to use (IL1 <= IL2 <= NLONS).
1000C     * LEV1      = first altitude level to use (LEV1 >=1).
1001C     * LEV2      = last altitude level to use (LEV1 < LEV2 <= NLEVS).
1002C     * NLONS     = number of longitudes.
1003C     * NLEVS     = number of vertical levels.
1004C     * NAZMTH    = azimuthal array dimension (NAZMTH >= NAZ).
1005C
1006C  Constants in DATA statements.
1007C
1008C     * COS45 = cosine of 45 degrees.           
1009C     * UMIN  = minimum allowable value for zonal or meridional
1010C     *         wind component (m/s).
1011C
1012C  Subroutine arguements.
1013C
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)
1018C
1019C  Internal variables.
1020C
1021      INTEGER  I, L
1022      REAL U, V, COS45, UMIN
1023C
1024      DATA  COS45 / 0.7071068 /
1025      DATA  UMIN / 0.001 /
1026C-----------------------------------------------------------------------     
1027C
1028C  Case with 4 azimuths.
1029C
1030
1031C      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
1046C
1047C  Case with 8 azimuths.
1048C
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
1067C
1068      RETURN
1069C-----------------------------------------------------------------------
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)
1076C
1077C  Calculate zonal and meridional components of the vertical flux
1078C  of horizontal momentum and corresponding wave drag (force per unit mass)
1079C  on a longitude by altitude grid for the Hines' Doppler spread
1080C  GWD parameterization scheme.
1081C  NOTE: only 4 or 8 azimuths can be used.
1082C
1083C  Aug. 6/95 - C. McLandress
1084C
1085C  Output arguements:
1086C
1087C     * FLUX_U = zonal component of vertical momentum flux (Pascals)
1088C     * FLUX_V = meridional component of vertical momentum flux (Pascals)
1089C     * DRAG_U = zonal component of drag (m/s^2).
1090C     * DRAG_V = meridional component of drag (m/s^2).
1091C
1092C  Input arguements:
1093C
1094C     * ALT       = altitudes (m).
1095C     * DENSITY   = background density (kg/m^3).
1096C     * DENSB     = background density at bottom level (kg/m^3).
1097C     * M_ALPHA   = cutoff vertical wavenumber (1/m).
1098C     * AK_ALPHA  = spectral amplitude factor (i.e., {AjKj} in m^4/s^2).
1099C     * K_ALPHA   = horizontal wavenumber (1/m).
1100C     * SLOPE     = slope of incident vertical wavenumber spectrum.
1101C     * NAZ       = actual number of horizontal azimuths used (must be 4 or 8).
1102C     * IL1       = first longitudinal index to use (IL1 >= 1).
1103C     * IL2       = last longitudinal index to use (IL1 <= IL2 <= NLONS).
1104C     * LEV1      = first altitude level to use (LEV1 >=1).
1105C     * LEV2      = last altitude level to use (LEV1 < LEV2 <= NLEVS).
1106C     * NLONS     = number of longitudes.
1107C     * NLEVS     = number of vertical levels.
1108C     * NAZMTH    = azimuthal array dimension (NAZMTH >= NAZ).
1109C
1110C     * LORMS       = .TRUE. for drag computation
1111C
1112C  Constant in DATA statement.
1113C
1114C     * COS45 = cosine of 45 degrees.           
1115C
1116C  Subroutine arguements.
1117C
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)
1126C
1127      LOGICAL LORMS(NLONS)
1128C
1129C  Internal variables.
1130C
1131      INTEGER  I, L, LEV1P, LEV2M
1132      REAL  COS45, PROD2, PROD4, PROD6, PROD8, DENDZ, DENDZ2
1133      DATA  COS45 / 0.7071068 /   
1134C-----------------------------------------------------------------------
1135C
1136      LEV1P = LEV1 + 1
1137      LEV2M = LEV2 - 1
1138      LEV2P = LEV2 + 1
1139C
1140C  Sum over azimuths for case where SLOPE = 1.
1141C
1142      IF (SLOPE.EQ.1.)  THEN
1143C
1144C  Case with 4 azimuths.
1145C
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
1156C
1157C  Case with 8 azimuths.
1158C
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
1177C
1178      END IF
1179C
1180C  Sum over azimuths for case where SLOPE not equal to 1.
1181C
1182      IF (SLOPE.NE.1.)  THEN
1183C
1184C  Case with 4 azimuths.
1185C
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
1198C
1199C  Case with 8 azimuths.
1200C
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
1219C
1220      END IF
1221C
1222C  Calculate flux from sum.
1223C
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
1230C
1231C  Calculate drag at intermediate levels using centered differences
1232C     
1233      DO 120 L = LEV1P,LEV2M
1234        DO 110 I = IL1,IL2
1235        IF (LORMS(I)) THEN
1236ccc       DENDZ2 = DENSITY(I,L) * ( ALT(I,L+1) - ALT(I,L-1) )
1237          DENDZ2 = DENSITY(I,L) * ( ALT(I,L-1) - ALT(I,L) )
1238ccc       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
1240ccc       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
1246C
1247C  Drag at first and last levels using one-side differences.
1248C
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
1270150   CONTINUE
1271      ENDIF
1272C
1273      RETURN
1274C-----------------------------------------------------------------------
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)
1281C
1282C  This routine calculates the gravity wave induced heating and
1283C  diffusion coefficient on a longitude by altitude grid for 
1284C  the Hines' Doppler spread gravity wave drag parameterization scheme.
1285C
1286C  Aug. 6/95 - C. McLandress
1287C
1288C  Output arguements:
1289C
1290C     * HEAT   = gravity wave heating (K/sec).
1291C     * DIFFCO = diffusion coefficient (m^2/sec)
1292C
1293C  Input arguements:
1294C
1295C     * M_ALPHA     = cutoff vertical wavenumber (1/m).
1296C     * DMDZ_ALPHA  = vertical derivative of cutoff wavenumber.
1297C     * AK_ALPHA    = spectral amplitude factor of each azimuth
1298C                     (i.e., {AjKj} in m^4/s^2).
1299C     * K_ALPHA     = horizontal wavenumber of each azimuth (1/m).
1300C     * BVFREQ      = background Brunt Vassala frequency (rad/sec).
1301C     * DENSITY     = background density (kg/m^3).
1302C     * DENSB       = background density at bottom level (kg/m^3).
1303C     * SIGMA_T     = total rms horizontal wind (m/s).
1304C     * VISC_MOL    = molecular viscosity (m^2/s).
1305C     * KSTAR       = typical gravity wave horizontal wavenumber (1/m).
1306C     * SLOPE       = slope of incident vertical wavenumber spectrum.
1307C     * F2,F3,F5,F6 = Hines's fudge factors.
1308C     * NAZ         = actual number of horizontal azimuths used.
1309C     * IL1         = first longitudinal index to use (IL1 >= 1).
1310C     * IL2         = last longitudinal index to use (IL1 <= IL2 <= NLONS).
1311C     * LEV1        = first altitude level to use (LEV1 >=1).
1312C     * LEV2        = last altitude level to use (LEV1 < LEV2 <= NLEVS).
1313C     * NLONS       = number of longitudes.
1314C     * NLEVS       = number of vertical levels.
1315C     * NAZMTH      = azimuthal array dimension (NAZMTH >= NAZ).
1316C
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)
1324C
1325C Internal variables.
1326C
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
1330C
1331C specific heat at constant pressure
1332C
1333      DATA  CPGAS / 1004. /
1334C             
1335C minimum permissible viscosity
1336C
1337      DATA  VISC_MIN / 1.E-10 /       
1338C-----------------------------------------------------------------------     
1339C
1340C  Initialize heating array.
1341C
1342      DO 20 L = 1,NLEVS
1343        DO 10 I = 1,NLONS
1344          HEAT(I,L) = 0.
1345  10    CONTINUE
1346  20  CONTINUE
1347C
1348C  Perform sum over azimuths for case where SLOPE = 1.
1349C
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
1360C
1361C  Perform sum over azimuths for case where SLOPE not 1.
1362C
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
1374C
1375C  Heating and diffusion.
1376C
1377      DO 100 L = LEV1,LEV2
1378        DO 90 I = IL1,IL2
1379C
1380C  Maximum permissible value of cutoff wavenumber is the smaller
1381C  of the instability-induced wavenumber (M_SUB_M_TURB) and
1382C  that imposed by molecular viscosity (M_SUB_M_MOL).
1383C
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 )
1388C
1389          HEATNG = - HEAT(I,L) * F5 * BVFREQ(I,L) / M_SUB_M
1390     ^               * DENSB(I) / DENSITY(I,L)
1391          DIFFCO(I,L) = F6 * HEATNG**0.33333333 / M_SUB_M**1.33333333
1392          HEAT(I,L)   = HEATNG / CPGAS
1393C
1394 90     CONTINUE
1395 100  CONTINUE
1396C
1397      RETURN
1398C-----------------------------------------------------------------------
1399      END
1400
1401      SUBROUTINE HINES_SIGMA (SIGMA_T,SIGMA_ALPHA,SIGSQH_ALPHA,
1402     1                        NAZ,LEV,IL1,IL2,NLONS,NLEVS,NAZMTH)
1403C
1404C  This routine calculates the total rms and azimuthal rms horizontal
1405C  velocities at a given level on a longitude by altitude grid for
1406C  the Hines' Doppler spread GWD parameterization scheme.
1407C  NOTE: only four or eight azimuths can be used.
1408C
1409C  Aug. 7/95 - C. McLandress
1410C
1411C  Output arguements:
1412C
1413C     * SIGMA_T      = total rms horizontal wind (m/s).
1414C     * SIGMA_ALPHA  = total rms wind in each azimuth (m/s).
1415C
1416C  Input arguements:
1417C
1418C     * SIGSQH_ALPHA = portion of wind variance from waves having wave
1419C     *                normals in the alpha azimuth (m/s).
1420C     * NAZ       = actual number of horizontal azimuths used (must be 4 or 8).
1421C     * LEV       = altitude level to process.
1422C     * IL1       = first longitudinal index to use (IL1 >= 1).
1423C     * IL2       = last longitudinal index to use (IL1 <= IL2 <= NLONS).
1424C     * NLONS     = number of longitudes.
1425C     * NLEVS     = number of vertical levels.
1426C     * NAZMTH    = azimuthal array dimension (NAZMTH >= NAZ).
1427C
1428C  Subroutine arguements.
1429C
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)
1435C
1436C  Internal variables.
1437C
1438      INTEGER  I, N
1439      REAL  SUM_EVEN, SUM_ODD
1440C-----------------------------------------------------------------------     
1441C
1442C  Calculate azimuthal rms velocity for the 4 azimuth case.
1443C
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
1454C
1455C  Calculate azimuthal rms velocity for the 8 azimuth case.
1456C
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
1481C
1482C  Calculate total rms velocity.
1483C
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
1495C
1496      RETURN
1497C-----------------------------------------------------------------------     
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)
1503C
1504C  This routine calculates the vertical wavenumber integral
1505C  for a single vertical level at each azimuth on a longitude grid
1506C  for the Hines' Doppler spread GWD parameterization scheme.
1507C  NOTE: (1) only spectral slopes of 1, 1.5 or 2 are permitted.
1508C        (2) the integral is written in terms of the product QM
1509C            which by construction is always less than 1. Series
1510C            solutions are used for small |QM| and analytical solutions
1511C            for remaining values.
1512C
1513C  Aug. 8/95 - C. McLandress
1514C
1515C  Output arguement:
1516C
1517C     * I_ALPHA = Hines' integral.
1518C
1519C  Input arguements:
1520C
1521C     * V_ALPHA = azimuthal wind component (m/s).
1522C     * M_ALPHA = azimuthal cutoff vertical wavenumber (1/m).
1523C     * BVFB    = background Brunt Vassala frequency at model bottom.
1524C     * SLOPE   = slope of initial vertical wavenumber spectrum
1525C     *           (must use SLOPE = 1., 1.5 or 2.)
1526C     * NAZ     = actual number of horizontal azimuths used.
1527C     * LEV     = altitude level to process.
1528C     * IL1     = first longitudinal index to use (IL1 >= 1).
1529C     * IL2     = last longitudinal index to use (IL1 <= IL2 <= NLONS).
1530C     * NLONS   = number of longitudes.
1531C     * NLEVS   = number of vertical levels.
1532C     * NAZMTH  = azimuthal array dimension (NAZMTH >= NAZ).
1533C
1534C     * LORMS       = .TRUE. for drag computation
1535C
1536C  Constants in DATA statements:
1537C
1538C     * QMIN = minimum value of Q_ALPHA (avoids indeterminant form of integral)
1539C     * QM_MIN = minimum value of Q_ALPHA * M_ALPHA (used to avoid numerical
1540C     *          problems).
1541C
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
1547C
1548      LOGICAL LORMS(NLONS)
1549C
1550C  Internal variables.
1551C
1552      INTEGER  I, N
1553      REAL  Q_ALPHA, QM, SQRTQM, Q_MIN, QM_MIN
1554C
1555      DATA  Q_MIN / 1.0 /, QM_MIN / 0.01 /
1556C-----------------------------------------------------------------------     
1557C
1558C  For integer value SLOPE = 1.
1559C
1560      IF (SLOPE .EQ. 1.)  THEN
1561C
1562        DO 20 N = 1,NAZ
1563          DO 10 I = IL1,IL2
1564          IF (LORMS(I)) THEN
1565C
1566            Q_ALPHA = V_ALPHA(I,LEV,N) / BVFB(I)
1567            QM = Q_ALPHA * M_ALPHA(I,LEV,N)
1568C
1569C  If |QM| is small then use first 4 terms series of Taylor series
1570C  expansion of integral in order to avoid indeterminate form of integral,
1571C  otherwise use analytical form of integral.
1572C
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) = ( QM**2/2. + QM**3/3. + QM**4/4.
1578     ^                           + QM**5/5. ) / Q_ALPHA**2
1579              END IF
1580            ELSE
1581              I_ALPHA(I,N) = - ( ALOG(1.-QM) + QM ) / Q_ALPHA**2
1582            END IF
1583C
1584          ENDIF
1585 10       CONTINUE
1586 20     CONTINUE
1587C
1588      END IF
1589C
1590C  For integer value SLOPE = 2.
1591C
1592      IF (SLOPE .EQ. 2.)  THEN
1593C
1594        DO 40 N = 1,NAZ
1595          DO 30 I = IL1,IL2
1596          IF (LORMS(I)) THEN
1597C
1598            Q_ALPHA = V_ALPHA(I,LEV,N) / BVFB(I)
1599            QM = Q_ALPHA * M_ALPHA(I,LEV,N)
1600C
1601C  If |QM| is small then use first 4 terms series of Taylor series
1602C  expansion of integral in order to avoid indeterminate form of integral,
1603C  otherwise use analytical form of integral.
1604C
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) = ( QM**3/3. + QM**4/4. + QM**5/5.
1610     ^                           + QM**6/6. ) / Q_ALPHA**3
1611              END IF
1612            ELSE
1613              I_ALPHA(I,N) = - ( ALOG(1.-QM) + QM + QM**2/2.)
1614     ^                         / Q_ALPHA**3
1615            END IF
1616C
1617          ENDIF
1618 30       CONTINUE
1619 40     CONTINUE
1620C
1621      END IF
1622C
1623C  For real value SLOPE = 1.5
1624C
1625      IF (SLOPE .EQ. 1.5)  THEN
1626C
1627        DO 60 N = 1,NAZ
1628          DO 50 I = IL1,IL2
1629          IF (LORMS(I)) THEN
1630C
1631            Q_ALPHA = V_ALPHA(I,LEV,N) / BVFB(I)
1632            QM = Q_ALPHA * M_ALPHA(I,LEV,N)       
1633C
1634C  If |QM| is small then use first 4 terms series of Taylor series
1635C  expansion of integral in order to avoid indeterminate form of integral,
1636C  otherwise use analytical form of integral.
1637C
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     ^                           + QM**3/4.5 + QM**4/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
1657C
1658          ENDIF
1659 50       CONTINUE
1660 60     CONTINUE
1661C
1662      END IF
1663C
1664C  If integral is negative (which in principal should not happen) then
1665C  print a message and some info since execution will abort when calculating
1666C  the variances.
1667C
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
1687C
1688      RETURN
1689C-----------------------------------------------------------------------
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)
1695C
1696C  This routine specifies various parameters needed for the
1697C  the Hines' Doppler spread gravity wave drag parameterization scheme.
1698C
1699C  Aug. 8/95 - C. McLandress
1700C
1701C  Output arguements:
1702C
1703C     * NAZ        = actual number of horizontal azimuths used
1704C     *              (code set up presently for only NAZ = 4 or 8).
1705C     * SLOPE      = slope of incident vertical wavenumber spectrum
1706C     *              (code set up presently for SLOPE 1., 1.5 or 2.).
1707C     * F1         = "fudge factor" used in calculation of trial value of
1708C     *              azimuthal cutoff wavenumber M_ALPHA (1.2 <= F1 <= 1.9).
1709C     * F2         = "fudge factor" used in calculation of maximum
1710C     *              permissible instabiliy-induced cutoff wavenumber
1711C     *              M_SUB_M_TURB (0.1 <= F2 <= 1.4).
1712C     * F3         = "fudge factor" used in calculation of maximum
1713C     *              permissible molecular viscosity-induced cutoff wavenumber
1714C     *              M_SUB_M_MOL (0.1 <= F2 <= 1.4).
1715C     * F5         = "fudge factor" used in calculation of heating rate
1716C     *              (1 <= F5 <= 3).
1717C     * F6         = "fudge factor" used in calculation of turbulent
1718C     *              diffusivity coefficient.
1719C     * KSTAR      = typical gravity wave horizontal wavenumber (1/m)
1720C     *              used in calculation of M_SUB_M_TURB.
1721C     * ICUTOFF    = 1 to exponentially damp off GWD, heating and diffusion
1722C     *              arrays above ALT_CUTOFF; otherwise arrays not modified.
1723C     * ALT_CUTOFF = altitude in meters above which exponential decay applied.
1724C     * SMCO       = smoother used to smooth cutoff vertical wavenumbers
1725C     *              and total rms winds before calculating drag or heating.
1726C     *              (==> a 1:SMCO:1 stencil used; SMCO >= 1.).
1727C     * NSMAX      = number of times smoother applied ( >= 1),
1728C     *            = 0 means no smoothing performed.
1729C     * IHEATCAL   = 1 to calculate heating rates and diffusion coefficient.
1730C     *            = 0 means only drag and flux calculated.
1731C     * K_ALPHA    = horizontal wavenumber of each azimuth (1/m) which
1732C     *              is set here to KSTAR.
1733C     * IERROR     = error flag.
1734C     *            = 0 no errors.
1735C     *            = 10 ==> NAZ > NAZMTH
1736C     *            = 20 ==> invalid number of azimuths (NAZ must be 4 or 8).
1737C     *            = 30 ==> invalid slope (SLOPE must be 1., 1.5 or 2.).
1738C     *            = 40 ==> invalid smoother (SMCO must be >= 1.)
1739C
1740C  Input arguements:
1741C
1742C     * NMESSG  = output unit number where messages to be printed.
1743C     * NLONS   = number of longitudes.
1744C     * NAZMTH  = azimuthal array dimension (NAZMTH >= NAZ).
1745C
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
1751C
1752C Internal variables.
1753C
1754      INTEGER  I, N
1755C-----------------------------------------------------------------------     
1756C
1757C  Specify constants.
1758C
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
1774c      SMCO       = 1.0
1775      NSMAX      = 5
1776c      NSMAX      = 2
1777      IHEATCAL   = 0
1778C
1779C  Print information to output file.
1780C
1781c      WRITE (NMESSG,6000)
1782c 6000 FORMAT (/' Subroutine HINES_SETUP:')
1783c      WRITE (NMESSG,*)  '  SLOPE = ', SLOPE
1784c      WRITE (NMESSG,*)  '  NAZ = ', NAZ
1785c      WRITE (NMESSG,*)  '  F1,F2,F3  = ', F1, F2, F3
1786c      WRITE (NMESSG,*)  '  F5,F6     = ', F5, F6
1787c      WRITE (NMESSG,*)  '  KSTAR     = ', KSTAR
1788c     >           ,'  COSLAT     = ', COSLAT
1789c      IF (ICUTOFF .EQ. 1)  THEN
1790c        WRITE (NMESSG,*) '  Drag exponentially damped above ',
1791c     &                       ALT_CUTOFF/1.E3
1792c     END IF
1793c      IF (NSMAX.LT.1 )  THEN
1794c        WRITE (NMESSG,*) '  No smoothing of cutoff wavenumbers, etc'
1795c      ELSE
1796c        WRITE (NMESSG,*) '  Cutoff wavenumbers and sig_t smoothed:'
1797c        WRITE (NMESSG,*) '    SMCO  =', SMCO
1798c        WRITE (NMESSG,*) '    NSMAX =', NSMAX
1799c     END IF
1800C
1801C  Check that things are setup correctly and log error if not
1802C
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
1808C
1809C  Use single value for azimuthal-dependent horizontal wavenumber.
1810C
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
1816C
1817      RETURN
1818C-----------------------------------------------------------------------
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)
1826C
1827C  Print out altitude profiles of various quantities from
1828C  Hines' Doppler spread gravity wave drag parameterization scheme.
1829C  (NOTE: only for NAZ = 4 or 8).
1830C
1831C  Aug. 8/95 - C. McLandress
1832C
1833C  Input arguements:
1834C
1835C     * IU_PRINT = 1 to print out values in east-west direction.
1836C     * IV_PRINT = 1 to print out values in north-south direction.
1837C     * NMESSG   = unit number for printed output.
1838C     * ILPRT1   = first longitudinal index to print.
1839C     * ILPRT2   = last longitudinal index to print.
1840C     * LEVPRT1  = first altitude level to print.
1841C     * LEVPRT2  = last altitude level to print.
1842C
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)
1851C
1852C  Internal variables.
1853C
1854      INTEGER  N_EAST, N_WEST, N_NORTH, N_SOUTH
1855      INTEGER  I, L
1856C-----------------------------------------------------------------------
1857C
1858C  Azimuthal indices of cardinal directions.
1859C
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
1870C
1871C  Print out values for range of longitudes.
1872C
1873      DO 100 I = ILPRT1,ILPRT2
1874C
1875C  Print east-west wind, sigmas, cutoff wavenumbers, flux and drag.
1876C
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
1893C
1894C  Print north-south winds, sigmas, cutoff wavenumbers, flux and drag.
1895C
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
1911C
1912 100  CONTINUE
1913C
1914      RETURN
1915C-----------------------------------------------------------------------
1916      END
1917
1918      SUBROUTINE HINES_EXP (DATA,DATA_ZMAX,ALT,ALT_EXP,IORDER,
1919     1                      IL1,IL2,LEV1,LEV2,NLONS,NLEVS)
1920C
1921C  This routine exponentially damps a longitude by altitude array
1922C  of data above a specified altitude.
1923C
1924C  Aug. 13/95 - C. McLandress
1925C
1926C  Output arguements:
1927C
1928C     * DATA = modified data array.
1929C
1930C  Input arguements:
1931C
1932C     * DATA    = original data array.
1933C     * ALT     = altitudes.
1934C     * ALT_EXP = altitude above which exponential decay applied.
1935C     * IORDER  = 1 means vertical levels are indexed from top down
1936C     *           (i.e., highest level indexed 1 and lowest level NLEVS);
1937C     *           .NE. 1 highest level is index NLEVS.
1938C     * IL1     = first longitudinal index to use (IL1 >= 1).
1939C     * IL2     = last longitudinal index to use (IL1 <= IL2 <= NLONS).
1940C     * LEV1    = first altitude level to use (LEV1 >=1).
1941C     * LEV2    = last altitude level to use (LEV1 < LEV2 <= NLEVS).
1942C     * NLONS   = number of longitudes.
1943C     * NLEVS   = number of vertical
1944C
1945C  Input work arrays:
1946C
1947C     * DATA_ZMAX = data values just above altitude ALT_EXP.
1948C
1949      INTEGER  IORDER, IL1, IL2, LEV1, LEV2, NLONS, NLEVS
1950      REAL  ALT_EXP
1951      REAL  DATA(NLONS,NLEVS), DATA_ZMAX(NLONS), ALT(NLONS,NLEVS)
1952C
1953C Internal variables.
1954C
1955      INTEGER  LEVBOT, LEVTOP, LINCR, I, L
1956      REAL  HSCALE
1957      DATA  HSCALE / 5.E3 /
1958C-----------------------------------------------------------------------     
1959C
1960C  Index of lowest altitude level (bottom of drag calculation).
1961C
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
1970C
1971C  Data values at first level above ALT_EXP.
1972C
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
1980C
1981C  Exponentially damp field above ALT_EXP to model top at L=1.
1982C
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
1990C
1991      RETURN
1992C-----------------------------------------------------------------------
1993      END
1994
1995      SUBROUTINE VERT_SMOOTH (DATA,WORK,COEFF,NSMOOTH,
1996     1                        IL1,IL2,LEV1,LEV2,NLONS,NLEVS)
1997C
1998C  Smooth a longitude by altitude array in the vertical over a
1999C  specified number of levels using a three point smoother.
2000C
2001C  NOTE: input array DATA is modified on output!
2002C
2003C  Aug. 3/95 - C. McLandress
2004C
2005C  Output arguement:
2006C
2007C     * DATA    = smoothed array (on output).
2008C
2009C  Input arguements:
2010C
2011C     * DATA    = unsmoothed array of data (on input).
2012C     * WORK    = work array of same dimension as DATA.
2013C     * COEFF   = smoothing coefficient for a 1:COEFF:1 stencil.
2014C     *           (e.g., COEFF = 2 will result in a smoother which
2015C     *           weights the level L gridpoint by two and the two
2016C     *           adjecent levels (L+1 and L-1) by one).
2017C     * NSMOOTH = number of times to smooth in vertical.
2018C     *           (e.g., NSMOOTH=1 means smoothed only once,
2019C     *           NSMOOTH=2 means smoothing repeated twice, etc.)
2020C     * IL1     = first longitudinal index to use (IL1 >= 1).
2021C     * IL2     = last longitudinal index to use (IL1 <= IL2 <= NLONS).
2022C     * LEV1    = first altitude level to use (LEV1 >=1).
2023C     * LEV2    = last altitude level to use (LEV1 < LEV2 <= NLEVS).
2024C     * NLONS   = number of longitudes.
2025C     * NLEVS   = number of vertical levels.
2026C
2027C  Subroutine arguements.
2028C
2029      INTEGER  NSMOOTH, IL1, IL2, LEV1, LEV2, NLONS, NLEVS
2030      REAL  COEFF
2031      REAL  DATA(NLONS,NLEVS), WORK(NLONS,NLEVS)
2032C
2033C  Internal variables.
2034C
2035      INTEGER  I, L, NS, LEV1P, LEV2M
2036      REAL  SUM_WTS
2037C-----------------------------------------------------------------------     
2038C
2039C  Calculate sum of weights.
2040C
2041      SUM_WTS = COEFF + 2.
2042C
2043      LEV1P = LEV1 + 1
2044      LEV2M = LEV2 - 1
2045C
2046C  Smooth NSMOOTH times
2047C
2048      DO 50 NS = 1,NSMOOTH
2049C
2050C  Copy data into work array.
2051C
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
2057C
2058C  Smooth array WORK in vertical direction and put into DATA.
2059C
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
2066C
2067 50   CONTINUE
2068C
2069      RETURN
2070C-----------------------------------------------------------------------
2071      END
2072
2073
2074
2075
2076     
2077     
Note: See TracBrowser for help on using the repository browser.