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

! $Id: hines_gwd.F90 2346 2015-08-21 15:13:46Z fairhead $

SUBROUTINE hines_gwd(nlon, nlev, dtime, paphm1x, papm1x, rlat, tx, ux, vx, &
    zustrhi, zvstrhi, d_t_hin, d_u_hin, d_v_hin)

  ! ########################################################################
  ! Parametrization of the momentum flux deposition due to a broad band
  ! spectrum of gravity waves, following Hines (1997a,b), as coded by
  ! McLANDRESS (1995). Modified by McFARLANE and MANZINI (1995-1997)
  ! MAECHAM model stand alone version
  ! ########################################################################


  USE dimphy
  IMPLICIT NONE

  include "YOEGWD.h"
  include "YOMCST.h"

  INTEGER nazmth
  PARAMETER (nazmth=8)

  ! INPUT ARGUMENTS.
  ! ----- ----------

  ! - 2D
  ! PAPHM1   : HALF LEVEL PRESSURE (T-DT)
  ! PAPM1    : FULL LEVEL PRESSURE (T-DT)
  ! PTM1     : TEMPERATURE (T-DT)
  ! PUM1     : ZONAL WIND (T-DT)
  ! PVM1     : MERIDIONAL WIND (T-DT)


  ! REFERENCE.
  ! ----------
  ! SEE MODEL DOCUMENTATION

  ! AUTHOR.
  ! -------

  ! N. MCFARLANE   DKRZ-HAMBURG   MAY 1995
  ! STAND ALONE E. MANZINI MPI-HAMBURG FEBRUARY 1997

  ! BASED ON A COMBINATION OF THE OROGRAPHIC SCHEME BY N.MCFARLANE 1987
  ! AND THE HINES SCHEME AS CODED BY C. MCLANDRESS 1995.



  ! ym      INTEGER KLEVM1

  REAL paphm1(klon, klev+1), papm1(klon, klev)
  REAL ptm1(klon, klev), pum1(klon, klev), pvm1(klon, klev)
  REAL prflux(klon)
  ! 1
  ! 1
  ! 1
  REAL rlat(klon), coslat(klon)

  REAL th(klon, klev), utendgw(klon, klev), vtendgw(klon, klev), &
    pressg(klon), uhs(klon, klev), vhs(klon, klev), zpr(klon)

  ! * VERTICAL POSITIONING ARRAYS.

  REAL sgj(klon, klev), shj(klon, klev), shxkj(klon, klev), dsgj(klon, klev)

  ! * LOGICAL SWITCHES TO CONTROL ROOF DRAG, ENVELOP GW DRAG AND
  ! * HINES' DOPPLER SPREADING EXTROWAVE GW DRAG.
  ! * LOZPR IS TRUE FOR ZPR ENHANCEMENT


  ! * WORK ARRAYS.

  REAL m_alpha(klon, klev, nazmth), v_alpha(klon, klev, nazmth), &
    sigma_alpha(klon, klev, nazmth), sigsqh_alpha(klon, klev, nazmth), &
    drag_u(klon, klev), drag_v(klon, klev), flux_u(klon, klev), &
    flux_v(klon, klev), heat(klon, klev), diffco(klon, klev), &
    bvfreq(klon, klev), density(klon, klev), sigma_t(klon, klev), &
    visc_mol(klon, klev), alt(klon, klev), sigsqmcw(klon, klev, nazmth), &
    sigmatm(klon, klev), ak_alpha(klon, nazmth), k_alpha(klon, nazmth), &
    mmin_alpha(klon, nazmth), i_alpha(klon, nazmth), rmswind(klon), &
    bvfbot(klon), densbot(klon)
  REAL smoothr1(klon, klev), smoothr2(klon, klev)
  REAL sigalpmc(klon, klev, nazmth)
  REAL f2mod(klon, klev)

  ! * THES ARE THE INPUT PARAMETERS FOR HINES ROUTINE AND
  ! * ARE SPECIFIED IN ROUTINE HINES_SETUP. SINCE THIS IS CALLED
  ! * ONLY AT FIRST CALL TO THIS ROUTINE THESE VARIABLES MUST BE SAVED
  ! * FOR USE AT SUBSEQUENT CALLS. THIS CAN BE AVOIDED BY CALLING
  ! * HINES_SETUP IN MAIN PROGRAM AND PASSING THE PARAMETERS AS
  ! * SUBROUTINE ARGUEMENTS.


  REAL rmscon
  INTEGER nmessg, iprint, ilrms
  INTEGER ifl

  INTEGER naz, icutoff, nsmax, iheatcal
  REAL slope, f1, f2, f3, f5, f6, kstar(klon), alt_cutoff, smco

  ! PROVIDED AS INPUT

  INTEGER nlon, nlev

  REAL dtime
  REAL paphm1x(nlon, nlev+1), papm1x(nlon, nlev)
  REAL ux(nlon, nlev), vx(nlon, nlev), tx(nlon, nlev)

  ! VARIABLES FOR OUTPUT


  REAL d_t_hin(nlon, nlev), d_u_hin(nlon, nlev), d_v_hin(nlon, nlev)
  REAL zustrhi(nlon), zvstrhi(nlon)


  ! * LOGICAL SWITCHES TO CONTROL PRECIP ENHANCEMENT AND
  ! * HINES' DOPPLER SPREADING EXTROWAVE GW DRAG.
  ! * LOZPR IS TRUE FOR ZPR ENHANCEMENT

  LOGICAL lozpr, lorms(klon)

  ! LOCAL PARAMETERS TO MAKE THINGS WORK (TEMPORARY VARIABLE)

  REAL rhoh2o, zpcons, rgocp, zlat, dttdsf, ratio, hscal
  INTEGER i, j, l, jl, jk, le, lref, lrefp, levbot

  ! DATA PARAMETERS NEEDED, EXPLAINED LATER

  REAL v0, vmin, dmpscal, taufac, hmin, apibt, cpart, fcrit
  REAL pcrit, pcons
  INTEGER iplev, ierror



  ! PRINT *,' IT IS STARTED HINES GOING ON...'




  ! *    COMPUTATIONAL CONSTANTS.
  ! ------------- ----------


  d_t_hin(:, :) = 0.

  rhoh2o = 1000.
  zpcons = (1000.*86400.)/rhoh2o
  ! ym      KLEVM1=KLEV-1


  DO jl = kidia, kfdia
    paphm1(jl, 1) = paphm1x(jl, klev+1)
    DO jk = 1, klev
      le = klev + 1 - jk
      paphm1(jl, jk+1) = paphm1x(jl, le)
      papm1(jl, jk) = papm1x(jl, le)
      ptm1(jl, jk) = tx(jl, le)
      pum1(jl, jk) = ux(jl, le)
      pvm1(jl, jk) = vx(jl, le)
    END DO
  END DO

  ! Define constants and arrays needed for the ccc/mam gwd scheme
  ! *Constants:

  rgocp = rd/rcpd
  lrefp = klev - 1
  lref = klev - 2
  ! 1
  ! 1    *Arrays
  ! 1
  DO jk = 1, klev
    DO jl = kidia, kfdia
      shj(jl, jk) = papm1(jl, jk)/paphm1(jl, klev+1)
      sgj(jl, jk) = papm1(jl, jk)/paphm1(jl, klev+1)
      dsgj(jl, jk) = (paphm1(jl,jk+1)-paphm1(jl,jk))/paphm1(jl, klev+1)
      shxkj(jl, jk) = (papm1(jl,jk)/paphm1(jl,klev+1))**rgocp
      th(jl, jk) = ptm1(jl, jk)
    END DO
  END DO

  ! C
  DO jl = kidia, kfdia
    pressg(jl) = paphm1(jl, klev+1)
  END DO


  DO jl = kidia, kfdia
    prflux(jl) = 0.0
    zpr(jl) = zpcons*prflux(jl)
    zlat = (rlat(jl)/180.)*rpi
    coslat(jl) = cos(zlat)
  END DO

  ! /#########################################################################
  ! /
  ! /

  ! * AUG. 14/95 - C. MCLANDRESS.
  ! * SEP.    95   N. MCFARLANE.

  ! * THIS ROUTINE CALCULATES THE HORIZONTAL WIND TENDENCIES
  ! * DUE TO MCFARLANE'S OROGRAPHIC GW DRAG SCHEME, HINES'
  ! * DOPPLER SPREAD SCHEME FOR "EXTROWAVES" AND ADDS ON
  ! * ROOF DRAG. IT IS BASED ON THE ROUTINE GWDFLX8.

  ! * LREFP IS THE INDEX OF THE MODEL LEVEL BELOW THE REFERENCE LEVEL
  ! * I/O ARRAYS PASSED FROM MAIN.
  ! * (PRESSG = SURFACE PRESSURE)




  ! * CONSTANTS VALUES DEFINED IN DATA STATEMENT ARE :
  ! * VMIN     = MIMINUM WIND IN THE DIRECTION OF REFERENCE LEVEL
  ! *            WIND BEFORE WE CONSIDER BREAKING TO HAVE OCCURED.
  ! * DMPSCAL  = DAMPING TIME FOR GW DRAG IN SECONDS.
  ! * TAUFAC   = 1/(LENGTH SCALE).
  ! * HMIN     = MIMINUM ENVELOPE HEIGHT REQUIRED TO PRODUCE GW DRAG.
  ! * V0       = VALUE OF WIND THAT APPROXIMATES ZERO.


  DATA vmin/5.0/, v0/1.E-10/, taufac/5.E-6/, hmin/40000./, dmpscal/6.5E+6/, &
    apibt/1.5708/, cpart/0.7/, fcrit/1./

  ! * HINES EXTROWAVE GWD CONSTANTS DEFINED IN DATA STATEMENT ARE:
  ! * RMSCON = ROOT MEAN SQUARE GRAVITY WAVE WIND AT LOWEST LEVEL (M/S).
  ! * NMESSG  = UNIT NUMBER FOR PRINTED MESSAGES.
  ! * IPRINT  = 1 TO DO PRINT OUT SOME HINES ARRAYS.
  ! * IFL     = FIRST CALL FLAG TO HINES_SETUP ("SAVE" IT)
  ! * PCRIT = CRITICAL VALUE OF ZPR (MM/D)
  ! * IPLEV = LEVEL OF APPLICATION OF PRCIT
  ! * PCONS = FACTOR OF ZPR ENHANCEMENT


  DATA pcrit/5./, pcons/4.75/

  iplev = lrefp - 1

  DATA rmscon/1.00/iprint/2/, nmessg/6/
  DATA ifl/0/

  lozpr = .FALSE.

  ! -----------------------------------------------------------------------



  ! * SET ERROR FLAG

  ierror = 0

  ! * SPECIFY VARIOUS PARAMETERS FOR HINES ROUTINE AT VERY FIRST CALL.
  ! * (NOTE THAT ARRAY K_ALPHA IS SPECIFIED SO MAKE SURE THAT
  ! * IT IS NOT OVERWRITTEN LATER ON).

  CALL hines_setup(naz, slope, f1, f2, f3, f5, f6, kstar, icutoff, &
    alt_cutoff, smco, nsmax, iheatcal, k_alpha, ierror, nmessg, klon, nazmth, &
    coslat)
  IF (ierror/=0) GO TO 999

  ! * START GWD CALCULATIONS.

  lref = lrefp - 1


  DO j = 1, nazmth
    DO l = 1, klev
      DO i = kidia, klon
        sigsqmcw(i, l, j) = 0.
      END DO
    END DO
  END DO



  ! * INITIALIZE NECESSARY ARRAYS.

  DO l = 1, klev
    DO i = kidia, kfdia
      utendgw(i, l) = 0.
      vtendgw(i, l) = 0.

      uhs(i, l) = 0.
      vhs(i, l) = 0.

    END DO
  END DO

  ! * IF USING HINES SCHEME THEN CALCULATE B V FREQUENCY AT ALL POINTS
  ! * AND SMOOTH BVFREQ.

  DO l = 2, klev
    DO i = kidia, kfdia
      dttdsf = (th(i,l)/shxkj(i,l)-th(i,l-1)/shxkj(i,l-1))/ &
        (shj(i,l)-shj(i,l-1))
      dttdsf = min(dttdsf, -5./sgj(i,l))
      bvfreq(i, l) = sqrt(-dttdsf*sgj(i,l)*(sgj(i,l)**rgocp)/rd)*rg/ptm1(i, l &
        )
    END DO
  END DO
  DO l = 1, klev
    DO i = kidia, kfdia
      IF (l==1) THEN
        bvfreq(i, l) = bvfreq(i, l+1)
      END IF
      IF (l>1) THEN
        ratio = 5.*log(sgj(i,l)/sgj(i,l-1))
        bvfreq(i, l) = (bvfreq(i,l-1)+ratio*bvfreq(i,l))/(1.+ratio)
      END IF
    END DO
  END DO

  ! * CALCULATE GW DRAG DUE TO HINES' EXTROWAVES
  ! * SET MOLECULAR VISCOSITY TO A VERY SMALL VALUE.
  ! * IF THE MODEL TOP IS GREATER THAN 100 KM THEN THE ACTUAL
  ! * VISCOSITY COEFFICIENT COULD BE SPECIFIED HERE.

  DO l = 1, klev
    DO i = kidia, kfdia
      visc_mol(i, l) = 1.5E-5
      drag_u(i, l) = 0.
      drag_v(i, l) = 0.
      flux_u(i, l) = 0.
      flux_v(i, l) = 0.
      heat(i, l) = 0.
      diffco(i, l) = 0.
    END DO
  END DO

  ! * ALTITUDE AND DENSITY AT BOTTOM.

  DO i = kidia, kfdia
    hscal = rd*ptm1(i, klev)/rg
    density(i, klev) = sgj(i, klev)*pressg(i)/(rg*hscal)
    alt(i, klev) = 0.
  END DO

  ! * ALTITUDE AND DENSITY AT REMAINING LEVELS.

  DO l = klev - 1, 1, -1
    DO i = kidia, kfdia
      hscal = rd*ptm1(i, l)/rg
      alt(i, l) = alt(i, l+1) + hscal*dsgj(i, l)/sgj(i, l)
      density(i, l) = sgj(i, l)*pressg(i)/(rg*hscal)
    END DO
  END DO


  ! * INITIALIZE SWITCHES FOR HINES GWD CALCULATION

  ilrms = 0

  DO i = kidia, kfdia
    lorms(i) = .FALSE.
  END DO


  ! * DEFILE BOTTOM LAUNCH LEVEL

  levbot = iplev

  ! * BACKGROUND WIND MINUS VALUE AT BOTTOM LAUNCH LEVEL.

  DO l = 1, levbot
    DO i = kidia, kfdia
      uhs(i, l) = pum1(i, l) - pum1(i, levbot)
      vhs(i, l) = pvm1(i, l) - pvm1(i, levbot)
    END DO
  END DO

  ! * SPECIFY ROOT MEAN SQUARE WIND AT BOTTOM LAUNCH LEVEL.

  DO i = kidia, kfdia
    rmswind(i) = rmscon
  END DO

  IF (lozpr) THEN
    DO i = kidia, kfdia
      IF (zpr(i)>pcrit) THEN
        rmswind(i) = rmscon + ((zpr(i)-pcrit)/zpr(i))*pcons
      END IF
    END DO
  END IF

  DO i = kidia, kfdia
    IF (rmswind(i)>0.0) THEN
      ilrms = ilrms + 1
      lorms(i) = .TRUE.
    END IF
  END DO

  ! * CALCULATE GWD (NOTE THAT DIFFUSION COEFFICIENT AND
  ! * HEATING RATE ONLY CALCULATED IF IHEATCAL = 1).

  IF (ilrms>0) THEN

    CALL hines_extro0(drag_u, drag_v, heat, diffco, flux_u, flux_v, uhs, vhs, &
      bvfreq, density, visc_mol, alt, rmswind, k_alpha, m_alpha, v_alpha, &
      sigma_alpha, sigsqh_alpha, ak_alpha, mmin_alpha, i_alpha, sigma_t, &
      densbot, bvfbot, 1, iheatcal, icutoff, iprint, nsmax, smco, alt_cutoff, &
      kstar, slope, f1, f2, f3, f5, f6, naz, sigsqmcw, sigmatm, kidia, klon, &
      1, levbot, klon, klev, nazmth, lorms, smoothr1, smoothr2, sigalpmc, &
      f2mod)

    ! * ADD ON HINES' GWD TENDENCIES TO OROGRAPHIC TENDENCIES AND
    ! * APPLY HINES' GW DRAG ON (UROW,VROW) WORK ARRAYS.

    DO l = 1, klev
      DO i = kidia, kfdia
        utendgw(i, l) = utendgw(i, l) + drag_u(i, l)
        vtendgw(i, l) = vtendgw(i, l) + drag_v(i, l)
      END DO
    END DO


    ! * END OF HINES CALCULATIONS.

  END IF

  ! -----------------------------------------------------------------------

  DO jl = kidia, kfdia
    zustrhi(jl) = flux_u(jl, 1)
    zvstrhi(jl) = flux_v(jl, 1)
    DO jk = 1, klev
      le = klev - jk + 1
      d_u_hin(jl, jk) = utendgw(jl, le)*dtime
      d_v_hin(jl, jk) = vtendgw(jl, le)*dtime
    END DO
  END DO

  ! PRINT *,'UTENDGW:',UTENDGW

  ! PRINT *,' HINES HAS BEEN COMPLETED (LONG ISNT IT...)'

  RETURN
