source: LMDZ5/trunk/libf/phylmd/hines_gwd.F90 @ 2306

Last change on this file since 2306 was 2197, checked in by Ehouarn Millour, 10 years ago

Added 'implicit none' statements and proper variable definitions where they were missing.
EM

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