MOM_geothermal.F90

!> Implemented geothermal heating at the ocean bottom.
module mom_geothermal

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

use mom_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
use mom_diag_mediator, only : register_static_field, time_type, diag_ctrl
use mom_domains,       only : pass_var
use mom_error_handler, only : mom_error, fatal, warning
use mom_file_parser,   only : get_param, log_param, log_version, param_file_type
use mom_io,            only : mom_read_data, slasher
use mom_grid,          only : ocean_grid_type
use mom_unit_scaling,  only : unit_scale_type
use mom_variables,     only : thermo_var_ptrs
use mom_verticalgrid,  only : verticalgrid_type, get_thickness_units
use mom_eos,           only : calculate_density, calculate_density_derivs

implicit none ; private

#include <MOM_memory.h>

public geothermal_entraining, geothermal_in_place, geothermal_init, geothermal_end

!> Control structure for geothermal heating
type, public :: geothermal_cs ; private
  real    :: drcv_dt_inplace !<   The value of dRcv_dT above which (dRcv_dT is
                             !! negative) the water is heated in place instead
                             !! of moving upward between layers [R degC-1 ~> kg m-3 degC-1].
  real, pointer :: geo_heat(:,:) => null() !< The geothermal heat flux [J m-2 T-1 ~> W m-2].
  real    :: geothermal_thick !< The thickness over which geothermal heating is
                             !! applied [H ~> m or kg m-2].
  logical :: apply_geothermal !< If true, geothermal heating will be applied
                             !! otherwise GEOTHERMAL_SCALE has been set to 0 and
                             !! there is no heat to apply.

  type(time_type), pointer :: time => null() !< A pointer to the ocean model's clock.
  type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
                                             !! regulate the timing of diagnostic output.
  integer :: id_internal_heat_heat_tendency = -1   !< ID for diagnostic of heat tendency
  integer :: id_internal_heat_temp_tendency = -1   !< ID for diagnostic of temperature tendency
  integer :: id_internal_heat_h_tendency = -1      !< ID for diagnostic of thickness tendency

end type geothermal_cs

contains

