MOM_controlled_forcing.F90

!> Use control-theory to adjust the surface heat flux and precipitation.
!!
!! Adjustments are based on the time-mean or periodically (seasonally) varying
!! anomalies from the observed state.
!!
!! The techniques behind this are described in Hallberg and Adcroft (2018, in prep.).
module mom_controlled_forcing

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

use mom_diag_mediator, only : post_data, query_averaging_enabled
use mom_diag_mediator, only : register_diag_field, diag_ctrl, safe_alloc_ptr
use mom_domains, only : pass_var, pass_vector, agrid, to_south, to_west, to_all
use mom_error_handler, only : mom_error, fatal, warning, mom_mesg, is_root_pe
use mom_file_parser, only : read_param, get_param, log_param, log_version, param_file_type
use mom_forcing_type, only : forcing
use mom_grid, only : ocean_grid_type
use mom_io, only : vardesc, var_desc
use mom_restart, only : register_restart_field, mom_restart_cs
use mom_time_manager, only : time_type, operator(+), operator(/), operator(-)
use mom_time_manager, only : get_date, set_date
use mom_time_manager, only : time_type_to_real, real_to_time
use mom_unit_scaling, only : unit_scale_type
use mom_variables, only : surface

implicit none ; private

#include <MOM_memory.h>

public apply_ctrl_forcing, register_ctrl_forcing_restarts
public controlled_forcing_init, controlled_forcing_end

!> Control structure for MOM_controlled_forcing
type, public :: ctrl_forcing_cs ; private
  logical :: use_temperature !< If true, temperature and salinity are used as
                             !! state variables.
  logical :: do_integrated   !< If true, use time-integrated anomalies to control
                             !! the surface state.
  integer :: num_cycle       !< The number of elements in the forcing cycle.
  real    :: heat_int_rate  !< The rate at which heating anomalies accumulate [s-1].
  real    :: prec_int_rate  !< The rate at which precipitation anomalies accumulate [s-1].
  real    :: heat_cyc_rate  !< The rate at which cyclical heating anomaliess
                            !! accumulate [s-1].
  real    :: prec_cyc_rate  !< The rate at which cyclical precipitation anomaliess
                            !! accumulate [s-1].
  real    :: len2           !< The square of the length scale over which the anomalies
                            !! are smoothed via a Laplacian filter [m2].
  real    :: lam_heat       !< A constant of proportionality between SST anomalies
                            !! and heat fluxes [W m-2 degC-1].
  real    :: lam_prec       !< A constant of proportionality between SSS anomalies
                            !! (normalised by mean SSS) and precipitation [kg m-2].
  real    :: lam_cyc_heat   !< A constant of proportionality between cyclical SST
                            !! anomalies and corrective heat fluxes [W m-2 degC-1].
  real    :: lam_cyc_prec   !< A constant of proportionality between cyclical SSS
                            !! anomalies (normalised by mean SSS) and corrective
                            !! precipitation [kg m-2].

  !>@{ Pointers for data.
  !! \todo Needs more complete documentation.
  real, pointer, dimension(:) :: &
    avg_time => null()
  real, pointer, dimension(:,:) :: &
    heat_0 => null(), &
    precip_0 => null()
  real, pointer, dimension(:,:,:) :: &
    heat_cyc => null(), &
    precip_cyc => null(), &
    avg_sst_anom => null(), &
    avg_sss_anom => null(), &
    avg_sss => null()
  !>@}
  type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
                            !! regulate the timing of diagnostic output.
  integer :: id_heat_0 = -1 !< Diagnostic handle
end type ctrl_forcing_cs

contains

!> This subroutine calls any of the other subroutines in this file
!! that are needed to specify the current surface forcing fields.
subroutine apply_ctrl_forcing(SST_anom, SSS_anom, SSS_mean, virt_heat, virt_precip, &
                              day_start, dt, G, US, CS)
  type(ocean_grid_type), intent(inout) :: g                    !< The ocean's grid structure.
  real, dimension(SZI_(G),SZJ_(G)), intent(in)    :: sst_anom  !< The sea surface temperature
                                                               !! anomalies [degC].
  real, dimension(SZI_(G),SZJ_(G)), intent(in)    :: sss_anom  !< The sea surface salinity
                                                               !! anomlies [ppt].
  real, dimension(SZI_(G),SZJ_(G)), intent(in)    :: sss_mean  !< The mean sea surface
                                                               !! salinity [ppt].
  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: virt_heat !< Virtual (corrective) heat
                                                               !! fluxes that are augmented
                                                               !! in this subroutine [W m-2].
  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: virt_precip !< Virtual (corrective)
                                                               !! precipitation fluxes that
                                                               !! are augmented in this
                                                               !! subroutine [kg m-2 s-1].
  type(time_type),       intent(in)    :: day_start      !< Start time of the fluxes.
  real,                  intent(in)    :: dt             !< Length of time over which these
                                                         !! fluxes will be applied [s].
  type(unit_scale_type), intent(in)    :: us             !< A dimensional unit scaling type
  type(ctrl_forcing_cs), pointer       :: cs             !< A pointer to the control structure
                                                         !! returned by a previous call to
                                                         !! ctrl_forcing_init.
