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