user_change_diffusivity.F90

!> Increments the diapycnal diffusivity in a specified band of latitudes and densities.
module user_change_diffusivity

! This file is part of MOM6. See LICENSE.md for the license.

use mom_diag_mediator, only : diag_ctrl, time_type
use mom_error_handler, only : mom_error, is_root_pe, fatal, warning, note
use mom_file_parser,   only : get_param, log_version, param_file_type
use mom_grid,          only : ocean_grid_type
use mom_unit_scaling, only : unit_scale_type
use mom_variables,     only : thermo_var_ptrs, vertvisc_type, p3d
use mom_verticalgrid,  only : verticalgrid_type
use mom_eos,           only : calculate_density, eos_domain

implicit none ; private

#include <MOM_memory.h>

public user_change_diff, user_change_diff_init, user_change_diff_end

! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]".  If the units
! vary with the Boussinesq approximation, the Boussinesq variant is given first.

!> Control structure for user_change_diffusivity
type, public :: user_change_diff_cs ; private
  real :: kd_add        !< The scale of a diffusivity that is added everywhere
                        !! without any filtering or scaling [Z2 T-1 ~> m2 s-1].
  real :: lat_range(4)  !< 4 values that define the latitude range over which
                        !! a diffusivity scaled by Kd_add is added [degLat].
  real :: rho_range(4)  !< 4 values that define the coordinate potential
                        !! density range over which a diffusivity scaled by
                        !! Kd_add is added [R ~> kg m-3].
  logical :: use_abs_lat  !< If true, use the absolute value of latitude when
                          !! setting lat_range.
  type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
                          !! regulate the timing of diagnostic output.
end type user_change_diff_cs

contains