999 CONTINUE

  ! * IF ERROR DETECTED THEN ABORT.

  WRITE (nmessg, 6000)
  WRITE (nmessg, 6010) ierror
6000 FORMAT (/' EXECUTION ABORTED IN GWDOREXV')
6010 FORMAT ('     ERROR FLAG =', I4)


  RETURN
END SUBROUTINE hines_gwd
! /
! /


SUBROUTINE hines_extro0(drag_u, drag_v, heat, diffco, flux_u, flux_v, vel_u, &
    vel_v, bvfreq, density, visc_mol, alt, rmswind, k_alpha, m_alpha, &
    v_alpha, sigma_alpha, sigsqh_alpha, ak_alpha, mmin_alpha, i_alpha, &
    sigma_t, densb, bvfb, iorder, iheatcal, icutoff, iprint, nsmax, smco, &
    alt_cutoff, kstar, slope, f1, f2, f3, f5, f6, naz, sigsqmcw, sigmatm, &
    il1, il2, lev1, lev2, nlons, nlevs, nazmth, lorms, smoothr1, smoothr2, &
    sigalpmc, f2mod)

  IMPLICIT NONE

  ! Main routine for Hines' "extrowave" gravity wave parameterization based
  ! on Hines' Doppler spread theory. This routine calculates zonal
  ! and meridional components of gravity wave drag, heating rates
  ! and diffusion coefficient on a longitude by altitude grid.
  ! No "mythical" lower boundary region calculation is made so it
  ! is assumed that lowest level winds are weak (i.e, approximately zero).

  ! Aug. 13/95 - C. McLandress
  ! SEPT. /95  - N.McFarlane

  ! Modifications:

  ! Output arguements:

  ! * DRAG_U = zonal component of gravity wave drag (m/s^2).
  ! * DRAG_V = meridional component of gravity wave drag (m/s^2).
  ! * HEAT   = gravity wave heating (K/sec).
  ! * DIFFCO = diffusion coefficient (m^2/sec)
  ! * FLUX_U = zonal component of vertical momentum flux (Pascals)
  ! * FLUX_V = meridional component of vertical momentum flux (Pascals)

  ! Input arguements:

  ! * VEL_U      = background zonal wind component (m/s).
  ! * VEL_V      = background meridional wind component (m/s).
  ! * BVFREQ     = background Brunt Vassala frequency (radians/sec).
  ! * DENSITY    = background density (kg/m^3)
  ! * VISC_MOL   = molecular viscosity (m^2/s)
  ! * ALT        = altitude of momentum, density, buoyancy levels (m)
  ! *              (NOTE: levels ordered so that ALT(I,1) > ALT(I,2), etc.)
  ! * RMSWIND   = root mean square gravity wave wind at lowest level (m/s).
  ! * K_ALPHA    = horizontal wavenumber of each azimuth (1/m).
  ! * IORDER       = 1 means vertical levels are indexed from top down
  ! *              (i.e., highest level indexed 1 and lowest level NLEVS);
  ! *           .NE. 1 highest level is index NLEVS.
  ! * IHEATCAL   = 1 to calculate heating rates and diffusion coefficient.
  ! * IPRINT     = 1 to print out various arrays.
  ! * ICUTOFF    = 1 to exponentially damp GWD, heating and diffusion
  ! *              arrays above ALT_CUTOFF; otherwise arrays not modified.
  ! * ALT_CUTOFF = altitude in meters above which exponential decay applied.
  ! * SMCO       = smoothing factor used to smooth cutoff vertical
  ! *              wavenumbers and total rms winds in vertical direction
  ! *              before calculating drag or heating
  ! *              (SMCO >= 1 ==> 1:SMCO:1 stencil used).
  ! * NSMAX      = number of times smoother applied ( >= 1),
  ! *            = 0 means no smoothing performed.
  ! * KSTAR      = typical gravity wave horizontal wavenumber (1/m).
  ! * SLOPE      = slope of incident vertical wavenumber spectrum
  ! *              (SLOPE must equal 1., 1.5 or 2.).
  ! * F1 to F6   = Hines's fudge factors (F4 not needed since used for
  ! *              vertical flux of vertical momentum).
  ! * NAZ        = actual number of horizontal azimuths used.
  ! * IL1        = first longitudinal index to use (IL1 >= 1).
  ! * IL2        = last longitudinal index to use (IL1 <= IL2 <= NLONS).
  ! * LEV1       = index of first level for drag calculation.
  ! * LEV2       = index of last level for drag calculation
  ! *              (i.e., LEV1 < LEV2 <= NLEVS).
  ! * NLONS      = number of longitudes.
  ! * NLEVS      = number of vertical levels.
  ! * NAZMTH     = azimuthal array dimension (NAZMTH >= NAZ).

  ! Work arrays.

  ! * M_ALPHA      = cutoff vertical wavenumber (1/m).
  ! * V_ALPHA      = wind component at each azimuth (m/s) and if IHEATCAL=1
  ! *                holds vertical derivative of cutoff wavenumber.
  ! * SIGMA_ALPHA  = total rms wind in each azimuth (m/s).
  ! * SIGSQH_ALPHA = portion of wind variance from waves having wave
  ! *                normals in the alpha azimuth (m/s).
  ! * SIGMA_T      = total rms horizontal wind (m/s).
  ! * AK_ALPHA     = spectral amplitude factor at each azimuth
  ! *                (i.e.,{AjKj}) in m^4/s^2.
  ! * I_ALPHA      = Hines' integral.
  ! * MMIN_ALPHA   = minimum value of cutoff wavenumber.
  ! * DENSB        = background density at bottom level.
  ! * BVFB         = buoyancy frequency at bottom level and
  ! *                work array for ICUTOFF = 1.

  ! * LORMS       = .TRUE. for drag computation

  INTEGER naz, nlons, nlevs, nazmth, il1, il2, lev1, lev2
  INTEGER icutoff, nsmax, iorder, iheatcal, iprint
  REAL kstar(nlons), f1, f2, f3, f5, f6, slope
  REAL alt_cutoff, smco
  REAL drag_u(nlons, nlevs), drag_v(nlons, nlevs)
  REAL heat(nlons, nlevs), diffco(nlons, nlevs)
  REAL flux_u(nlons, nlevs), flux_v(nlons, nlevs)
  REAL vel_u(nlons, nlevs), vel_v(nlons, nlevs)
  REAL bvfreq(nlons, nlevs), density(nlons, nlevs)
  REAL visc_mol(nlons, nlevs), alt(nlons, nlevs)
  REAL rmswind(nlons), bvfb(nlons), densb(nlons)
  REAL sigma_t(nlons, nlevs), sigsqmcw(nlons, nlevs, nazmth)
  REAL sigma_alpha(nlons, nlevs, nazmth), sigmatm(nlons, nlevs)
  REAL sigsqh_alpha(nlons, nlevs, nazmth)
  REAL m_alpha(nlons, nlevs, nazmth), v_alpha(nlons, nlevs, nazmth)
  REAL ak_alpha(nlons, nazmth), k_alpha(nlons, nazmth)
  REAL mmin_alpha(nlons, nazmth), i_alpha(nlons, nazmth)
  REAL smoothr1(nlons, nlevs), smoothr2(nlons, nlevs)
  REAL sigalpmc(nlons, nlevs, nazmth)
  REAL f2mod(nlons, nlevs)

  LOGICAL lorms(nlons)

  ! Internal variables.

  INTEGER levbot, levtop, i, n, l, lev1p, lev2m
  INTEGER ilprt1, ilprt2
  ! -----------------------------------------------------------------------

  ! PRINT *,' IN HINES_EXTRO0'
  lev1p = lev1 + 1
  lev2m = lev2 - 1

  ! Index of lowest altitude level (bottom of drag calculation).

  levbot = lev2
  levtop = lev1
  IF (iorder/=1) THEN
    WRITE (6, 1)
