MOM_marine_ice.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Routines incorporating the effects of marine ice (sea-ice and icebergs) into
6!! the ocean model dynamics and thermodynamics.
8
9use mom_constants, only : hlf
11use mom_domains, only : pass_var, pass_vector, agrid, bgrid_ne, cgrid_ne
12use mom_domains, only : to_all, omit_corners
13use mom_error_handler, only : mom_error, fatal, warning
14use mom_error_handler, only : mom_error, mom_mesg, fatal, warning
15use mom_file_parser, only : get_param, log_version, param_file_type
16use mom_forcing_type, only : allocate_forcing_type
18use mom_grid, only : ocean_grid_type
19use mom_time_manager, only : time_type
21use mom_variables, only : surface
22
23implicit none ; private
24
25#include <MOM_memory.h>
26
28
29!> Control structure for MOM_marine_ice
30type, public :: marine_ice_cs ; private
31 real :: kv_iceberg !< The viscosity of the icebergs [L4 Z-2 T-1 ~> m2 s-1] (for ice rigidity)
32 real :: berg_area_threshold !< Fraction of grid cell which iceberg must occupy
33 !! so that fluxes below are set to zero [nondim]. (0.5 is a
34 !! good value to use.) Not applied for negative values.
35 real :: latent_heat_fusion !< Latent heat of fusion [Q ~> J kg-1]
36 real :: density_iceberg !< A typical density of icebergs [R ~> kg m-3] (for ice rigidity)
37
38 type(time_type), pointer :: time !< A pointer to the ocean model's clock.
39 type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the timing of diagnostic output.
40end type marine_ice_cs
41
42contains
43
44!> add_berg_flux_to_shelf adds rigidity and ice-area coverage due to icebergs
45!! to the forces type fields, and adds ice-areal coverage and modifies various
46!! thermodynamic fluxes due to the presence of icebergs.
47subroutine iceberg_forces(G, forces, use_ice_shelf, sfc_state, time_step, CS)
48 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure
49 type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces
50 type(surface), intent(inout) :: sfc_state !< A structure containing fields that
51 !! describe the surface state of the ocean.
52 logical, intent(in) :: use_ice_shelf !< If true, this configuration uses ice shelves.
53 real, intent(in) :: time_step !< The coupling time step [T ~> s].
54 type(marine_ice_cs), pointer :: cs !< Pointer to the control structure for MOM_marine_ice
55
56 real :: kv_rho_ice ! The viscosity of ice divided by its density [L4 Z-2 T-1 R-1 ~> m5 kg-1 s-1].
57 integer :: i, j, is, ie, js, je
58 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
59 !This routine adds iceberg data to the ice shelf data (if ice shelf is used)
60 !which can then be used to change the top of ocean boundary condition used in
61 !the ocean model. This routine is taken from the add_shelf_flux subroutine
62 !within the ice shelf model.
63
64 if (.not.associated(cs)) return
65
66 if (.not.(associated(forces%area_berg) .and. associated(forces%mass_berg) ) ) return
67
68 if (.not.(associated(forces%frac_shelf_u) .and. associated(forces%frac_shelf_v) .and. &
69 associated(forces%rigidity_ice_u) .and. associated(forces%rigidity_ice_v)) ) return
70
71 ! This section sets or augments the values of fields in forces.
72 if (.not. use_ice_shelf) then
73 forces%frac_shelf_u(:,:) = 0.0 ; forces%frac_shelf_v(:,:) = 0.0
74 endif
75 if (.not. forces%accumulate_rigidity) then
76 forces%rigidity_ice_u(:,:) = 0.0 ; forces%rigidity_ice_v(:,:) = 0.0
77 endif
78
79 call pass_var(forces%area_berg, g%domain, to_all+omit_corners, halo=1, complete=.false.)
80 call pass_var(forces%mass_berg, g%domain, to_all+omit_corners, halo=1, complete=.true.)
81 kv_rho_ice = cs%kv_iceberg / cs%density_iceberg
82 do j=js,je ; do i=is-1,ie
83 if ((g%areaT(i,j) + g%areaT(i+1,j) > 0.0)) & ! .and. (G%dxdy_u(I,j) > 0.0)) &
84 forces%frac_shelf_u(i,j) = forces%frac_shelf_u(i,j) + &
85 ((forces%area_berg(i,j)*g%areaT(i,j)) + (forces%area_berg(i+1,j)*g%areaT(i+1,j))) / &
86 (g%areaT(i,j) + g%areaT(i+1,j))
87 forces%rigidity_ice_u(i,j) = forces%rigidity_ice_u(i,j) + kv_rho_ice * &
88 min(forces%mass_berg(i,j), forces%mass_berg(i+1,j))
89 enddo ; enddo
90 do j=js-1,je ; do i=is,ie
91 if ((g%areaT(i,j) + g%areaT(i,j+1) > 0.0)) & ! .and. (G%dxdy_v(i,J) > 0.0)) &
92 forces%frac_shelf_v(i,j) = forces%frac_shelf_v(i,j) + &
93 ((forces%area_berg(i,j)*g%areaT(i,j)) + (forces%area_berg(i,j+1)*g%areaT(i,j+1))) / &
94 (g%areaT(i,j) + g%areaT(i,j+1))
95 forces%rigidity_ice_v(i,j) = forces%rigidity_ice_v(i,j) + kv_rho_ice * &
96 min(forces%mass_berg(i,j), forces%mass_berg(i,j+1))
97 enddo ; enddo
98
99end subroutine iceberg_forces
100
101!> iceberg_fluxes adds ice-area-coverage and modifies various
102!! thermodynamic fluxes due to the presence of icebergs.
103subroutine iceberg_fluxes(G, US, fluxes, use_ice_shelf, sfc_state, time_step, CS)
104 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure
105 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
106 type(forcing), intent(inout) :: fluxes !< A structure with pointers to themodynamic,
107 !! tracer and mass exchange forcing fields
108 type(surface), intent(inout) :: sfc_state !< A structure containing fields that
109 !! describe the surface state of the ocean.
110 logical, intent(in) :: use_ice_shelf !< If true, this configuration uses ice shelves.
111 real, intent(in) :: time_step !< The coupling time step [T ~> s].
112 type(marine_ice_cs), pointer :: cs !< Pointer to the control structure for MOM_marine_ice
113
114 real :: fraz ! refreezing rate [R Z T-1 ~> kg m-2 s-1]
115 real :: i_dt_lhf ! The inverse of the timestep times the latent heat of fusion [Q-1 T-1 ~> kg J-1 s-1].
116 integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
117 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
118 isd = g%isd ; jsd = g%jsd ; ied = g%ied ; jed = g%jed
119 ! This routine adds iceberg data to the ice shelf data (if ice shelf is used)
120 ! which can then be used to change the top of ocean boundary condition used in
121 ! the ocean model. This routine is taken from the add_shelf_flux subroutine
122 ! within the ice shelf model.
123
124 if (.not.associated(cs)) return
125 if (.not.(associated(fluxes%area_berg) .and. associated(fluxes%ustar_berg) .and. &
126 associated(fluxes%mass_berg) ) ) return
127 if (.not.(associated(fluxes%frac_shelf_h) .and. associated(fluxes%ustar_shelf)) ) return
128
129
130 if (.not.(associated(fluxes%area_berg) .and. associated(fluxes%ustar_berg) .and. &
131 associated(fluxes%mass_berg) ) ) return
132 if (.not. use_ice_shelf) then
133 fluxes%frac_shelf_h(:,:) = 0.
134 fluxes%ustar_shelf(:,:) = 0.
135 endif
136 do j=jsd,jed ; do i=isd,ied ; if (g%areaT(i,j) > 0.0) then
137 fluxes%frac_shelf_h(i,j) = fluxes%frac_shelf_h(i,j) + fluxes%area_berg(i,j)
138 fluxes%ustar_shelf(i,j) = fluxes%ustar_shelf(i,j) + fluxes%ustar_berg(i,j)
139 endif ; enddo ; enddo
140
141 !Zero'ing out other fluxes under the tabular icebergs
142 if (cs%berg_area_threshold >= 0.) then
143 i_dt_lhf = 1.0 / (time_step * cs%latent_heat_fusion)
144 do j=jsd,jed ; do i=isd,ied
145 if (fluxes%frac_shelf_h(i,j) > cs%berg_area_threshold) then
146 ! Only applying for ice shelf covering most of cell.
147
148 if (associated(fluxes%sw)) fluxes%sw(i,j) = 0.0
149 if (associated(fluxes%lw)) fluxes%lw(i,j) = 0.0
150 if (associated(fluxes%latent)) fluxes%latent(i,j) = 0.0
151 if (associated(fluxes%evap)) fluxes%evap(i,j) = 0.0
152
153 ! Add frazil formation diagnosed by the ocean model [Q R Z ~> J m-2] in the
154 ! form of surface layer evaporation [R Z T-1 ~> kg m-2 s-1]. Update lprec in the
155 ! control structure for diagnostic purposes.
156
157 if (allocated(sfc_state%frazil)) then
158 fraz = sfc_state%frazil(i,j) * i_dt_lhf
159 if (associated(fluxes%evap)) fluxes%evap(i,j) = fluxes%evap(i,j) - fraz
160 ! if (associated(fluxes%lprec)) fluxes%lprec(i,j) = fluxes%lprec(i,j) - fraz
161 sfc_state%frazil(i,j) = 0.0
162 endif
163
164 !Alon: Should these be set to zero too?
165 if (associated(fluxes%sens)) fluxes%sens(i,j) = 0.0
166 if (associated(fluxes%salt_flux)) fluxes%salt_flux(i,j) = 0.0
167 if (associated(fluxes%lprec)) fluxes%lprec(i,j) = 0.0
168 endif
169 enddo ; enddo
170 endif
171
172end subroutine iceberg_fluxes
173
174!> Initialize control structure for MOM_marine_ice
175subroutine marine_ice_init(Time, G, param_file, diag, CS)
176 type(time_type), target, intent(in) :: time !< Current model time
177 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
178 type(param_file_type), intent(in) :: param_file !< Runtime parameter handles
179 type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure
180 type(marine_ice_cs), pointer :: cs !< Pointer to the control structure for MOM_marine_ice
181
182 ! This include declares and sets the variable "version".
183# include "version_variable.h"
184 character(len=40) :: mdl = "MOM_marine_ice" ! This module's name.
185
186 if (associated(cs)) then
187 call mom_error(warning, "marine_ice_init called with an associated control structure.")
188 return
189 else ; allocate(cs) ; endif
190
191 ! Write all relevant parameters to the model log.
192 call log_version(mdl, version)
193
194 call get_param(param_file, mdl, "KV_ICEBERG", cs%kv_iceberg, &
195 "The viscosity of the icebergs", &
196 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)
197 call get_param(param_file, mdl, "DENSITY_ICEBERGS", cs%density_iceberg, &
198 "A typical density of icebergs.", units="kg m-3", default=917.0, scale=g%US%kg_m3_to_R)
199 call get_param(param_file, mdl, "LATENT_HEAT_FUSION", cs%latent_heat_fusion, &
200 "The latent heat of fusion.", units="J/kg", default=hlf, scale=g%US%J_kg_to_Q)
201 call get_param(param_file, mdl, "BERG_AREA_THRESHOLD", cs%berg_area_threshold, &
202 "Fraction of grid cell which iceberg must occupy, so that fluxes "//&
203 "below berg are set to zero. Not applied for negative values.", &
204 units="nondim", default=-1.0)
205
206end subroutine marine_ice_init
207
208end module mom_marine_ice