LCOV - code coverage report
Current view: | top level - physics/cam - gw_common.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 202 | 230 | 87.8 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 9 | 13 | 69.2 %

              Line data    Source code
           1             : module gw_common
           2             :
           3             : !
           4             : ! This module contains code common to different gravity wave
           5             : ! parameterizations.
           6             : !
           7             : use gw_utils, only: r8
           8             : use coords_1d, only: Coords1D
           9             :
          10             :
          11             : implicit none
          12             : private
          13             : save
          14             :
          15             : ! Public interface.
          16             :
          17             : public :: GWBand
          18             :
          19             : public :: gw_common_init
          20             : public :: gw_prof
          21             : public :: gw_drag_prof
          22             : public :: qbo_hdepth_scaling
          23             : public :: calc_taucd, momentum_flux, momentum_fixer
          24             : public :: energy_change, energy_fixer
          25             : public :: coriolis_speed, adjust_inertial
          26             :
          27             : public :: pver
          28             : public :: west, east, north, south
          29             : public :: pi
          30             : public :: gravit
          31             : public :: rair
          32             :
          33             : ! Number of levels in the atmosphere.
          34             : integer, protected :: pver = 0
          35             :
          36             : ! Whether or not to enforce an upper boundary condition of tau = 0.
          37             : logical :: tau_0_ubc = .false.
          38             :
          39             : ! Index the cardinal directions.
          40             : integer, parameter :: west = 1
          41             : integer, parameter :: east = 2
          42             : integer, parameter :: south = 3
          43             : integer, parameter :: north = 4
          44             :
          45             : ! Scaling factor for generating QBO
          46             : real(r8), protected :: qbo_hdepth_scaling
          47             :
          48             : ! 3.14159...
          49             : real(r8), parameter :: pi = acos(-1._r8)
          50             :
          51             : ! Acceleration due to gravity.
          52             : real(r8), protected :: gravit = huge(1._r8)
          53             :
          54             : ! Gas constant for dry air.
          55             : real(r8), protected :: rair = huge(1._r8)
          56             :
          57             : !
          58             : ! Private variables
          59             : !
          60             :
          61             : ! Interface levels for gravity wave sources.
          62             : integer :: ktop = huge(1)
          63             :
          64             : ! Background diffusivity.
          65             : real(r8), parameter :: dback = 0.05_r8
          66             :
          67             : ! rair/gravit
          68             : real(r8) :: rog = huge(1._r8)
          69             :
          70             : ! Newtonian cooling coefficients.
          71             : real(r8), allocatable :: alpha(:)
          72             :
          73             : ! Inverse Prandtl number.
          74             : real(r8) :: prndl
          75             :
          76             : !
          77             : ! Limits to keep values reasonable.
          78             : !
          79             :
          80             : ! Minimum non-zero stress.
          81             : real(r8), parameter :: taumin = 1.e-10_r8
          82             : ! Maximum wind tendency from stress divergence (before efficiency applied).
          83             : ! 400 m/s/day
          84             : real(r8), parameter :: tndmax = 400._r8 / 86400._r8
          85             : ! Maximum allowed change in u-c (before efficiency applied).
          86             : real(r8), parameter :: umcfac = 0.5_r8
          87             : ! Minimum value of (u-c)**2.
          88             : real(r8), parameter :: ubmc2mn = 0.01_r8
          89             :
          90             : ! Type describing a band of wavelengths into which gravity waves can be
          91             : ! emitted.
          92             : ! Currently this has to have uniform spacing (i.e. adjacent elements of
          93             : ! cref are exactly dc apart).
          94             : type :: GWBand
          95             :    ! Dimension of the spectrum.
          96             :    integer :: ngwv
          97             :    ! Delta between nearest phase speeds [m/s].
          98             :    real(r8) :: dc
          99             :    ! Reference speeds [m/s].
         100             :    real(r8), allocatable :: cref(:)
         101             :    ! Critical Froude number, squared
         102             :    real(r8) :: fcrit2
         103             :    ! Horizontal wave number [1/m].
         104             :    real(r8) :: kwv
         105             :    ! Effective horizontal wave number [1/m] (fcrit2*kwv).
         106             :    real(r8) :: effkwv
         107             : end type GWBand
         108             :
         109             : interface GWBand
         110             :    module procedure new_GWBand
         111             : end interface
         112             :
         113             : contains
         114             :
         115             : !==========================================================================
         116             :
         117             : ! Constructor for a GWBand that calculates derived components.
         118        6144 : function new_GWBand(ngwv, dc, fcrit2, wavelength) result(band)
         119             :   ! Used directly to set the type's components.
         120             :   integer, intent(in) :: ngwv
         121             :   real(r8), intent(in) :: dc
         122             :   real(r8), intent(in) :: fcrit2
         123             :   ! Wavelength in meters.
         124             :   real(r8), intent(in) :: wavelength
         125             :
         126             :   ! Output.
         127             :   type(GWBand) :: band
         128             :
         129             :   ! Wavenumber index.
         130             :   integer :: l
         131             :
         132             :   ! Simple assignments.
         133        6144 :   band%ngwv = ngwv
         134        6144 :   band%dc = dc
         135        6144 :   band%fcrit2 = fcrit2
         136             :
         137             :   ! Uniform phase speed reference grid.
         138       18432 :   allocate(band%cref(-ngwv:ngwv))
         139      325632 :   band%cref = [( dc * l, l = -ngwv, ngwv )]
         140             :
         141             :   ! Wavenumber and effective wavenumber come from the wavelength.
         142        6144 :   band%kwv = 2._r8*pi / wavelength
         143        6144 :   band%effkwv = band%fcrit2 * band%kwv
         144             :
         145        6144 : end function new_GWBand
         146             :
         147             : !==========================================================================
         148             :
         149        1536 : subroutine gw_common_init(pver_in, &
         150        3072 :      tau_0_ubc_in, ktop_in, gravit_in, rair_in, alpha_in, &
         151             :      prndl_in, qbo_hdepth_scaling_in, errstring)
         152             :
         153             :   integer,  intent(in) :: pver_in
         154             :   logical,  intent(in) :: tau_0_ubc_in
         155             :   integer,  intent(in) :: ktop_in
         156             :   real(r8), intent(in) :: gravit_in
         157             :   real(r8), intent(in) :: rair_in
         158             :   real(r8), intent(in) :: alpha_in(:)
         159             :   real(r8), intent(in) :: prndl_in
         160             :   real(r8), intent(in) :: qbo_hdepth_scaling_in
         161             :   ! Report any errors from this routine.
         162             :   character(len=*), intent(out) :: errstring
         163             :
         164             :   integer :: ierr
         165             :
         166        1536 :   errstring = ""
         167             :
         168        1536 :   pver = pver_in
         169        1536 :   tau_0_ubc = tau_0_ubc_in
         170        1536 :   ktop = ktop_in
         171        1536 :   gravit = gravit_in
         172        1536 :   rair = rair_in
         173        4608 :   allocate(alpha(pver+1), stat=ierr, errmsg=errstring)
         174           0 :   if (ierr /= 0) return
         175      147456 :   alpha = alpha_in
         176        1536 :   prndl = prndl_in
         177        1536 :   qbo_hdepth_scaling = qbo_hdepth_scaling_in
         178             :
         179        1536 :   rog = rair/gravit
         180             :
         181        1536 : end subroutine gw_common_init
         182             :
         183             : !==========================================================================
         184             :
         185     1489176 : subroutine gw_prof (ncol, p, cpair, t, rhoi, nm, ni)
         186             :   !-----------------------------------------------------------------------
         187             :   ! Compute profiles of background state quantities for the multiple
         188             :   ! gravity wave drag parameterization.
         189             :   !
         190             :   ! The parameterization is assumed to operate only where water vapor
         191             :   ! concentrations are negligible in determining the density.
         192             :   !-----------------------------------------------------------------------
         193             :   use gw_utils, only: midpoint_interp
         194             :   !------------------------------Arguments--------------------------------
         195             :   ! Column dimension.
         196             :   integer, intent(in) :: ncol
         197             :   ! Pressure coordinates.
         198             :   type(Coords1D), intent(in) :: p
         199             :
         200             :   ! Specific heat of dry air, constant pressure.
         201             :   real(r8), intent(in) :: cpair
         202             :   ! Midpoint temperatures.
         203             :   real(r8), intent(in) :: t(ncol,pver)
         204             :
         205             :   ! Interface density.
         206             :   real(r8), intent(out) :: rhoi(ncol,pver+1)
         207             :   ! Midpoint and interface Brunt-Vaisalla frequencies.
         208             :   real(r8), intent(out) :: nm(ncol,pver), ni(ncol,pver+1)
         209             :
         210             :   !---------------------------Local Storage-------------------------------
         211             :   ! Column and level indices.
         212             :   integer :: i,k
         213             :
         214             :   ! dt/dp
         215             :   real(r8) :: dtdp
         216             :   ! Brunt-Vaisalla frequency squared.
         217             :   real(r8) :: n2
         218             :
         219             :   ! Interface temperature.
         220     2978352 :   real(r8) :: ti(ncol,pver+1)
         221             :
         222             :   ! Minimum value of Brunt-Vaisalla frequency squared.
         223             :   real(r8), parameter :: n2min = 5.e-5_r8
         224             :
         225             :   !------------------------------------------------------------------------
         226             :   ! Determine the interface densities and Brunt-Vaisala frequencies.
         227             :   !------------------------------------------------------------------------
         228             :
         229             :   ! The top interface values are calculated assuming an isothermal
         230             :   ! atmosphere above the top level.
         231     1489176 :   k = 1
         232    24865776 :   do i = 1, ncol
         233    23376600 :      ti(i,k) = t(i,k)
         234    23376600 :      rhoi(i,k) = p%ifc(i,k) / (rair*ti(i,k))
         235    24865776 :      ni(i,k) = sqrt(gravit*gravit / (cpair*ti(i,k)))
         236             :   end do
         237             :
         238             :   ! Interior points use centered differences.
         239     1489176 :   ti(:,2:pver) = midpoint_interp(t)
         240   138493368 :   do k = 2, pver
         241  2289140568 :      do i = 1, ncol
         242  2150647200 :         rhoi(i,k) = p%ifc(i,k) / (rair*ti(i,k))
         243  2150647200 :         dtdp = (t(i,k)-t(i,k-1)) * p%rdst(i,k-1)
         244  2150647200 :         n2 = gravit*gravit/ti(i,k) * (1._r8/cpair - rhoi(i,k)*dtdp)
         245  2287651392 :         ni(i,k) = sqrt(max(n2min, n2))
         246             :      end do
         247             :   end do
         248             :
         249             :   ! Bottom interface uses bottom level temperature, density; next interface
         250             :   ! B-V frequency.
         251     1489176 :   k = pver+1
         252    24865776 :   do i = 1, ncol
         253    23376600 :      ti(i,k) = t(i,k-1)
         254    23376600 :      rhoi(i,k) = p%ifc(i,k) / (rair*ti(i,k))
         255    24865776 :      ni(i,k) = ni(i,k-1)
         256             :   end do
         257             :
         258             :   !------------------------------------------------------------------------
         259             :   ! Determine the midpoint Brunt-Vaisala frequencies.
         260             :   !------------------------------------------------------------------------
         261     1489176 :   nm = midpoint_interp(ni)
         262             :
         263     1489176 : end subroutine gw_prof
         264             :
         265             : !==========================================================================
         266             :
         267    19359288 : subroutine gw_drag_prof(ncol, band, p, src_level, tend_level, dt, &
         268    19359288 :      t, vramp,   &
         269    77437152 :      piln, rhoi,    nm,   ni,  ubm,  ubi,  xv,    yv,   &
         270    96796440 :      effgw,      c, kvtt, q,   dse,  tau,  utgw,  vtgw, &
         271    38718576 :      ttgw, qtgw, egwdffi,   gwut, dttdf, dttke, ro_adjust, &
         272    29783520 :      kwvrdg, satfac_in, lapply_effgw_in, lapply_vdiff, tau_diag )
         273             :
         274             :   !-----------------------------------------------------------------------
         275             :   ! Solve for the drag profile from the multiple gravity wave drag
         276             :   ! parameterization.
         277             :   ! 1. scan up from the wave source to determine the stress profile
         278             :   ! 2. scan down the stress profile to determine the tendencies
         279             :   !     => apply bounds to the tendency
         280             :   !          a. from wkb solution
         281             :   !          b. from computational stability constraints
         282             :   !     => adjust stress on interface below to reflect actual bounded
         283             :   !        tendency
         284             :   !-----------------------------------------------------------------------
         285             :
         286             :   use gw_diffusion, only: gw_ediff, gw_diff_tend
         287             :   use linear_1d_operators, only: TriDiagDecomp
         288             :
         289             :   !------------------------------Arguments--------------------------------
         290             :   ! Column dimension.
         291             :   integer, intent(in) :: ncol
         292             :   ! Wavelengths.
         293             :   type(GWBand), intent(in) :: band
         294             :   ! Pressure coordinates.
         295             :   type(Coords1D), intent(in) :: p
         296             :   ! Level from which gravity waves are propagated upward.
         297             :   integer, intent(in) :: src_level(ncol)
         298             :   ! Lowest level where wind tendencies are calculated.
         299             :   integer, intent(in) :: tend_level(ncol)
         300             :   ! Using tend_level > src_level allows the orographic waves to prescribe
         301             :   ! wave propagation up to a certain level, but then allow wind tendencies
         302             :   ! and adjustments to tau below that level.
         303             :
         304             :   ! Time step.
         305             :   real(r8), intent(in) :: dt
         306             :
         307             :   ! Midpoint and interface temperatures.
         308             :   real(r8), intent(in) :: t(ncol,pver)
         309             :   ! Log of interface pressures.
         310             :   real(r8), intent(in) :: piln(ncol,pver+1)
         311             :   ! Interface densities.
         312             :   real(r8), intent(in) :: rhoi(ncol,pver+1)
         313             :   ! Midpoint and interface Brunt-Vaisalla frequencies.
         314             :   real(r8), intent(in) :: nm(ncol,pver), ni(ncol,pver+1)
         315             :   ! Projection of wind at midpoints and interfaces.
         316             :   real(r8), intent(in) :: ubm(ncol,pver), ubi(ncol,pver+1)
         317             :   ! Unit vectors of source wind (zonal and meridional components).
         318             :   real(r8), intent(in) :: xv(ncol), yv(ncol)
         319             :   ! Tendency efficiency.
         320             :   real(r8), intent(in) :: effgw(ncol)
         321             :   ! Wave phase speeds for each column.
         322             :   real(r8), intent(in) :: c(ncol,-band%ngwv:band%ngwv)
         323             :   ! Molecular thermal diffusivity.
         324             :   real(r8), intent(in) :: kvtt(ncol,pver+1)
         325             :   ! Constituent array.
         326             :   real(r8), intent(in) :: q(:,:,:)
         327             :   ! Dry static energy.
         328             :   real(r8), intent(in) :: dse(ncol,pver)
         329             :   ! Coefficient to ramp down diffusion coeff.
         330             :   real(r8), pointer, intent(in) :: vramp(:)
         331             :
         332             :   ! Wave Reynolds stress.
         333             :   real(r8), intent(inout) :: tau(ncol,-band%ngwv:band%ngwv,pver+1)
         334             :   ! Zonal/meridional wind tendencies.
         335             :   real(r8), intent(out) :: utgw(ncol,pver), vtgw(ncol,pver)
         336             :   ! Gravity wave heating tendency.
         337             :   real(r8), intent(out) :: ttgw(ncol,pver)
         338             :   ! Gravity wave constituent tendency.
         339             :   real(r8), intent(out) :: qtgw(:,:,:)
         340             :
         341             :   ! Effective gravity wave diffusivity at interfaces.
         342             :   real(r8), intent(out) :: egwdffi(ncol,pver+1)
         343             :
         344             :   ! Gravity wave wind tendency for each wave.
         345             :   real(r8), intent(out) :: gwut(ncol,pver,-band%ngwv:band%ngwv)
         346             :
         347             :   ! Temperature tendencies from diffusion and kinetic energy.
         348             :   real(r8), intent(out) :: dttdf(ncol,pver)
         349             :   real(r8), intent(out) :: dttke(ncol,pver)
         350             :
         351             :   ! Adjustment parameter for IGWs.
         352             :   real(r8), intent(in), optional :: &
         353             :        ro_adjust(ncol,-band%ngwv:band%ngwv,pver+1)
         354             :
         355             :   ! Diagnosed horizontal wavenumber for ridges.
         356             :   real(r8), intent(in), optional :: &
         357             :        kwvrdg(ncol)
         358             :
         359             :   ! Factor for saturation calculation. Here backwards
         360             :   ! compatibility. I believe it should be 1.0 (jtb).
         361             :   ! Looks like it has been 2.0 for a while in CAM.
         362             :   real(r8), intent(in), optional :: &
         363             :        satfac_in
         364             :
         365             :   logical, intent(in), optional :: lapply_effgw_in, lapply_vdiff
         366             :   ! Provisional Wave Reynolds stress.
         367             :   real(r8), intent(out), optional :: tau_diag(ncol,pver+1)
         368             :
         369             :   !---------------------------Local storage-------------------------------
         370             :
         371             :   ! Level, wavenumber, constituent and column loop indices.
         372             :   integer :: k, l, m, i
         373             :
         374             :   ! Lowest tendency and source levels.
         375             :   integer :: kbot_tend, kbot_src
         376             :
         377             :   ! "Total" and saturation diffusivity.
         378    38718576 :   real(r8) :: d(ncol)
         379             :   ! Imaginary part of vertical wavenumber.
         380    38718576 :   real(r8) :: mi(ncol)
         381             :   ! Stress after damping.
         382    38718576 :   real(r8) :: taudmp(ncol)
         383             :   ! Saturation stress.
         384    38718576 :   real(r8) :: tausat(ncol)
         385             :   ! (ub-c) and (ub-c)**2
         386    38718576 :   real(r8) :: ubmc(ncol), ubmc2(ncol)
         387             :   ! Temporary ubar tendencies (overall, and at wave l).
         388    38718576 :   real(r8) :: ubt(ncol,pver), ubtl(ncol)
         389    38718576 :   real(r8) :: wrk(ncol)
         390             :   ! Ratio used for ubt tndmax limiting.
         391    38718576 :   real(r8) :: ubt_lim_ratio(ncol)
         392             :
         393             :   ! saturation factor. Defaults to 2.0
         394             :   ! unless overidden by satfac_in
         395             :   real(r8) :: satfac
         396             :
         397             :   logical :: lapply_effgw,do_vertical_diffusion
         398             :
         399             :   ! LU decomposition.
         400    19359288 :   type(TriDiagDecomp) :: decomp
         401             :
         402             :   !------------------------------------------------------------------------
         403             :
         404    19359288 :   if (present(satfac_in)) then
         405    14891760 :      satfac = satfac_in
         406             :   else
         407             :      satfac = 2._r8
         408             :   endif
         409             :
         410             :   ! Default behavior is to apply vertical diffusion.
         411             :   ! The user has the option to turn off vert diffusion
         412    19359288 :   do_vertical_diffusion = .true.
         413    19359288 :   if (present(lapply_vdiff)) then
         414    14891760 :      do_vertical_diffusion = lapply_vdiff
         415             :   endif
         416             :
         417             :   ! Default behavior is to apply effgw and
         418             :   ! tendency limiters as designed by Sean
         419             :   ! Santos (lapply_effgw=.TRUE.). However,
         420             :   ! WACCM non-oro GW need to be retuned before
         421             :   ! this can done to them. --jtb 03/02/16
         422    19359288 :   if (present(lapply_effgw_in)) then
         423     4467528 :       lapply_effgw = lapply_effgw_in
         424             :   else
         425             :       lapply_effgw = .TRUE.
         426             :   endif
         427             :
         428             :
         429             :   ! Lowest levels that loops need to iterate over.
         430   323255088 :   kbot_tend = maxval(tend_level)
         431   323255088 :   kbot_src = maxval(src_level)
         432             :
         433             :   ! Initialize gravity wave drag tendencies to zero.
         434             :
         435 30082082472 :   utgw = 0._r8
         436 30082082472 :   vtgw = 0._r8
         437             :
         438 >32629*10^7 :   gwut = 0._r8
         439             :
         440 30082082472 :   dttke = 0._r8
         441 30082082472 :   ttgw = 0._r8
         442             :
         443 30082082472 :   dttdf = 0._r8
         444 >12333*10^8 :   qtgw = 0._r8
         445             :
         446             :   ! Workaround floating point exception issues on Intel by initializing
         447             :   ! everything that's first set in a where block.
         448   323255088 :   mi = 0._r8
         449   323255088 :   taudmp = 0._r8
         450   323255088 :   tausat = 0._r8
         451   323255088 :   ubmc = 0._r8
         452   323255088 :   ubmc2 = 0._r8
         453   323255088 :   wrk = 0._r8
         454             :
         455             :   !------------------------------------------------------------------------
         456             :   ! Compute the stress profiles and diffusivities
         457             :   !------------------------------------------------------------------------
         458             :
         459             :   ! Loop from bottom to top to get stress profiles.
         460             :   ! do k = kbot_src-1, ktop, -1 !++jtb I think this is right
         461   634212796 :   do k = kbot_src, ktop, -1  !++ but this is in model now
         462             :
         463             :      ! Determine the diffusivity for each column.
         464             :
         465 10270322476 :      d = dback + kvtt(:,k)
         466             :
         467 12861896704 :      do l = -band%ngwv, band%ngwv
         468             :
         469             :         ! Determine the absolute value of the saturation stress.
         470             :         ! Define critical levels where the sign of (u-c) changes between
         471             :         ! interfaces.
         472 >20424*10^7 :         ubmc = ubi(:,k) - c(:,l)
         473             :
         474 >20424*10^7 :         tausat = 0.0_r8
         475             :
         476 12227683908 :         if (present(kwvrdg)) then
         477 15601583979 :            where (src_level >= k)
         478             :               ! Test to see if u-c has the same sign here as the level below.
         479   311290601 :               where (ubmc > 0.0_r8 .eqv. ubi(:,k+1) > c(:,l))
         480   311290601 :                  tausat = abs(  kwvrdg  * rhoi(:,k) * ubmc**3 / &
         481             :                     (satfac*ni(:,k)))
         482             :               end where
         483             :            end where
         484             :         else
         485 >58520*10^7 :            where (src_level >= k)
         486             :               ! Test to see if u-c has the same sign here as the level below.
         487 11916393307 :               where (ubmc > 0.0_r8 .eqv. ubi(:,k+1) > c(:,l))
         488             :                  tausat = abs(band%effkwv * rhoi(:,k) * ubmc**3 / &
         489             :                     (satfac*ni(:,k)))
         490             :               end where
         491             :            end where
         492             :         end if
         493             :
         494 12227683908 :         if (present(ro_adjust)) then
         495           0 :            where (src_level >= k)
         496           0 :               tausat = tausat * sqrt(ro_adjust(:,l,k))
         497             :            end where
         498             :         end if
         499             :
         500 12842537416 :         if (present(kwvrdg)) then
         501 49826245723 :            where (src_level >= k)
         502             :               ! Compute stress for each wave. The stress at this level is the
         503             :               ! min of the saturation stress and the stress at the level below
         504             :               ! reduced by damping. The sign of the stress must be the same as
         505             :               ! at the level below.
         506             :
         507             :               ubmc2 = max(ubmc**2, ubmc2mn)
         508   311290601 :               mi = ni(:,k) / (2._r8 *   kwvrdg * ubmc2) * &  ! Is this 2._r8 related to satfac?
         509   311290601 :                  (alpha(k) + ni(:,k)**2/ubmc2 * d)
         510   311290601 :               wrk = -2._r8*mi*rog*t(:,k)*(piln(:,k+1) - piln(:,k))
         511             :
         512             :               taudmp = tau(:,l,k+1)
         513             :
         514             :               ! For some reason, PGI 14.1 loses bit-for-bit reproducibility if
         515             :               ! we limit tau, so instead limit the arrays used to set it.
         516             :               where (tausat <= taumin) tausat = 0._r8
         517             :               where (taudmp <= taumin) taudmp = 0._r8
         518             :
         519             :               tau(:,l,k) = min(taudmp, tausat)
         520             :            end where
         521             :
         522             :         else
         523             :
         524 >19069*10^8 :            where (src_level >= k)
         525             :
         526             :               ! Compute stress for each wave. The stress at this level is the
         527             :               ! min of the saturation stress and the stress at the level below
         528             :               ! reduced by damping. The sign of the stress must be the same as
         529             :               ! at the level below.
         530             :
         531             :               ubmc2 = max(ubmc**2, ubmc2mn)
         532             :               mi = ni(:,k) / (2._r8 * band%kwv * ubmc2) * &
         533 11916393307 :                  (alpha(k) + ni(:,k)**2/ubmc2 * d)
         534 11916393307 :               wrk = -2._r8*mi*rog*t(:,k)*(piln(:,k+1) - piln(:,k))
         535             :
         536             :               taudmp = tau(:,l,k+1) * exp(wrk)
         537             :
         538             :               ! For some reason, PGI 14.1 loses bit-for-bit reproducibility if
         539             :               ! we limit tau, so instead limit the arrays used to set it.
         540             :               where (tausat <= taumin) tausat = 0._r8
         541             :               where (taudmp <= taumin) taudmp = 0._r8
         542             :
         543             :               tau(:,l,k) = min(taudmp, tausat)
         544             :            end where
         545             :         endif
         546             :
         547             :      end do
         548             :   end do
         549             :
         550             :   ! Force tau at the top of the model to zero, if requested.
         551  3525433704 :   if (tau_0_ubc) tau(:,:,ktop) = 0._r8
         552             :
         553             :   ! Write out pre-adjustment tau profile for diagnostc purposes.
         554             :   ! Current implementation only makes sense for orographic waves.
         555             :   ! Fix later.
         556    19359288 :   if (PRESENT(tau_diag)) then
         557 23388721200 :      tau_diag(:,:) = tau(:,0,:)
         558             :   end if
         559             :
         560             :   ! Apply efficiency to completed stress profile.
         561    19359288 :   if (lapply_effgw) then
         562  1414717200 :      do k = ktop, kbot_tend+1
         563  2814542640 :         do l = -band%ngwv, band%ngwv
         564 24773654880 :            where (k-1 <= tend_level)
         565  1399825440 :               tau(:,l,k) = tau(:,l,k) * effgw
         566             :            end where
         567             :         end do
         568             :      end do
         569             :   end if
         570             :
         571             :   !------------------------------------------------------------------------
         572             :   ! Compute the tendencies from the stress divergence.
         573             :   !------------------------------------------------------------------------
         574             :
         575             :   ! Loop over levels from top to bottom
         576  1707855875 :   do k = ktop, kbot_tend
         577             :
         578             :      ! Accumulate the mean wind tendency over wavenumber.
         579 28194966163 :      ubt(:,k) = 0.0_r8
         580             :
         581 14989823574 :      do l = -band%ngwv, band%ngwv    ! loop over wave
         582             :
         583             :         ! Determine the wind tendency, including excess stress carried down
         584             :         ! from above.
         585 >22216*10^7 :         ubtl = gravit * (tau(:,l,k+1)-tau(:,l,k)) * p%rdel(:,k)
         586             :
         587             :         ! Apply first tendency limit to maintain numerical stability.
         588             :         ! Enforce du/dt < |c-u|/dt  so u-c cannot change sign
         589             :         !    (u^n+1 = u^n + du/dt * dt)
         590             :         ! The limiter is somewhat stricter, so that we don't come anywhere
         591             :         ! near reversing c-u.
         592 >22216*10^7 :         ubtl = min(ubtl, umcfac * abs(c(:,l)-ubm(:,k)) / dt)
         593             :
         594 >20042*10^7 :         if (.not. lapply_effgw) ubtl = min(ubtl, tndmax)
         595             :
         596 >85044*10^7 :         where (k <= tend_level)
         597             :
         598             :            ! Save tendency for each wave (for later computation of kzz).
         599             :            ! sign function returns magnitude of ubtl with sign of c-ubm
         600             :            ! Renders ubt/ubm check for mountain waves unecessary
         601             :            gwut(:,k,l) = sign(ubtl, c(:,l)-ubm(:,k))
         602             :            ubt(:,k) = ubt(:,k) + gwut(:,k,l)
         603             :
         604             :         end where
         605             :
         606             :      end do
         607             :
         608  1688496587 :      if (lapply_effgw) then
         609             :         ! Apply second tendency limit to maintain numerical stability.
         610             :         ! Enforce du/dt < tndmax so that ridicuously large tendencies are not
         611             :         ! permitted.
         612             :         ! This can only happen above tend_level, so don't bother checking the
         613             :         ! level explicitly.
         614 >11008*10^7 :         where (abs(ubt(:,k)) > tndmax)
         615             :            ubt_lim_ratio = tndmax/abs(ubt(:,k))
         616             :            ubt(:,k) = ubt_lim_ratio * ubt(:,k)
         617             :         elsewhere
         618             :            ubt_lim_ratio = 1._r8
         619             :         end where
         620             :      else
         621  5069794483 :         ubt_lim_ratio = 1._r8
         622             :      end if
         623             :
         624 14989823574 :      do l = -band%ngwv, band%ngwv
         625 >22216*10^7 :         gwut(:,k,l) = ubt_lim_ratio*gwut(:,k,l)
         626             :         ! Redetermine the effective stress on the interface below from the
         627             :         ! wind tendency. If the wind tendency was limited above, then the
         628             :         ! new stress will be smaller than the old stress, causing stress
         629             :         ! divergence in the next layer down. This smoothes large stress
         630             :         ! divergences downward while conserving total stress.
         631             :
         632             :         ! Protection on SMALL gwut to prevent floating point
         633             :         ! issues.
         634             :         !--------------------------------------------------
         635 >22216*10^7 :         where( abs(gwut(:,k,l)) < 1.e-15_r8 )
         636             :            gwut(:,k,l) = 0._r8
         637             :         endwhere
         638             :
         639 >22385*10^7 :         where (k <= tend_level)
         640 26602653974 :            tau(:,l,k+1) = tau(:,l,k) + &
         641 26602653974 :                 abs(gwut(:,k,l)) * p%del(:,k) / gravit
         642             :         end where
         643             :      end do
         644             :
         645             :      ! Project the mean wind tendency onto the components.
         646 81207905315 :      where (k <= tend_level)
         647             :         utgw(:,k) = ubt(:,k) * xv
         648             :         vtgw(:,k) = ubt(:,k) * yv
         649             :      end where
         650             :
         651  1707855875 :      if (associated(vramp)) then
         652           0 :         utgw(:,k) = utgw(:,k) * vramp(k)
         653           0 :         vtgw(:,k) = vtgw(:,k) * vramp(k)
         654             :      endif
         655             :
         656             :      ! End of level loop.
         657             :   end do
         658             :
         659             :
         660             :   ! Block to undo Sean Santos mods to effgw and limiters.
         661             :   ! Here because non-oro GW in WACCM need extensive re-tuning
         662             :   ! before Sean's mods can be adopted. --jtb 03/02/16
         663             :   !==========================================
         664    19359288 :   if (.not.(lapply_effgw)) then
         665   312497963 :      do k = ktop, kbot_tend+1
         666 12423973326 :         do l = -band%ngwv, band%ngwv
         667 >20260*10^7 :            where (k-1 <= tend_level)
         668 12111475363 :               tau(:,l,k) = tau(:,l,k) * effgw
         669             :            end where
         670             :         end do
         671             :      end do
         672   308030435 :      do k = ktop, kbot_tend
         673 12219956214 :         do l = -band%ngwv, band%ngwv
         674 >19934*10^7 :            gwut(:,k,l) = gwut(:,k,l) * effgw
         675             :         end do
         676  5069794483 :         utgw(:,k) = utgw(:,k) * effgw
         677  5074262011 :         vtgw(:,k) = vtgw(:,k) * effgw
         678             :      end do
         679             :   end if
         680             :   !===========================================
         681             :
         682    19359288 :   if (do_vertical_diffusion) then
         683             :
         684             :      ! Calculate effective diffusivity and LU decomposition for the
         685             :      ! vertical diffusion solver.
         686             :      call gw_ediff (ncol, pver, band%ngwv, kbot_tend, ktop, tend_level, &
         687             :           gwut, ubm, nm, rhoi, dt, prndl, gravit, p, c, vramp, &
         688    19359288 :           egwdffi, decomp, ro_adjust=ro_adjust)
         689             :
         690             :      ! Calculate tendency on each constituent.
         691   813090096 :      do m = 1, size(q,3)
         692             :
         693           0 :         call gw_diff_tend(ncol, pver, kbot_tend, ktop, q(:,:,m), &
         694   813090096 :              dt, decomp, qtgw(:,:,m))
         695             :
         696             :      enddo
         697             :
         698             :      ! Calculate tendency from diffusing dry static energy (dttdf).
         699    19359288 :      call gw_diff_tend(ncol, pver, kbot_tend, ktop, dse, dt, decomp, dttdf)
         700             :
         701             :   endif
         702             :
         703             :   ! Evaluate second temperature tendency term: Conversion of kinetic
         704             :   ! energy into thermal.
         705   229333104 :   do l = -band%ngwv, band%ngwv
         706 13530660091 :      do k = ktop, kbot_tend
         707 >22237*10^7 :         dttke(:,k) = dttke(:,k) - (ubm(:,k) - c(:,l)) * gwut(:,k,l)
         708             :      end do
         709             :   end do
         710             :
         711 30082082472 :   ttgw = dttke + dttdf
         712             :
         713    19359288 :   if (associated(vramp)) then
         714           0 :      do k = ktop, kbot_tend
         715           0 :         ttgw(:,k) = ttgw(:,k) * vramp(k)
         716             :      enddo
         717             :   endif
         718             :
         719             :   ! Deallocate decomp.
         720    19359288 :   call decomp%finalize()
         721             :
         722    53610336 : end subroutine gw_drag_prof
         723             :
         724             : !==========================================================================
         725             :
         726             : ! Calculate Reynolds stress for waves propagating in each cardinal
         727             : ! direction.
         728             :
         729     4467528 : function calc_taucd(ncol, ngwv, tend_level, tau, c, xv, yv, ubi) &
         730    22337640 :      result(taucd)
         731             :
         732             :   ! Column and gravity wave wavenumber dimensions.
         733             :   integer, intent(in) :: ncol, ngwv
         734             :   ! Lowest level where wind tendencies are calculated.
         735             :   integer, intent(in) :: tend_level(:)
         736             :   ! Wave Reynolds stress.
         737             :   real(r8), intent(in) :: tau(:,-ngwv:,:)
         738             :   ! Wave phase speeds for each column.
         739             :   real(r8), intent(in) :: c(:,-ngwv:)
         740             :   ! Unit vectors of source wind (zonal and meridional components).
         741             :   real(r8), intent(in) :: xv(:), yv(:)
         742             :   ! Projection of wind at interfaces.
         743             :   real(r8), intent(in) :: ubi(:,:)
         744             :
         745             :   real(r8) :: taucd(ncol,pver+1,4)
         746             :
         747             :   ! Indices.
         748             :   integer :: i, k, l
         749             :
         750             :   ! ubi at tend_level.
         751     8935056 :   real(r8) :: ubi_tend(ncol)
         752             :
         753             :   ! Signed wave Reynolds stress.
         754     8935056 :   real(r8) :: tausg(ncol)
         755             :
         756             :   ! Reynolds stress for waves propagating behind and forward of the wind.
         757     4467528 :   real(r8) :: taub(ncol)
         758     4467528 :   real(r8) :: tauf(ncol)
         759             :
         760 28070932968 :   taucd = 0._r8
         761    74597328 :   tausg = 0._r8
         762             :
         763   144727128 :   ubi_tend = (/ (ubi(i,tend_level(i)+1), i = 1, ncol) /)
         764             :
         765   382627763 :   do k = ktop, maxval(tend_level)+1
         766             :
         767  5144391811 :      taub = 0._r8
         768  5144391811 :      tauf = 0._r8
         769             :
         770 12419505798 :      do l = -ngwv, ngwv
         771 >98757*10^7 :         where (k-1 <= tend_level)
         772             :
         773 12111475363 :            tausg = sign(tau(:,l,k), c(:,l)-ubi(:,k))
         774             :
         775 12111475363 :            where ( c(:,l) < ubi_tend )
         776             :               taub = taub + tausg
         777             :            elsewhere
         778             :               tauf = tauf + tausg
         779             :            end where
         780             :
         781             :         end where
         782             :      end do
         783             :
         784 54128533969 :      where (k-1 <= tend_level)
         785             :         where (xv > 0._r8)
         786   308030435 :            taucd(:,k,east) = tauf * xv
         787             :            taucd(:,k,west) = taub * xv
         788             :         elsewhere
         789             :            taucd(:,k,east) = taub * xv
         790             :            taucd(:,k,west) = tauf * xv
         791             :         end where
         792             :
         793             :         where ( yv > 0._r8)
         794   308030435 :            taucd(:,k,north) = tauf * yv
         795             :            taucd(:,k,south) = taub * yv
         796             :         elsewhere
         797             :            taucd(:,k,north) = taub * yv
         798             :            taucd(:,k,south) = tauf * yv
         799             :         end where
         800             :      end where
         801             :
         802             :   end do
         803             :
         804     4467528 : end function calc_taucd
         805             :
         806             : !==========================================================================
         807             :
         808             : ! Calculate the amount of momentum conveyed from below the gravity wave
         809             : ! region, to the region where gravity waves are calculated.
         810     2978352 : subroutine momentum_flux(tend_level, taucd, um_flux, vm_flux)
         811             :
         812             :   ! Bottom stress level.
         813             :   integer, intent(in) :: tend_level(:)
         814             :   ! Projected stresses.
         815             :   real(r8), intent(in) :: taucd(:,:,:)
         816             :   ! Components of momentum change sourced from the bottom.
         817             :   real(r8), intent(out) :: um_flux(:), vm_flux(:)
         818             :
         819             :   integer :: i
         820             :
         821             :   ! Tendency for U & V below source level.
         822    49731552 :   do i = 1, size(tend_level)
         823    93506400 :      um_flux(i) = taucd(i,tend_level(i)+1, east) + &
         824   140259600 :                   taucd(i,tend_level(i)+1, west)
         825    93506400 :      vm_flux(i) = taucd(i,tend_level(i)+1,north) + &
         826    96484752 :                   taucd(i,tend_level(i)+1,south)
         827             :   end do
         828             :
         829     2978352 : end subroutine momentum_flux
         830             :
         831             : !==========================================================================
         832             :
         833             : ! Subtracts a change in momentum in the gravity wave levels from wind
         834             : ! tendencies in lower levels, ensuring momentum conservation.
         835     2978352 : subroutine momentum_fixer(tend_level, p, um_flux, vm_flux, utgw, vtgw)
         836             :
         837             :   ! Bottom stress level.
         838             :   integer, intent(in) :: tend_level(:)
         839             :   ! Pressure coordinates.
         840             :   type(Coords1D), intent(in) :: p
         841             :   ! Components of momentum change sourced from the bottom.
         842             :   real(r8), intent(in) :: um_flux(:), vm_flux(:)
         843             :   ! Wind tendencies.
         844             :   real(r8), intent(inout) :: utgw(:,:), vtgw(:,:)
         845             :
         846             :   ! Indices.
         847             :   integer :: i, k
         848             :   ! Reciprocal of total mass.
         849     5956704 :   real(r8) :: rdm(size(tend_level))
         850             :   ! Average changes in velocity from momentum change being spread over
         851             :   ! total mass.
         852     5956704 :   real(r8) :: du(size(tend_level)), dv(size(tend_level))
         853             :
         854             :   ! Total mass from ground to source level: rho*dz = dp/gravit
         855    49731552 :   do i = 1, size(tend_level)
         856    49731552 :      rdm(i) = gravit/(p%ifc(i,pver+1)-p%ifc(i,tend_level(i)+1))
         857             :   end do
         858             :
         859             :   ! Average velocity changes.
         860    52709904 :   du = -um_flux*rdm
         861    52709904 :   dv = -vm_flux*rdm
         862             :
         863   169533657 :   do k = minval(tend_level)+1, pver
         864 10245051657 :      where (k > tend_level)
         865   239604210 :         utgw(:,k) = utgw(:,k) + du
         866   239604210 :         vtgw(:,k) = vtgw(:,k) + dv
         867             :      end where
         868             :   end do
         869             :
         870     2978352 : end subroutine momentum_fixer
         871             :
         872             : !==========================================================================
         873             :
         874             : ! Calculate the change in total energy from tendencies up to this point.
         875     4467528 : subroutine energy_change(dt, p, u, v, dudt, dvdt, dsdt, de)
         876             :
         877             :   ! Time step.
         878             :   real(r8), intent(in) :: dt
         879             :   ! Pressure coordinates.
         880             :   type(Coords1D), intent(in) :: p
         881             :   ! Winds at start of time step.
         882             :   real(r8), intent(in) :: u(:,:), v(:,:)
         883             :   ! Wind tendencies.
         884             :   real(r8), intent(in) :: dudt(:,:), dvdt(:,:)
         885             :   ! Heating tendency.
         886             :   real(r8), intent(in) :: dsdt(:,:)
         887             :   ! Change in energy.
         888             :   real(r8), intent(out) :: de(:)
         889             :
         890             :   ! Level index.
         891             :   integer :: k
         892             :
         893             :   ! Net gain/loss of total energy in the column.
         894    74597328 :   de = 0.0_r8
         895   419947632 :   do k = 1, pver
         896   415480104 :      de = de + p%del(:,k)/gravit * (dsdt(:,k) + &
         897           0 :           dudt(:,k)*(u(:,k)+dudt(:,k)*0.5_r8*dt) + &
         898  7357499136 :           dvdt(:,k)*(v(:,k)+dvdt(:,k)*0.5_r8*dt) )
         899             :   end do
         900             :
         901     4467528 : end subroutine energy_change
         902             :
         903             : !==========================================================================
         904             :
         905             : ! Subtract change in energy from the heating tendency in the levels below
         906             : ! the gravity wave region.
         907     2978352 : subroutine energy_fixer(tend_level, p, de, ttgw)
         908             :
         909             :   ! Bottom stress level.
         910             :   integer, intent(in) :: tend_level(:)
         911             :   ! Pressure coordinates.
         912             :   type(Coords1D), intent(in) :: p
         913             :   ! Change in energy.
         914             :   real(r8), intent(in) :: de(:)
         915             :   ! Heating tendency.
         916             :   real(r8), intent(inout) :: ttgw(:,:)
         917             :
         918             :   ! Column/level indices.
         919             :   integer :: i, k
         920             :   ! Energy change to apply divided by all the mass it is spread across.
         921     5956704 :   real(r8) :: de_dm(size(tend_level))
         922             :
         923    49731552 :   do i = 1, size(tend_level)
         924    49731552 :      de_dm(i) = -de(i)*gravit/(p%ifc(i,pver+1)-p%ifc(i,tend_level(i)+1))
         925             :   end do
         926             :
         927             :   ! Subtract net gain/loss of total energy below tend_level.
         928   169533657 :   do k = minval(tend_level)+1, pver
         929  2123274276 :      where (k > tend_level)
         930   239604210 :         ttgw(:,k) = ttgw(:,k) + de_dm
         931             :      end where
         932             :   end do
         933             :
         934     2978352 : end subroutine energy_fixer
         935             :
         936             : !==========================================================================
         937             :
         938             : ! Calculates absolute value of the local Coriolis frequency divided by the
         939             : ! spatial frequency kwv, which gives a characteristic speed in m/s.
         940           0 : function coriolis_speed(band, lat)
         941             :   ! Inertial gravity wave lengths.
         942             :   type(GWBand), intent(in) :: band
         943             :   ! Latitude in radians.
         944             :   real(r8), intent(in) :: lat(:)
         945             :
         946             :   real(r8) :: coriolis_speed(size(lat))
         947             :
         948             :   ! 24*3600 = 86400 seconds in a day.
         949             :   real(r8), parameter :: omega_earth = 2._r8*pi/86400._r8
         950             :
         951           0 :   coriolis_speed = abs(sin(lat) * 2._r8 * omega_earth / band%kwv)
         952             :
         953             : end function coriolis_speed
         954             :
         955             : !==========================================================================
         956             :
         957           0 : subroutine adjust_inertial(band, tend_level, &
         958           0 :      u_coriolis, c, ubi, tau, ro_adjust)
         959             :   ! Inertial gravity wave lengths.
         960             :   type(GWBand), intent(in) :: band
         961             :   ! Levels above which tau is calculated.
         962             :   integer, intent(in) :: tend_level(:)
         963             :   ! Absolute value of the Coriolis frequency for each column,
         964             :   ! divided by kwv [m/s].
         965             :   real(r8), intent(in) :: u_coriolis(:)
         966             :   ! Wave propagation speed.
         967             :   real(r8), intent(in) :: c(:,-band%ngwv:)
         968             :   ! Wind speed in the direction of wave propagation.
         969             :   real(r8), intent(in) :: ubi(:,:)
         970             :
         971             :   ! Tau will be adjusted by blocking wave propagation through cells where
         972             :   ! the Coriolis effect prevents it.
         973             :   real(r8), intent(inout) :: tau(:,-band%ngwv:,:)
         974             :   ! Dimensionless Coriolis term used to reduce gravity wave strength.
         975             :   ! Equal to max(0, 1 - (1/ro)^2), where ro is the Rossby number of the
         976             :   ! wind with respect to inertial waves.
         977             :   real(r8), intent(out) :: ro_adjust(:,-band%ngwv:,:)
         978             :
         979             :   ! Column/level/wavenumber indices.
         980             :   integer :: i, k, l
         981             :
         982             :   ! For each column and wavenumber, are we clear of levels that block
         983             :   ! upward propagation?
         984           0 :   logical :: unblocked_mask(size(tend_level),-band%ngwv:band%ngwv)
         985             :
         986           0 :   unblocked_mask = .true.
         987           0 :   ro_adjust = 0._r8
         988             :
         989             :   ! Iterate from the bottom up, through every interface level where tau is
         990             :   ! set.
         991           0 :   do k = maxval(tend_level)+1, ktop, -1
         992           0 :      do l = -band%ngwv, band%ngwv
         993           0 :         do i = 1, size(tend_level)
         994             :            ! Only operate on valid levels for this column.
         995           0 :            if (k <= tend_level(i) + 1) then
         996             :               ! Block waves if Coriolis is too strong.
         997             :               ! By setting the mask in this way, we avoid division by zero.
         998           0 :               unblocked_mask(i,l) = unblocked_mask(i,l) .and. &
         999           0 :                    (abs(ubi(i,k) - c(i,l)) > u_coriolis(i))
        1000           0 :               if (unblocked_mask(i,l)) then
        1001           0 :                  ro_adjust(i,l,k) = &
        1002           0 :                       1._r8 - (u_coriolis(i)/(ubi(i,k)-c(i,l)))**2
        1003             :               else
        1004           0 :                  tau(i,l,k) = 0._r8
        1005             :               end if
        1006             :            end if
        1007             :         end do
        1008             :      end do
        1009             :   end do
        1010             :
        1011           0 : end subroutine adjust_inertial
        1012             :
        1013           0 : end module gw_common

Generated by: LCOV version 1.14