1   FORMAT (2X, ' error: IORDER NOT ONE! ')
  END IF

  ! Buoyancy and density at bottom level.

  DO i = il1, il2
    bvfb(i) = bvfreq(i, levbot)
    densb(i) = density(i, levbot)
  END DO

  ! initialize some variables

  DO n = 1, naz
    DO l = lev1, lev2
      DO i = il1, il2
        m_alpha(i, l, n) = 0.0
      END DO
    END DO
  END DO
  DO l = lev1, lev2
    DO i = il1, il2
      sigma_t(i, l) = 0.0
    END DO
  END DO
  DO n = 1, naz
    DO i = il1, il2
      i_alpha(i, n) = 0.0
    END DO
  END DO

  ! Compute azimuthal wind components from zonal and meridional winds.

  CALL hines_wind(v_alpha, vel_u, vel_v, naz, il1, il2, lev1, lev2, nlons, &
    nlevs, nazmth)

  ! Calculate cutoff vertical wavenumber and velocity variances.

  CALL hines_wavnum(m_alpha, sigma_alpha, sigsqh_alpha, sigma_t, ak_alpha, &
    v_alpha, visc_mol, density, densb, bvfreq, bvfb, rmswind, i_alpha, &
    mmin_alpha, kstar, slope, f1, f2, f3, naz, levbot, levtop, il1, il2, &
    nlons, nlevs, nazmth, sigsqmcw, sigmatm, lorms, sigalpmc, f2mod)
  ! Smooth cutoff wavenumbers and total rms velocity in the vertical
  ! direction NSMAX times, using FLUX_U as temporary work array.

  IF (nsmax>0) THEN
    DO n = 1, naz
      DO l = lev1, lev2
        DO i = il1, il2
          smoothr1(i, l) = m_alpha(i, l, n)
        END DO
      END DO
      CALL vert_smooth(smoothr1, smoothr2, smco, nsmax, il1, il2, lev1, lev2, &
        nlons, nlevs)
      DO l = lev1, lev2
        DO i = il1, il2
          m_alpha(i, l, n) = smoothr1(i, l)
        END DO
      END DO
    END DO
    CALL vert_smooth(sigma_t, smoothr2, smco, nsmax, il1, il2, lev1, lev2, &
      nlons, nlevs)
  END IF

  ! Calculate zonal and meridional components of the
  ! momentum flux and drag.

  CALL hines_flux(flux_u, flux_v, drag_u, drag_v, alt, density, densb, &
    m_alpha, ak_alpha, k_alpha, slope, naz, il1, il2, lev1, lev2, nlons, &
    nlevs, nazmth, lorms)

  ! Cutoff drag above ALT_CUTOFF, using BVFB as temporary work array.

  IF (icutoff==1) THEN
    CALL hines_exp(drag_u, bvfb, alt, alt_cutoff, iorder, il1, il2, lev1, &
      lev2, nlons, nlevs)
    CALL hines_exp(drag_v, bvfb, alt, alt_cutoff, iorder, il1, il2, lev1, &
      lev2, nlons, nlevs)
  END IF

  ! Print out various arrays for diagnostic purposes.

  IF (iprint==1) THEN
    ilprt1 = 15
    ilprt2 = 16
    CALL hines_print(flux_u, flux_v, drag_u, drag_v, alt, sigma_t, &
      sigma_alpha, v_alpha, m_alpha, 1, 1, 6, ilprt1, ilprt2, lev1, lev2, &
      naz, nlons, nlevs, nazmth)
  END IF

  ! If not calculating heating rate and diffusion coefficient then finished.

  IF (iheatcal/=1) RETURN

  ! Calculate vertical derivative of cutoff wavenumber (store
  ! in array V_ALPHA) using centered differences at interior gridpoints
  ! and one-sided differences at first and last levels.

  DO n = 1, naz
    DO l = lev1p, lev2m
      DO i = il1, il2
        v_alpha(i, l, n) = (m_alpha(i,l+1,n)-m_alpha(i,l-1,n))/ &
          (alt(i,l+1)-alt(i,l-1))
      END DO
    END DO
    DO i = il1, il2
      v_alpha(i, lev1, n) = (m_alpha(i,lev1p,n)-m_alpha(i,lev1,n))/ &
        (alt(i,lev1p)-alt(i,lev1))
    END DO
    DO i = il1, il2
      v_alpha(i, lev2, n) = (m_alpha(i,lev2,n)-m_alpha(i,lev2m,n))/ &
        (alt(i,lev2)-alt(i,lev2m))
    END DO
  END DO

  ! Heating rate and diffusion coefficient.

  CALL hines_heat(heat, diffco, m_alpha, v_alpha, ak_alpha, k_alpha, bvfreq, &
    density, densb, sigma_t, visc_mol, kstar, slope, f2, f3, f5, f6, naz, &
    il1, il2, lev1, lev2, nlons, nlevs, nazmth)

  ! Finished.

  RETURN
  ! -----------------------------------------------------------------------
END SUBROUTINE hines_extro0

