LCOV - code coverage report
Current view: | top level - physics/cam - gw_drag.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 318 | 1126 | 28.2 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 3 | 13 | 23.1 %

              Line data    Source code
           1             : module gw_drag
           2             :
           3             : !--------------------------------------------------------------------------
           4             : ! CAM and WACCM gravity wave parameterizations were merged by Sean Patrick
           5             : ! Santos in Summer 2013, and at the same time, gw_drag was split into
           6             : ! various modules. This is the CAM interface and driver module. The below
           7             : ! notes are for the old CAM and WACCM versions of gw_drag.
           8             : !--------------------------------------------------------------------------
           9             : ! This file came from wa17 and was modified by Fabrizio: 07-02-2004
          10             : ! Standard gw_drag with modification (6) of latitude profile of gw spectrum
          11             : !--------------------------------------------------------------------------
          12             : ! Purpose:
          13             : !
          14             : ! Module to compute the forcing due to parameterized gravity waves. Both an
          15             : ! orographic and an internal source spectrum are considered.
          16             : !
          17             : ! Author: Byron Boville
          18             : !
          19             : !--------------------------------------------------------------------------
          20             :   use shr_kind_mod,   only: r8=>shr_kind_r8, cl=>shr_kind_cl
          21             :   use shr_log_mod,    only: errMsg => shr_log_errMsg
          22             :   use shr_assert_mod, only: shr_assert
          23             :
          24             :   use ppgrid,         only: pcols, pver, begchunk, endchunk
          25             :   use constituents,   only: pcnst
          26             :   use physics_types,  only: physics_state, physics_ptend, physics_ptend_init
          27             :   use spmd_utils,     only: masterproc
          28             :   use cam_history,    only: outfld
          29             :   use cam_logfile,    only: iulog
          30             :   use cam_abortutils, only: endrun
          31             :   use error_messages, only: alloc_err
          32             :
          33             :
          34             :   use ref_pres,       only: do_molec_diff, nbot_molec, press_lim_idx
          35             :   use physconst,      only: cpair
          36             :
          37             :   ! These are the actual switches for different gravity wave sources.
          38             :   use phys_control,   only: use_gw_oro, use_gw_front, use_gw_front_igw, &
          39             :                             use_gw_convect_dp, use_gw_convect_sh,       &
          40             :                             use_simple_phys, use_gw_movmtn_pbl
          41             :
          42             :   use gw_common,      only: GWBand
          43             :   use gw_convect,     only: BeresSourceDesc
          44             :   use gw_movmtn,      only: MovMtnSourceDesc
          45             :   use gw_front,       only: CMSourceDesc
          46             :
          47             : ! Typical module header
          48             :   implicit none
          49             :   private
          50             :   save
          51             :
          52             : !
          53             : ! PUBLIC: interfaces
          54             : !
          55             :   public :: gw_drag_readnl           ! Read namelist
          56             :   public :: gw_init                  ! Initialization
          57             :   public :: gw_tend                  ! interface to actual parameterization
          58             :
          59             : !
          60             : ! PRIVATE: Rest of the data and interfaces are private to this module
          61             : !
          62             :   real(r8), parameter :: unset_r8 = huge(1._r8)
          63             :
          64             :   ! A mid-scale "band" with only stationary waves (l = 0).
          65             :   type(GWBand) :: band_oro
          66             :   ! Medium scale waves.
          67             :   type(GWBand) :: band_mid
          68             :   ! Long scale waves for IGWs.
          69             :   type(GWBand) :: band_long
          70             :   ! Medium scale waves for moving mountain
          71             :   type(GWBand) :: band_movmtn
          72             :
          73             :   ! Top level for gravity waves.
          74             :   integer, parameter :: ktop = 1
          75             :   ! Bottom level for frontal waves.
          76             :   integer :: kbot_front
          77             :
          78             :   ! Factor for SH orographic waves.
          79             :   real(r8) :: gw_oro_south_fac = 1._r8
          80             :
          81             :   ! Frontogenesis function critical threshold.
          82             :   real(r8) :: frontgfc = unset_r8
          83             :
          84             :   ! Tendency efficiencies.
          85             :
          86             :   ! Ridge scheme.
          87             :   logical  :: use_gw_rdg_beta      = .false.
          88             :   integer  :: n_rdg_beta           = 1
          89             :   real(r8) :: effgw_rdg_beta       = unset_r8
          90             :   real(r8) :: effgw_rdg_beta_max   = unset_r8
          91             :   real(r8) :: rdg_beta_cd_llb      = unset_r8  ! Low-level obstacle drag coefficient Ridge scheme.
          92             :   logical  :: trpd_leewv_rdg_beta  = .false.
          93             :
          94             :   logical  :: use_gw_rdg_gamma     = .false.
          95             :   integer  :: n_rdg_gamma          = -1
          96             :   real(r8) :: effgw_rdg_gamma      = unset_r8
          97             :   real(r8) :: effgw_rdg_gamma_max  = unset_r8
          98             :   real(r8) :: rdg_gamma_cd_llb     = unset_r8
          99             :   logical  :: trpd_leewv_rdg_gamma = .false.
         100             :   character(len=cl) :: bnd_rdggm   = 'bnd_rdggm' ! full pathname for meso-Gamma ridge dataset
         101             :
         102             :   ! Orography.
         103             :   real(r8) :: effgw_oro = unset_r8
         104             :   ! C&M scheme.
         105             :   real(r8) :: effgw_cm = unset_r8
         106             :   ! C&M scheme (inertial waves).
         107             :   real(r8) :: effgw_cm_igw = unset_r8
         108             :   ! Beres (deep convection).
         109             :   real(r8) :: effgw_beres_dp = unset_r8
         110             :   ! Beres (shallow convection).
         111             :   real(r8) :: effgw_beres_sh = unset_r8
         112             :
         113             :   ! Horzontal wavelengths [m].
         114             :   real(r8), parameter :: wavelength_mid = 1.e5_r8
         115             :   real(r8), parameter :: wavelength_long = 3.e5_r8
         116             :
         117             :   ! Background stress source strengths.
         118             :   real(r8) :: taubgnd = unset_r8
         119             :   real(r8) :: taubgnd_igw = unset_r8
         120             :
         121             :   ! Whether or not to use a polar taper for frontally generated waves.
         122             :   logical :: gw_polar_taper = .false.
         123             :
         124             :   ! Whether or not to enforce an upper boundary condition of tau = 0.
         125             :   ! (Like many variables, this is only here to hold the value between
         126             :   ! the readnl phase and the init phase of the CAM physics; only gw_common
         127             :   ! should actually use it.)
         128             :   logical :: tau_0_ubc = .false.
         129             :
         130             :   ! Whether or not to limit tau *before* applying any efficiency factors.
         131             :   logical :: gw_limit_tau_without_eff = .false.
         132             :
         133             :   ! Whether or not to apply tendency max
         134             :   logical :: gw_apply_tndmax = .true.
         135             :
         136             :   ! Files to read Beres source spectra from.
         137             :   character(len=cl) :: gw_drag_file = ""
         138             :   character(len=cl) :: gw_drag_file_sh = ""
         139             :   character(len=cl) :: gw_drag_file_mm = ""
         140             :
         141             :   ! Beres settings and table.
         142             :   type(BeresSourceDesc) :: beres_dp_desc
         143             :   type(BeresSourceDesc) :: beres_sh_desc
         144             :
         145             :   ! Moving mountain settings and table.
         146             :   type(MovMtnSourceDesc) :: movmtn_desc
         147             :
         148             :   ! Frontogenesis wave settings.
         149             :   type(CMSourceDesc) :: cm_desc
         150             :   type(CMSourceDesc) :: cm_igw_desc
         151             :
         152             :   ! Indices into pbuf
         153             :   integer :: kvt_idx      = -1
         154             :   integer :: ttend_dp_idx = -1
         155             :   integer :: ttend_sh_idx = -1
         156             :   integer :: frontgf_idx  = -1
         157             :   integer :: frontga_idx  = -1
         158             :   integer :: sgh_idx      = -1
         159             :
         160             :   ! From CLUBB
         161             :   integer :: ttend_clubb_idx  = -1
         162             :   integer :: upwp_clubb_gw_idx   = -1
         163             :   integer :: vpwp_clubb_gw_idx   = -1
         164             :   integer :: thlp2_clubb_gw_idx  = -1
         165             :   integer :: wpthlp_clubb_gw_idx  = -1
         166             :
         167             :   ! anisotropic ridge fields
         168             :   integer, parameter :: prdg = 16
         169             :
         170             :   real(r8), allocatable, dimension(:,:), target :: &
         171             :      rdg_gbxar
         172             :
         173             :      ! Meso Beta
         174             :   real(r8), allocatable, dimension(:,:,:), target :: &
         175             :      rdg_hwdth,  &
         176             :      rdg_clngt,  &
         177             :      rdg_mxdis,  &
         178             :      rdg_anixy,  &
         179             :      rdg_angll
         180             :
         181             :      ! Meso Gamma
         182             :   real(r8), allocatable, dimension(:,:,:), target :: &
         183             :      rdg_hwdthg, &
         184             :      rdg_clngtg, &
         185             :      rdg_mxdisg, &
         186             :      rdg_anixyg, &
         187             :      rdg_angllg
         188             :
         189             :   ! Water constituent indices for budget
         190             :   integer :: ixcldliq = -1
         191             :   integer :: ixcldice = -1
         192             :
         193             :   ! Prefixes for history field names
         194             :   character(len=1), parameter :: cm_pf = " "
         195             :   character(len=1), parameter :: cm_igw_pf = "I"
         196             :   character(len=1), parameter :: beres_dp_pf = "B"
         197             :   character(len=1), parameter :: beres_sh_pf = "S"
         198             :
         199             :   ! namelist
         200             :   logical          :: history_amwg                   ! output the variables used by the AMWG diag package
         201             :   logical  :: gw_lndscl_sgh = .true. ! scale SGH by land frac
         202             :   real(r8) :: gw_prndl = 0.25_r8
         203             :   real(r8) :: gw_qbo_hdepth_scaling = 1._r8 ! heating depth scaling factor
         204             :
         205             :   ! Width of gaussian used to create frontogenesis tau profile [m s-1].
         206             :   real(r8) :: front_gaussian_width = -huge(1._r8)
         207             :
         208             :   real(r8) :: alpha_gw_movmtn
         209             :
         210             :   logical :: gw_top_taper=.false.
         211             :   real(r8), pointer :: vramp(:)=>null()
         212             :
         213             : !==========================================================================
         214             : contains
         215             : !==========================================================================
         216             :
         217        1536 : subroutine gw_drag_readnl(nlfile)
         218             :
         219             :   use namelist_utils,  only: find_group_name
         220             :   use units,           only: getunit, freeunit
         221             :   use spmd_utils,      only: mpicom, mstrid=>masterprocid, mpi_real8, &
         222             :                              mpi_character, mpi_logical, mpi_integer
         223             :   use gw_rdg,          only: gw_rdg_readnl
         224             :
         225             :   ! File containing namelist input.
         226             :   character(len=*), intent(in) :: nlfile
         227             :
         228             :   ! Local variables
         229             :   integer :: unitn, ierr
         230             :   character(len=*), parameter :: sub = 'gw_drag_readnl'
         231             :
         232             :   ! Maximum wave number and width of spectrum bins.
         233             :   integer :: pgwv = -1
         234             :   real(r8) :: gw_dc = unset_r8
         235             :   integer :: pgwv_long = -1
         236             :   real(r8) :: gw_dc_long = unset_r8
         237             :
         238             :   namelist /gw_drag_nl/ pgwv, gw_dc, pgwv_long, gw_dc_long, tau_0_ubc, &
         239             :        effgw_beres_dp, effgw_beres_sh, effgw_cm, effgw_cm_igw, effgw_oro, &
         240             :        frontgfc, gw_drag_file, gw_drag_file_sh, gw_drag_file_mm, taubgnd, &
         241             :        taubgnd_igw, gw_polar_taper, &
         242             :        use_gw_rdg_beta, n_rdg_beta, effgw_rdg_beta, effgw_rdg_beta_max, &
         243             :        rdg_beta_cd_llb, trpd_leewv_rdg_beta, &
         244             :        use_gw_rdg_gamma, n_rdg_gamma, effgw_rdg_gamma, effgw_rdg_gamma_max, &
         245             :        rdg_gamma_cd_llb, trpd_leewv_rdg_gamma, bnd_rdggm, &
         246             :        gw_oro_south_fac, gw_limit_tau_without_eff, &
         247             :        gw_lndscl_sgh, gw_prndl, gw_apply_tndmax, gw_qbo_hdepth_scaling, &
         248             :        gw_top_taper, front_gaussian_width, alpha_gw_movmtn
         249             :   !----------------------------------------------------------------------
         250             :
         251        1536 :   if (use_simple_phys) return
         252             :
         253        1536 :   if (masterproc) then
         254           2 :      unitn = getunit()
         255           2 :      open( unitn, file=trim(nlfile), status='old' )
         256           2 :      call find_group_name(unitn, 'gw_drag_nl', status=ierr)
         257           2 :      if (ierr == 0) then
         258           2 :         read(unitn, gw_drag_nl, iostat=ierr)
         259           2 :         if (ierr /= 0) then
         260           0 :            call endrun(sub // ':: ERROR reading namelist')
         261             :         end if
         262             :      end if
         263           2 :      close(unitn)
         264           2 :      call freeunit(unitn)
         265             :   end if
         266             :
         267        1536 :   call mpi_bcast(pgwv, 1, mpi_integer, mstrid, mpicom, ierr)
         268        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: pgwv")
         269        1536 :   call mpi_bcast(gw_dc, 1, mpi_real8, mstrid, mpicom, ierr)
         270        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_dc")
         271        1536 :   call mpi_bcast(pgwv_long, 1, mpi_integer, mstrid, mpicom, ierr)
         272        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: pgwv_long")
         273        1536 :   call mpi_bcast(gw_dc_long, 1, mpi_real8, mstrid, mpicom, ierr)
         274        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_dc_long")
         275        1536 :   call mpi_bcast(tau_0_ubc, 1, mpi_logical, mstrid, mpicom, ierr)
         276        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: tau_0_ubc")
         277        1536 :   call mpi_bcast(effgw_beres_dp, 1, mpi_real8, mstrid, mpicom, ierr)
         278        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: effgw_beres_dp")
         279        1536 :   call mpi_bcast(effgw_beres_sh, 1, mpi_real8, mstrid, mpicom, ierr)
         280        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: effgw_beres_sh")
         281        1536 :   call mpi_bcast(effgw_cm, 1, mpi_real8, mstrid, mpicom, ierr)
         282        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: effgw_cm")
         283        1536 :   call mpi_bcast(effgw_cm_igw, 1, mpi_real8, mstrid, mpicom, ierr)
         284        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: effgw_cm_igw")
         285        1536 :   call mpi_bcast(effgw_oro, 1, mpi_real8, mstrid, mpicom, ierr)
         286        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: effgw_oro")
         287             :
         288        1536 :   call mpi_bcast(use_gw_rdg_beta, 1, mpi_logical, mstrid, mpicom, ierr)
         289        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: use_gw_rdg_beta")
         290        1536 :   call mpi_bcast(n_rdg_beta, 1, mpi_integer, mstrid, mpicom, ierr)
         291        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: n_rdg_beta")
         292        1536 :   call mpi_bcast(effgw_rdg_beta, 1, mpi_real8, mstrid, mpicom, ierr)
         293        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: effgw_rdg_beta")
         294        1536 :   call mpi_bcast(effgw_rdg_beta_max, 1, mpi_real8, mstrid, mpicom, ierr)
         295        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: effgw_rdg_beta_max")
         296        1536 :   call mpi_bcast(rdg_beta_cd_llb, 1, mpi_real8, mstrid, mpicom, ierr)
         297        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: rdg_beta_cd_llb")
         298        1536 :   call mpi_bcast(trpd_leewv_rdg_beta, 1, mpi_logical, mstrid, mpicom, ierr)
         299        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: trpd_leewv_rdg_beta")
         300             :
         301        1536 :   call mpi_bcast(use_gw_rdg_gamma, 1, mpi_logical, mstrid, mpicom, ierr)
         302        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: use_gw_rdg_gamma")
         303        1536 :   call mpi_bcast(n_rdg_gamma, 1, mpi_integer, mstrid, mpicom, ierr)
         304        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: n_rdg_gamma")
         305        1536 :   call mpi_bcast(effgw_rdg_gamma, 1, mpi_real8, mstrid, mpicom, ierr)
         306        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: effgw_rdg_gamma")
         307        1536 :   call mpi_bcast(effgw_rdg_gamma_max, 1, mpi_real8, mstrid, mpicom, ierr)
         308        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: effgw_rdg_gamma_max")
         309        1536 :   call mpi_bcast(rdg_gamma_cd_llb, 1, mpi_real8, mstrid, mpicom, ierr)
         310        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: rdg_gamma_cd_llb")
         311        1536 :   call mpi_bcast(trpd_leewv_rdg_gamma, 1, mpi_logical, mstrid, mpicom, ierr)
         312        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: trpd_leewv_rdg_gamma")
         313        1536 :   call mpi_bcast(bnd_rdggm, len(bnd_rdggm), mpi_character, mstrid, mpicom, ierr)
         314        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: bnd_rdggm")
         315             :
         316        1536 :   call mpi_bcast(gw_oro_south_fac, 1, mpi_real8, mstrid, mpicom, ierr)
         317        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_oro_south_fac")
         318        1536 :   call mpi_bcast(frontgfc, 1, mpi_real8, mstrid, mpicom, ierr)
         319        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: frontgfc")
         320        1536 :   call mpi_bcast(taubgnd, 1, mpi_real8, mstrid, mpicom, ierr)
         321        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: taubgnd")
         322        1536 :   call mpi_bcast(taubgnd_igw, 1, mpi_real8, mstrid, mpicom, ierr)
         323        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: taubgnd_igw")
         324             :
         325        1536 :   call mpi_bcast(gw_polar_taper, 1, mpi_logical, mstrid, mpicom, ierr)
         326        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_polar_taper")
         327        1536 :   call mpi_bcast(gw_limit_tau_without_eff, 1, mpi_logical, mstrid, mpicom, ierr)
         328        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_limit_tau_without_eff")
         329        1536 :   call mpi_bcast(gw_apply_tndmax, 1, mpi_logical, mstrid, mpicom, ierr)
         330        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_apply_tndmax")
         331             :
         332        1536 :   call mpi_bcast(gw_top_taper, 1, mpi_logical, mstrid, mpicom, ierr)
         333        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_top_taper")
         334             :
         335        1536 :   call mpi_bcast(gw_lndscl_sgh, 1, mpi_logical, mstrid, mpicom, ierr)
         336        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_lndscl_sgh")
         337        1536 :   call mpi_bcast(gw_prndl, 1, mpi_real8, mstrid, mpicom, ierr)
         338        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_prndl")
         339        1536 :   call mpi_bcast(gw_qbo_hdepth_scaling, 1, mpi_real8, mstrid, mpicom, ierr)
         340        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_qbo_hdepth_scaling")
         341             :
         342        1536 :   call mpi_bcast(gw_drag_file, len(gw_drag_file), mpi_character, mstrid, mpicom, ierr)
         343        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_drag_file")
         344        1536 :   call mpi_bcast(gw_drag_file_sh, len(gw_drag_file_sh), mpi_character, mstrid, mpicom, ierr)
         345        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_drag_file_sh")
         346        1536 :   call mpi_bcast(gw_drag_file_mm, len(gw_drag_file_mm), mpi_character, mstrid, mpicom, ierr)
         347        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_drag_file_mm")
         348             :
         349        1536 :   call mpi_bcast(front_gaussian_width, 1, mpi_real8, mstrid, mpicom, ierr)
         350        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: front_gaussian_width")
         351             :
         352        1536 :   call mpi_bcast(alpha_gw_movmtn, 1, mpi_real8, mstrid, mpicom, ierr)
         353        1536 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: alpha_gw_movmtn")
         354             :
         355             :   ! Check if pgwv was set.
         356             :   call shr_assert(pgwv >= 0, &
         357             :        "gw_drag_readnl: pgwv must be set via the namelist and &
         358             :        &non-negative."// &
         359        1536 :        errMsg(__FILE__, __LINE__))
         360             :
         361             :   ! Check if gw_dc was set.
         362             :   call shr_assert(gw_dc /= unset_r8, &
         363             :        "gw_drag_readnl: gw_dc must be set via the namelist."// &
         364        1536 :        errMsg(__FILE__, __LINE__))
         365             :
         366        1536 :   band_oro = GWBand(0, gw_dc, 1.0_r8, wavelength_mid)
         367        1536 :   band_mid = GWBand(pgwv, gw_dc, 1.0_r8, wavelength_mid)
         368        1536 :   band_long = GWBand(pgwv_long, gw_dc_long, 1.0_r8, wavelength_long)
         369        1536 :   band_movmtn = GWBand(0, gw_dc, 1.0_r8, wavelength_mid)
         370             :
         371        1536 :   if (use_gw_rdg_gamma .or. use_gw_rdg_beta) then
         372           0 :      call gw_rdg_readnl(nlfile)
         373             :   end if
         374             :
         375             : end subroutine gw_drag_readnl
         376             :
         377             : !==========================================================================
         378             :
         379        1536 : subroutine gw_init()
         380             :   !-----------------------------------------------------------------------
         381             :   ! Time independent initialization for multiple gravity wave
         382             :   ! parameterization.
         383             :   !-----------------------------------------------------------------------
         384             :
         385             :   use cam_history,      only: addfld, add_default, horiz_only
         386             :   use cam_history,      only: register_vector_field
         387             :   use interpolate_data, only: lininterp
         388             :   use phys_control,     only: phys_getopts
         389             :   use physics_buffer,   only: pbuf_get_index
         390             :   use constituents,     only: cnst_get_ind
         391             :
         392             :   use cam_initfiles,    only: topo_file_get_id
         393             :
         394             :   ! temporary for restart with ridge scheme
         395             :   use cam_initfiles,    only: bnd_topo
         396             :
         397             :   use cam_pio_utils,    only: cam_pio_openfile
         398             :   use cam_grid_support, only: cam_grid_check, cam_grid_id
         399             :   use cam_grid_support, only: cam_grid_get_dim_names
         400             :   use pio,              only: file_desc_t, pio_nowrite, pio_closefile
         401             :   use ncdio_atm,        only: infld
         402             :   use ioFileMod,        only: getfil
         403             :
         404             :   use ref_pres,   only: pref_edge, pref_mid
         405             :   use physconst,  only: gravit, rair, rearth
         406             :
         407             :   use gw_common,  only: gw_common_init
         408             :   use gw_front,   only: gaussian_cm_desc
         409             :
         410             :   !---------------------------Local storage-------------------------------
         411             :
         412             :   integer          :: i, l, k
         413             :   character(len=1) :: cn
         414             :
         415             :   ! Index for levels at specific pressures.
         416             :   integer :: kfront
         417             :
         418             :   ! output tendencies and state variables for CAM4 temperature,
         419             :   ! water vapor, cloud ice and cloud liquid budgets.
         420             :   logical :: history_budget
         421             :   ! output history file number for budget fields
         422             :   integer :: history_budget_histfile_num
         423             :   ! output variables of interest in WACCM runs
         424             :   logical :: history_waccm
         425             :
         426             :   ! Interpolated Newtonian cooling coefficients.
         427             :   real(r8) :: alpha(pver+1)
         428             :
         429             :   ! Levels of pre-calculated Newtonian cooling (1/day).
         430             :   ! The following profile is digitized from:
         431             :   ! Wehrbein and Leovy (JAS, 39, 1532-1544, 1982) figure 5
         432             :
         433             :   integer, parameter :: nalph = 71
         434             :   real(r8) :: alpha0(nalph) = [ &
         435             :        0.1_r8,         0.1_r8,         0.1_r8,         0.1_r8,         &
         436             :        0.1_r8,         0.1_r8,         0.1_r8,         0.1_r8,         &
         437             :        0.1_r8,         0.1_r8,         0.10133333_r8,  0.104_r8,       &
         438             :        0.108_r8,       0.112_r8,       0.116_r8,       0.12066667_r8,  &
         439             :        0.126_r8,       0.132_r8,       0.138_r8,       0.144_r8,       &
         440             :        0.15133333_r8,  0.16_r8,        0.17_r8,        0.18_r8,        &
         441             :        0.19_r8,        0.19933333_r8,  0.208_r8,       0.216_r8,       &
         442             :        0.224_r8,       0.232_r8,       0.23466667_r8,  0.232_r8,       &
         443             :        0.224_r8,       0.216_r8,       0.208_r8,       0.20133333_r8,  &
         444             :        0.196_r8,       0.192_r8,       0.188_r8,       0.184_r8,       &
         445             :        0.18266667_r8,  0.184_r8,       0.188_r8,       0.192_r8,       &
         446             :        0.196_r8,       0.19333333_r8,  0.184_r8,       0.168_r8,       &
         447             :        0.152_r8,       0.136_r8,       0.12133333_r8,  0.108_r8,       &
         448             :        0.096_r8,       0.084_r8,       0.072_r8,       0.061_r8,       &
         449             :        0.051_r8,       0.042_r8,       0.033_r8,       0.024_r8,       &
         450             :        0.017666667_r8, 0.014_r8,       0.013_r8,       0.012_r8,       &
         451             :        0.011_r8,       0.010333333_r8, 0.01_r8,        0.01_r8,        &
         452             :        0.01_r8,        0.01_r8,        0.01_r8                         &
         453             :        ]
         454             :
         455             :   ! Pressure levels that were used to calculate alpha0 (hPa).
         456             :   real(r8) :: palph(nalph) = [ &
         457             :        2.06115E-06_r8, 2.74280E-06_r8, 3.64988E-06_r8, 4.85694E-06_r8, &
         458             :        6.46319E-06_r8, 8.60065E-06_r8, 1.14450E-05_r8, 1.52300E-05_r8, &
         459             :        2.02667E-05_r8, 2.69692E-05_r8, 3.58882E-05_r8, 4.77568E-05_r8, &
         460             :        6.35507E-05_r8, 8.45676E-05_r8, 0.000112535_r8, 0.000149752_r8, &
         461             :        0.000199277_r8, 0.000265180_r8, 0.000352878_r8, 0.000469579_r8, &
         462             :        0.000624875_r8, 0.000831529_r8, 0.00110653_r8,  0.00147247_r8,  &
         463             :        0.00195943_r8,  0.00260744_r8,  0.00346975_r8,  0.00461724_r8,  &
         464             :        0.00614421_r8,  0.00817618_r8,  0.0108801_r8,   0.0144783_r8,   &
         465             :        0.0192665_r8,   0.0256382_r8,   0.0341170_r8,   0.0453999_r8,   &
         466             :        0.0604142_r8,   0.0803939_r8,   0.106981_r8,    0.142361_r8,    &
         467             :        0.189442_r8,    0.252093_r8,    0.335463_r8,    0.446404_r8,    &
         468             :        0.594036_r8,    0.790490_r8,    1.05192_r8,     1.39980_r8,     &
         469             :        1.86273_r8,     2.47875_r8,     3.29851_r8,     4.38936_r8,     &
         470             :        5.84098_r8,     7.77266_r8,     10.3432_r8,     13.7638_r8,     &
         471             :        18.3156_r8,     24.3728_r8,     32.4332_r8,     43.1593_r8,     &
         472             :        57.4326_r8,     76.4263_r8,     101.701_r8,     135.335_r8,     &
         473             :        180.092_r8,     239.651_r8,     318.907_r8,     424.373_r8,     &
         474             :        564.718_r8,     751.477_r8,     1000._r8                        &
         475             :        ]
         476             :
         477             :   ! Read data from file
         478             :   type(file_desc_t), pointer :: fh_topo
         479             :   type(file_desc_t)  :: fh_rdggm
         480             :   integer            :: grid_id
         481             :   character(len=8)   :: dim1name, dim2name
         482             :   logical            :: found
         483             :   character(len=cl) :: bnd_rdggm_loc   ! filepath of topo file on local disk
         484             :
         485             :   ! Allow reporting of error messages.
         486             :   character(len=128) :: errstring
         487             :   character(len=*), parameter :: sub = 'gw_init'
         488             :
         489             :   ! temporary workaround for restart w/ ridge scheme
         490             :   character(len=cl) :: bnd_topo_loc   ! filepath of topo file on local disk
         491             :
         492             :   integer :: botndx,topndx
         493             :
         494             :   !-----------------------------------------------------------------------
         495             :
         496        1536 :   if (do_molec_diff) then
         497           0 :      kvt_idx     = pbuf_get_index('kvt')
         498             :   end if
         499             :
         500        1536 :   if (masterproc) then
         501           2 :      write(iulog,*) ' '
         502           2 :      write(iulog,*) "GW_DRAG: band_mid%ngwv = ", band_mid%ngwv
         503           4 :      do l = -band_mid%ngwv, band_mid%ngwv
         504             :         write (iulog,'(A,I0,A,F7.2)') &
         505           4 :              "GW_DRAG: band_mid%cref(",l,") = ",band_mid%cref(l)
         506             :      enddo
         507           2 :      write(iulog,*) 'GW_DRAG: band_mid%kwv = ', band_mid%kwv
         508           2 :      write(iulog,*) 'GW_DRAG: band_mid%fcrit2 = ', band_mid%fcrit2
         509           2 :      write(iulog,*) ' '
         510           2 :      write(iulog,*) "GW_DRAG: band_long%ngwv = ", band_long%ngwv
         511           4 :      do l = -band_long%ngwv, band_long%ngwv
         512             :         write (iulog,'(A,I2,A,F7.2)') &
         513           4 :              "GW_DRAG: band_long%cref(",l,") = ",band_long%cref(l)
         514             :      enddo
         515           2 :      write(iulog,*) 'GW_DRAG: band_long%kwv = ', band_long%kwv
         516           2 :      write(iulog,*) 'GW_DRAG: band_long%fcrit2 = ', band_long%fcrit2
         517           2 :      write(iulog,*) ' '
         518             :   end if
         519             :
         520             :   ! pre-calculated newtonian damping:
         521             :   !     * convert to s-1
         522             :   !     * ensure it is not smaller than 1e-6
         523             :   !     * convert palph from hpa to pa
         524             :
         525      110592 :   do k=1,nalph
         526      109056 :      alpha0(k) = alpha0(k) / 86400._r8
         527      109056 :      alpha0(k) = max(alpha0(k), 1.e-6_r8)
         528      110592 :      palph(k) = palph(k)*1.e2_r8
         529             :   end do
         530             :
         531             :   ! interpolate to current vertical grid and obtain alpha
         532             :
         533        1536 :   call lininterp (alpha0  ,palph, nalph , alpha  , pref_edge , pver+1)
         534        1536 :   if (masterproc) then
         535           2 :      write (iulog,*) 'gw_init: newtonian damping (1/day):'
         536           2 :      write (iulog,fmt='(a4,a12,a10)') ' k  ','  pref_edge      ', &
         537           4 :           '  alpha   '
         538          56 :      do k = 1, pver+1
         539          54 :         write (iulog,fmt='(i4,1e12.5,1f10.2)') k,pref_edge(k), &
         540         110 :              alpha(k)*86400._r8
         541             :      end do
         542             :   end if
         543             :
         544        1536 :   if (masterproc) then
         545           2 :      write(iulog,*) 'KTOP        =',ktop
         546             :   end if
         547             :
         548             :   ! Used to decide whether temperature tendencies should be output.
         549             :   call phys_getopts( history_budget_out = history_budget, &
         550             :        history_budget_histfile_num_out = history_budget_histfile_num, &
         551             :        history_waccm_out = history_waccm, &
         552        1536 :        history_amwg_out   = history_amwg  )
         553             :
         554             :   ! Initialize subordinate modules.
         555             :   call gw_common_init(pver,&
         556             :        tau_0_ubc, ktop, gravit, rair, alpha, gw_prndl, &
         557        1536 :        gw_qbo_hdepth_scaling, errstring )
         558             :   call shr_assert(trim(errstring) == "", "gw_common_init: "//errstring// &
         559        1536 :        errMsg(__FILE__, __LINE__))
         560             :
         561        1536 :   if ( use_gw_oro ) then
         562             :
         563        1536 :      if (effgw_oro == unset_r8) then
         564             :         call endrun("gw_init: Orographic gravity waves enabled, &
         565           0 :              &but effgw_oro was not set.")
         566             :      end if
         567             :   end if
         568             :
         569        1536 :   if (use_gw_oro .or. use_gw_rdg_beta .or. use_gw_rdg_gamma) then
         570             :
         571        1536 :      sgh_idx = pbuf_get_index('SGH')
         572             :
         573             :      ! Declare history variables for orographic term
         574             :      call addfld ('TAUAORO',    (/ 'ilev' /), 'I','N m-2',  &
         575        3072 :           'Total stress from original OGW scheme')
         576             :      call addfld ('TTGWORO',    (/ 'lev' /), 'A','K s-1',  &
         577        3072 :           'T tendency - orographic gravity wave drag')
         578             :      call addfld ('TTGWSDFORO', (/ 'lev' /), 'A','K s-1',  &
         579        3072 :           'T tendency - orographic gravity wave, diffusion.')
         580             :      call addfld ('TTGWSKEORO', (/ 'lev' /), 'A','K s-1',  &
         581        3072 :           'T tendency - orographic gravity wave, breaking KE.')
         582             :      call addfld ('UTGWORO',    (/ 'lev' /), 'A','m s-2', &
         583        3072 :           'U tendency - orographic gravity wave drag')
         584             :      call addfld ('VTGWORO',    (/ 'lev' /), 'A','m s-2', &
         585        3072 :           'V tendency - orographic gravity wave drag')
         586        1536 :      call register_vector_field('UTGWORO', 'VTGWORO')
         587             :      call addfld ('TAUGWX',     horiz_only,  'A','N m-2', &
         588        1536 :           'Zonal gravity wave surface stress')
         589             :      call addfld ('TAUGWY',     horiz_only,  'A','N m-2', &
         590        1536 :           'Meridional gravity wave surface stress')
         591        1536 :      call register_vector_field('TAUGWX', 'TAUGWY')
         592             :
         593        1536 :      if (history_amwg) then
         594        1536 :         call add_default('TAUGWX  ', 1, ' ')
         595        1536 :         call add_default('TAUGWY  ', 1, ' ')
         596             :      end if
         597             :
         598        1536 :      if (history_waccm) then
         599           0 :         call add_default('UTGWORO ', 1, ' ')
         600           0 :         call add_default('VTGWORO ', 1, ' ')
         601           0 :         call add_default('TAUGWX  ', 1, ' ')
         602           0 :         call add_default('TAUGWY  ', 1, ' ')
         603             :      end if
         604             :
         605             :   end if
         606             :
         607        1536 :   if (use_gw_rdg_beta .or. use_gw_rdg_gamma) then
         608           0 :      grid_id = cam_grid_id('physgrid')
         609           0 :      if (.not. cam_grid_check(grid_id)) then
         610           0 :         call endrun(sub//': ERROR: no "physgrid" grid')
         611             :      end if
         612           0 :      call cam_grid_get_dim_names(grid_id, dim1name, dim2name)
         613             :   end if
         614             :
         615        1536 :   if (use_gw_rdg_beta) then
         616             :
         617           0 :      if (effgw_rdg_beta == unset_r8) then
         618             :         call endrun(sub//": ERROR: Anisotropic OGW enabled, &
         619           0 :                      &but effgw_rdg_beta was not set.")
         620             :      end if
         621             :
         622           0 :      fh_topo => topo_file_get_id()
         623           0 :      bnd_topo_loc = ' '
         624           0 :      if (.not. associated(fh_topo)) then
         625             :
         626             :         ! Try to open topo file here.  This workaround will not be needed
         627             :         ! once the refactored initialization sequence is on trunk.
         628             :
         629           0 :         allocate(fh_topo)
         630             :         ! Error exit is from getfil if file not found.
         631           0 :         call getfil(bnd_topo, bnd_topo_loc)
         632           0 :         call cam_pio_openfile(fh_topo, bnd_topo_loc, PIO_NOWRITE)
         633             :
         634             :      end if
         635             :
         636             :      ! Get beta ridge data
         637             :      allocate( &
         638           0 :         rdg_gbxar(pcols,begchunk:endchunk),      &
         639           0 :         rdg_hwdth(pcols,prdg,begchunk:endchunk), &
         640           0 :         rdg_clngt(pcols,prdg,begchunk:endchunk), &
         641           0 :         rdg_mxdis(pcols,prdg,begchunk:endchunk), &
         642           0 :         rdg_anixy(pcols,prdg,begchunk:endchunk), &
         643           0 :         rdg_angll(pcols,prdg,begchunk:endchunk)  )
         644             :
         645             :      call infld('GBXAR', fh_topo, dim1name, dim2name, 1, pcols, &
         646           0 :                          begchunk, endchunk, rdg_gbxar, found, gridname='physgrid')
         647           0 :      if (.not. found) call endrun(sub//': ERROR: GBXAR not found on topo file')
         648           0 :      rdg_gbxar = rdg_gbxar * (rearth/1000._r8)*(rearth/1000._r8) ! transform to km^2
         649             :
         650             :      call infld('HWDTH', fh_topo, dim1name, 'nrdg', dim2name, 1, pcols, &
         651           0 :                 1, prdg, begchunk, endchunk, rdg_hwdth, found, gridname='physgrid')
         652           0 :      if (.not. found) call endrun(sub//': ERROR: HWDTH not found on topo file')
         653             :
         654             :      call infld('CLNGT', fh_topo, dim1name, 'nrdg', dim2name, 1, pcols, &
         655           0 :                 1, prdg, begchunk, endchunk, rdg_clngt, found, gridname='physgrid')
         656           0 :      if (.not. found) call endrun(sub//': ERROR: CLNGT not found on topo file')
         657             :
         658             :      call infld('MXDIS', fh_topo, dim1name, 'nrdg', dim2name, 1, pcols, &
         659           0 :                 1, prdg, begchunk, endchunk, rdg_mxdis, found, gridname='physgrid')
         660           0 :      if (.not. found) call endrun(sub//': ERROR: MXDIS not found on topo file')
         661             :
         662             :      call infld('ANIXY', fh_topo, dim1name, 'nrdg', dim2name, 1, pcols, &
         663           0 :                 1, prdg, begchunk, endchunk, rdg_anixy, found, gridname='physgrid')
         664           0 :      if (.not. found) call endrun(sub//': ERROR: ANIXY not found on topo file')
         665             :
         666             :      call infld('ANGLL', fh_topo, dim1name, 'nrdg', dim2name, 1, pcols, &
         667           0 :                 1, prdg, begchunk, endchunk, rdg_angll, found, gridname='physgrid')
         668           0 :      if (.not. found) call endrun(sub//': ERROR: ANGLL not found on topo file')
         669             :
         670             :      ! close topo file only if it was opened here
         671           0 :      if (len_trim(bnd_topo_loc) > 0) then
         672           0 :         call pio_closefile(fh_topo)
         673             :      end if
         674             :
         675             :      call addfld ('WBR_HT1',     horiz_only,  'I','m', &
         676           0 :           'Wave breaking height for DSW')
         677             :      call addfld ('TLB_HT1',     horiz_only,  'I','m', &
         678           0 :           'Form drag layer height')
         679             :      call addfld ('BWV_HT1',     horiz_only,  'I','m', &
         680           0 :           'Bottom of freely-propagating OGW regime')
         681             :      call addfld ('TAUDSW1',     horiz_only,  'I','Nm-2', &
         682           0 :           'DSW enhanced drag')
         683             :      call addfld ('TAUORO1',     horiz_only,  'I','Nm-2', &
         684           0 :           'lower BC on propagating wave stress')
         685             :      call addfld ('UBMSRC1',     horiz_only,  'I','ms-1', &
         686           0 :           'below-peak-level on-ridge wind')
         687             :      call addfld ('USRC1',     horiz_only,  'I','ms-1', &
         688           0 :           'below-peak-level Zonal wind')
         689             :      call addfld ('VSRC1',     horiz_only,  'I','ms-1', &
         690           0 :           'below-peak-level Meridional wind')
         691             :      call addfld ('NSRC1',     horiz_only,  'I','s-1', &
         692           0 :           'below-peak-level stratification')
         693             :      call addfld ('MXDIS1',     horiz_only,  'I','m', &
         694           0 :           'Ridge/obstacle height')
         695             :      call addfld ('ANGLL1',     horiz_only,  'I','degrees', &
         696           0 :           'orientation clockwise w/resp north-south')
         697             :      call addfld ('ANIXY1',     horiz_only,  'I','1', &
         698           0 :           'Ridge quality')
         699             :      call addfld ('HWDTH1',     horiz_only,  'I','km', &
         700           0 :           'Ridge width')
         701             :      call addfld ('CLNGT1',     horiz_only,  'I','km', &
         702           0 :           'Ridge length')
         703             :      call addfld ('GBXAR1',     horiz_only,  'I','km+2', &
         704           0 :           'grid box area')
         705             :
         706             :      call addfld ('Fr1_DIAG',     horiz_only,  'I','1', &
         707           0 :           'Critical Froude number for linear waves')
         708             :      call addfld ('Fr2_DIAG',     horiz_only,  'I','1', &
         709           0 :           'Critical Froude number for blocked flow')
         710             :      call addfld ('Frx_DIAG',     horiz_only,  'I','1', &
         711           0 :           'Obstacle Froude Number')
         712             :
         713             :      call addfld('UEGW',  (/ 'lev' /) , 'A'  ,'s-1' ,  &
         714           0 :           'Zonal wind profile-entry to GW ' )
         715             :      call addfld('VEGW',  (/ 'lev' /) , 'A'  ,'s-1' ,  &
         716           0 :           'Merdional wind profile-entry to GW ' )
         717           0 :      call register_vector_field('UEGW','VEGW')
         718             :      call addfld('TEGW',  (/ 'lev' /) , 'A'  ,'K' ,  &
         719           0 :           'Temperature profile-entry to GW ' )
         720             :      call addfld('ZEGW',  (/ 'ilev' /) , 'A'  ,'m' ,  &
         721           0 :           'interface geopotential heights in GW code ' )
         722             :      call addfld('ZMGW',  (/ 'lev' /) , 'A'  ,'m' ,  &
         723           0 :           'midlayer geopotential heights in GW code ' )
         724             :
         725             :      call addfld('TAUM1_DIAG' , (/ 'ilev' /) , 'I'  ,'N m-2' , &
         726           0 :           'Ridge based momentum flux profile')
         727             :      call addfld('TAU1RDGBETAM' , (/ 'ilev' /) , 'I'  ,'N m-2' , &
         728           0 :           'Ridge based momentum flux profile')
         729             :      call addfld('UBM1BETA',  (/ 'lev' /) , 'A'  ,'s-1' ,  &
         730           0 :           'On-ridge wind profile           ' )
         731             :      call addfld('UBT1RDGBETA' , (/ 'lev' /) , 'I'  ,'m s-1' , &
         732           0 :           'On-ridge wind tendency from ridge 1     ')
         733             :
         734           0 :      do i = 1, 6
         735           0 :         write(cn, '(i1)') i
         736             :         call addfld('TAU'//cn//'RDGBETAY' , (/ 'ilev' /), 'I', 'N m-2', &
         737           0 :           'Ridge based momentum flux profile')
         738             :         call addfld('TAU'//cn//'RDGBETAX' , (/ 'ilev' /), 'I', 'N m-2', &
         739           0 :           'Ridge based momentum flux profile')
         740           0 :         call register_vector_field('TAU'//cn//'RDGBETAX','TAU'//cn//'RDGBETAY')
         741             :         call addfld('UT'//cn//'RDGBETA',    (/ 'lev' /),  'I', 'm s-1', &
         742           0 :           'U wind tendency from ridge '//cn)
         743             :         call addfld('VT'//cn//'RDGBETA',    (/ 'lev' /),  'I', 'm s-1', &
         744           0 :           'V wind tendency from ridge '//cn)
         745           0 :         call register_vector_field('UT'//cn//'RDGBETA','VT'//cn//'RDGBETA')
         746             :      end do
         747             :
         748             :      call addfld('TAUARDGBETAY' , (/ 'ilev' /) , 'I'  ,'N m-2' , &
         749           0 :           'Ridge based momentum flux profile')
         750             :      call addfld('TAUARDGBETAX' , (/ 'ilev' /) , 'I'  ,'N m-2' , &
         751           0 :           'Ridge based momentum flux profile')
         752           0 :      call register_vector_field('TAUARDGBETAX','TAUARDGBETAY')
         753             :
         754           0 :      if (history_waccm) then
         755           0 :         call add_default('TAUARDGBETAX', 1, ' ')
         756           0 :         call add_default('TAUARDGBETAY  ', 1, ' ')
         757             :      end if
         758             :
         759             :   end if
         760             :
         761        1536 :   if (use_gw_rdg_gamma) then
         762             :
         763           0 :      if (effgw_rdg_gamma == unset_r8) then
         764           0 :         call endrun(sub//": ERROR: Anisotropic OGW enabled, but effgw_rdg_gamma was not set.")
         765             :      end if
         766             :
         767           0 :      call getfil(bnd_rdggm, bnd_rdggm_loc, iflag=1, lexist=found)
         768           0 :      if (found) then
         769           0 :         call cam_pio_openfile(fh_rdggm, bnd_rdggm_loc, PIO_NOWRITE)
         770             :      else
         771             :         call endrun(sub//': ERROR: file for gamma ridges not found: bnd_rdggm='// &
         772           0 :                    trim(bnd_rdggm))
         773             :      end if
         774             :
         775           0 :      if (.not. allocated(rdg_gbxar)) then
         776           0 :         allocate(rdg_gbxar(pcols,begchunk:endchunk))
         777             :         call infld('GBXAR', fh_rdggm, dim1name, dim2name, 1, pcols, &
         778           0 :                             begchunk, endchunk, rdg_gbxar, found, gridname='physgrid')
         779           0 :         if (.not. found) call endrun(sub//': ERROR: GBXAR not found on bnd_rdggm')
         780           0 :         rdg_gbxar = rdg_gbxar * (rearth/1000._r8)*(rearth/1000._r8) ! transform to km^2
         781             :      end if
         782             :
         783             :      ! Get meso-gamma ridge data
         784             :      allocate( &
         785           0 :         rdg_hwdthg(pcols,prdg,begchunk:endchunk), &
         786           0 :         rdg_clngtg(pcols,prdg,begchunk:endchunk), &
         787           0 :         rdg_mxdisg(pcols,prdg,begchunk:endchunk), &
         788           0 :         rdg_anixyg(pcols,prdg,begchunk:endchunk), &
         789           0 :         rdg_angllg(pcols,prdg,begchunk:endchunk)  )
         790             :
         791             :      call infld('HWDTH', fh_rdggm, dim1name, 'nrdg', dim2name, 1, pcols, &
         792           0 :                 1, prdg, begchunk, endchunk, rdg_hwdthg, found, gridname='physgrid')
         793           0 :      if (.not. found) call endrun(sub//': ERROR: HWDTH not found on bnd_rdggm')
         794             :
         795             :      call infld('CLNGT', fh_rdggm, dim1name, 'nrdg', dim2name, 1, pcols, &
         796           0 :                 1, prdg, begchunk, endchunk, rdg_clngtg, found, gridname='physgrid')
         797           0 :      if (.not. found) call endrun(sub//': ERROR: CLNGT not found on bnd_rdggm')
         798             :
         799             :      call infld('MXDIS', fh_rdggm, dim1name, 'nrdg', dim2name, 1, pcols, &
         800           0 :                 1, prdg, begchunk, endchunk, rdg_mxdisg, found, gridname='physgrid')
         801           0 :      if (.not. found) call endrun(sub//': ERROR: MXDIS not found on bnd_rdggm')
         802             :
         803             :      call infld('ANIXY', fh_rdggm, dim1name, 'nrdg', dim2name, 1, pcols, &
         804           0 :                 1, prdg, begchunk, endchunk, rdg_anixyg, found, gridname='physgrid')
         805           0 :      if (.not. found) call endrun(sub//': ERROR: ANIXY not found on bnd_rdggm')
         806             :
         807             :      call infld('ANGLL', fh_rdggm, dim1name, 'nrdg', dim2name, 1, pcols, &
         808           0 :                 1, prdg, begchunk, endchunk, rdg_angllg, found, gridname='physgrid')
         809           0 :      if (.not. found) call endrun(sub//': ERROR: ANGLL not found on bnd_rdggm')
         810             :
         811           0 :      call pio_closefile(fh_rdggm)
         812             :
         813             :      call addfld ('TAU1RDGGAMMAM' , (/ 'ilev' /) , 'I'  ,'N m-2' , &
         814           0 :           'Ridge based momentum flux profile')
         815             :      call addfld ('UBM1GAMMA',  (/ 'lev' /) , 'A'  ,'s-1' ,  &
         816           0 :           'On-ridge wind profile           ' )
         817             :      call addfld ('UBT1RDGGAMMA' , (/ 'lev' /) , 'I'  ,'m s-1' , &
         818           0 :           'On-ridge wind tendency from ridge 1     ')
         819             :
         820           0 :      do i = 1, 6
         821           0 :         write(cn, '(i1)') i
         822             :         call addfld('TAU'//cn//'RDGGAMMAY', (/ 'ilev' /), 'I', 'N m-2', &
         823           0 :           'Ridge based momentum flux profile')
         824             :         call addfld('TAU'//cn//'RDGGAMMAX', (/ 'ilev' /), 'I', 'N m-2', &
         825           0 :           'Ridge based momentum flux profile')
         826             :         call addfld('UT'//cn//'RDGGAMMA' , (/ 'lev' /),  'I', 'm s-1', &
         827           0 :           'U wind tendency from ridge '//cn)
         828             :         call addfld('VT'//cn//'RDGGAMMA' , (/ 'lev' /),  'I', 'm s-1', &
         829           0 :           'V wind tendency from ridge '//cn)
         830           0 :         call register_vector_field('UT'//cn//'RDGGAMMA','VT'//cn//'RDGGAMMA')
         831             :      end do
         832             :
         833             :      call addfld ('TAUARDGGAMMAY' , (/ 'ilev' /) , 'I'  ,'N m-2' , &
         834           0 :           'Ridge based momentum flux profile')
         835             :      call addfld ('TAUARDGGAMMAX' , (/ 'ilev' /) , 'I'  ,'N m-2' , &
         836           0 :           'Ridge based momentum flux profile')
         837           0 :      call register_vector_field('TAUARDGGAMMAX','TAUARDGGAMMAY')
         838             :      call addfld ('TAURDGGMX',     horiz_only,  'A','N m-2', &
         839           0 :           'Zonal gravity wave surface stress')
         840             :      call addfld ('TAURDGGMY',     horiz_only,  'A','N m-2', &
         841           0 :           'Meridional gravity wave surface stress')
         842           0 :      call register_vector_field('TAURDGGMX','TAURDGGMY')
         843             :      call addfld ('UTRDGGM' , (/ 'lev' /) , 'I'  ,'m s-1' , &
         844           0 :           'U wind tendency from ridge 6     ')
         845             :      call addfld ('VTRDGGM' , (/ 'lev' /) , 'I'  ,'m s-1' , &
         846           0 :           'V wind tendency from ridge 6     ')
         847           0 :      call register_vector_field('UTRDGGM','VTRDGGM')
         848             :   end if
         849             :
         850        1536 :   if (use_gw_front .or. use_gw_front_igw) then
         851             :
         852           0 :      frontgf_idx = pbuf_get_index('FRONTGF')
         853           0 :      frontga_idx = pbuf_get_index('FRONTGA')
         854             :
         855             :      call shr_assert(unset_r8 /= frontgfc, &
         856             :           "gw_init: Frontogenesis enabled, but frontgfc was &
         857             :           & not set!"// &
         858           0 :           errMsg(__FILE__, __LINE__))
         859             :
         860           0 :      do k = 0, pver
         861             :         ! Check frontogenesis at 600 hPa.
         862           0 :         if (pref_edge(k+1) < 60000._r8) kfront = k+1
         863             :      end do
         864             :
         865             :      ! Source waves from 500 hPa.
         866           0 :      kbot_front = maxloc(pref_edge, 1, (pref_edge < 50000._r8)) - 1
         867             :
         868           0 :      if (masterproc) then
         869           0 :         write (iulog,*) 'KFRONT      =',kfront
         870           0 :         write (iulog,*) 'KBOT_FRONT  =',kbot_front
         871           0 :         write(iulog,*) ' '
         872             :      end if
         873             :
         874             :      call addfld ('FRONTGF', (/ 'lev' /), 'A', 'K^2/M^2/S', &
         875           0 :           'Frontogenesis function at gws src level')
         876             :      call addfld ('FRONTGFA', (/ 'lev' /), 'A', 'K^2/M^2/S', &
         877           0 :           'Frontogenesis function at gws src level')
         878             :
         879           0 :      if (history_waccm) then
         880           0 :         call add_default('FRONTGF', 1, ' ')
         881           0 :         call add_default('FRONTGFA', 1, ' ')
         882             :      end if
         883             :
         884             :   end if
         885             :
         886        1536 :   if (use_gw_front) then
         887             :
         888             :      call shr_assert(all(unset_r8 /= [ effgw_cm, taubgnd ]), &
         889             :           "gw_init: Frontogenesis mid-scale waves enabled, but not &
         890             :           &all required namelist variables were set!"// &
         891           0 :           errMsg(__FILE__, __LINE__))
         892             :
         893           0 :      if (masterproc) then
         894           0 :         write(iulog,*) 'gw_init: gw spectrum taubgnd, ', &
         895           0 :              'effgw_cm = ',taubgnd, effgw_cm
         896           0 :         write(iulog,*) ' '
         897             :      end if
         898             :
         899             :      cm_desc = gaussian_cm_desc(band_mid, kbot_front, kfront, frontgfc, &
         900           0 :           taubgnd, front_gaussian_width)
         901             :
         902             :      ! Output for gravity waves from frontogenesis.
         903             :      call gw_spec_addflds(prefix=cm_pf, scheme="C&M", band=band_mid, &
         904           0 :           history_defaults=history_waccm)
         905             :
         906             :   end if
         907             :
         908        1536 :   if (use_gw_front_igw) then
         909             :
         910             :      call shr_assert(all(unset_r8 /= [ effgw_cm_igw, taubgnd_igw ]), &
         911             :           "gw_init: Frontogenesis inertial waves enabled, but not &
         912             :           &all required namelist variables were set!"// &
         913           0 :           errMsg(__FILE__, __LINE__))
         914             :
         915           0 :      if (masterproc) then
         916           0 :         write(iulog,*) 'gw_init: gw spectrum taubgnd_igw, ', &
         917           0 :              'effgw_cm_igw = ',taubgnd_igw, effgw_cm_igw
         918           0 :         write(iulog,*) ' '
         919             :      end if
         920             :
         921             :      cm_igw_desc = gaussian_cm_desc(band_long, kbot_front, kfront, frontgfc, &
         922           0 :           taubgnd_igw, front_gaussian_width)
         923             :
         924             :      ! Output for gravity waves from frontogenesis.
         925             :      call gw_spec_addflds(prefix=cm_igw_pf, scheme="C&M IGW", &
         926           0 :           band=band_long, history_defaults=history_waccm)
         927             :
         928             :   end if
         929             :
         930             :   ! ========= Moving Mountain initialization! ==========================
         931        1536 :   if (use_gw_movmtn_pbl) then
         932             :
         933             :      ! get pbuf indices for CLUBB couplings
         934           0 :      ttend_clubb_idx     = pbuf_get_index('TTEND_CLUBB')
         935           0 :      thlp2_clubb_gw_idx  = pbuf_get_index('THLP2_CLUBB_GW')
         936           0 :      upwp_clubb_gw_idx   = pbuf_get_index('UPWP_CLUBB_GW')
         937           0 :      vpwp_clubb_gw_idx   = pbuf_get_index('VPWP_CLUBB_GW')
         938           0 :      wpthlp_clubb_gw_idx  = pbuf_get_index('WPTHLP_CLUBB_GW')
         939             :
         940           0 :      if (masterproc) then
         941           0 :         write (iulog,*) 'Moving Mountain development code call init_movmtn'
         942             :      end if
         943             :
         944             :
         945             :      ! Confirm moving mountain file is enabled
         946             :      call shr_assert(trim(gw_drag_file_mm) /= "", &
         947             :           "gw_init: No gw_drag_file provided for DP GW moving mountain lookup &
         948             :           &table. Set this via namelist."// &
         949           0 :           errMsg(__FILE__, __LINE__))
         950             :
         951           0 :      call gw_init_movmtn(gw_drag_file_mm, band_movmtn, movmtn_desc)
         952             :
         953           0 :      do k = 0, pver
         954             :         ! 950 hPa index
         955           0 :         if (pref_edge(k+1) < 95000._r8) movmtn_desc%k = k+1
         956             :      end do
         957             :
         958             :     ! Don't use deep convection heating depths below this limit.
         959           0 :      movmtn_desc%min_hdepth = 1._r8
         960           0 :      if (masterproc) then
         961           0 :         write (iulog,*) 'Moving mountain deep level =',movmtn_desc%k
         962             :      end if
         963             :
         964             :      call addfld ('GWUT_MOVMTN',(/ 'lev' /), 'I','m s-2', &
         965           0 :           'Mov Mtn dragforce - ubm component')
         966             :      call addfld ('UTGW_MOVMTN',(/ 'lev' /), 'I','m s-2', &
         967           0 :           'Mov Mtn dragforce - u component')
         968             :      call addfld ('VTGW_MOVMTN',(/ 'lev' /), 'I','m s-2', &
         969           0 :           'Mov Mtn dragforce - v component')
         970             :      call addfld('TAU_MOVMTN', (/ 'ilev' /), 'I', 'N m-2', &
         971           0 :           'Moving Mountain momentum flux profile')
         972             :      call addfld('U_MOVMTN_IN', (/ 'lev' /), 'I', 'm s-1', &
         973           0 :           'Moving Mountain - midpoint zonal input wind')
         974             :      call addfld('V_MOVMTN_IN', (/ 'lev' /), 'I', 'm s-1', &
         975           0 :           'Moving Mountain - midpoint meridional input wind')
         976             :      call addfld('UBI_MOVMTN', (/ 'ilev' /), 'I', 'm s-1', &
         977           0 :           'Moving Mountain - interface wind in direction of wave')
         978             :      call addfld('UBM_MOVMTN', (/ 'lev' /), 'I', 'm s-1', &
         979           0 :           'Moving Mountain - midpoint wind in direction of wave')
         980             :      call addfld ('HDEPTH_MOVMTN',horiz_only,'I','km', &
         981           0 :           'Heating Depth')
         982             :      call addfld ('UCELL_MOVMTN',horiz_only,'I','m s-1', &
         983           0 :           'Gravity Wave Moving Mountain - Source-level X-wind')
         984             :      call addfld ('VCELL_MOVMTN',horiz_only,'I','m s-1', &
         985           0 :           'Gravity Wave Moving Mountain - Source-level Y-wind')
         986             :      call addfld ('CS_MOVMTN',horiz_only,'I','m s-1', &
         987           0 :           'Gravity Wave Moving Mountain - phase speed in direction of wave')
         988             :      call addfld ('STEER_LEVEL_MOVMTN',horiz_only,'I','1', &
         989           0 :           'Gravity Wave Moving Mountain - steering level for movmtn GW')
         990             :      call addfld ('SRC_LEVEL_MOVMTN',horiz_only,'I','1', &
         991           0 :           'Gravity Wave Moving Mountain - launch level for movmtn GW')
         992             :      call addfld ('TND_LEVEL_MOVMTN',horiz_only,'I','1', &
         993           0 :           'Gravity Wave Moving Mountain - tendency lowest level for movmtn GW')
         994             :      call addfld ('NETDT_MOVMTN',(/ 'lev' /),'I','K s-1', &
         995           0 :           'Gravity Wave Moving Mountain - Net heating rate')
         996             :      call addfld ('TTEND_CLUBB',(/ 'lev' /),'A','K s-1', &
         997           0 :           'Gravity Wave Moving Mountain - CLUBB Net heating rate')
         998             :      call addfld ('THLP2_CLUBB_GW',(/ 'ilev' /),'A','K+2', &
         999           0 :           'Gravity Wave Moving Mountain - THLP variance from CLUBB to GW')
        1000             :      call addfld ('WPTHLP_CLUBB_GW',(/ 'ilev' /),'A','Km s-2', &
        1001           0 :           'Gravity Wave Moving Mountain - WPTHLP from CLUBB to GW')
        1002             :      call addfld ('UPWP_CLUBB_GW',(/ 'ilev' /),'A','m+2 s-2', &
        1003           0 :           'Gravity Wave Moving Mountain - X-momflux from CLUBB to GW')
        1004             :      call addfld ('VPWP_CLUBB_GW',(/ 'ilev' /),'A','m+2 s-2', &
        1005           0 :           'Gravity Wave Moving Mountain - Y-momflux from CLUBB to GW')
        1006             :      call addfld ('XPWP_SRC_MOVMTN',horiz_only,'I','m+2 s-2', &
        1007           0 :           'Gravity Wave Moving Mountain - flux source for moving mtn')
        1008             :
        1009             :   end if
        1010             :
        1011        1536 :   if (use_gw_convect_dp) then
        1012             :
        1013           0 :      ttend_dp_idx    = pbuf_get_index('TTEND_DP')
        1014             :
        1015             :      ! Set the deep scheme specification components.
        1016           0 :      beres_dp_desc%storm_shift = .true.
        1017             :
        1018           0 :      do k = 0, pver
        1019             :         ! 700 hPa index
        1020           0 :         if (pref_edge(k+1) < 70000._r8) beres_dp_desc%k = k+1
        1021             :      end do
        1022             :
        1023           0 :      if (masterproc) then
        1024           0 :         write (iulog,*) 'Beres deep level =',beres_dp_desc%k
        1025             :      end if
        1026             :
        1027             :      ! Don't use deep convection heating depths below this limit.
        1028             :      ! This is important for QBO. Bad result if left at 2.5 km.
        1029           0 :      beres_dp_desc%min_hdepth = 1000._r8
        1030             :
        1031             :      ! Read Beres file.
        1032             :
        1033             :      call shr_assert(trim(gw_drag_file) /= "", &
        1034             :           "gw_init: No gw_drag_file provided for Beres deep &
        1035             :           &scheme. Set this via namelist."// &
        1036           0 :           errMsg(__FILE__, __LINE__))
        1037             :
        1038           0 :      call gw_init_beres(gw_drag_file, band_mid, beres_dp_desc)
        1039             :
        1040             :      ! Output for gravity waves from the Beres scheme (deep).
        1041             :      call gw_spec_addflds(prefix=beres_dp_pf, scheme="Beres (deep)", &
        1042           0 :           band=band_mid, history_defaults=history_waccm)
        1043             :
        1044             :      call addfld ('NETDT',(/ 'lev' /), 'A','K s-1', &
        1045           0 :           'Net heating rate')
        1046             :      call addfld ('MAXQ0',horiz_only  ,  'A','K day-1', &
        1047           0 :           'Max column heating rate')
        1048             :      call addfld ('HDEPTH',horiz_only,    'A','km', &
        1049           0 :           'Heating Depth')
        1050             :
        1051           0 :      if (history_waccm) then
        1052           0 :         call add_default('NETDT    ', 1, ' ')
        1053           0 :         call add_default('HDEPTH   ', 1, ' ')
        1054           0 :         call add_default('MAXQ0    ', 1, ' ')
        1055             :      end if
        1056             :
        1057             :   end if
        1058             :
        1059        1536 :   if (use_gw_convect_sh) then
        1060             :
        1061           0 :      ttend_sh_idx    = pbuf_get_index('TTEND_SH')
        1062             :
        1063             :      ! Set the shallow scheme specification components.
        1064           0 :      beres_sh_desc%storm_shift = .false.
        1065             :
        1066           0 :      do k = 0, pver
        1067             :         ! 900 hPa index
        1068           0 :         if (pref_edge(k+1) < 90000._r8) beres_sh_desc%k = k+1
        1069             :      end do
        1070             :
        1071           0 :      if (masterproc) then
        1072           0 :         write (iulog,*) 'Beres shallow level =',beres_sh_desc%k
        1073             :      end if
        1074             :
        1075             :      ! Use all heating depths for shallow convection.
        1076           0 :      beres_sh_desc%min_hdepth = 0._r8
        1077             :
        1078             :      ! Read Beres file.
        1079             :
        1080             :      call shr_assert(trim(gw_drag_file_sh) /= "", &
        1081             :           "gw_init: No gw_drag_file_sh provided for Beres shallow &
        1082             :           &scheme. Set this via namelist."// &
        1083           0 :           errMsg(__FILE__, __LINE__))
        1084             :
        1085           0 :      call gw_init_beres(gw_drag_file_sh, band_mid, beres_sh_desc)
        1086             :
        1087             :      ! Output for gravity waves from the Beres scheme (shallow).
        1088             :      call gw_spec_addflds(prefix=beres_sh_pf, scheme="Beres (shallow)", &
        1089           0 :           band=band_mid, history_defaults=history_waccm)
        1090             :
        1091             :      call addfld ('SNETDT',(/ 'lev' /), 'A','K s-1', &
        1092           0 :           'Net heating rate')
        1093             :      call addfld ('SMAXQ0',horiz_only  ,  'A','K day-1', &
        1094           0 :           'Max column heating rate')
        1095             :      call addfld ('SHDEPTH',horiz_only,    'A','km', &
        1096           0 :           'Heating Depth')
        1097             :
        1098           0 :      if (history_waccm) then
        1099           0 :         call add_default('SNETDT   ', 1, ' ')
        1100           0 :         call add_default('SHDEPTH  ', 1, ' ')
        1101           0 :         call add_default('SMAXQ0   ', 1, ' ')
        1102             :      end if
        1103             :
        1104             :   end if
        1105             :
        1106             :   call addfld ('EKGW' ,(/ 'ilev' /), 'A','M2/S', &
        1107        3072 :        'Effective Kzz due to diffusion by gravity waves')
        1108             :
        1109        1536 :   if (history_waccm) then
        1110           0 :      call add_default('EKGW', 1, ' ')
        1111             :   end if
        1112             :
        1113             :   call addfld ('UTGW_TOTAL',    (/ 'lev' /), 'A','m s-2', &
        1114        3072 :        'Total U tendency due to gravity wave drag')
        1115             :   call addfld ('VTGW_TOTAL',    (/ 'lev' /), 'A','m s-2', &
        1116        3072 :        'Total V tendency due to gravity wave drag')
        1117        1536 :   call register_vector_field('UTGW_TOTAL', 'VTGW_TOTAL')
        1118             :
        1119             :   ! Total temperature tendency output.
        1120             :   call addfld ('TTGW', (/ 'lev' /), 'A', 'K s-1',  &
        1121        3072 :        'T tendency - gravity wave drag')
        1122             :
        1123             :   ! Water budget terms.
        1124             :   call addfld('QTGW',(/ 'lev' /), 'A','kg/kg/s', &
        1125        3072 :        'Q tendency - gravity wave drag')
        1126             :   call addfld('CLDLIQTGW',(/ 'lev' /), 'A','kg/kg/s', &
        1127        3072 :        'CLDLIQ tendency - gravity wave drag')
        1128             :   call addfld('CLDICETGW',(/ 'lev' /), 'A','kg/kg/s', &
        1129        3072 :        'CLDICE tendency - gravity wave drag')
        1130             :
        1131        1536 :   if ( history_budget ) then
        1132           0 :      call add_default('TTGW', history_budget_histfile_num, ' ')
        1133           0 :      call add_default('QTGW', history_budget_histfile_num, ' ')
        1134           0 :      call add_default('CLDLIQTGW', history_budget_histfile_num, ' ')
        1135           0 :      call add_default('CLDICETGW', history_budget_histfile_num, ' ')
        1136             :   end if
        1137             :
        1138             :   ! Get indices to actually output the above.
        1139        1536 :   call cnst_get_ind("CLDLIQ", ixcldliq)
        1140        1536 :   call cnst_get_ind("CLDICE", ixcldice)
        1141             :
        1142        1536 :   if (gw_top_taper) then
        1143           0 :      allocate(vramp(pver))
        1144           0 :      vramp(:) = 1._r8
        1145           0 :      topndx = 1
        1146           0 :      botndx = press_lim_idx( 0.6E-02_r8, top=.true. )
        1147           0 :      if (botndx>1) then
        1148           0 :         do k=botndx,topndx,-1
        1149           0 :            vramp(k) = vramp(k+1)/(pref_edge(k+1)/pref_edge(k))
        1150             :         end do
        1151           0 :         if (masterproc) then
        1152           0 :            write(iulog,'(A)') 'GW taper coef (vramp):'
        1153           0 :            do k=1,pver
        1154           0 :               write(iulog,"('k: ',I4,' taper coef,press(Pa): ',F12.8,E12.4)") k, vramp(k), pref_mid(k)
        1155             :            enddo
        1156             :         endif
        1157             :      endif
        1158             :   end if
        1159             :
        1160        3072 : end subroutine gw_init
        1161             :
        1162             : !==========================================================================
        1163             :
        1164           0 : subroutine gw_init_beres(file_name, band, desc)
        1165             :
        1166        1536 :   use ioFileMod, only: getfil
        1167             :   use pio, only: file_desc_t, pio_nowrite, pio_inq_varid, pio_get_var, &
        1168             :        pio_closefile
        1169             :   use cam_pio_utils, only: cam_pio_openfile
        1170             :
        1171             :   character(len=*), intent(in) :: file_name
        1172             :   type(GWBand), intent(in) :: band
        1173             :
        1174             :   type(BeresSourceDesc), intent(inout) :: desc
        1175             :
        1176             :   type(file_desc_t) :: gw_file_desc
        1177             :
        1178             :   ! PIO variable ids and error code.
        1179             :   integer :: mfccid, hdid, stat
        1180             :
        1181             :   ! Number of wavenumbers in the input file.
        1182             :   integer :: ngwv_file
        1183             :
        1184             :   ! Full path to gw_drag_file.
        1185             :   character(len=cl) :: file_path
        1186             :
        1187             :   character(len=cl) :: msg
        1188             :
        1189             :   !----------------------------------------------------------------------
        1190             :   ! read in look-up table for source spectra
        1191             :   !-----------------------------------------------------------------------
        1192             :
        1193           0 :   call getfil(file_name, file_path)
        1194           0 :   call cam_pio_openfile(gw_file_desc, file_path, pio_nowrite)
        1195             :
        1196             :   ! Get HD (heating depth) dimension.
        1197             :
        1198           0 :   desc%maxh = get_pio_dimlen(gw_file_desc, "HD", file_path)
        1199             :
        1200             :   ! Get MW (mean wind) dimension.
        1201             :
        1202           0 :   desc%maxuh = get_pio_dimlen(gw_file_desc, "MW", file_path)
        1203             :
        1204             :   ! Get PS (phase speed) dimension.
        1205             :
        1206           0 :   ngwv_file = get_pio_dimlen(gw_file_desc, "PS", file_path)
        1207             :
        1208             :   ! Number in each direction is half of total (and minus phase speed of 0).
        1209           0 :   desc%maxuh = (desc%maxuh-1)/2
        1210           0 :   ngwv_file = (ngwv_file-1)/2
        1211             :
        1212             :   call shr_assert(ngwv_file >= band%ngwv, &
        1213             :        "gw_init_beres: PhaseSpeed in lookup table file does not cover the whole &
        1214           0 :        &spectrum implied by the model's ngwv. ")
        1215             :
        1216             :   ! Allocate hd and get data.
        1217             :
        1218           0 :   allocate(desc%hd(desc%maxh), stat=stat, errmsg=msg)
        1219             :
        1220             :   call shr_assert(stat == 0, &
        1221             :        "gw_init_beres: Allocation error (hd): "//msg// &
        1222           0 :        errMsg(__FILE__, __LINE__))
        1223             :
        1224           0 :   stat = pio_inq_varid(gw_file_desc,'HD',hdid)
        1225             :
        1226             :   call handle_pio_error(stat, &
        1227           0 :        'Error finding HD in: '//trim(file_path))
        1228             :
        1229             :   stat = pio_get_var(gw_file_desc, hdid, start=[1], count=[desc%maxh], &
        1230           0 :        ival=desc%hd)
        1231             :
        1232             :   call handle_pio_error(stat, &
        1233           0 :        'Error reading HD from: '//trim(file_path))
        1234             :
        1235             :   ! While not currently documented in the file, it uses kilometers. Convert
        1236             :   ! to meters.
        1237           0 :   desc%hd = desc%hd*1000._r8
        1238             :
        1239             :   ! Allocate mfcc. "desc%maxh" and "desc%maxuh" are from the file, but the
        1240             :   ! model determines wavenumber dimension.
        1241             :
        1242             :   allocate(desc%mfcc(desc%maxh,-desc%maxuh:desc%maxuh,&
        1243           0 :        -band%ngwv:band%ngwv), stat=stat, errmsg=msg)
        1244             :
        1245             :   call shr_assert(stat == 0, &
        1246             :        "gw_init_beres: Allocation error (mfcc): "//msg// &
        1247           0 :        errMsg(__FILE__, __LINE__))
        1248             :
        1249             :   ! Get mfcc data.
        1250             :
        1251           0 :   stat = pio_inq_varid(gw_file_desc,'mfcc',mfccid)
        1252             :
        1253             :   call handle_pio_error(stat, &
        1254           0 :        'Error finding mfcc in: '//trim(file_path))
        1255             :
        1256             :   stat = pio_get_var(gw_file_desc, mfccid, &
        1257           0 :        start=[1,1,ngwv_file-band%ngwv+1], count=shape(desc%mfcc), &
        1258           0 :        ival=desc%mfcc)
        1259             :
        1260             :   call handle_pio_error(stat, &
        1261           0 :        'Error reading mfcc from: '//trim(file_path))
        1262             :
        1263           0 :   call pio_closefile(gw_file_desc)
        1264             :
        1265           0 :   if (masterproc) then
        1266             :
        1267           0 :      write(iulog,*) "Read in source spectra from file."
        1268           0 :      write(iulog,*) "mfcc max, min = ", &
        1269           0 :           maxval(desc%mfcc),", ",minval(desc%mfcc)
        1270             :
        1271             :   endif
        1272             :
        1273           0 : end subroutine gw_init_beres
        1274             :
        1275             : !==============================================================
        1276           0 : subroutine gw_init_movmtn(file_name, band, desc)
        1277             :
        1278           0 :   use ioFileMod, only: getfil
        1279             :   use pio, only: file_desc_t, pio_nowrite, pio_inq_varid, pio_get_var, &
        1280             :        pio_closefile
        1281             :   use cam_pio_utils, only: cam_pio_openfile
        1282             :
        1283             :   character(len=*), intent(in) :: file_name
        1284             :   type(GWBand), intent(in) :: band
        1285             :
        1286             :   type(MovMtnSourceDesc), intent(inout) :: desc
        1287             :
        1288             :   type(file_desc_t) :: gw_file_desc
        1289             :
        1290             :   ! PIO variable ids and error code.
        1291             :   integer :: mfccid, uhid, hdid, stat
        1292             :
        1293             :   ! Number of wavenumbers in the input file.
        1294             :   integer :: ngwv_file
        1295             :
        1296             :   ! Full path to gw_drag_file.
        1297             :   character(len=cl) :: file_path
        1298             :
        1299             :   character(len=cl) :: msg
        1300             :
        1301             :   !----------------------------------------------------------------------
        1302             :   ! read in look-up table for source spectra
        1303             :   !-----------------------------------------------------------------------
        1304             :
        1305           0 :   call getfil(file_name, file_path)
        1306             :
        1307           0 :   call cam_pio_openfile(gw_file_desc, file_path, pio_nowrite)
        1308             :
        1309             :   ! Get HD (heating depth) dimension.
        1310             :
        1311           0 :   desc%maxh = 15 !get_pio_dimlen(gw_file_desc, "HD", file_path)
        1312             :
        1313             :   ! Get MW (mean wind) dimension.
        1314             :
        1315           0 :  desc%maxuh = 241 ! get_pio_dimlen(gw_file_desc, "MW", file_path)
        1316             :
        1317             :   ! Get PS (phase speed) dimension.
        1318             :
        1319           0 :   ngwv_file = 0 !get_pio_dimlen(gw_file_desc, "PS", file_path)
        1320             :
        1321             :   ! Number in each direction is half of total (and minus phase speed of 0).
        1322           0 :   desc%maxuh = (desc%maxuh-1)/2
        1323           0 :   ngwv_file = (ngwv_file-1)/2
        1324             :
        1325             :   call shr_assert(ngwv_file >= band%ngwv, &
        1326           0 :        "gw_movmtn_init: PhaseSpeed in lookup table inconsistent with moving mountain")
        1327             :
        1328             :   ! Allocate hd and get data.
        1329             :
        1330           0 :   allocate(desc%hd(desc%maxh), stat=stat, errmsg=msg)
        1331             :
        1332             :   call shr_assert(stat == 0, &
        1333             :        "gw_init_movmtn: Allocation error (hd): "//msg// &
        1334           0 :        errMsg(__FILE__, __LINE__))
        1335             :
        1336           0 :   stat = pio_inq_varid(gw_file_desc,'HDEPTH',hdid)
        1337             :
        1338             :   call handle_pio_error(stat, &
        1339           0 :        'Error finding HD in: '//trim(file_path))
        1340             :
        1341             :   stat = pio_get_var(gw_file_desc, hdid, start=[1], count=[desc%maxh], &
        1342           0 :        ival=desc%hd)
        1343             :
        1344             :   call handle_pio_error(stat, &
        1345           0 :        'Error reading HD from: '//trim(file_path))
        1346             :
        1347             :   ! While not currently documented in the file, it uses kilometers. Convert
        1348             :   ! to meters.
        1349           0 :   desc%hd = desc%hd*1000._r8
        1350             :
        1351             :  ! Allocate wind and get data.
        1352             :
        1353           0 :   allocate(desc%uh(desc%maxuh), stat=stat, errmsg=msg)
        1354             :
        1355             :   call shr_assert(stat == 0, &
        1356             :        "gw_init_movmtn: Allocation error (uh): "//msg// &
        1357           0 :        errMsg(__FILE__, __LINE__))
        1358             :
        1359           0 :   stat = pio_inq_varid(gw_file_desc,'UARR',uhid)
        1360             :
        1361             :   call handle_pio_error(stat, &
        1362           0 :        'Error finding UH in: '//trim(file_path))
        1363             :
        1364             :   stat = pio_get_var(gw_file_desc, uhid, start=[1], count=[desc%maxuh], &
        1365           0 :        ival=desc%uh)
        1366             :
        1367             :   call handle_pio_error(stat, &
        1368           0 :        'Error reading UH from: '//trim(file_path))
        1369             :
        1370             :   ! Allocate mfcc. "desc%maxh" and "desc%maxuh" are from the file, but the
        1371             :   ! model determines wavenumber dimension.
        1372             :
        1373             :   allocate(desc%mfcc(desc%maxh,-desc%maxuh:desc%maxuh,&
        1374           0 :        -band%ngwv:band%ngwv), stat=stat, errmsg=msg)
        1375             :
        1376             :   call shr_assert(stat == 0, &
        1377             :        "gw_init_movmtn: Allocation error (mfcc): "//msg// &
        1378           0 :        errMsg(__FILE__, __LINE__))
        1379             :
        1380             :   ! Get mfcc data.
        1381             :
        1382           0 :   stat = pio_inq_varid(gw_file_desc,'NEWMF',mfccid)
        1383             :
        1384             :   call handle_pio_error(stat, &
        1385           0 :        'Error finding mfcc in: '//trim(file_path))
        1386             :
        1387             :   stat = pio_get_var(gw_file_desc, mfccid, &
        1388           0 :        start=[1,1], count=shape(desc%mfcc), &
        1389           0 :        ival=desc%mfcc)
        1390             :
        1391             :   call handle_pio_error(stat, &
        1392           0 :        'Error reading mfcc from: '//trim(file_path))
        1393             :
        1394           0 :   call pio_closefile(gw_file_desc)
        1395             :
        1396           0 :   if (masterproc) then
        1397             :
        1398           0 :      write(iulog,*) "Read in Mov Mountain source file."
        1399             :
        1400             :   endif
        1401             :
        1402           0 : end subroutine gw_init_movmtn
        1403             : !==========================================================================
        1404             :
        1405             : ! Utility to reduce the repetitiveness of reads during initialization.
        1406           0 : function get_pio_dimlen(file_desc, dim_name, file_path) result(dimlen)
        1407             :
        1408           0 :   use pio, only: file_desc_t, pio_inq_dimid, pio_inq_dimlen
        1409             :
        1410             :   type(file_desc_t), intent(in) :: file_desc
        1411             :   character(len=*), intent(in) :: dim_name
        1412             :
        1413             :   ! File path, for use in error messages only.
        1414             :   character(len=*), intent(in) :: file_path
        1415             :
        1416             :   integer :: dimlen
        1417             :
        1418             :   integer :: dimid, stat
        1419             :
        1420           0 :   stat = pio_inq_dimid(file_desc, dim_name, dimid)
        1421             :
        1422             :   call handle_pio_error(stat, &
        1423           0 :        "Error finding dimension "//dim_name//" in: "//file_path)
        1424             :
        1425           0 :   stat = pio_inq_dimlen(file_desc, dimid, dimlen)
        1426             :
        1427             :   call handle_pio_error(stat, &
        1428           0 :        "Error reading dimension "//dim_name//" from: "//file_path)
        1429             :
        1430           0 : end function get_pio_dimlen
        1431             :
        1432             : !==========================================================================
        1433             :
        1434             : ! In fact, we'd usually expect PIO errors to abort the run before you can
        1435             : ! even check the error code. But just in case, use this little assert.
        1436           0 : subroutine handle_pio_error(stat, message)
        1437             :   use pio, only: pio_noerr
        1438             :   integer, intent(in) :: stat
        1439             :   character(len=*) :: message
        1440             :
        1441             :   call shr_assert(stat == pio_noerr, &
        1442             :        "PIO error in gw_init_beres: "//trim(message)// &
        1443           0 :        errMsg(__FILE__, __LINE__))
        1444             :
        1445           0 : end subroutine handle_pio_error
        1446             :
        1447             : !==========================================================================
        1448             :
        1449     5956704 : subroutine gw_tend(state, pbuf, dt, ptend, cam_in, flx_heat)
        1450             :   !-----------------------------------------------------------------------
        1451             :   ! Interface for multiple gravity wave drag parameterization.
        1452             :   !-----------------------------------------------------------------------
        1453             :
        1454             :   use physics_types,   only: physics_state_copy, set_dry_to_wet
        1455             :   use constituents,    only: cnst_type
        1456             :   use physics_buffer,  only: physics_buffer_desc, pbuf_get_field
        1457             :   use camsrfexch,      only: cam_in_t
        1458             :   ! Location-dependent cpair
        1459             :   use air_composition, only: cpairv
        1460             :   use physconst,       only: pi
        1461             :   use coords_1d,       only: Coords1D
        1462             :   use gw_common,       only: gw_prof, gw_drag_prof, calc_taucd
        1463             :   use gw_common,       only: momentum_flux, momentum_fixer, energy_change
        1464             :   use gw_common,       only: energy_fixer, coriolis_speed, adjust_inertial
        1465             :   use gw_oro,          only: gw_oro_src
        1466             :   use gw_front,        only: gw_cm_src
        1467             :   use gw_convect,      only: gw_beres_src
        1468             :   use gw_movmtn,       only: gw_movmtn_src
        1469             :
        1470             :   !------------------------------Arguments--------------------------------
        1471             :   type(physics_state), intent(in) :: state   ! physics state structure
        1472             :   type(physics_buffer_desc), pointer :: pbuf(:) ! Physics buffer
        1473             :   real(r8), intent(in) :: dt                    ! time step
        1474             :   ! Parameterization net tendencies.
        1475             :   type(physics_ptend), intent(out):: ptend
        1476             :   type(cam_in_t), intent(in) :: cam_in
        1477             :   real(r8), intent(out) :: flx_heat(pcols)
        1478             :
        1479             :   !---------------------------Local storage-------------------------------
        1480             :
        1481     1489176 :   type(physics_state) :: state1     ! Local copy of state variable
        1482             :
        1483             :   integer :: lchnk                  ! chunk identifier
        1484             :   integer :: ncol                   ! number of atmospheric columns
        1485             :   integer :: istat
        1486             :
        1487             :   integer :: i, k                   ! loop indices
        1488             :
        1489     1489176 :   type(Coords1D) :: p               ! Pressure coordinates
        1490             :
        1491     2978352 :   real(r8) :: ttgw(state%ncol,pver) ! temperature tendency
        1492     2978352 :   real(r8) :: utgw(state%ncol,pver) ! zonal wind tendency
        1493     2978352 :   real(r8) :: vtgw(state%ncol,pver) ! meridional wind tendency
        1494             :
        1495     2978352 :   real(r8) :: ni(state%ncol,pver+1) ! interface Brunt-Vaisalla frequency
        1496     2978352 :   real(r8) :: nm(state%ncol,pver)   ! midpoint Brunt-Vaisalla frequency
        1497     2978352 :   real(r8) :: rhoi(state%ncol,pver+1)     ! interface density
        1498     1489176 :   real(r8), allocatable :: tau(:,:,:)  ! wave Reynolds stress
        1499     2978352 :   real(r8) :: tau0x(state%ncol)     ! c=0 sfc. stress (zonal)
        1500     2978352 :   real(r8) :: tau0y(state%ncol)     ! c=0 sfc. stress (meridional)
        1501     2978352 :   real(r8) :: ubi(state%ncol,pver+1)! projection of wind at interfaces
        1502     2978352 :   real(r8) :: ubm(state%ncol,pver)  ! projection of wind at midpoints
        1503     2978352 :   real(r8) :: xv(state%ncol)        ! unit vector of source wind (x)
        1504     2978352 :   real(r8) :: yv(state%ncol)        ! unit vector of source wind (y)
        1505             :
        1506             :   integer :: m                      ! dummy integers
        1507     2978352 :   real(r8) :: qtgw(state%ncol,pver,pcnst) ! constituents tendencies
        1508             :
        1509             :   ! Reynolds stress for waves propagating in each cardinal direction.
        1510     2978352 :   real(r8) :: taucd(state%ncol,pver+1,4)
        1511             :
        1512             :   ! gravity wave wind tendency for each wave
        1513     1489176 :   real(r8), allocatable :: gwut(:,:,:)
        1514             :
        1515             :   ! Temperature tendencies from diffusion and kinetic energy.
        1516     2978352 :   real(r8) :: dttdf(state%ncol,pver)
        1517     2978352 :   real(r8) :: dttke(state%ncol,pver)
        1518             :
        1519             :   ! Wave phase speeds for each column
        1520     1489176 :   real(r8), allocatable :: phase_speeds(:,:)
        1521             :
        1522             :   ! Efficiency for a gravity wave source.
        1523     2978352 :   real(r8) :: effgw(state%ncol)
        1524             :
        1525             :   ! Coriolis characteristic speed.
        1526     2978352 :   real(r8) :: u_coriolis(state%ncol)
        1527             :
        1528             :   ! Adjustment for inertial gravity waves.
        1529     1489176 :   real(r8), allocatable :: ro_adjust(:,:,:)
        1530             :
        1531             :   ! pbuf fields
        1532             :   ! Molecular diffusivity
        1533     1489176 :   real(r8), pointer :: kvt_in(:,:)
        1534     2978352 :   real(r8) :: kvtt(state%ncol,pver+1)
        1535             :
        1536             :   ! Frontogenesis
        1537     1489176 :   real(r8), pointer :: frontgf(:,:)
        1538     1489176 :   real(r8), pointer :: frontga(:,:)
        1539             :
        1540             :   ! Temperature change due to deep convection.
        1541     1489176 :   real(r8), pointer :: ttend_dp(:,:)
        1542             :   ! Temperature change due to shallow convection.
        1543     1489176 :   real(r8), pointer :: ttend_sh(:,:)
        1544             :
        1545             :   !  New couplings from CLUBB
        1546     1489176 :   real(r8), pointer :: ttend_clubb(:,:)
        1547     1489176 :   real(r8), pointer :: thlp2_clubb_gw(:,:)
        1548     1489176 :   real(r8), pointer :: wpthlp_clubb_gw(:,:)
        1549     1489176 :   real(r8), pointer :: upwp_clubb_gw(:,:)
        1550     1489176 :   real(r8), pointer :: vpwp_clubb_gw(:,:)
        1551     2978352 :   real(r8) :: xpwp_clubb(state%ncol,pver+1)
        1552             :
        1553             :
        1554             :   ! Standard deviation of orography.
        1555     1489176 :   real(r8), pointer :: sgh(:)
        1556             :
        1557             :   ! gridbox area
        1558     1489176 :   real(r8), pointer :: gbxar(:)
        1559             :
        1560             :      ! Beta ridges
        1561             :   ! width of ridges.
        1562     1489176 :   real(r8), pointer :: hwdth(:,:)
        1563             :   ! length of ridges.
        1564     1489176 :   real(r8), pointer :: clngt(:,:)
        1565             :   ! Maximum deviations of ridges.
        1566     1489176 :   real(r8), pointer :: mxdis(:,:)
        1567             :   ! orientation of ridges.
        1568     1489176 :   real(r8), pointer :: angll(:,:)
        1569             :   ! anisotropy of ridges.
        1570     1489176 :   real(r8), pointer :: anixy(:,:)
        1571             :
        1572             :      ! Gamma ridges
        1573             :   ! width of ridges.
        1574     1489176 :   real(r8), pointer :: hwdthg(:,:)
        1575             :   ! length of ridges.
        1576     1489176 :   real(r8), pointer :: clngtg(:,:)
        1577             :   ! Maximum deviations of ridges.
        1578     1489176 :   real(r8), pointer :: mxdisg(:,:)
        1579             :   ! orientation of ridges.
        1580     1489176 :   real(r8), pointer :: angllg(:,:)
        1581             :   ! anisotropy of ridges.
        1582     1489176 :   real(r8), pointer :: anixyg(:,:)
        1583             :
        1584             :   ! Indices of gravity wave source and lowest level where wind tendencies
        1585             :   ! are allowed.
        1586     2978352 :   integer :: src_level(state%ncol)
        1587     2978352 :   integer :: tend_level(state%ncol)
        1588             :
        1589             :   ! Convective source heating depth.
        1590             :   ! heating depth
        1591     2978352 :   real(r8) :: hdepth(state%ncol)
        1592             :   ! maximum heating rate
        1593     2978352 :   real(r8) :: maxq0(state%ncol)
        1594             :
        1595             :   ! Scale sgh to account for landfrac.
        1596     2978352 :   real(r8) :: sgh_scaled(state%ncol)
        1597             :
        1598             :   ! Parameters for the IGW polar taper.
        1599             :   real(r8), parameter :: degree2radian = pi/180._r8
        1600             :   real(r8), parameter :: al0 = 82.5_r8 * degree2radian
        1601             :   real(r8), parameter :: dlat0 = 5.0_r8 * degree2radian
        1602             :
        1603             :   ! effective gw diffusivity at interfaces needed for output
        1604     2978352 :   real(r8) :: egwdffi(state%ncol,pver+1)
        1605             :   ! sum from the two types of spectral GW
        1606     2978352 :   real(r8) :: egwdffi_tot(state%ncol,pver+1)
        1607             :
        1608             :   ! Momentum fluxes used by fixer.
        1609     2978352 :   real(r8) :: um_flux(state%ncol), vm_flux(state%ncol)
        1610             :   ! Energy change used by fixer.
        1611     2978352 :   real(r8) :: de(state%ncol)
        1612             :
        1613             :   ! Which constituents are being affected by diffusion.
        1614             :   logical  :: lq(pcnst)
        1615             :
        1616             :   ! Contiguous copies of state arrays.
        1617     2978352 :   real(r8) :: dse(state%ncol,pver)
        1618     2978352 :   real(r8) :: t(state%ncol,pver)
        1619     2978352 :   real(r8) :: u(state%ncol,pver)
        1620     2978352 :   real(r8) :: v(state%ncol,pver)
        1621     2978352 :   real(r8) :: q(state%ncol,pver,pcnst)
        1622     2978352 :   real(r8) :: piln(state%ncol,pver+1)
        1623     2978352 :   real(r8) :: zm(state%ncol,pver)
        1624     2978352 :   real(r8) :: zi(state%ncol,pver+1)
        1625             :
        1626             :   !------------------------------------------------------------------------
        1627             :
        1628             :   ! Make local copy of input state.
        1629     1489176 :   call physics_state_copy(state, state1)
        1630             :
        1631             :   ! constituents are all treated as wet mmr
        1632     1489176 :   call set_dry_to_wet(state1, convert_cnst_type='dry')
        1633             :
        1634     1489176 :   lchnk = state1%lchnk
        1635     1489176 :   ncol  = state1%ncol
        1636             :
        1637     1489176 :   p = Coords1D(state1%pint(:ncol,:))
        1638             :
        1639   647999352 :   dse = state1%s(:ncol,:)
        1640   647999352 :   t = state1%t(:ncol,:)
        1641   647999352 :   u = state1%u(:ncol,:)
        1642   647999352 :   v = state1%v(:ncol,:)
        1643  1945487232 :   q = state1%q(:ncol,:,:)
        1644   672865128 :   piln = state1%lnpint(:ncol,:)
        1645   647999352 :   zm = state1%zm(:ncol,:)
        1646   672865128 :   zi = state1%zi(:ncol,:)
        1647             :
        1648     5956704 :   lq = .true.
        1649             :   call physics_ptend_init(ptend, state1%psetcols, "Gravity wave drag", &
        1650     1489176 :        ls=.true., lu=.true., lv=.true., lq=lq)
        1651             :
        1652             :   ! Profiles of background state variables
        1653     1489176 :   call gw_prof(ncol, p, cpair, t, rhoi, nm, ni)
        1654             :
        1655     1489176 :   if (do_molec_diff) then
        1656             :      !--------------------------------------------------------
        1657             :      ! Initialize and calculate local molecular diffusivity
        1658             :      !--------------------------------------------------------
        1659             :
        1660           0 :      call pbuf_get_field(pbuf, kvt_idx, kvt_in)  ! kvt_in(1:pcols,1:pver+1)
        1661             :
        1662             :      ! Set kvtt from pbuf field; kvtt still needs a factor of 1/cpairv.
        1663           0 :      kvtt = kvt_in(:ncol,:)
        1664             :
        1665             :      ! Use linear extrapolation of cpairv to top interface.
        1666             :      kvtt(:,1) = kvtt(:,1) / &
        1667           0 :           (1.5_r8*cpairv(:ncol,1,lchnk) - &
        1668           0 :           0.5_r8*cpairv(:ncol,2,lchnk))
        1669             :
        1670             :      ! Interpolate cpairv to other interfaces.
        1671           0 :      do k = 2, nbot_molec
        1672           0 :         kvtt(:,k) = kvtt(:,k) / &
        1673           0 :              (cpairv(:ncol,k+1,lchnk)+cpairv(:ncol,k,lchnk)) * 2._r8
        1674             :      enddo
        1675             :
        1676             :   else
        1677             :
        1678   672865128 :      kvtt = 0._r8
        1679             :
        1680             :   end if
        1681             :
        1682     1489176 :   if (use_gw_front_igw) then
        1683           0 :      u_coriolis = coriolis_speed(band_long, state1%lat(:ncol))
        1684             :   end if
        1685             :
        1686             :   ! Totals that accumulate over different sources.
        1687   672865128 :   egwdffi_tot = 0._r8
        1688     1489176 :   flx_heat = 0._r8
        1689             :
        1690     1489176 :   if (use_gw_movmtn_pbl) then
        1691             :      !------------------------------------------------------------------
        1692             :      !Convective moving mountain gravity waves (Beres scheme).
        1693             :      !------------------------------------------------------------------
        1694             :
        1695           0 :      call outfld('U_MOVMTN_IN', u, ncol, lchnk)
        1696           0 :      call outfld('V_MOVMTN_IN', v, ncol, lchnk)
        1697             :
        1698             :      ! Allocate wavenumber fields.
        1699           0 :      allocate(tau(ncol,-band_movmtn%ngwv:band_movmtn%ngwv,pver+1),stat=istat)
        1700           0 :      call alloc_err(istat,'gw_tend','tau',ncol*(band_movmtn%ngwv**2+1)*(pver+1))
        1701           0 :      allocate(gwut(ncol,pver,-band_movmtn%ngwv:band_movmtn%ngwv),stat=istat)
        1702           0 :      call alloc_err(istat,'gw_tend','gwut',ncol*pver*band_movmtn%ngwv**2+1)
        1703           0 :      allocate(phase_speeds(ncol,-band_movmtn%ngwv:band_movmtn%ngwv),stat=istat)
        1704           0 :      call alloc_err(istat,'gw_tend','phase_speeds',ncol*band_movmtn%ngwv**2+1)
        1705             :
        1706             :      ! Set up heating
        1707           0 :      if (ttend_dp_idx > 0) then
        1708           0 :         call pbuf_get_field(pbuf, ttend_dp_idx, ttend_dp)
        1709             :      else
        1710           0 :         allocate(ttend_dp(pcols,pver), stat=istat)
        1711           0 :         call alloc_err(istat, 'gw_tend', 'ttend_dp', pcols*pver)
        1712           0 :         ttend_dp = 0.0_r8
        1713             :      end if
        1714             :
        1715             :      !   New couplings from CLUBB
        1716           0 :      call pbuf_get_field(pbuf, ttend_clubb_idx, ttend_clubb)
        1717           0 :      call pbuf_get_field(pbuf, thlp2_clubb_gw_idx, thlp2_clubb_gw)
        1718           0 :      call pbuf_get_field(pbuf, wpthlp_clubb_gw_idx, wpthlp_clubb_gw)
        1719           0 :      call pbuf_get_field(pbuf, upwp_clubb_gw_idx, upwp_clubb_gw)
        1720           0 :      call pbuf_get_field(pbuf, vpwp_clubb_gw_idx, vpwp_clubb_gw)
        1721             :
        1722           0 :      xpwp_clubb(:ncol,:) = sqrt( upwp_clubb_gw(:ncol,:)**2 + vpwp_clubb_gw(:ncol,:)**2 )
        1723             :
        1724           0 :      effgw = 1._r8
        1725             :      call gw_movmtn_src(ncol, lchnk, band_movmtn , movmtn_desc, &
        1726           0 :           u, v, ttend_dp(:ncol,:), ttend_clubb(:ncol,:), xpwp_clubb(:ncol,:)  , &
        1727             :           zm, alpha_gw_movmtn, src_level, tend_level, &
        1728             :           tau, ubm, ubi, xv, yv, &
        1729           0 :           phase_speeds, hdepth)
        1730             :      !-------------------------------------------------------------
        1731             :      ! gw_movmtn_src returns wave-relative wind profiles ubm,ubi
        1732             :      ! and unit vector components describing direction of wavevector
        1733             :      ! and application of wave-drag force. I believe correct setting
        1734             :      ! for c is c=0, since it is incorporated in ubm and (xv,yv)
        1735             :      !--------------------------------------------------------------
        1736             :
        1737           0 :      call outfld('SRC_LEVEL_MOVMTN', real(src_level,r8), ncol, lchnk)
        1738           0 :      call outfld('TND_LEVEL_MOVMTN', real(tend_level,r8), ncol, lchnk)
        1739           0 :      call outfld('UBI_MOVMTN', ubi, ncol, lchnk)
        1740           0 :      call outfld('UBM_MOVMTN', ubm, ncol, lchnk)
        1741             :
        1742             :      call gw_drag_prof(ncol, band_movmtn, p, src_level, tend_level, dt, &
        1743             :           t, vramp,    &
        1744             :           piln, rhoi,       nm,   ni, ubm,  ubi,  xv,    yv,   &
        1745             :           effgw,   phase_speeds,       kvtt, q,  dse,  tau,  utgw,  vtgw, &
        1746             :           ttgw, qtgw, egwdffi,  gwut, dttdf, dttke,            &
        1747           0 :           lapply_effgw_in=gw_apply_tndmax )
        1748             :
        1749             :      ! Project stress into directional components.
        1750           0 :      taucd = calc_taucd(ncol, band_movmtn%ngwv, tend_level, tau, phase_speeds, xv, yv, ubi)
        1751             :
        1752             :      !  add the diffusion coefficients
        1753           0 :      do k = 1, pver+1
        1754           0 :         egwdffi_tot(:,k) = egwdffi_tot(:,k) + egwdffi(:,k)
        1755             :      end do
        1756             :
        1757             :      ! Store constituents tendencies
        1758           0 :      do m=1, pcnst
        1759           0 :         do k = 1, pver
        1760           0 :            ptend%q(:ncol,k,m) = ptend%q(:ncol,k,m) + qtgw(:,k,m)
        1761             :         end do
        1762             :      end do
        1763             :
        1764             :       ! Add the momentum tendencies to the output tendency arrays.
        1765           0 :      do k = 1, pver
        1766           0 :         ptend%u(:ncol,k) = ptend%u(:ncol,k) + utgw(:,k)
        1767           0 :         ptend%v(:ncol,k) = ptend%v(:ncol,k) + vtgw(:,k)
        1768             :      end do
        1769             :
        1770           0 :      do k = 1, pver
        1771           0 :         ptend%s(:ncol,k) = ptend%s(:ncol,k) + ttgw(:,k)
        1772             :      end do
        1773             :
        1774           0 :      call outfld('TAU_MOVMTN', tau(:,0,:), ncol, lchnk)
        1775           0 :      call outfld('GWUT_MOVMTN', gwut(:,:,0), ncol, lchnk)
        1776           0 :      call outfld('VTGW_MOVMTN', vtgw, ncol, lchnk)
        1777           0 :      call outfld('UTGW_MOVMTN', utgw, ncol, lchnk)
        1778           0 :      call outfld('HDEPTH_MOVMTN', hdepth/1000._r8, ncol, lchnk)
        1779           0 :      call outfld('NETDT_MOVMTN', ttend_dp, pcols, lchnk)
        1780           0 :      call outfld('TTEND_CLUBB', ttend_clubb, pcols, lchnk)
        1781           0 :      call outfld('THLP2_CLUBB_GW', thlp2_clubb_gw, pcols, lchnk)
        1782           0 :      call outfld('WPTHLP_CLUBB_GW', wpthlp_clubb_gw, pcols, lchnk)
        1783           0 :      call outfld('UPWP_CLUBB_GW', upwp_clubb_gw, pcols, lchnk)
        1784           0 :      call outfld('VPWP_CLUBB_GW', vpwp_clubb_gw, pcols, lchnk)
        1785             :
        1786             :      !Deallocate variables that are no longer used:
        1787           0 :      deallocate(tau, gwut, phase_speeds)
        1788             :
        1789             :      !Deallocate/nullify ttend_dp if not a pbuf variable:
        1790           0 :      if (ttend_dp_idx <= 0) then
        1791           0 :        deallocate(ttend_dp)
        1792             :        nullify(ttend_dp)
        1793             :      end if
        1794             :
        1795             :   end if
        1796             :
        1797     1489176 :   if (use_gw_convect_dp) then
        1798             :      !------------------------------------------------------------------
        1799             :      ! Convective gravity waves (Beres scheme, deep).
        1800             :      !------------------------------------------------------------------
        1801             :
        1802             :      ! Allocate wavenumber fields.
        1803           0 :      allocate(tau(ncol,-band_mid%ngwv:band_mid%ngwv,pver+1),stat=istat)
        1804           0 :      call alloc_err(istat,'gw_tend','tau',ncol*(band_mid%ngwv**2+1)*(pver+1))
        1805           0 :      allocate(gwut(ncol,pver,-band_mid%ngwv:band_mid%ngwv),stat=istat)
        1806           0 :      call alloc_err(istat,'gw_tend','gwut',ncol*pver*(band_mid%ngwv**2+1))
        1807           0 :      allocate(phase_speeds(ncol,-band_mid%ngwv:band_mid%ngwv),stat=istat)
        1808           0 :      call alloc_err(istat,'gw_tend','tau',ncol*(band_mid%ngwv**2+1))
        1809             :
        1810             :      ! Set up heating
        1811           0 :      call pbuf_get_field(pbuf, ttend_dp_idx, ttend_dp)
        1812             :
        1813             :      ! Efficiency of gravity wave momentum transfer.
        1814             :      ! This is really only to remove the pole points.
        1815           0 :      where (pi/2._r8 - abs(state1%lat(:ncol)) >= 4*epsilon(1._r8))
        1816             :         effgw = effgw_beres_dp
        1817             :      elsewhere
        1818             :         effgw = 0._r8
        1819             :      end where
        1820             :
        1821             :      ! Determine wave sources for Beres deep scheme
        1822             :      call gw_beres_src(ncol, band_mid, beres_dp_desc, &
        1823           0 :           u, v, ttend_dp(:ncol,:), zm, src_level, tend_level, tau, &
        1824           0 :           ubm, ubi, xv, yv, phase_speeds, hdepth, maxq0)
        1825             :
        1826             :      ! Solve for the drag profile with Beres source spectrum.
        1827             :      call gw_drag_prof(ncol, band_mid, p, src_level, tend_level, dt, &
        1828             :           t, vramp,    &
        1829             :           piln, rhoi,       nm,   ni, ubm,  ubi,  xv,    yv,   &
        1830             :           effgw,   phase_speeds,       kvtt, q,  dse,  tau,  utgw,  vtgw, &
        1831             :           ttgw, qtgw, egwdffi,  gwut, dttdf, dttke,            &
        1832           0 :           lapply_effgw_in=gw_apply_tndmax)
        1833             :
        1834             :      ! Project stress into directional components.
        1835           0 :      taucd = calc_taucd(ncol, band_mid%ngwv, tend_level, tau, phase_speeds, xv, yv, ubi)
        1836             :
        1837             :      !  add the diffusion coefficients
        1838           0 :      do k = 1, pver+1
        1839           0 :         egwdffi_tot(:,k) = egwdffi_tot(:,k) + egwdffi(:,k)
        1840             :      end do
        1841             :
        1842             :      ! Store constituents tendencies
        1843           0 :      do m=1, pcnst
        1844           0 :         do k = 1, pver
        1845           0 :            ptend%q(:ncol,k,m) = ptend%q(:ncol,k,m) + qtgw(:,k,m)
        1846             :         end do
        1847             :      end do
        1848             :
        1849             :      ! Find momentum flux, and use it to fix the wind tendencies below
        1850             :      ! the gravity wave region.
        1851           0 :      call momentum_flux(tend_level, taucd, um_flux, vm_flux)
        1852           0 :      call momentum_fixer(tend_level, p, um_flux, vm_flux, utgw, vtgw)
        1853             :
        1854             :      ! Add the momentum tendencies to the output tendency arrays.
        1855           0 :      do k = 1, pver
        1856           0 :         ptend%u(:ncol,k) = ptend%u(:ncol,k) + utgw(:,k)
        1857           0 :         ptend%v(:ncol,k) = ptend%v(:ncol,k) + vtgw(:,k)
        1858             :      end do
        1859             :
        1860             :      ! Find energy change in the current state, and use fixer to apply
        1861             :      ! the difference in lower levels.
        1862           0 :      call energy_change(dt, p, u, v, ptend%u(:ncol,:), &
        1863           0 :           ptend%v(:ncol,:), ptend%s(:ncol,:)+ttgw, de)
        1864           0 :      call energy_fixer(tend_level, p, de-flx_heat(:ncol), ttgw)
        1865             :
        1866           0 :      do k = 1, pver
        1867           0 :         ptend%s(:ncol,k) = ptend%s(:ncol,k) + ttgw(:,k)
        1868             :      end do
        1869             :
        1870             :      ! Change ttgw to a temperature tendency before outputing it.
        1871           0 :      ttgw = ttgw / cpair
        1872             :      call gw_spec_outflds(beres_dp_pf, lchnk, ncol, band_mid, phase_speeds, u, v, &
        1873           0 :           xv, yv, gwut, dttdf, dttke, tau(:,:,2:), utgw, vtgw, ttgw, &
        1874           0 :           taucd)
        1875             :
        1876             :      ! Diagnostic outputs (convert hdepth to km).
        1877           0 :      call outfld('NETDT', ttend_dp, pcols, lchnk)
        1878           0 :      call outfld('HDEPTH', hdepth/1000._r8, ncol, lchnk)
        1879           0 :      call outfld('MAXQ0', maxq0, ncol, lchnk)
        1880             :
        1881           0 :      deallocate(tau, gwut, phase_speeds)
        1882             :
        1883             :   end if
        1884             :
        1885     1489176 :   if (use_gw_convect_sh) then
        1886             :      !------------------------------------------------------------------
        1887             :      ! Convective gravity waves (Beres scheme, shallow).
        1888             :      !------------------------------------------------------------------
        1889             :
        1890             :      ! Allocate wavenumber fields.
        1891           0 :      allocate(tau(ncol,-band_mid%ngwv:band_mid%ngwv,pver+1),stat=istat)
        1892           0 :      call alloc_err(istat,'gw_tend','tau',ncol*(band_mid%ngwv**2+1)*(pver+1))
        1893           0 :      allocate(gwut(ncol,pver,-band_mid%ngwv:band_mid%ngwv),stat=istat)
        1894           0 :      call alloc_err(istat,'gw_tend','gwut',ncol*pver*(band_mid%ngwv**2+1))
        1895           0 :      allocate(phase_speeds(ncol,-band_mid%ngwv:band_mid%ngwv),stat=istat)
        1896           0 :      call alloc_err(istat,'gw_tend','phase_speeds',ncol*(band_mid%ngwv**2+1))
        1897             :
        1898             :      ! Set up heating
        1899           0 :      call pbuf_get_field(pbuf, ttend_sh_idx, ttend_sh)
        1900             :
        1901             :      ! Efficiency of gravity wave momentum transfer.
        1902             :      ! This is really only to remove the pole points.
        1903           0 :      where (pi/2._r8 - abs(state1%lat(:ncol)) >= 4*epsilon(1._r8))
        1904             :         effgw = effgw_beres_sh
        1905             :      elsewhere
        1906             :         effgw = 0._r8
        1907             :      end where
        1908             :
        1909             :      ! Determine wave sources for Beres shallow scheme
        1910             :      call gw_beres_src(ncol, band_mid, beres_sh_desc, &
        1911           0 :           u, v, ttend_sh(:ncol,:), zm, src_level, tend_level, tau, &
        1912           0 :           ubm, ubi, xv, yv, phase_speeds, hdepth, maxq0)
        1913             :
        1914             :      ! Solve for the drag profile with Beres source spectrum.
        1915             :      call gw_drag_prof(ncol, band_mid, p, src_level, tend_level,  dt, &
        1916             :           t, vramp,    &
        1917             :           piln, rhoi,       nm,   ni, ubm,  ubi,  xv,    yv,   &
        1918             :           effgw,   phase_speeds,       kvtt, q,  dse,  tau,  utgw,  vtgw, &
        1919             :           ttgw, qtgw, egwdffi,  gwut, dttdf, dttke,            &
        1920           0 :           lapply_effgw_in=gw_apply_tndmax)
        1921             :
        1922             :      ! Project stress into directional components.
        1923           0 :      taucd = calc_taucd(ncol, band_mid%ngwv, tend_level, tau, phase_speeds, xv, yv, ubi)
        1924             :
        1925             :      !  add the diffusion coefficients
        1926           0 :      do k = 1, pver+1
        1927           0 :         egwdffi_tot(:,k) = egwdffi_tot(:,k) + egwdffi(:,k)
        1928             :      end do
        1929             :
        1930             :      ! Store constituents tendencies
        1931           0 :      do m=1, pcnst
        1932           0 :         do k = 1, pver
        1933           0 :            ptend%q(:ncol,k,m) = ptend%q(:ncol,k,m) + qtgw(:,k,m)
        1934             :         end do
        1935             :      end do
        1936             :
        1937             :      ! Add the momentum tendencies to the output tendency arrays.
        1938             :      ! Don't calculate fixers, since we are too close to the ground to
        1939             :      ! spread momentum/energy differences across low layers.
        1940           0 :      do k = 1, pver
        1941           0 :         ptend%u(:ncol,k) = ptend%u(:ncol,k) + utgw(:,k)
        1942           0 :         ptend%v(:ncol,k) = ptend%v(:ncol,k) + vtgw(:,k)
        1943           0 :         ptend%s(:ncol,k) = ptend%s(:ncol,k) + ttgw(:,k)
        1944             :      end do
        1945             :
        1946             :      ! Calculate energy change for output to CAM's energy checker.
        1947             :      ! This is sort of cheating; we don't have a good a priori idea of the
        1948             :      ! energy coming from surface stress, so we just integrate what we and
        1949             :      ! actually have so far and overwrite flx_heat with that.
        1950           0 :      call energy_change(dt, p, u, v, ptend%u(:ncol,:), &
        1951           0 :           ptend%v(:ncol,:), ptend%s(:ncol,:), de)
        1952           0 :      flx_heat(:ncol) = de
        1953             :
        1954             :      ! Change ttgw to a temperature tendency before outputing it.
        1955           0 :      ttgw = ttgw / cpair
        1956             :      call gw_spec_outflds(beres_sh_pf, lchnk, ncol, band_mid, phase_speeds, u, v, &
        1957           0 :           xv, yv, gwut, dttdf, dttke, tau(:,:,2:), utgw, vtgw, ttgw, &
        1958           0 :           taucd)
        1959             :
        1960             :      ! Diagnostic outputs (convert SHDEPTH to km).
        1961           0 :      call outfld ('SNETDT', ttend_sh, pcols, lchnk)
        1962           0 :      call outfld ('SHDEPTH', hdepth/1000._r8, ncol, lchnk)
        1963           0 :      call outfld ('SMAXQ0', maxq0, ncol, lchnk)
        1964             :
        1965           0 :      deallocate(tau, gwut, phase_speeds)
        1966             :
        1967             :   end if
        1968             :
        1969     1489176 :   if (use_gw_front .or. use_gw_front_igw) then
        1970             :      ! Get frontogenesis physics buffer fields set by dynamics.
        1971           0 :      call pbuf_get_field(pbuf, frontgf_idx, frontgf)
        1972           0 :      call pbuf_get_field(pbuf, frontga_idx, frontga)
        1973             :
        1974             :      ! Output for diagnostics.
        1975           0 :      call outfld ('FRONTGF', frontgf, pcols, lchnk)
        1976           0 :      call outfld ('FRONTGFA', frontga, pcols, lchnk)
        1977             :   end if
        1978             :
        1979     1489176 :   if (use_gw_front) then
        1980             :      !------------------------------------------------------------------
        1981             :      ! Frontally generated gravity waves
        1982             :      !------------------------------------------------------------------
        1983             :
        1984             :      ! Allocate wavenumber fields.
        1985           0 :      allocate(tau(ncol,-band_mid%ngwv:band_mid%ngwv,pver+1),stat=istat)
        1986           0 :      call alloc_err(istat,'gw_tend','tau',ncol*(band_mid%ngwv**2+1)*(pver+1))
        1987           0 :      allocate(gwut(ncol,pver,-band_mid%ngwv:band_mid%ngwv),stat=istat)
        1988           0 :      call alloc_err(istat,'gw_tend','gwut',ncol*pver*(band_mid%ngwv**2+1))
        1989           0 :      allocate(phase_speeds(ncol,-band_mid%ngwv:band_mid%ngwv),stat=istat)
        1990           0 :      call alloc_err(istat,'gw_tend','tau',ncol*(band_mid%ngwv**2+1))
        1991             :
        1992             :      ! Efficiency of gravity wave momentum transfer.
        1993           0 :      effgw = effgw_cm
        1994             :      ! Frontogenesis is too high at the poles (at least for the FV
        1995             :      ! dycore), so introduce a polar taper.
        1996           0 :      if (gw_polar_taper) effgw = effgw * cos(state1%lat(:ncol))
        1997             :
        1998             :      ! Determine the wave source for C&M background spectrum
        1999           0 :      call gw_cm_src(ncol, band_mid, cm_desc, u, v, frontgf(:ncol,:), &
        2000           0 :           src_level, tend_level, tau, ubm, ubi, xv, yv, phase_speeds)
        2001             :
        2002             :      ! Solve for the drag profile with C&M source spectrum.
        2003             :      call gw_drag_prof(ncol, band_mid, p, src_level, tend_level,  dt, &
        2004             :           t, vramp,   &
        2005             :           piln, rhoi,       nm,   ni, ubm,  ubi,  xv,    yv,   &
        2006             :           effgw,   phase_speeds,       kvtt, q,  dse,  tau,  utgw,  vtgw, &
        2007             :           ttgw, qtgw, egwdffi,  gwut, dttdf, dttke,            &
        2008           0 :           lapply_effgw_in=gw_apply_tndmax)
        2009             :
        2010             :      ! Project stress into directional components.
        2011           0 :      taucd = calc_taucd(ncol, band_mid%ngwv, tend_level, tau, phase_speeds, xv, yv, ubi)
        2012             :
        2013             :      !  add the diffusion coefficients
        2014           0 :      do k = 1, pver+1
        2015           0 :         egwdffi_tot(:,k) = egwdffi_tot(:,k) + egwdffi(:,k)
        2016             :      end do
        2017             :
        2018             :      !Add the constituent tendencies
        2019           0 :      do m=1, pcnst
        2020           0 :         do k = 1, pver
        2021           0 :            ptend%q(:ncol,k,m) = ptend%q(:ncol,k,m) + qtgw(:,k,m)
        2022             :         end do
        2023             :      end do
        2024             :
        2025             :      ! Find momentum flux, and use it to fix the wind tendencies below
        2026             :      ! the gravity wave region.
        2027           0 :      call momentum_flux(tend_level, taucd, um_flux, vm_flux)
        2028           0 :      call momentum_fixer(tend_level, p, um_flux, vm_flux, utgw, vtgw)
        2029             :
        2030             :      ! add the momentum tendencies to the output tendency arrays
        2031           0 :      do k = 1, pver
        2032           0 :         ptend%u(:ncol,k) = ptend%u(:ncol,k) + utgw(:,k)
        2033           0 :         ptend%v(:ncol,k) = ptend%v(:ncol,k) + vtgw(:,k)
        2034             :      end do
        2035             :
        2036             :      ! Find energy change in the current state, and use fixer to apply
        2037             :      ! the difference in lower levels.
        2038           0 :      call energy_change(dt, p, u, v, ptend%u(:ncol,:), &
        2039           0 :           ptend%v(:ncol,:), ptend%s(:ncol,:)+ttgw, de)
        2040           0 :      call energy_fixer(tend_level, p, de-flx_heat(:ncol), ttgw)
        2041             :
        2042           0 :      do k = 1, pver
        2043           0 :         ptend%s(:ncol,k) = ptend%s(:ncol,k) + ttgw(:,k)
        2044             :      end do
        2045             :
        2046             :      ! Change ttgw to a temperature tendency before outputing it.
        2047           0 :      ttgw = ttgw / cpair
        2048             :      call gw_spec_outflds(cm_pf, lchnk, ncol, band_mid, phase_speeds, u, v, &
        2049           0 :           xv, yv, gwut, dttdf, dttke, tau(:,:,2:), utgw, vtgw, ttgw, &
        2050           0 :           taucd)
        2051             :
        2052           0 :      deallocate(tau, gwut, phase_speeds)
        2053             :
        2054             :   end if
        2055             :
        2056     1489176 :   if (use_gw_front_igw) then
        2057             :      !------------------------------------------------------------------
        2058             :      ! Frontally generated inertial gravity waves
        2059             :      !------------------------------------------------------------------
        2060             :
        2061             :      ! Allocate wavenumber fields.
        2062           0 :      allocate(tau(ncol,-band_long%ngwv:band_long%ngwv,pver+1),stat=istat)
        2063           0 :      call alloc_err(istat,'gw_tend','tau',ncol*(band_long%ngwv**2+1)*(pver+1))
        2064           0 :      allocate(gwut(ncol,pver,-band_long%ngwv:band_long%ngwv),stat=istat)
        2065           0 :      call alloc_err(istat,'gw_tend','gwut',ncol*pver*(band_long%ngwv**2+1))
        2066           0 :      allocate(phase_speeds(ncol,-band_long%ngwv:band_long%ngwv),stat=istat)
        2067           0 :      call alloc_err(istat,'gw_tend','phase_speeds',ncol*(band_long%ngwv**2+1))
        2068           0 :      allocate(ro_adjust(ncol,-band_long%ngwv:band_long%ngwv,pver+1),stat=istat)
        2069           0 :      call alloc_err(istat,'gw_tend','ro_adjust',ncol*(band_long%ngwv**2+1)*(pver+1))
        2070             :
        2071             :      ! Efficiency of gravity wave momentum transfer.
        2072           0 :      effgw = effgw_cm_igw
        2073             :
        2074             :      ! Frontogenesis is too high at the poles (at least for the FV
        2075             :      ! dycore), so introduce a polar taper.
        2076           0 :      if (gw_polar_taper) then
        2077           0 :         where (abs(state1%lat(:ncol)) <= 89._r8*degree2radian)
        2078             :            effgw = effgw * 0.25_r8 * &
        2079           0 :                  (1._r8+tanh((state1%lat(:ncol)+al0)/dlat0)) * &
        2080             :                  (1._r8-tanh((state1%lat(:ncol)-al0)/dlat0))
        2081             :         elsewhere
        2082             :            effgw = 0._r8
        2083             :         end where
        2084             :      end if
        2085             :
        2086             :      ! Determine the wave source for C&M background spectrum
        2087           0 :      call gw_cm_src(ncol, band_long, cm_igw_desc, u, v, frontgf(:ncol,:), &
        2088           0 :           src_level, tend_level, tau, ubm, ubi, xv, yv, phase_speeds)
        2089             :
        2090             :      call adjust_inertial(band_long, tend_level, u_coriolis, phase_speeds, ubi, &
        2091           0 :           tau, ro_adjust)
        2092             :
        2093             :      ! Solve for the drag profile with C&M source spectrum.
        2094             :      call gw_drag_prof(ncol, band_long, p, src_level, tend_level,  dt, &
        2095             :           t, vramp,    &
        2096             :           piln, rhoi,       nm,   ni, ubm,  ubi,  xv,    yv,   &
        2097             :           effgw,   phase_speeds,       kvtt, q,  dse,  tau,  utgw,  vtgw, &
        2098             :           ttgw, qtgw, egwdffi,  gwut, dttdf, dttke, ro_adjust=ro_adjust, &
        2099           0 :           lapply_effgw_in=gw_apply_tndmax)
        2100             :
        2101             :      ! Project stress into directional components.
        2102           0 :      taucd = calc_taucd(ncol, band_long%ngwv, tend_level, tau, phase_speeds, xv, yv, ubi)
        2103             :
        2104             :      !  add the diffusion coefficients
        2105           0 :      do k = 1, pver+1
        2106           0 :         egwdffi_tot(:,k) = egwdffi_tot(:,k) + egwdffi(:,k)
        2107             :      end do
        2108             :
        2109             :      !Add the constituent tendencies
        2110           0 :      do m=1, pcnst
        2111           0 :         do k = 1, pver
        2112           0 :            ptend%q(:ncol,k,m) = ptend%q(:ncol,k,m) + qtgw(:,k,m)
        2113             :         end do
        2114             :      end do
        2115             :
        2116             :      ! Find momentum flux, and use it to fix the wind tendencies below
        2117             :      ! the gravity wave region.
        2118           0 :      call momentum_flux(tend_level, taucd, um_flux, vm_flux)
        2119           0 :      call momentum_fixer(tend_level, p, um_flux, vm_flux, utgw, vtgw)
        2120             :
        2121             :      ! add the momentum tendencies to the output tendency arrays
        2122           0 :      do k = 1, pver
        2123           0 :         ptend%u(:ncol,k) = ptend%u(:ncol,k) + utgw(:,k)
        2124           0 :         ptend%v(:ncol,k) = ptend%v(:ncol,k) + vtgw(:,k)
        2125             :      end do
        2126             :
        2127             :      ! Find energy change in the current state, and use fixer to apply
        2128             :      ! the difference in lower levels.
        2129           0 :      call energy_change(dt, p, u, v, ptend%u(:ncol,:), &
        2130           0 :           ptend%v(:ncol,:), ptend%s(:ncol,:)+ttgw, de)
        2131           0 :      call energy_fixer(tend_level, p, de-flx_heat(:ncol), ttgw)
        2132             :
        2133           0 :      do k = 1, pver
        2134           0 :         ptend%s(:ncol,k) = ptend%s(:ncol,k) + ttgw(:,k)
        2135             :      end do
        2136             :
        2137             :      ! Change ttgw to a temperature tendency before outputing it.
        2138           0 :      ttgw = ttgw / cpair
        2139             :      call gw_spec_outflds(cm_igw_pf, lchnk, ncol, band_long, phase_speeds, u, v, &
        2140           0 :           xv, yv, gwut, dttdf, dttke, tau(:,:,2:), utgw, vtgw, ttgw, &
        2141           0 :           taucd)
        2142             :
        2143           0 :      deallocate(tau, gwut, phase_speeds, ro_adjust)
        2144             :
        2145             :   end if
        2146             :
        2147     1489176 :   if (use_gw_oro) then
        2148             :      !---------------------------------------------------------------------
        2149             :      ! Orographic stationary gravity waves
        2150             :      !---------------------------------------------------------------------
        2151             :
        2152             :      ! Allocate wavenumber fields.
        2153     5956704 :      allocate(tau(ncol,band_oro%ngwv:band_oro%ngwv,pver+1),stat=istat)
        2154     1489176 :      call alloc_err(istat,'gw_tend','tau',ncol*(band_oro%ngwv**2+1)*(pver+1))
        2155     4467528 :      allocate(gwut(ncol,pver,band_oro%ngwv:band_oro%ngwv),stat=istat)
        2156     1489176 :      call alloc_err(istat,'gw_tend','gwut',ncol*pver*(band_oro%ngwv**2+1))
        2157     4467528 :      allocate(phase_speeds(ncol,band_oro%ngwv:band_oro%ngwv),stat=istat)
        2158     1489176 :      call alloc_err(istat,'gw_tend','phase_speeds',ncol*(band_oro%ngwv**2+1))
        2159             :
        2160             :      ! Efficiency of gravity wave momentum transfer.
        2161             :      ! Take into account that wave sources are only over land.
        2162     1489176 :      call pbuf_get_field(pbuf, sgh_idx, sgh)
        2163             :
        2164     1489176 :      if (gw_lndscl_sgh) then
        2165   121350528 :         where (cam_in%landfrac(:ncol) >= epsilon(1._r8))
        2166     1489176 :            effgw = effgw_oro * cam_in%landfrac(:ncol)
        2167     1489176 :            sgh_scaled = sgh(:ncol) / sqrt(cam_in%landfrac(:ncol))
        2168             :         elsewhere
        2169             :            effgw = 0._r8
        2170             :            sgh_scaled = 0._r8
        2171             :         end where
        2172             :
        2173             :         ! Determine the orographic wave source
        2174             :         call gw_oro_src(ncol, band_oro, p, &
        2175             :            u, v, t, sgh_scaled, zm, nm, &
        2176     1489176 :            src_level, tend_level, tau, ubm, ubi, xv, yv, phase_speeds)
        2177             :      else
        2178           0 :         effgw = effgw_oro
        2179             :
        2180             :         ! Determine the orographic wave source
        2181             :         call gw_oro_src(ncol, band_oro, p, &
        2182           0 :              u, v, t, sgh(:ncol), zm, nm, &
        2183           0 :              src_level, tend_level, tau, ubm, ubi, xv, yv, phase_speeds)
        2184             :      endif
        2185    24865776 :      do i = 1, ncol
        2186    24865776 :         if (state1%lat(i) < 0._r8) then
        2187   642856500 :            tau(i,:,:) = tau(i,:,:) * gw_oro_south_fac
        2188             :         end if
        2189             :      end do
        2190             :
        2191             :      ! Solve for the drag profile with orographic sources.
        2192             :      call gw_drag_prof(ncol, band_oro, p, src_level, tend_level,   dt,   &
        2193             :           t, vramp,   &
        2194             :           piln, rhoi,       nm,   ni, ubm,  ubi,  xv,    yv,   &
        2195             :           effgw, phase_speeds, kvtt, q,  dse,  tau,  utgw,  vtgw, &
        2196             :           ttgw, qtgw, egwdffi,  gwut, dttdf, dttke,            &
        2197     1489176 :           lapply_effgw_in=gw_apply_tndmax)
        2198             :
        2199             :      ! For orographic waves, don't bother with taucd, since there are no
        2200             :      ! momentum conservation routines or directional diagnostics.
        2201             :
        2202             :      !  add the diffusion coefficients
        2203    41696928 :      do k = 1, pver+1
        2204   672865128 :         egwdffi_tot(:,k) = egwdffi_tot(:,k) + egwdffi(:,k)
        2205             :      end do
        2206             :
        2207             :      ! Add the orographic tendencies to the spectrum tendencies.
        2208             :      ! Don't calculate fixers, since we are too close to the ground to
        2209             :      ! spread momentum/energy differences across low layers.
        2210    40207752 :      do k = 1, pver
        2211   646510176 :         ptend%u(:ncol,k) = ptend%u(:ncol,k) + utgw(:,k)
        2212   646510176 :         ptend%v(:ncol,k) = ptend%v(:ncol,k) + vtgw(:,k)
        2213   646510176 :         ptend%s(:ncol,k) = ptend%s(:ncol,k) + ttgw(:,k)
        2214             :         ! Convert to temperature tendency for output.
        2215   647999352 :         ttgw(:,k) = ttgw(:,k) / cpairv(:ncol, k, lchnk)
        2216             :      end do
        2217             :
        2218             :      ! Calculate energy change for output to CAM's energy checker.
        2219             :      ! This is sort of cheating; we don't have a good a priori idea of the
        2220             :      ! energy coming from surface stress, so we just integrate what we and
        2221             :      ! actually have so far and overwrite flx_heat with that.
        2222           0 :      call energy_change(dt, p, u, v, ptend%u(:ncol,:), &
        2223     1489176 :           ptend%v(:ncol,:), ptend%s(:ncol,:), de)
        2224    24865776 :      flx_heat(:ncol) = de
        2225             :
        2226     5956704 :      do m = 1, pcnst
        2227   122112432 :         do k = 1, pver
        2228  1943998056 :            ptend%q(:ncol,k,m) = ptend%q(:ncol,k,m) + qtgw(:,k,m)
        2229             :         end do
        2230             :      end do
        2231             :
        2232             :      ! Write output fields to history file
        2233     1489176 :      call outfld('TAUAORO', tau(:,0,:),  ncol, lchnk)
        2234     1489176 :      call outfld('UTGWORO', utgw,  ncol, lchnk)
        2235     1489176 :      call outfld('VTGWORO', vtgw,  ncol, lchnk)
        2236     1489176 :      call outfld('TTGWORO', ttgw,  ncol, lchnk)
        2237   647999352 :      call outfld('TTGWSDFORO', dttdf / cpair,  ncol, lchnk)
        2238   647999352 :      call outfld('TTGWSKEORO', dttke / cpair,  ncol, lchnk)
        2239    24865776 :      tau0x = tau(:,0,pver+1) * xv
        2240    24865776 :      tau0y = tau(:,0,pver+1) * yv
        2241     1489176 :      call outfld('TAUGWX', tau0x, ncol, lchnk)
        2242     1489176 :      call outfld('TAUGWY', tau0y, ncol, lchnk)
        2243             :
        2244     1489176 :      deallocate(tau, gwut, phase_speeds)
        2245             :
        2246             :   end if
        2247             :
        2248     1489176 :   if (use_gw_rdg_beta) then
        2249             :      !---------------------------------------------------------------------
        2250             :      ! Orographic stationary gravity waves
        2251             :      !---------------------------------------------------------------------
        2252             :
        2253             :      ! Efficiency of gravity wave momentum transfer.
        2254             :      ! Take into account that wave sources are only over land.
        2255             :
        2256           0 :      hwdth => rdg_hwdth(:ncol,:,lchnk)
        2257           0 :      clngt => rdg_clngt(:ncol,:,lchnk)
        2258           0 :      gbxar => rdg_gbxar(:ncol,lchnk)
        2259           0 :      mxdis => rdg_mxdis(:ncol,:,lchnk)
        2260           0 :      angll => rdg_angll(:ncol,:,lchnk)
        2261           0 :      anixy => rdg_anixy(:ncol,:,lchnk)
        2262             :
        2263           0 :      where(mxdis < 0._r8)
        2264             :         mxdis = 0._r8
        2265             :      end where
        2266             :
        2267             :      ! Save state at top of routine
        2268             :      ! Useful for unit testing checks
        2269           0 :      call outfld('UEGW', u ,  ncol, lchnk)
        2270           0 :      call outfld('VEGW', v ,  ncol, lchnk)
        2271           0 :      call outfld('TEGW', t ,  ncol, lchnk)
        2272           0 :      call outfld('ZEGW', zi , ncol, lchnk)
        2273           0 :      call outfld('ZMGW', zm , ncol, lchnk)
        2274             :
        2275             :      call gw_rdg_calc(&
        2276             :         'BETA ', ncol, lchnk, n_rdg_beta, dt,     &
        2277             :         u, v, t, p, piln, zm, zi,                 &
        2278             :         nm, ni, rhoi, kvtt, q, dse,               &
        2279             :         effgw_rdg_beta, effgw_rdg_beta_max,       &
        2280             :         hwdth, clngt, gbxar, mxdis, angll, anixy, &
        2281             :         rdg_beta_cd_llb, trpd_leewv_rdg_beta,     &
        2282           0 :         ptend, flx_heat)
        2283             :
        2284             :   end if
        2285             :
        2286     1489176 :   if (use_gw_rdg_gamma) then
        2287             :      !---------------------------------------------------------------------
        2288             :      ! Orographic stationary gravity waves
        2289             :      !---------------------------------------------------------------------
        2290             :
        2291             :      ! Efficiency of gravity wave momentum transfer.
        2292             :      ! Take into account that wave sources are only over land.
        2293             :
        2294           0 :      hwdthg => rdg_hwdthg(:ncol,:,lchnk)
        2295           0 :      clngtg => rdg_clngtg(:ncol,:,lchnk)
        2296           0 :      gbxar  => rdg_gbxar(:ncol,lchnk)
        2297           0 :      mxdisg => rdg_mxdisg(:ncol,:,lchnk)
        2298           0 :      angllg => rdg_angllg(:ncol,:,lchnk)
        2299           0 :      anixyg => rdg_anixyg(:ncol,:,lchnk)
        2300             :
        2301           0 :      where(mxdisg < 0._r8)
        2302             :         mxdisg = 0._r8
        2303             :      end where
        2304             :
        2305             :      call gw_rdg_calc(&
        2306             :         'GAMMA', ncol, lchnk, n_rdg_gamma, dt,         &
        2307             :         u, v, t, p, piln, zm, zi,                      &
        2308             :         nm, ni, rhoi, kvtt, q, dse,                    &
        2309             :         effgw_rdg_gamma, effgw_rdg_gamma_max,          &
        2310             :         hwdthg, clngtg, gbxar, mxdisg, angllg, anixyg, &
        2311             :         rdg_gamma_cd_llb, trpd_leewv_rdg_gamma,        &
        2312           0 :         ptend, flx_heat)
        2313             :
        2314             :   endif
        2315             :
        2316             :   ! Convert the tendencies for the dry constituents to dry air basis.
        2317     5956704 :   do m = 1, pcnst
        2318     5956704 :      if (cnst_type(m).eq.'dry') then
        2319           0 :         do k = 1, pver
        2320           0 :            do i = 1, ncol
        2321           0 :               ptend%q(i,k,m) = ptend%q(i,k,m)*state1%pdel(i,k)/state1%pdeldry(i,k)
        2322             :            end do
        2323             :         end do
        2324             :      end if
        2325             :   end do
        2326             :
        2327             :   ! Write totals to history file.
        2328     1489176 :   call outfld('EKGW', egwdffi_tot , ncol, lchnk)
        2329   659704968 :   call outfld('TTGW', ptend%s/cpairv(:,:,lchnk),  pcols, lchnk)
        2330             :
        2331     1489176 :   call outfld('UTGW_TOTAL', ptend%u, pcols, lchnk)
        2332     1489176 :   call outfld('VTGW_TOTAL', ptend%v, pcols, lchnk)
        2333             :
        2334     1489176 :   call outfld('QTGW', ptend%q(:,:,1), pcols, lchnk)
        2335     1489176 :   call outfld('CLDLIQTGW', ptend%q(:,:,ixcldliq), pcols, lchnk)
        2336     1489176 :   call outfld('CLDICETGW', ptend%q(:,:,ixcldice), pcols, lchnk)
        2337             :
        2338             :   ! Destroy objects.
        2339     1489176 :   call p%finalize()
        2340             :
        2341     2978352 : end subroutine gw_tend
        2342             :
        2343             : !==========================================================================
        2344             :
        2345           0 : subroutine gw_rdg_calc( &
        2346           0 :    type, ncol, lchnk, n_rdg, dt, &
        2347           0 :    u, v, t, p, piln, zm, zi, &
        2348           0 :    nm, ni, rhoi, kvtt, q, dse, &
        2349             :    effgw_rdg, effgw_rdg_max, &
        2350           0 :    hwdth, clngt, gbxar, &
        2351           0 :    mxdis, angll, anixy, &
        2352             :    rdg_cd_llb, trpd_leewv, &
        2353             :    ptend, flx_heat)
        2354             :
        2355     1489176 :    use coords_1d,  only: Coords1D
        2356             :    use gw_rdg,     only: gw_rdg_src, gw_rdg_belowpeak, gw_rdg_break_trap, gw_rdg_do_vdiff
        2357             :    use gw_common,  only: gw_drag_prof, energy_change
        2358             :
        2359             :    character(len=5), intent(in) :: type         ! BETA or GAMMA
        2360             :    integer,          intent(in) :: ncol         ! number of atmospheric columns
        2361             :    integer,          intent(in) :: lchnk        ! chunk identifier
        2362             :    integer,          intent(in) :: n_rdg
        2363             :    real(r8),         intent(in) :: dt           ! Time step.
        2364             :
        2365             :    real(r8),         intent(in) :: u(ncol,pver)    ! Midpoint zonal winds. ( m s-1)
        2366             :    real(r8),         intent(in) :: v(ncol,pver)    ! Midpoint meridional winds. ( m s-1)
        2367             :    real(r8),         intent(in) :: t(ncol,pver)    ! Midpoint temperatures. (K)
        2368             :    type(Coords1D),   intent(in) :: p               ! Pressure coordinates.
        2369             :    real(r8),         intent(in) :: piln(ncol,pver+1)  ! Log of interface pressures.
        2370             :    real(r8),         intent(in) :: zm(ncol,pver)   ! Midpoint altitudes above ground (m).
        2371             :    real(r8),         intent(in) :: zi(ncol,pver+1) ! Interface altitudes above ground (m).
        2372             :    real(r8),         intent(in) :: nm(ncol,pver)   ! Midpoint Brunt-Vaisalla frequencies (s-1).
        2373             :    real(r8),         intent(in) :: ni(ncol,pver+1) ! Interface Brunt-Vaisalla frequencies (s-1).
        2374             :    real(r8),         intent(in) :: rhoi(ncol,pver+1) ! Interface density (kg m-3).
        2375             :    real(r8),         intent(in) :: kvtt(ncol,pver+1) ! Molecular thermal diffusivity.
        2376             :    real(r8),         intent(in) :: q(:,:,:)        ! Constituent array.
        2377             :    real(r8),         intent(in) :: dse(ncol,pver)  ! Dry static energy.
        2378             :
        2379             :
        2380             :    real(r8),         intent(in) :: effgw_rdg       ! Tendency efficiency.
        2381             :    real(r8),         intent(in) :: effgw_rdg_max
        2382             :    real(r8),         intent(in) :: hwdth(ncol,prdg) ! width of ridges.
        2383             :    real(r8),         intent(in) :: clngt(ncol,prdg) ! length of ridges.
        2384             :    real(r8),         intent(in) :: gbxar(ncol)      ! gridbox area
        2385             :
        2386             :    real(r8),         intent(in) :: mxdis(ncol,prdg) ! Height estimate for ridge (m).
        2387             :    real(r8),         intent(in) :: angll(ncol,prdg) ! orientation of ridges.
        2388             :    real(r8),         intent(in) :: anixy(ncol,prdg) ! Anisotropy parameter.
        2389             :
        2390             :    real(r8),         intent(in) :: rdg_cd_llb      ! Drag coefficient for low-level flow
        2391             :    logical,          intent(in) :: trpd_leewv
        2392             :
        2393             :    type(physics_ptend), intent(inout):: ptend   ! Parameterization net tendencies.
        2394             :
        2395             :    real(r8),        intent(out) :: flx_heat(pcols)
        2396             :
        2397             :    !---------------------------Local storage-------------------------------
        2398             :
        2399             :    integer :: k, m, nn, istat
        2400             :
        2401           0 :    real(r8), allocatable :: tau(:,:,:)  ! wave Reynolds stress
        2402             :    ! gravity wave wind tendency for each wave
        2403           0 :    real(r8), allocatable :: gwut(:,:,:)
        2404             :    ! Wave phase speeds for each column
        2405           0 :    real(r8), allocatable :: phase_speeds(:,:)
        2406             :
        2407             :    ! Isotropic source flag [anisotropic orography].
        2408           0 :    integer  :: isoflag(ncol)
        2409             :
        2410             :    ! horiz wavenumber [anisotropic orography].
        2411           0 :    real(r8) :: kwvrdg(ncol)
        2412             :
        2413             :    ! Efficiency for a gravity wave source.
        2414           0 :    real(r8) :: effgw(ncol)
        2415             :
        2416             :    ! Indices of top gravity wave source level and lowest level where wind
        2417             :    ! tendencies are allowed.
        2418           0 :    integer :: src_level(ncol)
        2419           0 :    integer :: tend_level(ncol)
        2420           0 :    integer :: bwv_level(ncol)
        2421           0 :    integer :: tlb_level(ncol)
        2422             :
        2423             :    ! Projection of wind at midpoints and interfaces.
        2424           0 :    real(r8) :: ubm(ncol,pver)
        2425           0 :    real(r8) :: ubi(ncol,pver+1)
        2426             :
        2427             :    ! Unit vectors of source wind (zonal and meridional components).
        2428           0 :    real(r8) :: xv(ncol)
        2429           0 :    real(r8) :: yv(ncol)
        2430             :
        2431             :    ! Averages over source region.
        2432           0 :    real(r8) :: ubmsrc(ncol) ! On-ridge wind.
        2433           0 :    real(r8) :: usrc(ncol)   ! Zonal wind.
        2434           0 :    real(r8) :: vsrc(ncol)   ! Meridional wind.
        2435           0 :    real(r8) :: nsrc(ncol)   ! B-V frequency.
        2436           0 :    real(r8) :: rsrc(ncol)   ! Density.
        2437             :
        2438             :    ! normalized wavenumber
        2439           0 :    real(r8) :: m2src(ncol)
        2440             :
        2441             :    ! Top of low-level flow layer.
        2442           0 :    real(r8) :: tlb(ncol)
        2443             :
        2444             :    ! Bottom of linear wave region.
        2445           0 :    real(r8) :: bwv(ncol)
        2446             :
        2447             :    ! Froude numbers for flow/drag regimes
        2448           0 :    real(r8) :: Fr1(ncol)
        2449           0 :    real(r8) :: Fr2(ncol)
        2450           0 :    real(r8) :: Frx(ncol)
        2451             :
        2452             :    ! Wave Reynolds stresses at source level
        2453           0 :    real(r8) :: tauoro(ncol)
        2454           0 :    real(r8) :: taudsw(ncol)
        2455             :
        2456             :    ! Surface streamline displacement height for linear waves.
        2457           0 :    real(r8) :: hdspwv(ncol)
        2458             :
        2459             :    ! Surface streamline displacement height for downslope wind regime.
        2460           0 :    real(r8) :: hdspdw(ncol)
        2461             :
        2462             :    ! Wave breaking level
        2463           0 :    real(r8) :: wbr(ncol)
        2464             :
        2465           0 :    real(r8) :: utgw(ncol,pver)       ! zonal wind tendency
        2466           0 :    real(r8) :: vtgw(ncol,pver)       ! meridional wind tendency
        2467           0 :    real(r8) :: ttgw(ncol,pver)       ! temperature tendency
        2468           0 :    real(r8) :: qtgw(ncol,pver,pcnst) ! constituents tendencies
        2469             :
        2470             :    ! Effective gravity wave diffusivity at interfaces.
        2471           0 :    real(r8) :: egwdffi(ncol,pver+1)
        2472             :
        2473             :    ! Temperature tendencies from diffusion and kinetic energy.
        2474           0 :    real(r8) :: dttdf(ncol,pver)
        2475           0 :    real(r8) :: dttke(ncol,pver)
        2476             :
        2477             :    ! Wave stress in zonal/meridional direction
        2478           0 :    real(r8) :: taurx(ncol,pver+1)
        2479           0 :    real(r8) :: taurx0(ncol,pver+1)
        2480           0 :    real(r8) :: taury(ncol,pver+1)
        2481           0 :    real(r8) :: taury0(ncol,pver+1)
        2482             :    ! Provisional absolute wave stress from gw_drag_prof
        2483           0 :    real(r8) :: tau_diag(ncol,pver+1)
        2484             :
        2485             :    ! U,V tendency accumulators
        2486           0 :    real(r8) :: utrdg(ncol,pver)
        2487           0 :    real(r8) :: vtrdg(ncol,pver)
        2488           0 :    real(r8) :: ttrdg(ncol,pver)
        2489             :
        2490             :    ! Energy change used by fixer.
        2491           0 :    real(r8) :: de(ncol)
        2492             :
        2493             :    character(len=1) :: cn
        2494             :    character(len=9) :: fname(4)
        2495             :    !----------------------------------------------------------------------------
        2496             :
        2497             :    ! Allocate wavenumber fields.
        2498           0 :    allocate(tau(ncol,band_oro%ngwv:band_oro%ngwv,pver+1),stat=istat)
        2499           0 :    call alloc_err(istat,'gw_rdg_calc','tau',ncol*(band_oro%ngwv**2+1)*(pver+1))
        2500           0 :    allocate(gwut(ncol,pver,band_oro%ngwv:band_oro%ngwv),stat=istat)
        2501           0 :    call alloc_err(istat,'rdg_calc','gwut',ncol*pver*(band_oro%ngwv**2+1))
        2502           0 :    allocate(phase_speeds(ncol,band_oro%ngwv:band_oro%ngwv),stat=istat)
        2503           0 :    call alloc_err(istat,'rdg_calc','phase_speeds',ncol*(band_oro%ngwv**2+1))
        2504             :
        2505             :    ! initialize accumulated momentum fluxes and tendencies
        2506           0 :    taurx = 0._r8
        2507           0 :    taury = 0._r8
        2508           0 :    ttrdg = 0._r8
        2509           0 :    utrdg = 0._r8
        2510           0 :    vtrdg = 0._r8
        2511           0 :    tau_diag = -9999._r8
        2512             :
        2513           0 :    do nn = 1, n_rdg
        2514           0 :       kwvrdg  = 0.001_r8 / ( hwdth(:,nn) + 0.001_r8 ) ! this cant be done every time step !!!
        2515           0 :       isoflag = 0
        2516           0 :       effgw   = effgw_rdg * ( hwdth(1:ncol,nn)* clngt(1:ncol,nn) ) / gbxar(1:ncol)
        2517           0 :       effgw   = min( effgw_rdg_max , effgw )
        2518             :
        2519             :       call gw_rdg_src(ncol, band_oro, p, &
        2520             :          u, v, t, mxdis(:,nn), angll(:,nn), anixy(:,nn), kwvrdg, isoflag, zi, nm, &
        2521             :          src_level, tend_level, bwv_level, tlb_level, tau, ubm, ubi, xv, yv,  &
        2522           0 :          ubmsrc, usrc, vsrc, nsrc, rsrc, m2src, tlb, bwv, Fr1, Fr2, Frx, phase_speeds)
        2523             :
        2524             :       call gw_rdg_belowpeak(ncol, band_oro, rdg_cd_llb, &
        2525             :          t, mxdis(:,nn), anixy(:,nn), kwvrdg, &
        2526             :          zi, nm, ni, rhoi, &
        2527             :          src_level, tau, &
        2528             :          ubmsrc, nsrc, rsrc, m2src, tlb, bwv, Fr1, Fr2, Frx, &
        2529           0 :          tauoro, taudsw, hdspwv, hdspdw)
        2530             :
        2531             :
        2532             :       call gw_rdg_break_trap(ncol, band_oro, &
        2533             :          zi, nm, ni, ubm, ubi, rhoi, kwvrdg , bwv, tlb, wbr, &
        2534             :          src_level, tlb_level, hdspwv, hdspdw,  mxdis(:,nn), &
        2535             :          tauoro, taudsw, tau, &
        2536           0 :          ldo_trapped_waves=trpd_leewv)
        2537             :
        2538             :
        2539             :       call gw_drag_prof(ncol, band_oro, p, src_level, tend_level, dt, &
        2540             :          t, vramp,    &
        2541             :          piln, rhoi, nm, ni, ubm, ubi, xv, yv,   &
        2542             :          effgw, phase_speeds, kvtt, q, dse, tau, utgw, vtgw, &
        2543             :          ttgw, qtgw, egwdffi,   gwut, dttdf, dttke, &
        2544             :          kwvrdg=kwvrdg, &
        2545           0 :          satfac_in = 1._r8, lapply_vdiff=gw_rdg_do_vdiff , tau_diag=tau_diag )
        2546             :
        2547             :       ! Add the tendencies from each ridge to the totals.
        2548           0 :       do k = 1, pver
        2549             :          ! diagnostics
        2550           0 :          utrdg(:,k) = utrdg(:,k) + utgw(:,k)
        2551           0 :          vtrdg(:,k) = vtrdg(:,k) + vtgw(:,k)
        2552           0 :          ttrdg(:,k) = ttrdg(:,k) + ttgw(:,k)
        2553             :          ! physics tendencies
        2554           0 :          ptend%u(:ncol,k) = ptend%u(:ncol,k) + utgw(:,k)
        2555           0 :          ptend%v(:ncol,k) = ptend%v(:ncol,k) + vtgw(:,k)
        2556           0 :          ptend%s(:ncol,k) = ptend%s(:ncol,k) + ttgw(:,k)
        2557             :       end do
        2558             :
        2559           0 :       do m = 1, pcnst
        2560           0 :          do k = 1, pver
        2561           0 :             ptend%q(:ncol,k,m) = ptend%q(:ncol,k,m) + qtgw(:,k,m)
        2562             :          end do
        2563             :       end do
        2564             :
        2565           0 :       do k = 1, pver+1
        2566           0 :          taurx0(:,k) =  tau(:,0,k)*xv
        2567           0 :          taury0(:,k) =  tau(:,0,k)*yv
        2568           0 :          taurx(:,k)  =  taurx(:,k) + taurx0(:,k)
        2569           0 :          taury(:,k)  =  taury(:,k) + taury0(:,k)
        2570             :       end do
        2571             :
        2572           0 :       if (nn == 1) then
        2573           0 :          call outfld('BWV_HT1', bwv,     ncol, lchnk)
        2574           0 :          call outfld('TLB_HT1', tlb,     ncol, lchnk)
        2575           0 :          call outfld('WBR_HT1', wbr,     ncol, lchnk)
        2576           0 :          call outfld('TAUDSW1', taudsw,  ncol, lchnk)
        2577           0 :          call outfld('TAUORO1', tauoro,  ncol, lchnk)
        2578           0 :          call outfld('UBMSRC1', ubmsrc,  ncol, lchnk)
        2579           0 :          call outfld('USRC1',   usrc,    ncol, lchnk)
        2580           0 :          call outfld('VSRC1',   vsrc,    ncol, lchnk)
        2581           0 :          call outfld('NSRC1'  , nsrc,    ncol, lchnk)
        2582             :          ! Froude numbers
        2583           0 :          call outfld('Fr1_DIAG' , Fr1,    ncol, lchnk)
        2584           0 :          call outfld('Fr2_DIAG' , Fr2,    ncol, lchnk)
        2585           0 :          call outfld('Frx_DIAG' , Frx,    ncol, lchnk)
        2586             :          ! Ridge quantities - don't change.  Written for convenience
        2587           0 :          call outfld('MXDIS1' , mxdis(:,nn) ,  ncol, lchnk)
        2588           0 :          call outfld('ANGLL1' , angll(:,nn) ,  ncol, lchnk)
        2589           0 :          call outfld('ANIXY1' , anixy(:,nn) ,  ncol, lchnk)
        2590           0 :          call outfld('HWDTH1' , hwdth(:,nn) ,  ncol, lchnk)
        2591           0 :          call outfld('CLNGT1' , clngt(:,nn) ,  ncol, lchnk)
        2592           0 :          call outfld('GBXAR1' , gbxar ,        ncol, lchnk)
        2593           0 :          call outfld('TAUM1_DIAG' , tau_diag ,  ncol, lchnk)
        2594           0 :          call outfld('TAU1RDG'//trim(type)//'M', tau(:,0,:),  ncol, lchnk)
        2595           0 :          call outfld('UBM1'//trim(type),         ubm,         ncol, lchnk)
        2596           0 :          call outfld('UBT1RDG'//trim(type),      gwut,        ncol, lchnk)
        2597             :       end if
        2598             :
        2599           0 :       if (nn <= 6) then
        2600           0 :          write(cn, '(i1)') nn
        2601           0 :          call outfld('TAU'//cn//'RDG'//trim(type)//'X', taurx0,  ncol, lchnk)
        2602           0 :          call outfld('TAU'//cn//'RDG'//trim(type)//'Y', taury0,  ncol, lchnk)
        2603           0 :          call outfld('UT'//cn//'RDG'//trim(type),       utgw,    ncol, lchnk)
        2604           0 :          call outfld('VT'//cn//'RDG'//trim(type),       vtgw,    ncol, lchnk)
        2605             :       end if
        2606             :
        2607             :    end do ! end of loop over multiple ridges
        2608             :
        2609             :    ! Calculate energy change for output to CAM's energy checker.
        2610           0 :    call energy_change(dt, p, u, v, ptend%u(:ncol,:), &
        2611           0 :           ptend%v(:ncol,:), ptend%s(:ncol,:), de)
        2612           0 :    flx_heat(:ncol) = de
        2613             :
        2614           0 :    call outfld('TAUARDG'//trim(type)//'X', taurx,  ncol, lchnk)
        2615           0 :    call outfld('TAUARDG'//trim(type)//'Y', taury,  ncol, lchnk)
        2616             :
        2617           0 :    if (trim(type) == 'BETA') then
        2618           0 :       fname(1) = 'TAUGWX'
        2619           0 :       fname(2) = 'TAUGWY'
        2620           0 :       fname(3) = 'UTGWORO'
        2621           0 :       fname(4) = 'VTGWORO'
        2622           0 :    else if (trim(type) == 'GAMMA') then
        2623           0 :       fname(1) = 'TAURDGGMX'
        2624           0 :       fname(2) = 'TAURDGGMY'
        2625           0 :       fname(3) = 'UTRDGGM'
        2626           0 :       fname(4) = 'VTRDGGM'
        2627             :    else
        2628             :       call endrun('gw_rdg_calc: FATAL: type must be either BETA or GAMMA'&
        2629           0 :                   //' type= '//type)
        2630             :    end if
        2631             :
        2632           0 :    call outfld(fname(1), taurx(:,pver+1), ncol, lchnk)
        2633           0 :    call outfld(fname(2), taury(:,pver+1), ncol, lchnk)
        2634           0 :    call outfld(fname(3), utrdg,  ncol, lchnk)
        2635           0 :    call outfld(fname(4), vtrdg,  ncol, lchnk)
        2636           0 :    call outfld('TTGWORO', ttrdg / cpair,  ncol, lchnk)
        2637             :
        2638           0 :    deallocate(tau, gwut, phase_speeds)
        2639             :
        2640           0 : end subroutine gw_rdg_calc
        2641             :
        2642             : !==========================================================================
        2643             :
        2644             : ! Add all history fields for a gravity wave spectrum source.
        2645           0 : subroutine gw_spec_addflds(prefix, scheme, band, history_defaults)
        2646             :   use cam_history, only: addfld, add_default, register_vector_field
        2647             :
        2648             :   !------------------------------Arguments--------------------------------
        2649             :
        2650             :   ! One character prefix prepended to output fields.
        2651             :   character(len=1), intent(in) :: prefix
        2652             :   ! Gravity wave scheme name prepended to output field descriptions.
        2653             :   character(len=*), intent(in) :: scheme
        2654             :   ! Wave speeds.
        2655             :   type(GWBand), intent(in) :: band
        2656             :   ! Whether or not to call add_default for fields output by WACCM.
        2657             :   logical, intent(in) :: history_defaults
        2658             :
        2659             :   !---------------------------Local storage-------------------------------
        2660             :
        2661             :   integer :: l
        2662             :   ! 7 chars is enough for "-100.00"
        2663             :   character(len=7)  :: fnum
        2664             :   ! 10 chars is enough for "BTAUXSn32"
        2665             :   character(len=10) :: dumc1x, dumc1y
        2666             :   ! Allow 80 chars for description
        2667             :   character(len=80) dumc2
        2668             :
        2669             :   !-----------------------------------------------------------------------
        2670             :
        2671             :   ! Overall wind tendencies.
        2672             :   call addfld (trim(prefix)//'UTGWSPEC',(/ 'lev' /), 'A','m s-2', &
        2673           0 :        trim(scheme)//' U tendency - gravity wave spectrum')
        2674             :   call addfld (trim(prefix)//'VTGWSPEC',(/ 'lev' /), 'A','m s-2', &
        2675           0 :        trim(scheme)//' V tendency - gravity wave spectrum')
        2676           0 :   call register_vector_field(trim(prefix)//'UTGWSPEC',trim(prefix)//'VTGWSPEC')
        2677             :
        2678             :   call addfld (trim(prefix)//'TTGWSPEC',(/ 'lev' /), 'A','K s-1', &
        2679           0 :        trim(scheme)//' T tendency - gravity wave spectrum')
        2680             :
        2681             :   ! Wind tendencies broken across five spectral bins.
        2682             :   call addfld (trim(prefix)//'UTEND1',  (/ 'lev' /), 'A','m s-2', &
        2683           0 :        trim(scheme)//' U tendency   c < -40')
        2684             :   call addfld (trim(prefix)//'UTEND2',  (/ 'lev' /), 'A','m s-2', &
        2685           0 :        trim(scheme)//' U tendency  -40 < c < -15')
        2686             :   call addfld (trim(prefix)//'UTEND3',  (/ 'lev' /), 'A','m s-2', &
        2687           0 :        trim(scheme)//' U tendency  -15 < c <  15')
        2688             :   call addfld (trim(prefix)//'UTEND4',  (/ 'lev' /), 'A','m s-2', &
        2689           0 :        trim(scheme)//' U tendency   15 < c <  40')
        2690             :   call addfld (trim(prefix)//'UTEND5',  (/ 'lev' /), 'A','m s-2', &
        2691           0 :        trim(scheme)//' U tendency   40 < c ')
        2692             :
        2693             :   ! Reynold's stress toward each cardinal direction, and net zonal stress.
        2694             :   call addfld (trim(prefix)//'TAUE' ,   (/ 'ilev' /), 'A','Pa', &
        2695           0 :        trim(scheme)//' Eastward Reynolds stress')
        2696             :   call addfld (trim(prefix)//'TAUW' ,   (/ 'ilev' /), 'A','Pa', &
        2697           0 :        trim(scheme)//' Westward Reynolds stress')
        2698             :   call addfld (trim(prefix)//'TAUNET' , (/ 'ilev' /), 'A','Pa', &
        2699           0 :        trim(scheme)//' E+W Reynolds stress')
        2700             :   call addfld (trim(prefix)//'TAUN' ,   (/ 'ilev' /), 'A','Pa', &
        2701           0 :        trim(scheme)//' Northward Reynolds stress')
        2702             :   call addfld (trim(prefix)//'TAUS' ,   (/ 'ilev' /), 'A','Pa', &
        2703           0 :        trim(scheme)//' Southward Reynolds stress')
        2704             :
        2705             :   ! Momentum flux in each direction.
        2706             :   call addfld (trim(prefix)//'EMF',       (/ 'lev' /), 'A','Pa', &
        2707           0 :        trim(scheme)//' Eastward MF')
        2708             :   call addfld (trim(prefix)//'WMF',       (/ 'lev' /), 'A','Pa', &
        2709           0 :        trim(scheme)//' Westward MF')
        2710             :   call addfld (trim(prefix)//'NMF',       (/ 'lev' /), 'A','Pa', &
        2711           0 :        trim(scheme)//' Northward MF')
        2712             :   call addfld (trim(prefix)//'SMF',       (/ 'lev' /), 'A','Pa', &
        2713           0 :        trim(scheme)//' Southward MF')
        2714             :
        2715             :   ! Temperature tendency terms.
        2716             :   call addfld (trim(prefix)//'TTGWSDF' , (/ 'lev' /), 'A','K s-1', &
        2717           0 :        trim(scheme)//' t tendency - diffusion term')
        2718             :   call addfld (trim(prefix)//'TTGWSKE' , (/ 'lev' /), 'A','K s-1', &
        2719           0 :        trim(scheme)//' t tendency - kinetic energy conversion term')
        2720             :
        2721             :   ! Gravity wave source spectra by wave number.
        2722           0 :   do l=-band%ngwv,band%ngwv
        2723             :      ! String containing reference speed.
        2724           0 :      write (fnum,fmt='(f7.2)') band%cref(l)
        2725             :
        2726           0 :      dumc1x = tau_fld_name(l, prefix, x_not_y=.true.)
        2727           0 :      dumc1y = tau_fld_name(l, prefix, x_not_y=.false.)
        2728           0 :      dumc2 = trim(scheme)//" tau at c= "//trim(fnum)//" m s-1"
        2729           0 :      call addfld (trim(dumc1x),(/ 'lev' /), 'A','Pa',dumc2)
        2730           0 :      call addfld (trim(dumc1y),(/ 'lev' /), 'A','Pa',dumc2)
        2731             :
        2732             :   end do
        2733             :
        2734           0 :   if (history_defaults) then
        2735           0 :      call add_default(trim(prefix)//'UTGWSPEC', 1, ' ')
        2736           0 :      call add_default(trim(prefix)//'VTGWSPEC', 1, ' ')
        2737           0 :      call add_default(trim(prefix)//'TTGWSPEC', 1, ' ')
        2738           0 :      call add_default(trim(prefix)//'TAUE', 1, ' ')
        2739           0 :      call add_default(trim(prefix)//'TAUW', 1, ' ')
        2740           0 :      call add_default(trim(prefix)//'TAUNET', 1, ' ')
        2741           0 :      call add_default(trim(prefix)//'TAUN', 1, ' ')
        2742           0 :      call add_default(trim(prefix)//'TAUS', 1, ' ')
        2743             :   end if
        2744             :
        2745           0 : end subroutine gw_spec_addflds
        2746             :
        2747             : !==========================================================================
        2748             :
        2749             : ! Outputs for spectral waves.
        2750           0 : subroutine gw_spec_outflds(prefix, lchnk, ncol, band, phase_speeds, u, v, xv, yv, &
        2751           0 :      gwut, dttdf, dttke, tau, utgw, vtgw, ttgw, taucd)
        2752             :
        2753           0 :   use gw_common, only: west, east, south, north
        2754             :
        2755             :   ! One-character prefix prepended to output fields.
        2756             :   character(len=1), intent(in) :: prefix
        2757             :   ! Chunk and number of columns in the chunk.
        2758             :   integer, intent(in) :: lchnk
        2759             :   integer, intent(in) :: ncol
        2760             :   ! Wave speeds.
        2761             :   type(GWBand), intent(in) :: band
        2762             :   ! Wave phase speeds for each column.
        2763             :   real(r8), intent(in) :: phase_speeds(ncol,-band%ngwv:band%ngwv)
        2764             :   ! Winds at cell midpoints.
        2765             :   real(r8), intent(in) :: u(ncol,pver)
        2766             :   real(r8), intent(in) :: v(ncol,pver)
        2767             :   ! Unit vector in the direction of wind at source level.
        2768             :   real(r8), intent(in) :: xv(ncol)
        2769             :   real(r8), intent(in) :: yv(ncol)
        2770             :   ! Wind tendency for each wave.
        2771             :   real(r8), intent(in) :: gwut(ncol,pver,-band%ngwv:band%ngwv)
        2772             :   ! Temperature tendencies from diffusion and kinetic energy.
        2773             :   real(r8) :: dttdf(ncol,pver)
        2774             :   real(r8) :: dttke(ncol,pver)
        2775             :   ! Wave Reynolds stress.
        2776             :   real(r8), intent(in) :: tau(ncol,-band%ngwv:band%ngwv,pver)
        2777             :   ! Zonal and meridional total wind tendencies.
        2778             :   real(r8), intent(in) :: utgw(ncol,pver)
        2779             :   real(r8), intent(in) :: vtgw(ncol,pver)
        2780             :   ! Temperature tendencies.
        2781             :   real(r8), intent(in) :: ttgw(ncol,pver)
        2782             :   ! Reynolds stress for waves propagating in each cardinal direction.
        2783             :   real(r8), intent(in) :: taucd(ncol,pver+1,4)
        2784             :
        2785             :   ! Indices
        2786             :   integer :: i, k, l
        2787           0 :   integer :: ix(ncol, -band%ngwv:band%ngwv), iy(ncol, -band%ngwv:band%ngwv)
        2788           0 :   integer :: iu(ncol), iv(ncol)
        2789             :
        2790             :   ! Zonal wind tendency, broken up into five bins.
        2791           0 :   real(r8) :: utb(ncol, pver, 5)
        2792             :   ! Definition of the bin boundaries.
        2793             :   real(r8), parameter :: bounds(4) = (/ -40._r8, -15._r8, &
        2794             :        15._r8, 40._r8 /)
        2795             :
        2796             :   ! Momentum flux in the four cardinal directions.
        2797           0 :   real(r8) :: mf(ncol, pver, 4)
        2798             :
        2799             :   ! Wave stress in zonal/meridional direction
        2800           0 :   real(r8) :: taux(ncol,-band%ngwv:band%ngwv,pver)
        2801           0 :   real(r8) :: tauy(ncol,-band%ngwv:band%ngwv,pver)
        2802             :
        2803             :   ! Temporaries for output
        2804           0 :   real(r8) :: dummyx(ncol,pver)
        2805           0 :   real(r8) :: dummyy(ncol,pver)
        2806             :   ! Variable names
        2807             :   character(len=10) :: dumc1x, dumc1y
        2808             :
        2809             :
        2810             :   ! Accumulate wind tendencies binned according to phase speed.
        2811             :
        2812           0 :   utb = 0._r8
        2813             :
        2814             :   ! Find which output bin the phase speed corresponds to.
        2815           0 :   ix = find_bin(phase_speeds)
        2816             :
        2817             :   ! Put the wind tendency in that bin.
        2818           0 :   do l = -band%ngwv, band%ngwv
        2819           0 :      do k = 1, pver
        2820           0 :         do i = 1, ncol
        2821           0 :            utb(i,k,ix(i,l)) = utb(i,k,ix(i,l)) + gwut(i,k,l)
        2822             :         end do
        2823             :      end do
        2824             :   end do
        2825             :
        2826             :   ! Find just the zonal part.
        2827           0 :   do l = 1, 5
        2828           0 :      do k = 1, pver
        2829           0 :         utb(:, k, l) = utb(:, k, l) * xv
        2830             :      end do
        2831             :   end do
        2832             :
        2833           0 :   call outfld(trim(prefix)//'UTEND1', utb(:,:,1), ncol, lchnk)
        2834           0 :   call outfld(trim(prefix)//'UTEND2', utb(:,:,2), ncol, lchnk)
        2835           0 :   call outfld(trim(prefix)//'UTEND3', utb(:,:,3), ncol, lchnk)
        2836           0 :   call outfld(trim(prefix)//'UTEND4', utb(:,:,4), ncol, lchnk)
        2837           0 :   call outfld(trim(prefix)//'UTEND5', utb(:,:,5), ncol, lchnk)
        2838             :
        2839             :   ! Output temperature tendencies due to diffusion and from kinetic energy.
        2840           0 :   call outfld(trim(prefix)//'TTGWSDF', dttdf / cpair, ncol, lchnk)
        2841           0 :   call outfld(trim(prefix)//'TTGWSKE', dttke / cpair, ncol, lchnk)
        2842             :
        2843             :
        2844             :   ! Output tau broken down into zonal and meridional components.
        2845             :
        2846           0 :   taux = 0._r8
        2847           0 :   tauy = 0._r8
        2848             :
        2849             :   ! Project phase_speeds, and convert each component to a wavenumber index.
        2850             :   ! These are mappings from the wavenumber index of tau to those of taux
        2851             :   ! and tauy, respectively.
        2852           0 :   do l=-band%ngwv,band%ngwv
        2853           0 :      ix(:,l) = c_to_l(phase_speeds(:,l)*xv)
        2854           0 :      iy(:,l) = c_to_l(phase_speeds(:,l)*yv)
        2855             :   end do
        2856             :
        2857             :   ! Find projection of tau.
        2858           0 :   do k = 1, pver
        2859           0 :      do l = -band%ngwv,band%ngwv
        2860           0 :         do i = 1, ncol
        2861           0 :            taux(i,ix(i,l),k) = taux(i,ix(i,l),k) &
        2862           0 :                 + abs(tau(i,l,k)*xv(i))
        2863           0 :            tauy(i,iy(i,l),k) = tauy(i,iy(i,l),k) &
        2864           0 :                 + abs(tau(i,l,k)*yv(i))
        2865             :         end do
        2866             :      end do
        2867             :   end do
        2868             :
        2869           0 :   do l=-band%ngwv,band%ngwv
        2870             :
        2871           0 :      dummyx = taux(:,l,:)
        2872           0 :      dummyy = tauy(:,l,:)
        2873             :
        2874           0 :      dumc1x = tau_fld_name(l, prefix, x_not_y=.true.)
        2875           0 :      dumc1y = tau_fld_name(l, prefix, x_not_y=.false.)
        2876             :
        2877           0 :      call outfld(dumc1x,dummyx,ncol,lchnk)
        2878           0 :      call outfld(dumc1y,dummyy,ncol,lchnk)
        2879             :
        2880             :   enddo
        2881             :
        2882             :
        2883             :   ! Output momentum flux in each cardinal direction.
        2884           0 :   mf = 0._r8
        2885             :
        2886           0 :   do k = 1, pver
        2887             :
        2888             :      ! Convert wind speed components to wavenumber indices.
        2889           0 :      iu = c_to_l(u(:,k))
        2890           0 :      iv = c_to_l(v(:,k))
        2891             :
        2892             :      ! Sum tau components in each cardinal direction.
        2893             :      ! Split west/east and north/south based on whether wave speed exceeds
        2894             :      ! wind speed.
        2895           0 :      do l = -band%ngwv, band%ngwv
        2896             :
        2897           0 :         where (iu > l)
        2898           0 :            mf(:,k,west) = mf(:,k,west) + taux(:,l,k)
        2899             :         elsewhere
        2900           0 :            mf(:,k,east) = mf(:,k,east) + taux(:,l,k)
        2901             :         end where
        2902             :
        2903           0 :         where (iv > l)
        2904           0 :            mf(:,k,south) = mf(:,k,south) + tauy(:,l,k)
        2905             :         elsewhere
        2906           0 :            mf(:,k,north) = mf(:,k,north) + tauy(:,l,k)
        2907             :         end where
        2908             :
        2909             :      end do
        2910             :
        2911             :   end do
        2912             :
        2913           0 :   call outfld(trim(prefix)//'WMF',mf(:,:,west),ncol,lchnk)
        2914           0 :   call outfld(trim(prefix)//'EMF',mf(:,:,east),ncol,lchnk)
        2915           0 :   call outfld(trim(prefix)//'SMF',mf(:,:,south),ncol,lchnk)
        2916           0 :   call outfld(trim(prefix)//'NMF',mf(:,:,north),ncol,lchnk)
        2917             :
        2918             :   ! Simple output fields written to history file.
        2919             :   ! Total wind tendencies.
        2920           0 :   call outfld (trim(prefix)//'UTGWSPEC', utgw , ncol, lchnk)
        2921           0 :   call outfld (trim(prefix)//'VTGWSPEC', vtgw , ncol, lchnk)
        2922           0 :   call outfld (trim(prefix)//'TTGWSPEC', ttgw , ncol, lchnk)
        2923             :
        2924             :   ! Tau in each direction.
        2925           0 :   call outfld (trim(prefix)//'TAUE', taucd(:,:,east), ncol, lchnk)
        2926           0 :   call outfld (trim(prefix)//'TAUW', taucd(:,:,west), ncol, lchnk)
        2927           0 :   call outfld (trim(prefix)//'TAUN', taucd(:,:,north), ncol, lchnk)
        2928           0 :   call outfld (trim(prefix)//'TAUS', taucd(:,:,south), ncol, lchnk)
        2929             :
        2930             :   call outfld (trim(prefix)//'TAUNET', taucd(:,:,east)+taucd(:,:,west), &
        2931           0 :        ncol, lchnk)
        2932             :
        2933             : contains
        2934             :
        2935             :   ! Given a value, finds which bin marked by "bounds" the value falls
        2936             :   ! into.
        2937           0 :   elemental function find_bin(val) result(idx)
        2938             :     real(r8), intent(in) :: val
        2939             :
        2940             :     integer :: idx
        2941             :
        2942             :     ! We just have to count how many bounds are exceeded.
        2943           0 :     if (val >= 0._r8) then
        2944           0 :        idx = count(val > bounds) + 1
        2945             :     else
        2946           0 :        idx = count(val >= bounds) + 1
        2947             :     end if
        2948             :
        2949           0 :   end function find_bin
        2950             :
        2951             :   ! Convert a speed to a wavenumber between -ngwv and ngwv.
        2952           0 :   elemental function c_to_l(c) result(l)
        2953             :     real(r8), intent(in) :: c
        2954             :
        2955             :     integer :: l
        2956             :
        2957           0 :     l = min( max(int(c/band%dc),-band%ngwv), band%ngwv )
        2958             :
        2959           0 :   end function c_to_l
        2960             :
        2961             : end subroutine gw_spec_outflds
        2962             :
        2963             : !==========================================================================
        2964             :
        2965             : ! Generates names for tau output across the wave spectrum (e.g.
        2966             : ! BTAUXSn01 or TAUYSp05).
        2967             : ! Probably this should use a wavenumber dimension on one field rather
        2968             : ! than creating a ton of numbered fields.
        2969           0 : character(len=9) pure function tau_fld_name(l, prefix, x_not_y)
        2970             :   ! Wavenumber
        2971             :   integer, intent(in) :: l
        2972             :   ! Single-character prefix for output
        2973             :   character(len=1), intent(in) :: prefix
        2974             :   ! X or Y?
        2975             :   logical, intent(in) :: x_not_y
        2976             :
        2977             :   character(len=2) :: num_str
        2978             :
        2979           0 :   tau_fld_name = trim(prefix)
        2980             :
        2981           0 :   tau_fld_name = trim(tau_fld_name)//"TAU"
        2982             :
        2983           0 :   if (x_not_y) then
        2984           0 :      tau_fld_name = trim(tau_fld_name)//"XS"
        2985             :   else
        2986           0 :      tau_fld_name = trim(tau_fld_name)//"YS"
        2987             :   end if
        2988             :
        2989           0 :   if (l < 0) then
        2990           0 :      tau_fld_name = trim(tau_fld_name)//"n"
        2991             :   else
        2992           0 :      tau_fld_name = trim(tau_fld_name)//"p"
        2993             :   end if
        2994             :
        2995           0 :   write(num_str,'(I2.2)') abs(l)
        2996             :
        2997           0 :   tau_fld_name = trim(tau_fld_name)//num_str
        2998             :
        2999           0 : end function tau_fld_name
        3000             :
        3001             : !==========================================================================
        3002             :
        3003             : end module gw_drag

Generated by: LCOV version 1.14