!> This subroutine provides an interface for a user to use to modify the
!! main code to alter the diffusivities as needed.  The specific example
!! implemented here augments the diffusivity for a specified range of latitude
!! and coordinate potential density.
subroutine user_change_diff(h, tv, G, GV, US, CS, Kd_lay, Kd_int, T_f, S_f, Kd_int_add)
  type(ocean_grid_type),                    intent(in)    :: G   !< The ocean's grid structure.
  type(verticalgrid_type),                  intent(in)    :: GV  !< The ocean's vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in)    :: h   !< Layer thickness [H ~> m or kg m-2].
  type(thermo_var_ptrs),                    intent(in)    :: tv  !< A structure containing pointers
                                                                 !! to any available thermodynamic
                                                                 !! fields. Absent fields have NULL ptrs.
  type(unit_scale_type),                    intent(in)    :: US  !< A dimensional unit scaling type
  type(user_change_diff_cs),                pointer       :: CS  !< This module's control structure.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   optional, intent(inout) :: Kd_lay !< The diapycnal diffusivity of
                                                                  !! each layer [Z2 T-1 ~> m2 s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), optional, intent(inout) :: Kd_int !< The diapycnal diffusivity
                                                                  !! at each interface [Z2 T-1 ~> m2 s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   optional, intent(in)    :: T_f !< Temperature with massless
                                                                  !! layers filled in vertically [degC].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   optional, intent(in)    :: S_f !< Salinity with massless
                                                                  !! layers filled in vertically [ppt].
  real, dimension(:,:,:),                     optional, pointer       :: Kd_int_add !< The diapycnal
                                                                  !! diffusivity that is being added at
                                                                  !! each interface [Z2 T-1 ~> m2 s-1].
  ! Local variables
  real :: Rcv(szi_(g),szk_(g)) ! The coordinate density in layers [R ~> kg m-3].
  real :: p_ref(szi_(g))       ! An array of tv%P_Ref pressures [R L2 T-2 ~> Pa].
  real :: rho_fn      ! The density dependence of the input function, 0-1 [nondim].
  real :: lat_fn      ! The latitude dependence of the input function, 0-1 [nondim].
  logical :: use_EOS  ! If true, density is calculated from T & S using an
                      ! equation of state.
  logical :: store_Kd_add  ! Save the added diffusivity as a diagnostic if true.
  integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
  integer :: i, j, k, is, ie, js, je, nz
  integer :: isd, ied, jsd, jed

  real :: kappa_fill  ! diffusivity used to fill massless layers
  real :: dt_fill     ! timestep used to fill massless layers
  character(len=200) :: mesg

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed

  if (.not.associated(cs)) call mom_error(fatal,"user_set_diffusivity: "//&
         "Module must be initialized before it is used.")

  use_eos = associated(tv%eqn_of_state)
  if (.not.use_eos) return
  store_kd_add = .false.
  if (present(kd_int_add)) store_kd_add = associated(kd_int_add)

  if (.not.range_ok(cs%lat_range)) then
    write(mesg, '(4(1pe15.6))') cs%lat_range(1:4)
    call mom_error(fatal, "user_set_diffusivity: bad latitude range: \n  "//&
                    trim(mesg))
  endif
  if (.not.range_ok(cs%rho_range)) then
    write(mesg, '(4(1pe15.6))') cs%rho_range(1:4)
    call mom_error(fatal, "user_set_diffusivity: bad density range: \n  "//&
                    trim(mesg))
  endif

  if (store_kd_add) kd_int_add(:,:,:) = 0.0

  do i=is,ie ; p_ref(i) = tv%P_Ref ; enddo
  eosdom(:) = eos_domain(g%HI)
  do j=js,je
    if (present(t_f) .and. present(s_f)) then
      do k=1,nz
        call calculate_density(t_f(:,j,k), s_f(:,j,k), p_ref, rcv(:,k), tv%eqn_of_state, eosdom)
      enddo
    else
      do k=1,nz
        call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_ref, rcv(:,k), tv%eqn_of_state, eosdom)
      enddo
    endif

    if (present(kd_lay)) then
      do k=1,nz ; do i=is,ie
        if (cs%use_abs_lat) then
          lat_fn = val_weights(abs(g%geoLatT(i,j)), cs%lat_range)
        else
          lat_fn = val_weights(g%geoLatT(i,j), cs%lat_range)
        endif
        rho_fn = val_weights(rcv(i,k), cs%rho_range)
        if (rho_fn * lat_fn > 0.0) &
          kd_lay(i,j,k) = kd_lay(i,j,k) + cs%Kd_add * rho_fn * lat_fn
      enddo ; enddo
    endif
    if (present(kd_int)) then
      do k=2,nz ; do i=is,ie
        if (cs%use_abs_lat) then
          lat_fn = val_weights(abs(g%geoLatT(i,j)), cs%lat_range)
        else
          lat_fn = val_weights(g%geoLatT(i,j), cs%lat_range)
        endif
        rho_fn = val_weights( 0.5*(rcv(i,k-1) + rcv(i,k)), cs%rho_range)
        if (rho_fn * lat_fn > 0.0) then
          kd_int(i,j,k) = kd_int(i,j,k) + cs%Kd_add * rho_fn * lat_fn
          if (store_kd_add) kd_int_add(i,j,k) = cs%Kd_add * rho_fn * lat_fn
        endif
      enddo ; enddo
    endif
  enddo

end subroutine user_change_diff

!> This subroutine checks whether the 4 values of range are in ascending order.
function range_ok(range) result(OK)
  real, dimension(4), intent(in) :: range  !< Four values to check.
  logical                        :: OK     !< Return value.

  ok = ((range(1) <= range(2)) .and. (range(2) <= range(3)) .and. &
        (range(3) <= range(4)))

end function range_ok

!> This subroutine returns a value that goes smoothly from 0 to 1, stays
!! at 1, and then goes smoothly back to 0 at the four values of range.  The
!! transitions are cubic, and have zero first derivatives where the curves
!! hit 0 and 1.  The values in range must be in ascending order, as can be
!! checked by calling range_OK.
function val_weights(val, range) result(ans)
  real,               intent(in) :: val    !< Value for which we need an answer [arbitrary units].
  real, dimension(4), intent(in) :: range  !< Range over which the answer is non-zero [arbitrary units].
  real                           :: ans    !< Return value [nondim].
  ! Local variables
  real :: x   ! A nondimensional number between 0 and 1.

  ans = 0.0
  if ((val > range(1)) .and. (val < range(4))) then
    if (val < range(2)) then
      ! x goes from 0 to 1; ans goes from 0 to 1, with 0 derivatives at the ends.
      x = (val - range(1)) / (range(2) - range(1))
      ans = x**2 * (3.0 - 2.0 * x)
    elseif (val > range(3)) then
      ! x goes from 0 to 1; ans goes from 0 to 1, with 0 derivatives at the ends.
      x = (range(4) - val) / (range(4) - range(3))
      ans = x**2 * (3.0 - 2.0 * x)
    else
      ans = 1.0
    endif
  endif