SUBROUTINE hines_wavnum(m_alpha, sigma_alpha, sigsqh_alpha, sigma_t, &
    ak_alpha, v_alpha, visc_mol, density, densb, bvfreq, bvfb, rms_wind, &
    i_alpha, mmin_alpha, kstar, slope, f1, f2, f3, naz, levbot, levtop, il1, &
    il2, nlons, nlevs, nazmth, sigsqmcw, sigmatm, lorms, sigalpmc, f2mod)
  IMPLICIT NONE
  ! This routine calculates the cutoff vertical wavenumber and velocity
  ! variances on a longitude by altitude grid for the Hines' Doppler
  ! spread gravity wave drag parameterization scheme.
  ! NOTE: (1) only values of four or eight can be used for # azimuths (NAZ).
  ! (2) only values of 1.0, 1.5 or 2.0 can be used for slope (SLOPE).

  ! Aug. 10/95 - C. McLandress

  ! Output arguements:

  ! * M_ALPHA      = cutoff wavenumber at each azimuth (1/m).
  ! * SIGMA_ALPHA  = total rms wind in each azimuth (m/s).
  ! * SIGSQH_ALPHA = portion of wind variance from waves having wave
  ! *                normals in the alpha azimuth (m/s).
  ! * SIGMA_T      = total rms horizontal wind (m/s).
  ! * AK_ALPHA     = spectral amplitude factor at each azimuth
  ! *                (i.e.,{AjKj}) in m^4/s^2.

  ! Input arguements:

  ! * V_ALPHA  = wind component at each azimuth (m/s).
  ! * VISC_MOL = molecular viscosity (m^2/s)
  ! * DENSITY  = background density (kg/m^3).
  ! * DENSB    = background density at model bottom (kg/m^3).
  ! * BVFREQ   = background Brunt Vassala frequency (radians/sec).
  ! * BVFB     = background Brunt Vassala frequency at model bottom.
  ! * RMS_WIND = root mean square gravity wave wind at lowest level (m/s).
  ! * KSTAR    = typical gravity wave horizontal wavenumber (1/m).
  ! * SLOPE    = slope of incident vertical wavenumber spectrum
  ! *            (SLOPE = 1., 1.5 or 2.).
  ! * F1,F2,F3 = Hines's fudge factors.
  ! * NAZ      = actual number of horizontal azimuths used (4 or 8).
  ! * LEVBOT   = index of lowest vertical level.
  ! * LEVTOP   = index of highest vertical level
  ! *            (NOTE: if LEVTOP < LEVBOT then level index
  ! *             increases from top down).
  ! * IL1      = first longitudinal index to use (IL1 >= 1).
  ! * IL2      = last longitudinal index to use (IL1 <= IL2 <= NLONS).
  ! * NLONS    = number of longitudes.
  ! * NLEVS    = number of vertical levels.
  ! * NAZMTH   = azimuthal array dimension (NAZMTH >= NAZ).

  ! * LORMS       = .TRUE. for drag computation

  ! Input work arrays:

  ! * I_ALPHA    = Hines' integral at a single level.
  ! * MMIN_ALPHA = minimum value of cutoff wavenumber.

  INTEGER naz, levbot, levtop, il1, il2, nlons, nlevs, nazmth
  REAL slope, kstar(nlons), f1, f2, f3, f2mfac
  REAL m_alpha(nlons, nlevs, nazmth)
  REAL sigma_alpha(nlons, nlevs, nazmth)
  REAL sigalpmc(nlons, nlevs, nazmth)
  REAL sigsqh_alpha(nlons, nlevs, nazmth)
  REAL sigsqmcw(nlons, nlevs, nazmth)
  REAL sigma_t(nlons, nlevs)
  REAL sigmatm(nlons, nlevs)
  REAL ak_alpha(nlons, nazmth)
  REAL v_alpha(nlons, nlevs, nazmth)
  REAL visc_mol(nlons, nlevs)
  REAL f2mod(nlons, nlevs)
  REAL density(nlons, nlevs), densb(nlons)
  REAL bvfreq(nlons, nlevs), bvfb(nlons), rms_wind(nlons)
  REAL i_alpha(nlons, nazmth), mmin_alpha(nlons, nazmth)

  LOGICAL lorms(nlons)

  ! Internal variables.

  INTEGER i, l, n, lstart, lend, lincr, lbelow
  REAL m_sub_m_turb, m_sub_m_mol, m_trial
  REAL visc, visc_min, azfac, sp1

  ! c      REAL  N_OVER_M(1000), SIGFAC(1000)

  REAL n_over_m(nlons), sigfac(nlons)
  DATA visc_min/1.E-10/
  ! -----------------------------------------------------------------------


  ! PRINT *,'IN HINES_WAVNUM'
  sp1 = slope + 1.

  ! Indices of levels to process.

  IF (levbot>levtop) THEN
    lstart = levbot - 1
    lend = levtop
    lincr = -1
  ELSE
    WRITE (6, 1)
1   FORMAT (2X, ' error: IORDER NOT ONE! ')
  END IF

  ! Use horizontal isotropy to calculate azimuthal variances at bottom level.

  azfac = 1./real(naz)
  DO n = 1, naz
    DO i = il1, il2
      sigsqh_alpha(i, levbot, n) = azfac*rms_wind(i)**2
    END DO
  END DO

  ! Velocity variances at bottom level.

  CALL hines_sigma(sigma_t, sigma_alpha, sigsqh_alpha, naz, levbot, il1, il2, &
    nlons, nlevs, nazmth)

  CALL hines_sigma(sigmatm, sigalpmc, sigsqmcw, naz, levbot, il1, il2, nlons, &
    nlevs, nazmth)

  ! Calculate cutoff wavenumber and spectral amplitude factor
  ! at bottom level where it is assumed that background winds vanish
  ! and also initialize minimum value of cutoff wavnumber.

  DO n = 1, naz
    DO i = il1, il2
      IF (lorms(i)) THEN
        m_alpha(i, levbot, n) = bvfb(i)/(f1*sigma_alpha(i,levbot,n)+f2* &
          sigma_t(i,levbot))
        ak_alpha(i, n) = sigsqh_alpha(i, levbot, n)/ &
          (m_alpha(i,levbot,n)**sp1/sp1)
        mmin_alpha(i, n) = m_alpha(i, levbot, n)
      END IF
    END DO
  END DO

  ! Calculate quantities from the bottom upwards,
  ! starting one level above bottom.

  DO l = lstart, lend, lincr

    ! Level beneath present level.

    lbelow = l - lincr

    ! Calculate N/m_M where m_M is maximum permissible value of the vertical
    ! wavenumber (i.e., m > m_M are obliterated) and N is buoyancy frequency.
    ! m_M is taken as the smaller of the instability-induced
    ! wavenumber (M_SUB_M_TURB) and that imposed by molecular viscosity
    ! (M_SUB_M_MOL). Since variance at this level is not yet known
    ! use value at level below.

    DO i = il1, il2
      IF (lorms(i)) THEN

        f2mfac = sigmatm(i, lbelow)**2
        f2mod(i, lbelow) = 1. + 2.*f2mfac/(f2mfac+sigma_t(i,lbelow)**2)

        visc = amax1(visc_mol(i,l), visc_min)
        m_sub_m_turb = bvfreq(i, l)/(f2*f2mod(i,lbelow)*sigma_t(i,lbelow))
        m_sub_m_mol = (bvfreq(i,l)*kstar(i)/visc)**0.33333333/f3
        IF (m_sub_m_turb<m_sub_m_mol) THEN
          n_over_m(i) = f2*f2mod(i, lbelow)*sigma_t(i, lbelow)
        ELSE
          n_over_m(i) = bvfreq(i, l)/m_sub_m_mol
        END IF
      END IF
    END DO

    ! Calculate cutoff wavenumber at this level.

    DO n = 1, naz
      DO i = il1, il2
        IF (lorms(i)) THEN

          ! Calculate trial value (since variance at this level is not yet
          ! known
          ! use value at level below). If trial value is negative or if it
          ! exceeds
          ! minimum value (not permitted) then set it to minimum value.

          m_trial = bvfb(i)/(f1*(sigma_alpha(i,lbelow,n)+sigalpmc(i,lbelow, &
            n))+n_over_m(i)+v_alpha(i,l,n))
          IF (m_trial<=0. .OR. m_trial>mmin_alpha(i,n)) THEN
            m_trial = mmin_alpha(i, n)
          END IF
          m_alpha(i, l, n) = m_trial

          ! Reset minimum value of cutoff wavenumber if necessary.

          IF (m_alpha(i,l,n)<mmin_alpha(i,n)) THEN
            mmin_alpha(i, n) = m_alpha(i, l, n)
          END IF

        END IF
      END DO
    END DO

    ! Calculate the Hines integral at this level.

    CALL hines_intgrl(i_alpha, v_alpha, m_alpha, bvfb, slope, naz, l, il1, &
      il2, nlons, nlevs, nazmth, lorms)


    ! Calculate the velocity variances at this level.

    DO i = il1, il2
      sigfac(i) = densb(i)/density(i, l)*bvfreq(i, l)/bvfb(i)
    END DO
    DO n = 1, naz
      DO i = il1, il2
        sigsqh_alpha(i, l, n) = sigfac(i)*ak_alpha(i, n)*i_alpha(i, n)
      END DO
    END DO
    CALL hines_sigma(sigma_t, sigma_alpha, sigsqh_alpha, naz, l, il1, il2, &
      nlons, nlevs, nazmth)

    CALL hines_sigma(sigmatm, sigalpmc, sigsqmcw, naz, l, il1, il2, nlons, &
      nlevs, nazmth)

    ! End of level loop.

  END DO

  RETURN
  ! -----------------------------------------------------------------------
END SUBROUTINE hines_wavnum

SUBROUTINE hines_wind(v_alpha, vel_u, vel_v, naz, il1, il2, lev1, lev2, &
    nlons, nlevs, nazmth)
    IMPLICIT NONE
  ! This routine calculates the azimuthal horizontal background wind
  ! components
  ! on a longitude by altitude grid for the case of 4 or 8 azimuths for
  ! the Hines' Doppler spread GWD parameterization scheme.

  ! Aug. 7/95 - C. McLandress

  ! Output arguement:

  ! * V_ALPHA   = background wind component at each azimuth (m/s).
  ! *             (note: first azimuth is in eastward direction
  ! *              and rotate in counterclockwise direction.)

  ! Input arguements:

  ! * VEL_U     = background zonal wind component (m/s).
  ! * VEL_V     = background meridional wind component (m/s).
  ! * NAZ       = actual number of horizontal azimuths used (must be 4 or 8).
  ! * IL1       = first longitudinal index to use (IL1 >= 1).
  ! * IL2       = last longitudinal index to use (IL1 <= IL2 <= NLONS).
  ! * LEV1      = first altitude level to use (LEV1 >=1).
  ! * LEV2      = last altitude level to use (LEV1 < LEV2 <= NLEVS).
  ! * NLONS     = number of longitudes.
  ! * NLEVS     = number of vertical levels.
  ! * NAZMTH    = azimuthal array dimension (NAZMTH >= NAZ).

  ! Constants in DATA statements.

  ! * COS45 = cosine of 45 degrees.
  ! * UMIN  = minimum allowable value for zonal or meridional
  ! *         wind component (m/s).

  ! Subroutine arguements.

  INTEGER naz, il1, il2, lev1, lev2
  INTEGER nlons, nlevs, nazmth
  REAL v_alpha(nlons, nlevs, nazmth)
  REAL vel_u(nlons, nlevs), vel_v(nlons, nlevs)

  ! Internal variables.

  INTEGER i, l
  REAL u, v, cos45, umin

  DATA cos45/0.7071068/
  DATA umin/0.001/
  ! -----------------------------------------------------------------------

  ! Case with 4 azimuths.


  ! PRINT *,'IN HINES_WIND'
  IF (naz==4) THEN
    DO l = lev1, lev2
      DO i = il1, il2
        u = vel_u(i, l)
        v = vel_v(i, l)
        IF (abs(u)<umin) u = umin
        IF (abs(v)<umin) v = umin
        v_alpha(i, l, 1) = u
        v_alpha(i, l, 2) = v
        v_alpha(i, l, 3) = -u
        v_alpha(i, l, 4) = -v
      END DO
    END DO
  END IF

  ! Case with 8 azimuths.

  IF (naz==8) THEN
    DO l = lev1, lev2
      DO i = il1, il2
        u = vel_u(i, l)
        v = vel_v(i, l)
        IF (abs(u)<umin) u = umin
        IF (abs(v)<umin) v = umin
        v_alpha(i, l, 1) = u
        v_alpha(i, l, 2) = cos45*(v+u)
        v_alpha(i, l, 3) = v
        v_alpha(i, l, 4) = cos45*(v-u)
        v_alpha(i, l, 5) = -u
        v_alpha(i, l, 6) = -v_alpha(i, l, 2)
        v_alpha(i, l, 7) = -v
        v_alpha(i, l, 8) = -v_alpha(i, l, 4)
      END DO
    END DO
  END IF

  RETURN
  ! -----------------------------------------------------------------------
