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