LCOV - code coverage report
Current view: | top level - physics/cam - rayleigh_friction.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 34 | 55 | 61.8 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 3 | 3 | 100.0 %

              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

Generated by: LCOV version 1.14