END SUBROUTINE hines_wind

SUBROUTINE hines_flux(flux_u, flux_v, drag_u, drag_v, alt, density, densb, &
    m_alpha, ak_alpha, k_alpha, slope, naz, il1, il2, lev1, lev2, nlons, &
    nlevs, nazmth, lorms)
    IMPLICIT NONE
  ! Calculate zonal and meridional components of the vertical flux
  ! of horizontal momentum and corresponding wave drag (force per unit mass)
  ! on a longitude by altitude grid for the Hines' Doppler spread
  ! GWD parameterization scheme.
  ! NOTE: only 4 or 8 azimuths can be used.

  ! Aug. 6/95 - C. McLandress

  ! Output arguements:

  ! * FLUX_U = zonal component of vertical momentum flux (Pascals)
  ! * FLUX_V = meridional component of vertical momentum flux (Pascals)
  ! * DRAG_U = zonal component of drag (m/s^2).
  ! * DRAG_V = meridional component of drag (m/s^2).

  ! Input arguements:

  ! * ALT       = altitudes (m).
  ! * DENSITY   = background density (kg/m^3).
  ! * DENSB     = background density at bottom level (kg/m^3).
  ! * M_ALPHA   = cutoff vertical wavenumber (1/m).
  ! * AK_ALPHA  = spectral amplitude factor (i.e., {AjKj} in m^4/s^2).
  ! * K_ALPHA   = horizontal wavenumber (1/m).
  ! * SLOPE     = slope of incident vertical wavenumber spectrum.
  ! * NAZ       = actual number of horizontal azimuths used (must be 4 or 8).
  ! * IL1       = first longitudinal index to use (IL1 >= 1).
  ! * IL2       = last longitudinal index to use (IL1 <= IL2 <= NLONS).
  ! * LEV1      = first altitude level to use (LEV1 >=1).
  ! * LEV2      = last altitude level to use (LEV1 < LEV2 <= NLEVS).
  ! * NLONS     = number of longitudes.
  ! * NLEVS     = number of vertical levels.
  ! * NAZMTH    = azimuthal array dimension (NAZMTH >= NAZ).

  ! * LORMS       = .TRUE. for drag computation

  ! Constant in DATA statement.

  ! * COS45 = cosine of 45 degrees.

  ! Subroutine arguements.

  INTEGER naz, il1, il2, lev1, lev2
  INTEGER nlons, nlevs, nazmth
  REAL slope
  REAL flux_u(nlons, nlevs), flux_v(nlons, nlevs)
  REAL drag_u(nlons, nlevs), drag_v(nlons, nlevs)
  REAL alt(nlons, nlevs), density(nlons, nlevs), densb(nlons)
  REAL m_alpha(nlons, nlevs, nazmth)
  REAL ak_alpha(nlons, nazmth), k_alpha(nlons, nazmth)

  LOGICAL lorms(nlons)

  ! Internal variables.

  INTEGER i, l, lev1p, lev2m, lev2p
  REAL cos45, prod2, prod4, prod6, prod8, dendz, dendz2
  DATA cos45/0.7071068/
  ! -----------------------------------------------------------------------

  lev1p = lev1 + 1
  lev2m = lev2 - 1
  lev2p = lev2 + 1

  ! Sum over azimuths for case where SLOPE = 1.

  IF (slope==1.) THEN

    ! Case with 4 azimuths.

    IF (naz==4) THEN
      DO l = lev1, lev2
        DO i = il1, il2
          flux_u(i, l) = ak_alpha(i, 1)*k_alpha(i, 1)*m_alpha(i, l, 1) - &
            ak_alpha(i, 3)*k_alpha(i, 3)*m_alpha(i, l, 3)
          flux_v(i, l) = ak_alpha(i, 2)*k_alpha(i, 2)*m_alpha(i, l, 2) - &
            ak_alpha(i, 4)*k_alpha(i, 4)*m_alpha(i, l, 4)
        END DO
      END DO
    END IF

    ! Case with 8 azimuths.

    IF (naz==8) THEN
      DO l = lev1, lev2
        DO i = il1, il2
          prod2 = ak_alpha(i, 2)*k_alpha(i, 2)*m_alpha(i, l, 2)
          prod4 = ak_alpha(i, 4)*k_alpha(i, 4)*m_alpha(i, l, 4)
          prod6 = ak_alpha(i, 6)*k_alpha(i, 6)*m_alpha(i, l, 6)
          prod8 = ak_alpha(i, 8)*k_alpha(i, 8)*m_alpha(i, l, 8)
          flux_u(i, l) = ak_alpha(i, 1)*k_alpha(i, 1)*m_alpha(i, l, 1) - &
            ak_alpha(i, 5)*k_alpha(i, 5)*m_alpha(i, l, 5) + &
            cos45*(prod2-prod4-prod6+prod8)
          flux_v(i, l) = ak_alpha(i, 3)*k_alpha(i, 3)*m_alpha(i, l, 3) - &
            ak_alpha(i, 7)*k_alpha(i, 7)*m_alpha(i, l, 7) + &
            cos45*(prod2+prod4-prod6-prod8)
        END DO
      END DO
    END IF

  END IF

  ! Sum over azimuths for case where SLOPE not equal to 1.

  IF (slope/=1.) THEN

    ! Case with 4 azimuths.

    IF (naz==4) THEN
      DO l = lev1, lev2
        DO i = il1, il2
          flux_u(i, l) = ak_alpha(i, 1)*k_alpha(i, 1)* &
            m_alpha(i, l, 1)**slope - ak_alpha(i, 3)*k_alpha(i, 3)*m_alpha(i, &
            l, 3)**slope
          flux_v(i, l) = ak_alpha(i, 2)*k_alpha(i, 2)* &
            m_alpha(i, l, 2)**slope - ak_alpha(i, 4)*k_alpha(i, 4)*m_alpha(i, &
            l, 4)**slope
        END DO
      END DO
    END IF

    ! Case with 8 azimuths.

    IF (naz==8) THEN
      DO l = lev1, lev2
        DO i = il1, il2
          prod2 = ak_alpha(i, 2)*k_alpha(i, 2)*m_alpha(i, l, 2)**slope
          prod4 = ak_alpha(i, 4)*k_alpha(i, 4)*m_alpha(i, l, 4)**slope
          prod6 = ak_alpha(i, 6)*k_alpha(i, 6)*m_alpha(i, l, 6)**slope
          prod8 = ak_alpha(i, 8)*k_alpha(i, 8)*m_alpha(i, l, 8)**slope
          flux_u(i, l) = ak_alpha(i, 1)*k_alpha(i, 1)* &
            m_alpha(i, l, 1)**slope - ak_alpha(i, 5)*k_alpha(i, 5)*m_alpha(i, &
            l, 5)**slope + cos45*(prod2-prod4-prod6+prod8)
          flux_v(i, l) = ak_alpha(i, 3)*k_alpha(i, 3)* &
            m_alpha(i, l, 3)**slope - ak_alpha(i, 7)*k_alpha(i, 7)*m_alpha(i, &
            l, 7)**slope + cos45*(prod2+prod4-prod6-prod8)
        END DO
      END DO
    END IF

  END IF

  ! Calculate flux from sum.

  DO l = lev1, lev2
    DO i = il1, il2
      flux_u(i, l) = flux_u(i, l)*densb(i)/slope
      flux_v(i, l) = flux_v(i, l)*densb(i)/slope
    END DO
  END DO

  ! Calculate drag at intermediate levels using centered differences

  DO l = lev1p, lev2m
    DO i = il1, il2
      IF (lorms(i)) THEN
        ! cc       DENDZ2 = DENSITY(I,L) * ( ALT(I,L+1) - ALT(I,L-1) )
        dendz2 = density(i, l)*(alt(i,l-1)-alt(i,l))
        ! cc       DRAG_U(I,L) = - ( FLUX_U(I,L+1) - FLUX_U(I,L-1) ) / DENDZ2
        drag_u(i, l) = -(flux_u(i,l-1)-flux_u(i,l))/dendz2
        ! cc       DRAG_V(I,L) = - ( FLUX_V(I,L+1) - FLUX_V(I,L-1) ) / DENDZ2
        drag_v(i, l) = -(flux_v(i,l-1)-flux_v(i,l))/dendz2

      END IF
    END DO
  END DO

  ! Drag at first and last levels using one-side differences.

  DO i = il1, il2
    IF (lorms(i)) THEN
      dendz = density(i, lev1)*(alt(i,lev1)-alt(i,lev1p))
      drag_u(i, lev1) = flux_u(i, lev1)/dendz
      drag_v(i, lev1) = flux_v(i, lev1)/dendz
    END IF
  END DO
  DO i = il1, il2
    IF (lorms(i)) THEN
      dendz = density(i, lev2)*(alt(i,lev2m)-alt(i,lev2))
      drag_u(i, lev2) = -(flux_u(i,lev2m)-flux_u(i,lev2))/dendz
      drag_v(i, lev2) = -(flux_v(i,lev2m)-flux_v(i,lev2))/dendz
    END IF
  END DO
  IF (nlevs>lev2) THEN
    DO i = il1, il2
      IF (lorms(i)) THEN
        dendz = density(i, lev2p)*(alt(i,lev2)-alt(i,lev2p))
        drag_u(i, lev2p) = -flux_u(i, lev2)/dendz
        drag_v(i, lev2p) = -flux_v(i, lev2)/dendz
      END IF
    END DO
  END IF

  RETURN
  ! -----------------------------------------------------------------------
END SUBROUTINE hines_flux

