diff --git a/schemes/rayleigh_friction/rayleigh_friction.F90 b/schemes/rayleigh_friction/rayleigh_friction.F90 new file mode 100644 index 00000000..6d7a543b --- /dev/null +++ b/schemes/rayleigh_friction/rayleigh_friction.F90 @@ -0,0 +1,160 @@ +module rayleigh_friction + +!--------------------------------------------------------------------------------- +! Module to apply rayleigh friction in region of model top. +! We specify a decay rate profile that is largest at the model top and +! drops off vertically using a hyperbolic tangent profile. +! We compute the tendencies in u and v using an Euler backward scheme. +! We then apply the negative of the kinetic energy tendency to "s", the dry +! static energy. +!---------------------------Code history-------------------------------- +! This is a new routine written by Art Mirin in collaboration with Phil Rasch. +! Initial coding for this version: Art Mirin, May 2007. +! First CCPP conversion by Andrew Gettelman, Sep 2021 +! Final CCPP completion by Kate T-C, Oct 2024 +!--------------------------------------------------------------------------------- + + use ccpp_kinds, only: kind_phys + + implicit none + + save + private ! Make default type private to the module +! +! PUBLIC: interfaces +! + public rayleigh_friction_init + public rayleigh_friction_run + +! PRIVATE: module variables +integer :: rayk0 = 2 ! vertical level at which rayleigh friction term is centered +real(kind_phys) :: raykrange = 0._kind_phys ! range of rayleigh friction profile + ! if 0, range is set to satisfy x=2 (see below) +real(kind_phys) :: raytau0 = 0._kind_phys ! approximate value of decay time at model top (days) + ! if 0., no rayleigh friction is applied +! Local +real (kind_phys) :: krange ! range of rayleigh friction profile +real (kind_phys) :: tau0 ! approximate value of decay time at model top +real (kind_phys) :: otau0 ! inverse of tau0 +real (kind_phys), allocatable :: otau(:) ! inverse decay time versus vertical level + +! We apply a profile of the form otau0 * [1 + tanh (x)] / 2 , where +! x = (k0 - k) / krange. The default is for x to equal 2 at k=1, meaning +! krange = (k0 - 1) / 2. The default is applied when raykrange is set to 0. +! If otau0 = 0, no term is applied. + +contains + + !> \section arg_table_rayleigh_friction_init Argument Table + !! \htmlinclude rayleigh_friction_init.html + subroutine rayleigh_friction_init(pver, raytau0_nl, raykrange_nl, rayk0_nl, masterproc, iulog, errmsg, errflg) + + integer, intent(in) :: pver + real (kind_phys), intent(in) :: raytau0_nl + real (kind_phys), intent(in) :: raykrange_nl + integer, intent(in) :: rayk0_nl + logical, intent(in) :: masterproc + integer, intent(in) :: iulog + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! local variables + character(len=*), parameter :: subname = 'rayleigh_friction_init' + real (kind_phys) x + integer k, ierr + + errmsg = '' + errflg = 0 + raytau0 = raytau0_nl + raykrange = raykrange_nl + rayk0 = rayk0_nl + + if ( raytau0 > 0._kind_phys) then ! Rayleigh Friction enabled + ! Allocate module data + allocate(otau(pver), stat=ierr) + if (ierr /= 0) then + errflg = ierr + errmsg = trim(subname)//': Allocate of otau failed' + return + end if + + !----------------------------------------------------------------------- + ! Compute tau array + !----------------------------------------------------------------------- + + krange = raykrange + if (raykrange .eq. 0._kind_phys) krange = (rayk0 - 1) / 2._kind_phys + + tau0 = (86400._kind_phys) * raytau0 ! convert to seconds + otau0 = 0._kind_phys + if (tau0 .ne. 0._kind_phys) otau0 = 1._kind_phys/tau0 + + do k = 1, pver + x = (rayk0 - k) / krange + otau(k) = otau0 * (1 + tanh(x)) / (2._kind_phys) + enddo + end if ! Rayleigh Friction enabled + + if ( masterproc ) then + if (raytau0 > 0._kind_phys) then + write(iulog,*) 'tuning parameters rayleigh_friction_init: rayk0',rayk0 + write(iulog,*) 'tuning parameters rayleigh_friction_init: raykrange',raykrange + write(iulog,*) 'tuning parameters rayleigh_friction_init: raytau0',raytau0 + else + write(iulog,*) 'Rayleigh Friction not enabled.' + endif + endif + + end subroutine rayleigh_friction_init + + !========================================================================================= + !> \section arg_table_rayleigh_friction_run Argument Table + !! \htmlinclude rayleigh_friction_run.html + subroutine rayleigh_friction_run(ncol, pver, ztodt, u, v, dudt, dvdt, dsdt, errmsg, errflg) + + !------------------------------Arguments-------------------------------- + integer, intent(in) :: ncol + integer, intent(in) :: pver + real(kind_phys), intent(in) :: ztodt !physics timestep + real(kind_phys), intent(in) :: u(:,:) + real(kind_phys), intent(in) :: v(:,:) + real(kind_phys), intent(out) :: dudt(:,:) !tendency_of_eastward_wind + real(kind_phys), intent(out) :: dvdt(:,:) !tendency_of_northward_wind + real(kind_phys), intent(out) :: dsdt(:,:) !heating_rate + + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errflg + + !---------------------------Local storage------------------------------- + character(len=*), parameter :: subname = 'rayleigh_friction_run' + integer :: k ! level + real(kind_phys) :: rztodt ! 1./ztodt + real(kind_phys) :: c1, c2, c3 ! temporary variables + !----------------------------------------------------------------------- + + ! initialize values + errmsg = '' + errflg = 0 + dudt(:,:)=0._kind_phys + dvdt(:,:)=0._kind_phys + dsdt(:,:) =0._kind_phys + + if (raytau0 .eq. 0._kind_phys) return + + rztodt = 1._kind_phys/ztodt + + ! u, v and s are modified by rayleigh friction + + do k = 1, pver + c2 = 1._kind_phys / (1._kind_phys + otau(k)*ztodt) + c1 = -otau(k) * c2 + c3 = 0.5_kind_phys * (1._kind_phys - c2*c2) * rztodt + dudt(:ncol,k) = c1 * u(:ncol,k) + dvdt(:ncol,k) = c1 * v(:ncol,k) + dsdt(:ncol,k) = c3 * (u(:ncol,k)**2 + v(:ncol,k)**2) + enddo + + end subroutine rayleigh_friction_run + +end module rayleigh_friction + diff --git a/schemes/rayleigh_friction/rayleigh_friction.meta b/schemes/rayleigh_friction/rayleigh_friction.meta new file mode 100644 index 00000000..7c9eae12 --- /dev/null +++ b/schemes/rayleigh_friction/rayleigh_friction.meta @@ -0,0 +1,126 @@ +[ccpp-table-properties] + name = rayleigh_friction + type = scheme + +[ccpp-arg-table] + name = rayleigh_friction_init + type = scheme +[ pver ] + standard_name = vertical_layer_dimension + units = count + type = integer + dimensions = () + intent = in +[ raytau0_nl ] + standard_name = model_top_decay_time_for_rayleigh_friction + long_name = Approx decay time at model top for RF + units = days + type = real | kind = kind_phys + dimensions = () + intent = in +[ raykrange_nl ] + standard_name = number_of_vertical_layers_for_rayleigh_friction + long_name = Range of levels with RF applied + units = count + type = real | kind = kind_phys + dimensions = () + intent = in +[ rayk0_nl ] + standard_name = center_vertical_layer_for_rayleigh_friction + long_name = Vertical level at which RF profile is centered + units = index + type = integer + dimensions = () + intent = in +[ masterproc ] + standard_name = flag_for_mpi_root + units = flag + type = logical + dimensions = () + intent = in +[ iulog ] + standard_name = log_output_unit + units = 1 + type = integer + dimensions = () + intent = in +[ errmsg ] + standard_name = ccpp_error_message + long_name = Error message for error handling in CCPP + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + long_name = Error flag for error handling in CCPP + units = 1 + type = integer + dimensions = () + intent = out + +[ccpp-arg-table] + name = rayleigh_friction_run + type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ pver ] + standard_name = vertical_layer_dimension + units = count + type = integer + dimensions = () + intent = in +[ ztodt ] + standard_name = timestep_for_physics + units = s + type = real | kind = kind_phys + dimensions = () + intent = in +[ u ] + standard_name = eastward_wind + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ v ] + standard_name = northward_wind + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ dudt ] + standard_name = tendency_of_eastward_wind + units = m s-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ dvdt ] + standard_name = tendency_of_northward_wind + units = m s-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ dsdt ] + standard_name = tendency_of_dry_air_enthalpy_at_constant_pressure + units = J kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ errmsg ] + standard_name = ccpp_error_message + long_name = Error message for error handling in CCPP + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + long_name = Error flag for error handling in CCPP + units = 1 + type = integer + dimensions = () + intent = out diff --git a/schemes/rayleigh_friction/rayleigh_friction_namelist.xml b/schemes/rayleigh_friction/rayleigh_friction_namelist.xml new file mode 100644 index 00000000..279e344f --- /dev/null +++ b/schemes/rayleigh_friction/rayleigh_friction_namelist.xml @@ -0,0 +1,56 @@ + + + + + + + + + + + integer + rayleigh_friction + rayleigh_friction_nl + center_vertical_layer_for_rayleigh_friction + index + + Variable to specify the vertical index at which the + Rayleigh friction term is centered (the peak value). + + + 2 + + + + + real + kind_phys + rayleigh_friction + rayleigh_friction_nl + number_of_vertical_layers_for_rayleigh_friction + count + + Range (number of layers) of the Rayleigh Friction profile. If 0, range is set to satisfy center=2. + + + 0 + + + + + real + kind_phys + rayleigh_friction + rayleigh_friction_nl + model_top_decay_time_for_rayleigh_friction + days + + Approximate value of max decay time (in units of days) at model top. If 0, no Rayleigh Friction is computed or applied. + + + 0 + + + + + diff --git a/schemes/sima_diagnostics/rayleigh_friction_diagnostics.F90 b/schemes/sima_diagnostics/rayleigh_friction_diagnostics.F90 new file mode 100644 index 00000000..1b4ae55d --- /dev/null +++ b/schemes/sima_diagnostics/rayleigh_friction_diagnostics.F90 @@ -0,0 +1,64 @@ +module rayleigh_friction_diagnostics + + use ccpp_kinds, only: kind_phys + + implicit none + private + save + + public :: rayleigh_friction_diagnostics_init ! init routine + public :: rayleigh_friction_diagnostics_run ! main routine + +CONTAINS + + !> \section arg_table_rayleigh_friction_diagnostics_init Argument Table + !! \htmlinclude rayleigh_friction_diagnostics_init.html + subroutine rayleigh_friction_diagnostics_init(errmsg, errflg) + use cam_history, only: history_add_field + use cam_history_support, only: horiz_only + + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! Local variables: + + errmsg = '' + errflg = 0 + + ! History add field calls + call history_add_field('UTEND_RAYLEIGH', 'Zonal wind tendency due to Rayleigh Friction', 'lev', 'avg', 'm s-2') + call history_add_field('VTEND_RAYLEIGH', 'Meridional wind tendency due to Rayleigh Friction', 'lev', 'avg', 'm s-2') + call history_add_field('STEND_RAYLEIGH', 'Dry air enthalpy tendency due to Rayleigh Friction', 'lev', 'avg', 'J kg-1') + + end subroutine rayleigh_friction_diagnostics_init + + !> \section arg_table_rayleigh_friction_diagnostics_run Argument Table + !! \htmlinclude rayleigh_friction_diagnostics_run.html + subroutine rayleigh_friction_diagnostics_run(dudt, dvdt, dsdt, errmsg, errflg) + + use cam_history, only: history_out_field + !------------------------------------------------ + ! Input / output parameters + !------------------------------------------------ + ! State variables + real(kind_phys), intent(in) :: dudt(:,:) !tendency_of_eastward_wind due to RF + real(kind_phys), intent(in) :: dvdt(:,:) !tendency_of_northward_wind due to RF + real(kind_phys), intent(in) :: dsdt(:,:) !tendency_of_dry_air_enthalpy_at_constant_pressure due to RF + + ! CCPP error handling variables + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errflg + + errmsg = '' + errflg = 0 + + ! History out field calls + call history_out_field('UTEND_RAYLEIGH', dudt) + call history_out_field('VTEND_RAYLEIGH', dvdt) + call history_out_field('STEND_RAYLEIGH', dsdt) + + end subroutine rayleigh_friction_diagnostics_run + + !======================================================================= + +end module rayleigh_friction_diagnostics diff --git a/schemes/sima_diagnostics/rayleigh_friction_diagnostics.meta b/schemes/sima_diagnostics/rayleigh_friction_diagnostics.meta new file mode 100644 index 00000000..b481b852 --- /dev/null +++ b/schemes/sima_diagnostics/rayleigh_friction_diagnostics.meta @@ -0,0 +1,53 @@ +[ccpp-table-properties] + name = rayleigh_friction_diagnostics + type = scheme + +[ccpp-arg-table] + name = rayleigh_friction_diagnostics_init + type = scheme +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out + +[ccpp-arg-table] + name = rayleigh_friction_diagnostics_run + type = scheme +[ dudt ] + standard_name = tendency_of_eastward_wind + units = m s-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ dvdt ] + standard_name = tendency_of_northward_wind + units = m s-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ dsdt ] + standard_name = tendency_of_dry_air_enthalpy_at_constant_pressure + units = J kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out diff --git a/test/test_suites/suite_rayleigh_friction.xml b/test/test_suites/suite_rayleigh_friction.xml new file mode 100644 index 00000000..3dc15696 --- /dev/null +++ b/test/test_suites/suite_rayleigh_friction.xml @@ -0,0 +1,14 @@ + + + + + + rayleigh_friction + rayleigh_friction_diagnostics + apply_heating_rate + apply_tendency_of_eastward_wind + apply_tendency_of_northward_wind + geopotential_temp + + +