mom_geothermal Module Reference

Data Types
type	geothermal_cs

Functions/Subroutines
subroutine, public	geothermal (h, tv, dt, ea, eb, G, GV, CS)
	This subroutine 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)? More...
subroutine, public	geothermal_init (Time, G, param_file, diag, CS)
subroutine, public	geothermal_end (CS)

Function/Subroutine Documentation

geothermal()

subroutine, public mom_geothermal::geothermal	(	real, dimension(szi_(g),szj_(g),szk_(g)), intent(inout)	h,
		type(thermo_var_ptrs), intent(inout)	tv,
		real, intent(in)	dt,
		real, dimension(szi_(g),szj_(g),szk_(g)), intent(inout)	ea,
		real, dimension(szi_(g),szj_(g),szk_(g)), intent(inout)	eb,
		type(ocean_grid_type), intent(inout)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(geothermal_cs), pointer	CS
	)

This subroutine 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)?

Parameters

[in,out]	g	The ocean's grid structure.
[in]	gv	The ocean's vertical grid structure.
[in,out]	h	Layer thicknesses, in H (usually m or kg m-2).
[in,out]	tv	A structure containing pointers to any available thermodynamic fields. Absent fields have NULL ptrs.
[in]	dt	Time increment, in s.
[in,out]	ea	The amount of fluid moved downward into a layer; this should be increased due to mixed layer detrainment, in the same units as h - usually m or kg m-2 (i.e., H).
[in,out]	eb	The amount of fluid moved upward into a layer; this should be increased due to mixed layer entrainment, in the same units as h - usually m or kg m-2 (i.e., H)
	cs	The control structure returned by a previous call to geothermal_init.

Definition at line 95 of file MOM_geothermal.F90.

References mom_eos::calculate_density_derivs(), and mom_error_handler::mom_error().

Referenced by mom_diabatic_driver::diabatic().

  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, in H
                                                                !! (usually 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, in 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, in the same
                                                                !! units as h - usually m or kg m-2
                                                                !! (i.e., H).
  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, in the same units as
                                                                !! h - usually m or kg m-2 (i.e., H)
  type(geothermal_cs),                      pointer       :: cs !< The control structure returned by
                                                                !! a previous call to
                                                                !! geothermal_init.

!   This subroutine 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)?

! Arguments: h - Layer thickness, in m or kg m-2. (Intent in/out)
!                The units of h are referred to as H below.
!  (in/out)  tv - A structure containing pointers to any available
!                 thermodynamic fields. Absent fields have NULL ptrs.
!  (in)      dt - Time increment, in s.
!  (in/out)  ea - The amount of fluid moved downward into a layer; this should
!                 be increased due to mixed layer detrainment, in the same units
!                 as h - usually m or kg m-2 (i.e., H).
!  (in/out)  eb - The amount of fluid moved upward into a layer; this should
!                 be increased due to mixed layer entrainment, in the same units
!                 as h - usually m or kg m-2 (i.e., H).
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS - The control structure returned by a previous call to
!                 geothermal_init.

!  real :: resid(SZI_(G),SZJ_(G)) !z1l: never been used.

  real, dimension(SZI_(G)) :: &
    heat_rem,  & ! remaining heat (H * degC)
    h_geo_rem, & ! remaining thickness to apply geothermal heating (units of H)
    rcv_bl,    & ! coordinate density in the deepest variable density layer (kg/m3)
    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 K-1)
    drcv_ds_    ! partial derivative of coordinate density wrt saln (kg m-3 ppt-1)

  real :: angstrom, h_neglect  ! small thicknesses in H
  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 K-1); usually negative
  real :: h_heated      ! thickness that is being heated (units of H)
  real :: heat_avail    ! heating available for the present layer (units of Kelvin * H)
  real :: heat_in_place ! heating to warm present layer w/o movement between layers (K * H)
  real :: heat_trans    ! heating available to move water from present layer to target layer (K * H)
  real :: heating       ! heating used to move water from present layer to target layer (K * H)
                        ! 0 <= heating <= heat_trans
  real :: h_transfer    ! thickness moved between layers (units of H)
  real :: wt_in_place   ! relative weighting that goes from 0 to 1 (non-dim)
  real :: i_h           ! inverse thickness (units of 1/H)
  real :: dtemp         ! temperature increase in a layer (Kelvin)
  real :: irho_cp       ! inverse of heat capacity per unit layer volume (units K H m2 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 (.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_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

Here is the call graph for this function:

Here is the caller graph for this function:

geothermal_end()

subroutine, public mom_geothermal::geothermal_end

(

type(geothermal_cs), pointer

CS

)

Definition at line 473 of file MOM_geothermal.F90.

  type(geothermal_cs), pointer :: cs

  if (associated(cs%geo_heat)) deallocate(cs%geo_heat)
  if (associated(cs)) deallocate(cs)

geothermal_init()

subroutine, public mom_geothermal::geothermal_init	(	type(time_type), intent(in), target	Time,
		type(ocean_grid_type), intent(in)	G,
		type(param_file_type), intent(in)	param_file,
		type(diag_ctrl), intent(inout), target	diag,
		type(geothermal_cs), pointer	CS
	)

Parameters

[in]	time	Current model time.
[in]	g	The ocean's grid structure.
[in]	param_file	A structure to parse for run-time parameters.
[in,out]	diag	Structure used to regulate diagnostic output.
	cs	Pointer pointing to the module control structure.

Definition at line 383 of file MOM_geothermal.F90.

References mom_error_handler::mom_error(), and mom_diag_mediator::register_static_field().

Referenced by mom_diabatic_driver::diabatic_driver_init().

  type(time_type), target, intent(in)    :: time !< 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 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.

! Arguments:
!  (in)      Time       - current model time
!  (in)      G          - ocean grid structure
!  (in)      param_file - structure indicating the open file to parse for
!                         model parameter values
!  (in)      diag       - structure used to regulate diagnostic output
!  (in/out)  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
  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 \n"//&
                 "factor for the heat flux read from GEOTHERMAL_FILE, or \n"//&
                 "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 \n"//&
                 "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 \n"//&
                 "simply heats water in place instead of moving it between \n"//&
                 "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 \n"//&
                 "GEOTHERMAL_FILE.", default="geo_heat")
    call read_data(filename, trim(geotherm_var), cs%geo_heat, &
                   domain=g%Domain%mpp_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

  ! 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')
  if (id > 0) call post_data(id, cs%geo_heat, diag, .true.)

Here is the call graph for this function:

Here is the caller graph for this function: