Line data    Source code

       1             : 
       2             : module rayleigh_friction
       3             : 
       4             : !---------------------------------------------------------------------------------
       5             : ! Module to apply rayleigh friction in region of model top.
       6             : ! We specify a decay rate profile that is largest at the model top and
       7             : ! drops off vertically using a hyperbolic tangent profile.
       8             : ! We compute the tendencies in u and v using an Euler backward scheme.
       9             : ! We then apply the negative of the kinetic energy tendency to "s", the dry
      10             : ! static energy.
      11             : !
      12             : ! calling sequence:
      13             : !
      14             : !  rayleigh_friction_init          initializes rayleigh friction constants
      15             : !  rayleigh_friction_tend          computes rayleigh friction tendencies
      16             : !
      17             : !---------------------------Code history--------------------------------
      18             : ! This is a new routine written by Art Mirin in collaboration with Phil Rasch.
      19             : ! Initial coding for this version:  Art Mirin, May 2007.
      20             : !---------------------------------------------------------------------------------
      21             : 
      22             : use shr_kind_mod,     only: r8 => shr_kind_r8
      23             : use ppgrid,           only: pver
      24             : use spmd_utils,       only: masterproc
      25             : use phys_control,     only: use_simple_phys
      26             : use cam_logfile,      only: iulog
      27             : use cam_abortutils,   only: endrun
      28             : 
      29             : implicit none
      30             : private
      31             : save
      32             :   
      33             : ! Public interfaces
      34             : public :: &
      35             :    rayleigh_friction_readnl, &! read namelist
      36             :    rayleigh_friction_init,   &! Initialization
      37             :    rayleigh_friction_tend     ! Computation of tendencies
      38             : 
      39             : ! Namelist variables
      40             : integer  :: rayk0 = 2           ! vertical level at which rayleigh friction term is centered
      41             : real(r8) :: raykrange = 0._r8   ! range of rayleigh friction profile 
      42             :                                 ! if 0, range is set to satisfy x=2 (see below)
      43             : real(r8) :: raytau0 = 0._r8     ! approximate value of decay time at model top (days)
      44             :                                 ! if 0., no rayleigh friction is applied
      45             : ! Local
      46             : real (r8) :: krange         ! range of rayleigh friction profile 
      47             : real (r8) :: tau0           ! approximate value of decay time at model top
      48             : real (r8) :: otau0          ! inverse of tau0
      49             : real (r8) :: otau(pver)     ! inverse decay time versus vertical level
      50             : 
      51             : ! We apply a profile of the form otau0 * [1 + tanh (x)] / 2 , where
      52             : ! x = (k0 - k) / krange. The default is for x to equal 2 at k=1, meaning
      53             : ! krange = (k0 - 1) / 2. The default is applied when raykrange is set to 0.
      54             : ! If otau0 = 0, no term is applied.
      55             : 
      56             : !===============================================================================
      57             : contains
      58             : !===============================================================================
      59             : 
      60     1490712 : subroutine rayleigh_friction_readnl(nlfile)
      61             : 
      62             :    use namelist_utils,  only: find_group_name
      63             :    use units,           only: getunit, freeunit
      64             :    use spmd_utils,      only: mpicom, mstrid=>masterprocid, mpi_integer, mpi_real8
      65             : 
      66             :    character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input
      67             : 
      68             :    ! Local variables
      69             :    integer :: unitn, ierr
      70             :    character(len=*), parameter :: sub = 'rayleigh_friction_readnl'
      71             : 
      72             :    namelist /rayleigh_friction_nl/ rayk0, raykrange, raytau0
      73             :    !-----------------------------------------------------------------------------
      74             : 
      75        1536 :    if (use_simple_phys) return
      76             : 
      77        1536 :    if (masterproc) then
      78           2 :       unitn = getunit()
      79           2 :       open( unitn, file=trim(nlfile), status='old' )
      80           2 :       call find_group_name(unitn, 'rayleigh_friction_nl', status=ierr)
      81           2 :       if (ierr == 0) then
      82           0 :          read(unitn, rayleigh_friction_nl, iostat=ierr)
      83           0 :          if (ierr /= 0) then
      84           0 :             call endrun(sub//': FATAL: reading namelist')
      85             :          end if
      86             :       end if
      87           2 :       close(unitn)
      88           2 :       call freeunit(unitn)
      89             :    end if
      90             : 
      91        1536 :    call mpi_bcast(rayk0, 1, mpi_integer, mstrid, mpicom, ierr)
      92        1536 :    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: rayk0")
      93        1536 :    call mpi_bcast(raykrange, 1, mpi_real8, mstrid, mpicom, ierr)
      94        1536 :    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: raykrange")
      95        1536 :    call mpi_bcast(raytau0, 1, mpi_real8, mstrid, mpicom, ierr)
      96        1536 :    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: raytau0")
      97             : 
      98        1536 :    if (masterproc) then
      99           2 :       if (raytau0 > 0._r8) then
     100           0 :          write (iulog,*) 'Rayleigh friction options: '
     101           0 :          write (iulog,*) '  rayk0     = ', rayk0
     102           0 :          write (iulog,*) '  raykrange = ', raykrange
     103           0 :          write (iulog,*) '  raytau0   = ', raytau0
     104             :       else
     105           2 :          write (iulog,*) 'Rayleigh friction not enabled.'
     106             :       end if
     107             :    end if
     108             : 
     109             : end subroutine rayleigh_friction_readnl
     110             : 
     111             : !===============================================================================
     112             : 
     113        1536 : subroutine rayleigh_friction_init()
     114             : 
     115             :    !---------------------------Local storage-------------------------------
     116             :    real (r8) x
     117             :    integer k
     118             : 
     119             :    !-----------------------------------------------------------------------
     120             :    ! Compute tau array
     121             :    !-----------------------------------------------------------------------
     122             : 
     123        1536 :    krange = raykrange
     124        1536 :    if (raykrange .eq. 0._r8) krange = (rayk0 - 1) / 2._r8
     125             : 
     126        1536 :    tau0 = (86400._r8) * raytau0   ! convert to seconds
     127        1536 :    otau0 = 0._r8
     128        1536 :    if (tau0 .ne. 0._r8) otau0 = 1._r8/tau0
     129             : 
     130      144384 :    do k = 1, pver
     131      142848 :       x = (rayk0 - k) / krange
     132      144384 :       otau(k) = otau0 * (1 + tanh(x)) / (2._r8)
     133             :    enddo
     134             : 
     135        1536 :    if (masterproc) then
     136           2 :       if (tau0 > 0._r8) then
     137           0 :          write (iulog,*) 'Rayleigh friction - krange = ', krange
     138           0 :          write (iulog,*) 'Rayleigh friction - otau0 = ', otau0
     139           0 :          write (iulog,*) 'Rayleigh friction decay rate profile'
     140           0 :          do k = 1, pver
     141           0 :             write (iulog,*) '   k = ', k, '   otau = ', otau(k)
     142             :          enddo
     143             :       end if
     144             :    end if
     145             : 
     146        1536 : end subroutine rayleigh_friction_init
     147             :   
     148             : !=========================================================================================
     149             : 
     150    62545392 : subroutine rayleigh_friction_tend(                                     &
     151             :    ztodt    ,state    ,ptend    )
     152             : 
     153             :    !-----------------------------------------------------------------------
     154             :    ! compute tendencies for rayleigh friction
     155             :    !-----------------------------------------------------------------------
     156             :    use physics_types, only: physics_state, physics_ptend, physics_ptend_init
     157             : 
     158             :    !------------------------------Arguments--------------------------------
     159             :    real(r8),            intent(in) :: ztodt      ! physics timestep
     160             :    type(physics_state), intent(in) :: state      ! physics state variables
     161             :     
     162             :    type(physics_ptend), intent(out):: ptend      ! individual parameterization tendencies
     163             : 
     164             :    !---------------------------Local storage-------------------------------
     165             :    integer :: ncol                                ! number of atmospheric columns
     166             :    integer :: k                                   ! level
     167             :    real(r8) :: rztodt                             ! 1./ztodt
     168             :    real(r8) :: c1, c2, c3                         ! temporary variables
     169             :    !-----------------------------------------------------------------------
     170             : 
     171     1489176 :    call physics_ptend_init(ptend, state%psetcols, 'rayleigh friction', ls=.true., lu=.true., lv=.true.)
     172             : 
     173     1489176 :    if (otau0 .eq. 0._r8) return
     174             : 
     175           0 :    rztodt = 1._r8/ztodt
     176           0 :    ncol  = state%ncol
     177             : 
     178             :    ! u, v and s are modified by rayleigh friction
     179             : 
     180           0 :    do k = 1, pver
     181           0 :       c2 = 1._r8 / (1._r8 + otau(k)*ztodt)
     182           0 :       c1 = -otau(k) * c2
     183           0 :       c3 = 0.5_r8 * (1._r8 - c2*c2) * rztodt
     184           0 :       ptend%u(:ncol,k) = c1 * state%u(:ncol,k)
     185           0 :       ptend%v(:ncol,k) = c1 * state%v(:ncol,k)
     186           0 :       ptend%s(:ncol,k) = c3 * (state%u(:ncol,k)**2 + state%v(:ncol,k)**2)
     187             :    enddo
     188             : 
     189     1489176 : end subroutine rayleigh_friction_tend
     190             : 
     191             : end module rayleigh_friction