SUBROUTINE hines_heat(heat, diffco, m_alpha, dmdz_alpha, ak_alpha, k_alpha, &
    bvfreq, density, densb, sigma_t, visc_mol, kstar, slope, f2, f3, f5, f6, &
    naz, il1, il2, lev1, lev2, nlons, nlevs, nazmth)
  IMPLICIT NONE
  ! This routine calculates the gravity wave induced heating and
  ! diffusion coefficient on a longitude by altitude grid for
  ! the Hines' Doppler spread gravity wave drag parameterization scheme.

  ! Aug. 6/95 - C. McLandress

  ! Output arguements:

  ! * HEAT   = gravity wave heating (K/sec).
  ! * DIFFCO = diffusion coefficient (m^2/sec)

  ! Input arguements:

  ! * M_ALPHA     = cutoff vertical wavenumber (1/m).
  ! * DMDZ_ALPHA  = vertical derivative of cutoff wavenumber.
  ! * AK_ALPHA    = spectral amplitude factor of each azimuth
  ! (i.e., {AjKj} in m^4/s^2).
  ! * K_ALPHA     = horizontal wavenumber of each azimuth (1/m).
  ! * BVFREQ      = background Brunt Vassala frequency (rad/sec).
  ! * DENSITY     = background density (kg/m^3).
  ! * DENSB       = background density at bottom level (kg/m^3).
  ! * SIGMA_T     = total rms horizontal wind (m/s).
  ! * VISC_MOL    = molecular viscosity (m^2/s).
  ! * KSTAR       = typical gravity wave horizontal wavenumber (1/m).
  ! * SLOPE       = slope of incident vertical wavenumber spectrum.
  ! * F2,F3,F5,F6 = Hines's fudge factors.
  ! * NAZ         = actual number of horizontal azimuths used.
  ! * IL1         = first longitudinal index to use (IL1 >= 1).
  ! * IL2         = last longitudinal index to use (IL1 <= IL2 <= NLONS).
  ! * LEV1        = first altitude level to use (LEV1 >=1).
  ! * LEV2        = last altitude level to use (LEV1 < LEV2 <= NLEVS).
  ! * NLONS       = number of longitudes.
  ! * NLEVS       = number of vertical levels.
  ! * NAZMTH      = azimuthal array dimension (NAZMTH >= NAZ).

  INTEGER naz, il1, il2, lev1, lev2, nlons, nlevs, nazmth
  REAL kstar(nlons), slope, f2, f3, f5, f6
  REAL heat(nlons, nlevs), diffco(nlons, nlevs)
  REAL m_alpha(nlons, nlevs, nazmth), dmdz_alpha(nlons, nlevs, nazmth)
  REAL ak_alpha(nlons, nazmth), k_alpha(nlons, nazmth)
  REAL bvfreq(nlons, nlevs), density(nlons, nlevs), densb(nlons)
  REAL sigma_t(nlons, nlevs), visc_mol(nlons, nlevs)

  ! Internal variables.

  INTEGER i, l, n
  REAL m_sub_m_turb, m_sub_m_mol, m_sub_m, heatng
  REAL visc, visc_min, cpgas, sm1

  ! specific heat at constant pressure

  DATA cpgas/1004./

  ! minimum permissible viscosity

  DATA visc_min/1.E-10/
  ! -----------------------------------------------------------------------

  ! Initialize heating array.

  DO l = 1, nlevs
    DO i = 1, nlons
      heat(i, l) = 0.
    END DO
  END DO

  ! Perform sum over azimuths for case where SLOPE = 1.

  IF (slope==1.) THEN
    DO n = 1, naz
      DO l = lev1, lev2
        DO i = il1, il2
          heat(i, l) = heat(i, l) + ak_alpha(i, n)*k_alpha(i, n)*dmdz_alpha(i &
            , l, n)
        END DO
      END DO
    END DO
  END IF

  ! Perform sum over azimuths for case where SLOPE not 1.

  IF (slope/=1.) THEN
    sm1 = slope - 1.
    DO n = 1, naz
      DO l = lev1, lev2
        DO i = il1, il2
          heat(i, l) = heat(i, l) + ak_alpha(i, n)*k_alpha(i, n)*m_alpha(i, l &
            , n)**sm1*dmdz_alpha(i, l, n)
        END DO
      END DO
    END DO
  END IF

  ! Heating and diffusion.

  DO l = lev1, lev2
    DO i = il1, il2

      ! Maximum permissible value of cutoff wavenumber is the smaller
      ! of the instability-induced wavenumber (M_SUB_M_TURB) and
      ! that imposed by molecular viscosity (M_SUB_M_MOL).

      visc = amax1(visc_mol(i,l), visc_min)
      m_sub_m_turb = bvfreq(i, l)/(f2*sigma_t(i,l))
      m_sub_m_mol = (bvfreq(i,l)*kstar(i)/visc)**0.33333333/f3
      m_sub_m = amin1(m_sub_m_turb, m_sub_m_mol)

      heatng = -heat(i, l)*f5*bvfreq(i, l)/m_sub_m*densb(i)/density(i, l)
      diffco(i, l) = f6*heatng**0.33333333/m_sub_m**1.33333333
      heat(i, l) = heatng/cpgas

    END DO
  END DO

  RETURN
  ! -----------------------------------------------------------------------
END SUBROUTINE hines_heat

SUBROUTINE hines_sigma(sigma_t, sigma_alpha, sigsqh_alpha, naz, lev, il1, &
    il2, nlons, nlevs, nazmth)
  IMPLICIT NONE
  ! This routine calculates the total rms and azimuthal rms horizontal
  ! velocities at a given level on a longitude by altitude grid for
  ! the Hines' Doppler spread GWD parameterization scheme.
  ! NOTE: only four or eight azimuths can be used.

  ! Aug. 7/95 - C. McLandress

  ! Output arguements:

  ! * SIGMA_T      = total rms horizontal wind (m/s).
  ! * SIGMA_ALPHA  = total rms wind in each azimuth (m/s).

  ! Input arguements:

  ! * SIGSQH_ALPHA = portion of wind variance from waves having wave
  ! *                normals in the alpha azimuth (m/s).
  ! * NAZ       = actual number of horizontal azimuths used (must be 4 or 8).
  ! * LEV       = altitude level to process.
  ! * IL1       = first longitudinal index to use (IL1 >= 1).
  ! * IL2       = last longitudinal index to use (IL1 <= IL2 <= NLONS).
  ! * NLONS     = number of longitudes.
  ! * NLEVS     = number of vertical levels.
  ! * NAZMTH    = azimuthal array dimension (NAZMTH >= NAZ).

  ! Subroutine arguements.

  INTEGER lev, naz, il1, il2
  INTEGER nlons, nlevs, nazmth
  REAL sigma_t(nlons, nlevs)
  REAL sigma_alpha(nlons, nlevs, nazmth)
  REAL sigsqh_alpha(nlons, nlevs, nazmth)

  ! Internal variables.

  INTEGER i, n
  REAL sum_even, sum_odd
  ! -----------------------------------------------------------------------

  ! Calculate azimuthal rms velocity for the 4 azimuth case.

  IF (naz==4) THEN
    DO i = il1, il2
      sigma_alpha(i, lev, 1) = sqrt(sigsqh_alpha(i,lev,1)+sigsqh_alpha(i,lev, &
        3))
      sigma_alpha(i, lev, 2) = sqrt(sigsqh_alpha(i,lev,2)+sigsqh_alpha(i,lev, &
        4))
      sigma_alpha(i, lev, 3) = sigma_alpha(i, lev, 1)
      sigma_alpha(i, lev, 4) = sigma_alpha(i, lev, 2)
    END DO
  END IF

  ! Calculate azimuthal rms velocity for the 8 azimuth case.

  IF (naz==8) THEN
    DO i = il1, il2
      sum_odd = (sigsqh_alpha(i,lev,1)+sigsqh_alpha(i,lev,3)+ &
        sigsqh_alpha(i,lev,5)+sigsqh_alpha(i,lev,7))/2.
      sum_even = (sigsqh_alpha(i,lev,2)+sigsqh_alpha(i,lev,4)+ &
        sigsqh_alpha(i,lev,6)+sigsqh_alpha(i,lev,8))/2.
      sigma_alpha(i, lev, 1) = sqrt(sigsqh_alpha(i,lev,1)+sigsqh_alpha(i,lev, &
        5)+sum_even)
      sigma_alpha(i, lev, 2) = sqrt(sigsqh_alpha(i,lev,2)+sigsqh_alpha(i,lev, &
        6)+sum_odd)
      sigma_alpha(i, lev, 3) = sqrt(sigsqh_alpha(i,lev,3)+sigsqh_alpha(i,lev, &
        7)+sum_even)
      sigma_alpha(i, lev, 4) = sqrt(sigsqh_alpha(i,lev,4)+sigsqh_alpha(i,lev, &
        8)+sum_odd)
      sigma_alpha(i, lev, 5) = sigma_alpha(i, lev, 1)
      sigma_alpha(i, lev, 6) = sigma_alpha(i, lev, 2)
      sigma_alpha(i, lev, 7) = sigma_alpha(i, lev, 3)
      sigma_alpha(i, lev, 8) = sigma_alpha(i, lev, 4)
    END DO
  END IF

  ! Calculate total rms velocity.

  DO i = il1, il2
    sigma_t(i, lev) = 0.
  END DO
  DO n = 1, naz
    DO i = il1, il2
      sigma_t(i, lev) = sigma_t(i, lev) + sigsqh_alpha(i, lev, n)
    END DO
  END DO
  DO i = il1, il2
    sigma_t(i, lev) = sqrt(sigma_t(i,lev))
  END DO

  RETURN
  ! -----------------------------------------------------------------------
END SUBROUTINE hines_sigma

