Ignore:
Timestamp:
Mar 5, 2014, 2:19:12 PM (11 years ago)
Author:
lguez
Message:

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

-- replaced #include by include.

File:
1 moved

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/hines_gwd.F90

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