user_surface_forcing.F90

! This file is part of MOM6, the Modular Ocean Model version 6.

! See the LICENSE file for licensing information.

! SPDX-License-Identifier: Apache-2.0

!> Template for user to code up surface forcing.

module user_surface_forcing

use mom_diag_mediator, only : post_data, query_averaging_enabled

use mom_diag_mediator, only : register_diag_field, diag_ctrl, safe_alloc_ptr

use mom_domains, only : pass_var, pass_vector, agrid

use mom_error_handler, only : mom_error, fatal, warning, is_root_pe

use mom_file_parser, only : get_param, param_file_type, log_version

use mom_forcing_type, only : forcing, mech_forcing

use mom_forcing_type, only : allocate_forcing_type, allocate_mech_forcing

use mom_grid, only : ocean_grid_type

use mom_time_manager, only : time_type, operator(+), operator(/)

use mom_tracer_flow_control, only : call_tracer_set_forcing

use mom_tracer_flow_control, only : tracer_flow_control_cs

use mom_unit_scaling, only : unit_scale_type

use mom_variables, only : surface

implicit none ; private

public user_wind_forcing, user_buoyancy_forcing, user_surface_forcing_init

!>   This control structure should be used to store any run-time variables

!! associated with the user-specified forcing.

!!

!! It can be readily modified for a specific case, and because it is private there

!! will be no changes needed in other code (although they will have to be recompiled).

type, public :: user_surface_forcing_cs ; private

  !   The variables in the canonical example are used for some common

  ! cases, but do not need to be used.

  logical :: use_temperature !< If true, temperature and salinity are used as state variables.

  logical :: restorebuoy     !< If true, use restoring surface buoyancy forcing.

  real :: rho0               !< The density used in the Boussinesq approximation [R ~> kg m-3].

  real :: g_earth            !< The gravitational acceleration [L2 Z-1 T-2 ~> m s-2].

  real :: flux_const         !< The restoring rate at the surface [Z T-1 ~> m s-1].

  real :: rho_restore        !< The density that is used to convert piston velocities into salt

                             !! or heat fluxes with salinity or temperature restoring [R ~> kg m-3]

  real :: gust_const         !< A constant unresolved background gustiness

                             !! that contributes to ustar [R Z2 T-2 ~> Pa].

  type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the

                             !! timing of diagnostic output.

end type user_surface_forcing_cs

contains

!> This subroutine sets the surface wind stresses, forces%taux and forces%tauy, in [R Z L T-2 ~> Pa].

!! These are the stresses in the direction of the model grid (i.e. the same

!! direction as the u- and v- velocities).

subroutine user_wind_forcing(sfc_state, forces, day, G, US, CS)

  type(surface),                 intent(inout) :: sfc_state !< A structure containing fields that

                                                       !! describe the surface state of the ocean.

  type(mech_forcing),            intent(inout) :: forces !< A structure with the driving mechanical forces

  type(time_type),               intent(in)    :: day  !< The time of the fluxes

  type(ocean_grid_type),         intent(inout) :: g    !< The ocean's grid structure

  type(unit_scale_type),         intent(in)    :: us   !< A dimensional unit scaling type

  type(user_surface_forcing_cs), pointer       :: cs   !< A pointer to the control structure returned

                                                       !! by a previous call to user_surface_forcing_init

  ! Local variables

  integer :: i, j, is, ie, js, je, isq, ieq, jsq, jeq

  !   When modifying the code, comment out this error message.  It is here

  ! so that the original (unmodified) version is not accidentally used.

  call mom_error(fatal, "User_wind_surface_forcing: " // &

     "User forcing routine called without modification." )

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

  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB

  ! Allocate the forcing arrays, if necessary.

  call allocate_mech_forcing(g, forces, stress=.true., ustar=.true., tau_mag=.true.)

  !  Set the surface wind stresses, in units of [R L Z T-2 ~> Pa].  A positive taux

  !  accelerates the ocean to the (pseudo-)east.

  !  The i-loop extends to is-1 so that taux can be used later in the

  ! calculation of ustar - otherwise the lower bound would be Isq.

  do j=js,je ; do i=is-1,ieq

    ! Change this to the desired expression.

    forces%taux(i,j) = g%mask2dCu(i,j) * 0.0*us%Pa_to_RLZ_T2

  enddo ; enddo

  do j=js-1,jeq ; do i=is,ie

    forces%tauy(i,j) = g%mask2dCv(i,j) * 0.0  ! Change this to the desired expression.

  enddo ; enddo

  !    Set the surface friction velocity, in units of [Z T-1 ~> m s-1].  ustar

  !  is always positive.

  if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie

    !  This expression can be changed if desired, but need not be.

    forces%tau_mag(i,j) = g%mask2dT(i,j) * (cs%gust_const + &

            us%L_to_Z*sqrt(0.5*((forces%taux(i-1,j)**2) + (forces%taux(i,j)**2)) + &

                           0.5*((forces%tauy(i,j-1)**2) + (forces%tauy(i,j)**2))))

    if (associated(forces%ustar)) &

      forces%ustar(i,j) = g%mask2dT(i,j) * sqrt(forces%tau_mag(i,j) * (1.0/cs%Rho0))

  enddo ; enddo ; endif

