MOM_marine_ice.F90

!> Routines incorporating the effects of marine ice (sea-ice and icebergs) into
!! the ocean model dynamics and thermodynamics.
module mom_marine_ice

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

use mom_constants,     only : hlf
use mom_diag_mediator, only : post_data, query_averaging_enabled, diag_ctrl
use mom_domains,       only : pass_var, pass_vector, agrid, bgrid_ne, cgrid_ne
use mom_domains,       only : to_all, omit_corners
use mom_error_handler, only : mom_error, fatal, warning
use mom_error_handler, only : mom_error, mom_mesg, fatal, warning
use mom_file_parser,   only : get_param, log_version, param_file_type
use mom_forcing_type,  only : allocate_forcing_type
use mom_forcing_type,  only : forcing, mech_forcing
use mom_grid,          only : ocean_grid_type
use mom_time_manager,  only : time_type
use mom_unit_scaling,  only : unit_scale_type
use mom_variables,     only : surface

implicit none ; private

#include <MOM_memory.h>

public iceberg_forces, iceberg_fluxes, marine_ice_init

!> Control structure for MOM_marine_ice
type, public :: marine_ice_cs ; private
  real :: kv_iceberg          !< The viscosity of the icebergs [L4 Z-2 T-1 ~> m2 s-1] (for ice rigidity)
  real :: berg_area_threshold !< Fraction of grid cell which iceberg must occupy
                              !! so that fluxes below are set to zero. (0.5 is a
                              !! good value to use.) Not applied for negative values.
  real :: latent_heat_fusion  !< Latent heat of fusion [Q ~> J kg-1]
  real :: density_iceberg     !< A typical density of icebergs [R ~> kg m-3] (for ice rigidity)

  type(time_type), pointer :: time !< A pointer to the ocean model's clock.
  type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the timing of diagnostic output.
end type marine_ice_cs

contains

!> add_berg_flux_to_shelf adds rigidity and ice-area coverage due to icebergs
!! to the forces type fields, and adds ice-areal coverage and modifies various
!! thermodynamic fluxes due to the presence of icebergs.
subroutine iceberg_forces(G, forces, use_ice_shelf, sfc_state, time_step, CS)
  type(ocean_grid_type), intent(inout) :: G       !< The ocean's grid structure
  type(mech_forcing),    intent(inout) :: forces  !< A structure with the driving mechanical forces
  type(surface),         intent(inout) :: sfc_state !< A structure containing fields that
                                                    !! describe the surface state of the ocean.
  logical,               intent(in)    :: use_ice_shelf  !< If true, this configuration uses ice shelves.
  real,                  intent(in)    :: time_step  !< The coupling time step [s].
  type(marine_ice_cs),   pointer       :: CS      !< Pointer to the control structure for MOM_marine_ice

  real :: kv_rho_ice ! The viscosity of ice divided by its density [L4 Z-2 T-1 R-1 ~> m5 kg-1 s-1].
  integer :: i, j, is, ie, js, je
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
  !This routine adds iceberg data to the ice shelf data (if ice shelf is used)
  !which can then be used to change the top of ocean boundary condition used in
  !the ocean model. This routine is taken from the add_shelf_flux subroutine
  !within the ice shelf model.

  if (.not.associated(cs)) return

  if (.not.(associated(forces%area_berg) .and.  associated(forces%mass_berg) ) ) return

  if (.not.(associated(forces%frac_shelf_u) .and. associated(forces%frac_shelf_v) .and. &
            associated(forces%rigidity_ice_u) .and. associated(forces%rigidity_ice_v)) ) return

  ! This section sets or augments the values of fields in forces.
  if (.not. use_ice_shelf) then
    forces%frac_shelf_u(:,:) = 0.0 ; forces%frac_shelf_v(:,:) = 0.0
  endif
  if (.not. forces%accumulate_rigidity) then
    forces%rigidity_ice_u(:,:) = 0.0 ; forces%rigidity_ice_v(:,:) = 0.0
  endif

  call pass_var(forces%area_berg, g%domain, to_all+omit_corners, halo=1, complete=.false.)
  call pass_var(forces%mass_berg, g%domain, to_all+omit_corners, halo=1, complete=.true.)
  kv_rho_ice = cs%kv_iceberg / cs%density_iceberg
  do j=js,je ; do i=is-1,ie
    if ((g%areaT(i,j) + g%areaT(i+1,j) > 0.0)) & ! .and. (G%dxdy_u(I,j) > 0.0)) &
      forces%frac_shelf_u(i,j) = forces%frac_shelf_u(i,j) + &
           (forces%area_berg(i,j)*g%areaT(i,j) + forces%area_berg(i+1,j)*g%areaT(i+1,j)) / &
           (g%areaT(i,j) + g%areaT(i+1,j))
    forces%rigidity_ice_u(i,j) = forces%rigidity_ice_u(i,j) + kv_rho_ice * &
                        min(forces%mass_berg(i,j), forces%mass_berg(i+1,j))
  enddo ; enddo
  do j=js-1,je ; do i=is,ie
    if ((g%areaT(i,j) + g%areaT(i,j+1) > 0.0)) & ! .and. (G%dxdy_v(i,J) > 0.0)) &
      forces%frac_shelf_v(i,j) = forces%frac_shelf_v(i,j) + &
           (forces%area_berg(i,j)*g%areaT(i,j) + forces%area_berg(i,j+1)*g%areaT(i,j+1)) / &
           (g%areaT(i,j) + g%areaT(i,j+1))
    forces%rigidity_ice_v(i,j) = forces%rigidity_ice_v(i,j) + kv_rho_ice * &
                         min(forces%mass_berg(i,j), forces%mass_berg(i,j+1))
  enddo ; enddo

