Complete CCPPization of Beljaars orographic form drag scheme

schemes/beljaars_drag/beljaars_drag.F90

@@ -0,0 +1,101 @@
+! Beljaars Sub-Grid Orographic (SGO) Form Drag Parameterization
+! Returns drag profile and integrated stress associated with subgrid mountains
+! with horizontal length scales nominally below 3 km.
+! Based on Beljaars et al. (2004, QJRMS) https://doi.org/10.1256/qj.03.73.
+!
+! Original author: J. Bacmeister, March 2016, based on TMS
+module beljaars_drag
+  implicit none
+  private
+
+  public :: beljaars_drag_run
+
+contains
+!> \section arg_table_beljaars_drag_run Argument Table
+!! \htmlinclude beljaars_drag_run.html
+  ! Compute Beljaars SGO form drag profile and surface stresses
+  pure subroutine beljaars_drag_run( &
+    do_beljaars, &
+    ncol, pver, &
+    u, v, delp, zm, sgh30, &
+    gravit, &
+    ! below output:
+    drag, taux, tauy, &
+    errmsg, errflg)
+    use ccpp_kinds, only: kind_phys
+
+    logical,            intent(in)  :: do_beljaars  ! is Beljaars active?
+    integer,            intent(in)  :: ncol
+    integer,            intent(in)  :: pver
+    real(kind_phys),    intent(in)  :: u(:, :)      ! zonal wind [m s-1]
+    real(kind_phys),    intent(in)  :: v(:, :)      ! meridional wind [m s-1]
+    real(kind_phys),    intent(in)  :: delp(:, :)   ! air pressure thickness [Pa]
+    real(kind_phys),    intent(in)  :: zm(:, :)     ! geopotential height wrt surface [m]
+    real(kind_phys),    intent(in)  :: sgh30(:)     ! standard deviation of subgrid orography [m]
+    real(kind_phys),    intent(in)  :: gravit       ! gravitational acceleration [m s-2]
+
+    real(kind_phys),    intent(out) :: drag(:, :)   ! SGO drag profile [s-1]
+    real(kind_phys),    intent(out) :: taux(:)      ! surface zonal wind stress [N m-2]
+    real(kind_phys),    intent(out) :: tauy(:)      ! surface meridional wind stress [N m-2]
+    character(len=*),   intent(out) :: errmsg
+    integer,            intent(out) :: errflg
+
+    ! Local variables
+    integer  :: i, k
+
+    real(kind_phys) :: vmag                         ! velocity magnitude [m s-1]
+
+    real(kind_phys), parameter :: alpha = 12.0_kind_phys
+    real(kind_phys), parameter :: beta  = 1.0_kind_phys
+    real(kind_phys), parameter :: n1    = -1.9_kind_phys
+    real(kind_phys), parameter :: n2    = -2.8_kind_phys
+    real(kind_phys), parameter :: Cmd   = 0.005_kind_phys
+    real(kind_phys), parameter :: Ccorr = 0.6_kind_phys * 5.0_kind_phys
+    real(kind_phys), parameter :: kflt  = 0.00035_kind_phys   ! m-1
+    real(kind_phys), parameter :: k1    = 0.003_kind_phys     ! m-1
+    real(kind_phys), parameter :: IH    = 0.00102_kind_phys   ! m-1
+    real(kind_phys) :: a1(ncol), a2(ncol)
+
+    errmsg = ''
+    errflg = 0
+
+    if(.not. do_beljaars) then
+      ! if not doing Beljaars, return zero drag and stresses from routine.
+      drag(:,:) = 0._kind_phys
+      taux(:)   = 0._kind_phys
+      tauy(:)   = 0._kind_phys
+      return
+    end if
+
+    a1(1:ncol) = (sgh30(1:ncol) * sgh30(1:ncol)) / (IH * (kflt**n1))
+    a2(1:ncol) = a1(1:ncol) * k1**(n1 - n2)
+
+    do k = 1, pver
+      do i = 1, ncol
+        vmag      = sqrt(u(i, k)**2 + v(i, k)**2)
+        drag(i, k) = -alpha * beta * Cmd * Ccorr * vmag * 2.109_kind_phys * &
+                     exp(-(zm(i, k) / 1500.0_kind_phys) * sqrt(zm(i, k) / 1500.0_kind_phys)) * &
+                     (zm(i, k)**(-1.2_kind_phys)) * a2(i)
+      end do
+    end do
+
+    ! Diagnose effective surface drag in X and Y by integrating in the vertical
+    !
+    ! taux, tauy stresses generated by Beljaars drag here
+    ! uses the pre-vertical diffusion winds (and not the updated winds)
+    ! however, at this point only these winds are available because the diffusion
+    ! solver runs after the orographic drag.
+    ! the actual provisionally updated winds are actually used to recompute taubljx/y
+    ! which are added to the updated residual stress. see the physics scheme
+    ! beljaars_add_updated_residual_stress.
+    taux(:) = 0.0_kind_phys
+    tauy(:) = 0.0_kind_phys
+    do k = 1, pver
+      do i = 1, ncol
+        taux(i) = taux(i) + drag(i, k) * u(i, k) * delp(i, k) / gravit
+        tauy(i) = tauy(i) + drag(i, k) * v(i, k) * delp(i, k) / gravit
+      end do
+    end do
+
+  end subroutine beljaars_drag_run
+end module beljaars_drag

