! $Id: hines_gwd.F90 5158 2024-08-02 12:12:03Z abarral $ SUBROUTINE hines_gwd(nlon, nlev, dtime, paphm1x, papm1x, rlat, tx, ux, vx, & zustrhi, zvstrhi, d_t_hin, d_u_hin, d_v_hin) ! ######################################################################## ! Parametrization of the momentum flux deposition due to a broad band ! spectrum of gravity waves, following Hines (1997a,b), as coded by ! McLANDRESS (1995). Modified by McFARLANE and MANZINI (1995-1997) ! MAECHAM model stand alone version ! ######################################################################## USE dimphy USE lmdz_YOEGWD, ONLY: GFRCRIT, GKWAKE, GRCRIT, GVCRIT, GKDRAG, GKLIFT, GHMAX, GRAHILO, GSIGCR, NKTOPG, NSTRA, GSSEC, GTSEC, GVSEC, & GWD_RANDO_RUWMAX, gwd_rando_sat, GWD_FRONT_RUWMAX, gwd_front_sat USE lmdz_yomcst IMPLICIT NONE INTEGER nazmth PARAMETER (nazmth = 8) ! INPUT ARGUMENTS. ! ----- ---------- ! - 2D ! PAPHM1 : HALF LEVEL PRESSURE (T-DT) ! PAPM1 : FULL LEVEL PRESSURE (T-DT) ! PTM1 : TEMPERATURE (T-DT) ! PUM1 : ZONAL WIND (T-DT) ! PVM1 : MERIDIONAL WIND (T-DT) ! REFERENCE. ! ---------- ! SEE MODEL DOCUMENTATION ! AUTHOR. ! ------- ! N. MCFARLANE DKRZ-HAMBURG MAY 1995 ! STAND ALONE E. MANZINI MPI-HAMBURG FEBRUARY 1997 ! BASED ON A COMBINATION OF THE OROGRAPHIC SCHEME BY N.MCFARLANE 1987 ! AND THE HINES SCHEME AS CODED BY C. MCLANDRESS 1995. ! ym INTEGER KLEVM1 REAL paphm1(klon, klev + 1), papm1(klon, klev) REAL ptm1(klon, klev), pum1(klon, klev), pvm1(klon, klev) REAL prflux(klon) ! 1 ! 1 ! 1 REAL rlat(klon), coslat(klon) REAL th(klon, klev), utendgw(klon, klev), vtendgw(klon, klev), & pressg(klon), uhs(klon, klev), vhs(klon, klev), zpr(klon) ! * VERTICAL POSITIONING ARRAYS. REAL sgj(klon, klev), shj(klon, klev), shxkj(klon, klev), dsgj(klon, klev) ! * LOGICAL SWITCHES TO CONTROL ROOF DRAG, ENVELOP GW DRAG AND ! * HINES' DOPPLER SPREADING EXTROWAVE GW DRAG. ! * LOZPR IS TRUE FOR ZPR ENHANCEMENT ! * WORK ARRAYS. REAL m_alpha(klon, klev, nazmth), v_alpha(klon, klev, nazmth), & sigma_alpha(klon, klev, nazmth), sigsqh_alpha(klon, klev, nazmth), & drag_u(klon, klev), drag_v(klon, klev), flux_u(klon, klev), & flux_v(klon, klev), heat(klon, klev), diffco(klon, klev), & bvfreq(klon, klev), density(klon, klev), sigma_t(klon, klev), & visc_mol(klon, klev), alt(klon, klev), sigsqmcw(klon, klev, nazmth), & sigmatm(klon, klev), ak_alpha(klon, nazmth), k_alpha(klon, nazmth), & mmin_alpha(klon, nazmth), i_alpha(klon, nazmth), rmswind(klon), & bvfbot(klon), densbot(klon) REAL smoothr1(klon, klev), smoothr2(klon, klev) REAL sigalpmc(klon, klev, nazmth) REAL f2mod(klon, klev) ! * THES ARE THE INPUT PARAMETERS FOR HINES ROUTINE AND ! * ARE SPECIFIED IN ROUTINE HINES_SETUP. SINCE THIS IS CALLED ! * ONLY AT FIRST CALL TO THIS ROUTINE THESE VARIABLES MUST BE SAVED ! * FOR USE AT SUBSEQUENT CALLS. THIS CAN BE AVOIDED BY CALLING ! * HINES_SETUP IN MAIN PROGRAM AND PASSING THE PARAMETERS AS ! * SUBROUTINE ARGUEMENTS. REAL rmscon INTEGER nmessg, iprint, ilrms INTEGER ifl INTEGER naz, icutoff, nsmax, iheatcal REAL slope, f1, f2, f3, f5, f6, kstar(klon), alt_cutoff, smco ! PROVIDED AS INPUT INTEGER nlon, nlev REAL dtime REAL paphm1x(nlon, nlev + 1), papm1x(nlon, nlev) REAL ux(nlon, nlev), vx(nlon, nlev), tx(nlon, nlev) ! VARIABLES FOR OUTPUT REAL d_t_hin(nlon, nlev), d_u_hin(nlon, nlev), d_v_hin(nlon, nlev) REAL zustrhi(nlon), zvstrhi(nlon) ! * LOGICAL SWITCHES TO CONTROL PRECIP ENHANCEMENT AND ! * HINES' DOPPLER SPREADING EXTROWAVE GW DRAG. ! * LOZPR IS TRUE FOR ZPR ENHANCEMENT LOGICAL lozpr, lorms(klon) ! LOCAL PARAMETERS TO MAKE THINGS WORK (TEMPORARY VARIABLE) REAL rhoh2o, zpcons, rgocp, zlat, dttdsf, ratio, hscal INTEGER i, j, l, jl, jk, le, lref, lrefp, levbot ! DATA PARAMETERS NEEDED, EXPLAINED LATER REAL v0, vmin, dmpscal, taufac, hmin, apibt, cpart, fcrit REAL pcrit, pcons INTEGER iplev, ierror ! /######################################################################### ! / ! / ! * AUG. 14/95 - C. MCLANDRESS. ! * SEP. 95 N. MCFARLANE. ! * THIS ROUTINE CALCULATES THE HORIZONTAL WIND TENDENCIES ! * DUE TO MCFARLANE'S OROGRAPHIC GW DRAG SCHEME, HINES' ! * DOPPLER SPREAD SCHEME FOR "EXTROWAVES" AND ADDS ON ! * ROOF DRAG. IT IS BASED ON THE ROUTINE GWDFLX8. ! * LREFP IS THE INDEX OF THE MODEL LEVEL BELOW THE REFERENCE LEVEL ! * I/O ARRAYS PASSED FROM MAIN. ! * (PRESSG = SURFACE PRESSURE) ! * CONSTANTS VALUES DEFINED IN DATA STATEMENT ARE : ! * VMIN = MIMINUM WIND IN THE DIRECTION OF REFERENCE LEVEL ! * WIND BEFORE WE CONSIDER BREAKING TO HAVE OCCURED. ! * DMPSCAL = DAMPING TIME FOR GW DRAG IN SECONDS. ! * TAUFAC = 1/(LENGTH SCALE). ! * HMIN = MIMINUM ENVELOPE HEIGHT REQUIRED TO PRODUCE GW DRAG. ! * V0 = VALUE OF WIND THAT APPROXIMATES ZERO. DATA vmin/5.0/, v0/1.E-10/, taufac/5.E-6/, hmin/40000./, dmpscal/6.5E+6/, & apibt/1.5708/, cpart/0.7/, fcrit/1./ ! * HINES EXTROWAVE GWD CONSTANTS DEFINED IN DATA STATEMENT ARE: ! * RMSCON = ROOT MEAN SQUARE GRAVITY WAVE WIND AT LOWEST LEVEL (M/S). ! * NMESSG = UNIT NUMBER FOR PRINTED MESSAGES. ! * IPRINT = 1 TO DO PRINT OUT SOME HINES ARRAYS. ! * IFL = FIRST CALL FLAG TO HINES_SETUP ("SAVE" IT) ! * PCRIT = CRITICAL VALUE OF ZPR (MM/D) ! * IPLEV = LEVEL OF APPLICATION OF PRCIT ! * PCONS = FACTOR OF ZPR ENHANCEMENT DATA pcrit/5./, pcons/4.75/ DATA rmscon/1.00/iprint/2/, nmessg/6/ DATA ifl/0/ iplev = lrefp - 1 lozpr = .FALSE. ! * COMPUTATIONAL CONSTANTS. ! ------------- ---------- d_t_hin(:, :) = 0. rhoh2o = 1000. zpcons = (1000. * 86400.) / rhoh2o ! ym KLEVM1=KLEV-1 DO jl = kidia, kfdia paphm1(jl, 1) = paphm1x(jl, klev + 1) DO jk = 1, klev le = klev + 1 - jk paphm1(jl, jk + 1) = paphm1x(jl, le) papm1(jl, jk) = papm1x(jl, le) ptm1(jl, jk) = tx(jl, le) pum1(jl, jk) = ux(jl, le) pvm1(jl, jk) = vx(jl, le) END DO END DO ! Define constants and arrays needed for the ccc/mam gwd scheme ! *Constants: rgocp = rd / rcpd lrefp = klev - 1 lref = klev - 2 ! 1 ! 1 *Arrays ! 1 DO jk = 1, klev DO jl = kidia, kfdia shj(jl, jk) = papm1(jl, jk) / paphm1(jl, klev + 1) sgj(jl, jk) = papm1(jl, jk) / paphm1(jl, klev + 1) dsgj(jl, jk) = (paphm1(jl, jk + 1) - paphm1(jl, jk)) / paphm1(jl, klev + 1) shxkj(jl, jk) = (papm1(jl, jk) / paphm1(jl, klev + 1))**rgocp th(jl, jk) = ptm1(jl, jk) END DO END DO ! C DO jl = kidia, kfdia pressg(jl) = paphm1(jl, klev + 1) END DO DO jl = kidia, kfdia prflux(jl) = 0.0 zpr(jl) = zpcons * prflux(jl) zlat = (rlat(jl) / 180.) * rpi coslat(jl) = cos(zlat) END DO ! ----------------------------------------------------------------------- ! * SET ERROR FLAG ierror = 0 ! * SPECIFY VARIOUS PARAMETERS FOR HINES ROUTINE AT VERY FIRST CALL. ! * (NOTE THAT ARRAY K_ALPHA IS SPECIFIED SO MAKE SURE THAT ! * IT IS NOT OVERWRITTEN LATER ON). CALL hines_setup(naz, slope, f1, f2, f3, f5, f6, kstar, icutoff, & alt_cutoff, smco, nsmax, iheatcal, k_alpha, ierror, nmessg, klon, nazmth, & coslat) IF (ierror/=0) GO TO 999 ! * START GWD CALCULATIONS. lref = lrefp - 1 DO j = 1, nazmth DO l = 1, klev DO i = kidia, klon sigsqmcw(i, l, j) = 0. END DO END DO END DO ! * INITIALIZE NECESSARY ARRAYS. DO l = 1, klev DO i = kidia, kfdia utendgw(i, l) = 0. vtendgw(i, l) = 0. uhs(i, l) = 0. vhs(i, l) = 0. END DO END DO ! * IF USING HINES SCHEME THEN CALCULATE B V FREQUENCY AT ALL POINTS ! * AND SMOOTH BVFREQ. DO l = 2, klev DO i = kidia, kfdia dttdsf = (th(i, l) / shxkj(i, l) - th(i, l - 1) / shxkj(i, l - 1)) / & (shj(i, l) - shj(i, l - 1)) dttdsf = min(dttdsf, -5. / sgj(i, l)) bvfreq(i, l) = sqrt(-dttdsf * sgj(i, l) * (sgj(i, l)**rgocp) / rd) * rg / ptm1(i, l & ) END DO END DO DO l = 1, klev DO i = kidia, kfdia IF (l==1) THEN bvfreq(i, l) = bvfreq(i, l + 1) END IF IF (l>1) THEN ratio = 5. * log(sgj(i, l) / sgj(i, l - 1)) bvfreq(i, l) = (bvfreq(i, l - 1) + ratio * bvfreq(i, l)) / (1. + ratio) END IF END DO END DO ! * CALCULATE GW DRAG DUE TO HINES' EXTROWAVES ! * SET MOLECULAR VISCOSITY TO A VERY SMALL VALUE. ! * IF THE MODEL TOP IS GREATER THAN 100 KM THEN THE ACTUAL ! * VISCOSITY COEFFICIENT COULD BE SPECIFIED HERE. DO l = 1, klev DO i = kidia, kfdia visc_mol(i, l) = 1.5E-5 drag_u(i, l) = 0. drag_v(i, l) = 0. flux_u(i, l) = 0. flux_v(i, l) = 0. heat(i, l) = 0. diffco(i, l) = 0. END DO END DO ! * ALTITUDE AND DENSITY AT BOTTOM. DO i = kidia, kfdia hscal = rd * ptm1(i, klev) / rg density(i, klev) = sgj(i, klev) * pressg(i) / (rg * hscal) alt(i, klev) = 0. END DO ! * ALTITUDE AND DENSITY AT REMAINING LEVELS. DO l = klev - 1, 1, -1 DO i = kidia, kfdia hscal = rd * ptm1(i, l) / rg alt(i, l) = alt(i, l + 1) + hscal * dsgj(i, l) / sgj(i, l) density(i, l) = sgj(i, l) * pressg(i) / (rg * hscal) END DO END DO ! * INITIALIZE SWITCHES FOR HINES GWD CALCULATION ilrms = 0 DO i = kidia, kfdia lorms(i) = .FALSE. END DO ! * DEFILE BOTTOM LAUNCH LEVEL levbot = iplev ! * BACKGROUND WIND MINUS VALUE AT BOTTOM LAUNCH LEVEL. DO l = 1, levbot DO i = kidia, kfdia uhs(i, l) = pum1(i, l) - pum1(i, levbot) vhs(i, l) = pvm1(i, l) - pvm1(i, levbot) END DO END DO ! * SPECIFY ROOT MEAN SQUARE WIND AT BOTTOM LAUNCH LEVEL. DO i = kidia, kfdia rmswind(i) = rmscon END DO IF (lozpr) THEN DO i = kidia, kfdia IF (zpr(i)>pcrit) THEN rmswind(i) = rmscon + ((zpr(i) - pcrit) / zpr(i)) * pcons END IF END DO END IF DO i = kidia, kfdia IF (rmswind(i)>0.0) THEN ilrms = ilrms + 1 lorms(i) = .TRUE. END IF END DO ! * CALCULATE GWD (NOTE THAT DIFFUSION COEFFICIENT AND ! * HEATING RATE ONLY CALCULATED IF IHEATCAL = 1). IF (ilrms>0) THEN CALL hines_extro0(drag_u, drag_v, heat, diffco, flux_u, flux_v, uhs, vhs, & bvfreq, density, visc_mol, alt, rmswind, k_alpha, m_alpha, v_alpha, & sigma_alpha, sigsqh_alpha, ak_alpha, mmin_alpha, i_alpha, sigma_t, & densbot, bvfbot, 1, iheatcal, icutoff, iprint, nsmax, smco, alt_cutoff, & kstar, slope, f1, f2, f3, f5, f6, naz, sigsqmcw, sigmatm, kidia, klon, & 1, levbot, klon, klev, nazmth, lorms, smoothr1, smoothr2, sigalpmc, & f2mod) ! * ADD ON HINES' GWD TENDENCIES TO OROGRAPHIC TENDENCIES AND ! * APPLY HINES' GW DRAG ON (UROW,VROW) WORK ARRAYS. DO l = 1, klev DO i = kidia, kfdia utendgw(i, l) = utendgw(i, l) + drag_u(i, l) vtendgw(i, l) = vtendgw(i, l) + drag_v(i, l) END DO END DO ! * END OF HINES CALCULATIONS. END IF ! ----------------------------------------------------------------------- DO jl = kidia, kfdia zustrhi(jl) = flux_u(jl, 1) zvstrhi(jl) = flux_v(jl, 1) DO jk = 1, klev le = klev - jk + 1 d_u_hin(jl, jk) = utendgw(jl, le) * dtime d_v_hin(jl, jk) = vtendgw(jl, le) * dtime END DO END DO ! PRINT *,'UTENDGW:',UTENDGW ! PRINT *,' HINES HAS BEEN COMPLETED (LONG ISNT IT...)' RETURN 999 CONTINUE ! * IF ERROR DETECTED THEN ABORT. WRITE (nmessg, 6000) WRITE (nmessg, 6010) ierror 6000 FORMAT (/' EXECUTION ABORTED IN GWDOREXV') 6010 FORMAT (' ERROR FLAG =', I4) END SUBROUTINE hines_gwd ! / ! / SUBROUTINE hines_extro0(drag_u, drag_v, heat, diffco, flux_u, flux_v, vel_u, & vel_v, bvfreq, density, visc_mol, alt, rmswind, k_alpha, m_alpha, & v_alpha, sigma_alpha, sigsqh_alpha, ak_alpha, mmin_alpha, i_alpha, & sigma_t, densb, bvfb, iorder, iheatcal, icutoff, iprint, nsmax, smco, & alt_cutoff, kstar, slope, f1, f2, f3, f5, f6, naz, sigsqmcw, sigmatm, & il1, il2, lev1, lev2, nlons, nlevs, nazmth, lorms, smoothr1, smoothr2, & sigalpmc, f2mod) IMPLICIT NONE ! Main routine for Hines' "extrowave" gravity wave parameterization based ! on Hines' Doppler spread theory. This routine calculates zonal ! and meridional components of gravity wave drag, heating rates ! and diffusion coefficient on a longitude by altitude grid. ! No "mythical" lower boundary region calculation is made so it ! is assumed that lowest level winds are weak (i.e, approximately zero). ! Aug. 13/95 - C. McLandress ! SEPT. /95 - N.McFarlane ! Modifications: ! Output arguements: ! * DRAG_U = zonal component of gravity wave drag (m/s^2). ! * DRAG_V = meridional component of gravity wave drag (m/s^2). ! * HEAT = gravity wave heating (K/sec). ! * DIFFCO = diffusion coefficient (m^2/sec) ! * FLUX_U = zonal component of vertical momentum flux (Pascals) ! * FLUX_V = meridional component of vertical momentum flux (Pascals) ! Input arguements: ! * VEL_U = background zonal wind component (m/s). ! * VEL_V = background meridional wind component (m/s). ! * BVFREQ = background Brunt Vassala frequency (radians/sec). ! * DENSITY = background density (kg/m^3) ! * VISC_MOL = molecular viscosity (m^2/s) ! * ALT = altitude of momentum, density, buoyancy levels (m) ! * (NOTE: levels ordered so that ALT(I,1) > ALT(I,2), etc.) ! * RMSWIND = root mean square gravity wave wind at lowest level (m/s). ! * K_ALPHA = horizontal wavenumber of each azimuth (1/m). ! * IORDER = 1 means vertical levels are indexed from top down ! * (i.e., highest level indexed 1 and lowest level NLEVS); ! * .NE. 1 highest level is index NLEVS. ! * IHEATCAL = 1 to calculate heating rates and diffusion coefficient. ! * IPRINT = 1 to print out various arrays. ! * ICUTOFF = 1 to exponentially damp GWD, heating and diffusion ! * arrays above ALT_CUTOFF; otherwise arrays not modified. ! * ALT_CUTOFF = altitude in meters above which exponential decay applied. ! * SMCO = smoothing factor used to smooth cutoff vertical ! * wavenumbers and total rms winds in vertical direction ! * before calculating drag or heating ! * (SMCO >= 1 ==> 1:SMCO:1 stencil used). ! * NSMAX = number of times smoother applied ( >= 1), ! * = 0 means no smoothing performed. ! * KSTAR = typical gravity wave horizontal wavenumber (1/m). ! * SLOPE = slope of incident vertical wavenumber spectrum ! * (SLOPE must equal 1., 1.5 or 2.). ! * F1 to F6 = Hines's fudge factors (F4 not needed since used for ! * vertical flux of vertical momentum). ! * NAZ = actual number of horizontal azimuths used. ! * IL1 = first longitudinal index to use (IL1 >= 1). ! * IL2 = last longitudinal index to use (IL1 <= IL2 <= NLONS). ! * LEV1 = index of first level for drag calculation. ! * LEV2 = index of last level for drag calculation ! * (i.e., LEV1 < LEV2 <= NLEVS). ! * NLONS = number of longitudes. ! * NLEVS = number of vertical levels. ! * NAZMTH = azimuthal array dimension (NAZMTH >= NAZ). ! Work arrays. ! * M_ALPHA = cutoff vertical wavenumber (1/m). ! * V_ALPHA = wind component at each azimuth (m/s) and if IHEATCAL=1 ! * holds vertical derivative of cutoff wavenumber. ! * SIGMA_ALPHA = total rms wind in each azimuth (m/s). ! * SIGSQH_ALPHA = portion of wind variance from waves having wave ! * normals in the alpha azimuth (m/s). ! * SIGMA_T = total rms horizontal wind (m/s). ! * AK_ALPHA = spectral amplitude factor at each azimuth ! * (i.e.,{AjKj}) in m^4/s^2. ! * I_ALPHA = Hines' integral. ! * MMIN_ALPHA = minimum value of cutoff wavenumber. ! * DENSB = background density at bottom level. ! * BVFB = buoyancy frequency at bottom level and ! * work array for ICUTOFF = 1. ! * LORMS = .TRUE. for drag computation INTEGER naz, nlons, nlevs, nazmth, il1, il2, lev1, lev2 INTEGER icutoff, nsmax, iorder, iheatcal, iprint REAL kstar(nlons), f1, f2, f3, f5, f6, slope REAL alt_cutoff, smco REAL drag_u(nlons, nlevs), drag_v(nlons, nlevs) REAL heat(nlons, nlevs), diffco(nlons, nlevs) REAL flux_u(nlons, nlevs), flux_v(nlons, nlevs) REAL vel_u(nlons, nlevs), vel_v(nlons, nlevs) REAL bvfreq(nlons, nlevs), density(nlons, nlevs) REAL visc_mol(nlons, nlevs), alt(nlons, nlevs) REAL rmswind(nlons), bvfb(nlons), densb(nlons) REAL sigma_t(nlons, nlevs), sigsqmcw(nlons, nlevs, nazmth) REAL sigma_alpha(nlons, nlevs, nazmth), sigmatm(nlons, nlevs) REAL sigsqh_alpha(nlons, nlevs, nazmth) REAL m_alpha(nlons, nlevs, nazmth), v_alpha(nlons, nlevs, nazmth) REAL ak_alpha(nlons, nazmth), k_alpha(nlons, nazmth) REAL mmin_alpha(nlons, nazmth), i_alpha(nlons, nazmth) REAL smoothr1(nlons, nlevs), smoothr2(nlons, nlevs) REAL sigalpmc(nlons, nlevs, nazmth) REAL f2mod(nlons, nlevs) LOGICAL lorms(nlons) ! Internal variables. INTEGER levbot, levtop, i, n, l, lev1p, lev2m INTEGER ilprt1, ilprt2 ! ----------------------------------------------------------------------- ! PRINT *,' IN HINES_EXTRO0' lev1p = lev1 + 1 lev2m = lev2 - 1 ! Index of lowest altitude level (bottom of drag calculation). levbot = lev2 levtop = lev1 IF (iorder/=1) THEN WRITE (6, 1) 1 FORMAT (2X, ' error: IORDER NOT ONE! ') END IF ! Buoyancy and density at bottom level. DO i = il1, il2 bvfb(i) = bvfreq(i, levbot) densb(i) = density(i, levbot) END DO ! initialize some variables DO n = 1, naz DO l = lev1, lev2 DO i = il1, il2 m_alpha(i, l, n) = 0.0 END DO END DO END DO DO l = lev1, lev2 DO i = il1, il2 sigma_t(i, l) = 0.0 END DO END DO DO n = 1, naz DO i = il1, il2 i_alpha(i, n) = 0.0 END DO END DO ! Compute azimuthal wind components from zonal and meridional winds. CALL hines_wind(v_alpha, vel_u, vel_v, naz, il1, il2, lev1, lev2, nlons, & nlevs, nazmth) ! Calculate cutoff vertical wavenumber and velocity variances. CALL hines_wavnum(m_alpha, sigma_alpha, sigsqh_alpha, sigma_t, ak_alpha, & v_alpha, visc_mol, density, densb, bvfreq, bvfb, rmswind, i_alpha, & mmin_alpha, kstar, slope, f1, f2, f3, naz, levbot, levtop, il1, il2, & nlons, nlevs, nazmth, sigsqmcw, sigmatm, lorms, sigalpmc, f2mod) ! Smooth cutoff wavenumbers and total rms velocity in the vertical ! direction NSMAX times, using FLUX_U as temporary work array. IF (nsmax>0) THEN DO n = 1, naz DO l = lev1, lev2 DO i = il1, il2 smoothr1(i, l) = m_alpha(i, l, n) END DO END DO CALL vert_smooth(smoothr1, smoothr2, smco, nsmax, il1, il2, lev1, lev2, & nlons, nlevs) DO l = lev1, lev2 DO i = il1, il2 m_alpha(i, l, n) = smoothr1(i, l) END DO END DO END DO CALL vert_smooth(sigma_t, smoothr2, smco, nsmax, il1, il2, lev1, lev2, & nlons, nlevs) END IF ! Calculate zonal and meridional components of the ! momentum flux and drag. CALL hines_flux(flux_u, flux_v, drag_u, drag_v, alt, density, densb, & m_alpha, ak_alpha, k_alpha, slope, naz, il1, il2, lev1, lev2, nlons, & nlevs, nazmth, lorms) ! Cutoff drag above ALT_CUTOFF, using BVFB as temporary work array. IF (icutoff==1) THEN CALL hines_exp(drag_u, bvfb, alt, alt_cutoff, iorder, il1, il2, lev1, & lev2, nlons, nlevs) CALL hines_exp(drag_v, bvfb, alt, alt_cutoff, iorder, il1, il2, lev1, & lev2, nlons, nlevs) END IF ! Print out various arrays for diagnostic purposes. IF (iprint==1) THEN ilprt1 = 15 ilprt2 = 16 CALL hines_print(flux_u, flux_v, drag_u, drag_v, alt, sigma_t, & sigma_alpha, v_alpha, m_alpha, 1, 1, 6, ilprt1, ilprt2, lev1, lev2, & naz, nlons, nlevs, nazmth) END IF ! If not calculating heating rate and diffusion coefficient then finished. IF (iheatcal/=1) RETURN ! Calculate vertical derivative of cutoff wavenumber (store ! in array V_ALPHA) using centered differences at interior gridpoints ! and one-sided differences at first and last levels. DO n = 1, naz DO l = lev1p, lev2m DO i = il1, il2 v_alpha(i, l, n) = (m_alpha(i, l + 1, n) - m_alpha(i, l - 1, n)) / & (alt(i, l + 1) - alt(i, l - 1)) END DO END DO DO i = il1, il2 v_alpha(i, lev1, n) = (m_alpha(i, lev1p, n) - m_alpha(i, lev1, n)) / & (alt(i, lev1p) - alt(i, lev1)) END DO DO i = il1, il2 v_alpha(i, lev2, n) = (m_alpha(i, lev2, n) - m_alpha(i, lev2m, n)) / & (alt(i, lev2) - alt(i, lev2m)) END DO END DO ! Heating rate and diffusion coefficient. CALL hines_heat(heat, diffco, m_alpha, v_alpha, ak_alpha, k_alpha, bvfreq, & density, densb, sigma_t, visc_mol, kstar, slope, f2, f3, f5, f6, naz, & il1, il2, lev1, lev2, nlons, nlevs, nazmth) ! Finished. RETURN ! ----------------------------------------------------------------------- END SUBROUTINE hines_extro0 SUBROUTINE hines_wavnum(m_alpha, sigma_alpha, sigsqh_alpha, sigma_t, & ak_alpha, v_alpha, visc_mol, density, densb, bvfreq, bvfb, rms_wind, & i_alpha, mmin_alpha, kstar, slope, f1, f2, f3, naz, levbot, levtop, il1, & il2, nlons, nlevs, nazmth, sigsqmcw, sigmatm, lorms, sigalpmc, f2mod) IMPLICIT NONE ! This routine calculates the cutoff vertical wavenumber and velocity ! variances on a longitude by altitude grid for the Hines' Doppler ! spread gravity wave drag parameterization scheme. ! NOTE: (1) only values of four or eight can be used for # azimuths (NAZ). ! (2) only values of 1.0, 1.5 or 2.0 can be used for slope (SLOPE). ! Aug. 10/95 - C. McLandress ! Output arguements: ! * M_ALPHA = cutoff wavenumber at each azimuth (1/m). ! * SIGMA_ALPHA = total rms wind in each azimuth (m/s). ! * SIGSQH_ALPHA = portion of wind variance from waves having wave ! * normals in the alpha azimuth (m/s). ! * SIGMA_T = total rms horizontal wind (m/s). ! * AK_ALPHA = spectral amplitude factor at each azimuth ! * (i.e.,{AjKj}) in m^4/s^2. ! Input arguements: ! * V_ALPHA = wind component at each azimuth (m/s). ! * VISC_MOL = molecular viscosity (m^2/s) ! * DENSITY = background density (kg/m^3). ! * DENSB = background density at model bottom (kg/m^3). ! * BVFREQ = background Brunt Vassala frequency (radians/sec). ! * BVFB = background Brunt Vassala frequency at model bottom. ! * RMS_WIND = root mean square gravity wave wind at lowest level (m/s). ! * KSTAR = typical gravity wave horizontal wavenumber (1/m). ! * SLOPE = slope of incident vertical wavenumber spectrum ! * (SLOPE = 1., 1.5 or 2.). ! * F1,F2,F3 = Hines's fudge factors. ! * NAZ = actual number of horizontal azimuths used (4 or 8). ! * LEVBOT = index of lowest vertical level. ! * LEVTOP = index of highest vertical level ! * (NOTE: if LEVTOP < LEVBOT then level index ! * increases from top down). ! * IL1 = first longitudinal index to use (IL1 >= 1). ! * IL2 = last longitudinal index to use (IL1 <= IL2 <= NLONS). ! * NLONS = number of longitudes. ! * NLEVS = number of vertical levels. ! * NAZMTH = azimuthal array dimension (NAZMTH >= NAZ). ! * LORMS = .TRUE. for drag computation ! Input work arrays: ! * I_ALPHA = Hines' integral at a single level. ! * MMIN_ALPHA = minimum value of cutoff wavenumber. INTEGER naz, levbot, levtop, il1, il2, nlons, nlevs, nazmth REAL slope, kstar(nlons), f1, f2, f3, f2mfac REAL m_alpha(nlons, nlevs, nazmth) REAL sigma_alpha(nlons, nlevs, nazmth) REAL sigalpmc(nlons, nlevs, nazmth) REAL sigsqh_alpha(nlons, nlevs, nazmth) REAL sigsqmcw(nlons, nlevs, nazmth) REAL sigma_t(nlons, nlevs) REAL sigmatm(nlons, nlevs) REAL ak_alpha(nlons, nazmth) REAL v_alpha(nlons, nlevs, nazmth) REAL visc_mol(nlons, nlevs) REAL f2mod(nlons, nlevs) REAL density(nlons, nlevs), densb(nlons) REAL bvfreq(nlons, nlevs), bvfb(nlons), rms_wind(nlons) REAL i_alpha(nlons, nazmth), mmin_alpha(nlons, nazmth) LOGICAL lorms(nlons) ! Internal variables. INTEGER i, l, n, lstart, lend, lincr, lbelow REAL m_sub_m_turb, m_sub_m_mol, m_trial REAL visc, visc_min, azfac, sp1 ! c REAL N_OVER_M(1000), SIGFAC(1000) REAL n_over_m(nlons), sigfac(nlons) DATA visc_min/1.E-10/ ! ----------------------------------------------------------------------- ! PRINT *,'IN HINES_WAVNUM' sp1 = slope + 1. ! Indices of levels to process. IF (levbot>levtop) THEN lstart = levbot - 1 lend = levtop lincr = -1 ELSE WRITE (6, 1) 1 FORMAT (2X, ' error: IORDER NOT ONE! ') END IF ! Use horizontal isotropy to calculate azimuthal variances at bottom level. azfac = 1. / real(naz) DO n = 1, naz DO i = il1, il2 sigsqh_alpha(i, levbot, n) = azfac * rms_wind(i)**2 END DO END DO ! Velocity variances at bottom level. CALL hines_sigma(sigma_t, sigma_alpha, sigsqh_alpha, naz, levbot, il1, il2, & nlons, nlevs, nazmth) CALL hines_sigma(sigmatm, sigalpmc, sigsqmcw, naz, levbot, il1, il2, nlons, & nlevs, nazmth) ! Calculate cutoff wavenumber and spectral amplitude factor ! at bottom level where it is assumed that background winds vanish ! and also initialize minimum value of cutoff wavnumber. DO n = 1, naz DO i = il1, il2 IF (lorms(i)) THEN m_alpha(i, levbot, n) = bvfb(i) / (f1 * sigma_alpha(i, levbot, n) + f2 * & sigma_t(i, levbot)) ak_alpha(i, n) = sigsqh_alpha(i, levbot, n) / & (m_alpha(i, levbot, n)**sp1 / sp1) mmin_alpha(i, n) = m_alpha(i, levbot, n) END IF END DO END DO ! Calculate quantities from the bottom upwards, ! starting one level above bottom. DO l = lstart, lend, lincr ! Level beneath present level. lbelow = l - lincr ! Calculate N/m_M where m_M is maximum permissible value of the vertical ! wavenumber (i.e., m > m_M are obliterated) and N is buoyancy frequency. ! m_M is taken as the smaller of the instability-induced ! wavenumber (M_SUB_M_TURB) and that imposed by molecular viscosity ! (M_SUB_M_MOL). Since variance at this level is not yet known ! use value at level below. DO i = il1, il2 IF (lorms(i)) THEN f2mfac = sigmatm(i, lbelow)**2 f2mod(i, lbelow) = 1. + 2. * f2mfac / (f2mfac + sigma_t(i, lbelow)**2) visc = amax1(visc_mol(i, l), visc_min) m_sub_m_turb = bvfreq(i, l) / (f2 * f2mod(i, lbelow) * sigma_t(i, lbelow)) m_sub_m_mol = (bvfreq(i, l) * kstar(i) / visc)**0.33333333 / f3 IF (m_sub_m_turbmmin_alpha(i, n)) THEN m_trial = mmin_alpha(i, n) END IF m_alpha(i, l, n) = m_trial ! Reset minimum value of cutoff wavenumber if necessary. IF (m_alpha(i, l, n)= 1). ! * IL2 = last longitudinal index to use (IL1 <= IL2 <= NLONS). ! * LEV1 = first altitude level to use (LEV1 >=1). ! * LEV2 = last altitude level to use (LEV1 < LEV2 <= NLEVS). ! * NLONS = number of longitudes. ! * NLEVS = number of vertical levels. ! * NAZMTH = azimuthal array dimension (NAZMTH >= NAZ). ! Constants in DATA statements. ! * COS45 = cosine of 45 degrees. ! * UMIN = minimum allowable value for zonal or meridional ! * wind component (m/s). ! Subroutine arguements. INTEGER naz, il1, il2, lev1, lev2 INTEGER nlons, nlevs, nazmth REAL v_alpha(nlons, nlevs, nazmth) REAL vel_u(nlons, nlevs), vel_v(nlons, nlevs) ! Internal variables. INTEGER i, l REAL u, v, cos45, umin DATA cos45/0.7071068/ DATA umin/0.001/ ! ----------------------------------------------------------------------- ! Case with 4 azimuths. ! PRINT *,'IN HINES_WIND' IF (naz==4) THEN DO l = lev1, lev2 DO i = il1, il2 u = vel_u(i, l) v = vel_v(i, l) IF (abs(u)= 1). ! * IL2 = last longitudinal index to use (IL1 <= IL2 <= NLONS). ! * LEV1 = first altitude level to use (LEV1 >=1). ! * LEV2 = last altitude level to use (LEV1 < LEV2 <= NLEVS). ! * NLONS = number of longitudes. ! * NLEVS = number of vertical levels. ! * NAZMTH = azimuthal array dimension (NAZMTH >= NAZ). ! * LORMS = .TRUE. for drag computation ! Constant in DATA statement. ! * COS45 = cosine of 45 degrees. ! Subroutine arguements. INTEGER naz, il1, il2, lev1, lev2 INTEGER nlons, nlevs, nazmth REAL slope REAL flux_u(nlons, nlevs), flux_v(nlons, nlevs) REAL drag_u(nlons, nlevs), drag_v(nlons, nlevs) REAL alt(nlons, nlevs), density(nlons, nlevs), densb(nlons) REAL m_alpha(nlons, nlevs, nazmth) REAL ak_alpha(nlons, nazmth), k_alpha(nlons, nazmth) LOGICAL lorms(nlons) ! Internal variables. INTEGER i, l, lev1p, lev2m, lev2p REAL cos45, prod2, prod4, prod6, prod8, dendz, dendz2 DATA cos45/0.7071068/ ! ----------------------------------------------------------------------- lev1p = lev1 + 1 lev2m = lev2 - 1 lev2p = lev2 + 1 ! Sum over azimuths for case where SLOPE = 1. IF (slope==1.) THEN ! Case with 4 azimuths. IF (naz==4) THEN DO l = lev1, lev2 DO i = il1, il2 flux_u(i, l) = ak_alpha(i, 1) * k_alpha(i, 1) * m_alpha(i, l, 1) - & ak_alpha(i, 3) * k_alpha(i, 3) * m_alpha(i, l, 3) flux_v(i, l) = ak_alpha(i, 2) * k_alpha(i, 2) * m_alpha(i, l, 2) - & ak_alpha(i, 4) * k_alpha(i, 4) * m_alpha(i, l, 4) END DO END DO END IF ! Case with 8 azimuths. IF (naz==8) THEN DO l = lev1, lev2 DO i = il1, il2 prod2 = ak_alpha(i, 2) * k_alpha(i, 2) * m_alpha(i, l, 2) prod4 = ak_alpha(i, 4) * k_alpha(i, 4) * m_alpha(i, l, 4) prod6 = ak_alpha(i, 6) * k_alpha(i, 6) * m_alpha(i, l, 6) prod8 = ak_alpha(i, 8) * k_alpha(i, 8) * m_alpha(i, l, 8) flux_u(i, l) = ak_alpha(i, 1) * k_alpha(i, 1) * m_alpha(i, l, 1) - & ak_alpha(i, 5) * k_alpha(i, 5) * m_alpha(i, l, 5) + & cos45 * (prod2 - prod4 - prod6 + prod8) flux_v(i, l) = ak_alpha(i, 3) * k_alpha(i, 3) * m_alpha(i, l, 3) - & ak_alpha(i, 7) * k_alpha(i, 7) * m_alpha(i, l, 7) + & cos45 * (prod2 + prod4 - prod6 - prod8) END DO END DO END IF END IF ! Sum over azimuths for case where SLOPE not equal to 1. IF (slope/=1.) THEN ! Case with 4 azimuths. IF (naz==4) THEN DO l = lev1, lev2 DO i = il1, il2 flux_u(i, l) = ak_alpha(i, 1) * k_alpha(i, 1) * & m_alpha(i, l, 1)**slope - ak_alpha(i, 3) * k_alpha(i, 3) * m_alpha(i, & l, 3)**slope flux_v(i, l) = ak_alpha(i, 2) * k_alpha(i, 2) * & m_alpha(i, l, 2)**slope - ak_alpha(i, 4) * k_alpha(i, 4) * m_alpha(i, & l, 4)**slope END DO END DO END IF ! Case with 8 azimuths. IF (naz==8) THEN DO l = lev1, lev2 DO i = il1, il2 prod2 = ak_alpha(i, 2) * k_alpha(i, 2) * m_alpha(i, l, 2)**slope prod4 = ak_alpha(i, 4) * k_alpha(i, 4) * m_alpha(i, l, 4)**slope prod6 = ak_alpha(i, 6) * k_alpha(i, 6) * m_alpha(i, l, 6)**slope prod8 = ak_alpha(i, 8) * k_alpha(i, 8) * m_alpha(i, l, 8)**slope flux_u(i, l) = ak_alpha(i, 1) * k_alpha(i, 1) * & m_alpha(i, l, 1)**slope - ak_alpha(i, 5) * k_alpha(i, 5) * m_alpha(i, & l, 5)**slope + cos45 * (prod2 - prod4 - prod6 + prod8) flux_v(i, l) = ak_alpha(i, 3) * k_alpha(i, 3) * & m_alpha(i, l, 3)**slope - ak_alpha(i, 7) * k_alpha(i, 7) * m_alpha(i, & l, 7)**slope + cos45 * (prod2 + prod4 - prod6 - prod8) END DO END DO END IF END IF ! Calculate flux from sum. DO l = lev1, lev2 DO i = il1, il2 flux_u(i, l) = flux_u(i, l) * densb(i) / slope flux_v(i, l) = flux_v(i, l) * densb(i) / slope END DO END DO ! Calculate drag at intermediate levels using centered differences DO l = lev1p, lev2m DO i = il1, il2 IF (lorms(i)) THEN ! cc DENDZ2 = DENSITY(I,L) * ( ALT(I,L+1) - ALT(I,L-1) ) dendz2 = density(i, l) * (alt(i, l - 1) - alt(i, l)) ! cc DRAG_U(I,L) = - ( FLUX_U(I,L+1) - FLUX_U(I,L-1) ) / DENDZ2 drag_u(i, l) = -(flux_u(i, l - 1) - flux_u(i, l)) / dendz2 ! cc DRAG_V(I,L) = - ( FLUX_V(I,L+1) - FLUX_V(I,L-1) ) / DENDZ2 drag_v(i, l) = -(flux_v(i, l - 1) - flux_v(i, l)) / dendz2 END IF END DO END DO ! Drag at first and last levels using one-side differences. DO i = il1, il2 IF (lorms(i)) THEN dendz = density(i, lev1) * (alt(i, lev1) - alt(i, lev1p)) drag_u(i, lev1) = flux_u(i, lev1) / dendz drag_v(i, lev1) = flux_v(i, lev1) / dendz END IF END DO DO i = il1, il2 IF (lorms(i)) THEN dendz = density(i, lev2) * (alt(i, lev2m) - alt(i, lev2)) drag_u(i, lev2) = -(flux_u(i, lev2m) - flux_u(i, lev2)) / dendz drag_v(i, lev2) = -(flux_v(i, lev2m) - flux_v(i, lev2)) / dendz END IF END DO IF (nlevs>lev2) THEN DO i = il1, il2 IF (lorms(i)) THEN dendz = density(i, lev2p) * (alt(i, lev2) - alt(i, lev2p)) drag_u(i, lev2p) = -flux_u(i, lev2) / dendz drag_v(i, lev2p) = -flux_v(i, lev2) / dendz END IF END DO END IF RETURN ! ----------------------------------------------------------------------- END SUBROUTINE hines_flux SUBROUTINE hines_heat(heat, diffco, m_alpha, dmdz_alpha, ak_alpha, k_alpha, & bvfreq, density, densb, sigma_t, visc_mol, kstar, slope, f2, f3, f5, f6, & naz, il1, il2, lev1, lev2, nlons, nlevs, nazmth) IMPLICIT NONE ! This routine calculates the gravity wave induced heating and ! diffusion coefficient on a longitude by altitude grid for ! the Hines' Doppler spread gravity wave drag parameterization scheme. ! Aug. 6/95 - C. McLandress ! Output arguements: ! * HEAT = gravity wave heating (K/sec). ! * DIFFCO = diffusion coefficient (m^2/sec) ! Input arguements: ! * M_ALPHA = cutoff vertical wavenumber (1/m). ! * DMDZ_ALPHA = vertical derivative of cutoff wavenumber. ! * AK_ALPHA = spectral amplitude factor of each azimuth ! (i.e., {AjKj} in m^4/s^2). ! * K_ALPHA = horizontal wavenumber of each azimuth (1/m). ! * BVFREQ = background Brunt Vassala frequency (rad/sec). ! * DENSITY = background density (kg/m^3). ! * DENSB = background density at bottom level (kg/m^3). ! * SIGMA_T = total rms horizontal wind (m/s). ! * VISC_MOL = molecular viscosity (m^2/s). ! * KSTAR = typical gravity wave horizontal wavenumber (1/m). ! * SLOPE = slope of incident vertical wavenumber spectrum. ! * F2,F3,F5,F6 = Hines's fudge factors. ! * NAZ = actual number of horizontal azimuths used. ! * IL1 = first longitudinal index to use (IL1 >= 1). ! * IL2 = last longitudinal index to use (IL1 <= IL2 <= NLONS). ! * LEV1 = first altitude level to use (LEV1 >=1). ! * LEV2 = last altitude level to use (LEV1 < LEV2 <= NLEVS). ! * NLONS = number of longitudes. ! * NLEVS = number of vertical levels. ! * NAZMTH = azimuthal array dimension (NAZMTH >= NAZ). INTEGER naz, il1, il2, lev1, lev2, nlons, nlevs, nazmth REAL kstar(nlons), slope, f2, f3, f5, f6 REAL heat(nlons, nlevs), diffco(nlons, nlevs) REAL m_alpha(nlons, nlevs, nazmth), dmdz_alpha(nlons, nlevs, nazmth) REAL ak_alpha(nlons, nazmth), k_alpha(nlons, nazmth) REAL bvfreq(nlons, nlevs), density(nlons, nlevs), densb(nlons) REAL sigma_t(nlons, nlevs), visc_mol(nlons, nlevs) ! Internal variables. INTEGER i, l, n REAL m_sub_m_turb, m_sub_m_mol, m_sub_m, heatng REAL visc, visc_min, cpgas, sm1 ! specific heat at constant pressure DATA cpgas/1004./ ! minimum permissible viscosity DATA visc_min/1.E-10/ ! ----------------------------------------------------------------------- ! Initialize heating array. DO l = 1, nlevs DO i = 1, nlons heat(i, l) = 0. END DO END DO ! Perform sum over azimuths for case where SLOPE = 1. IF (slope==1.) THEN DO n = 1, naz DO l = lev1, lev2 DO i = il1, il2 heat(i, l) = heat(i, l) + ak_alpha(i, n) * k_alpha(i, n) * dmdz_alpha(i & , l, n) END DO END DO END DO END IF ! Perform sum over azimuths for case where SLOPE not 1. IF (slope/=1.) THEN sm1 = slope - 1. DO n = 1, naz DO l = lev1, lev2 DO i = il1, il2 heat(i, l) = heat(i, l) + ak_alpha(i, n) * k_alpha(i, n) * m_alpha(i, l & , n)**sm1 * dmdz_alpha(i, l, n) END DO END DO END DO END IF ! Heating and diffusion. DO l = lev1, lev2 DO i = il1, il2 ! Maximum permissible value of cutoff wavenumber is the smaller ! of the instability-induced wavenumber (M_SUB_M_TURB) and ! that imposed by molecular viscosity (M_SUB_M_MOL). visc = amax1(visc_mol(i, l), visc_min) m_sub_m_turb = bvfreq(i, l) / (f2 * sigma_t(i, l)) m_sub_m_mol = (bvfreq(i, l) * kstar(i) / visc)**0.33333333 / f3 m_sub_m = amin1(m_sub_m_turb, m_sub_m_mol) heatng = -heat(i, l) * f5 * bvfreq(i, l) / m_sub_m * densb(i) / density(i, l) diffco(i, l) = f6 * heatng**0.33333333 / m_sub_m**1.33333333 heat(i, l) = heatng / cpgas END DO END DO RETURN ! ----------------------------------------------------------------------- END SUBROUTINE hines_heat SUBROUTINE hines_sigma(sigma_t, sigma_alpha, sigsqh_alpha, naz, lev, il1, & il2, nlons, nlevs, nazmth) IMPLICIT NONE ! This routine calculates the total rms and azimuthal rms horizontal ! velocities at a given level on a longitude by altitude grid for ! the Hines' Doppler spread GWD parameterization scheme. ! NOTE: only four or eight azimuths can be used. ! Aug. 7/95 - C. McLandress ! Output arguements: ! * SIGMA_T = total rms horizontal wind (m/s). ! * SIGMA_ALPHA = total rms wind in each azimuth (m/s). ! Input arguements: ! * SIGSQH_ALPHA = portion of wind variance from waves having wave ! * normals in the alpha azimuth (m/s). ! * NAZ = actual number of horizontal azimuths used (must be 4 or 8). ! * LEV = altitude level to process. ! * IL1 = first longitudinal index to use (IL1 >= 1). ! * IL2 = last longitudinal index to use (IL1 <= IL2 <= NLONS). ! * NLONS = number of longitudes. ! * NLEVS = number of vertical levels. ! * NAZMTH = azimuthal array dimension (NAZMTH >= NAZ). ! Subroutine arguements. INTEGER lev, naz, il1, il2 INTEGER nlons, nlevs, nazmth REAL sigma_t(nlons, nlevs) REAL sigma_alpha(nlons, nlevs, nazmth) REAL sigsqh_alpha(nlons, nlevs, nazmth) ! Internal variables. INTEGER i, n REAL sum_even, sum_odd ! ----------------------------------------------------------------------- ! Calculate azimuthal rms velocity for the 4 azimuth case. IF (naz==4) THEN DO i = il1, il2 sigma_alpha(i, lev, 1) = sqrt(sigsqh_alpha(i, lev, 1) + sigsqh_alpha(i, lev, & 3)) sigma_alpha(i, lev, 2) = sqrt(sigsqh_alpha(i, lev, 2) + sigsqh_alpha(i, lev, & 4)) sigma_alpha(i, lev, 3) = sigma_alpha(i, lev, 1) sigma_alpha(i, lev, 4) = sigma_alpha(i, lev, 2) END DO END IF ! Calculate azimuthal rms velocity for the 8 azimuth case. IF (naz==8) THEN DO i = il1, il2 sum_odd = (sigsqh_alpha(i, lev, 1) + sigsqh_alpha(i, lev, 3) + & sigsqh_alpha(i, lev, 5) + sigsqh_alpha(i, lev, 7)) / 2. sum_even = (sigsqh_alpha(i, lev, 2) + sigsqh_alpha(i, lev, 4) + & sigsqh_alpha(i, lev, 6) + sigsqh_alpha(i, lev, 8)) / 2. sigma_alpha(i, lev, 1) = sqrt(sigsqh_alpha(i, lev, 1) + sigsqh_alpha(i, lev, & 5) + sum_even) sigma_alpha(i, lev, 2) = sqrt(sigsqh_alpha(i, lev, 2) + sigsqh_alpha(i, lev, & 6) + sum_odd) sigma_alpha(i, lev, 3) = sqrt(sigsqh_alpha(i, lev, 3) + sigsqh_alpha(i, lev, & 7) + sum_even) sigma_alpha(i, lev, 4) = sqrt(sigsqh_alpha(i, lev, 4) + sigsqh_alpha(i, lev, & 8) + sum_odd) sigma_alpha(i, lev, 5) = sigma_alpha(i, lev, 1) sigma_alpha(i, lev, 6) = sigma_alpha(i, lev, 2) sigma_alpha(i, lev, 7) = sigma_alpha(i, lev, 3) sigma_alpha(i, lev, 8) = sigma_alpha(i, lev, 4) END DO END IF ! Calculate total rms velocity. DO i = il1, il2 sigma_t(i, lev) = 0. END DO DO n = 1, naz DO i = il1, il2 sigma_t(i, lev) = sigma_t(i, lev) + sigsqh_alpha(i, lev, n) END DO END DO DO i = il1, il2 sigma_t(i, lev) = sqrt(sigma_t(i, lev)) END DO RETURN ! ----------------------------------------------------------------------- END SUBROUTINE hines_sigma SUBROUTINE hines_intgrl(i_alpha, v_alpha, m_alpha, bvfb, slope, naz, lev, & il1, il2, nlons, nlevs, nazmth, lorms) IMPLICIT NONE ! This routine calculates the vertical wavenumber integral ! for a single vertical level at each azimuth on a longitude grid ! for the Hines' Doppler spread GWD parameterization scheme. ! NOTE: (1) only spectral slopes of 1, 1.5 or 2 are permitted. ! (2) the integral is written in terms of the product QM ! which by construction is always less than 1. Series ! solutions are used for small |QM| and analytical solutions ! for remaining values. ! Aug. 8/95 - C. McLandress ! Output arguement: ! * I_ALPHA = Hines' integral. ! Input arguements: ! * V_ALPHA = azimuthal wind component (m/s). ! * M_ALPHA = azimuthal cutoff vertical wavenumber (1/m). ! * BVFB = background Brunt Vassala frequency at model bottom. ! * SLOPE = slope of initial vertical wavenumber spectrum ! * (must use SLOPE = 1., 1.5 or 2.) ! * NAZ = actual number of horizontal azimuths used. ! * LEV = altitude level to process. ! * IL1 = first longitudinal index to use (IL1 >= 1). ! * IL2 = last longitudinal index to use (IL1 <= IL2 <= NLONS). ! * NLONS = number of longitudes. ! * NLEVS = number of vertical levels. ! * NAZMTH = azimuthal array dimension (NAZMTH >= NAZ). ! * LORMS = .TRUE. for drag computation ! Constants in DATA statements: ! * QMIN = minimum value of Q_ALPHA (avoids indeterminant form of integral) ! * QM_MIN = minimum value of Q_ALPHA * M_ALPHA (used to avoid numerical ! * problems). INTEGER lev, naz, il1, il2, nlons, nlevs, nazmth REAL i_alpha(nlons, nazmth) REAL v_alpha(nlons, nlevs, nazmth) REAL m_alpha(nlons, nlevs, nazmth) REAL bvfb(nlons), slope LOGICAL lorms(nlons) ! Internal variables. INTEGER i, n REAL q_alpha, qm, sqrtqm, q_min, qm_min DATA q_min/1.0/, qm_min/0.01/ ! ----------------------------------------------------------------------- ! For integer value SLOPE = 1. IF (slope==1.) THEN DO n = 1, naz DO i = il1, il2 IF (lorms(i)) THEN q_alpha = v_alpha(i, lev, n) / bvfb(i) qm = q_alpha * m_alpha(i, lev, n) ! If |QM| is small then use first 4 terms series of Taylor series ! expansion of integral in order to avoid indeterminate form of ! integral, ! otherwise use analytical form of integral. IF (abs(q_alpha)=0.) THEN i_alpha(i, n) = (alog((1. + sqrtqm) / (1. - sqrtqm)) - 2. * sqrtqm * (1. + qm & / 3.)) / q_alpha**2.5 ELSE i_alpha(i, n) = 2. * (atan(sqrtqm) + sqrtqm * (qm / 3. - 1.)) / & abs(q_alpha)**2.5 END IF END IF END IF END DO END DO END IF ! If integral is negative (which in principal should not happen) THEN ! print a message and some info since execution will abort when calculating ! the variances. ! DO 80 N = 1,NAZ ! DO 70 I = IL1,IL2 ! IF (I_ALPHA(I,N).LT.0.) THEN ! WRITE (6,*) ! WRITE (6,*) '******************************' ! WRITE (6,*) 'Hines integral I_ALPHA < 0 ' ! WRITE (6,*) ' longitude I=',I ! WRITE (6,*) ' azimuth N=',N ! WRITE (6,*) ' level LEV=',LEV ! WRITE (6,*) ' I_ALPHA =',I_ALPHA(I,N) ! WRITE (6,*) ' V_ALPHA =',V_ALPHA(I,LEV,N) ! WRITE (6,*) ' M_ALPHA =',M_ALPHA(I,LEV,N) ! WRITE (6,*) ' Q_ALPHA =',V_ALPHA(I,LEV,N) / BVFB(I) ! WRITE (6,*) ' QM =',V_ALPHA(I,LEV,N) / BVFB(I) ! ^ * M_ALPHA(I,LEV,N) ! WRITE (6,*) '******************************' ! END IF ! 70 CONTINUE ! 80 CONTINUE RETURN ! ----------------------------------------------------------------------- END SUBROUTINE hines_intgrl SUBROUTINE hines_setup(naz, slope, f1, f2, f3, f5, f6, kstar, icutoff, & alt_cutoff, smco, nsmax, iheatcal, k_alpha, ierror, nmessg, nlons, & nazmth, coslat) IMPLICIT NONE ! This routine specifies various parameters needed for the ! the Hines' Doppler spread gravity wave drag parameterization scheme. ! Aug. 8/95 - C. McLandress ! Output arguements: ! * NAZ = actual number of horizontal azimuths used ! * (code set up presently for only NAZ = 4 or 8). ! * SLOPE = slope of incident vertical wavenumber spectrum ! * (code set up presently for SLOPE 1., 1.5 or 2.). ! * F1 = "fudge factor" used in calculation of trial value of ! * azimuthal cutoff wavenumber M_ALPHA (1.2 <= F1 <= 1.9). ! * F2 = "fudge factor" used in calculation of maximum ! * permissible instabiliy-induced cutoff wavenumber ! * M_SUB_M_TURB (0.1 <= F2 <= 1.4). ! * F3 = "fudge factor" used in calculation of maximum ! * permissible molecular viscosity-induced cutoff wavenumber ! * M_SUB_M_MOL (0.1 <= F2 <= 1.4). ! * F5 = "fudge factor" used in calculation of heating rate ! * (1 <= F5 <= 3). ! * F6 = "fudge factor" used in calculation of turbulent ! * diffusivity coefficient. ! * KSTAR = typical gravity wave horizontal wavenumber (1/m) ! * used in calculation of M_SUB_M_TURB. ! * ICUTOFF = 1 to exponentially damp off GWD, heating and diffusion ! * arrays above ALT_CUTOFF; otherwise arrays not modified. ! * ALT_CUTOFF = altitude in meters above which exponential decay applied. ! * SMCO = smoother used to smooth cutoff vertical wavenumbers ! * and total rms winds before calculating drag or heating. ! * (==> a 1:SMCO:1 stencil used; SMCO >= 1.). ! * NSMAX = number of times smoother applied ( >= 1), ! * = 0 means no smoothing performed. ! * IHEATCAL = 1 to calculate heating rates and diffusion coefficient. ! * = 0 means only drag and flux calculated. ! * K_ALPHA = horizontal wavenumber of each azimuth (1/m) which ! * is set here to KSTAR. ! * IERROR = error flag. ! * = 0 no errors. ! * = 10 ==> NAZ > NAZMTH ! * = 20 ==> invalid number of azimuths (NAZ must be 4 or 8). ! * = 30 ==> invalid slope (SLOPE must be 1., 1.5 or 2.). ! * = 40 ==> invalid smoother (SMCO must be >= 1.) ! Input arguements: ! * NMESSG = output unit number where messages to be printed. ! * NLONS = number of longitudes. ! * NAZMTH = azimuthal array dimension (NAZMTH >= NAZ). INTEGER naz, nlons, nazmth, iheatcal, icutoff INTEGER nmessg, nsmax, ierror REAL kstar(nlons), slope, f1, f2, f3, f5, f6, alt_cutoff, smco REAL k_alpha(nlons, nazmth), coslat(nlons) REAL ksmin, ksmax ! Internal variables. INTEGER i, n ! ----------------------------------------------------------------------- ! Specify constants. naz = 8 slope = 1. f1 = 1.5 f2 = 0.3 f3 = 1.0 f5 = 3.0 f6 = 1.0 ksmin = 1.E-5 ksmax = 1.E-4 DO i = 1, nlons kstar(i) = ksmin / (coslat(i) + (ksmin / ksmax)) END DO icutoff = 1 alt_cutoff = 105.E3 smco = 2.0 ! SMCO = 1.0 nsmax = 5 ! NSMAX = 2 iheatcal = 0 ! Print information to output file. ! WRITE (NMESSG,6000) ! 6000 FORMAT (/' Subroutine HINES_SETUP:') ! WRITE (NMESSG,*) ' SLOPE = ', SLOPE ! WRITE (NMESSG,*) ' NAZ = ', NAZ ! WRITE (NMESSG,*) ' F1,F2,F3 = ', F1, F2, F3 ! WRITE (NMESSG,*) ' F5,F6 = ', F5, F6 ! WRITE (NMESSG,*) ' KSTAR = ', KSTAR ! > ,' COSLAT = ', COSLAT ! IF (ICUTOFF .EQ. 1) THEN ! WRITE (NMESSG,*) ' Drag exponentially damped above ', ! & ALT_CUTOFF/1.E3 ! END IF ! IF (NSMAX.LT.1 ) THEN ! WRITE (NMESSG,*) ' No smoothing of cutoff wavenumbers, etc' ! ELSE ! WRITE (NMESSG,*) ' Cutoff wavenumbers and sig_t smoothed:' ! WRITE (NMESSG,*) ' SMCO =', SMCO ! WRITE (NMESSG,*) ' NSMAX =', NSMAX ! END IF ! Check that things are setup correctly and log error if not ierror = 0 IF (naz>nazmth) ierror = 10 IF (naz/=4 .AND. naz/=8) ierror = 20 IF (slope/=1. .AND. slope/=1.5 .AND. slope/=2.) ierror = 30 IF (smco<1.) ierror = 40 ! Use single value for azimuthal-dependent horizontal wavenumber. DO n = 1, naz DO i = 1, nlons k_alpha(i, n) = kstar(i) END DO END DO RETURN ! ----------------------------------------------------------------------- END SUBROUTINE hines_setup SUBROUTINE hines_print(flux_u, flux_v, drag_u, drag_v, alt, sigma_t, & sigma_alpha, v_alpha, m_alpha, iu_print, iv_print, nmessg, ilprt1, & ilprt2, levprt1, levprt2, naz, nlons, nlevs, nazmth) IMPLICIT NONE ! Print out altitude profiles of various quantities from ! Hines' Doppler spread gravity wave drag parameterization scheme. ! (NOTE: only for NAZ = 4 or 8). ! Aug. 8/95 - C. McLandress ! Input arguements: ! * IU_PRINT = 1 to print out values in east-west direction. ! * IV_PRINT = 1 to print out values in north-south direction. ! * NMESSG = unit number for printed output. ! * ILPRT1 = first longitudinal index to print. ! * ILPRT2 = last longitudinal index to print. ! * LEVPRT1 = first altitude level to print. ! * LEVPRT2 = last altitude level to print. INTEGER naz, ilprt1, ilprt2, levprt1, levprt2 INTEGER nlons, nlevs, nazmth INTEGER iu_print, iv_print, nmessg REAL flux_u(nlons, nlevs), flux_v(nlons, nlevs) REAL drag_u(nlons, nlevs), drag_v(nlons, nlevs) REAL alt(nlons, nlevs), sigma_t(nlons, nlevs) REAL sigma_alpha(nlons, nlevs, nazmth) REAL v_alpha(nlons, nlevs, nazmth), m_alpha(nlons, nlevs, nazmth) ! Internal variables. INTEGER n_east, n_west, n_north, n_south INTEGER i, l ! ----------------------------------------------------------------------- ! Azimuthal indices of cardinal directions. n_east = 1 IF (naz==4) THEN n_west = 3 n_north = 2 n_south = 4 ELSE IF (naz==8) THEN n_west = 5 n_north = 3 n_south = 7 END IF ! Print out values for range of longitudes. DO i = ilprt1, ilprt2 ! Print east-west wind, sigmas, cutoff wavenumbers, flux and drag. IF (iu_print==1) THEN WRITE (nmessg, *) WRITE (nmessg, 6001) i WRITE (nmessg, 6005) 6001 FORMAT ('Hines GW (east-west) at longitude I =', I3) 6005 FORMAT (15X, ' U ', 2X, 'sig_E', 2X, 'sig_T', 3X, 'm_E', 4X, 'm_W', 4X, & 'fluxU', 5X, 'gwdU') DO l = levprt1, levprt2 WRITE (nmessg, 6701) alt(i, l) / 1.E3, v_alpha(i, l, n_east), & sigma_alpha(i, l, n_east), sigma_t(i, l), & m_alpha(i, l, n_east) * 1.E3, m_alpha(i, l, n_west) * 1.E3, & flux_u(i, l) * 1.E5, drag_u(i, l) * 24. * 3600. END DO 6701 FORMAT (' z=', F7.2, 1X, 3F7.1, 2F7.3, F9.4, F9.3) END IF ! Print north-south winds, sigmas, cutoff wavenumbers, flux and drag. IF (iv_print==1) THEN WRITE (nmessg, *) WRITE (nmessg, 6002) i 6002 FORMAT ('Hines GW (north-south) at longitude I =', I3) WRITE (nmessg, 6006) 6006 FORMAT (15X, ' V ', 2X, 'sig_N', 2X, 'sig_T', 3X, 'm_N', 4X, 'm_S', 4X, & 'fluxV', 5X, 'gwdV') DO l = levprt1, levprt2 WRITE (nmessg, 6701) alt(i, l) / 1.E3, v_alpha(i, l, n_north), & sigma_alpha(i, l, n_north), sigma_t(i, l), & m_alpha(i, l, n_north) * 1.E3, m_alpha(i, l, n_south) * 1.E3, & flux_v(i, l) * 1.E5, drag_v(i, l) * 24. * 3600. END DO END IF END DO RETURN ! ----------------------------------------------------------------------- END SUBROUTINE hines_print SUBROUTINE hines_exp(data, data_zmax, alt, alt_exp, iorder, il1, il2, lev1, & lev2, nlons, nlevs) IMPLICIT NONE ! This routine exponentially damps a longitude by altitude array ! of data above a specified altitude. ! Aug. 13/95 - C. McLandress ! Output arguements: ! * DATA = modified data array. ! Input arguements: ! * DATA = original data array. ! * ALT = altitudes. ! * ALT_EXP = altitude above which exponential decay applied. ! * IORDER = 1 means vertical levels are indexed from top down ! * (i.e., highest level indexed 1 and lowest level NLEVS); ! * .NE. 1 highest level is index NLEVS. ! * IL1 = first longitudinal index to use (IL1 >= 1). ! * IL2 = last longitudinal index to use (IL1 <= IL2 <= NLONS). ! * LEV1 = first altitude level to use (LEV1 >=1). ! * LEV2 = last altitude level to use (LEV1 < LEV2 <= NLEVS). ! * NLONS = number of longitudes. ! * NLEVS = number of vertical ! Input work arrays: ! * DATA_ZMAX = data values just above altitude ALT_EXP. INTEGER iorder, il1, il2, lev1, lev2, nlons, nlevs REAL alt_exp REAL data(nlons, nlevs), data_zmax(nlons), alt(nlons, nlevs) ! Internal variables. INTEGER levbot, levtop, lincr, i, l REAL hscale DATA hscale/5.E3/ ! ----------------------------------------------------------------------- ! Index of lowest altitude level (bottom of drag calculation). levbot = lev2 levtop = lev1 lincr = 1 IF (iorder/=1) THEN levbot = lev1 levtop = lev2 lincr = -1 END IF ! Data values at first level above ALT_EXP. DO i = il1, il2 DO l = levtop, levbot, lincr IF (alt(i, l)>=alt_exp) THEN data_zmax(i) = data(i, l) END IF END DO END DO ! Exponentially damp field above ALT_EXP to model top at L=1. DO l = 1, lev2 DO i = il1, il2 IF (alt(i, l)>=alt_exp) THEN data(i, l) = data_zmax(i) * exp((alt_exp - alt(i, l)) / hscale) END IF END DO END DO RETURN ! ----------------------------------------------------------------------- END SUBROUTINE hines_exp SUBROUTINE vert_smooth(data, work, coeff, nsmooth, il1, il2, lev1, lev2, & nlons, nlevs) IMPLICIT NONE ! Smooth a longitude by altitude array in the vertical over a ! specified number of levels using a three point smoother. ! NOTE: input array DATA is modified on output! ! Aug. 3/95 - C. McLandress ! Output arguement: ! * DATA = smoothed array (on output). ! Input arguements: ! * DATA = unsmoothed array of data (on input). ! * WORK = work array of same dimension as DATA. ! * COEFF = smoothing coefficient for a 1:COEFF:1 stencil. ! * (e.g., COEFF = 2 will result in a smoother which ! * weights the level L gridpoint by two and the two ! * adjecent levels (L+1 and L-1) by one). ! * NSMOOTH = number of times to smooth in vertical. ! * (e.g., NSMOOTH=1 means smoothed only once, ! * NSMOOTH=2 means smoothing repeated twice, etc.) ! * IL1 = first longitudinal index to use (IL1 >= 1). ! * IL2 = last longitudinal index to use (IL1 <= IL2 <= NLONS). ! * LEV1 = first altitude level to use (LEV1 >=1). ! * LEV2 = last altitude level to use (LEV1 < LEV2 <= NLEVS). ! * NLONS = number of longitudes. ! * NLEVS = number of vertical levels. ! Subroutine arguements. INTEGER nsmooth, il1, il2, lev1, lev2, nlons, nlevs REAL coeff REAL data(nlons, nlevs), work(nlons, nlevs) ! Internal variables. INTEGER i, l, ns, lev1p, lev2m REAL sum_wts ! ----------------------------------------------------------------------- ! Calculate sum of weights. sum_wts = coeff + 2. lev1p = lev1 + 1 lev2m = lev2 - 1 ! Smooth NSMOOTH times DO ns = 1, nsmooth ! Copy data into work array. DO l = lev1, lev2 DO i = il1, il2 work(i, l) = data(i, l) END DO END DO ! Smooth array WORK in vertical direction and put into DATA. DO l = lev1p, lev2m DO i = il1, il2 data(i, l) = (work(i, l + 1) + coeff * work(i, l) + work(i, l - 1)) / sum_wts END DO END DO END DO RETURN ! ----------------------------------------------------------------------- END SUBROUTINE vert_smooth