!> Applies geothermal heating, including the movement of water
!! between isopycnal layers to match the target densities.  The heating is
!! applied to the bottommost layers that occur within GEOTHERMAL_THICKNESS of the bottom. If
!! the partial derivative of the coordinate density with temperature is positive
!! or very small, the layers are simply heated in place.  Any heat that can not
!! be applied to the ocean is returned (WHERE)?
subroutine geothermal_entraining(h, tv, dt, ea, eb, G, GV, US, CS, halo)
  type(ocean_grid_type),                    intent(inout) :: g  !< The ocean's grid structure.
  type(verticalgrid_type),                  intent(in)    :: gv !< The ocean's vertical grid structure.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h  !< Layer thicknesses [H ~> m or kg m-2]
  type(thermo_var_ptrs),                    intent(inout) :: tv !< A structure containing pointers
                                                                !! to any available thermodynamic
                                                                !! fields. Absent fields have NULL
                                                                !! ptrs.
  real,                                     intent(in)    :: dt !< Time increment [T ~> s].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: ea !< The amount of fluid moved
                                                                !! downward into a layer; this
                                                                !! should be increased due to mixed
                                                                !! layer detrainment [H ~> m or kg m-2]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: eb !< The amount of fluid moved upward
                                                                !! into a layer; this should be
                                                                !! increased due to mixed layer
                                                                !! entrainment [H ~> m or kg m-2].
  type(unit_scale_type),                    intent(in)    :: us !< A dimensional unit scaling type
  type(geothermal_cs),                      pointer       :: cs !< The control structure returned by
                                                                !! a previous call to
                                                                !! geothermal_init.
  integer,                        optional, intent(in)    :: halo !< Halo width over which to work
  ! Local variables
  real, dimension(SZI_(G)) :: &
    heat_rem,  & ! remaining heat [H degC ~> m degC or kg degC m-2]
    h_geo_rem, & ! remaining thickness to apply geothermal heating [H ~> m or kg m-2]
    rcv_bl,    & ! coordinate density in the deepest variable density layer [R ~> kg m-3]
    p_ref        ! coordinate densities reference pressure [R L2 T-2 ~> Pa]

  real, dimension(2) :: &
    t2, s2, &   ! temp and saln in the present and target layers [degC] and [ppt]
    drcv_dt_, & ! partial derivative of coordinate density wrt temp [R degC-1 ~> kg m-3 degC-1]
    drcv_ds_    ! partial derivative of coordinate density wrt saln [R ppt-1 ~> kg m-3 ppt-1]

  real :: angstrom, h_neglect  ! small thicknesses [H ~> m or kg m-2]
  real :: rcv           ! coordinate density of present layer [R ~> kg m-3]
  real :: rcv_tgt       ! coordinate density of target layer [R ~> kg m-3]
  real :: drcv          ! difference between Rcv and Rcv_tgt [R ~> kg m-3]
  real :: drcv_dt       ! partial derivative of coordinate density wrt temp
                        ! in the present layer [R degC-1 ~> kg m-3 degC-1]; usually negative
  real :: h_heated      ! thickness that is being heated [H ~> m or kg m-2]
  real :: heat_avail    ! heating available for the present layer [degC H ~> degC m or degC kg m-2]
  real :: heat_in_place ! heating to warm present layer w/o movement between layers
                        ! [degC H ~> degC m or degC kg m-2]
  real :: heat_trans    ! heating available to move water from present layer to target
                        ! layer [degC H ~> degC m or degC kg m-2]
  real :: heating       ! heating used to move water from present layer to target layer
                        ! [degC H ~> degC m or degC kg m-2]
                        ! 0 <= heating <= heat_trans
  real :: h_transfer    ! thickness moved between layers [H ~> m or kg m-2]
  real :: wt_in_place   ! relative weighting that goes from 0 to 1 [nondim]
  real :: i_h           ! inverse thickness [H-1 ~> m-1 or m2 kg-1]
  real :: dtemp         ! temperature increase in a layer [degC]
  real :: irho_cp       ! inverse of heat capacity per unit layer volume
                        ! [degC H Q-1 R-1 Z-1 ~> degC m3 J-1 or degC kg J-1]

  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: t_old   ! Temperature of each layer
                                                      ! before any heat is added,
                                                      ! for diagnostics [degC]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_old   ! Thickness of each layer
                                                      ! before any heat is added,
                                                      ! for diagnostics [H ~> m or kg m-2]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: work_3d ! Scratch variable used to
                                                      ! calculate change in heat
                                                      ! due to geothermal
  real :: idt           ! inverse of the timestep [T-1 ~> s-1]

  logical :: do_i(szi_(g))
  logical :: compute_h_old, compute_t_old
  integer :: i, j, k, is, ie, js, je, nz, k2, i2
  integer :: isj, iej, num_left, nkmb, k_tgt

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  if (present(halo)) then
    is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
  endif

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

  nkmb      = gv%nk_rho_varies
  irho_cp   = 1.0 / (gv%H_to_RZ * tv%C_p)
  angstrom  = gv%Angstrom_H
  h_neglect = gv%H_subroundoff
  p_ref(:)  = tv%P_Ref
  idt       = 1.0 / dt

  if (.not.associated(tv%T)) call mom_error(fatal, "MOM geothermal_entraining: "//&
      "Geothermal heating can only be applied if T & S are state variables.")

!  do j=js,je ; do i=is,ie
!    resid(i,j) = tv%internal_heat(i,j)
!  enddo ; enddo

  ! Conditionals for tracking diagnostic depdendencies
  compute_h_old = cs%id_internal_heat_h_tendency > 0 &
                  .or. cs%id_internal_heat_heat_tendency > 0 &
                  .or. cs%id_internal_heat_temp_tendency > 0

  compute_t_old = cs%id_internal_heat_heat_tendency > 0 &
                  .or. cs%id_internal_heat_temp_tendency > 0

  if (cs%id_internal_heat_heat_tendency > 0) work_3d(:,:,:) = 0.0

  if (compute_h_old .or. compute_t_old) then ; do k=1,nz ; do j=js,je ; do i=is,ie
    ! Save temperature and thickness before any changes are made (for diagnostics)
    h_old(i,j,k) = h(i,j,k)
    t_old(i,j,k) = tv%T(i,j,k)
  enddo ; enddo ; enddo ; endif