SUBROUTINE hines_intgrl(i_alpha, v_alpha, m_alpha, bvfb, slope, naz, lev, &
    il1, il2, nlons, nlevs, nazmth, lorms)
  IMPLICIT NONE
  ! This routine calculates the vertical wavenumber integral
  ! for a single vertical level at each azimuth on a longitude grid
  ! for the Hines' Doppler spread GWD parameterization scheme.
  ! NOTE: (1) only spectral slopes of 1, 1.5 or 2 are permitted.
  ! (2) the integral is written in terms of the product QM
  ! which by construction is always less than 1. Series
  ! solutions are used for small |QM| and analytical solutions
  ! for remaining values.

  ! Aug. 8/95 - C. McLandress

  ! Output arguement:

  ! * I_ALPHA = Hines' integral.

  ! Input arguements:

  ! * V_ALPHA = azimuthal wind component (m/s).
  ! * M_ALPHA = azimuthal cutoff vertical wavenumber (1/m).
  ! * BVFB    = background Brunt Vassala frequency at model bottom.
  ! * SLOPE   = slope of initial vertical wavenumber spectrum
  ! *           (must use SLOPE = 1., 1.5 or 2.)
  ! * NAZ     = actual number of horizontal azimuths used.
  ! * LEV     = altitude level to process.
  ! * IL1     = first longitudinal index to use (IL1 >= 1).
  ! * IL2     = last longitudinal index to use (IL1 <= IL2 <= NLONS).
  ! * NLONS   = number of longitudes.
  ! * NLEVS   = number of vertical levels.
  ! * NAZMTH  = azimuthal array dimension (NAZMTH >= NAZ).

  ! * LORMS       = .TRUE. for drag computation

  ! Constants in DATA statements:

  ! * QMIN = minimum value of Q_ALPHA (avoids indeterminant form of integral)
  ! * QM_MIN = minimum value of Q_ALPHA * M_ALPHA (used to avoid numerical
  ! *          problems).

  INTEGER lev, naz, il1, il2, nlons, nlevs, nazmth
  REAL i_alpha(nlons, nazmth)
  REAL v_alpha(nlons, nlevs, nazmth)
  REAL m_alpha(nlons, nlevs, nazmth)
  REAL bvfb(nlons), slope

  LOGICAL lorms(nlons)

  ! Internal variables.

  INTEGER i, n
  REAL q_alpha, qm, sqrtqm, q_min, qm_min

  DATA q_min/1.0/, qm_min/0.01/
  ! -----------------------------------------------------------------------

  ! For integer value SLOPE = 1.

  IF (slope==1.) THEN

    DO n = 1, naz
      DO i = il1, il2
        IF (lorms(i)) THEN

          q_alpha = v_alpha(i, lev, n)/bvfb(i)
          qm = q_alpha*m_alpha(i, lev, n)

          ! If |QM| is small then use first 4 terms series of Taylor series
          ! expansion of integral in order to avoid indeterminate form of
          ! integral,
          ! otherwise use analytical form of integral.

          IF (abs(q_alpha)<q_min .OR. abs(qm)<qm_min) THEN
            IF (q_alpha==0.) THEN
              i_alpha(i, n) = m_alpha(i, lev, n)**2/2.
            ELSE
              i_alpha(i, n) = (qm**2/2.+qm**3/3.+qm**4/4.+qm**5/5.)/ &
                q_alpha**2
            END IF
          ELSE
            i_alpha(i, n) = -(alog(1.-qm)+qm)/q_alpha**2
          END IF

        END IF
      END DO
    END DO

  END IF

  ! For integer value SLOPE = 2.

  IF (slope==2.) THEN

    DO n = 1, naz
      DO i = il1, il2
        IF (lorms(i)) THEN

          q_alpha = v_alpha(i, lev, n)/bvfb(i)
          qm = q_alpha*m_alpha(i, lev, n)

          ! If |QM| is small then use first 4 terms series of Taylor series
          ! expansion of integral in order to avoid indeterminate form of
          ! integral,
          ! otherwise use analytical form of integral.

          IF (abs(q_alpha)<q_min .OR. abs(qm)<qm_min) THEN
            IF (q_alpha==0.) THEN
              i_alpha(i, n) = m_alpha(i, lev, n)**3/3.
            ELSE
              i_alpha(i, n) = (qm**3/3.+qm**4/4.+qm**5/5.+qm**6/6.)/ &
                q_alpha**3
            END IF
          ELSE
            i_alpha(i, n) = -(alog(1.-qm)+qm+qm**2/2.)/q_alpha**3
          END IF

        END IF
      END DO
    END DO

  END IF

  ! For real value SLOPE = 1.5

  IF (slope==1.5) THEN

    DO n = 1, naz
      DO i = il1, il2
        IF (lorms(i)) THEN

          q_alpha = v_alpha(i, lev, n)/bvfb(i)
          qm = q_alpha*m_alpha(i, lev, n)

          ! If |QM| is small then use first 4 terms series of Taylor series
          ! expansion of integral in order to avoid indeterminate form of
          ! integral,
          ! otherwise use analytical form of integral.

          IF (abs(q_alpha)<q_min .OR. abs(qm)<qm_min) THEN
            IF (q_alpha==0.) THEN
              i_alpha(i, n) = m_alpha(i, lev, n)**2.5/2.5
            ELSE
              i_alpha(i, n) = (qm/2.5+qm**2/3.5+qm**3/4.5+qm**4/5.5)* &
                m_alpha(i, lev, n)**1.5/q_alpha
            END IF
          ELSE
            qm = abs(qm)
            sqrtqm = sqrt(qm)
            IF (q_alpha>=0.) THEN
              i_alpha(i, n) = (alog((1.+sqrtqm)/(1.-sqrtqm))-2.*sqrtqm*(1.+qm &
                /3.))/q_alpha**2.5
            ELSE
              i_alpha(i, n) = 2.*(atan(sqrtqm)+sqrtqm*(qm/3.-1.))/ &
                abs(q_alpha)**2.5
            END IF
          END IF

        END IF
      END DO
    END DO

  END IF

  ! If integral is negative (which in principal should not happen) then
  ! print a message and some info since execution will abort when calculating
  ! the variances.

  ! DO 80 N = 1,NAZ
  ! DO 70 I = IL1,IL2
  ! IF (I_ALPHA(I,N).LT.0.)  THEN
  ! WRITE (6,*)
  ! WRITE (6,*) '******************************'
  ! WRITE (6,*) 'Hines integral I_ALPHA < 0 '
  ! WRITE (6,*) '  longitude I=',I
  ! WRITE (6,*) '  azimuth   N=',N
  ! WRITE (6,*) '  level   LEV=',LEV
  ! WRITE (6,*) '  I_ALPHA =',I_ALPHA(I,N)
  ! WRITE (6,*) '  V_ALPHA =',V_ALPHA(I,LEV,N)
  ! WRITE (6,*) '  M_ALPHA =',M_ALPHA(I,LEV,N)
  ! WRITE (6,*) '  Q_ALPHA =',V_ALPHA(I,LEV,N) / BVFB(I)
  ! WRITE (6,*) '  QM      =',V_ALPHA(I,LEV,N) / BVFB(I)
  ! ^                                * M_ALPHA(I,LEV,N)
  ! WRITE (6,*) '******************************'
  ! END IF
  ! 70     CONTINUE
  ! 80   CONTINUE

  RETURN
  ! -----------------------------------------------------------------------
END SUBROUTINE hines_intgrl

SUBROUTINE hines_setup(naz, slope, f1, f2, f3, f5, f6, kstar, icutoff, &
    alt_cutoff, smco, nsmax, iheatcal, k_alpha, ierror, nmessg, nlons, &
    nazmth, coslat)
  IMPLICIT NONE
  ! This routine specifies various parameters needed for the
  ! the Hines' Doppler spread gravity wave drag parameterization scheme.

  ! Aug. 8/95 - C. McLandress

  ! Output arguements:

  ! * NAZ        = actual number of horizontal azimuths used
  ! *              (code set up presently for only NAZ = 4 or 8).
  ! * SLOPE      = slope of incident vertical wavenumber spectrum
  ! *              (code set up presently for SLOPE 1., 1.5 or 2.).
  ! * F1         = "fudge factor" used in calculation of trial value of
  ! *              azimuthal cutoff wavenumber M_ALPHA (1.2 <= F1 <= 1.9).
  ! * F2         = "fudge factor" used in calculation of maximum
  ! *              permissible instabiliy-induced cutoff wavenumber
  ! *              M_SUB_M_TURB (0.1 <= F2 <= 1.4).
  ! * F3         = "fudge factor" used in calculation of maximum
  ! *              permissible molecular viscosity-induced cutoff wavenumber
  ! *              M_SUB_M_MOL (0.1 <= F2 <= 1.4).
  ! * F5         = "fudge factor" used in calculation of heating rate
  ! *              (1 <= F5 <= 3).
  ! * F6         = "fudge factor" used in calculation of turbulent
  ! *              diffusivity coefficient.
  ! * KSTAR      = typical gravity wave horizontal wavenumber (1/m)
  ! *              used in calculation of M_SUB_M_TURB.
  ! * ICUTOFF    = 1 to exponentially damp off GWD, heating and diffusion
  ! *              arrays above ALT_CUTOFF; otherwise arrays not modified.
  ! * ALT_CUTOFF = altitude in meters above which exponential decay applied.
  ! * SMCO       = smoother used to smooth cutoff vertical wavenumbers
  ! *              and total rms winds before calculating drag or heating.
  ! *              (==> a 1:SMCO:1 stencil used; SMCO >= 1.).
  ! * NSMAX      = number of times smoother applied ( >= 1),
  ! *            = 0 means no smoothing performed.
  ! * IHEATCAL   = 1 to calculate heating rates and diffusion coefficient.
  ! *            = 0 means only drag and flux calculated.
  ! * K_ALPHA    = horizontal wavenumber of each azimuth (1/m) which
  ! *              is set here to KSTAR.
  ! * IERROR     = error flag.
  ! *            = 0 no errors.
  ! *            = 10 ==> NAZ > NAZMTH
  ! *            = 20 ==> invalid number of azimuths (NAZ must be 4 or 8).
  ! *            = 30 ==> invalid slope (SLOPE must be 1., 1.5 or 2.).
  ! *            = 40 ==> invalid smoother (SMCO must be >= 1.)

  ! Input arguements:

  ! * NMESSG  = output unit number where messages to be printed.
  ! * NLONS   = number of longitudes.
  ! * NAZMTH  = azimuthal array dimension (NAZMTH >= NAZ).

  INTEGER naz, nlons, nazmth, iheatcal, icutoff
  INTEGER nmessg, nsmax, ierror
  REAL kstar(nlons), slope, f1, f2, f3, f5, f6, alt_cutoff, smco
  REAL k_alpha(nlons, nazmth), coslat(nlons)
  REAL ksmin, ksmax

  ! Internal variables.

  INTEGER i, n
  ! -----------------------------------------------------------------------

  ! Specify constants.

  naz = 8
  slope = 1.
  f1 = 1.5
  f2 = 0.3
  f3 = 1.0
  f5 = 3.0
  f6 = 1.0
  ksmin = 1.E-5
  ksmax = 1.E-4
  DO i = 1, nlons
    kstar(i) = ksmin/(coslat(i)+(ksmin/ksmax))
  END DO
  icutoff = 1
  alt_cutoff = 105.E3
  smco = 2.0
  ! SMCO       = 1.0
  nsmax = 5
  ! NSMAX      = 2
  iheatcal = 0

  ! Print information to output file.

  ! WRITE (NMESSG,6000)
  ! 6000 FORMAT (/' Subroutine HINES_SETUP:')
  ! WRITE (NMESSG,*)  '  SLOPE = ', SLOPE
  ! WRITE (NMESSG,*)  '  NAZ = ', NAZ
  ! WRITE (NMESSG,*)  '  F1,F2,F3  = ', F1, F2, F3
  ! WRITE (NMESSG,*)  '  F5,F6     = ', F5, F6
  ! WRITE (NMESSG,*)  '  KSTAR     = ', KSTAR
  ! >           ,'  COSLAT     = ', COSLAT
  ! IF (ICUTOFF .EQ. 1)  THEN
  ! WRITE (NMESSG,*) '  Drag exponentially damped above ',
  ! &                       ALT_CUTOFF/1.E3
  ! END IF
  ! IF (NSMAX.LT.1 )  THEN
  ! WRITE (NMESSG,*) '  No smoothing of cutoff wavenumbers, etc'
  ! ELSE
  ! WRITE (NMESSG,*) '  Cutoff wavenumbers and sig_t smoothed:'
  ! WRITE (NMESSG,*) '    SMCO  =', SMCO
  ! WRITE (NMESSG,*) '    NSMAX =', NSMAX
  ! END IF

  ! Check that things are setup correctly and log error if not

  ierror = 0
  IF (naz>nazmth) ierror = 10
  IF (naz/=4 .AND. naz/=8) ierror = 20
  IF (slope/=1. .AND. slope/=1.5 .AND. slope/=2.) ierror = 30
  IF (smco<1.) ierror = 40

  ! Use single value for azimuthal-dependent horizontal wavenumber.

  DO n = 1, naz
    DO i = 1, nlons
      k_alpha(i, n) = kstar(i)
    END DO
  END DO

  RETURN
  ! -----------------------------------------------------------------------
