! $Id: hines_gwd.F90 3102 2017-12-03 20:27:42Z oboucher $

  ! ########################################################################

  ! 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

  ! 1

  ! 1

  ! 1

  ! * VERTICAL POSITIONING ARRAYS.

  ! * LOGICAL SWITCHES TO CONTROL ROOF DRAG, ENVELOP GW DRAG AND

  ! * HINES' DOPPLER SPREADING EXTROWAVE GW DRAG.

  ! * LOZPR IS TRUE FOR ZPR ENHANCEMENT

  ! * WORK ARRAYS.

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

  ! 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

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

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

  rhoh2o = 1000.

  zpcons = (1000.*86400.)/rhoh2o

  ! ym      KLEVM1=KLEV-1

    END DO

  END DO

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

  ! *Constants:

  ! 1

  ! 1    *Arrays

  ! 1

    END DO

  END DO

  ! C

  END DO

  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

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

      END DO

    END DO

  END DO

  ! * INITIALIZE NECESSARY ARRAYS.

    END DO

  END DO

  ! * IF USING HINES SCHEME THEN CALCULATE B V FREQUENCY AT ALL POINTS

  ! * AND SMOOTH BVFREQ.

      dttdsf = (th(i,l)/shxkj(i,l)-th(i,l-1)/shxkj(i,l-1))/ &

      bvfreq(i, l) = sqrt(-dttdsf*sgj(i,l)*(sgj(i,l)**rgocp)/rd)*rg/ptm1(i, l &

    END DO

  END DO

      END IF

      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.

    END DO

  END DO

  ! * ALTITUDE AND DENSITY AT BOTTOM.

  END DO

  ! * ALTITUDE AND DENSITY AT REMAINING LEVELS.

    END DO

  END DO

  ! * INITIALIZE SWITCHES FOR HINES GWD CALCULATION

  ilrms = 0

  END DO

  ! * DEFILE BOTTOM LAUNCH LEVEL

  ! * BACKGROUND WIND MINUS VALUE AT BOTTOM LAUNCH LEVEL.

    END DO

  END DO

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

  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

    END IF

  END DO

  ! * CALCULATE GWD (NOTE THAT DIFFUSION COEFFICIENT AND

  ! * HEATING RATE ONLY CALCULATED IF IHEATCAL = 1).

    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, &

    ! * ADD ON HINES' GWD TENDENCIES TO OROGRAPHIC TENDENCIES AND

    ! * APPLY HINES' GW DRAG ON (UROW,VROW) WORK ARRAYS.

      END DO

    END DO

    ! * END OF HINES CALCULATIONS.

  END IF

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

    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

! /

! /

    v_alpha, sigma_alpha, sigsqh_alpha, ak_alpha, mmin_alpha, i_alpha, &

    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'

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

1   FORMAT (2X, ' error: IORDER NOT ONE! ')

  END IF

  ! Buoyancy and density at bottom level.

  END DO

  ! initialize some variables

      END DO

    END DO

  END DO

    END DO

  END DO

    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, &

  ! 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, &

  ! Smooth cutoff wavenumbers and total rms velocity in the vertical

  ! direction NSMAX times, using FLUX_U as temporary work array.

        END DO

      END DO

      CALL vert_smooth(smoothr1, smoothr2, smco, nsmax, il1, il2, lev1, lev2, &

        END DO

      END DO

    END DO

    CALL vert_smooth(sigma_t, smoothr2, smco, nsmax, il1, il2, lev1, lev2, &

  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, &

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

    CALL hines_exp(drag_u, bvfb, alt, alt_cutoff, iorder, il1, il2, lev1, &

    CALL hines_exp(drag_v, bvfb, alt, alt_cutoff, iorder, il1, il2, lev1, &

  END IF

  ! Print out various arrays for diagnostic purposes.

    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, &

  END IF

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

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

        v_alpha(i, l, n) = (m_alpha(i,l+1,n)-m_alpha(i,l-1,n))/ &

      END DO

    END DO

      v_alpha(i, lev1, n) = (m_alpha(i,lev1p,n)-m_alpha(i,lev1,n))/ &

    END DO

      v_alpha(i, lev2, n) = (m_alpha(i,lev2,n)-m_alpha(i,lev2m,n))/ &

    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, &

  ! Finished.

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

END SUBROUTINE hines_extro0

    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)

  DATA visc_min/1.E-10/

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

  ! PRINT *,'IN HINES_WAVNUM'

  ! Indices of levels to process.

    lend = levtop

    lincr = -1

  ELSE

1   FORMAT (2X, ' error: IORDER NOT ONE! ')

  END IF

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

    END DO

  END DO

  ! Velocity variances at bottom level.

  CALL hines_sigma(sigma_t, sigma_alpha, sigsqh_alpha, naz, levbot, il1, il2, &

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

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

        m_alpha(i, levbot, n) = bvfb(i)/(f1*sigma_alpha(i,levbot,n)+f2* &

        ak_alpha(i, n) = sigsqh_alpha(i, levbot, n)/ &

      END IF

    END DO

  END DO

  ! Calculate quantities from the bottom upwards,

  ! starting one level above bottom.

    ! Level beneath present level.

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

        ELSE

        END IF

      END IF

    END DO

    ! Calculate cutoff wavenumber at this level.

          ! 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, &

          END IF

          ! Reset minimum value of cutoff wavenumber if necessary.

          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, &

    ! Calculate the velocity variances at this level.

    END DO

      END DO

    END DO

    CALL hines_sigma(sigma_t, sigma_alpha, sigsqh_alpha, naz, l, il1, il2, &

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

    ! End of level loop.

  END DO

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

END SUBROUTINE hines_wavnum

    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'