end subroutine user_wind_forcing

!>    This subroutine specifies the current surface fluxes of buoyancy or

!!  temperature and fresh water.  It may also be modified to add

!!  surface fluxes of user provided tracers.

subroutine user_buoyancy_forcing(sfc_state, fluxes, day, dt, G, US, CS)

  type(surface),                 intent(inout) :: sfc_state !< A structure containing fields that

                                                       !! describe the surface state of the ocean.

  type(forcing),                 intent(inout) :: fluxes !< A structure containing thermodynamic forcing fields

  type(time_type),               intent(in)    :: day  !< The time of the fluxes

  real,                          intent(in)    :: dt   !< The amount of time over which

                                                       !! the fluxes apply [T ~> s]

  type(ocean_grid_type),         intent(in)    :: g    !< The ocean's grid structure

  type(unit_scale_type),         intent(in)    :: us   !< A dimensional unit scaling type

  type(user_surface_forcing_cs), pointer       :: cs   !< A pointer to the control structure returned

                                                       !! by a previous call to user_surface_forcing_init

!    This subroutine specifies the current surface fluxes of buoyancy or

!  temperature and fresh water.  It may also be modified to add

!  surface fluxes of user provided tracers.

!    When temperature is used, there are long list of fluxes that need to be

!  set - essentially the same as for a full coupled model, but most of these

!  can be simply set to zero.  The net fresh water flux should probably be

!  set in fluxes%evap and fluxes%lprec, with any salinity restoring

!  appearing in fluxes%vprec, and the other water flux components

!  (fprec, lrunoff and frunoff) left as arrays full of zeros.

!  Evap is usually negative and precip is usually positive.  All heat fluxes

!  are in W m-2 and positive for heat going into the ocean.  All fresh water

