Skip to content

Commit

Permalink
+Create module MOM_marine_ice
Browse files Browse the repository at this point in the history
  Created a new module, MOM_marine_ice, to set dynamic properties for the ocean
due to icebergs and sea-ice.  The existing subroutine add_berg_flux_to_sHelf
has been moved into this new module, and there are also an init routine and a
control structure for this routine.  All answers are bitwise identical.
  • Loading branch information
Hallberg-NOAA committed May 3, 2018
1 parent ef4cd0b commit a1c5679
Show file tree
Hide file tree
Showing 2 changed files with 174 additions and 114 deletions.
118 changes: 4 additions & 114 deletions config_src/coupled_driver/ocean_model_MOM.F90
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,7 @@ module ocean_model_mod
use MOM_get_input, only : Get_MOM_Input, directories
use MOM_grid, only : ocean_grid_type
use MOM_io, only : close_file, file_exists, read_data, write_version_number
use MOM_marine_ice, only : add_berg_flux_to_shelf, marine_ice_init, marine_ice_CS
use MOM_restart, only : MOM_restart_CS, save_restart
use MOM_string_functions, only : uppercase
use MOM_surface_forcing, only : surface_forcing_init, convert_IOB_to_fluxes
Expand Down Expand Up @@ -210,6 +211,9 @@ module ocean_model_mod
Ice_shelf_CSp => NULL() !< A pointer to the control structure for the
!! ice shelf model that couples with MOM6. This
!! is null if there is no ice shelf.
type(marine_ice_CS), pointer :: &
marine_ice_CSp => NULL() !< A pointer to the control structure for the
!! marine ice effects module.
type(wave_parameters_cs), pointer :: &
Waves !< A structure containing pointers to the surface wave fields
type(surface_forcing_CS), pointer :: &
Expand Down Expand Up @@ -388,8 +392,6 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn)
OS%diag, OS%forces, OS%fluxes)
endif
if (OS%icebergs_apply_rigid_boundary) then
!call allocate_forcing_type(OS%grid, OS%fluxes, iceberg=.true.)
!This assumes that the iceshelf and ocean are on the same grid. I hope this is true
if (.not. OS%use_ice_shelf) &
call allocate_forcing_type(OS%grid, OS%fluxes, shelf=.true.)
endif
Expand Down Expand Up @@ -687,118 +689,6 @@ end subroutine update_ocean_model
! </SUBROUTINE> NAME="update_ocean_model"


subroutine add_berg_flux_to_shelf(G, forces, fluxes, use_ice_shelf, density_ice, kv_ice, &
latent_heat_fusion, sfc_state, time_step, berg_area_threshold)
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(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) :: kv_ice !< The viscosity of ice, in m2 s-1.
real, intent(in) :: density_ice !< A typical density of ice, in kg m-3.
real, intent(in) :: latent_heat_fusion !< The latent heat of fusion, in J kg-1.
real, intent(in) :: time_step !< The coupling time step, in s.
real, intent(in) :: berg_area_threshold !< Area threshold for zeroing fluxes below iceberg
! Arguments:
! (in) fluxes - A structure of surface fluxes that may be used.
! (in) G - The ocean's grid structure.
real :: fraz ! refreezing rate in kg m-2 s-1
real :: I_dt_LHF ! The inverse of the timestep times the latent heat of fusion, in kg J-1 s-1.
real :: kv_rho_ice ! The viscosity of ice divided by its density, in m5 kg-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(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

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

! 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
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 = kv_ice / density_ice
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
!### This halo update may be unnecessary. Test it. -RWH
call pass_vector(forces%frac_shelf_u, forces%frac_shelf_v, G%domain, TO_ALL, CGRID_NE)

! The remaining code sets or augments the values of fields in fluxes.

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 (berg_area_threshold >= 0.) then
I_dt_LHF = 1.0 / (time_step * latent_heat_fusion)
do j=jsd,jed ; do i=isd,ied
if (fluxes%frac_shelf_h(i,j) > 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 (J m-2) in the
! form of surface layer evaporation (kg m-2 s-1). Update lprec in the
! control structure for diagnostic purposes.

if (associated(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
!CS%lprec(i,j)=CS%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 add_berg_flux_to_shelf

!=======================================================================
! <SUBROUTINE NAME="ocean_model_restart">
!
Expand Down
170 changes: 170 additions & 0 deletions src/ice_shelf/MOM_marine_ice.F90
Original file line number Diff line number Diff line change
@@ -0,0 +1,170 @@
!> 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_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_variables, only : surface

implicit none ; private

#include <MOM_memory.h>

public add_berg_flux_to_shelf, marine_ice_init

!> Control structure for MOM_marine_ice
type, public :: marine_ice_CS ; private
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 add_berg_flux_to_shelf(G, forces, fluxes, use_ice_shelf, density_ice, kv_ice, &
latent_heat_fusion, sfc_state, time_step, berg_area_threshold)
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(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) :: kv_ice !< The viscosity of ice, in m2 s-1.
real, intent(in) :: density_ice !< A typical density of ice, in kg m-3.
real, intent(in) :: latent_heat_fusion !< The latent heat of fusion, in J kg-1.
real, intent(in) :: time_step !< The coupling time step, in s.
real, intent(in) :: berg_area_threshold !< Area threshold for zeroing fluxes below iceberg
! Arguments:
! (in) fluxes - A structure of surface fluxes that may be used.
! (in) G - The ocean's grid structure.
real :: fraz ! refreezing rate in kg m-2 s-1
real :: I_dt_LHF ! The inverse of the timestep times the latent heat of fusion, in kg J-1 s-1.
real :: kv_rho_ice ! The viscosity of ice divided by its density, in m5 kg-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(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

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

! 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
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 = kv_ice / density_ice
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
!### This halo update may be unnecessary. Test it. -RWH
call pass_vector(forces%frac_shelf_u, forces%frac_shelf_v, G%domain, TO_ALL, CGRID_NE)

! The remaining code sets or augments the values of fields in fluxes.

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 (berg_area_threshold >= 0.) then
I_dt_LHF = 1.0 / (time_step * latent_heat_fusion)
do j=jsd,jed ; do i=isd,ied
if (fluxes%frac_shelf_h(i,j) > 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 (J m-2) in the
! form of surface layer evaporation (kg m-2 s-1). Update lprec in the
! control structure for diagnostic purposes.

if (associated(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
!CS%lprec(i,j)=CS%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 add_berg_flux_to_shelf

!> 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 !< 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)

end subroutine marine_ice_init

end module MOM_marine_ice

0 comments on commit a1c5679

Please sign in to comment.