MOM_surface_forcing.F90

!> Functions that calculate the surface wind stresses and fluxes of buoyancy
!! or temperature/salinity andfresh water, in ocean-only (solo) mode.
!!
!! These functions are called every time step, even if the wind stresses
!! or buoyancy fluxes are constant in time - in that case these routines
!! return quickly without doing anything.  In addition, any I/O of forcing
!! fields is controlled by surface_forcing_init, located in this file.
module mom_surface_forcing

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

use mom_constants,           only : hlv, hlf
use mom_cpu_clock,           only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
use mom_cpu_clock,           only : clock_module
use mom_diag_mediator,       only : post_data, query_averaging_enabled
use mom_diag_mediator,       only : diag_ctrl, safe_alloc_ptr
use mom_domains,             only : pass_var, pass_vector, agrid, to_south, to_west, to_all
use mom_domains,             only : fill_symmetric_edges, cgrid_ne
use mom_error_handler,       only : mom_error, fatal, warning, mom_mesg, is_root_pe
use mom_error_handler,       only : calltree_enter, calltree_leave
use mom_file_parser,         only : get_param, log_version, param_file_type
use mom_string_functions,    only : uppercase
use mom_forcing_type,        only : forcing, mech_forcing
use mom_forcing_type,        only : set_net_mass_forcing, copy_common_forcing_fields
use mom_forcing_type,        only : set_derived_forcing_fields
use mom_forcing_type,        only : forcing_diags, mech_forcing_diags, register_forcing_type_diags
use mom_forcing_type,        only : allocate_forcing_type, deallocate_forcing_type
use mom_forcing_type,        only : allocate_mech_forcing, deallocate_mech_forcing
use mom_grid,                only : ocean_grid_type
use mom_get_input,           only : get_mom_input, directories
use mom_io,                  only : file_exists, mom_read_data, mom_read_vector, slasher
use mom_io,                  only : east_face, north_face, num_timelevels
use mom_restart,             only : register_restart_field, restart_init, mom_restart_cs
use mom_restart,             only : restart_init_end, save_restart, restore_state
use mom_time_manager,        only : time_type, operator(+), operator(/), get_time, time_type_to_real
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
use meso_surface_forcing,    only : meso_buoyancy_forcing
use meso_surface_forcing,    only : meso_surface_forcing_init, meso_surface_forcing_cs
use user_surface_forcing,    only : user_wind_forcing, user_buoyancy_forcing
use user_surface_forcing,    only : user_surface_forcing_init, user_surface_forcing_cs
use user_revise_forcing,     only : user_alter_forcing, user_revise_forcing_init
use user_revise_forcing,     only : user_revise_forcing_cs
use idealized_hurricane, only : idealized_hurricane_wind_init
use idealized_hurricane, only : idealized_hurricane_wind_forcing, scm_idealized_hurricane_wind_forcing
use idealized_hurricane, only : idealized_hurricane_cs
use scm_cvmix_tests,         only : scm_cvmix_tests_surface_forcing_init
use scm_cvmix_tests,         only : scm_cvmix_tests_wind_forcing
use scm_cvmix_tests,         only : scm_cvmix_tests_buoyancy_forcing
use scm_cvmix_tests,         only : scm_cvmix_tests_cs
use bfb_surface_forcing,    only : bfb_buoyancy_forcing
use bfb_surface_forcing,    only : bfb_surface_forcing_init, bfb_surface_forcing_cs
use dumbbell_surface_forcing,    only : dumbbell_surface_forcing_init, dumbbell_surface_forcing_cs
use dumbbell_surface_forcing, only    : dumbbell_buoyancy_forcing
use data_override_mod, only : data_override_init, data_override

implicit none ; private

#include <MOM_memory.h>

public set_forcing
public surface_forcing_init
public forcing_save_restart

!> Structure containing pointers to the forcing fields that may be used to drive MOM.
!!  All fluxes are positive into the ocean.
type, public :: surface_forcing_cs ; private

  logical :: use_temperature    !< if true, temp & salinity used as state variables
  logical :: restorebuoy        !< if true, use restoring surface buoyancy forcing
  logical :: adiabatic          !< if true, no diapycnal mass fluxes or surface buoyancy forcing
  logical :: variable_winds     !< if true, wind stresses vary with time
  logical :: variable_buoyforce !< if true, buoyancy forcing varies with time.
  real    :: south_lat          !< southern latitude of the domain
  real    :: len_lat            !< domain length in latitude

  real :: rho0                  !< Boussinesq reference density [R ~> kg m-3]
  real :: g_earth               !< gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
  real :: flux_const            !< piston velocity for surface restoring [Z T-1 ~> m s-1]
  real :: flux_const_t          !< piston velocity for surface temperature restoring [m s-1]
  real :: flux_const_s          !< piston velocity for surface salinity restoring [Z T-1 ~> m s-1]
  real :: latent_heat_fusion    !< latent heat of fusion times [Q ~> J kg-1]
  real :: latent_heat_vapor     !< latent heat of vaporization [Q ~> J kg-1]
  real :: tau_x0                !< Constant zonal wind stress used in the WIND_CONFIG="const" forcing
  real :: tau_y0                !< Constant meridional wind stress used in the WIND_CONFIG="const" forcing

  real    :: gust_const                 !< constant unresolved background gustiness for ustar [R L Z T-1 ~> Pa]
  logical :: read_gust_2d               !< if true, use 2-dimensional gustiness supplied from a file
  real, pointer :: gust(:,:) => null()  !< spatially varying unresolved background gustiness [R L Z T-1 ~> Pa]
                                        !! gust is used when read_gust_2d is true.

  real, pointer :: t_restore(:,:)    => null()  !< temperature to damp (restore) the SST to [degC]
  real, pointer :: s_restore(:,:)    => null()  !< salinity to damp (restore) the SSS [ppt]
  real, pointer :: dens_restore(:,:) => null()  !< density to damp (restore) surface density [R ~> kg m-3]

  integer :: buoy_last_lev_read = -1 !< The last time level read from buoyancy input files

  ! if WIND_CONFIG=='gyres' then use the following as  = A, B, C and n respectively for
  ! taux = A + B*sin(n*pi*y/L) + C*cos(n*pi*y/L)
  real :: gyres_taux_const   !< A constant wind stress [Pa].
  real :: gyres_taux_sin_amp !< The amplitude of cosine wind stress gyres [Pa], if WIND_CONFIG=='gyres'.
  real :: gyres_taux_cos_amp !< The amplitude of cosine wind stress gyres [Pa], if WIND_CONFIG=='gyres'.
  real :: gyres_taux_n_pis   !< The number of sine lobes in the basin if  if WIND_CONFIG=='gyres'
  logical :: answers_2018    !< If true, use the order of arithmetic and expressions that recover
                             !! the answers from the end of 2018.  Otherwise, use a form of the gyre
                             !! wind stresses that are rotationally invariant and more likely to be
                             !! the same between compilers.
  logical :: fix_ustar_gustless_bug         !< If true correct a bug in the time-averaging of the
                                            !! gustless wind friction velocity.
  ! if WIND_CONFIG=='scurves' then use the following to define a piecwise scurve profile
  real :: scurves_ydata(20) = 90. !< Latitudes of scurve nodes [degreesN]
  real :: scurves_taux(20) = 0.   !< Zonal wind stress values at scurve nodes [Pa]

  real :: t_north   !< target temperatures at north used in buoyancy_forcing_linear
  real :: t_south   !< target temperatures at south used in buoyancy_forcing_linear
  real :: s_north   !< target salinity at north used in buoyancy_forcing_linear
  real :: s_south   !< target salinity at south used in buoyancy_forcing_linear

  logical :: first_call_set_forcing = .true. !< True until after the first call to set_forcing
  logical :: archaic_omip_file = .true. !< If true use the variable names and data fields from
                                        !! a very old version of the OMIP forcing
  logical :: dataoverrideisinitialized = .false. !< If true, data override has been initialized

  real :: wind_scale          !< value by which wind-stresses are scaled, ND.
  real :: constantheatforcing !< value used for sensible heat flux when buoy_config="const" [Q R Z T-1 ~> W m-2]

  character(len=8)   :: wind_stagger !< A character indicating how the wind stress components
                              !! are staggered in WIND_FILE.  Valid values are A or C for now.
  type(tracer_flow_control_cs), pointer :: tracer_flow_csp => null() !< A pointer to the structure
                              !! that is used to orchestrate the calling of tracer packages
!#CTRL#  type(ctrl_forcing_CS), pointer :: ctrl_forcing_CSp => NULL()
  type(mom_restart_cs), pointer :: restart_csp => null() !< A pointer to the restart control structure

  type(diag_ctrl), pointer :: diag !< structure used to regulate timing of diagnostic output

  character(len=200) :: inputdir    !< directory where NetCDF input files are.
  character(len=200) :: wind_config !< indicator for wind forcing type (2gyre, USER, FILE..)
  character(len=200) :: wind_file   !< if wind_config is "file", file to use
  character(len=200) :: buoy_config !< indicator for buoyancy forcing type

  character(len=200) :: longwave_file     = '' !< The file from which the longwave heat flux is read
  character(len=200) :: shortwave_file    = '' !< The file from which the shortwave heat flux is read
  character(len=200) :: evaporation_file  = '' !< The file from which the evaporation is read
  character(len=200) :: sensibleheat_file = '' !< The file from which the sensible heat flux is read
  character(len=200) :: latentheat_file   = '' !< The file from which the latent heat flux is read

  character(len=200) :: rain_file   = '' !< The file from which the rainfall is read
  character(len=200) :: snow_file   = '' !< The file from which the snowfall is read
  character(len=200) :: runoff_file = '' !< The file from which the runoff is read

  character(len=200) :: longwaveup_file  = '' !< The file from which the upward longwave heat flux is read
  character(len=200) :: shortwaveup_file = '' !< The file from which the upward shorwave heat flux is read

  character(len=200) :: sstrestore_file      = '' !< The file from which to read the sea surface
                                                  !! temperature to restore toward
  character(len=200) :: salinityrestore_file = '' !< The file from which to read the sea surface
                                                  !! salinity to restore toward

  character(len=80)  :: stress_x_var  = '' !< X-windstress variable name in the input file
  character(len=80)  :: stress_y_var  = '' !< Y-windstress variable name in the input file
  character(len=80)  :: ustar_var     = '' !< ustar variable name in the input file
  character(len=80)  :: lw_var        = '' !< lonngwave heat flux variable name in the input file
  character(len=80)  :: sw_var        = '' !< shortwave heat flux variable name in the input file
  character(len=80)  :: latent_var    = '' !< latent heat flux variable name in the input file
  character(len=80)  :: sens_var      = '' !< sensible heat flux variable name in the input file
  character(len=80)  :: evap_var      = '' !< evaporation variable name in the input file
  character(len=80)  :: rain_var      = '' !< rainfall variable name in the input file
  character(len=80)  :: snow_var      = '' !< snowfall variable name in the input file
  character(len=80)  :: lrunoff_var   = '' !< liquid runoff variable name in the input file
  character(len=80)  :: frunoff_var   = '' !< frozen runoff variable name in the input file
  character(len=80)  :: sst_restore_var = '' !< target sea surface temeperature variable name in the input file
  character(len=80)  :: sss_restore_var = '' !< target sea surface salinity variable name in the input file

  ! These variables give the number of time levels in the various forcing files.
  integer :: wind_nlev = -1   !< The number of time levels in the file of wind stress
  integer :: sw_nlev   = -1   !< The number of time levels in the file of shortwave heat flux
  integer :: lw_nlev = -1     !< The number of time levels in the file of longwave heat flux
  integer :: latent_nlev = -1 !< The number of time levels in the file of latent heat flux
  integer :: sens_nlev = -1   !< The number of time levels in the file of sensible heat flux
  integer :: evap_nlev = -1   !< The number of time levels in the file of evaporation
  integer :: precip_nlev = -1 !< The number of time levels in the file of precipitation
  integer :: runoff_nlev = -1 !< The number of time levels in the file of runoff
  integer :: sst_nlev  = -1   !< The number of time levels in the file of target SST
  integer :: sss_nlev = -1    !< The number of time levels in the file of target SSS

  ! These variables give the last time level read for the various forcing files.
  integer :: wind_last_lev = -1   !< The last time level read of wind stress
  integer :: sw_last_lev   = -1   !< The last time level read of shortwave heat flux
  integer :: lw_last_lev = -1     !< The last time level read of longwave heat flux
  integer :: latent_last_lev = -1 !< The last time level read of latent heat flux
  integer :: sens_last_lev = -1   !< The last time level read of sensible heat flux
  integer :: evap_last_lev = -1   !< The last time level read of evaporation
  integer :: precip_last_lev = -1 !< The last time level read of precipitation
  integer :: runoff_last_lev = -1 !< The last time level read of runoff
  integer :: sst_last_lev  = -1   !< The last time level read of target SST
  integer :: sss_last_lev = -1    !< The last time level read of target SSS

  type(forcing_diags), public :: handles !< A structure with diagnostics handles

  !>@{ Control structures for named forcing packages
  type(user_revise_forcing_cs),  pointer :: urf_cs => null()
  type(user_surface_forcing_cs), pointer :: user_forcing_csp => null()
  type(bfb_surface_forcing_cs), pointer :: bfb_forcing_csp => null()
  type(dumbbell_surface_forcing_cs), pointer :: dumbbell_forcing_csp => null()
  type(meso_surface_forcing_cs), pointer :: meso_forcing_csp => null()
  type(idealized_hurricane_cs), pointer :: idealized_hurricane_csp => null()
  type(scm_cvmix_tests_cs),      pointer :: scm_cvmix_tests_csp => null()
  !>@}