!$OMP parallel do default(none) shared(is,ie,js,je,G,GV,US,CS,dt,Irho_cp,nkmb,tv, &
!$OMP                                  p_Ref,h,Angstrom,nz,H_neglect,eb,          &
!$OMP                                  h_old,T_old,work_3d,Idt)                   &
!$OMP                          private(heat_rem,do_i,h_geo_rem,num_left,          &
!$OMP                                  isj,iej,Rcv_BL,h_heated,heat_avail,k_tgt,  &
!$OMP                                  Rcv_tgt,Rcv,dRcv_dT,T2,S2,dRcv_dT_,        &
!$OMP                                  dRcv_dS_,heat_in_place,heat_trans,         &
!$OMP                                  wt_in_place,dTemp,dRcv,h_transfer,heating, &
!$OMP                                  I_h)

  do j=js,je
    ! 1. Only work on columns that are being heated.
    ! 2. Find the deepest layer with any mass.
    ! 3. Find the partial derivative of locally referenced potential density
    !  and coordinate density with temperature, and the density of the layer
    !  and the layer above.
    ! 4. Heat a portion of the bottommost layer until it matches the target
    !    density of the layer above, and move it.
    ! 4a. In the case of variable density layers, heat but do not move.
    ! 5. If there is still heat left over, repeat for the next layer up.
    ! This subroutine updates thickness, T & S, and increments eb accordingly.

    ! 6. If there is not enough mass in the ocean, pass some of the heat up
    !    from the ocean via the frazil field?

    num_left = 0
    do i=is,ie
      heat_rem(i) = g%mask2dT(i,j) * (cs%geo_heat(i,j) * (dt*irho_cp))
      do_i(i) = .true. ; if (heat_rem(i) <= 0.0) do_i(i) = .false.
      if (do_i(i)) num_left = num_left + 1
      h_geo_rem(i) = cs%Geothermal_thick
    enddo
    if (num_left == 0) cycle

    ! Find the first and last columns that need to be worked on.
    isj = ie+1 ; do i=is,ie ; if (do_i(i)) then ; isj = i ; exit ; endif ; enddo
    iej = is-1 ; do i=ie,is,-1 ; if (do_i(i)) then ; iej = i ; exit ; endif ; enddo

    if (nkmb > 0) then
      call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref(:), rcv_bl(:), &
                             tv%eqn_of_state, (/isj-(g%isd-1),iej-(g%isd-1)/) )
    else
      rcv_bl(:) = -1.0
    endif

    do k=nz,1,-1
      do i=isj,iej ; if (do_i(i)) then

        if (h(i,j,k) > angstrom) then
          if ((h(i,j,k)-angstrom) >= h_geo_rem(i)) then
            h_heated = h_geo_rem(i)
            heat_avail = heat_rem(i)
            h_geo_rem(i) = 0.0
          else
            h_heated = (h(i,j,k)-angstrom)
            heat_avail = heat_rem(i) * (h_heated / &
                                        (h_geo_rem(i) + h_neglect))
            h_geo_rem(i) = h_geo_rem(i) - h_heated
          endif

          if (k<=nkmb .or. nkmb<=0) then
            ! Simply heat the layer; convective adjustment occurs later
            ! if necessary.
            k_tgt = k
          elseif ((k==nkmb+1) .or. (gv%Rlay(k-1) < rcv_bl(i))) then
            ! Add enough heat to match the lowest buffer layer density.
            k_tgt = nkmb
            rcv_tgt = rcv_bl(i)
          else
            ! Add enough heat to match the target density of layer k-1.
            k_tgt = k-1
            rcv_tgt = gv%Rlay(k-1)
          endif

          if (k<=nkmb .or. nkmb<=0) then
            rcv = 0.0 ; drcv_dt = 0.0 ! Is this OK?
          else
            call calculate_density(tv%T(i,j,k), tv%S(i,j,k), tv%P_Ref, &
                         rcv, tv%eqn_of_state)
            t2(1) = tv%T(i,j,k) ; s2(1) = tv%S(i,j,k)
            t2(2) = tv%T(i,j,k_tgt) ; s2(2) = tv%S(i,j,k_tgt)
            call calculate_density_derivs(t2(:), s2(:), p_ref(:), drcv_dt_, drcv_ds_, &
                         tv%eqn_of_state, (/1,2/) )
            drcv_dt = 0.5*(drcv_dt_(1) + drcv_dt_(2))
          endif

          if ((drcv_dt >= 0.0) .or. (k<=nkmb .or. nkmb<=0)) then
            ! This applies to variable density layers.
            heat_in_place = heat_avail
            heat_trans = 0.0
          elseif (drcv_dt <= cs%dRcv_dT_inplace) then
            ! This is the option that usually applies in isopycnal coordinates.
            heat_in_place = min(heat_avail, max(0.0, h(i,j,k) * &
                                            ((gv%Rlay(k)-rcv) / drcv_dt)))
            heat_trans = heat_avail - heat_in_place
          else
            ! wt_in_place should go from 0 to 1.
            wt_in_place = (cs%dRcv_dT_inplace - drcv_dt) / cs%dRcv_dT_inplace
            heat_in_place = max(wt_in_place*heat_avail, &
                                h(i,j,k) * ((gv%Rlay(k)-rcv) / drcv_dt) )
            heat_trans = heat_avail - heat_in_place
          endif

          if (heat_in_place > 0.0) then
            ! This applies to variable density layers. In isopycnal coordinates
            ! this only arises for relatively fresh water near the freezing
            ! point, in which case heating in place will eventually cause things
            ! to sort themselves out, if only because the water will warm to
            ! the temperature of maximum density.
            dtemp = heat_in_place / (h(i,j,k) + h_neglect)
            tv%T(i,j,k) = tv%T(i,j,k) + dtemp
            heat_rem(i) = heat_rem(i) - heat_in_place
            rcv = rcv + drcv_dt * dtemp
          endif

          if (heat_trans > 0.0) then
            ! The second expression might never be used, but will avoid
            ! division by 0.
            drcv = max(rcv - rcv_tgt, 0.0)

            !   dTemp = -dRcv / dRcv_dT
            !   h_transfer = min(heat_rem(i) / dTemp, h(i,j,k)-Angstrom)
            if ((-drcv_dt * heat_trans) >= drcv * (h(i,j,k)-angstrom)) then
              h_transfer = h(i,j,k) - angstrom
              heating = (h_transfer * drcv) / (-drcv_dt)
              ! Since not all the heat has been applied, return the fraction
              ! of the layer thickness that has not yet been fully heated to
              ! the h_geo_rem.
              h_geo_rem(i) = h_geo_rem(i) + h_heated * &
                      ((heat_avail - (heating + heat_in_place)) / heat_avail)
            else
              h_transfer = (-drcv_dt * heat_trans) / drcv
              heating = heat_trans
            endif
            heat_rem(i) = heat_rem(i) - heating

            i_h = 1.0 / ((h(i,j,k_tgt) + h_neglect) + h_transfer)
            tv%T(i,j,k_tgt) = ((h(i,j,k_tgt) + h_neglect) * tv%T(i,j,k_tgt) + &
                               (h_transfer * tv%T(i,j,k) + heating)) * i_h
            tv%S(i,j,k_tgt) = ((h(i,j,k_tgt) + h_neglect) * tv%S(i,j,k_tgt) + &
                               h_transfer * tv%S(i,j,k)) * i_h

            h(i,j,k) = h(i,j,k) - h_transfer
            h(i,j,k_tgt) = h(i,j,k_tgt) + h_transfer
            eb(i,j,k_tgt) = eb(i,j,k_tgt) + h_transfer
            if (k_tgt < k-1) then
              do k2 = k_tgt+1,k-1
                eb(i,j,k2) = eb(i,j,k2) + h_transfer
              enddo
            endif
          endif

          if (heat_rem(i) <= 0.0) then
            do_i(i) = .false. ; num_left = num_left-1
            ! For efficiency, uncomment these?
            ! if ((i==isj) .and. (num_left > 0)) then ; do i2=isj+1,iej ; if (do_i(i2)) then
            !   isj = i2 ; exit ! Set the new starting value.
            ! endif ; enddo ; endif
            ! if ((i==iej) .and. (num_left > 0)) then ; do i2=iej-1,isj,-1 ; if (do_i(i2)) then
            !   iej = i2 ; exit ! Set the new ending value.
            ! endif ; enddo ; endif
          endif
        endif

        ! Calculate heat tendency due to addition and transfer of internal heat
        if (cs%id_internal_heat_heat_tendency > 0) then
          work_3d(i,j,k) = ((gv%H_to_RZ*tv%C_p) * idt) * (h(i,j,k) * tv%T(i,j,k) - h_old(i,j,k) * t_old(i,j,k))
        endif

      endif ; enddo
      if (num_left <= 0) exit
    enddo ! k-loop

    if (associated(tv%internal_heat)) then ; do i=is,ie
      tv%internal_heat(i,j) = tv%internal_heat(i,j) + gv%H_to_RZ * &
           (g%mask2dT(i,j) * (cs%geo_heat(i,j) * (dt*irho_cp)) - heat_rem(i))
    enddo ; endif
  enddo ! j-loop

  ! Post diagnostic of 3D tendencies (heat, temperature, and thickness) due to internal heat
  if (cs%id_internal_heat_heat_tendency > 0) then
    call post_data(cs%id_internal_heat_heat_tendency, work_3d, cs%diag, alt_h=h_old)
  endif
  if (cs%id_internal_heat_temp_tendency > 0) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      work_3d(i,j,k) = idt * (tv%T(i,j,k) - t_old(i,j,k))
    enddo; enddo; enddo
    call post_data(cs%id_internal_heat_temp_tendency, work_3d, cs%diag, alt_h=h_old)
  endif
  if (cs%id_internal_heat_h_tendency > 0) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      work_3d(i,j,k) = idt * (h(i,j,k) - h_old(i,j,k))
    enddo ; enddo ; enddo
    call post_data(cs%id_internal_heat_h_tendency, work_3d, cs%diag, alt_h=h_old)
  endif