END SUBROUTINE hines_setup

SUBROUTINE hines_print(flux_u, flux_v, drag_u, drag_v, alt, sigma_t, &
    sigma_alpha, v_alpha, m_alpha, iu_print, iv_print, nmessg, ilprt1, &
    ilprt2, levprt1, levprt2, naz, nlons, nlevs, nazmth)
  IMPLICIT NONE
  ! Print out altitude profiles of various quantities from
  ! Hines' Doppler spread gravity wave drag parameterization scheme.
  ! (NOTE: only for NAZ = 4 or 8).

  ! Aug. 8/95 - C. McLandress

  ! Input arguements:

  ! * IU_PRINT = 1 to print out values in east-west direction.
  ! * IV_PRINT = 1 to print out values in north-south direction.
  ! * NMESSG   = unit number for printed output.
  ! * ILPRT1   = first longitudinal index to print.
  ! * ILPRT2   = last longitudinal index to print.
  ! * LEVPRT1  = first altitude level to print.
  ! * LEVPRT2  = last altitude level to print.

  INTEGER naz, ilprt1, ilprt2, levprt1, levprt2
  INTEGER nlons, nlevs, nazmth
  INTEGER iu_print, iv_print, nmessg
  REAL flux_u(nlons, nlevs), flux_v(nlons, nlevs)
  REAL drag_u(nlons, nlevs), drag_v(nlons, nlevs)
  REAL alt(nlons, nlevs), sigma_t(nlons, nlevs)
  REAL sigma_alpha(nlons, nlevs, nazmth)
  REAL v_alpha(nlons, nlevs, nazmth), m_alpha(nlons, nlevs, nazmth)

  ! Internal variables.

  INTEGER n_east, n_west, n_north, n_south
  INTEGER i, l
  ! -----------------------------------------------------------------------

  ! Azimuthal indices of cardinal directions.

  n_east = 1
  IF (naz==4) THEN
    n_west = 3
    n_north = 2
    n_south = 4
  ELSE IF (naz==8) THEN
    n_west = 5
    n_north = 3
    n_south = 7
  END IF

  ! Print out values for range of longitudes.

  DO i = ilprt1, ilprt2

    ! Print east-west wind, sigmas, cutoff wavenumbers, flux and drag.

    IF (iu_print==1) THEN
      WRITE (nmessg, *)
      WRITE (nmessg, 6001) i
      WRITE (nmessg, 6005)
6001  FORMAT ('Hines GW (east-west) at longitude I =', I3)
6005  FORMAT (15X, ' U ', 2X, 'sig_E', 2X, 'sig_T', 3X, 'm_E', 4X, 'm_W', 4X, &
        'fluxU', 5X, 'gwdU')
      DO l = levprt1, levprt2
        WRITE (nmessg, 6701) alt(i, l)/1.E3, v_alpha(i, l, n_east), &
          sigma_alpha(i, l, n_east), sigma_t(i, l), &
          m_alpha(i, l, n_east)*1.E3, m_alpha(i, l, n_west)*1.E3, &
          flux_u(i, l)*1.E5, drag_u(i, l)*24.*3600.
      END DO
6701  FORMAT (' z=', F7.2, 1X, 3F7.1, 2F7.3, F9.4, F9.3)
    END IF

    ! Print north-south winds, sigmas, cutoff wavenumbers, flux and drag.

    IF (iv_print==1) THEN
      WRITE (nmessg, *)
      WRITE (nmessg, 6002) i
6002  FORMAT ('Hines GW (north-south) at longitude I =', I3)
      WRITE (nmessg, 6006)
6006  FORMAT (15X, ' V ', 2X, 'sig_N', 2X, 'sig_T', 3X, 'm_N', 4X, 'm_S', 4X, &
        'fluxV', 5X, 'gwdV')
      DO l = levprt1, levprt2
        WRITE (nmessg, 6701) alt(i, l)/1.E3, v_alpha(i, l, n_north), &
          sigma_alpha(i, l, n_north), sigma_t(i, l), &
          m_alpha(i, l, n_north)*1.E3, m_alpha(i, l, n_south)*1.E3, &
          flux_v(i, l)*1.E5, drag_v(i, l)*24.*3600.
      END DO
    END IF

  END DO

  RETURN
  ! -----------------------------------------------------------------------
END SUBROUTINE hines_print

SUBROUTINE hines_exp(data, data_zmax, alt, alt_exp, iorder, il1, il2, lev1, &
    lev2, nlons, nlevs)
  IMPLICIT NONE
  ! This routine exponentially damps a longitude by altitude array
  ! of data above a specified altitude.

  ! Aug. 13/95 - C. McLandress

  ! Output arguements:

  ! * DATA = modified data array.

  ! Input arguements:

  ! * DATA    = original data array.
  ! * ALT     = altitudes.
  ! * ALT_EXP = altitude above which exponential decay applied.
  ! * IORDER    = 1 means vertical levels are indexed from top down
  ! *           (i.e., highest level indexed 1 and lowest level NLEVS);
  ! *           .NE. 1 highest level is index NLEVS.
  ! * IL1     = first longitudinal index to use (IL1 >= 1).
  ! * IL2     = last longitudinal index to use (IL1 <= IL2 <= NLONS).
  ! * LEV1    = first altitude level to use (LEV1 >=1).
  ! * LEV2    = last altitude level to use (LEV1 < LEV2 <= NLEVS).
  ! * NLONS   = number of longitudes.
  ! * NLEVS   = number of vertical

  ! Input work arrays:

  ! * DATA_ZMAX = data values just above altitude ALT_EXP.

  INTEGER iorder, il1, il2, lev1, lev2, nlons, nlevs
  REAL alt_exp
  REAL data(nlons, nlevs), data_zmax(nlons), alt(nlons, nlevs)

  ! Internal variables.

  INTEGER levbot, levtop, lincr, i, l
  REAL hscale
  DATA hscale/5.E3/
  ! -----------------------------------------------------------------------

  ! Index of lowest altitude level (bottom of drag calculation).

  levbot = lev2
  levtop = lev1
  lincr = 1
  IF (iorder/=1) THEN
    levbot = lev1
    levtop = lev2
    lincr = -1
  END IF

  ! Data values at first level above ALT_EXP.

  DO i = il1, il2
    DO l = levtop, levbot, lincr
      IF (alt(i,l)>=alt_exp) THEN
        data_zmax(i) = data(i, l)
      END IF
    END DO
  END DO

  ! Exponentially damp field above ALT_EXP to model top at L=1.

  DO l = 1, lev2
    DO i = il1, il2
      IF (alt(i,l)>=alt_exp) THEN
        data(i, l) = data_zmax(i)*exp((alt_exp-alt(i,l))/hscale)
      END IF
    END DO
  END DO

  RETURN
  ! -----------------------------------------------------------------------
END SUBROUTINE hines_exp

SUBROUTINE vert_smooth(data, work, coeff, nsmooth, il1, il2, lev1, lev2, &
    nlons, nlevs)
  IMPLICIT NONE
  ! Smooth a longitude by altitude array in the vertical over a
  ! specified number of levels using a three point smoother.

  ! NOTE: input array DATA is modified on output!

  ! Aug. 3/95 - C. McLandress

  ! Output arguement:

  ! * DATA    = smoothed array (on output).

  ! Input arguements:

  ! * DATA    = unsmoothed array of data (on input).
  ! * WORK    = work array of same dimension as DATA.
  ! * COEFF   = smoothing coefficient for a 1:COEFF:1 stencil.
  ! *           (e.g., COEFF = 2 will result in a smoother which
  ! *           weights the level L gridpoint by two and the two
  ! *           adjecent levels (L+1 and L-1) by one).
  ! * NSMOOTH = number of times to smooth in vertical.
  ! *           (e.g., NSMOOTH=1 means smoothed only once,
  ! *           NSMOOTH=2 means smoothing repeated twice, etc.)
  ! * IL1     = first longitudinal index to use (IL1 >= 1).
  ! * IL2     = last longitudinal index to use (IL1 <= IL2 <= NLONS).
  ! * LEV1    = first altitude level to use (LEV1 >=1).
  ! * LEV2    = last altitude level to use (LEV1 < LEV2 <= NLEVS).
  ! * NLONS   = number of longitudes.
  ! * NLEVS   = number of vertical levels.

  ! Subroutine arguements.

  INTEGER nsmooth, il1, il2, lev1, lev2, nlons, nlevs
  REAL coeff
  REAL data(nlons, nlevs), work(nlons, nlevs)

  ! Internal variables.

  INTEGER i, l, ns, lev1p, lev2m
  REAL sum_wts
  ! -----------------------------------------------------------------------

  ! Calculate sum of weights.

  sum_wts = coeff + 2.

  lev1p = lev1 + 1
  lev2m = lev2 - 1

  ! Smooth NSMOOTH times

  DO ns = 1, nsmooth

    ! Copy data into work array.

    DO l = lev1, lev2
      DO i = il1, il2
        work(i, l) = data(i, l)
      END DO
    END DO

    ! Smooth array WORK in vertical direction and put into DATA.

    DO l = lev1p, lev2m
      DO i = il1, il2
        data(i, l) = (work(i,l+1)+coeff*work(i,l)+work(i,l-1))/sum_wts
      END DO
    END DO

  END DO

  RETURN
  ! -----------------------------------------------------------------------
END SUBROUTINE vert_smooth