end type surface_forcing_cs

integer :: id_clock_forcing !< A CPU time clock

contains

!> Calls subroutines in this file to get surface forcing fields.
!!
!! It also allocates and initializes the fields in the forcing and mech_forcing types
!! the first time it is called.
subroutine set_forcing(sfc_state, forces, fluxes, day_start, day_interval, 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(forcing),         intent(inout) :: fluxes !< A structure containing thermodynamic forcing fields
  type(time_type),       intent(in)    :: day_start !< The start time of the fluxes
  type(time_type),       intent(in)    :: day_interval !< Length of time over which these fluxes applied
  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(surface_forcing_cs), pointer    :: cs   !< pointer to control struct returned by
                                               !! a previous surface_forcing_init call
  ! Local variables
  real :: dt                     ! length of time over which fluxes applied [s]
  type(time_type) :: day_center  ! central time of the fluxes.
  integer :: isd, ied, jsd, jed
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed

  call cpu_clock_begin(id_clock_forcing)
  call calltree_enter("set_forcing, MOM_surface_forcing.F90")

  day_center = day_start + day_interval/2
  dt = time_type_to_real(day_interval)

  if (cs%first_call_set_forcing) then
    ! Allocate memory for the mechanical and thermodyanmic forcing fields.
    call allocate_mech_forcing(g, forces, stress=.true., ustar=.true., press=.true.)

    call allocate_forcing_type(g, fluxes, ustar=.true., fix_accum_bug=cs%fix_ustar_gustless_bug)
    if (trim(cs%buoy_config) /= "NONE") then
      if ( cs%use_temperature ) then
        call allocate_forcing_type(g, fluxes, water=.true., heat=.true., press=.true.)
        if (cs%restorebuoy) then
          call safe_alloc_ptr(cs%T_Restore,isd, ied, jsd, jed)
          call safe_alloc_ptr(fluxes%heat_added, isd, ied, jsd, jed)
          call safe_alloc_ptr(cs%S_Restore, isd, ied, jsd, jed)
        endif
      else ! CS%use_temperature false.
        call safe_alloc_ptr(fluxes%buoy, isd, ied, jsd, jed)

        if (cs%restorebuoy) call safe_alloc_ptr(cs%Dens_Restore, isd, ied, jsd, jed)
      endif  ! endif for  CS%use_temperature
    endif
  endif

  ! calls to various wind options
  if (cs%variable_winds .or. cs%first_call_set_forcing) then

    if (trim(cs%wind_config) == "file") then
      call wind_forcing_from_file(sfc_state, forces, day_center, g, us, cs)
    elseif (trim(cs%wind_config) == "data_override") then
      call wind_forcing_by_data_override(sfc_state, forces, day_center, g, us, cs)
    elseif (trim(cs%wind_config) == "2gyre") then
      call wind_forcing_2gyre(sfc_state, forces, day_center, g, us, cs)
    elseif (trim(cs%wind_config) == "1gyre") then
      call wind_forcing_1gyre(sfc_state, forces, day_center, g, us, cs)
    elseif (trim(cs%wind_config) == "gyres") then
      call wind_forcing_gyres(sfc_state, forces, day_center, g, us, cs)
    elseif (trim(cs%wind_config) == "zero") then
      call wind_forcing_const(sfc_state, forces, 0., 0., day_center, g, us, cs)
    elseif (trim(cs%wind_config) == "const") then
      call wind_forcing_const(sfc_state, forces, cs%tau_x0, cs%tau_y0, day_center, g, us, cs)
    elseif (trim(cs%wind_config) == "Neverworld" .or. trim(cs%wind_config) == "Neverland") then
      call neverworld_wind_forcing(sfc_state, forces, day_center, g, us, cs)
    elseif (trim(cs%wind_config) == "scurves") then
      call scurve_wind_forcing(sfc_state, forces, day_center, g, us, cs)
    elseif (trim(cs%wind_config) == "ideal_hurr") then
      call idealized_hurricane_wind_forcing(sfc_state, forces, day_center, g, us, cs%idealized_hurricane_CSp)
    elseif (trim(cs%wind_config) == "SCM_ideal_hurr") then
      call scm_idealized_hurricane_wind_forcing(sfc_state, forces, day_center, g, us, cs%idealized_hurricane_CSp)
    elseif (trim(cs%wind_config) == "SCM_CVmix_tests") then
      call scm_cvmix_tests_wind_forcing(sfc_state, forces, day_center, g, us, cs%SCM_CVmix_tests_CSp)
    elseif (trim(cs%wind_config) == "USER") then
      call user_wind_forcing(sfc_state, forces, day_center, g, us, cs%user_forcing_CSp)
    elseif (cs%variable_winds .and. .not.cs%first_call_set_forcing) then
      call mom_error(fatal, &
       "MOM_surface_forcing: Variable winds defined with no wind config")
    else
       call mom_error(fatal, &
       "MOM_surface_forcing:Unrecognized wind config "//trim(cs%wind_config))
    endif
  endif

  ! calls to various buoyancy forcing options
  if (cs%restorebuoy .and. .not.cs%variable_buoyforce) then
    call mom_error(fatal, "With RESTOREBUOY = True, VARIABLE_BUOYFORCE = True should be used. "//&
                          "Otherwise, this can lead to diverging solutions when a simulation "//&
                          "is continued using a restart file.")
  endif

  if ((cs%variable_buoyforce .or. cs%first_call_set_forcing) .and. &
      (.not.cs%adiabatic)) then
    if (trim(cs%buoy_config) == "file") then
      call buoyancy_forcing_from_files(sfc_state, fluxes, day_center, dt, g, us, cs)
    elseif (trim(cs%buoy_config) == "data_override") then
      call buoyancy_forcing_from_data_override(sfc_state, fluxes, day_center, dt, g, us, cs)
    elseif (trim(cs%buoy_config) == "zero") then
      call buoyancy_forcing_zero(sfc_state, fluxes, day_center, dt, g, cs)
    elseif (trim(cs%buoy_config) == "const") then
      call buoyancy_forcing_const(sfc_state, fluxes, day_center, dt, g, us, cs)
    elseif (trim(cs%buoy_config) == "linear") then
      call buoyancy_forcing_linear(sfc_state, fluxes, day_center, dt, g, us, cs)
    elseif (trim(cs%buoy_config) == "MESO") then
      call meso_buoyancy_forcing(sfc_state, fluxes, day_center, dt, g, us, cs%MESO_forcing_CSp)
    elseif (trim(cs%buoy_config) == "SCM_CVmix_tests") then
      call scm_cvmix_tests_buoyancy_forcing(sfc_state, fluxes, day_center, g, us, cs%SCM_CVmix_tests_CSp)
    elseif (trim(cs%buoy_config) == "USER") then
      call user_buoyancy_forcing(sfc_state, fluxes, day_center, dt, g, us, cs%user_forcing_CSp)
    elseif (trim(cs%buoy_config) == "BFB") then
      call bfb_buoyancy_forcing(sfc_state, fluxes, day_center, dt, g, us, cs%BFB_forcing_CSp)
    elseif (trim(cs%buoy_config) == "dumbbell") then
      call dumbbell_buoyancy_forcing(sfc_state, fluxes, day_center, dt, g, us, cs%dumbbell_forcing_CSp)
    elseif (trim(cs%buoy_config) == "NONE") then
      call mom_mesg("MOM_surface_forcing: buoyancy forcing has been set to omitted.")
    elseif (cs%variable_buoyforce .and. .not.cs%first_call_set_forcing) then
      call mom_error(fatal, &
       "MOM_surface_forcing: Variable buoy defined with no buoy config.")
    else
       call mom_error(fatal, &
       "MOM_surface_forcing: Unrecognized buoy config "//trim(cs%buoy_config))
    endif
  endif

  if (associated(cs%tracer_flow_CSp)) then
    call call_tracer_set_forcing(sfc_state, fluxes, day_start, day_interval, g, cs%tracer_flow_CSp)
  endif

  ! Allow for user-written code to alter the fluxes after all the above
  call user_alter_forcing(sfc_state, fluxes, day_center, g, cs%urf_CS)

  ! Fields that exist in both the forcing and mech_forcing types must be copied.
  if (cs%variable_winds .or. cs%first_call_set_forcing) then
    call copy_common_forcing_fields(forces, fluxes, g)
    call set_derived_forcing_fields(forces, fluxes, g, us, cs%Rho0)
  endif

  if ((cs%variable_buoyforce .or. cs%first_call_set_forcing) .and. &
      (.not.cs%adiabatic)) then
    call set_net_mass_forcing(fluxes, forces, g, us)
  endif

  cs%first_call_set_forcing = .false.

  call cpu_clock_end(id_clock_forcing)
  call calltree_leave("set_forcing")

end subroutine set_forcing

!> Sets the surface wind stresses to constant values
subroutine wind_forcing_const(sfc_state, forces, tau_x0, tau_y0, 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
  real,                     intent(in)    :: tau_x0 !< The zonal wind stress [Pa]
  real,                     intent(in)    :: tau_y0 !< The meridional wind stress [Pa]
  type(time_type),          intent(in)    :: day  !< The time of the fluxes
  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(surface_forcing_CS), pointer       :: CS   !< pointer to control struct returned by
                                                  !! a previous surface_forcing_init call
  ! Local variables
  real :: Pa_conversion ! A unit conversion factor from Pa to the internal units [R Z L T-2 Pa-1 ~> 1]
  real :: mag_tau
  integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq

  call calltree_enter("wind_forcing_const, MOM_surface_forcing.F90")
  is   = g%isc  ; ie   = g%iec  ; js   = g%jsc  ; je   = g%jec
  isq  = g%IscB ; ieq  = g%IecB ; jsq  = g%JscB ; jeq  = g%JecB
  pa_conversion = us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z

  !set steady surface wind stresses, in units of Pa.
  mag_tau = pa_conversion * sqrt( tau_x0**2 + tau_y0**2)

  do j=js,je ; do i=is-1,ieq
    forces%taux(i,j) = tau_x0 * pa_conversion
  enddo ; enddo

  do j=js-1,jeq ; do i=is,ie
    forces%tauy(i,j) = tau_y0 * pa_conversion
  enddo ; enddo

  if (cs%read_gust_2d) then
    if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie
      forces%ustar(i,j) = sqrt( us%L_to_Z * ( mag_tau + cs%gust(i,j) ) / cs%Rho0 )
    enddo ; enddo ; endif
  else
    if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie
      forces%ustar(i,j) = sqrt( us%L_to_Z * ( mag_tau + cs%gust_const ) / cs%Rho0 )
    enddo ; enddo ; endif
  endif

  call calltree_leave("wind_forcing_const")
end subroutine wind_forcing_const


!> Sets the surface wind stresses to set up two idealized gyres.
subroutine wind_forcing_2gyre(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(in)    :: G    !< The ocean's grid structure
  type(unit_scale_type),    intent(in)    :: US   !< A dimensional unit scaling type
  type(surface_forcing_CS), pointer       :: CS   !< pointer to control struct returned by
                                                  !! a previous surface_forcing_init call
  ! Local variables
  real :: PI
  integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq

  call calltree_enter("wind_forcing_2gyre, MOM_surface_forcing.F90")
  is   = g%isc  ; ie   = g%iec  ; js   = g%jsc  ; je   = g%jec
  isq  = g%IscB ; ieq  = g%IecB ; jsq  = g%JscB ; jeq  = g%JecB

  !set the steady surface wind stresses, in units of Pa.
  pi = 4.0*atan(1.0)

  do j=js,je ; do i=is-1,ieq
    forces%taux(i,j) = 0.1*us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z * &
                      (1.0 - cos(2.0*pi*(g%geoLatCu(i,j)-cs%South_lat) / cs%len_lat))
  enddo ; enddo

  do j=js-1,jeq ; do i=is,ie
    forces%tauy(i,j) = 0.0
  enddo ; enddo

  call calltree_leave("wind_forcing_2gyre")
end subroutine wind_forcing_2gyre


!> Sets the surface wind stresses to set up a single idealized gyre.
subroutine wind_forcing_1gyre(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(in)    :: G    !< The ocean's grid structure
  type(unit_scale_type),    intent(in)    :: US   !< A dimensional unit scaling type
  type(surface_forcing_CS), pointer       :: CS   !< pointer to control struct returned by
                                                  !! a previous surface_forcing_init call
  ! Local variables
  real :: PI
  integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq

  call calltree_enter("wind_forcing_1gyre, MOM_surface_forcing.F90")
  is   = g%isc  ; ie   = g%iec  ; js   = g%jsc  ; je   = g%jec
  isq  = g%IscB ; ieq  = g%IecB ; jsq  = g%JscB ; jeq  = g%JecB

  ! set the steady surface wind stresses, in units of Pa.
  pi = 4.0*atan(1.0)

  do j=js,je ; do i=is-1,ieq
    forces%taux(i,j) = -0.2*us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z * &
                       cos(pi*(g%geoLatCu(i,j)-cs%South_lat)/cs%len_lat)
  enddo ; enddo

  do j=js-1,jeq ; do i=is,ie
    forces%tauy(i,j) = 0.0
  enddo ; enddo

  call calltree_leave("wind_forcing_1gyre")
end subroutine wind_forcing_1gyre

!> Sets the surface wind stresses to set up idealized gyres.
subroutine wind_forcing_gyres(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(in)    :: G    !< The ocean's grid structure
  type(unit_scale_type),    intent(in)    :: US   !< A dimensional unit scaling type
  type(surface_forcing_CS), pointer       :: CS   !< pointer to control struct returned by
                                                  !! a previous surface_forcing_init call
  ! Local variables
  real :: PI, y, I_rho
  integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq

  call calltree_enter("wind_forcing_gyres, MOM_surface_forcing.F90")
  is   = g%isc  ; ie   = g%iec  ; js   = g%jsc  ; je   = g%jec
  isq  = g%IscB ; ieq  = g%IecB ; jsq  = g%JscB ; jeq  = g%JecB

  ! steady surface wind stresses [Pa]
  pi = 4.0*atan(1.0)

  do j=js-1,je+1 ; do i=is-1,ieq
    y = (g%geoLatCu(i,j)-cs%South_lat) / cs%len_lat
    forces%taux(i,j) = us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z * &
            (cs%gyres_taux_const +                            &
             (   cs%gyres_taux_sin_amp*sin(cs%gyres_taux_n_pis*pi*y)    &
               + cs%gyres_taux_cos_amp*cos(cs%gyres_taux_n_pis*pi*y) ))
  enddo ; enddo

  do j=js-1,jeq ; do i=is-1,ie+1
    forces%tauy(i,j) = 0.0
  enddo ; enddo

  ! set the friction velocity
  if (cs%answers_2018) then
    do j=js,je ; do i=is,ie
      forces%ustar(i,j) = sqrt(us%L_to_Z * ((cs%gust_const/cs%Rho0) + &
              sqrt(0.5*(forces%tauy(i,j-1)*forces%tauy(i,j-1) + forces%tauy(i,j)*forces%tauy(i,j) + &
                        forces%taux(i-1,j)*forces%taux(i-1,j) + forces%taux(i,j)*forces%taux(i,j)))/cs%Rho0) )
    enddo ; enddo
  else
    i_rho = us%L_to_Z / cs%Rho0
    do j=js,je ; do i=is,ie
      forces%ustar(i,j) = sqrt( (cs%gust_const + &
            sqrt(0.5*((forces%tauy(i,j-1)**2 + forces%tauy(i,j)**2) + &
                      (forces%taux(i-1,j)**2 + forces%taux(i,j)**2))) ) * i_rho )
    enddo ; enddo
  endif

  call calltree_leave("wind_forcing_gyres")
end subroutine wind_forcing_gyres

!> Sets the surface wind stresses, forces%taux and forces%tauy for the
!! Neverworld forcing configuration.
subroutine neverworld_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    !< Time used for determining the fluxes.
  type(ocean_grid_type),    intent(inout) :: G      !< Grid structure.
  type(unit_scale_type),    intent(in)    :: US     !< A dimensional unit scaling type
  type(surface_forcing_CS), pointer       :: CS     !< pointer to control struct returned by
                                                    !! a previous surface_forcing_init call
  ! Local variables
  integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
  real :: PI, I_rho, y
  real :: tau_max  ! The magnitude of the wind stress [R Z L T-2 ~> Pa]
  real :: off

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

  ! Allocate the forcing arrays, if necessary.
  call allocate_mech_forcing(g, forces, stress=.true.)

  !  Set the surface wind stresses, in units of 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.
  pi = 4.0*atan(1.0)
  forces%taux(:,:) = 0.0
  tau_max = 0.2 * us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z
  off = 0.02
  do j=js,je ; do i=is-1,ieq
    y = (g%geoLatT(i,j)-g%south_lat)/g%len_lat

    if (y <= 0.29) then
      forces%taux(i,j) = forces%taux(i,j) + tau_max * ( (1/0.29)*y - ( 1/(2*pi) )*sin( (2*pi*y) / 0.29 ) )
    endif
    if ((y > 0.29) .and. (y <= (0.8-off))) then
      forces%taux(i,j) = forces%taux(i,j) + tau_max *(0.35+0.65*cos(pi*(y-0.29)/(0.51-off))  )
    endif
    if ((y > (0.8-off)) .and. (y <= (1-off))) then
      forces%taux(i,j) = forces%taux(i,j) + tau_max *( 1.5*( (y-1+off) - (0.1/pi)*sin(10.0*pi*(y-0.8+off)) ) )
    endif
    forces%taux(i,j) = g%mask2dCu(i,j) * forces%taux(i,j)
  enddo ; enddo

  do j=js-1,jeq ; do i=is,ie
    forces%tauy(i,j) = g%mask2dCv(i,j) * 0.0
  enddo ; enddo

  ! Set the surface friction velocity, in units of m s-1.  ustar is always positive.
  if (associated(forces%ustar)) then
    i_rho = us%L_to_Z / cs%Rho0
    do j=js,je ; do i=is,ie
      forces%ustar(i,j) = sqrt( (cs%gust_const + &
            sqrt(0.5*((forces%tauy(i,j-1)**2 + forces%tauy(i,j)**2) + &
                      (forces%taux(i-1,j)**2 + forces%taux(i,j)**2))) ) * i_rho )
    enddo ; enddo
  endif

end subroutine neverworld_wind_forcing

!> Sets the zonal wind stresses to a piecewise series of s-curves.
subroutine scurve_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    !< Time used for determining the fluxes.
  type(ocean_grid_type),    intent(inout) :: G      !< Grid structure.
  type(unit_scale_type),    intent(in)    :: US     !< A dimensional unit scaling type
  type(surface_forcing_CS), pointer       :: CS     !< pointer to control struct returned by
                                                    !! a previous surface_forcing_init call
  ! Local variables
  integer :: i, j, kseg
  real :: lon, lat, I_rho, y, L
! real :: ydata(7) = (/ -70., -45., -15., 0., 15., 45., 70. /)
! real :: taudt(7) = (/ 0., 0.2, -0.1, -0.02, -0.1, 0.1, 0. /)

  ! Allocate the forcing arrays, if necessary.
  call allocate_mech_forcing(g, forces, stress=.true.)

  kseg = 1
  do j=g%jsd,g%jed ; do i=g%IsdB,g%IedB
    lon = g%geoLonCu(i,j)
    lat = g%geoLatCu(i,j)

    ! Find segment k s.t. ydata(k)<= lat < ydata(k+1)
    do while (lat>=cs%scurves_ydata(kseg+1) .and. kseg<6)
      kseg = kseg+1
    enddo
    do while (lat<cs%scurves_ydata(kseg) .and. kseg>1)
      kseg = kseg-1
    enddo

    y = lat - cs%scurves_ydata(kseg)
    l = cs%scurves_ydata(kseg+1) - cs%scurves_ydata(kseg)
    forces%taux(i,j) = cs%scurves_taux(kseg) +  &
                       ( cs%scurves_taux(kseg+1) - cs%scurves_taux(kseg) ) * scurve(y, l)
    forces%taux(i,j) = g%mask2dCu(i,j) * forces%taux(i,j)
  enddo ; enddo

  do j=g%JsdB,g%JedB ; do i=g%isd,g%ied
    forces%tauy(i,j) = g%mask2dCv(i,j) * 0.0
  enddo ; enddo

  ! Set the surface friction velocity, in units of m s-1.  ustar is always positive.
  if (associated(forces%ustar)) then
    i_rho = us%L_to_Z / cs%Rho0
    do j=g%jsc,g%jec ; do i=g%isc,g%iec
      forces%ustar(i,j) = sqrt( (cs%gust_const + &
            sqrt(0.5*((forces%tauy(i,j-1)**2 + forces%tauy(i,j)**2) + &
                      (forces%taux(i-1,j)**2 + forces%taux(i,j)**2))) ) * i_rho )
    enddo ; enddo
  endif

end subroutine scurve_wind_forcing

!> Returns the value of a cosine-bell function evaluated at x/L
real function scurve(x,L)
  real , intent(in) :: x       !< non-dimensional position
  real , intent(in) :: l       !< non-dimensional width
  real :: s

  s = x/l
  scurve = (3. - 2.*s) * (s*s)
end function scurve

! Sets the surface wind stresses from input files.
subroutine wind_forcing_from_file(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(surface_forcing_CS), pointer       :: CS   !< pointer to control struct returned by
                                                  !! a previous surface_forcing_init call
  ! Local variables
  character(len=200) :: filename  ! The name of the input file.
  real    :: temp_x(SZI_(G),SZJ_(G)) ! Pseudo-zonal and psuedo-meridional
  real    :: temp_y(SZI_(G),SZJ_(G)) ! wind stresses at h-points [R L Z T-1 ~> Pa].
  real    :: Pa_conversion           ! A unit conversion factor from Pa to the internal wind stress
                                     ! units [R Z L T-2 Pa-1 ~> 1]
  integer :: time_lev_daily          ! The time levels to read for fields with
  integer :: time_lev_monthly        ! daily and montly cycles.
  integer :: time_lev                ! The time level that is used for a field.
  integer :: days, seconds
  integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
  logical :: read_Ustar

  call calltree_enter("wind_forcing_from_file, MOM_surface_forcing.F90")
  is   = g%isc  ; ie   = g%iec  ; js   = g%jsc  ; je   = g%jec
  isq  = g%IscB ; ieq  = g%IecB ; jsq  = g%JscB ; jeq  = g%JecB
  pa_conversion = us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z

  call get_time(day, seconds, days)
  time_lev_daily = days - 365*floor(real(days) / 365.0)

  if (time_lev_daily < 31) then ; time_lev_monthly = 0
  elseif (time_lev_daily < 59)  then ; time_lev_monthly = 1
  elseif (time_lev_daily < 90)  then ; time_lev_monthly = 2
  elseif (time_lev_daily < 120) then ; time_lev_monthly = 3
  elseif (time_lev_daily < 151) then ; time_lev_monthly = 4
  elseif (time_lev_daily < 181) then ; time_lev_monthly = 5
  elseif (time_lev_daily < 212) then ; time_lev_monthly = 6
  elseif (time_lev_daily < 243) then ; time_lev_monthly = 7
  elseif (time_lev_daily < 273) then ; time_lev_monthly = 8
  elseif (time_lev_daily < 304) then ; time_lev_monthly = 9
  elseif (time_lev_daily < 334) then ; time_lev_monthly = 10
  else ; time_lev_monthly = 11
  endif

  time_lev_daily = time_lev_daily+1
  time_lev_monthly = time_lev_monthly+1

  select case (cs%wind_nlev)
    case (12) ; time_lev = time_lev_monthly
    case (365) ; time_lev = time_lev_daily
    case default ; time_lev = 1
  end select

  if (time_lev /= cs%wind_last_lev) then
    filename = trim(cs%wind_file)
    read_ustar = (len_trim(cs%ustar_var) > 0)
!    if (is_root_pe()) &
!      write(*,'("Wind_forcing Reading time level ",I," last was ",I,".")')&
!           time_lev-1,CS%wind_last_lev-1
    select case ( uppercase(cs%wind_stagger(1:1)) )
    case ("A")
      temp_x(:,:) = 0.0 ; temp_y(:,:) = 0.0
      call mom_read_vector(filename, cs%stress_x_var, cs%stress_y_var, &
                           temp_x(:,:), temp_y(:,:), g%Domain, stagger=agrid, &
                           timelevel=time_lev, scale=pa_conversion)

      call pass_vector(temp_x, temp_y, g%Domain, to_all, agrid)
      do j=js,je ; do i=is-1,ieq
        forces%taux(i,j) = 0.5 * cs%wind_scale * (temp_x(i,j) + temp_x(i+1,j))
      enddo ; enddo
      do j=js-1,jeq ; do i=is,ie
        forces%tauy(i,j) = 0.5 * cs%wind_scale * (temp_y(i,j) + temp_y(i,j+1))
      enddo ; enddo

      if (.not.read_ustar) then
        if (cs%read_gust_2d) then
          do j=js,je ; do i=is,ie
            forces%ustar(i,j) = sqrt((cs%gust(i,j) + &
                    sqrt(temp_x(i,j)*temp_x(i,j) + temp_y(i,j)*temp_y(i,j))) * us%L_to_Z / cs%Rho0)
          enddo ; enddo
        else
          do j=js,je ; do i=is,ie
            forces%ustar(i,j) = sqrt(us%L_to_Z * (cs%gust_const/cs%Rho0 + &
                    sqrt(temp_x(i,j)*temp_x(i,j) + temp_y(i,j)*temp_y(i,j)) / cs%Rho0) )
          enddo ; enddo
        endif
      endif
    case ("C")
      if (g%symmetric) then
        if (.not.associated(g%Domain_aux)) call mom_error(fatal, &
          " wind_forcing_from_file with C-grid input and symmetric memory "//&
          " called with a non-associated auxiliary domain in the grid type.")
        !   Read the data as though symmetric memory were not being used, and
        ! then translate it appropriately.
        temp_x(:,:) = 0.0 ; temp_y(:,:) = 0.0
        call mom_read_vector(filename, cs%stress_x_var, cs%stress_y_var, &
                             temp_x(:,:), temp_y(:,:), &
                             g%Domain_aux, stagger=cgrid_ne, timelevel=time_lev, &
                             scale=pa_conversion)
        do j=js,je ; do i=is,ie
          forces%taux(i,j) = cs%wind_scale * temp_x(i,j)
          forces%tauy(i,j) = cs%wind_scale * temp_y(i,j)
        enddo ; enddo
        call fill_symmetric_edges(forces%taux, forces%tauy, g%Domain, stagger=cgrid_ne)
      else
        call mom_read_vector(filename, cs%stress_x_var, cs%stress_y_var, &
                             forces%taux(:,:), forces%tauy(:,:), &
                             g%Domain, stagger=cgrid_ne, timelevel=time_lev, &
                             scale=pa_conversion)

        if (cs%wind_scale /= 1.0) then
          do j=js,je ; do i=isq,ieq
            forces%taux(i,j) = cs%wind_scale * forces%taux(i,j)
          enddo ; enddo
          do j=jsq,jeq ; do i=is,ie
            forces%tauy(i,j) = cs%wind_scale * forces%tauy(i,j)
          enddo ; enddo
        endif
      endif

      call pass_vector(forces%taux, forces%tauy, g%Domain, to_all)
      if (.not.read_ustar) then
        if (cs%read_gust_2d) then
          do j=js, je ; do i=is, ie
            forces%ustar(i,j) = sqrt((cs%gust(i,j) + &
                    sqrt(0.5*((forces%tauy(i,j-1)**2 + forces%tauy(i,j)**2) + &
                              (forces%taux(i-1,j)**2 + forces%taux(i,j)**2))) ) * us%L_to_Z / cs%Rho0 )
          enddo ; enddo
        else
          do j=js, je ; do i=is, ie
            forces%ustar(i,j) = sqrt(us%L_to_Z * ( (cs%gust_const/cs%Rho0) + &
                    sqrt(0.5*((forces%tauy(i,j-1)**2 + forces%tauy(i,j)**2) + &
                              (forces%taux(i-1,j)**2 + forces%taux(i,j)**2)))/cs%Rho0))
          enddo ; enddo
        endif
      endif
    case default
      call mom_error(fatal, "wind_forcing_from_file: Unrecognized stagger "//&
                      trim(cs%wind_stagger)//" is not 'A' or 'C'.")
    end select

    if (read_ustar) then
      call mom_read_data(filename, cs%Ustar_var, forces%ustar(:,:), &
                         g%Domain, timelevel=time_lev, scale=us%m_to_Z*us%T_to_s)
    endif

    cs%wind_last_lev = time_lev

  endif ! time_lev /= CS%wind_last_lev

  call calltree_leave("wind_forcing_from_file")
end subroutine wind_forcing_from_file


! Sets the surface wind stresses via the data override facility.
subroutine wind_forcing_by_data_override(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(surface_forcing_CS), pointer       :: CS   !< pointer to control struct returned by
                                                  !! a previous surface_forcing_init call
  ! Local variables
  real :: temp_x(SZI_(G),SZJ_(G)) ! Pseudo-zonal and psuedo-meridional
  real :: temp_y(SZI_(G),SZJ_(G)) ! wind stresses at h-points [Pa].
  real :: temp_ustar(SZI_(G),SZJ_(G)) ! ustar [m s-1] (not rescaled).
  real :: Pa_conversion ! A unit conversion factor from Pa to the internal units [R Z L T-2 Pa-1 ~> 1]
  integer :: i, j, is_in, ie_in, js_in, je_in
  logical :: read_uStar

  call calltree_enter("wind_forcing_by_data_override, MOM_surface_forcing.F90")

  if (.not.cs%dataOverrideIsInitialized) then
    call allocate_mech_forcing(g, forces, stress=.true., ustar=.true., press=.true.)
    call data_override_init(ocean_domain_in=g%Domain%mpp_domain)
    cs%dataOverrideIsInitialized = .true.
  endif

  is_in = g%isc - g%isd + 1 ; ie_in = g%iec - g%isd + 1
  js_in = g%jsc - g%jsd + 1 ; je_in = g%jec - g%jsd + 1
  pa_conversion = us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z

  temp_x(:,:) = 0.0 ; temp_y(:,:) = 0.0
  call data_override('OCN', 'taux', temp_x, day, is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)
  call data_override('OCN', 'tauy', temp_y, day, is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)
  call pass_vector(temp_x, temp_y, g%Domain, to_all, agrid)
  ! Ignore CS%wind_scale when using data_override ?????
  do j=g%jsc,g%jec ; do i=g%isc-1,g%IecB
    forces%taux(i,j) = pa_conversion * 0.5 * (temp_x(i,j) + temp_x(i+1,j))
  enddo ; enddo
  do j=g%jsc-1,g%JecB ; do i=g%isc,g%iec
    forces%tauy(i,j) = pa_conversion * 0.5 * (temp_y(i,j) + temp_y(i,j+1))
  enddo ; enddo

  read_ustar = (len_trim(cs%ustar_var) > 0) ! Need better control higher up ????
  if (read_ustar) then
    do j=g%jsc,g%jec ; do i=g%isc,g%iec ; temp_ustar(i,j) = us%Z_to_m*us%s_to_T*forces%ustar(i,j) ; enddo ; enddo
    call data_override('OCN', 'ustar', temp_ustar, day, is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)
    do j=g%jsc,g%jec ; do i=g%isc,g%iec ; forces%ustar(i,j) = us%m_to_Z*us%T_to_s*temp_ustar(i,j) ; enddo ; enddo
  else
    if (cs%read_gust_2d) then
      call data_override('OCN', 'gust', cs%gust, day, is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)
      do j=g%jsc,g%jec ; do i=g%isc,g%iec
        forces%ustar(i,j) = sqrt((pa_conversion * sqrt(temp_x(i,j)*temp_x(i,j) + &
            temp_y(i,j)*temp_y(i,j)) + cs%gust(i,j)) * us%L_to_Z / cs%Rho0)
      enddo ; enddo
    else
      do j=g%jsc,g%jec ; do i=g%isc,g%iec
        forces%ustar(i,j) = sqrt(us%L_to_Z * (pa_conversion*sqrt(temp_x(i,j)*temp_x(i,j) + &
            temp_y(i,j)*temp_y(i,j))/cs%Rho0 + cs%gust_const/cs%Rho0 ))
      enddo ; enddo
    endif
  endif

  call pass_vector(forces%taux, forces%tauy, g%Domain, to_all)
! call pass_var(forces%ustar, G%Domain, To_All)     Not needed  ?????

  call calltree_leave("wind_forcing_by_data_override")
end subroutine wind_forcing_by_data_override


!> Specifies zero surface bouyancy fluxes from input files.
subroutine buoyancy_forcing_from_files(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 [s]
  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(surface_forcing_CS), pointer    :: CS   !< pointer to control struct returned by
                                               !! a previous surface_forcing_init call
  ! Local variables
  real, dimension(SZI_(G),SZJ_(G)) :: &
    temp, &       ! A 2-d temporary work array with various units.
    SST_anom, &   ! Instantaneous sea surface temperature anomalies from a
                  ! target (observed) value [degC].
    sss_anom, &   ! Instantaneous sea surface salinity anomalies from a target
                  ! (observed) value [ppt].
    sss_mean      ! A (mean?) salinity about which to normalize local salinity
                  ! anomalies when calculating restorative precipitation
                  ! anomalies [ppt].

  real :: kg_m2_s_conversion  ! A combination of unit conversion factors for rescaling
                              ! mass fluxes [R Z s m2 kg-1 T-1 ~> 1].
  real :: rhoXcp ! reference density times heat capacity [Q R degC-1 ~> J m-3 degC-1]

  integer :: time_lev_daily     ! time levels to read for fields with daily cycle
  integer :: time_lev_monthly   ! time levels to read for fields with monthly cycle
  integer :: time_lev           ! time level that for a field

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

  call calltree_enter("buoyancy_forcing_from_files, MOM_surface_forcing.F90")

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

  if (cs%use_temperature) rhoxcp = cs%Rho0 * fluxes%C_p

  ! Read the buoyancy forcing file
  call get_time(day, seconds, days)

  time_lev_daily = days - 365*floor(real(days) / 365.0)

  if (time_lev_daily < 31) then ; time_lev_monthly = 0
  elseif (time_lev_daily < 59)  then ; time_lev_monthly = 1
  elseif (time_lev_daily < 90)  then ; time_lev_monthly = 2
  elseif (time_lev_daily < 120) then ; time_lev_monthly = 3
  elseif (time_lev_daily < 151) then ; time_lev_monthly = 4
  elseif (time_lev_daily < 181) then ; time_lev_monthly = 5
  elseif (time_lev_daily < 212) then ; time_lev_monthly = 6
  elseif (time_lev_daily < 243) then ; time_lev_monthly = 7
  elseif (time_lev_daily < 273) then ; time_lev_monthly = 8
  elseif (time_lev_daily < 304) then ; time_lev_monthly = 9
  elseif (time_lev_daily < 334) then ; time_lev_monthly = 10
  else ; time_lev_monthly = 11
  endif

  time_lev_daily   = time_lev_daily  +1
  time_lev_monthly = time_lev_monthly+1

  if (time_lev_daily /= cs%buoy_last_lev_read) then

    ! longwave
    select case (cs%LW_nlev)
      case (12)    ; time_lev = time_lev_monthly
      case (365)   ; time_lev = time_lev_daily
      case default ; time_lev = 1
    end select
    call mom_read_data(cs%longwave_file, cs%LW_var, fluxes%lw(:,:), &
                       g%Domain, timelevel=time_lev, scale=us%W_m2_to_QRZ_T)
    if (cs%archaic_OMIP_file) then
      call mom_read_data(cs%longwaveup_file, "lwup_sfc", temp(:,:), g%Domain, &
                         timelevel=time_lev, scale=us%W_m2_to_QRZ_T)
      do j=js,je ; do i=is,ie ; fluxes%LW(i,j) = fluxes%LW(i,j) - temp(i,j) ; enddo ; enddo
    endif
    cs%LW_last_lev = time_lev

    ! evaporation
    select case (cs%evap_nlev)
      case (12)    ; time_lev = time_lev_monthly
      case (365)   ; time_lev = time_lev_daily
      case default ; time_lev = 1
    end select
    if (cs%archaic_OMIP_file) then
      call mom_read_data(cs%evaporation_file, cs%evap_var, fluxes%evap(:,:), &
                     g%Domain, timelevel=time_lev, scale=-kg_m2_s_conversion)
      do j=js,je ; do i=is,ie
        fluxes%latent(i,j)           = cs%latent_heat_vapor*fluxes%evap(i,j)
        fluxes%latent_evap_diag(i,j) = fluxes%latent(i,j)
      enddo ; enddo
    else
      call mom_read_data(cs%evaporation_file, cs%evap_var, fluxes%evap(:,:), &
                     g%Domain, timelevel=time_lev, scale=kg_m2_s_conversion)
    endif
    cs%evap_last_lev = time_lev

    select case (cs%latent_nlev)
      case (12)    ; time_lev = time_lev_monthly
      case (365)   ; time_lev = time_lev_daily
      case default ; time_lev = 1
    end select
    if (.not.cs%archaic_OMIP_file) then
      call mom_read_data(cs%latentheat_file, cs%latent_var, fluxes%latent(:,:), &
                     g%Domain, timelevel=time_lev, scale=us%W_m2_to_QRZ_T)
      do j=js,je ; do i=is,ie
        fluxes%latent_evap_diag(i,j) = fluxes%latent(i,j)
      enddo ; enddo
    endif
    cs%latent_last_lev = time_lev

    select case (cs%sens_nlev)
      case (12)    ; time_lev = time_lev_monthly
      case (365)   ; time_lev = time_lev_daily
      case default ; time_lev = 1
    end select
    if (cs%archaic_OMIP_file) then
      call mom_read_data(cs%sensibleheat_file, cs%sens_var, fluxes%sens(:,:), &
                     g%Domain, timelevel=time_lev, scale=-us%W_m2_to_QRZ_T)
    else
      call mom_read_data(cs%sensibleheat_file, cs%sens_var, fluxes%sens(:,:), &
                     g%Domain, timelevel=time_lev, scale=us%W_m2_to_QRZ_T)
    endif
    cs%sens_last_lev = time_lev

    select case (cs%SW_nlev)
      case (12)    ; time_lev = time_lev_monthly
      case (365)   ; time_lev = time_lev_daily
      case default ; time_lev = 1
    end select
    call mom_read_data(cs%shortwave_file, cs%SW_var, fluxes%sw(:,:), g%Domain, &
                       timelevel=time_lev, scale=us%W_m2_to_QRZ_T)
    if (cs%archaic_OMIP_file) then
      call mom_read_data(cs%shortwaveup_file, "swup_sfc", temp(:,:), g%Domain, &
                         timelevel=time_lev, scale=us%W_m2_to_QRZ_T)
      do j=js,je ; do i=is,ie
        fluxes%sw(i,j) = fluxes%sw(i,j) - temp(i,j)
      enddo ; enddo
    endif
    cs%SW_last_lev = time_lev

    select case (cs%precip_nlev)
      case (12)    ; time_lev = time_lev_monthly
      case (365)   ; time_lev = time_lev_daily
      case default ; time_lev = 1
    end select
    call mom_read_data(cs%snow_file, cs%snow_var, &
             fluxes%fprec(:,:), g%Domain, timelevel=time_lev, scale=kg_m2_s_conversion)
    call mom_read_data(cs%rain_file, cs%rain_var, &
             fluxes%lprec(:,:), g%Domain, timelevel=time_lev, scale=kg_m2_s_conversion)
    if (cs%archaic_OMIP_file) then
      do j=js,je ; do i=is,ie
        fluxes%lprec(i,j) = fluxes%lprec(i,j) - fluxes%fprec(i,j)
      enddo ; enddo
    endif
    cs%precip_last_lev = time_lev

    select case (cs%runoff_nlev)
      case (12)    ; time_lev = time_lev_monthly
      case (365)   ; time_lev = time_lev_daily
      case default ; time_lev = 1
    end select
    if (cs%archaic_OMIP_file) then
      call mom_read_data(cs%runoff_file, cs%lrunoff_var, temp(:,:), &
                     g%Domain, timelevel=time_lev, scale=kg_m2_s_conversion)
      do j=js,je ; do i=is,ie
        fluxes%lrunoff(i,j) = temp(i,j)*us%m_to_L**2*g%IareaT(i,j)
      enddo ; enddo
      call mom_read_data(cs%runoff_file, cs%frunoff_var, temp(:,:), &
                     g%Domain, timelevel=time_lev, scale=kg_m2_s_conversion)
      do j=js,je ; do i=is,ie
        fluxes%frunoff(i,j) = temp(i,j)*us%m_to_L**2*g%IareaT(i,j)
      enddo ; enddo
    else
      call mom_read_data(cs%runoff_file, cs%lrunoff_var, fluxes%lrunoff(:,:), &
                     g%Domain, timelevel=time_lev, scale=kg_m2_s_conversion)
      call mom_read_data(cs%runoff_file, cs%frunoff_var, fluxes%frunoff(:,:), &
                     g%Domain, timelevel=time_lev, scale=kg_m2_s_conversion)
    endif
    cs%runoff_last_lev = time_lev

!     Read the SST and SSS fields for damping.
    if (cs%restorebuoy) then !#CTRL# .or. associated(CS%ctrl_forcing_CSp)) then
      select case (cs%SST_nlev)
        case (12)    ; time_lev = time_lev_monthly
        case (365)   ; time_lev = time_lev_daily
        case default ; time_lev = 1
      end select
      call mom_read_data(cs%SSTrestore_file, cs%SST_restore_var, &
               cs%T_Restore(:,:), g%Domain, timelevel=time_lev)
      cs%SST_last_lev = time_lev

      select case (cs%SSS_nlev)
        case (12)    ; time_lev = time_lev_monthly
        case (365)   ; time_lev = time_lev_daily
        case default ; time_lev = 1
      end select
      call mom_read_data(cs%salinityrestore_file, cs%SSS_restore_var, &
               cs%S_Restore(:,:), g%Domain, timelevel=time_lev)
      cs%SSS_last_lev = time_lev
    endif
    cs%buoy_last_lev_read = time_lev_daily

    ! mask out land points and compute heat content of water fluxes
    ! assume liquid precip enters ocean at SST
    ! assume frozen precip enters ocean at 0degC
    ! assume liquid runoff enters ocean at SST
    ! assume solid runoff (calving) enters ocean at 0degC
    ! mass leaving the ocean has heat_content determined in MOM_diabatic_driver.F90
    do j=js,je ; do i=is,ie
      fluxes%evap(i,j)    = fluxes%evap(i,j)    * g%mask2dT(i,j)
      fluxes%lprec(i,j)   = fluxes%lprec(i,j)   * g%mask2dT(i,j)
      fluxes%fprec(i,j)   = fluxes%fprec(i,j)   * g%mask2dT(i,j)
      fluxes%lrunoff(i,j) = fluxes%lrunoff(i,j) * g%mask2dT(i,j)
      fluxes%frunoff(i,j) = fluxes%frunoff(i,j) * g%mask2dT(i,j)
      fluxes%lw(i,j)      = fluxes%lw(i,j)      * g%mask2dT(i,j)
      fluxes%sens(i,j)    = fluxes%sens(i,j)    * g%mask2dT(i,j)
      fluxes%sw(i,j)      = fluxes%sw(i,j)      * g%mask2dT(i,j)
      fluxes%latent(i,j)  = fluxes%latent(i,j)  * g%mask2dT(i,j)

      fluxes%latent_evap_diag(i,j)     = fluxes%latent_evap_diag(i,j) * g%mask2dT(i,j)
      fluxes%latent_fprec_diag(i,j)    = -fluxes%fprec(i,j)*cs%latent_heat_fusion
      fluxes%latent_frunoff_diag(i,j)  = -fluxes%frunoff(i,j)*cs%latent_heat_fusion
    enddo ; enddo

  endif ! time_lev /= CS%buoy_last_lev_read


  ! restoring surface boundary fluxes
  if (cs%restorebuoy) then

    if (cs%use_temperature) then
      do j=js,je ; do i=is,ie
        if (g%mask2dT(i,j) > 0) then
          fluxes%heat_added(i,j) = g%mask2dT(i,j) * &
              ((cs%T_Restore(i,j) - sfc_state%SST(i,j)) * rhoxcp * cs%Flux_const_T)
          fluxes%vprec(i,j) = - (cs%Rho0*cs%Flux_const_S) * &
              (cs%S_Restore(i,j) - sfc_state%SSS(i,j)) / &
              (0.5*(sfc_state%SSS(i,j) + cs%S_Restore(i,j)))
        else
          fluxes%heat_added(i,j) = 0.0
          fluxes%vprec(i,j)      = 0.0
        endif
      enddo ; enddo
    else
      do j=js,je ; do i=is,ie
        if (g%mask2dT(i,j) > 0) then
          fluxes%buoy(i,j) = (cs%Dens_Restore(i,j) - sfc_state%sfc_density(i,j)) * &
                             (cs%G_Earth * cs%Flux_const / cs%Rho0)
        else
          fluxes%buoy(i,j) = 0.0
        endif
      enddo ; enddo
    endif

  else                                              ! not RESTOREBUOY
    if (.not.cs%use_temperature) then
      call mom_error(fatal, "buoyancy_forcing in MOM_surface_forcing: "// &
                     "The fluxes need to be defined without RESTOREBUOY.")
    endif

  endif                                             ! end RESTOREBUOY

!#CTRL# if (associated(CS%ctrl_forcing_CSp)) then
!#CTRL#   do j=js,je ; do i=is,ie
!#CTRL#     SST_anom(i,j) = sfc_state%SST(i,j) - CS%T_Restore(i,j)
!#CTRL#     SSS_anom(i,j) = sfc_state%SSS(i,j) - CS%S_Restore(i,j)
!#CTRL#     SSS_mean(i,j) = 0.5*(sfc_state%SSS(i,j) + CS%S_Restore(i,j))
!#CTRL#   enddo ; enddo
!#CTRL#   call apply_ctrl_forcing(SST_anom, SSS_anom, SSS_mean, fluxes%heat_added, &
!#CTRL#                           fluxes%vprec, day, dt, G, CS%ctrl_forcing_CSp)
!#CTRL# endif

  call calltree_leave("buoyancy_forcing_from_files")
end subroutine buoyancy_forcing_from_files

!> Specifies zero surface bouyancy fluxes from data over-ride.
subroutine buoyancy_forcing_from_data_override(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 [s]
  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(surface_forcing_CS), pointer       :: CS   !< pointer to control struct returned by
                                                  !! a previous surface_forcing_init call
  ! Local variables
  real, dimension(SZI_(G),SZJ_(G)) :: &
    temp, &       ! A 2-d temporary work array with various units.
    SST_anom, &   ! Instantaneous sea surface temperature anomalies from a
                  ! target (observed) value [degC].
    sss_anom, &   ! Instantaneous sea surface salinity anomalies from a target
                  ! (observed) value [ppt].
    sss_mean      ! A (mean?) salinity about which to normalize local salinity
                  ! anomalies when calculating restorative precipitation
                  ! anomalies [ppt].
  real :: kg_m2_s_conversion  ! A combination of unit conversion factors for rescaling
                              ! mass fluxes [R Z s m2 kg-1 T-1 ~> 1].
  real :: rhoXcp ! The mean density times the heat capacity [Q R degC-1 ~> J m-3 degC-1].

  integer :: time_lev_daily     ! The time levels to read for fields with
  integer :: time_lev_monthly   ! daily and montly cycles.
  integer :: itime_lev           ! The time level that is used for a field.

  integer :: days, seconds
  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
  integer :: is_in, ie_in, js_in, je_in

  call calltree_enter("buoyancy_forcing_from_data_override, MOM_surface_forcing.F90")

  is  = g%isc ; ie  = g%iec ; js  = g%jsc ; je  = g%jec
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  kg_m2_s_conversion = us%kg_m2s_to_RZ_T

  if (cs%use_temperature) rhoxcp = cs%Rho0 * fluxes%C_p

  if (.not.cs%dataOverrideIsInitialized) then
    call data_override_init(ocean_domain_in=g%Domain%mpp_domain)
    cs%dataOverrideIsInitialized = .true.
  endif

  is_in = g%isc - g%isd + 1
  ie_in = g%iec - g%isd + 1
  js_in = g%jsc - g%jsd + 1
  je_in = g%jec - g%jsd + 1

  call data_override('OCN', 'lw', fluxes%lw(:,:), day, &
       is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in) ! scale=US%W_m2_to_QRZ_T
  if (us%QRZ_T_to_W_m2 /= 1.0) then ; do j=js,je ; do i=is,ie
    fluxes%lw(i,j) = fluxes%lw(i,j) * us%W_m2_to_QRZ_T
  enddo ; enddo ; endif
  call data_override('OCN', 'evap', fluxes%evap(:,:), day, &
       is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)

  ! note the sign convention
  do j=js,je ; do i=is,ie
    ! The normal convention is that fluxes%evap positive into the ocean
    ! but evap is normally a positive quantity in the files
    ! This conversion is dangerous because it is not clear whether the data files have been read!
    fluxes%evap(i,j) = -kg_m2_s_conversion*fluxes%evap(i,j)
    fluxes%latent(i,j)           = cs%latent_heat_vapor*fluxes%evap(i,j)
    fluxes%latent_evap_diag(i,j) = fluxes%latent(i,j)
  enddo ; enddo

  call data_override('OCN', 'sens', fluxes%sens(:,:), day, &
       is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)

  ! note the sign convention
  do j=js,je ; do i=is,ie
     fluxes%sens(i,j) = -us%W_m2_to_QRZ_T * fluxes%sens(i,j)  ! Normal convention is positive into the ocean
                                           ! but sensible is normally a positive quantity in the files
  enddo ; enddo

  call data_override('OCN', 'sw', fluxes%sw(:,:), day, &
       is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in) ! scale=US%W_m2_to_QRZ_T
  if (us%QRZ_T_to_W_m2 /= 1.0) then ; do j=js,je ; do i=is,ie
    fluxes%sw(i,j) = fluxes%sw(i,j) * us%W_m2_to_QRZ_T
  enddo ; enddo ; endif

  call data_override('OCN', 'snow', fluxes%fprec(:,:), day, &
       is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in) ! scale=kg_m2_s_conversion

  call data_override('OCN', 'rain', fluxes%lprec(:,:), day, &
       is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in) ! scale=kg_m2_s_conversion

  call data_override('OCN', 'runoff', fluxes%lrunoff(:,:), day, &
       is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in) ! scale=kg_m2_s_conversion

  call data_override('OCN', 'calving', fluxes%frunoff(:,:), day, &
       is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in) ! scale=kg_m2_s_conversion

  if (kg_m2_s_conversion /= 1.0) then ; do j=js,je ; do i=is,ie
    fluxes%lprec(i,j) = fluxes%lprec(i,j) * kg_m2_s_conversion
    fluxes%fprec(i,j) = fluxes%fprec(i,j) * kg_m2_s_conversion
    fluxes%lrunoff(i,j) = fluxes%lrunoff(i,j) * kg_m2_s_conversion
    fluxes%frunoff(i,j) = fluxes%frunoff(i,j) * kg_m2_s_conversion
  enddo ; enddo ; endif

!     Read the SST and SSS fields for damping.
  if (cs%restorebuoy) then !#CTRL# .or. associated(CS%ctrl_forcing_CSp)) then
     call data_override('OCN', 'SST_restore', cs%T_restore(:,:), day, &
          is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)

     call data_override('OCN', 'SSS_restore', cs%S_restore(:,:), day, &
          is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)

  endif

  ! restoring boundary fluxes
  if (cs%restorebuoy) then
    if (cs%use_temperature) then
      do j=js,je ; do i=is,ie
        if (g%mask2dT(i,j) > 0) then
          fluxes%heat_added(i,j) = g%mask2dT(i,j) * &
              ((cs%T_Restore(i,j) - sfc_state%SST(i,j)) * rhoxcp * cs%Flux_const_T)
          fluxes%vprec(i,j) = - (cs%Rho0*cs%Flux_const_S) * &
              (cs%S_Restore(i,j) - sfc_state%SSS(i,j)) / &
              (0.5*(sfc_state%SSS(i,j) + cs%S_Restore(i,j)))
        else
          fluxes%heat_added(i,j) = 0.0
          fluxes%vprec(i,j)      = 0.0
        endif
      enddo ; enddo
    else
      do j=js,je ; do i=is,ie
        if (g%mask2dT(i,j) > 0) then
          fluxes%buoy(i,j) = (cs%Dens_Restore(i,j) - sfc_state%sfc_density(i,j)) * &
                             (cs%G_Earth * cs%Flux_const / cs%Rho0)
        else
          fluxes%buoy(i,j) = 0.0
        endif
      enddo ; enddo
    endif
  else                                              ! not RESTOREBUOY
    if (.not.cs%use_temperature) then
      call mom_error(fatal, "buoyancy_forcing in MOM_surface_forcing: "// &
                     "The fluxes need to be defined without RESTOREBUOY.")
    endif
  endif                                             ! end RESTOREBUOY


  ! mask out land points and compute heat content of water fluxes
  ! assume liquid precip enters ocean at SST
  ! assume frozen precip enters ocean at 0degC
  ! assume liquid runoff enters ocean at SST
  ! assume solid runoff (calving) enters ocean at 0degC
  ! mass leaving ocean has heat_content determined in MOM_diabatic_driver.F90
  do j=js,je ; do i=is,ie
    fluxes%evap(i,j)    = fluxes%evap(i,j)    * g%mask2dT(i,j)
    fluxes%lprec(i,j)   = fluxes%lprec(i,j)   * g%mask2dT(i,j)
    fluxes%fprec(i,j)   = fluxes%fprec(i,j)   * g%mask2dT(i,j)
    fluxes%lrunoff(i,j) = fluxes%lrunoff(i,j) * g%mask2dT(i,j)
    fluxes%frunoff(i,j) = fluxes%frunoff(i,j) * g%mask2dT(i,j)
    fluxes%lw(i,j)      = fluxes%lw(i,j)      * g%mask2dT(i,j)
    fluxes%latent(i,j)  = fluxes%latent(i,j)  * g%mask2dT(i,j)
    fluxes%sens(i,j)    = fluxes%sens(i,j)    * g%mask2dT(i,j)
    fluxes%sw(i,j)      = fluxes%sw(i,j)      * g%mask2dT(i,j)

    fluxes%latent_evap_diag(i,j)     = fluxes%latent_evap_diag(i,j) * g%mask2dT(i,j)
    fluxes%latent_fprec_diag(i,j)    = -fluxes%fprec(i,j)*cs%latent_heat_fusion
    fluxes%latent_frunoff_diag(i,j)  = -fluxes%frunoff(i,j)*cs%latent_heat_fusion
  enddo ; enddo


!#CTRL# if (associated(CS%ctrl_forcing_CSp)) then
!#CTRL#   do j=js,je ; do i=is,ie
!#CTRL#     SST_anom(i,j) = sfc_state%SST(i,j) - CS%T_Restore(i,j)
!#CTRL#     SSS_anom(i,j) = sfc_state%SSS(i,j) - CS%S_Restore(i,j)
!#CTRL#     SSS_mean(i,j) = 0.5*(sfc_state%SSS(i,j) + CS%S_Restore(i,j))
!#CTRL#   enddo ; enddo
!#CTRL#   call apply_ctrl_forcing(SST_anom, SSS_anom, SSS_mean, fluxes%heat_added, &
!#CTRL#                           fluxes%vprec, day, dt, G, CS%ctrl_forcing_CSp)
!#CTRL# endif

  call calltree_leave("buoyancy_forcing_from_data_override")
end subroutine buoyancy_forcing_from_data_override

!> This subroutine specifies zero surface bouyancy fluxes
subroutine buoyancy_forcing_zero(sfc_state, fluxes, day, dt, G, 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 [s]
  type(ocean_grid_type), intent(in)    :: G    !< The ocean's grid structure
  type(surface_forcing_CS), pointer    :: CS   !< pointer to control struct returned by
                                               !! a previous surface_forcing_init call
  ! Local variables
  integer :: i, j, is, ie, js, je

  call calltree_enter("buoyancy_forcing_zero, MOM_surface_forcing.F90")
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec

  if (cs%use_temperature) then
    do j=js,je ; do i=is,ie
      fluxes%evap(i,j)                 = 0.0
      fluxes%lprec(i,j)                = 0.0
      fluxes%fprec(i,j)                = 0.0
      fluxes%vprec(i,j)                = 0.0
      fluxes%lrunoff(i,j)              = 0.0
      fluxes%frunoff(i,j)              = 0.0
      fluxes%lw(i,j)                   = 0.0
      fluxes%latent(i,j)               = 0.0
      fluxes%sens(i,j)                 = 0.0
      fluxes%sw(i,j)                   = 0.0
      fluxes%latent_evap_diag(i,j)     = 0.0
      fluxes%latent_fprec_diag(i,j)    = 0.0
      fluxes%latent_frunoff_diag(i,j)  = 0.0
    enddo ; enddo
  else
    do j=js,je ; do i=is,ie
      fluxes%buoy(i,j) = 0.0
    enddo ; enddo
  endif

  call calltree_leave("buoyancy_forcing_zero")
end subroutine buoyancy_forcing_zero


!> Sets up spatially and temporally constant surface heat fluxes.
subroutine buoyancy_forcing_const(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 [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(surface_forcing_CS), pointer    :: CS   !< pointer to control struct returned by
                                               !! a previous surface_forcing_init call
  ! Local variables
  integer :: i, j, is, ie, js, je
  call calltree_enter("buoyancy_forcing_const, MOM_surface_forcing.F90")
  is  = g%isc ; ie  = g%iec ; js  = g%jsc ; je  = g%jec

  if (cs%use_temperature) then
    do j=js,je ; do i=is,ie
      fluxes%evap(i,j)                 = 0.0
      fluxes%lprec(i,j)                = 0.0
      fluxes%fprec(i,j)                = 0.0
      fluxes%vprec(i,j)                = 0.0
      fluxes%lrunoff(i,j)              = 0.0
      fluxes%frunoff(i,j)              = 0.0
      fluxes%lw(i,j)                   = 0.0
      fluxes%latent(i,j)               = 0.0
      fluxes%sens(i,j)                 = cs%constantHeatForcing * g%mask2dT(i,j)
      fluxes%sw(i,j)                   = 0.0
      fluxes%latent_evap_diag(i,j)     = 0.0
      fluxes%latent_fprec_diag(i,j)    = 0.0
      fluxes%latent_frunoff_diag(i,j)  = 0.0
    enddo ; enddo
  else
    do j=js,je ; do i=is,ie
      fluxes%buoy(i,j) = 0.0
    enddo ; enddo
  endif

  call calltree_leave("buoyancy_forcing_const")
end subroutine buoyancy_forcing_const

!> Sets surface fluxes of heat and salinity by restoring to temperature and
!! salinity profiles that vary linearly with latitude.
subroutine buoyancy_forcing_linear(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 [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(surface_forcing_CS), pointer    :: CS   !< pointer to control struct returned by
                                               !! a previous surface_forcing_init call
  ! Local variables
  real :: y, T_restore, S_restore
  integer :: i, j, is, ie, js, je

  call calltree_enter("buoyancy_forcing_linear, MOM_surface_forcing.F90")
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec

  ! This case has no surface buoyancy forcing.
  if (cs%use_temperature) then
    do j=js,je ; do i=is,ie
      fluxes%evap(i,j)                 = 0.0
      fluxes%lprec(i,j)                = 0.0
      fluxes%fprec(i,j)                = 0.0
      fluxes%vprec(i,j)                = 0.0
      fluxes%lrunoff(i,j)              = 0.0
      fluxes%frunoff(i,j)              = 0.0
      fluxes%lw(i,j)                   = 0.0
      fluxes%latent(i,j)               = 0.0
      fluxes%sens(i,j)                 = 0.0
      fluxes%sw(i,j)                   = 0.0
      fluxes%latent_evap_diag(i,j)     = 0.0
      fluxes%latent_fprec_diag(i,j)    = 0.0
      fluxes%latent_frunoff_diag(i,j)  = 0.0
    enddo ; enddo
  else
    do j=js,je ; do i=is,ie
      fluxes%buoy(i,j) = 0.0
    enddo ; enddo
  endif

  if (cs%restorebuoy) then
    if (cs%use_temperature) then
      do j=js,je ; do i=is,ie
        y = (g%geoLatCu(i,j)-cs%South_lat)/cs%len_lat
        t_restore = cs%T_south + (cs%T_north-cs%T_south)*y
        s_restore = cs%S_south + (cs%S_north-cs%S_south)*y
        if (g%mask2dT(i,j) > 0) then
          fluxes%heat_added(i,j) = g%mask2dT(i,j) * &
              ((t_restore - sfc_state%SST(i,j)) * ((cs%Rho0 * fluxes%C_p) * cs%Flux_const))
          fluxes%vprec(i,j) = - (cs%Rho0*cs%Flux_const) * &
              (s_restore - sfc_state%SSS(i,j)) / &
              (0.5*(sfc_state%SSS(i,j) + s_restore))
        else
          fluxes%heat_added(i,j) = 0.0
          fluxes%vprec(i,j)      = 0.0
        endif
      enddo ; enddo
    else
      call mom_error(fatal, "buoyancy_forcing_linear in MOM_surface_forcing: "// &
                     "RESTOREBUOY to linear not written yet.")
     !do j=js,je ; do i=is,ie
     !  if (G%mask2dT(i,j) > 0) then
     !    fluxes%buoy(i,j) = (CS%Dens_Restore(i,j) - sfc_state%sfc_density(i,j)) * &
     !                       (CS%G_Earth * CS%Flux_const / CS%Rho0)
     !  else
     !    fluxes%buoy(i,j) = 0.0
     !  endif
     !enddo ; enddo
    endif
  else                                              ! not RESTOREBUOY
    if (.not.cs%use_temperature) then
      call mom_error(fatal, "buoyancy_forcing_linear in MOM_surface_forcing: "// &
                     "The fluxes need to be defined without RESTOREBUOY.")
    endif
  endif                                             ! end RESTOREBUOY

  call calltree_leave("buoyancy_forcing_linear")
end subroutine buoyancy_forcing_linear

!> Save a restart file for the forcing fields
subroutine forcing_save_restart(CS, G, Time, directory, time_stamped, &
                                filename_suffix)
  type(surface_forcing_cs),   pointer       :: cs   !< pointer to control struct returned by
                                                    !! a previous surface_forcing_init call
  type(ocean_grid_type),      intent(inout) :: g    !< The ocean's grid structure
  type(time_type),            intent(in)    :: time !< model time at this call; needed for mpp_write calls
  character(len=*),           intent(in)    :: directory !< directory into which to write these restart files
  logical,          optional, intent(in)    :: time_stamped !< If true, the restart file names
                                                    !! include a unique time stamp; the  default is false.
  character(len=*), optional, intent(in)    :: filename_suffix !< optional suffix (e.g., a time-stamp)
                                                    !! to append to the restart fname

  if (.not.associated(cs)) return
  if (.not.associated(cs%restart_CSp)) return

  call save_restart(directory, time, g, cs%restart_CSp, time_stamped)

end subroutine forcing_save_restart

!> Initialize the surface forcing module
subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, tracer_flow_CSp)
  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(inout) :: diag !< structure used to regulate diagnostic output
  type(surface_forcing_cs),     pointer       :: cs   !< pointer to control struct returned by
                                                      !! a previous surface_forcing_init call
  type(tracer_flow_control_cs), pointer       :: tracer_flow_csp !< Forcing for tracers?

  ! Local variables
  type(directories)  :: dirs
  logical            :: new_sim
  type(time_type)    :: time_frc
  ! This include declares and sets the variable "version".
# include "version_variable.h"
  real :: flux_const_default ! The unscaled value of FLUXCONST [m day-1]
  logical :: default_2018_answers
  character(len=40)  :: mdl = "MOM_surface_forcing" ! This module's name.
  character(len=200) :: filename, gust_file ! The name of the gustiness input file.

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

  id_clock_forcing=cpu_clock_id('(Ocean surface forcing)', grain=clock_module)
  call cpu_clock_begin(id_clock_forcing)

  cs%diag => diag
  if (associated(tracer_flow_csp)) cs%tracer_flow_CSp => tracer_flow_csp

  ! 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, "INPUTDIR", cs%inputdir, &
                 "The directory in which all input files are found.", &
                 default=".")
  cs%inputdir = slasher(cs%inputdir)

  call get_param(param_file, mdl, "ADIABATIC", cs%adiabatic, &
                 "There are no diapycnal mass fluxes if ADIABATIC is "//&
                 "true. This assumes that KD = KDML = 0.0 and that "//&
                 "there is no buoyancy forcing, but makes the model "//&
                 "faster by eliminating subroutine calls.", default=.false.)
  call get_param(param_file, mdl, "VARIABLE_WINDS", cs%variable_winds, &
                 "If true, the winds vary in time after the initialization.", &
                 default=.true.)
  call get_param(param_file, mdl, "VARIABLE_BUOYFORCE", cs%variable_buoyforce, &
                 "If true, the buoyancy forcing varies in time after the "//&
                 "initialization of the model.", default=.true.)

  call get_param(param_file, mdl, "BUOY_CONFIG", cs%buoy_config, &
                 "The character string that indicates how buoyancy forcing "//&
                 "is specified. Valid options include (file), (zero), "//&
                 "(linear), (USER), (BFB) and (NONE).", default="zero")
  if (trim(cs%buoy_config) == "file") then
    call get_param(param_file, mdl, "ARCHAIC_OMIP_FORCING_FILE", cs%archaic_OMIP_file, &
                 "If true, use the forcing variable decomposition from "//&
                 "the old German OMIP prescription that predated CORE. If "//&
                 "false, use the variable groupings available from MOM "//&
                 "output diagnostics of forcing variables.", default=.true.)
    if (cs%archaic_OMIP_file) then
      call get_param(param_file, mdl, "LONGWAVEDOWN_FILE", cs%longwave_file, &
                 "The file with the downward longwave heat flux, in "//&
                 "variable lwdn_sfc.", fail_if_missing=.true.)
      call get_param(param_file, mdl, "LONGWAVEUP_FILE", cs%longwaveup_file, &
                 "The file with the upward longwave heat flux, in "//&
                 "variable lwup_sfc.", fail_if_missing=.true.)
      call get_param(param_file, mdl, "EVAPORATION_FILE", cs%evaporation_file, &
                 "The file with the evaporative moisture flux, in "//&
                 "variable evap.", fail_if_missing=.true.)
      call get_param(param_file, mdl, "SENSIBLEHEAT_FILE", cs%sensibleheat_file, &
                 "The file with the sensible heat flux, in "//&
                 "variable shflx.", fail_if_missing=.true.)
      call get_param(param_file, mdl, "SHORTWAVEUP_FILE", cs%shortwaveup_file, &
                 "The file with the upward shortwave heat flux.", &
                 fail_if_missing=.true.)
      call get_param(param_file, mdl, "SHORTWAVEDOWN_FILE", cs%shortwave_file, &
                 "The file with the downward shortwave heat flux.", &
                 fail_if_missing=.true.)
      call get_param(param_file, mdl, "SNOW_FILE", cs%snow_file, &
                 "The file with the downward frozen precip flux, in "//&
                 "variable snow.", fail_if_missing=.true.)
      call get_param(param_file, mdl, "PRECIP_FILE", cs%rain_file, &
                 "The file with the downward total precip flux, in "//&
                 "variable precip.", fail_if_missing=.true.)
      call get_param(param_file, mdl, "FRESHDISCHARGE_FILE", cs%runoff_file, &
                 "The file with the fresh and frozen runoff/calving fluxes, "//&
                 "invariables disch_w and disch_s.", fail_if_missing=.true.)

      ! These variable names are hard-coded, per the archaic OMIP conventions.
      cs%latentheat_file = cs%evaporation_file ; cs%latent_var = "evap"
      cs%LW_var = "lwdn_sfc"; cs%SW_var = "swdn_sfc"; cs%sens_var = "shflx"
      cs%evap_var = "evap"; cs%rain_var = "precip"; cs%snow_var = "snow"
      cs%lrunoff_var = "disch_w"; cs%frunoff_var = "disch_s"

    else
      call get_param(param_file, mdl, "LONGWAVE_FILE", cs%longwave_file, &
                 "The file with the longwave heat flux, in the variable "//&
                 "given by LONGWAVE_FORCING_VAR.", fail_if_missing=.true.)
      call get_param(param_file, mdl, "LONGWAVE_FORCING_VAR", cs%LW_var, &
                 "The variable with the longwave forcing field.", default="LW")

      call get_param(param_file, mdl, "SHORTWAVE_FILE", cs%shortwave_file, &
                 "The file with the shortwave heat flux, in the variable "//&
                 "given by SHORTWAVE_FORCING_VAR.", fail_if_missing=.true.)
      call get_param(param_file, mdl, "SHORTWAVE_FORCING_VAR", cs%SW_var, &
                 "The variable with the shortwave forcing field.", default="SW")

      call get_param(param_file, mdl, "EVAPORATION_FILE", cs%evaporation_file, &
                 "The file with the evaporative moisture flux, in the "//&
                 "variable given by EVAP_FORCING_VAR.", fail_if_missing=.true.)
      call get_param(param_file, mdl, "EVAP_FORCING_VAR", cs%evap_var, &
                 "The variable with the evaporative moisture flux.", &
                 default="evap")

      call get_param(param_file, mdl, "LATENTHEAT_FILE", cs%latentheat_file, &
                 "The file with the latent heat flux, in the variable "//&
                 "given by LATENT_FORCING_VAR.", fail_if_missing=.true.)
      call get_param(param_file, mdl, "LATENT_FORCING_VAR", cs%latent_var, &
                 "The variable with the latent heat flux.", default="latent")

      call get_param(param_file, mdl, "SENSIBLEHEAT_FILE", cs%sensibleheat_file, &
                 "The file with the sensible heat flux, in the variable "//&
                 "given by SENSIBLE_FORCING_VAR.", fail_if_missing=.true.)
      call get_param(param_file, mdl, "SENSIBLE_FORCING_VAR", cs%sens_var, &
                 "The variable with the sensible heat flux.", default="sensible")

      call get_param(param_file, mdl, "RAIN_FILE", cs%rain_file, &
                 "The file with the liquid precipitation flux, in the "//&
                 "variable given by RAIN_FORCING_VAR.", fail_if_missing=.true.)
      call get_param(param_file, mdl, "RAIN_FORCING_VAR", cs%rain_var, &
                 "The variable with the liquid precipitation flux.", &
                 default="liq_precip")
      call get_param(param_file, mdl, "SNOW_FILE", cs%snow_file, &
                 "The file with the frozen precipitation flux, in the "//&
                 "variable given by SNOW_FORCING_VAR.", fail_if_missing=.true.)
      call get_param(param_file, mdl, "SNOW_FORCING_VAR", cs%snow_var, &
                 "The variable with the frozen precipitation flux.", &
                 default="froz_precip")

      call get_param(param_file, mdl, "RUNOFF_FILE", cs%runoff_file, &
                 "The file with the fresh and frozen runoff/calving "//&
                 "fluxes, in variables given by LIQ_RUNOFF_FORCING_VAR "//&
                 "and FROZ_RUNOFF_FORCING_VAR.", fail_if_missing=.true.)
      call get_param(param_file, mdl, "LIQ_RUNOFF_FORCING_VAR", cs%lrunoff_var, &
                 "The variable with the liquid runoff flux.", &
                 default="liq_runoff")
      call get_param(param_file, mdl, "FROZ_RUNOFF_FORCING_VAR", cs%frunoff_var, &
                 "The variable with the frozen runoff flux.", &
                 default="froz_runoff")
    endif

    call get_param(param_file, mdl, "SSTRESTORE_FILE", cs%SSTrestore_file, &
                 "The file with the SST toward which to restore in the "//&
                 "variable given by SST_RESTORE_VAR.", fail_if_missing=.true.)
    call get_param(param_file, mdl, "SALINITYRESTORE_FILE", cs%salinityrestore_file, &
                 "The file with the surface salinity toward which to "//&
                 "restore in the variable given by SSS_RESTORE_VAR.", &
                 fail_if_missing=.true.)

    if (cs%archaic_OMIP_file) then
      cs%SST_restore_var = "TEMP" ; cs%SSS_restore_var = "SALT"
    else
      call get_param(param_file, mdl, "SST_RESTORE_VAR", cs%SST_restore_var, &
                 "The variable with the SST toward which to restore.", &
                 default="SST")
      call get_param(param_file, mdl, "SSS_RESTORE_VAR", cs%SSS_restore_var, &
                 "The variable with the SSS toward which to restore.", &
                 default="SSS")
    endif

    ! Add inputdir to the file names.
    cs%shortwave_file = trim(cs%inputdir)//trim(cs%shortwave_file)
    cs%longwave_file = trim(cs%inputdir)//trim(cs%longwave_file)
    cs%sensibleheat_file = trim(cs%inputdir)//trim(cs%sensibleheat_file)
    cs%latentheat_file = trim(cs%inputdir)//trim(cs%latentheat_file)
    cs%evaporation_file = trim(cs%inputdir)//trim(cs%evaporation_file)
    cs%snow_file = trim(cs%inputdir)//trim(cs%snow_file)
    cs%rain_file = trim(cs%inputdir)//trim(cs%rain_file)
    cs%runoff_file = trim(cs%inputdir)//trim(cs%runoff_file)

    cs%shortwaveup_file = trim(cs%inputdir)//trim(cs%shortwaveup_file)
    cs%longwaveup_file = trim(cs%inputdir)//trim(cs%longwaveup_file)

    cs%SSTrestore_file = trim(cs%inputdir)//trim(cs%SSTrestore_file)
    cs%salinityrestore_file = trim(cs%inputdir)//trim(cs%salinityrestore_file)
  elseif (trim(cs%buoy_config) == "const") then
    call get_param(param_file, mdl, "SENSIBLE_HEAT_FLUX", cs%constantHeatForcing, &
                 "A constant heat forcing (positive into ocean) applied "//&
                 "through the sensible heat flux field. ", &
                 units='W/m2', scale=us%W_m2_to_QRZ_T, fail_if_missing=.true.)
  endif
  call get_param(param_file, mdl, "WIND_CONFIG", cs%wind_config, &
                 "The character string that indicates how wind forcing "//&
                 "is specified. Valid options include (file), (2gyre), "//&
                 "(1gyre), (gyres), (zero), and (USER).", default="zero")
  if (trim(cs%wind_config) == "file") then
    call get_param(param_file, mdl, "WIND_FILE", cs%wind_file, &
                 "The file in which the wind stresses are found in "//&
                 "variables STRESS_X and STRESS_Y.", fail_if_missing=.true.)
    call get_param(param_file, mdl, "WINDSTRESS_X_VAR",cs%stress_x_var, &
                 "The name of the x-wind stress variable in WIND_FILE.", &
                 default="STRESS_X")
    call get_param(param_file, mdl, "WINDSTRESS_Y_VAR", cs%stress_y_var, &
                 "The name of the y-wind stress variable in WIND_FILE.", &
                 default="STRESS_Y")
    call get_param(param_file, mdl, "WIND_STAGGER",cs%wind_stagger, &
                 "A character indicating how the wind stress components "//&
                 "are staggered in WIND_FILE.  This may be A or C for now.", &
                 default="C")
    call get_param(param_file, mdl, "WINDSTRESS_SCALE", cs%wind_scale, &
                 "A value by which the wind stresses in WIND_FILE are rescaled.", &
                 default=1.0, units="nondim")
    call get_param(param_file, mdl, "USTAR_FORCING_VAR", cs%ustar_var, &
                 "The name of the friction velocity variable in WIND_FILE "//&
                 "or blank to get ustar from the wind stresses plus the "//&
                 "gustiness.", default=" ", units="nondim")
    cs%wind_file = trim(cs%inputdir) // trim(cs%wind_file)
  endif
  if (trim(cs%wind_config) == "gyres") then
    call get_param(param_file, mdl, "TAUX_CONST", cs%gyres_taux_const, &
                 "With the gyres wind_config, the constant offset in the "//&
                 "zonal wind stress profile: "//&
                 "  A in taux = A + B*sin(n*pi*y/L) + C*cos(n*pi*y/L).", &
                 units="Pa", default=0.0)
    call get_param(param_file, mdl, "TAUX_SIN_AMP",cs%gyres_taux_sin_amp, &
                 "With the gyres wind_config, the sine amplitude in the "//&
                 "zonal wind stress profile: "//&
                 "  B in taux = A + B*sin(n*pi*y/L) + C*cos(n*pi*y/L).", &
                 units="Pa", default=0.0)
    call get_param(param_file, mdl, "TAUX_COS_AMP",cs%gyres_taux_cos_amp, &
                 "With the gyres wind_config, the cosine amplitude in "//&
                 "the zonal wind stress profile: "//&
                 "  C in taux = A + B*sin(n*pi*y/L) + C*cos(n*pi*y/L).", &
                 units="Pa", default=0.0)
    call get_param(param_file, mdl, "TAUX_N_PIS",cs%gyres_taux_n_pis, &
                 "With the gyres wind_config, the number of gyres in "//&
                 "the zonal wind stress profile: "//&
                 "  n in taux = A + B*sin(n*pi*y/L) + C*cos(n*pi*y/L).", &
                 units="nondim", default=0.0)
    call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
                 "This sets the default value for the various _2018_ANSWERS parameters.", &
                 default=.false.)
    call get_param(param_file, mdl, "WIND_GYRES_2018_ANSWERS", cs%answers_2018, &
                 "If true, use the order of arithmetic and expressions that recover the answers "//&
                 "from the end of 2018.  Otherwise, use expressions for the gyre friction velocities "//&
                 "that are rotationally invariant and more likely to be the same between compilers.", &
                 default=default_2018_answers)
  else
    cs%answers_2018 = .false.
  endif
  if (trim(cs%wind_config) == "scurves") then
    call get_param(param_file, mdl, "WIND_SCURVES_LATS", cs%scurves_ydata, &
                 "A list of latitudes defining a piecewise scurve profile "//&
                 "for zonal wind stress.", &
                 units="degrees N", fail_if_missing=.true.)
    call get_param(param_file, mdl, "WIND_SCURVES_TAUX", cs%scurves_taux, &
                 "A list of zonal wind stress values at latitudes "//&
                 "WIND_SCURVES_LATS defining a piecewise scurve profile.", &
                 units="Pa", fail_if_missing=.true.)
  endif
  if ((trim(cs%wind_config) == "2gyre") .or. &
      (trim(cs%wind_config) == "1gyre") .or. &
      (trim(cs%wind_config) == "gyres") .or. &
      (trim(cs%buoy_config) == "linear")) then
    cs%south_lat = g%south_lat
    cs%len_lat = g%len_lat
  endif

  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 BOUSSINSEQ 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, "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.)
  call get_param(param_file, mdl, "LATENT_HEAT_FUSION", cs%latent_heat_fusion, &
                 "The latent heat of fusion.", default=hlf, &
                 units="J/kg", scale=us%J_kg_to_Q)
  call get_param(param_file, mdl, "LATENT_HEAT_VAPORIZATION", cs%latent_heat_vapor, &
                 "The latent heat of fusion.", default=hlv, units="J/kg", scale=us%J_kg_to_Q)
  if (cs%restorebuoy) then
    ! These three variables use non-standard time units, but are rescaled as they are read.
    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*us%T_to_s/86400.0, &
                 unscaled=flux_const_default)

    if (cs%use_temperature) then
      call get_param(param_file, mdl, "FLUXCONST_T", cs%Flux_const_T, &
           "The constant that relates the restoring surface temperature "//&
           "flux to the relative surface anomaly (akin to a piston "//&
           "velocity).  Note the non-MKS units.", &
           units="m day-1", scale=us%m_to_Z*us%T_to_s/86400.0, &
           default=flux_const_default)
      call get_param(param_file, mdl, "FLUXCONST_S", cs%Flux_const_S, &
           "The constant that relates the restoring surface salinity "//&
           "flux to the relative surface anomaly (akin to a piston "//&
           "velocity).  Note the non-MKS units.", &
           units="m day-1", scale=us%m_to_Z*us%T_to_s/86400.0, &
           default=flux_const_default)
    endif

    if (trim(cs%buoy_config) == "linear") then
      call get_param(param_file, mdl, "SST_NORTH", cs%T_north, &
                 "With buoy_config linear, the sea surface temperature "//&
                 "at the northern end of the domain toward which to "//&
                 "to restore.", units="deg C", default=0.0)
      call get_param(param_file, mdl, "SST_SOUTH", cs%T_south, &
                 "With buoy_config linear, the sea surface temperature "//&
                 "at the southern end of the domain toward which to "//&
                 "to restore.", units="deg C", default=0.0)
      call get_param(param_file, mdl, "SSS_NORTH", cs%S_north, &
                 "With buoy_config linear, the sea surface salinity "//&
                 "at the northern end of the domain toward which to "//&
                 "to restore.", units="PSU", default=35.0)
      call get_param(param_file, mdl, "SSS_SOUTH", cs%S_south, &
                 "With buoy_config linear, the sea surface salinity "//&
                 "at the southern end of the domain toward which to "//&
                 "to restore.", units="PSU", default=35.0)
    endif
  endif
  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, "GUST_CONST", cs%gust_const, &
                 "The background gustiness in the winds.", &
                 units="Pa", default=0.0, scale=us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z)
  call get_param(param_file, mdl, "FIX_USTAR_GUSTLESS_BUG", cs%fix_ustar_gustless_bug, &
                 "If true correct a bug in the time-averaging of the gustless wind friction velocity", &
                 default=.true.)
  call get_param(param_file, mdl, "READ_GUST_2D", cs%read_gust_2d, &
                 "If true, use a 2-dimensional gustiness supplied from "//&
                 "an input file", default=.false.)
  if (cs%read_gust_2d) then
    call get_param(param_file, mdl, "GUST_2D_FILE", gust_file, &
                 "The file in which the wind gustiness is found in "//&
                 "variable gustiness.", fail_if_missing=.true.)
    call safe_alloc_ptr(cs%gust,g%isd,g%ied,g%jsd,g%jed)
    filename = trim(cs%inputdir) // trim(gust_file)
    call mom_read_data(filename,'gustiness',cs%gust,g%domain, timelevel=1, &
                   scale=us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z) ! units in file should be Pa
  endif

!  All parameter settings are now known.

  if (trim(cs%wind_config) == "USER" .or. trim(cs%buoy_config) == "USER" ) then
    call user_surface_forcing_init(time, g, us, param_file, diag, cs%user_forcing_CSp)
  elseif (trim(cs%buoy_config) == "BFB" ) then
    call bfb_surface_forcing_init(time, g, us, param_file, diag, cs%BFB_forcing_CSp)
  elseif (trim(cs%buoy_config) == "dumbbell" ) then
    call dumbbell_surface_forcing_init(time, g, us, param_file, diag, cs%dumbbell_forcing_CSp)
  elseif (trim(cs%wind_config) == "MESO" .or. trim(cs%buoy_config) == "MESO" ) then
    call meso_surface_forcing_init(time, g, us, param_file, diag, cs%MESO_forcing_CSp)
  elseif (trim(cs%wind_config) == "ideal_hurr" .or.&
          trim(cs%wind_config) == "SCM_ideal_hurr") then
    call idealized_hurricane_wind_init(time, g, us, param_file, cs%idealized_hurricane_CSp)
  elseif (trim(cs%wind_config) == "const") then
    call get_param(param_file, mdl, "CONST_WIND_TAUX", cs%tau_x0, &
                 "With wind_config const, this is the constant zonal "//&
                 "wind-stress", units="Pa", fail_if_missing=.true.)
    call get_param(param_file, mdl, "CONST_WIND_TAUY", cs%tau_y0, &
                 "With wind_config const, this is the constant meridional "//&
                 "wind-stress", units="Pa", fail_if_missing=.true.)
  elseif (trim(cs%wind_config) == "SCM_CVmix_tests" .or. &
          trim(cs%buoy_config) == "SCM_CVmix_tests") then
    call scm_cvmix_tests_surface_forcing_init(time, g, param_file, cs%SCM_CVmix_tests_CSp)
  endif

  call register_forcing_type_diags(time, diag, us, cs%use_temperature, cs%handles)

  ! Set up any restart fields associated with the forcing.
  call restart_init(param_file, cs%restart_CSp, "MOM_forcing.res")
!#CTRL#  call register_ctrl_forcing_restarts(G, param_file, CS%ctrl_forcing_CSp, &
!#CTRL#                                      CS%restart_CSp)
  call restart_init_end(cs%restart_CSp)

  if (associated(cs%restart_CSp)) then
    call get_mom_input(dirs=dirs)

    new_sim = .false.
    if ((dirs%input_filename(1:1) == 'n') .and. &
        (len_trim(dirs%input_filename) == 1)) new_sim = .true.
    if (.not.new_sim) then
      call restore_state(dirs%input_filename, dirs%restart_input_dir, time_frc, &
                         g, cs%restart_CSp)
    endif
  endif

  ! Determine how many time levels are in each forcing variable.
  if (trim(cs%buoy_config) == "file") then
    cs%SW_nlev = num_timelevels(cs%shortwave_file, cs%SW_var, min_dims=3)
    cs%LW_nlev = num_timelevels(cs%longwave_file, cs%LW_var, min_dims=3)
    cs%latent_nlev = num_timelevels(cs%latentheat_file, cs%latent_var, 3)
    cs%sens_nlev = num_timelevels(cs%sensibleheat_file, cs%sens_var, min_dims=3)

    cs%evap_nlev = num_timelevels(cs%evaporation_file, cs%evap_var, min_dims=3)
    cs%precip_nlev = num_timelevels(cs%rain_file, cs%rain_var, min_dims=3)
    cs%runoff_nlev = num_timelevels(cs%runoff_file, cs%lrunoff_var, 3)

    cs%SST_nlev = num_timelevels(cs%SSTrestore_file, cs%SST_restore_var, 3)
    cs%SSS_nlev = num_timelevels(cs%salinityrestore_file, cs%SSS_restore_var, 3)
  endif

  if (trim(cs%wind_config) == "file") &
    cs%wind_nlev = num_timelevels(cs%wind_file, cs%stress_x_var, min_dims=3)

!#CTRL#  call controlled_forcing_init(Time, G, param_file, diag, CS%ctrl_forcing_CSp)

  call user_revise_forcing_init(param_file, cs%urf_CS)

  call cpu_clock_end(id_clock_forcing)
end subroutine surface_forcing_init


!> Deallocate memory associated with the surface forcing module
subroutine surface_forcing_end(CS, fluxes)
  type(surface_forcing_CS), pointer    :: CS   !< pointer to control struct returned by
                                               !! a previous surface_forcing_init call
  type(forcing), optional,  intent(inout) :: fluxes !< A structure containing thermodynamic forcing fields
! Arguments:  CS - A pointer to the control structure returned by a previous
!                  call to surface_forcing_init, it will be deallocated here.
!  (inout)    fluxes - A structure containing pointers to any possible
!                     forcing fields.  Unused fields have NULL ptrs.

  if (present(fluxes)) call deallocate_forcing_type(fluxes)

!#CTRL#  call controlled_forcing_end(CS%ctrl_forcing_CSp)

  if (associated(cs)) deallocate(cs)
  cs => null()

end subroutine surface_forcing_end

end module mom_surface_forcing