schemes/beljaars_drag/beljaars_drag.meta

@@ -0,0 +1,91 @@
+[ccpp-table-properties]
+  name = beljaars_drag
+  type = scheme
+
+[ccpp-arg-table]
+  name  = beljaars_drag_run
+  type  = scheme
+[ do_beljaars ]
+  standard_name = do_beljaars_drag_in_vertical_diffusion
+  units = flag
+  type = logical
+  dimensions = ()
+  intent = in
+[ 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
+[ 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
+[ delp ]
+  standard_name = air_pressure_thickness
+  units = Pa
+  type = real | kind = kind_phys
+  dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+  intent = in
+[ zm ]
+  standard_name = geopotential_height_wrt_surface
+  units = m
+  type = real | kind = kind_phys
+  dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+  intent = in
+[ sgh30 ]
+  standard_name = standard_deviation_of_subgrid_orography_for_turbulent_orographic_form_drag
+  units = m
+  type = real | kind = kind_phys
+  dimensions = (horizontal_loop_extent)
+  intent = in
+[ gravit ]
+  standard_name = standard_gravitational_acceleration
+  units = m s-2
+  type = real | kind = kind_phys
+  dimensions = ()
+  intent = in
+[ drag ]
+  standard_name = turbulent_orographic_form_drag_coefficent
+  units = s-1
+  type = real | kind = kind_phys
+  dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+  intent = out
+[ taux ]
+  standard_name = eastward_beljaars_surface_stress_tbd
+  units = N m-2
+  type = real | kind = kind_phys
+  dimensions = (horizontal_loop_extent)
+  intent = out
+[ tauy ]
+  standard_name = northward_beljaars_surface_stress_tbd
+  units = N m-2
+  type = real | kind = kind_phys
+  dimensions = (horizontal_loop_extent)
+  intent = out
+[ errmsg ]
+  standard_name = ccpp_error_message
+  units = none
+  type = character | kind = len=*
+  dimensions = ()
+  intent = out
+[ errflg ]
+  standard_name = ccpp_error_code
+  units = 1
+  type = integer
+  dimensions = ()
+  intent = out

schemes/beljaars_drag/beljaars_drag_interstitials.F90

@@ -0,0 +1,118 @@
+! Interstitial schemes for Beljaars drag.
+!
+! Put in a separate file than beljaars_drag.F90 because
+! scheme files with multiple schemes cannot have a namelist file attached.
+module beljaars_drag_interstitials
+
+  implicit none
+  private
+
+  ! public CCPP-compliant subroutines
+  public :: beljaars_add_wind_damping_rate_run
+  public :: beljaars_add_updated_residual_stress_run
+contains
+
+  ! Add Beljaars drag to the wind damping rate for vertical diffusion.
+  ! Has to run after vertical_diffusion_wind_damping_rate
+!> \section arg_table_beljaars_add_wind_damping_rate_run Argument Table
+!! \htmlinclude arg_table_beljaars_add_wind_damping_rate_run.html
+  pure subroutine beljaars_add_wind_damping_rate_run( &
+    ncol, pver, &
+    dragblj, &
+    tau_damp_rate, &
+    errmsg, errflg)
+
+    use ccpp_kinds, only: kind_phys
+
+    ! Input arguments
+    integer,            intent(in)    :: ncol
+    integer,            intent(in)    :: pver
+    real(kind_phys),    intent(in)    :: dragblj(:, :)    ! Drag profile from Beljaars SGO form drag > 0. [s-1]
+
+    ! Input/output arguments
+    real(kind_phys),    intent(inout) :: tau_damp_rate(:,:) ! Rate at which external (surface) stress damps wind speeds [s-1]
+
+    ! Output arguments
+    character(len=*),   intent(out)   :: errmsg  ! error message
+    integer,            intent(out)   :: errflg  ! error flag
+
+    errmsg = ''
+    errflg = 0
+
+    ! Beljaars et al SGO scheme incorporated here. It
+    ! appears as a "3D" tau_damp_rate specification.
+    tau_damp_rate(:ncol,:pver) = tau_damp_rate(:ncol,:pver) + dragblj(:ncol,:pver)
+
+  end subroutine beljaars_add_wind_damping_rate_run
+
+  ! Add Beljaars stress using updated provisional winds to surface stresses used
+  ! for horizontal momentum diffusion - kinetic energy dissipation.
+!> \section arg_table_beljaars_add_updated_residual_stress_run Argument Table
+!! \htmlinclude arg_table_beljaars_add_updated_residual_stress_run.html
+  subroutine beljaars_add_updated_residual_stress_run( &
+    ncol, pver, &
+    gravit, &
+    p, &
+    do_iss, &
+    itaures, &
+    dragblj, &
+    u1, v1, & ! provisionally updated winds
+    ! input/output
+    tauresx, tauresy, &
+    errmsg, errflg)
+
+    use ccpp_kinds, only: kind_phys
+    use coords_1d, only: Coords1D
+
+    ! Input arguments
+    integer,         intent(in)       :: ncol
+    integer,         intent(in)       :: pver
+    real(kind_phys), intent(in)       :: gravit
+    type(coords1d),  intent(in)       :: p                ! Pressure coordinates [Pa]
+    logical,         intent(in)       :: do_iss           ! Flag for implicit surface stress (namelist from vdiff)
+    logical,         intent(in)       :: itaures          ! Flag for updating tauresx tauresy in this subroutine.
+    real(kind_phys), intent(in)       :: dragblj(:, :)    ! Drag profile from Beljaars SGO form drag > 0. [s-1]
+    real(kind_phys), intent(in)       :: u1(:,:)          ! After vertical diffusion u-wind [m s-1]
+    real(kind_phys), intent(in)       :: v1(:,:)          ! After vertical diffusion v-wind [m s-1]
+
+    ! Input/output arguments
+    real(kind_phys), intent(inout)    :: tauresx(:)       ! Partially updated residual surface stress using provisional winds
+    real(kind_phys), intent(inout)    :: tauresy(:)
+    character(len=*),   intent(out)   :: errmsg           ! error message
+    integer,            intent(out)   :: errflg           ! error flag
+
+    integer :: i, k
+    real(kind_phys) :: taubljx(ncol)                      ! recomputed explicit/residual beljaars stress
+    real(kind_phys) :: taubljy(ncol)                      ! recomputed explicit/residual beljaars stress
+
+    errmsg = ''
+    errflg = 0
+
+    if(.not. do_iss) then
+      ! Not added to residual stress if implicit surface stress is off
+      return
+    end if
+
+    if(.not. itaures) then
+      ! Not added to residual stress if residual stress is not requested to be updated.
+      return
+    end if
+
+    do i = 1, ncol
+      ! Add vertically-integrated Beljaars drag to residual stress
+      ! these are calculated using provisionally-updated winds.
+      ! These are used only locally and not written back to model state.
+      taubljx(i) = 0._kind_phys
+      taubljy(i) = 0._kind_phys
+      do k = 1, pver
+        taubljx(i) = taubljx(i) + (1._kind_phys/gravit)*dragblj(i, k)*u1(i, k)*p%del(i, k)
+        taubljy(i) = taubljy(i) + (1._kind_phys/gravit)*dragblj(i, k)*v1(i, k)*p%del(i, k)
+      end do
+    end do
+
+    tauresx(:ncol) = tauresx(:ncol) + taubljx(:ncol)
+    tauresy(:ncol) = tauresy(:ncol) + taubljy(:ncol)
+
+  end subroutine beljaars_add_updated_residual_stress_run
+
+end module beljaars_drag_interstitials