!
  real, dimension(SZIB_(G),SZJ_(G)) :: &
    flux_heat_x, &
    flux_prec_x
  real, dimension(SZI_(G),SZJB_(G)) :: &
    flux_heat_y, &
    flux_prec_y
  type(time_type) :: day_end
  real    :: coef    ! A heat-flux coefficient [m2].
  real    :: mr_st, mr_end, mr_mid, mr_prev, mr_next
  real    :: dt_wt, dt_heat_rate, dt_prec_rate
  real    :: dt1_heat_rate, dt1_prec_rate, dt2_heat_rate, dt2_prec_rate
  real    :: wt_per1, wt_st, wt_end, wt_mid
  integer :: m_st, m_end, m_mid, m_u1, m_u2, m_u3
  integer :: yr, mon, day, hr, min, sec
  integer :: i, j, is, ie, js, je

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

  if (.not.associated(cs)) return
  if ((cs%num_cycle <= 0) .and. (.not.cs%do_integrated)) return

  day_end = day_start + real_to_time(dt)

  do j=js,je ; do i=is,ie
    virt_heat(i,j) = 0.0 ; virt_precip(i,j) = 0.0
  enddo ; enddo

  if (cs%do_integrated) then
    dt_heat_rate = dt * cs%heat_int_rate
    dt_prec_rate = dt * cs%prec_int_rate
    call pass_var(cs%heat_0, g%Domain, complete=.false.)
    call pass_var(cs%precip_0, g%Domain)

    do j=js,je ; do i=is-1,ie
      coef = cs%Len2 * (g%dy_Cu(i,j)*g%IdxCu(i,j))
      flux_heat_x(i,j) = coef * (cs%heat_0(i,j) - cs%heat_0(i+1,j))
      flux_prec_x(i,j) = coef * (cs%precip_0(i,j) - cs%precip_0(i+1,j))
    enddo ; enddo
    do j=js-1,je ; do i=is,ie
      coef = cs%Len2 * (g%dx_Cv(i,j)*g%IdyCv(i,j))
      flux_heat_y(i,j) = coef * (cs%heat_0(i,j) - cs%heat_0(i,j+1))
      flux_prec_y(i,j) = coef * (cs%precip_0(i,j) - cs%precip_0(i,j+1))
    enddo ; enddo
    do j=js,je ; do i=is,ie
      cs%heat_0(i,j) = cs%heat_0(i,j) + dt_heat_rate * ( &
         -cs%lam_heat*g%mask2dT(i,j)*sst_anom(i,j) + &
        (us%m_to_L**2*g%IareaT(i,j) * ((flux_heat_x(i-1,j) - flux_heat_x(i,j)) + &
                          (flux_heat_y(i,j-1) - flux_heat_y(i,j))) ) )

      cs%precip_0(i,j) = cs%precip_0(i,j) + dt_prec_rate * ( &
         cs%lam_prec * g%mask2dT(i,j)*(sss_anom(i,j) / sss_mean(i,j)) + &
        (us%m_to_L**2*g%IareaT(i,j) * ((flux_prec_x(i-1,j) - flux_prec_x(i,j)) + &
                          (flux_prec_y(i,j-1) - flux_prec_y(i,j))) ) )

      virt_heat(i,j) = virt_heat(i,j) + cs%heat_0(i,j)
      virt_precip(i,j) = virt_precip(i,j) + cs%precip_0(i,j)
    enddo ; enddo
  endif

  if (cs%num_cycle > 0) then
    ! Determine the current period, with values that run from 0 to CS%num_cycle.
    call get_date(day_start, yr, mon, day, hr, min, sec)
    mr_st = cs%num_cycle * (time_type_to_real(day_start - set_date(yr, 1, 1)) / &
                   time_type_to_real(set_date(yr+1, 1, 1) - set_date(yr, 1, 1)))

    call get_date(day_end, yr, mon, day, hr, min, sec)
    mr_end = cs%num_cycle * (time_type_to_real(day_end - set_date(yr, 1, 1)) / &
                   time_type_to_real(set_date(yr+1, 1, 1) - set_date(yr, 1, 1)))

    ! The Chapeau functions are centered at whole integer values that are nominally
    ! the end of the month to enable simple conversion from the fractional-years times
    ! CS%num_cycle.

    ! The month-average temperatures have as an index the month number.

    m_end = periodic_int(real(ceiling(mr_end)), cs%num_cycle)
    m_mid = periodic_int(real(ceiling(mr_st)), cs%num_cycle)
    m_st = periodic_int(mr_st, cs%num_cycle)

    mr_st = periodic_real(mr_st, cs%num_cycle)
    mr_end = periodic_real(mr_end, cs%num_cycle)
      !  mr_mid = periodic_real(ceiling(mr_st), CS%num_cycle)
    mr_prev = periodic_real(real(floor(mr_st)), cs%num_cycle)
    mr_next = periodic_real(real(m_end), cs%num_cycle)
    if (m_mid == m_end) then ; mr_mid = mr_end ! There is only one cell.
    else ; mr_mid = periodic_real(real(m_mid), cs%num_cycle) ; endif

    ! There may be two cells that run from mr_st to mr_mid and mr_mid to mr_end.

    ! The values of m for weights are all calculated relative to mr_prev, so
    ! check whether mr_mid, etc., need to be shifted by CS%num_cycle, so that these
    ! values satisfiy  mr_prev <= mr_st < mr_mid <= mr_end <= mr_next.
    if (mr_st < mr_prev) mr_prev = mr_prev - cs%num_cycle
    if (mr_mid < mr_st) mr_mid = mr_mid + cs%num_cycle
    if (mr_end < mr_st) mr_end = mr_end + cs%num_cycle
    if (mr_next < mr_prev) mr_next = mr_next + cs%num_cycle

    !### These might be removed later - they are to check the coding.
    if ((mr_mid < mr_st) .or. (mr_mid > mr_prev + 1.)) call mom_error(fatal, &
          "apply ctrl_forcing: m_mid interpolation out of bounds; fix the code.")
    if ((mr_end < mr_st) .or. (mr_end > mr_prev + 2.)) call mom_error(fatal, &
          "apply ctrl_forcing: m_end interpolation out of bounds; fix the code.")
    if (mr_end > mr_next) call mom_error(fatal, &
          "apply ctrl_forcing: mr_next interpolation out of bounds; fix the code.")

    wt_per1 = 1.0
    if (mr_mid < mr_end) wt_per1 = (mr_mid - mr_st) / (mr_end - mr_st)

    ! Find the 3 Chapeau-function weights, bearing in mind that m_end may be m_mid.
    wt_st = wt_per1 * (1. + (mr_prev - 0.5*(mr_st + mr_mid)))
    wt_end = (1.0-wt_per1) * (1. + (0.5*(mr_end + mr_mid) - mr_next))
    wt_mid = 1.0 - (wt_st + wt_end)
    if ((wt_st < 0.0) .or. (wt_end < 0.0) .or. (wt_mid < 0.0)) &
      call mom_error(fatal, "apply_ctrl_forcing: Negative m weights")
    if ((wt_st > 1.0) .or. (wt_end > 1.0) .or. (wt_mid > 1.0)) &
      call mom_error(fatal, "apply_ctrl_forcing: Excessive m weights")

    ! Add to vert_heat and vert_precip.
    do j=js,je ; do i=is,ie
      virt_heat(i,j) = virt_heat(i,j) + (wt_st * cs%heat_cyc(i,j,m_st) + &
                        (wt_mid * cs%heat_cyc(i,j,m_mid) + &
                         wt_end * cs%heat_cyc(i,j,m_end)))
      virt_precip(i,j) = virt_precip(i,j) + (wt_st * cs%precip_cyc(i,j,m_st) + &
                        (wt_mid * cs%precip_cyc(i,j,m_mid) + &
                         wt_end * cs%precip_cyc(i,j,m_end)))
    enddo ; enddo

    ! If different from the last period, take the average and determine the
    ! chapeau weighting

    ! The Chapeau functions are centered at whole integer values that are nominally
    ! the end of the month to enable simple conversion from the fractional-years times
    ! CS%num_cycle.

    ! The month-average temperatures have as an index the month number, so the averages
    ! apply to indicies m_end and m_mid.

    if (cs%avg_time(m_end) <= 0.0) then ! zero out the averages.
      cs%avg_time(m_end) = 0.0
      do j=js,je ; do i=is,ie
        cs%avg_SST_anom(i,j,m_end) = 0.0
        cs%avg_SSS_anom(i,j,m_end) = 0.0 ; cs%avg_SSS(i,j,m_end) = 0.0
      enddo ; enddo
    endif
    if (cs%avg_time(m_mid) <= 0.0) then ! zero out the averages.
      cs%avg_time(m_mid) = 0.0
      do j=js,je ; do i=is,ie
        cs%avg_SST_anom(i,j,m_mid) = 0.0
        cs%avg_SSS_anom(i,j,m_mid) = 0.0 ; cs%avg_SSS(i,j,m_mid) = 0.0
      enddo ; enddo
    endif

    ! Accumulate the average anomalies for this period.
    dt_wt = wt_per1 * dt
    cs%avg_time(m_mid) = cs%avg_time(m_mid) + dt_wt
    do j=js,je ; do i=is,ie
      cs%avg_SST_anom(i,j,m_mid) = cs%avg_SST_anom(i,j,m_mid) + &
                                   dt_wt * g%mask2dT(i,j) * sst_anom(i,j)
      cs%avg_SSS_anom(i,j,m_mid) = cs%avg_SSS_anom(i,j,m_mid) + &
                                   dt_wt * g%mask2dT(i,j) * sss_anom(i,j)
      cs%avg_SSS(i,j,m_mid) = cs%avg_SSS(i,j,m_mid) + dt_wt * sss_mean(i,j)
    enddo ; enddo
    if (wt_per1 < 1.0) then
      dt_wt = (1.0-wt_per1) * dt
      cs%avg_time(m_end) = cs%avg_time(m_end) + dt_wt
      do j=js,je ; do i=is,ie
        cs%avg_SST_anom(i,j,m_end) = cs%avg_SST_anom(i,j,m_end) + &
                                     dt_wt * g%mask2dT(i,j) * sst_anom(i,j)
        cs%avg_SSS_anom(i,j,m_end) = cs%avg_SSS_anom(i,j,m_end) + &
                                     dt_wt * g%mask2dT(i,j) * sss_anom(i,j)
        cs%avg_SSS(i,j,m_end) = cs%avg_SSS(i,j,m_end) + dt_wt * sss_mean(i,j)
      enddo ; enddo
    endif

    ! Update the Chapeau magnitudes for 4 cycles ago.
    m_u1 = periodic_int(m_st - 4.0, cs%num_cycle)
    m_u2 = periodic_int(m_st - 3.0, cs%num_cycle)
    m_u3 = periodic_int(m_st - 2.0, cs%num_cycle)

    if (cs%avg_time(m_u1) > 0.0) then
      do j=js,je ; do i=is,ie
        cs%avg_SST_anom(i,j,m_u1) = cs%avg_SST_anom(i,j,m_u1) / cs%avg_time(m_u1)
        cs%avg_SSS_anom(i,j,m_u1) = cs%avg_SSS_anom(i,j,m_u1) / cs%avg_time(m_u1)
        cs%avg_SSS(i,j,m_u1) = cs%avg_SSS(i,j,m_u1) / cs%avg_time(m_u1)
      enddo ; enddo
      cs%avg_time(m_u1) = -1.0
    endif
    if (cs%avg_time(m_u2) > 0.0) then
      do j=js,je ; do i=is,ie
        cs%avg_SST_anom(i,j,m_u2) = cs%avg_SST_anom(i,j,m_u2) / cs%avg_time(m_u2)
        cs%avg_SSS_anom(i,j,m_u2) = cs%avg_SSS_anom(i,j,m_u2) / cs%avg_time(m_u2)
        cs%avg_SSS(i,j,m_u2) = cs%avg_SSS(i,j,m_u2) / cs%avg_time(m_u2)
      enddo ; enddo
      cs%avg_time(m_u2) = -1.0
    endif
    if (cs%avg_time(m_u3) > 0.0) then
      do j=js,je ; do i=is,ie
        cs%avg_SST_anom(i,j,m_u3) = cs%avg_SST_anom(i,j,m_u3) / cs%avg_time(m_u3)
        cs%avg_SSS_anom(i,j,m_u3) = cs%avg_SSS_anom(i,j,m_u3) / cs%avg_time(m_u3)
        cs%avg_SSS(i,j,m_u3) = cs%avg_SSS(i,j,m_u3) / cs%avg_time(m_u3)
      enddo ; enddo
      cs%avg_time(m_u3) = -1.0
    endif

    dt1_heat_rate = wt_per1 * dt * cs%heat_cyc_rate
    dt1_prec_rate = wt_per1 * dt * cs%prec_cyc_rate
    dt2_heat_rate = (1.0-wt_per1) * dt * cs%heat_cyc_rate
    dt2_prec_rate = (1.0-wt_per1) * dt * cs%prec_cyc_rate

    if (wt_per1 < 1.0) then
      call pass_var(cs%heat_cyc(:,:,m_u2), g%Domain, complete=.false.)
      call pass_var(cs%precip_cyc(:,:,m_u2), g%Domain, complete=.false.)
    endif
    call pass_var(cs%heat_cyc(:,:,m_u1), g%Domain, complete=.false.)
    call pass_var(cs%precip_cyc(:,:,m_u1), g%Domain)

    if ((cs%avg_time(m_u1) == -1.0) .and. (cs%avg_time(m_u2) == -1.0)) then
      do j=js,je ; do i=is-1,ie
        coef = cs%Len2 * (g%dy_Cu(i,j)*g%IdxCu(i,j))
        flux_heat_x(i,j) = coef * (cs%heat_cyc(i,j,m_u1) - cs%heat_cyc(i+1,j,m_u1))
        flux_prec_x(i,j) = coef * (cs%precip_cyc(i,j,m_u1) - cs%precip_cyc(i+1,j,m_u1))
      enddo ; enddo
      do j=js-1,je ; do i=is,ie
        coef = cs%Len2 * (g%dx_Cv(i,j)*g%IdyCv(i,j))
        flux_heat_y(i,j) = coef * (cs%heat_cyc(i,j,m_u1) - cs%heat_cyc(i,j+1,m_u1))
        flux_prec_y(i,j) = coef * (cs%precip_cyc(i,j,m_u1) - cs%precip_cyc(i,j+1,m_u1))
      enddo ; enddo
      do j=js,je ; do i=is,ie
        cs%heat_cyc(i,j,m_u1) = cs%heat_cyc(i,j,m_u1) + dt1_heat_rate * ( &
           -cs%lam_cyc_heat*(cs%avg_SST_anom(i,j,m_u2) - cs%avg_SST_anom(i,j,m_u1)) + &
          (us%m_to_L**2*g%IareaT(i,j) * ((flux_heat_x(i-1,j) - flux_heat_x(i,j)) + &
                            (flux_heat_y(i,j-1) - flux_heat_y(i,j))) ) )

        cs%precip_cyc(i,j,m_u1) = cs%precip_cyc(i,j,m_u1) + dt1_prec_rate * ( &
          cs%lam_cyc_prec * (cs%avg_SSS_anom(i,j,m_u2) - cs%avg_SSS_anom(i,j,m_u1)) / &
                            (0.5*(cs%avg_SSS(i,j,m_u2) + cs%avg_SSS(i,j,m_u1))) + &
          (us%m_to_L**2*g%IareaT(i,j) * ((flux_prec_x(i-1,j) - flux_prec_x(i,j)) + &
                            (flux_prec_y(i,j-1) - flux_prec_y(i,j))) ) )
      enddo ; enddo
    endif

    if ((wt_per1 < 1.0) .and. (cs%avg_time(m_u1) == -1.0) .and. (cs%avg_time(m_u2) == -1.0))  then
      do j=js,je ; do i=is-1,ie
        coef = cs%Len2 * (g%dy_Cu(i,j)*g%IdxCu(i,j))
        flux_heat_x(i,j) = coef * (cs%heat_cyc(i,j,m_u2) - cs%heat_cyc(i+1,j,m_u2))
        flux_prec_x(i,j) = coef * (cs%precip_cyc(i,j,m_u2) - cs%precip_cyc(i+1,j,m_u2))
      enddo ; enddo
      do j=js-1,je ; do i=is,ie
        coef = cs%Len2 * (g%dx_Cv(i,j)*g%IdyCv(i,j))
        flux_heat_y(i,j) = coef * (cs%heat_cyc(i,j,m_u2) - cs%heat_cyc(i,j+1,m_u2))
        flux_prec_y(i,j) = coef * (cs%precip_cyc(i,j,m_u2) - cs%precip_cyc(i,j+1,m_u2))
      enddo ; enddo
      do j=js,je ; do i=is,ie
        cs%heat_cyc(i,j,m_u2) = cs%heat_cyc(i,j,m_u2) + dt1_heat_rate * ( &
         -cs%lam_cyc_heat*(cs%avg_SST_anom(i,j,m_u3) - cs%avg_SST_anom(i,j,m_u2)) + &
          (us%m_to_L**2*g%IareaT(i,j) * ((flux_heat_x(i-1,j) - flux_heat_x(i,j)) + &
                            (flux_heat_y(i,j-1) - flux_heat_y(i,j))) ) )

        cs%precip_cyc(i,j,m_u2) = cs%precip_cyc(i,j,m_u2) + dt1_prec_rate * ( &
          cs%lam_cyc_prec * (cs%avg_SSS_anom(i,j,m_u3) - cs%avg_SSS_anom(i,j,m_u2)) / &
                             (0.5*(cs%avg_SSS(i,j,m_u3) + cs%avg_SSS(i,j,m_u2))) + &
          (us%m_to_L**2*g%IareaT(i,j) * ((flux_prec_x(i-1,j) - flux_prec_x(i,j)) + &
                            (flux_prec_y(i,j-1) - flux_prec_y(i,j))) ) )
      enddo ; enddo
    endif

  endif ! (CS%num_cycle > 0)

