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_variables, only : thermo_var_ptrs
use mom_verticalgrid, only : verticalgrid_type
use mom_eos, only : calculate_density, calculate_density_derivs
 
implicit none ; private
 
#include <MOM_memory.h>
 
public geothermal, 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 [kg m-3 degC-1].
  real, pointer :: geo_heat(:,:) => null() !< The geothermal heat flux [W m-2].
  real    :: geothermal_thick !< The thickness over which geothermal heating is
                             !! applied [m] (not [H]).
  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.
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 ### 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(h, tv, dt, ea, eb, G, GV, 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 [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(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 [kg m-3]
    p_ref        ! coordiante densities reference pressure [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 [kg m-3 degC-1]
    drcv_ds_    ! partial derivative of coordinate density wrt saln [kg m-3 ppt-1]
 
  real :: angstrom, h_neglect  ! small thicknesses [H ~> m or kg m-2]
  real :: rcv           ! coordinate density of present layer [kg m-3]
  real :: rcv_tgt       ! coordinate density of target layer [kg m-3]
  real :: drcv          ! difference between Rcv and Rcv_tgt [kg m-3]
  real :: drcv_dt       ! partial derivative of coordinate density wrt temp
                        ! in the present layer [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 m2 J-1 ~> degC m3 J-1 or degC kg J-1]
 
  logical :: do_i(szi_(g))
  integer :: i, j, k, is, ie, js, je, nz, k2, i2
  integer :: isj, iej, num_start, 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_kg_m2 * tv%C_p)
  angstrom  = gv%Angstrom_H
  h_neglect = gv%H_subroundoff
  p_ref(:)  = tv%P_Ref
 
  if (.not.associated(tv%T)) call mom_error(fatal, "MOM geothermal: "//&
      "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
 
!$OMP parallel do default(none) shared(is,ie,js,je,G,GV,CS,dt,Irho_cp,nkmb,tv,    &
!$OMP                                  p_Ref,h,Angstrom,nz,H_neglect,eb)          &
!$OMP                          private(num_start,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_start = 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_start = num_start + 1
      h_geo_rem(i) = cs%Geothermal_thick * gv%m_to_H
    enddo
    if (num_start == 0) cycle
    num_left = num_start
 
    ! 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(:), isj, iej-isj+1, tv%eqn_of_state)
    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_, 1, 2, tv%eqn_of_state)
            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_transfer + h_neglect)
            tv%T(i,j,k_tgt) = (h(i,j,k_tgt) * 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) * 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
      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_kg_m2 * &
           (g%mask2dT(i,j) * (cs%geo_heat(i,j) * (dt*irho_cp)) - heat_rem(i))
    enddo ; endif
  enddo ! j-loop
 
!  do i=is,ie ; do j=js,je
!    resid(i,j) = tv%internal_heat(i,j) - resid(i,j) - GV%H_to_kg_m2 * &
!           (G%mask2dT(i,j) * (CS%geo_heat(i,j) * (dt*Irho_cp)))
!  enddo ; enddo
 
end subroutine geothermal
 
!> Initialize parameters and allocate memory associated with the geothermal heating module.
subroutine geothermal_init(Time, G, param_file, diag, CS)
  type(time_type), target, intent(in)    :: time !< Current model time.
  type(ocean_grid_type),   intent(inout) :: g    !< The ocean's grid structure.
  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.
! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_geothermal"  ! module name
  ! Local variables
  character(len=200) :: inputdir, geo_file, filename, geotherm_var
  real :: scale
  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", 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)
  cs%apply_geothermal = .not.(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)
  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",  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) * 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) * scale
    enddo ; enddo
  endif
  call pass_var(cs%geo_heat, g%domain)
 
  ! post the static geothermal heating field
  id = register_static_field('ocean_model', 'geo_heat', diag%axesT1,   &
        'Geothermal heat flux into ocean', 'W m-2',                    &
        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.)
 
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