LCOV - code coverage report
Current view: | top level - physics/cam - gw_rdg.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 303 | 0.0 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 0 | 4 | 0.0 %

              Line data    Source code
           1             : module gw_rdg
           2             :
           3             : !
           4             : ! This module handles gravity waves from orographic sources, and was
           5             : ! extracted from gw_drag in May 2013.
           6             : !
           7             : use shr_const_mod, only: pii => shr_const_pi
           8             : use gw_utils, only: r8
           9             : use gw_common, only: pver
          10             : use coords_1d, only: Coords1D
          11             : use spmd_utils,only: masterproc
          12             : use cam_abortutils, only: endrun
          13             :
          14             :
          15             : implicit none
          16             : private
          17             : save
          18             :
          19             : ! Public interface
          20             : public :: gw_rdg_readnl
          21             : public :: gw_rdg_src
          22             : public :: gw_rdg_belowpeak
          23             : public :: gw_rdg_break_trap
          24             : public :: gw_rdg_do_vdiff
          25             :
          26             : ! Tunable Parameters
          27             : !--------------------
          28             : logical            :: do_divstream
          29             :
          30             : !===========================================
          31             : ! Parameters for DS2017 (do_divstream=.T.)
          32             : !===========================================
          33             : ! Amplification factor - 1.0 for
          34             : ! high-drag/windstorm regime
          35             : real(r8), protected :: C_BetaMax_DS
          36             :
          37             : ! Max Ratio Fr2:Fr1 - 1.0
          38             : real(r8), protected :: C_GammaMax
          39             :
          40             : ! Normalized limits  for Fr2(Frx) function
          41             : real(r8), protected :: Frx0
          42             : real(r8), protected :: Frx1
          43             :
          44             :
          45             : !===========================================
          46             : ! Parameters for SM2000
          47             : !===========================================
          48             : ! Amplification factor - 1.0 for
          49             : ! high-drag/windstorm regime
          50             : real(r8), protected :: C_BetaMax_SM
          51             :
          52             :
          53             :
          54             : ! NOTE: Critical inverse Froude number Fr_c is
          55             : ! 1./(SQRT(2.)~0.707 in SM2000
          56             : ! (should be <= 1)
          57             : real(r8), protected :: Fr_c
          58             :
          59             : logical, protected :: gw_rdg_do_vdiff=.true.
          60             :
          61             : logical :: do_smooth_regimes
          62             : logical :: do_adjust_tauoro
          63             : logical :: do_backward_compat
          64             :
          65             :
          66             : ! Limiters (min/max values)
          67             : ! min surface displacement height for orographic waves (m)
          68             : real(r8), protected :: orohmin
          69             : ! min wind speed for orographic waves
          70             : real(r8), protected :: orovmin
          71             : ! min stratification allowing wave behavior
          72             : real(r8), protected :: orostratmin
          73             : ! min stratification allowing wave behavior
          74             : real(r8), protected :: orom2min
          75             :
          76             : !==========================================================================
          77             : contains
          78             : !==========================================================================
          79             :
          80           0 : subroutine gw_rdg_readnl(nlfile)
          81             :   use namelist_utils,  only: find_group_name
          82             :   use units,           only: getunit, freeunit
          83             :   use spmd_utils,      only: mpicom, mstrid=>masterprocid, mpi_real8, mpi_logical
          84             :
          85             :   ! File containing namelist input.
          86             :   character(len=*), intent(in) :: nlfile
          87             :
          88             :   ! Local variables
          89             :   integer :: unitn, ierr
          90             :   character(len=*), parameter :: sub = 'gw_rdg_readnl'
          91             :
          92             :   logical ::  gw_rdg_do_divstream, gw_rdg_do_smooth_regimes, gw_rdg_do_adjust_tauoro, &
          93             :               gw_rdg_do_backward_compat
          94             :
          95             :
          96             :   real(r8) :: gw_rdg_C_BetaMax_DS, gw_rdg_C_GammaMax, &
          97             :               gw_rdg_Frx0, gw_rdg_Frx1, gw_rdg_C_BetaMax_SM, gw_rdg_Fr_c, &
          98             :               gw_rdg_orohmin, gw_rdg_orovmin, gw_rdg_orostratmin, gw_rdg_orom2min
          99             :
         100             :   namelist /gw_rdg_nl/ gw_rdg_do_divstream, gw_rdg_C_BetaMax_DS, gw_rdg_C_GammaMax, &
         101             :                        gw_rdg_Frx0, gw_rdg_Frx1, gw_rdg_C_BetaMax_SM, gw_rdg_Fr_c, &
         102             :                        gw_rdg_do_smooth_regimes, gw_rdg_do_adjust_tauoro, &
         103             :                        gw_rdg_do_backward_compat, gw_rdg_orohmin, gw_rdg_orovmin, &
         104             :                        gw_rdg_orostratmin, gw_rdg_orom2min, gw_rdg_do_vdiff
         105             :
         106             :   !----------------------------------------------------------------------
         107             :
         108           0 :   if (masterproc) then
         109           0 :      unitn = getunit()
         110           0 :      open( unitn, file=trim(nlfile), status='old' )
         111           0 :      call find_group_name(unitn, 'gw_rdg_nl', status=ierr)
         112           0 :      if (ierr == 0) then
         113           0 :         read(unitn, gw_rdg_nl, iostat=ierr)
         114           0 :         if (ierr /= 0) then
         115           0 :            call endrun(sub // ':: ERROR reading namelist')
         116             :         end if
         117             :      end if
         118           0 :      close(unitn)
         119           0 :      call freeunit(unitn)
         120             :
         121             :      ! Set the local variables
         122           0 :      do_divstream        = gw_rdg_do_divstream
         123           0 :      C_BetaMax_DS        = gw_rdg_C_BetaMax_DS
         124           0 :      C_GammaMax          = gw_rdg_C_GammaMax
         125           0 :      Frx0                = gw_rdg_Frx0
         126           0 :      Frx1                = gw_rdg_Frx1
         127           0 :      C_BetaMax_SM        = gw_rdg_C_BetaMax_SM
         128           0 :      Fr_c                = gw_rdg_Fr_c
         129           0 :      do_smooth_regimes   = gw_rdg_do_smooth_regimes
         130           0 :      do_adjust_tauoro    = gw_rdg_do_adjust_tauoro
         131           0 :      do_backward_compat  = gw_rdg_do_backward_compat
         132           0 :      orohmin             = gw_rdg_orohmin
         133           0 :      orovmin             = gw_rdg_orovmin
         134           0 :      orostratmin         = gw_rdg_orostratmin
         135           0 :      orom2min            = gw_rdg_orom2min
         136             :   end if
         137             :
         138             :   ! Broadcast the local variables
         139             :
         140           0 :   call mpi_bcast(do_divstream, 1, mpi_logical, mstrid, mpicom, ierr)
         141           0 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: do_divstream")
         142           0 :   call mpi_bcast(do_smooth_regimes, 1, mpi_logical, mstrid, mpicom, ierr)
         143           0 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: do_smooth_regimes")
         144           0 :   call mpi_bcast(do_adjust_tauoro, 1, mpi_logical, mstrid, mpicom, ierr)
         145           0 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: do_adjust_tauoro")
         146           0 :   call mpi_bcast(do_backward_compat, 1, mpi_logical, mstrid, mpicom, ierr)
         147           0 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: do_backward_compat")
         148             :
         149           0 :   call mpi_bcast(C_BetaMax_DS, 1, mpi_real8, mstrid, mpicom, ierr)
         150           0 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: C_BetaMax_DS")
         151           0 :   call mpi_bcast(C_GammaMax, 1, mpi_real8, mstrid, mpicom, ierr)
         152           0 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: C_GammaMax")
         153           0 :   call mpi_bcast(Frx0, 1, mpi_real8, mstrid, mpicom, ierr)
         154           0 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: Frx0")
         155           0 :   call mpi_bcast(Frx1, 1, mpi_real8, mstrid, mpicom, ierr)
         156           0 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: Frx1")
         157           0 :   call mpi_bcast(C_BetaMax_SM, 1, mpi_real8, mstrid, mpicom, ierr)
         158           0 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: C_BetaMax_SM")
         159           0 :   call mpi_bcast(Fr_c, 1, mpi_real8, mstrid, mpicom, ierr)
         160           0 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: Fr_c")
         161           0 :   call mpi_bcast(orohmin, 1, mpi_real8, mstrid, mpicom, ierr)
         162           0 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: orohmin")
         163           0 :   call mpi_bcast(orovmin, 1, mpi_real8, mstrid, mpicom, ierr)
         164           0 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: orovmin")
         165           0 :   call mpi_bcast(orostratmin, 1, mpi_real8, mstrid, mpicom, ierr)
         166           0 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: orostratmin")
         167           0 :   call mpi_bcast(orom2min, 1, mpi_real8, mstrid, mpicom, ierr)
         168           0 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: orom2min")
         169             :
         170           0 :   call mpi_bcast(gw_rdg_do_vdiff, 1, mpi_logical, mstrid, mpicom, ierr)
         171           0 :   if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_rdg_do_vdiff")
         172             :
         173           0 :   if (Fr_c > 1.0_r8) call endrun(sub//": FATAL: Fr_c must be <= 1")
         174             :
         175           0 : end subroutine gw_rdg_readnl
         176             :
         177             :
         178           0 : subroutine gw_rdg_src(ncol, band, p, &
         179           0 :      u, v, t, mxdis, angxy, anixy, kwvrdg, iso, zi, nm, &
         180           0 :      src_level, tend_level, bwv_level ,tlb_level , tau, ubm, ubi, xv, yv,  &
         181           0 :      ubmsrc, usrc, vsrc, nsrc, rsrc, m2src, tlb, bwv, Fr1, Fr2, Frx, c)
         182             :   use gw_common, only: rair, GWBand
         183             :   use gw_utils, only:  dot_2d, midpoint_interp
         184             :   !-----------------------------------------------------------------------
         185             :   ! Orographic source for multiple gravity wave drag parameterization.
         186             :   !
         187             :   ! The stress is returned for a single wave with c=0, over orography.
         188             :   ! For points where the orographic variance is small (including ocean),
         189             :   ! the returned stress is zero.
         190             :   !------------------------------Arguments--------------------------------
         191             :   ! Column dimension.
         192             :   integer, intent(in) :: ncol
         193             :
         194             :   ! Band to emit orographic waves in.
         195             :   ! Regardless, we will only ever emit into l = 0.
         196             :   type(GWBand), intent(in) :: band
         197             :   ! Pressure coordinates.
         198             :   type(Coords1D), intent(in) :: p
         199             :
         200             :
         201             :   ! Midpoint zonal/meridional winds. ( m s-1)
         202             :   real(r8), intent(in) :: u(ncol,pver), v(ncol,pver)
         203             :   ! Midpoint temperatures. (K)
         204             :   real(r8), intent(in) :: t(ncol,pver)
         205             :   ! Height estimate for ridge (m) [anisotropic orography].
         206             :   real(r8), intent(in) :: mxdis(ncol)
         207             :   ! Angle of ridge axis w/resp to north (degrees) [anisotropic orography].
         208             :   real(r8), intent(in) :: angxy(ncol)
         209             :   ! Anisotropy parameter [anisotropic orography].
         210             :   real(r8), intent(in) :: anixy(ncol)
         211             :   ! horiz wavenumber [anisotropic orography].
         212             :   real(r8), intent(in) :: kwvrdg(ncol)
         213             :   ! Isotropic source flag [anisotropic orography].
         214             :   integer, intent(in)  :: iso(ncol)
         215             :   ! Interface altitudes above ground (m).
         216             :   real(r8), intent(in) :: zi(ncol,pver+1)
         217             :   ! Midpoint Brunt-Vaisalla frequencies (s-1).
         218             :   real(r8), intent(in) :: nm(ncol,pver)
         219             :
         220             :   ! Indices of top gravity wave source level and lowest level where wind
         221             :   ! tendencies are allowed.
         222             :   integer, intent(out) :: src_level(ncol)
         223             :   integer, intent(out) :: tend_level(ncol)
         224             :   integer, intent(out) :: bwv_level(ncol),tlb_level(ncol)
         225             :
         226             :   ! Averages over source region.
         227             :   real(r8), intent(out) :: nsrc(ncol) ! B-V frequency.
         228             :   real(r8), intent(out) :: rsrc(ncol) ! Density.
         229             :   real(r8), intent(out) :: usrc(ncol) ! Zonal wind.
         230             :   real(r8), intent(out) :: vsrc(ncol) ! Meridional wind.
         231             :   real(r8), intent(out) :: ubmsrc(ncol) ! On-ridge wind.
         232             :   ! Top of low-level flow layer.
         233             :   real(r8), intent(out) :: tlb(ncol)
         234             :   ! Bottom of linear wave region.
         235             :   real(r8), intent(out) :: bwv(ncol)
         236             :   ! normalized wavenumber
         237             :   real(r8), intent(out) :: m2src(ncol)
         238             :
         239             :
         240             :   ! Wave Reynolds stress.
         241             :   real(r8), intent(out) :: tau(ncol,-band%ngwv:band%ngwv,pver+1)
         242             :   ! Projection of wind at midpoints and interfaces.
         243             :   real(r8), intent(out) :: ubm(ncol,pver), ubi(ncol,pver+1)
         244             :   ! Unit vectors of source wind (zonal and meridional components).
         245             :   real(r8), intent(out) :: xv(ncol), yv(ncol)
         246             :   ! Phase speeds.
         247             :   real(r8), intent(out) :: c(ncol,-band%ngwv:band%ngwv)
         248             :   ! Froude numbers for flow/drag regimes
         249             :   real(r8), intent(out) :: Fr1(ncol), Fr2(ncol), Frx(ncol)
         250             :
         251             :   !---------------------------Local Storage-------------------------------
         252             :   ! Column and level indices.
         253             :   integer :: i, k
         254             :
         255             :   ! Surface streamline displacement height (2*sgh).
         256           0 :   real(r8) :: hdsp(ncol)
         257             :
         258             :   ! Difference in interface pressure across source region.
         259           0 :   real(r8) :: dpsrc(ncol)
         260             :   ! Thickness of downslope wind region.
         261           0 :   real(r8) :: ddw(ncol)
         262             :   ! Thickness of linear wave region.
         263           0 :   real(r8) :: dwv(ncol)
         264             :   ! Wind speed in source region.
         265           0 :   real(r8) :: wmsrc(ncol)
         266             :
         267           0 :   real(r8) :: ragl(ncol)
         268             :
         269             : !--------------------------------------------------------------------------
         270             : ! Check that ngwav is equal to zero, otherwise end the job
         271             : !--------------------------------------------------------------------------
         272           0 :   if (band%ngwv /= 0) call endrun(' gw_rdg_src :: ERROR - band%ngwv must be zero and it is not')
         273             :
         274             : !--------------------------------------------------------------------------
         275             : ! Average the basic state variables for the wave source over the depth of
         276             : ! the orographic standard deviation. Here we assume that the appropiate
         277             : ! values of wind, stability, etc. for determining the wave source are
         278             : ! averages over the depth of the atmosphere penterated by the typical
         279             : ! mountain.
         280             : ! Reduces to the bottom midpoint values when mxdis=0, such as over ocean.
         281             : !--------------------------------------------------------------------------
         282             :
         283           0 :   hdsp      = mxdis ! no longer multipied by 2
         284           0 :   src_level = pver+1
         285           0 :   bwv_level = -1
         286           0 :   tlb_level = -1
         287             :
         288           0 :   tau(:,0,:) = 0.0_r8
         289             :
         290             :   ! Find depth of "source layer" for mountain waves
         291             :   ! i.e., between ground and mountain top
         292           0 :   do k = pver, 1, -1
         293           0 :      do i = 1, ncol
         294             :         ! Need to have h >= z(k+1) here or code will bomb when h=0.
         295           0 :         if ( (hdsp(i) >= zi(i,k+1)) .and. (hdsp(i) < zi(i,k))   ) then
         296           0 :            src_level(i) = k
         297             :         end if
         298             :      end do
         299             :   end do
         300             :
         301           0 :   rsrc = 0._r8
         302           0 :   usrc = 0._r8
         303           0 :   vsrc = 0._r8
         304           0 :   nsrc = 0._r8
         305           0 :   do i = 1, ncol
         306           0 :       do k = pver, src_level(i), -1
         307           0 :            rsrc(i) = rsrc(i) + p%mid(i,k) / (rair*t(i,k))* p%del(i,k)
         308           0 :            usrc(i) = usrc(i) + u(i,k) * p%del(i,k)
         309           0 :            vsrc(i) = vsrc(i) + v(i,k) * p%del(i,k)
         310           0 :            nsrc(i) = nsrc(i) + nm(i,k)* p%del(i,k)
         311             :      end do
         312             :   end do
         313             :
         314             :
         315           0 :   do i = 1, ncol
         316           0 :      dpsrc(i) = p%ifc(i,pver+1) - p%ifc(i,src_level(i))
         317             :   end do
         318             :
         319           0 :   rsrc = rsrc / dpsrc
         320           0 :   usrc = usrc / dpsrc
         321           0 :   vsrc = vsrc / dpsrc
         322           0 :   nsrc = nsrc / dpsrc
         323             :
         324           0 :   wmsrc = sqrt( usrc**2 + vsrc**2 )
         325             :
         326             :
         327             :   ! Get the unit vector components
         328             :   ! Want agl=0 with U>0 to give xv=1
         329             :
         330           0 :   ragl = angxy * pii/180._r8
         331             :
         332             :   ! protect from wierd "bad" angles
         333             :   ! that may occur if hdsp is zero
         334           0 :   where( hdsp <= orohmin )
         335             :      ragl = 0._r8
         336             :   end where
         337             :
         338           0 :   yv   =-sin( ragl )
         339           0 :   xv   = cos( ragl )
         340             :
         341             :
         342             :   ! Kluge in possible "isotropic" obstacle.
         343           0 :   where( ( iso == 1 ) .and. (wmsrc > orovmin) )
         344             :        xv = usrc/wmsrc
         345             :        yv = vsrc/wmsrc
         346             :   end where
         347             :
         348             :
         349             :   ! Project the local wind at midpoints into the on-ridge direction
         350           0 :   do k = 1, pver
         351           0 :      ubm(:,k) = dot_2d(u(:,k), v(:,k), xv, yv)
         352             :   end do
         353           0 :   ubmsrc = dot_2d(usrc , vsrc , xv, yv)
         354             :
         355             :   ! Ensure on-ridge wind is positive at source level
         356           0 :   do k = 1, pver
         357           0 :      ubm(:,k) = sign( ubmsrc*0._r8+1._r8 , ubmsrc ) *  ubm(:,k)
         358             :   end do
         359             :
         360             :                   ! Sean says just use 1._r8 as
         361             :                   ! first argument
         362           0 :   xv  = sign( ubmsrc*0._r8+1._r8 , ubmsrc ) *  xv
         363           0 :   yv  = sign( ubmsrc*0._r8+1._r8 , ubmsrc ) *  yv
         364             :
         365             :   ! Now make ubmsrc positive and protect
         366             :   ! against zero
         367           0 :   ubmsrc = abs(ubmsrc)
         368           0 :   ubmsrc = max( 0.01_r8 , ubmsrc )
         369             :
         370             :
         371             :   ! The minimum stratification allowing GW behavior
         372             :   ! should really depend on horizontal scale since
         373             :   !
         374             :   !      m^2 ~ (N/U)^2 - k^2
         375             :   !
         376             :   ! Should also think about parameterizing
         377             :   ! trapped lee-waves.
         378             :
         379             :
         380             :   ! This needs to be made constistent with later
         381             :   ! treatment of nonhydrostatic effects.
         382           0 :   m2src = ( (nsrc/(ubmsrc+0.01_r8))**2 - kwvrdg**2 ) /((nsrc/(ubmsrc+0.01_r8))**2)
         383             :
         384             :
         385             :   !-------------------------------------------------------------
         386             :   ! Calculate provisional limits (in Z [m]) for 3 regimes. This
         387             :   ! will modified later if wave breaking or trapping are
         388             :   ! diagnosed
         389             :   !
         390             :   !                                            ^
         391             :   !                                            | *** linear propagation ***
         392             :   !  (H) -------- mountain top -------------   | *** or wave breaking  ****
         393             :   !                                            | *** regimes  *************
         394             :   ! (BWV)------ bottom of linear waves ----    |
         395             :   !                    :                       |
         396             :   !                 *******                    |
         397             :   !                    :                       |
         398             :   ! (TLB)--- top of flow diversion layer---    '
         399             :   !                   :
         400             :   !        **** flow diversion *****
         401             :   !                    :
         402             :   !============================================
         403             :
         404             :   !============================================
         405             :   ! For Dividing streamline para (DS2017)
         406             :   !--------------------------------------------
         407             :   ! High-drag downslope wind regime exists
         408             :   ! between bottom of linear waves and top of
         409             :   ! flow diversion. Linear waves can only
         410             :   ! attain vertical displacment of f1*U/N. So,
         411             :   ! bottom of linear waves is given by
         412             :   !
         413             :   !        BWV = H - Fr1*U/N
         414             :   !
         415             :   ! Downslope wind layer begins at BWV and
         416             :   ! extends below it until some maximum high
         417             :   ! drag obstacle height Fr2*U/N is attained
         418             :   ! (where Fr2 >= f1).  Below downslope wind
         419             :   ! there is flow diversion, so top of
         420             :   ! diversion layer (TLB) is equivalent to
         421             :   ! bottom of downslope wind layer and is;
         422             :   !
         423             :   !       TLB = H - Fr2*U/N
         424             :   !
         425             :   !-----------------------------------------
         426             :
         427             :   ! Critical inverse Froude number
         428             :   !-----------------------------------------------
         429           0 :   Fr1(:) = Fr_c * 1.00_r8
         430           0 :   Frx(:) = hdsp(:)*nsrc(:)/abs( ubmsrc(:) ) / Fr_c
         431             :
         432           0 :   if ( do_divstream ) then
         433             :      !------------------------------------------------
         434             :      ! Calculate Fr2(Frx) for DS2017
         435             :      !------------------------------------------------
         436           0 :      where(Frx <= Frx0)
         437             :           Fr2(:) = Fr1(:) + Fr1(:)* C_GammaMax * anixy(:)
         438             :      elsewhere((Frx > Frx0).and.(Frx <= Frx1) )
         439             :           Fr2(:) = Fr1(:) + Fr1(:)* C_GammaMax * anixy(:) &
         440             :                    * (Frx1 - Frx(:))/(Frx1-Frx0)
         441             :      elsewhere(Frx > Frx1)
         442             :           Fr2(:)=Fr1(:)
         443             :      endwhere
         444             :   else
         445             :   !------------------------------------------
         446             :   ! Regime distinctions entirely carried by
         447             :   ! amplification of taudsw (next subr)
         448             :   !------------------------------------------
         449           0 :      Fr2(:)=Fr1(:)
         450             :   end if
         451             :
         452             :
         453             :
         454           0 :   where( m2src > orom2min )
         455             :      ddw  = Fr2 * ( abs(ubmsrc) )/nsrc
         456             :   elsewhere
         457             :      ddw  = 0._r8
         458             :   endwhere
         459             :
         460             :
         461             :   ! If TLB is less than zero then obstacle is not
         462             :   ! high enough to produce an low-level diversion layer
         463           0 :   tlb = mxdis - ddw
         464           0 :   where( tlb < 0._r8)
         465             :      tlb = 0._r8
         466             :   endwhere
         467           0 :   do k = pver, pver/2, -1
         468           0 :      do i = 1, ncol
         469           0 :          if ( (tlb(i) > zi(i,k+1)) .and. (tlb(i) <= zi(i,k))   ) then
         470           0 :            tlb_level(i) = k
         471             :         end if
         472             :      end do
         473             :   end do
         474             :
         475             :
         476             :   ! Find *BOTTOM* of linear wave layer (BWV)
         477             :   !where ( nsrc > orostratmin )
         478           0 :   where( m2src > orom2min )
         479             :       dwv  = Fr1 * ( abs(ubmsrc) )/nsrc
         480             :   elsewhere
         481             :      dwv  = -9.999e9_r8 ! if weak strat - no waves
         482             :   endwhere
         483             :
         484           0 :   bwv = mxdis - dwv
         485           0 :   where(( bwv < 0._r8) .or. (dwv < 0._r8) )
         486             :      bwv = 0._r8
         487             :   endwhere
         488           0 :   do k = pver,1, -1
         489           0 :      do i = 1, ncol
         490           0 :         if ( (bwv(i) > zi(i,k+1)) .and. (bwv(i) <= zi(i,k))   ) then
         491           0 :            bwv_level(i) = k+1
         492             :         end if
         493             :      end do
         494             :   end do
         495             :
         496             :
         497             :
         498             :   ! Compute the interface wind projection by averaging the midpoint winds.
         499             :   ! Use the top level wind at the top interface.
         500           0 :   ubi(:,1) = ubm(:,1)
         501           0 :   ubi(:,2:pver) = midpoint_interp(ubm)
         502           0 :   ubi(:,pver+1) = ubm(:,pver)
         503             :
         504             :   ! Allow wind tendencies all the way to the model bottom.
         505           0 :   tend_level = pver
         506             :
         507             :   ! No spectrum; phase speed is just 0.
         508           0 :   c = 0._r8
         509             :
         510           0 :   where( m2src < orom2min )
         511             :      tlb = mxdis
         512             :      tlb_level = src_level
         513             :   endwhere
         514             :
         515             :
         516           0 : end subroutine gw_rdg_src
         517             :
         518             :
         519             : !==========================================================================
         520             :
         521           0 : subroutine gw_rdg_belowpeak(ncol, band, rdg_cd_llb, &
         522           0 :      t, mxdis, anixy, kwvrdg, zi, nm, ni, rhoi, &
         523           0 :      src_level , tau,  &
         524           0 :      ubmsrc, nsrc, rsrc, m2src,tlb,bwv,Fr1,Fr2,Frx, &
         525           0 :      tauoro,taudsw, hdspwv,hdspdw  )
         526             :
         527             :   use gw_common, only: GWBand
         528             :   !-----------------------------------------------------------------------
         529             :   ! Orographic source for multiple gravity wave drag parameterization.
         530             :   !
         531             :   ! The stress is returned for a single wave with c=0, over orography.
         532             :   ! For points where the orographic variance is small (including ocean),
         533             :   ! the returned stress is zero.
         534             :   !------------------------------Arguments--------------------------------
         535             :   ! Column dimension.
         536             :   integer, intent(in) :: ncol
         537             :   ! Band to emit orographic waves in.
         538             :   ! Regardless, we will only ever emit into l = 0.
         539             :   type(GWBand), intent(in) :: band
         540             :   ! Drag coefficient for low-level flow
         541             :   real(r8), intent(in) :: rdg_cd_llb
         542             :
         543             :
         544             :   ! Midpoint temperatures. (K)
         545             :   real(r8), intent(in) :: t(ncol,pver)
         546             :   ! Height estimate for ridge (m) [anisotropic orography].
         547             :   real(r8), intent(in) :: mxdis(ncol)
         548             :   ! Anisotropy parameter [0-1] [anisotropic orography].
         549             :   real(r8), intent(in) :: anixy(ncol)
         550             :   ! Inverse cross-ridge lengthscale (m-1) [anisotropic orography].
         551             :   real(r8), intent(inout) :: kwvrdg(ncol)
         552             :   ! Interface altitudes above ground (m).
         553             :   real(r8), intent(in) :: zi(ncol,pver+1)
         554             :   ! Midpoint Brunt-Vaisalla frequencies (s-1).
         555             :   real(r8), intent(in) :: nm(ncol,pver)
         556             :   ! Interface Brunt-Vaisalla frequencies (s-1).
         557             :   real(r8), intent(in) :: ni(ncol,pver+1)
         558             :   ! Interface density (kg m-3).
         559             :   real(r8), intent(in) :: rhoi(ncol,pver+1)
         560             :
         561             :   ! Indices of top gravity wave source level
         562             :   integer, intent(inout) :: src_level(ncol)
         563             :
         564             :   ! Wave Reynolds stress.
         565             :   real(r8), intent(inout) :: tau(ncol,-band%ngwv:band%ngwv,pver+1)
         566             :   ! Top of low-level flow layer.
         567             :   real(r8), intent(in)    :: tlb(ncol)
         568             :   ! Bottom of linear wave region.
         569             :   real(r8), intent(in)    :: bwv(ncol)
         570             :   ! surface stress from linear waves.
         571             :   real(r8), intent(out) :: tauoro(ncol)
         572             :   ! surface stress for downslope wind regime.
         573             :   real(r8), intent(out) :: taudsw(ncol)
         574             :
         575             :   ! Surface streamline displacement height for linear waves.
         576             :   real(r8), intent(out) :: hdspwv(ncol)
         577             :   ! Surface streamline displacement height for downslope wind regime.
         578             :   real(r8), intent(out) :: hdspdw(ncol)
         579             :
         580             :
         581             :
         582             :   ! Froude numbers for flow/drag regimes
         583             :   real(r8), intent(in) :: Fr1(ncol), Fr2(ncol),Frx(ncol)
         584             :
         585             :   ! Averages over source region.
         586             :   real(r8), intent(in) :: m2src(ncol) ! normalized non-hydro wavenumber
         587             :   real(r8), intent(in) :: nsrc(ncol)  ! B-V frequency.
         588             :   real(r8), intent(in) :: rsrc(ncol)  ! Density.
         589             :   real(r8), intent(in) :: ubmsrc(ncol) ! On-ridge wind.
         590             :
         591             :
         592             :   !logical, intent(in), optional :: forcetlb
         593             :
         594             :   !---------------------------Local Storage-------------------------------
         595             :   ! Column and level indices.
         596             :   integer :: i, k
         597             :
         598           0 :   real(r8) :: Coeff_LB(ncol),tausat,ubsrcx(ncol),dswamp
         599           0 :   real(r8) :: taulin(ncol),BetaMax
         600             :
         601             :   ! ubsrcx introduced to account for situations with high shear, strong strat.
         602           0 :   do i = 1, ncol
         603           0 :         ubsrcx(i)    = max( ubmsrc(i)  , 0._r8 )
         604             :   end do
         605             :
         606           0 :   do i = 1, ncol
         607           0 :      if ( m2src(i) > orom2min )   then
         608           0 :         hdspwv(i) = min( mxdis(i) , Fr1(i) * ubsrcx(i) / nsrc(i) )
         609             :      else
         610           0 :         hdspwv(i) = 0._r8
         611             :      end if
         612             :   end do
         613             :
         614           0 :   if (do_divstream) then
         615           0 :      do i = 1, ncol
         616           0 :         if ( m2src(i) > orom2min )   then
         617           0 :            hdspdw(i) = min( mxdis(i) , Fr2(i) * ubsrcx(i) / nsrc(i) )
         618             :         else
         619           0 :            hdspdw(i) = 0._r8
         620             :         end if
         621             :      end do
         622             :   else
         623           0 :      do i = 1, ncol
         624             :         ! Needed only to mark where a DSW occurs
         625           0 :         if ( m2src(i) > orom2min )   then
         626           0 :            hdspdw(i) = mxdis(i)
         627             :         else
         628           0 :            hdspdw(i) = 0._r8
         629             :         end if
         630             :      end do
         631             :   end if
         632             :
         633             :   ! Calculate form drag coefficient ("CD")
         634             :   !--------------------------------------
         635           0 :   Coeff_LB = rdg_cd_llb*anixy
         636             :
         637             :   ! Determine the orographic c=0 source term following McFarlane (1987).
         638             :   ! Set the source top interface index to pver, if the orographic term is
         639             :   ! zero.
         640             :   !
         641             :   ! This formula is basically from
         642             :   !
         643             :   !      tau(src) = rho * u' * w'
         644             :   ! where
         645             :   !      u' ~ N*h'  and w' ~ U*h'/b  (b="breite")
         646             :   !
         647             :   ! and 1/b has been replaced with k (kwvrdg)
         648             :   !
         649           0 :   do i = 1, ncol
         650           0 :      if ( ( src_level(i) > 0 ) .and. ( m2src(i) > orom2min ) ) then
         651             :         tauoro(i) = kwvrdg(i) * ( hdspwv(i)**2 ) * rsrc(i) * nsrc(i) &
         652           0 :              * ubsrcx(i)
         653             :         taudsw(i) = kwvrdg(i) * ( hdspdw(i)**2 ) * rsrc(i) * nsrc(i) &
         654           0 :              * ubsrcx(i)
         655             :      else
         656           0 :         tauoro(i) = 0._r8
         657           0 :         taudsw(i) = 0._r8
         658             :      end if
         659             :   end do
         660             :
         661           0 :   if (do_divstream) then
         662           0 :      do i = 1, ncol
         663           0 :            taulin(i) = 0._r8
         664             :      end do
         665             :   !---------------------------------------
         666             :   ! Need linear drag when divstream is not used is used
         667             :   !---------------------------------------
         668             :   else
         669           0 :      do i = 1, ncol
         670           0 :         if ( ( src_level(i) > 0 ) .and. ( m2src(i) > orom2min ) ) then
         671             :            taulin(i) = kwvrdg(i) * ( mxdis(i)**2 ) * rsrc(i) * nsrc(i) &
         672           0 :                 * ubsrcx(i)
         673             :         else
         674           0 :            taulin(i) = 0._r8
         675             :         end if
         676             :      end do
         677             :   end if
         678             :
         679           0 :   if ( do_divstream ) then
         680             :   ! Amplify DSW between Frx=1. and Frx=Frx1
         681           0 :      do i = 1,ncol
         682           0 :         dswamp=0._r8
         683           0 :         BetaMax   = C_BetaMax_DS * anixy(i)
         684           0 :         if ( (Frx(i)>1._r8).and.(Frx(i)<=Frx1)) then
         685           0 :            dswamp = (Frx(i)-1._r8)*(Frx1-Frx(i))/(0.25_r8*(Frx1-1._r8)**2)
         686             :         end if
         687           0 :         taudsw(i) = (1._r8 + BetaMax*dswamp)*taudsw(i)
         688             :      end do
         689             :   else
         690             :   !-------------------
         691             :   ! Scinocca&McFarlane
         692             :   !--------------------
         693           0 :      do i = 1, ncol
         694           0 :         BetaMax   = C_BetaMax_SM * anixy(i)
         695           0 :         if ( (Frx(i) >=1._r8) .and. (Frx(i) < 1.5_r8) ) then
         696           0 :            dswamp = 2._r8 * BetaMax * (Frx(i) -1._r8)
         697           0 :         else if ( ( Frx(i) >= 1.5_r8 ) .and. (Frx(i) < 3._r8 ) ) then
         698             :            dswamp = ( 1._r8 + BetaMax - (0.666_r8**2) ) * ( 0.666_r8*(3._r8 - Frx(i) ))**2  &
         699           0 :                       + ( 1._r8 / Frx(i) )**2  -1._r8
         700             :         else
         701             :            dswamp    = 0._r8
         702             :         end if
         703           0 :         if ( (Frx(i) >=1._r8) .and. (Frx(i) < 3._r8) ) then
         704           0 :           taudsw(i) = (1._r8 + dswamp )*taulin(i) - tauoro(i)
         705             :         else
         706           0 :           taudsw(i) = 0._r8
         707             :         endif
         708             :         ! This code defines "taudsw" as SUM of freely-propagating
         709             :         ! DSW enhancement. Different than in SM2000
         710           0 :         taudsw(i) = taudsw(i) + tauoro(i)
         711             :      end do
         712             :  !----------------------------------------------------
         713             :   end if
         714             :
         715             :
         716           0 :   do i = 1, ncol
         717           0 :      if ( m2src(i) > orom2min )   then
         718           0 :         where ( ( zi(i,:) < mxdis(i) ) .and. ( zi(i,:) >= bwv(i) ) )
         719           0 :              tau(i,0,:) =  tauoro(i)
         720             :         else where ( ( zi(i,:) < bwv(i) ) .and. ( zi(i,:) >= tlb(i) ) )
         721             :              tau(i,0,:) =  tauoro(i) +( taudsw(i)-tauoro(i) )* &
         722             :                                          ( bwv(i) - zi(i,:) ) / &
         723             :                                          ( bwv(i) - tlb(i) )
         724             :         endwhere
         725             :         ! low-level form drag on obstacle. Quantity kwvrdg (~1/b) appears for consistency
         726             :         ! with tauoro and taudsw forms. Should be weighted by L*b/A_g before applied to flow.
         727           0 :         where ( ( zi(i,:) < tlb(i) ) .and. ( zi(i,:) >= 0._r8 ) )
         728             :              tau(i,0,:) =  taudsw(i) +  &
         729             :                            Coeff_LB(i) * kwvrdg(i) * rsrc(i) * 0.5_r8 * (ubsrcx(i)**2) * ( tlb(i) - zi(i,:) )
         730             :         endwhere
         731             :
         732           0 :         if (do_smooth_regimes) then
         733             :         !  This blocks accounts for case where both mxdis and tlb fall
         734             :         !  between adjacent edges
         735           0 :            do k=1,pver
         736           0 :               if ( (zi(i,k) >= tlb(i)).and.(zi(i,k+1) < tlb(i)).and. &
         737           0 :                    (zi(i,k) >= mxdis(i)).and.(zi(i,k+1) < mxdis(i)) ) then
         738           0 :                  src_level(i) = src_level(i)-1
         739           0 :                  tau(i,0,k) = tauoro(i)
         740             :               end if
         741             :            end do
         742             :         end if
         743             :
         744             :      else     !----------------------------------------------
         745             :              ! This block allows low-level dynamics to occur
         746             :              ! even if m2 is less than orom2min
         747           0 :         where ( ( zi(i,:) < tlb(i) ) .and. ( zi(i,:) >= 0._r8 ) )
         748           0 :                tau(i,0,:) =  taudsw(i) +  &
         749             :                    Coeff_LB(i) * kwvrdg(i) * rsrc(i) * 0.5_r8 * &
         750             :                    (ubsrcx(i)**2) * ( tlb(i) - zi(i,:) )
         751             :         endwhere
         752             :      endif
         753             :   end do
         754             :
         755             :   ! This may be redundant with newest version of gw_drag_prof.
         756             :   ! That code reaches down to level k=src_level+1. (jtb 1/5/16)
         757           0 :   do i = 1, ncol
         758           0 :      k=src_level(i)
         759           0 :      if ( ni(i,k) > orostratmin ) then
         760             :          tausat    =  (Fr_c**2) * kwvrdg(i) * rhoi(i,k) * ubsrcx(i)**3 / &
         761           0 :               (1._r8*ni(i,k))
         762             :      else
         763             :          tausat = 0._r8
         764             :      endif
         765           0 :      tau(i,0,src_level(i)) = min( tauoro(i), tausat )
         766             :   end do
         767             :
         768             :
         769             :
         770             :   ! Final clean-up. Do nothing if obstacle less than orohmin
         771           0 :   do i = 1, ncol
         772           0 :      if ( mxdis(i) < orohmin ) then
         773           0 :         tau(i,0,:) = 0._r8
         774           0 :         tauoro(i)  = 0._r8
         775           0 :         taudsw(i)  = 0._r8
         776             :      endif
         777             :   end do
         778             :
         779             :           ! Disable vertical propagation if Scorer param is
         780             :           ! too small.
         781           0 :   do i = 1, ncol
         782           0 :      if ( m2src(i) <= orom2min ) then
         783           0 :         src_level(i)=1
         784             :      endif
         785             :   end do
         786             :
         787             :
         788             :
         789           0 : end subroutine gw_rdg_belowpeak
         790             :
         791             : !==========================================================================
         792           0 : subroutine gw_rdg_break_trap(ncol, band, &
         793           0 :      zi, nm, ni, ubm, ubi, rhoi, kwvrdg, bwv, tlb, wbr, &
         794           0 :      src_level, tlb_level, &
         795           0 :      hdspwv, hdspdw, mxdis, &
         796           0 :      tauoro, taudsw,  tau, &
         797             :      ldo_trapped_waves, wdth_kwv_scale_in )
         798             :   use gw_common, only: GWBand
         799             :   !-----------------------------------------------------------------------
         800             :   ! Parameterization of high-drag regimes and trapped lee-waves for CAM
         801             :   !
         802             :   !------------------------------Arguments--------------------------------
         803             :   ! Column dimension.
         804             :   integer, intent(in) :: ncol
         805             :   ! Band to emit orographic waves in.
         806             :   ! Regardless, we will only ever emit into l = 0.
         807             :   type(GWBand), intent(in) :: band
         808             :
         809             :
         810             :   ! Height estimate for ridge (m) [anisotropic orography].
         811             :   !real(r8), intent(in) :: mxdis(ncol)
         812             :   ! Horz wavenumber for ridge (1/m) [anisotropic orography].
         813             :   real(r8), intent(in) :: kwvrdg(ncol)
         814             :   ! Interface altitudes above ground (m).
         815             :   real(r8), intent(in) :: zi(ncol,pver+1)
         816             :   ! Midpoint Brunt-Vaisalla frequencies (s-1).
         817             :   real(r8), intent(in) :: nm(ncol,pver)
         818             :   ! Interface Brunt-Vaisalla frequencies (s-1).
         819             :   real(r8), intent(in) :: ni(ncol,pver+1)
         820             :
         821             :   ! Indices of gravity wave sources.
         822             :   integer, intent(inout) :: src_level(ncol), tlb_level(ncol)
         823             :
         824             :   ! Wave Reynolds stress.
         825             :   real(r8), intent(inout) :: tau(ncol,-band%ngwv:band%ngwv,pver+1)
         826             :   ! Wave Reynolds stresses at source.
         827             :   real(r8), intent(in)    :: taudsw(ncol)
         828             :   real(r8), intent(inout) :: tauoro(ncol)
         829             :   ! Projection of wind at midpoints and interfaces.
         830             :   real(r8), intent(in) :: ubm(ncol,pver)
         831             :   real(r8), intent(in) :: ubi(ncol,pver+1)
         832             :   ! Interface density (kg m-3).
         833             :   real(r8), intent(in) :: rhoi(ncol,pver+1)
         834             :
         835             :   ! Top of low-level flow layer.
         836             :   real(r8), intent(in) :: tlb(ncol)
         837             :   ! Bottom of linear wave region.
         838             :   real(r8), intent(in) :: bwv(ncol)
         839             :
         840             :   ! Surface streamline displacement height for linear waves.
         841             :   real(r8), intent(in) :: hdspwv(ncol)
         842             :   ! Surface streamline displacement height for downslope wind regime.
         843             :   real(r8), intent(in) :: hdspdw(ncol)
         844             :   ! Ridge height.
         845             :   real(r8), intent(in) :: mxdis(ncol)
         846             :
         847             :
         848             :   ! Wave breaking level
         849             :   real(r8), intent(out) :: wbr(ncol)
         850             :
         851             :   logical, intent(in), optional :: ldo_trapped_waves
         852             :   real(r8), intent(in), optional :: wdth_kwv_scale_in
         853             :
         854             :   !---------------------------Local Storage-------------------------------
         855             :   ! Column and level indices.
         856             :   integer :: i, k, kp1, non_hydro
         857           0 :   real(r8):: m2(ncol,pver),delz(ncol),tausat(ncol),trn(ncol)
         858           0 :   real(r8):: wbrx(ncol)
         859           0 :   real(r8):: phswkb(ncol,pver+1)
         860             :   logical :: lldo_trapped_waves
         861             :   real(r8):: wdth_kwv_scale
         862             :   ! Indices of important levels.
         863           0 :   integer :: trn_level(ncol)
         864             :
         865           0 :   if (present(ldo_trapped_waves)) then
         866           0 :      lldo_trapped_waves = ldo_trapped_waves
         867           0 :      if(lldo_trapped_waves) then
         868             :        non_hydro = 1
         869             :      else
         870           0 :        non_hydro = 0
         871             :      endif
         872             :   else
         873             :      lldo_trapped_waves = .false.
         874             :      non_hydro = 0
         875             :   endif
         876             :
         877           0 :   if (present(wdth_kwv_scale_in)) then
         878           0 :      wdth_kwv_scale = wdth_kwv_scale_in
         879             :   else
         880             :      wdth_kwv_scale = 1._r8
         881             :   endif
         882             :
         883             :   ! Calculate vertical wavenumber**2
         884             :   !---------------------------------
         885           0 :   m2 = (nm  / (abs(ubm)+.01_r8))**2
         886           0 :   do k=pver,1,-1
         887           0 :      m2(:,k) = m2(:,k) - non_hydro*(wdth_kwv_scale*kwvrdg)**2
         888             :      ! sweeping up, zero out m2 above first occurence
         889             :      ! of m2(:,k)<=0
         890           0 :      kp1=min( k+1, pver )
         891           0 :      where( (m2(:,k) <= 0.0_r8 ).or.(m2(:,kp1) <= 0.0_r8 ) )
         892             :         m2(:,k) = 0._r8
         893             :      endwhere
         894             :   end do
         895             :
         896             :   ! Take square root of m**2 and
         897             :   ! do vertical integral to find
         898             :   ! WKB phase.
         899             :   !-----------------------------
         900           0 :   m2 = SQRT( m2 )
         901           0 :   phswkb(:,:)=0
         902           0 :   do k=pver,1,-1
         903           0 :      where( zi(:,k) > tlb(:) )
         904           0 :         delz(:) = min( zi(:,k)-zi(:,k+1) , zi(:,k)-tlb(:) )
         905           0 :         phswkb(:,k) = phswkb(:,k+1) + m2(:,k)*delz(:)
         906             :      endwhere
         907             :   end do
         908             :
         909             :   ! Identify top edge of layer in which phswkb reaches 3*pi/2
         910             :   ! - approximately the "breaking level"
         911             :   !----------------------------------------------------------
         912           0 :   wbr(:)=0._r8
         913           0 :   wbrx(:)=0._r8
         914           0 :   if (do_smooth_regimes) then
         915           0 :      do k=pver,1,-1
         916           0 :      where( (phswkb(:,k+1)<1.5_r8*pii).and.(phswkb(:,k)>=1.5_r8*pii) &
         917           0 :             .and.(hdspdw(:)>hdspwv(:)) )
         918             :         wbr(:)  = zi(:,k)
         919             :         ! Extrapolation to make regime
         920             :         ! transitions smoother
         921             :         wbrx(:) = zi(:,k)   - ( phswkb(:,k) -  1.5_r8*pii ) &
         922           0 :                             / ( m2(:,k) + 1.e-6_r8 )
         923             :         src_level(:) = k-1
         924             :      endwhere
         925             :      end do
         926             :   else
         927           0 :      do k=pver,1,-1
         928           0 :      where( (phswkb(:,k+1)<1.5_r8*pii).and.(phswkb(:,k)>=1.5_r8*pii) &
         929           0 :             .and.(hdspdw(:)>hdspwv(:)) )
         930             :         wbr(:)  = zi(:,k)
         931             :         src_level(:) = k
         932             :      endwhere
         933             :      end do
         934             :   end if
         935             :
         936             :   ! Adjust tauoro at new source levels if needed.
         937             :   ! This is problematic if Fr_c<1.0. Not sure why.
         938             :   !----------------------------------------------------------
         939           0 :   if (do_adjust_tauoro) then
         940           0 :      do i = 1,ncol
         941           0 :         if (wbr(i) > 0._r8 ) then
         942           0 :             tausat(i) = (Fr_c**2) * kwvrdg(i)  * rhoi( i, src_level(i) ) &
         943             :                       * abs(ubi(i , src_level(i) ))**3  &
         944           0 :                       / ni( i , src_level(i) )
         945           0 :             tauoro(i) = min( tauoro(i), tausat(i) )
         946             :         end if
         947             :      end do
         948             :   end if
         949             :
         950           0 :   if (do_smooth_regimes) then
         951           0 :      do i = 1, ncol
         952           0 :      do k=1,pver+1
         953           0 :         if ( ( zi(i,k) <= wbr(i) ) .and. ( zi(i,k) > tlb(i) ) ) then
         954           0 :            tau(i,0,k) =  tauoro(i) + (taudsw(i)-tauoro(i)) * &
         955             :                           ( wbrx(i) - zi(i,k) ) / &
         956           0 :                           ( wbrx(i) - tlb(i)  )
         957           0 :            tau(i,0,k) = max( tau(i,0,k), tauoro(i) )
         958             :         endif
         959             :      end do
         960             :      end do
         961             :   else
         962             :   ! Following is for backwards B4B compatibility with earlier versions
         963             :   ! ("N1" and "N5" -- Note: "N5" used do_backward_compat=.true.)
         964           0 :      if (.not.do_backward_compat) then
         965           0 :         do i = 1, ncol
         966           0 :         do k=1,pver+1
         967           0 :            if ( ( zi(i,k) <  wbr(i) ) .and. ( zi(i,k) >= tlb(i) ) ) then
         968           0 :               tau(i,0,k) =  tauoro(i) + (taudsw(i)-tauoro(i)) * &
         969             :                             ( wbr(i) - zi(i,k) ) / &
         970           0 :                             ( wbr(i) - tlb(i)  )
         971             :            endif
         972             :         end do
         973             :         end do
         974             :      else
         975           0 :         do i = 1, ncol
         976           0 :         do k=1,pver+1
         977           0 :            if ( ( zi(i,k) <= wbr(i) ) .and. ( zi(i,k) > tlb(i) ) ) then
         978           0 :               tau(i,0,k) =  tauoro(i) + (taudsw(i)-tauoro(i)) * &
         979             :                             ( wbr(i) - zi(i,k) ) / &
         980           0 :                             ( wbr(i) - tlb(i)  )
         981             :            endif
         982             :         end do
         983             :         end do
         984             :      end if
         985             :   end if
         986             :
         987           0 :   if (lldo_trapped_waves) then
         988             :
         989             :   ! Identify top edge of layer in which Scorer param drops below 0
         990             :   ! - approximately the "turning level"
         991             :   !----------------------------------------------------------
         992           0 :      trn(:)=1.e8_r8
         993           0 :      trn_level(:) = 0 ! pver+1
         994           0 :      where( m2(:,pver)<= 0._r8 )
         995           0 :          trn(:) = zi(:,pver)
         996             :          trn_level(:) = pver
         997             :      endwhere
         998           0 :      do k=pver-1,1,-1
         999           0 :         where( (m2(:,k+1)> 0._r8).and.(m2(:,k)<= 0._r8) )
        1000           0 :            trn(:) = zi(:,k)
        1001             :            trn_level(:) = k
        1002             :         endwhere
        1003             :      end do
        1004             :
        1005           0 :      do i = 1,ncol
        1006             :      ! Case: Turning below mountain top
        1007           0 :         if ( (trn(i) < mxdis(i)).and.(trn_level(i)>=1) ) then
        1008           0 :             tau(i,0,:) =  tau(i,0,:) - max( tauoro(i),taudsw(i) )
        1009           0 :             tau(i,0,:) =  max( tau(i,0,:) , 0._r8 )
        1010           0 :             tau(i,0,1:tlb_level(i))=0._r8
        1011           0 :             src_level(i) = 1 ! disable any more tau calculation
        1012             :         end if
        1013             :         ! Case: Turning but no breaking
        1014           0 :         if ( (wbr(i) == 0._r8 ).and.(trn(i)>mxdis(i)).and.(trn_level(i)>=1) ) then
        1015           0 :            where ( ( zi(i,:) <= trn(i) ) .and. ( zi(i,:) >= bwv(i) ) )
        1016           0 :                tau(i,0,:) =  tauoro(i) * &
        1017             :                              ( trn(i) - zi(i,:) ) / &
        1018             :                              ( trn(i) - bwv(i)  )
        1019             :            end where
        1020           0 :            src_level(i) = 1 ! disable any more tau calculation
        1021             :         end if
        1022             :         ! Case: Turning AND breaking. Turning ABOVE breaking
        1023           0 :         if ( (wbr(i) > 0._r8 ).and.(trn(i) >= wbr(i)).and.(trn_level(i)>=1) ) then
        1024           0 :            where ( ( zi(i,:) <= trn(i) ) .and. ( zi(i,:) >= wbr(i) ) )
        1025           0 :                tau(i,0,:) =   tauoro(i) * &
        1026             :                              ( trn(i) - zi(i,:) ) / &
        1027             :                              ( trn(i) - wbr(i)  )
        1028             :            endwhere
        1029           0 :            src_level(i) = 1 ! disable any more tau calculation
        1030             :         end if
        1031             :         ! Case: Turning AND breaking. Turning BELOW breaking
        1032           0 :         if ( (wbr(i) > 0._r8 ).and.(trn(i) < wbr(i)).and.(trn_level(i)>=1) ) then
        1033           0 :            tauoro(i) = 0._r8
        1034           0 :            where ( ( zi(i,:) < wbr(i) ) .and. ( zi(i,:) >= tlb(i) ) )
        1035           0 :                tau(i,0,:) =  tauoro(i) + (taudsw(i)-tauoro(i)) * &
        1036             :                              ( wbr(i) - zi(i,:) ) / &
        1037             :                              ( wbr(i) - tlb(i)  )
        1038             :            endwhere
        1039           0 :            src_level(i) = 1 ! disable any more tau calculation
        1040             :         end if
        1041             :      end do
        1042             :   end if
        1043             :
        1044           0 :   end subroutine gw_rdg_break_trap
        1045             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        1046             :
        1047             :
        1048             :
        1049             : end module gw_rdg

Generated by: LCOV version 1.14