!  do j=js,je ; do i=is,ie
!    resid(i,j) = tv%internal_heat(i,j) - resid(i,j) - GV%H_to_RZ * &
!           (G%mask2dT(i,j) * (CS%geo_heat(i,j) * (dt*Irho_cp)))
!  enddo ; enddo

end subroutine geothermal_entraining

!> Applies geothermal heating to the bottommost layers that occur within GEOTHERMAL_THICKNESS of
!! the bottom, by simply heating the water in place.  Any heat that can not be applied to the ocean
!! is returned (WHERE)?
subroutine geothermal_in_place(h, tv, dt, G, GV, US, CS, halo)
  type(ocean_grid_type),                    intent(inout) :: g  !< The ocean's grid structure.
  type(verticalgrid_type),                  intent(in)    :: gv !< The ocean's vertical grid structure.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in)    :: h  !< Layer thicknesses [H ~> m or kg m-2]
  type(thermo_var_ptrs),                    intent(inout) :: tv !< A structure containing pointers
                                                                !! to any available thermodynamic
                                                                !! fields. Absent fields have NULL
                                                                !! ptrs.
  real,                                     intent(in)    :: dt !< Time increment [T ~> s].
  type(unit_scale_type),                    intent(in)    :: us !< A dimensional unit scaling type
  type(geothermal_cs),                      pointer       :: cs !< The control structure returned by
                                                                !! a previous call to geothermal_init.
  integer,                        optional, intent(in)    :: halo !< Halo width over which to work

  ! Local variables
  real, dimension(SZI_(G)) :: &
    heat_rem,  & ! remaining heat [H degC ~> m degC or kg degC m-2]
    h_geo_rem    ! remaining thickness to apply geothermal heating [H ~> m or kg m-2]

  real :: angstrom, h_neglect  ! small thicknesses [H ~> m or kg m-2]
  real :: heat_here     ! heating applied to the present layer [degC H ~> degC m or degC kg m-2]
  real :: dtemp         ! temperature increase in a layer [degC]
  real :: irho_cp       ! inverse of heat capacity per unit layer volume
                        ! [degC H Q-1 R-1 Z-1 ~> degC m3 J-1 or degC kg J-1]

  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
    dtdt_diag           ! Diagnostic of temperature tendency [degC T-1 ~> degC s-1] which might be
                        ! converted into a layer-integrated heat tendency [Q R Z T-1 ~> W m-2]
  real :: idt           ! inverse of the timestep [T-1 ~> s-1]
  logical :: do_any     ! True if there is more to be done on the current j-row.
  logical :: calc_diags ! True if diagnostic tendencies are needed.
  integer :: i, j, k, is, ie, js, je, nz, i2, isj, iej

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  if (present(halo)) then
    is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
  endif

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

  irho_cp   = 1.0 / (gv%H_to_RZ * tv%C_p)
  angstrom  = gv%Angstrom_H
  h_neglect = gv%H_subroundoff
  idt       = 1.0 / dt

  if (.not.associated(tv%T)) call mom_error(fatal, "MOM geothermal_in_place: "//&
      "Geothermal heating can only be applied if T & S are state variables.")