end function val_weights

!> Set up the module control structure.
subroutine user_change_diff_init(Time, G, GV, US, param_file, diag, CS)
  type(time_type),           intent(in)    :: Time       !< The current model time.
  type(ocean_grid_type),     intent(in)    :: G          !< The ocean's grid structure.
  type(verticalgrid_type),   intent(in)    :: GV         !< The ocean's vertical grid structure
  type(unit_scale_type),     intent(in)    :: US         !< A dimensional unit scaling type
  type(param_file_type),     intent(in)    :: param_file !< A structure indicating the
                                                         !! open file to parse for
                                                         !! model parameter values.
  type(diag_ctrl), target,   intent(inout) :: diag       !< A structure that is used to
                                                         !! regulate diagnostic output.
  type(user_change_diff_cs), pointer       :: CS         !< A pointer that is set to
                                                         !! point to the control
                                                         !! structure for this module.

! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "user_set_diffusivity"  ! This module's name.
  character(len=200) :: mesg
  integer :: i, j, is, ie, js, je

  if (associated(cs)) then
    call mom_error(warning, "diabatic_entrain_init called with an associated "// &
                            "control structure.")
    return
  endif
  allocate(cs)

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec

  cs%diag => diag

  ! Read all relevant parameters and write them to the model log.
  call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl, "USER_KD_ADD", cs%Kd_add, &
                 "A user-specified additional diffusivity over a range of "//&
                 "latitude and density.", default=0.0, units="m2 s-1", &
                 scale=us%m2_s_to_Z2_T)
  if (cs%Kd_add /= 0.0) then
    call get_param(param_file, mdl, "USER_KD_ADD_LAT_RANGE", cs%lat_range(:), &
                 "Four successive values that define a range of latitudes "//&
                 "over which the user-specified extra diffusivity is "//&
                 "applied.  The four values specify the latitudes at "//&
                 "which the extra diffusivity starts to increase from 0, "//&
                 "hits its full value, starts to decrease again, and is "//&
                 "back to 0.", units="degree", default=-1.0e9)
    call get_param(param_file, mdl, "USER_KD_ADD_RHO_RANGE", cs%rho_range(:), &
                 "Four successive values that define a range of potential "//&
                 "densities over which the user-given extra diffusivity "//&
                 "is applied.  The four values specify the density at "//&
                 "which the extra diffusivity starts to increase from 0, "//&
                 "hits its full value, starts to decrease again, and is "//&
                 "back to 0.", units="kg m-3", default=-1.0e9, scale=us%kg_m3_to_R)
    call get_param(param_file, mdl, "USER_KD_ADD_USE_ABS_LAT", cs%use_abs_lat, &
                 "If true, use the absolute value of latitude when "//&
                 "checking whether a point fits into range of latitudes.", &
                 default=.false.)
  endif

  if (.not.range_ok(cs%lat_range)) then
    write(mesg, '(4(1pe15.6))') cs%lat_range(1:4)
    call mom_error(fatal, "user_set_diffusivity: bad latitude range: \n  "//&
                    trim(mesg))
  endif
  if (.not.range_ok(cs%rho_range)) then
    write(mesg, '(4(1pe15.6))') cs%rho_range(1:4)
    call mom_error(fatal, "user_set_diffusivity: bad density range: \n  "//&
                    trim(mesg))
  endif

end subroutine user_change_diff_init

!> Clean up the module control structure.
subroutine user_change_diff_end(CS)
  type(user_change_diff_cs), pointer :: CS !< A pointer that is set to point to the control
                                           !! structure for this module.

  if (associated(cs)) deallocate(cs)

end subroutine user_change_diff_end

end module user_change_diffusivity