end subroutine apply_ctrl_forcing

!> This function maps rval into an integer in the range from 1 to num_period.
function periodic_int(rval, num_period) result (m)
  real,    intent(in) :: rval       !< Input for mapping.
  integer, intent(in) :: num_period !< Maximum output.
  integer             :: m          !< Return value.

  m = floor(rval)
  if (m <= 0) then
    m = m + num_period * (1 + (abs(m) / num_period))
  elseif (m > num_period) then
    m = m - num_period * ((m-1) / num_period)
  endif
end function

!> This function shifts rval by an integer multiple of num_period so that
!! 0 <= val_out < num_period.
function periodic_real(rval, num_period) result(val_out)
  real,    intent(in) :: rval       !< Input to be shifted into valid range.
  integer, intent(in) :: num_period !< Maximum valid value.
  real                :: val_out    !< Return value.
  integer :: nshft

  if (rval < 0) then ; nshft = floor(abs(rval) / num_period) + 1
  elseif (rval < num_period) then ; nshft = 0
  else ; nshft = -1*floor(rval / num_period) ; endif

  val_out = rval + nshft * num_period
end function


!> This subroutine is used to allocate and register any fields in this module
!! that should be written to or read from the restart file.
subroutine register_ctrl_forcing_restarts(G, param_file, CS, restart_CS)
  type(ocean_grid_type), intent(in) :: g          !< The ocean's grid structure.
  type(param_file_type), intent(in) :: param_file !< A structure indicating the
                                                  !! open file to parse for model
                                                  !! parameter values.
  type(ctrl_forcing_cs), pointer :: cs            !< A pointer that is set to point to the
                                                  !! control structure for this module.
  type(mom_restart_cs),  pointer :: restart_cs    !< A pointer to the restart control structure.

  logical :: controlled, use_temperature
  character (len=8) :: period_str
  type(vardesc) :: vd
  integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

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

  controlled = .false.
  call read_param(param_file, "CONTROLLED_FORCING", controlled)
  if (.not.controlled) return

  use_temperature = .true.
  call read_param(param_file, "ENABLE_THERMODYNAMICS", use_temperature)
  if (.not.use_temperature) call mom_error(fatal, &
    "register_ctrl_forcing_restarts: CONTROLLED_FORCING only works with "//&
    "ENABLE_THERMODYNAMICS defined.")

  allocate(cs)

  cs%do_integrated = .true. ; cs%num_cycle = 0
  call read_param(param_file, "CTRL_FORCE_INTEGRATED", cs%do_integrated)
  call read_param(param_file, "CTRL_FORCE_NUM_CYCLE", cs%num_cycle)

  if (cs%do_integrated) then
    call safe_alloc_ptr(cs%heat_0,isd,ied,jsd,jed) ; cs%heat_0(:,:) = 0.0
    call safe_alloc_ptr(cs%precip_0,isd,ied,jsd,jed) ; cs%precip_0(:,:) = 0.0
    vd = var_desc("Ctrl_heat","W m-2","Control Integrative Heating",z_grid='1')
    call register_restart_field(cs%heat_0, vd, .false., restart_cs)
    vd = var_desc("Ctrl_precip","kg m-2 s-1","Control Integrative Precipitation",z_grid='1')
    call register_restart_field(cs%precip_0, vd, .false., restart_cs)
  endif

  if (cs%num_cycle > 0) then
    write (period_str, '(i8)') cs%num_cycle
    period_str = trim('p ')//trim(adjustl(period_str))
    call safe_alloc_ptr(cs%heat_cyc,isd,ied,jsd,jed,cs%num_cycle) ; cs%heat_cyc(:,:,:) = 0.0
    call safe_alloc_ptr(cs%precip_cyc,isd,ied,jsd,jed,cs%num_cycle) ; cs%precip_cyc(:,:,:) = 0.0
    vd = var_desc("Ctrl_heat_cycle", "W m-2","Cyclical Control Heating",&
                  z_grid='1', t_grid=period_str)
    call register_restart_field(cs%heat_cyc, vd, .false., restart_cs)
    vd = var_desc("Ctrl_precip_cycle","kg m-2 s-1","Cyclical Control Precipitation", &
                  z_grid='1', t_grid=period_str)
    call register_restart_field(cs%precip_cyc, vd, .false., restart_cs)

    call safe_alloc_ptr(cs%avg_time,cs%num_cycle) ; cs%avg_time(:) = 0.0
    vd = var_desc("avg_time","sec","Cyclical accumulated averaging time", &
                  '1',z_grid='1',t_grid=period_str)
    call register_restart_field(cs%avg_time, vd, .false., restart_cs)

    call safe_alloc_ptr(cs%avg_SST_anom,isd,ied,jsd,jed,cs%num_cycle) ; cs%avg_SST_anom(:,:,:) = 0.0
    call safe_alloc_ptr(cs%avg_SSS_anom,isd,ied,jsd,jed,cs%num_cycle) ; cs%avg_SSS_anom(:,:,:) = 0.0
    vd = var_desc("avg_SST_anom","deg C","Cyclical average SST Anomaly", &
                  z_grid='1',t_grid=period_str)
    call register_restart_field(cs%avg_SST_anom, vd, .false., restart_cs)
    vd = var_desc("avg_SSS_anom","g kg-1","Cyclical average SSS Anomaly", &
                  z_grid='1',t_grid=period_str)
    call register_restart_field(cs%avg_SSS_anom, vd, .false., restart_cs)
  endif