!  do i=is,ie ; do j=js,je
!    resid(i,j) = tv%internal_heat(i,j)
!  enddo ; enddo

  ! Conditionals for tracking diagnostic depdendencies
  calc_diags = (cs%id_internal_heat_heat_tendency > 0) .or. (cs%id_internal_heat_temp_tendency > 0)

  if (calc_diags) dtdt_diag(:,:,:) = 0.0

  !$OMP parallel do default(shared) private(heat_rem,do_any,h_geo_rem,isj,iej,heat_here,dTemp)
  do j=js,je
    ! Only work on columns that are being heated, and heat the near-bottom water.

    ! If there is not enough mass in the ocean, pass some of the heat up
    ! from the ocean via the frazil field?

    do_any = .false.
    do i=is,ie
      heat_rem(i) = g%mask2dT(i,j) * (cs%geo_heat(i,j) * (dt*irho_cp))
      if (heat_rem(i) > 0.0) do_any = .true.
      h_geo_rem(i) = cs%Geothermal_thick
    enddo
    if (.not.do_any) cycle

    ! Find the first and last columns that need to be worked on.
    isj = ie+1 ; do i=is,ie ; if (heat_rem(i) > 0.0) then ; isj = i ; exit ; endif ; enddo
    iej = is-1 ; do i=ie,is,-1 ; if (heat_rem(i) > 0.0) then ; iej = i ; exit ; endif ; enddo

    do k=nz,1,-1
      do_any = .false.
      do i=isj,iej
        if ((heat_rem(i) > 0.0) .and. (h(i,j,k) > angstrom)) then
          ! Apply some or all of the remaining heat to this layer.
          ! Convective adjustment occurs outside of this module if necessary.
          if ((h(i,j,k)-angstrom) >= h_geo_rem(i)) then
            heat_here = heat_rem(i)
            h_geo_rem(i) = 0.0
            heat_rem(i) = 0.0
          else
            heat_here = heat_rem(i) * ((h(i,j,k)-angstrom) / (h_geo_rem(i) + h_neglect))
            h_geo_rem(i) = h_geo_rem(i) - (h(i,j,k)-angstrom)
            heat_rem(i) = heat_rem(i) - heat_here
          endif

          dtemp = heat_here / (h(i,j,k) + h_neglect)
          tv%T(i,j,k) = tv%T(i,j,k) + dtemp
          if (calc_diags) dtdt_diag(i,j,k) = dtemp * idt
        endif

        if (heat_rem(i) > 0.0) do_any= .true.
      enddo

      if (.not.do_any) exit
    enddo ! k-loop

    if (associated(tv%internal_heat)) then ; do i=is,ie
      tv%internal_heat(i,j) = tv%internal_heat(i,j) + gv%H_to_RZ * &
           (g%mask2dT(i,j) * (cs%geo_heat(i,j) * (dt*irho_cp)) - heat_rem(i))
    enddo ; endif
  enddo ! j-loop

  ! Post diagnostics of 3D tendencies of heat and temperature due to geothermal heat
  if (cs%id_internal_heat_temp_tendency > 0) then
    call post_data(cs%id_internal_heat_temp_tendency, dtdt_diag, cs%diag, alt_h=h)
  endif
  if (cs%id_internal_heat_heat_tendency > 0) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      ! Dangerously reuse dTdt_diag for a related variable with different units, going from
      ! units of [degC T-1 ~> degC s-1] to units of [Q R Z T-1 ~> W m-2]
      dtdt_diag(i,j,k) = (gv%H_to_RZ*tv%C_p) * (h(i,j,k) * dtdt_diag(i,j,k))
    enddo ; enddo ; enddo
    call post_data(cs%id_internal_heat_heat_tendency, dtdt_diag, cs%diag, alt_h=h)
  endif

