source: LMDZ6/branches/Amaury_dev/libf/phylmd/hines_gwd.F90 @ 5218

! $Id: hines_gwd.F90 5158 2024-08-02 12:12:03Z abarral $

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
  USE lmdz_YOEGWD, ONLY: GFRCRIT, GKWAKE, GRCRIT, GVCRIT, GKDRAG, GKLIFT, GHMAX, GRAHILO, GSIGCR, NKTOPG, NSTRA, GSSEC, GTSEC, GVSEC, &
          GWD_RANDO_RUWMAX, gwd_rando_sat, GWD_FRONT_RUWMAX, gwd_front_sat
  USE lmdz_yomcst

  IMPLICIT NONE

  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

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

  ! * 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/

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

  iplev = lrefp - 1

  lozpr = .FALSE.



  ! *    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

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



  ! * 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)

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