end subroutine iceberg_forces

!> iceberg_fluxes adds ice-area-coverage and modifies various
!! thermodynamic fluxes due to the presence of icebergs.
subroutine iceberg_fluxes(G, US, fluxes, use_ice_shelf, sfc_state, time_step, CS)
  type(ocean_grid_type), intent(inout) :: G       !< The ocean's grid structure
  type(unit_scale_type), intent(in)    :: US      !< A dimensional unit scaling type
  type(forcing),         intent(inout) :: fluxes  !< A structure with pointers to themodynamic,
                                                  !! tracer and mass exchange forcing fields
  type(surface),         intent(inout) :: sfc_state !< A structure containing fields that
                                                    !! describe the surface state of the ocean.
  logical,               intent(in)    :: use_ice_shelf  !< If true, this configuration uses ice shelves.
  real,                  intent(in)    :: time_step   !< The coupling time step [s].
  type(marine_ice_cs),   pointer       :: CS      !< Pointer to the control structure for MOM_marine_ice

  real :: fraz      ! refreezing rate [R Z T-1 ~> kg m-2 s-1]
  real :: I_dt_LHF  ! The inverse of the timestep times the latent heat of fusion times [Q-1 T-1 ~> kg J-1 s-1].
  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
  isd = g%isd ; jsd = g%jsd ; ied = g%ied ; jed = g%jed
  ! This routine adds iceberg data to the ice shelf data (if ice shelf is used)
  ! which can then be used to change the top of ocean boundary condition used in
  ! the ocean model. This routine is taken from the add_shelf_flux subroutine
  ! within the ice shelf model.

  if (.not.associated(cs)) return
  if (.not.(associated(fluxes%area_berg) .and. associated(fluxes%ustar_berg) .and. &
            associated(fluxes%mass_berg) ) ) return
  if (.not.(associated(fluxes%frac_shelf_h) .and. associated(fluxes%ustar_shelf)) ) return


  if (.not.(associated(fluxes%area_berg) .and. associated(fluxes%ustar_berg) .and. &
            associated(fluxes%mass_berg) ) ) return
  if (.not. use_ice_shelf) then
    fluxes%frac_shelf_h(:,:) = 0.
    fluxes%ustar_shelf(:,:) = 0.
  endif
  do j=jsd,jed ; do i=isd,ied ; if (g%areaT(i,j) > 0.0) then
    fluxes%frac_shelf_h(i,j) = fluxes%frac_shelf_h(i,j) + fluxes%area_berg(i,j)
    fluxes%ustar_shelf(i,j)  = fluxes%ustar_shelf(i,j)  + fluxes%ustar_berg(i,j)
  endif ; enddo ; enddo

  !Zero'ing out other fluxes under the tabular icebergs
  if (cs%berg_area_threshold >= 0.) then
    i_dt_lhf = 1.0 / (us%s_to_T*time_step * cs%latent_heat_fusion)
    do j=jsd,jed ; do i=isd,ied
      if (fluxes%frac_shelf_h(i,j) > cs%berg_area_threshold) then
        ! Only applying for ice shelf covering most of cell.

        if (associated(fluxes%sw)) fluxes%sw(i,j) = 0.0
        if (associated(fluxes%lw)) fluxes%lw(i,j) = 0.0
        if (associated(fluxes%latent)) fluxes%latent(i,j) = 0.0
        if (associated(fluxes%evap)) fluxes%evap(i,j) = 0.0

        ! Add frazil formation diagnosed by the ocean model [Q R Z ~> J m-2] in the
        ! form of surface layer evaporation [R Z T-1 ~> kg m-2 s-1]. Update lprec in the
        ! control structure for diagnostic purposes.

        if (allocated(sfc_state%frazil)) then
          fraz = sfc_state%frazil(i,j) * i_dt_lhf
          if (associated(fluxes%evap))  fluxes%evap(i,j)  = fluxes%evap(i,j)  - fraz
        ! if (associated(fluxes%lprec)) fluxes%lprec(i,j) = fluxes%lprec(i,j) - fraz
          sfc_state%frazil(i,j) = 0.0
        endif

        !Alon: Should these be set to zero too?
        if (associated(fluxes%sens)) fluxes%sens(i,j) = 0.0
        if (associated(fluxes%salt_flux)) fluxes%salt_flux(i,j) = 0.0
        if (associated(fluxes%lprec)) fluxes%lprec(i,j) = 0.0
      endif
    enddo ; enddo
  endif

