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

              Line data    Source code
           1             : module trb_mtn_stress
           2             :
           3             :   implicit none
           4             :   private
           5             :   save
           6             :
           7             :   public init_tms                             ! Initialization
           8             :   public compute_tms                          ! Full routine
           9             :
          10             :   ! ------------ !
          11             :   ! Private data !
          12             :   ! ------------ !
          13             :
          14             :   integer,  parameter :: r8 = selected_real_kind(12) ! 8 byte real
          15             :
          16             :   real(r8), parameter :: horomin= 1._r8       ! Minimum value of subgrid orographic height for mountain stress [ m ]
          17             :   real(r8), parameter :: z0max  = 100._r8     ! Maximum value of z_0 for orography [ m ]
          18             :   real(r8), parameter :: dv2min = 0.01_r8     ! Minimum shear squared [ m2/s2 ]
          19             :   real(r8)            :: orocnst              ! Converts from standard deviation to height [ no unit ]
          20             :   real(r8)            :: z0fac                ! Factor determining z_0 from orographic standard deviation [ no unit ]
          21             :   real(r8)            :: karman               ! von Karman constant
          22             :   real(r8)            :: gravit               ! Acceleration due to gravity
          23             :   real(r8)            :: rair                 ! Gas constant for dry air
          24             :
          25             : contains
          26             :
          27             :   !============================================================================ !
          28             :   !                                                                             !
          29             :   !============================================================================ !
          30             :
          31           0 :   subroutine init_tms( kind, oro_in, z0fac_in, karman_in, gravit_in, rair_in, &
          32             :        errstring)
          33             :
          34             :     integer, intent(in) :: kind
          35             :
          36             :     real(r8), intent(in) :: oro_in, z0fac_in, karman_in, gravit_in, rair_in
          37             :
          38             :     character(len=*), intent(out) :: errstring
          39             :
          40           0 :     errstring = ' '
          41             :
          42           0 :     if ( kind /= r8 ) then
          43           0 :        errstring = 'inconsistent KIND of reals passed to init_tms'
          44           0 :        return
          45             :     endif
          46             :
          47           0 :     orocnst  = oro_in
          48           0 :     z0fac    = z0fac_in
          49           0 :     karman   = karman_in
          50           0 :     gravit   = gravit_in
          51           0 :     rair     = rair_in
          52             :
          53           0 :   end subroutine init_tms
          54             :
          55             :   !============================================================================ !
          56             :   !                                                                             !
          57             :   !============================================================================ !
          58             :
          59           0 :   subroutine compute_tms( pcols    , pver    , ncol    ,                     &
          60           0 :                           u        , v       , t       , pmid    , exner   , &
          61           0 :                           zm       , sgh     , ksrf    , taux    , tauy    , &
          62           0 :                           landfrac )
          63             :
          64             :     !------------------------------------------------------------------------------ !
          65             :     ! Turbulent mountain stress parameterization                                    !
          66             :     !                                                                               !
          67             :     ! Returns surface drag coefficient and stress associated with subgrid mountains !
          68             :     ! For points where the orographic variance is small ( including ocean ),        !
          69             :     ! the returned surface drag coefficient and stress is zero.                     !
          70             :     !                                                                               !
          71             :     ! Lastly arranged : Sungsu Park. Jan. 2010.                                     !
          72             :     !------------------------------------------------------------------------------ !
          73             :
          74             :     ! ---------------------- !
          75             :     ! Input-Output Arguments !
          76             :     ! ---------------------- !
          77             :
          78             :     integer,  intent(in)  :: pcols                 ! Number of columns dimensioned
          79             :     integer,  intent(in)  :: pver                  ! Number of model layers
          80             :     integer,  intent(in)  :: ncol                  ! Number of columns actually used
          81             :
          82             :     real(r8), intent(in)  :: u(pcols,pver)         ! Layer mid-point zonal wind [ m/s ]
          83             :     real(r8), intent(in)  :: v(pcols,pver)         ! Layer mid-point meridional wind [ m/s ]
          84             :     real(r8), intent(in)  :: t(pcols,pver)         ! Layer mid-point temperature [ K ]
          85             :     real(r8), intent(in)  :: pmid(pcols,pver)      ! Layer mid-point pressure [ Pa ]
          86             :     real(r8), intent(in)  :: exner(pcols,pver)     ! Layer mid-point exner function [ no unit ]
          87             :     real(r8), intent(in)  :: zm(pcols,pver)        ! Layer mid-point height [ m ]
          88             :     real(r8), intent(in)  :: sgh(pcols)            ! Standard deviation of orography [ m ]
          89             :     real(r8), intent(in)  :: landfrac(pcols)       ! Land fraction [ fraction ]
          90             :
          91             :     real(r8), intent(out) :: ksrf(pcols)           ! Surface drag coefficient [ kg/s/m2 ]
          92             :     real(r8), intent(out) :: taux(pcols)           ! Surface zonal      wind stress [ N/m2 ]
          93             :     real(r8), intent(out) :: tauy(pcols)           ! Surface meridional wind stress [ N/m2 ]
          94             :
          95             :     ! --------------- !
          96             :     ! Local Variables !
          97             :     ! --------------- !
          98             :
          99             :     integer  :: i                                  ! Loop index
         100             :     integer  :: kb, kt                             ! Bottom and top of source region
         101             :
         102             :     real(r8) :: horo                               ! Orographic height [ m ]
         103             :     real(r8) :: z0oro                              ! Orographic z0 for momentum [ m ]
         104             :     real(r8) :: dv2                                ! (delta v)**2 [ m2/s2 ]
         105             :     real(r8) :: ri                                 ! Richardson number [ no unit ]
         106             :     real(r8) :: stabfri                            ! Instability function of Richardson number [ no unit ]
         107             :     real(r8) :: rho                                ! Density [ kg/m3 ]
         108             :     real(r8) :: cd                                 ! Drag coefficient [ no unit ]
         109             :     real(r8) :: vmag                               ! Velocity magnitude [ m /s ]
         110             :
         111             :     ! ----------------------- !
         112             :     ! Main Computation Begins !
         113             :     ! ----------------------- !
         114             :
         115           0 :     do i = 1, ncol
         116             :
         117             :      ! determine subgrid orgraphic height ( mean to peak )
         118             :
         119           0 :        horo = orocnst * sgh(i)
         120             :
         121             :      ! No mountain stress if horo is too small
         122             :
         123           0 :        if( horo < horomin ) then
         124             :
         125           0 :            ksrf(i) = 0._r8
         126           0 :            taux(i) = 0._r8
         127           0 :            tauy(i) = 0._r8
         128             :
         129             :        else
         130             :
         131             :          ! Determine z0m for orography
         132             :
         133           0 :            z0oro = min( z0fac * horo, z0max )
         134             :
         135             :          ! Calculate neutral drag coefficient
         136             :
         137           0 :            cd = ( karman / log( ( zm(i,pver) + z0oro ) / z0oro) )**2
         138             :
         139             :          ! Calculate the Richardson number over the lowest 2 layers
         140             :
         141           0 :            kt  = pver - 1
         142           0 :            kb  = pver
         143           0 :            dv2 = max( ( u(i,kt) - u(i,kb) )**2 + ( v(i,kt) - v(i,kb) )**2, dv2min )
         144             :
         145             :          ! Modification : Below computation of Ri is wrong. Note that 'Exner' function here is
         146             :          !                inverse exner function. Here, exner function is not multiplied in
         147             :          !                the denominator. Also, we should use moist Ri not dry Ri.
         148             :          !                Also, this approach using the two lowest model layers can be potentially
         149             :          !                sensitive to the vertical resolution.
         150             :          ! OK. I only modified the part associated with exner function.
         151             :
         152             :            ri  = 2._r8 * gravit * ( t(i,kt) * exner(i,kt) - t(i,kb) * exner(i,kb) ) * ( zm(i,kt) - zm(i,kb) ) &
         153           0 :                                 / ( ( t(i,kt) * exner(i,kt) + t(i,kb) * exner(i,kb) ) * dv2 )
         154             :
         155             :          ! ri  = 2._r8 * gravit * ( t(i,kt) * exner(i,kt) - t(i,kb) * exner(i,kb) ) * ( zm(i,kt) - zm(i,kb) ) &
         156             :          !                      / ( ( t(i,kt) + t(i,kb) ) * dv2 )
         157             :
         158             :          ! Calculate the instability function and modify the neutral drag cofficient.
         159             :          ! We should probably follow more elegant approach like Louis et al (1982) or Bretherton and Park (2009)
         160             :          ! but for now we use very crude approach : just 1 for ri < 0, 0 for ri > 1, and linear ramping.
         161             :
         162           0 :            stabfri = max( 0._r8, min( 1._r8, 1._r8 - ri ) )
         163           0 :            cd      = cd * stabfri
         164             :
         165             :          ! Compute density, velocity magnitude and stress using bottom level properties
         166             :
         167           0 :            rho     = pmid(i,pver) / ( rair * t(i,pver) )
         168           0 :            vmag    = sqrt( u(i,pver)**2 + v(i,pver)**2 )
         169           0 :            ksrf(i) = rho * cd * vmag * landfrac(i)
         170           0 :            taux(i) = -ksrf(i) * u(i,pver)
         171           0 :            tauy(i) = -ksrf(i) * v(i,pver)
         172             :
         173             :        end if
         174             :
         175             :     end do
         176             :
         177           0 :     return
         178             :   end subroutine compute_tms
         179             :
         180             : end module trb_mtn_stress

Generated by: LCOV version 1.14