phys/module_cam_trb_mtn_stress.F

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