!  fluxes are in [R Z T-1 ~> kg m-2 s-1] and positive for water moving into the ocean.

  ! Local variables

  real :: temp_restore   ! The temperature that is being restored toward [C ~> degC].

  real :: salin_restore  ! The salinity that is being restored toward [S ~> ppt]

  real :: density_restore  ! The potential density that is being restored

                         ! toward [R ~> kg m-3].

  real :: rhoxcp         ! The mean density times the heat capacity [Q R C-1 ~> J m-3 degC-1].

  real :: buoy_rest_const  ! A constant relating density anomalies to the

                           ! restoring buoyancy flux [L2 T-3 R-1 ~> m5 s-3 kg-1].

  integer :: i, j, is, ie, js, je

  integer :: isd, ied, jsd, jed

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

  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed

  !   When modifying the code, comment out this error message.  It is here

  ! so that the original (unmodified) version is not accidentally used.

  call mom_error(fatal, "User_buoyancy_surface_forcing: " // &

    "User forcing routine called without modification." )

  ! Allocate and zero out the forcing arrays, as necessary.  This portion is

  ! usually not changed.

  if (cs%use_temperature) then

    call safe_alloc_ptr(fluxes%evap, isd, ied, jsd, jed)

    call safe_alloc_ptr(fluxes%lprec, isd, ied, jsd, jed)

    call safe_alloc_ptr(fluxes%fprec, isd, ied, jsd, jed)

    call safe_alloc_ptr(fluxes%lrunoff, isd, ied, jsd, jed)

    call safe_alloc_ptr(fluxes%frunoff, isd, ied, jsd, jed)

    call safe_alloc_ptr(fluxes%vprec, isd, ied, jsd, jed)

    call safe_alloc_ptr(fluxes%sw, isd, ied, jsd, jed)

    call safe_alloc_ptr(fluxes%lw, isd, ied, jsd, jed)

    call safe_alloc_ptr(fluxes%latent, isd, ied, jsd, jed)

    call safe_alloc_ptr(fluxes%sens, isd, ied, jsd, jed)

  else ! This is the buoyancy only mode.

    call safe_alloc_ptr(fluxes%buoy, isd, ied, jsd, jed)

  endif

  ! MODIFY THE CODE IN THE FOLLOWING LOOPS TO SET THE BUOYANCY FORCING TERMS.

  if ( cs%use_temperature ) then

    ! Set whichever fluxes are to be used here.  Any fluxes that

    ! are always zero do not need to be changed here.

    do j=js,je ; do i=is,ie

      ! Fluxes of fresh water through the surface are in units of [R Z T-1 ~> kg m-2 s-1]

      ! and are positive downward - i.e. evaporation should be negative.

      fluxes%evap(i,j) = -0.0 * g%mask2dT(i,j)

      fluxes%lprec(i,j) = 0.0 * g%mask2dT(i,j)

      ! vprec will be set later, if it is needed for salinity restoring.

      fluxes%vprec(i,j) = 0.0

      !   Heat fluxes are in units of [Q R Z T-1 ~> W m-2] and are positive into the ocean.

      fluxes%lw(i,j) = 0.0 * g%mask2dT(i,j)

      fluxes%latent(i,j) = 0.0 * g%mask2dT(i,j)

      fluxes%sens(i,j) = 0.0 * g%mask2dT(i,j)

      fluxes%sw(i,j) = 0.0 * g%mask2dT(i,j)

    enddo ; enddo

  else ! This is the buoyancy only mode.

    do j=js,je ; do i=is,ie

      !   fluxes%buoy is the buoyancy flux into the ocean [L2 T-3 ~> m2 s-3].  A positive

      ! buoyancy flux is of the same sign as heating the ocean.

      fluxes%buoy(i,j) = 0.0 * g%mask2dT(i,j)

    enddo ; enddo

  endif

  if (cs%restorebuoy) then

    if (cs%use_temperature) then

      call safe_alloc_ptr(fluxes%heat_added, isd, ied, jsd, jed)

      !   When modifying the code, comment out this error message.  It is here

      ! so that the original (unmodified) version is not accidentally used.

      call mom_error(fatal, "User_buoyancy_surface_forcing: " // &

        "Temperature and salinity restoring used without modification." )

      rhoxcp = cs%rho_restore * fluxes%C_p

      do j=js,je ; do i=is,ie

        !   Set Temp_restore and Salin_restore to the temperature (in [C ~> degC]) and

        ! salinity (in [S ~> ppt]) that are being restored toward.

        temp_restore = 0.0

        salin_restore = 0.0

        fluxes%heat_added(i,j) = (g%mask2dT(i,j) * (rhoxcp * cs%Flux_const)) * &

            (temp_restore - sfc_state%SST(i,j))

        fluxes%vprec(i,j) = - (g%mask2dT(i,j) * (cs%rho_restore*cs%Flux_const)) * &

            ((salin_restore - sfc_state%SSS(i,j)) / (0.5 * (salin_restore + sfc_state%SSS(i,j))))

      enddo ; enddo

    else

      !   When modifying the code, comment out this error message.  It is here

      ! so that the original (unmodified) version is not accidentally used.

      call mom_error(fatal, "User_buoyancy_surface_forcing: " // &

        "Buoyancy restoring used without modification." )

      ! The -1 is because density has the opposite sign to buoyancy.

      buoy_rest_const = -1.0 * (cs%G_Earth * cs%Flux_const) / cs%rho_restore

      do j=js,je ; do i=is,ie

       !   Set density_restore to an expression for the surface potential

       ! density [R ~> kg m-3] that is being restored toward.

        density_restore = 1030.0*us%kg_m3_to_R

        fluxes%buoy(i,j) = g%mask2dT(i,j) * buoy_rest_const * &

                          (density_restore - sfc_state%sfc_density(i,j))

      enddo ; enddo

    endif

  endif                                             ! end RESTOREBUOY