end subroutine iceberg_fluxes

!> Initialize control structure for MOM_marine_ice
subroutine marine_ice_init(Time, G, param_file, diag, CS)
  type(time_type), target, intent(in)    :: Time !< Current model time
  type(ocean_grid_type),   intent(in)    :: G !< Ocean grid structure
  type(param_file_type),   intent(in)    :: param_file !< Runtime parameter handles
  type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure
  type(marine_ice_cs),     pointer       :: CS   !< Pointer to the control structure for MOM_marine_ice
! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_marine_ice"  ! This module's name.

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

  ! Write all relevant parameters to the model log.
  call log_version(mdl, version)

  call get_param(param_file, mdl, "KV_ICEBERG",  cs%kv_iceberg, &
                 "The viscosity of the icebergs",  &
                 units="m2 s-1", default=1.0e10, scale=g%US%Z_to_L**2*g%US%m_to_L**2*g%US%T_to_s)
  call get_param(param_file, mdl, "DENSITY_ICEBERGS",  cs%density_iceberg, &
                 "A typical density of icebergs.", units="kg m-3", default=917.0, scale=g%US%kg_m3_to_R)
  call get_param(param_file, mdl, "LATENT_HEAT_FUSION", cs%latent_heat_fusion, &
                 "The latent heat of fusion.", units="J/kg", default=hlf, scale=g%US%J_kg_to_Q)
  call get_param(param_file, mdl, "BERG_AREA_THRESHOLD", cs%berg_area_threshold, &
                 "Fraction of grid cell which iceberg must occupy, so that fluxes "//&
                 "below berg are set to zero. Not applied for negative "//&
                 "values.", units="non-dim", default=-1.0)

end subroutine marine_ice_init

end module mom_marine_ice