end subroutine register_ctrl_forcing_restarts

!> Set up this modules control structure.
subroutine controlled_forcing_init(Time, G, param_file, diag, CS)
  type(time_type),           intent(in) :: time       !< The current model time.
  type(ocean_grid_type),     intent(in) :: g          !< The ocean's grid structure.
  type(param_file_type),     intent(in) :: param_file !< A structure indicating the
                                                      !! open file to parse for model
                                                      !! parameter values.
  type(diag_ctrl), target,   intent(in) :: diag       !< A structure that is used to regulate
                                                      !! diagnostic output.
  type(ctrl_forcing_cs),     pointer    :: cs         !< A pointer that is set to point to the
                                                      !! control structure for this module.
  real :: smooth_len
  logical :: do_integrated
  integer :: num_cycle
! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_controlled_forcing" ! This module's name.

  ! These should have already been called.
  ! call read_param(param_file, "CTRL_FORCE_INTEGRATED", CS%do_integrated)
  ! call read_param(param_file, "CTRL_FORCE_NUM_CYCLE", CS%num_cycle)

  if (associated(cs)) then
    do_integrated = cs%do_integrated ; num_cycle = cs%num_cycle
  else
    do_integrated = .false. ; num_cycle = 0
  endif

  ! Read all relevant parameters and write them to the model log.
  call log_version(param_file, mdl, version, "")
  call log_param(param_file, mdl, "CTRL_FORCE_INTEGRATED", do_integrated, &
                 "If true, use a PI controller to determine the surface "//&
                 "forcing that is consistent with the observed mean properties.", &
                 default=.false.)
  call log_param(param_file, mdl, "CTRL_FORCE_NUM_CYCLE", num_cycle, &
                 "The number of cycles per year in the controlled forcing, "//&
                 "or 0 for no cyclic forcing.", default=0)

  if (.not.associated(cs)) return

  cs%diag => diag

  call get_param(param_file, mdl, "CTRL_FORCE_HEAT_INT_RATE", cs%heat_int_rate, &
                 "The integrated rate at which heat flux anomalies are "//&
                 "accumulated.", units="s-1", default=0.0)
  call get_param(param_file, mdl, "CTRL_FORCE_PREC_INT_RATE", cs%prec_int_rate, &
                 "The integrated rate at which precipitation anomalies "//&
                 "are accumulated.", units="s-1", default=0.0)
  call get_param(param_file, mdl, "CTRL_FORCE_HEAT_CYC_RATE", cs%heat_cyc_rate, &
                 "The integrated rate at which cyclical heat flux "//&
                 "anomalies are accumulated.", units="s-1", default=0.0)
  call get_param(param_file, mdl, "CTRL_FORCE_PREC_CYC_RATE", cs%prec_cyc_rate, &
                 "The integrated rate at which cyclical precipitation "//&
                 "anomalies are accumulated.", units="s-1", default=0.0)
  call get_param(param_file, mdl, "CTRL_FORCE_SMOOTH_LENGTH", smooth_len, &
                 "The length scales over which controlled forcing "//&
                 "anomalies are smoothed.", units="m", default=0.0)
  call get_param(param_file, mdl, "CTRL_FORCE_LAMDA_HEAT", cs%lam_heat, &
                 "A constant of proportionality between SST anomalies "//&
                 "and controlling heat fluxes", "W m-2 K-1", default=0.0)
  call get_param(param_file, mdl, "CTRL_FORCE_LAMDA_PREC", cs%lam_prec, &
                 "A constant of proportionality between SSS anomalies "//&
                 "(normalised by mean SSS) and controlling precipitation.", &
                 "kg m-2", default=0.0)
  call get_param(param_file, mdl, "CTRL_FORCE_LAMDA_CYC_HEAT", cs%lam_cyc_heat, &
                 "A constant of proportionality between SST anomalies "//&
                 "and cyclical controlling heat fluxes", "W m-2 K-1", default=0.0)
  call get_param(param_file, mdl, "CTRL_FORCE_LAMDA_CYC_PREC", cs%lam_cyc_prec, &
                 "A constant of proportionality between SSS anomalies "//&
                 "(normalised by mean SSS) and cyclical controlling "//&
                 "precipitation.", "kg m-2", default=0.0)

  cs%Len2 = smooth_len**2