!  do j=js,je ; do i=is,ie
!    resid(i,j) = tv%internal_heat(i,j) - resid(i,j) - GV%H_to_RZ * &
!           (G%mask2dT(i,j) * (CS%geo_heat(i,j) * (dt*Irho_cp)))
!  enddo ; enddo

end subroutine geothermal_in_place

!> Initialize parameters and allocate memory associated with the geothermal heating module.
subroutine geothermal_init(Time, G, GV, US, param_file, diag, CS, useALEalgorithm)
  type(time_type), target, intent(in)    :: time !< Current model time.
  type(ocean_grid_type),   intent(inout) :: g    !< The ocean's grid structure.
  type(verticalgrid_type), intent(in)    :: gv   !< The ocean's vertical grid structure.
  type(unit_scale_type),   intent(in)    :: us   !< A dimensional unit scaling type
  type(param_file_type),   intent(in)    :: param_file !< A structure to parse for run-time
                                                 !! parameters.
  type(diag_ctrl), target, intent(inout) :: diag !< Structure used to regulate diagnostic output.
  type(geothermal_cs),     pointer       :: cs   !< Pointer pointing to the module control
                                                 !! structure.
  logical,       optional, intent(in)    :: usealealgorithm  !< logical for whether to use ALE remapping

! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_geothermal"  ! module name
  character(len=48)  :: thickness_units
  ! Local variables
  character(len=200) :: inputdir, geo_file, filename, geotherm_var
  real :: geo_scale  ! A constant heat flux or dimensionally rescaled geothermal flux scaling factor
                     ! [Q R Z T-1 ~> W m-2] or [Q R Z m2 s J-1 T-1 ~> 1]
  integer :: i, j, isd, ied, jsd, jed, id
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed

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

  cs%diag => diag
  cs%Time => time

  ! write parameters to the model log.
  call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl, "GEOTHERMAL_SCALE", geo_scale, &
                 "The constant geothermal heat flux, a rescaling "//&
                 "factor for the heat flux read from GEOTHERMAL_FILE, or "//&
                 "0 to disable the geothermal heating.", &
                 units="W m-2 or various", default=0.0, scale=us%W_m2_to_QRZ_T)
  cs%apply_geothermal = .not.(geo_scale == 0.0)
  if (.not.cs%apply_geothermal) return

  call safe_alloc_ptr(cs%geo_heat, isd, ied, jsd, jed) ; cs%geo_heat(:,:) = 0.0

  call get_param(param_file, mdl, "GEOTHERMAL_FILE", geo_file, &
                 "The file from which the geothermal heating is to be "//&
                 "read, or blank to use a constant heating rate.", default=" ")
  call get_param(param_file, mdl, "GEOTHERMAL_THICKNESS", cs%geothermal_thick, &
                 "The thickness over which to apply geothermal heating.", &
                 units="m", default=0.1, scale=gv%m_to_H)
  call get_param(param_file, mdl, "GEOTHERMAL_DRHO_DT_INPLACE", cs%dRcv_dT_inplace, &
                 "The value of drho_dT above which geothermal heating "//&
                 "simply heats water in place instead of moving it between "//&
                 "isopycnal layers.  This must be negative.", &
                 units="kg m-3 K-1", scale=us%kg_m3_to_R, default=-0.01)
  if (cs%dRcv_dT_inplace >= 0.0) call mom_error(fatal, "geothermal_init: "//&
         "GEOTHERMAL_DRHO_DT_INPLACE must be negative.")

  if (len_trim(geo_file) >= 1) then
    call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
    inputdir = slasher(inputdir)
    filename = trim(inputdir)//trim(geo_file)
    call log_param(param_file, mdl, "INPUTDIR/GEOTHERMAL_FILE", filename)
    call get_param(param_file, mdl, "GEOTHERMAL_VARNAME", geotherm_var, &
                 "The name of the geothermal heating variable in "//&
                 "GEOTHERMAL_FILE.", default="geo_heat")
    call mom_read_data(filename, trim(geotherm_var), cs%geo_heat, g%Domain)
    do j=jsd,jed ; do i=isd,ied
      cs%geo_heat(i,j) = (g%mask2dT(i,j) * geo_scale) * cs%geo_heat(i,j)
    enddo ; enddo
  else
    do j=jsd,jed ; do i=isd,ied
      cs%geo_heat(i,j) = g%mask2dT(i,j) * geo_scale
    enddo ; enddo
  endif
  call pass_var(cs%geo_heat, g%domain)

  thickness_units = get_thickness_units(gv)

  ! post the static geothermal heating field
  id = register_static_field('ocean_model', 'geo_heat', diag%axesT1,   &
        'Geothermal heat flux into ocean', 'W m-2', conversion=us%QRZ_T_to_W_m2, &
        cmor_field_name='hfgeou', cmor_units='W m-2',                  &
        cmor_standard_name='upward_geothermal_heat_flux_at_sea_floor', &
        cmor_long_name='Upward geothermal heat flux at sea floor', &
        x_cell_method='mean', y_cell_method='mean', area_cell_method='mean')
  if (id > 0) call post_data(id, cs%geo_heat, diag, .true.)

  ! Diagnostic for tendencies due to internal heat (in 3d)
  cs%id_internal_heat_heat_tendency=register_diag_field('ocean_model', &
        'internal_heat_heat_tendency', diag%axesTL, time,              &
        'Heat tendency (in 3D) due to internal (geothermal) sources',  &
        'W m-2', conversion=us%QRZ_T_to_W_m2, v_extensive=.true.)
  cs%id_internal_heat_temp_tendency=register_diag_field('ocean_model', &
        'internal_heat_temp_tendency', diag%axesTL, time,              &
        'Temperature tendency (in 3D) due to internal (geothermal) sources', &
        'degC s-1', conversion=us%s_to_T, v_extensive=.true.)
  if (present(usealealgorithm)) then ; if (.not.usealealgorithm) then
    ! Do not offer this diagnostic if heating will be in place.
    cs%id_internal_heat_h_tendency=register_diag_field('ocean_model',    &
        'internal_heat_h_tendency', diag%axesTL, time,                &
        'Thickness tendency (in 3D) due to internal (geothermal) sources', &
        trim(thickness_units), conversion=gv%H_to_MKS*us%s_to_T, v_extensive=.true.)
  endif ; endif

end subroutine geothermal_init

!> Clean up and deallocate memory associated with the geothermal heating module.
subroutine geothermal_end(CS)
  type(geothermal_cs), pointer :: cs !< Geothermal heating control structure that
                                     !! will be deallocated in this subroutine.

  if (associated(cs%geo_heat)) deallocate(cs%geo_heat)
  if (associated(cs)) deallocate(cs)
end subroutine geothermal_end

!> \namespace mom_geothermal
!!
!! Geothermal heating can be added either in a layered isopycnal mode, in which the heating raises the density
!! of the layer to the target density of the layer above, and then moves the water into that layer, or in a
!! simple Eulerian mode, in which the bottommost GEOTHERMAL_THICKNESS are heated.  Geothermal heating will also
!! provide a buoyant source of bottom TKE that can be used to further mix the near-bottom water. In cold fresh
!! water lakes where heating increases density, water should be moved into deeper layers, but this is not
!! implemented yet.

end module mom_geothermal