end subroutine user_buoyancy_forcing

!> This subroutine initializes the USER_surface_forcing module

subroutine user_surface_forcing_init(Time, G, 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(unit_scale_type),         intent(in) :: us   !< A dimensional unit scaling type

  type(param_file_type),         intent(in) :: param_file !< A structure to parse for run-time parameters

  type(diag_ctrl), target,       intent(in) :: diag !< A structure that is used to regulate diagnostic output.

  type(user_surface_forcing_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_surface_forcing" ! This module's name.

  if (associated(cs)) then

    call mom_error(warning, "USER_surface_forcing_init called with an associated "// &

                             "control structure.")

    return

  endif

  allocate(cs)

  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, "ENABLE_THERMODYNAMICS", cs%use_temperature, &

                 "If true, Temperature and salinity are used as state "//&

                 "variables.", default=.true.)

  call get_param(param_file, mdl, "G_EARTH", cs%G_Earth, &

                 "The gravitational acceleration of the Earth.", &

                 units="m s-2", default = 9.80, scale=us%m_to_L**2*us%Z_to_m*us%T_to_s**2)

  call get_param(param_file, mdl, "RHO_0", cs%Rho0, &

                 "The mean ocean density used with BOUSSINESQ true to "//&

                 "calculate accelerations and the mass for conservation "//&

                 "properties, or with BOUSSINESQ false to convert some "//&

                 "parameters from vertical units of m to kg m-2.", &

                 units="kg m-3", default=1035.0, scale=us%kg_m3_to_R)

  call get_param(param_file, mdl, "GUST_CONST", cs%gust_const, &

                 "The background gustiness in the winds.", &

                 units="Pa", default=0.0, scale=us%Pa_to_RLZ_T2*us%L_to_Z)

  call get_param(param_file, mdl, "RESTOREBUOY", cs%restorebuoy, &

                 "If true, the buoyancy fluxes drive the model back "//&

                 "toward some specified surface state with a rate "//&

                 "given by FLUXCONST.", default= .false.)

  if (cs%restorebuoy) then

    call get_param(param_file, mdl, "FLUXCONST", cs%Flux_const, &

                 "The constant that relates the restoring surface fluxes to the relative "//&

                 "surface anomalies (akin to a piston velocity).  Note the non-MKS units.", &

                 default=0.0, units="m day-1", scale=us%m_to_Z/(86400.0*us%s_to_T))

  endif

  call get_param(param_file, mdl, "RESTORE_FLUX_RHO", cs%rho_restore, &

                 "The density that is used to convert piston velocities into salt or heat "//&

                 "fluxes with RESTORE_SALINITY or RESTORE_TEMPERATURE.", &

                 units="kg m-3", default=cs%Rho0*us%R_to_kg_m3, scale=us%kg_m3_to_R, &

                 do_not_log=(cs%Flux_const==0.0).or.(.not.cs%restorebuoy))

end subroutine user_surface_forcing_init

!! \namespace user_surface_forcing

!!

!! This file contains the subroutines that a user should modify to

!! to set the surface wind stresses and fluxes of buoyancy or

!! temperature and fresh water.  They are called when the run-time

!! parameters WIND_CONFIG or BUOY_CONFIG are set to "USER".  The

!! standard version has simple examples, along with run-time error

!! messages that will cause the model to abort if this code has no

!! been modified.  This code is intended for use with relatively

!! simple specifications of the forcing.  For more complicated forms,

!! it is probably a good idea to read the forcing from input files

!! using "file" for WIND_CONFIG and BUOY_CONFIG.

!!

!! USER_wind_forcing() should set the surface wind stresses (taux and

!! tauy) perhaps along with the surface friction velocity (ustar).

!!

!! USER_buoyancy() forcing is used to set the surface buoyancy

!! forcing, which may include a number of fresh water flux fields

!! (evap, lprec, fprec, lrunoff, frunoff, and

!! vprec) and the surface heat fluxes (sw, lw, latent and sens)

!! if temperature and salinity are state variables, or it may simply

!! be the buoyancy flux if it is not.  This routine also has coded a

!! restoring to surface values of temperature and salinity.

end module user_surface_forcing