! ### REPLACE THIS WITH ANY DIAGNOSTICS FROM THIS MODULE.
!  CS%id_taux = register_diag_field('ocean_model', 'taux', diag%axesu1, Time, &
!        'Zonal Wind Stress', 'Pascal')

end subroutine controlled_forcing_init

!> Clean up this modules control structure.
subroutine controlled_forcing_end(CS)
  type(ctrl_forcing_cs),    pointer :: cs !< A pointer to the control structure
                                          !! returned by a previous call to
                                          !! controlled_forcing_init, it will be
                                          !! deallocated here.

  if (associated(cs)) then
    if (associated(cs%heat_0))       deallocate(cs%heat_0)
    if (associated(cs%precip_0))     deallocate(cs%precip_0)
    if (associated(cs%heat_cyc))     deallocate(cs%heat_cyc)
    if (associated(cs%precip_cyc))   deallocate(cs%precip_cyc)
    if (associated(cs%avg_SST_anom)) deallocate(cs%avg_SST_anom)
    if (associated(cs%avg_SSS_anom)) deallocate(cs%avg_SSS_anom)
    if (associated(cs%avg_SSS))      deallocate(cs%avg_SSS)

    deallocate(cs)
  endif
  cs => null()

end subroutine controlled_forcing_end

!> \namespace mom_controlled_forcing
!!                                                                     *
!!  By Robert Hallberg, July 2011                                      *
!!                                                                     *
!!    This program contains the subroutines that use control-theory    *
!!  to adjust the surface heat flux and precipitation, based on the    *
!!  time-mean or periodically (seasonally) varying anomalies from the  *
!!  observed state.  The techniques behind this are described in       *
!!  Hallberg and Adcroft (2011, in prep.).                             *
!!                                                                     *
!!  Macros written all in capital letters are defined in MOM_memory.h. *
!!                                                                     *
!!     A small fragment of the grid is shown below:                    *
!!                                                                     *
!!    j+1  x ^ x ^ x   At x:  q                                        *
!!    j+1  > o > o >   At ^:  v, tauy                                  *
!!    j    x ^ x ^ x   At >:  u, taux                                  *
!!    j    > o > o >   At o:  h, fluxes.                               *
!!    j-1  x ^ x ^ x                                                   *
!!        i-1  i  i+1  At x & ^:                                       *
!!           i  i+1    At > & o:                                       *
!!                                                                     *
!!  The boundaries always run through q grid points (x).               *
end module mom_controlled_forcing