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

Last change on this file was 2346, checked in by Ehouarn Millour, 9 years ago

Physics/dynamics separation:

  • remove all references to dimensions.h from physics. nbp_lon (==iim) , nbp_lat (==jjm+1) and nbp_lev (==llm) from mod_grid_phy_lmdz should be used instead.
  • added module regular_lonlat_mod in phy_common to store information about the global (lon-lat) grid cell boundaries and centers.

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