user_shelf_init.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!> This module specifies the initial values and evolving properties of the
6!! MOM6 ice shelf, using user-provided code.
8
9use mom_error_handler, only : mom_mesg, mom_error, fatal, is_root_pe
10use mom_file_parser, only : get_param, log_version, param_file_type
11use mom_grid, only : ocean_grid_type
12use mom_time_manager, only : time_type, set_time, time_type_to_real
14
15implicit none ; private
16
17#include <MOM_memory.h>
18
21
22! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
23! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
24! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
25! vary with the Boussinesq approximation, the Boussinesq variant is given first.
26
27!> The control structure for the user_ice_shelf module
28type, public :: user_ice_shelf_cs ; private
29 real :: rho_ocean !< The ocean's typical density [R ~> kg m-3].
30 real :: max_draft !< The maximum ocean draft of the ice shelf [Z ~> m].
31 real :: min_draft !< The minimum ocean draft of the ice shelf [Z ~> m].
32 real :: flat_shelf_width !< The range over which the shelf is min_draft thick [km].
33 real :: shelf_slope_scale !< The range over which the shelf slopes [km].
34 real :: pos_shelf_edge_0 !< The x-position of the shelf edge at time 0 [km].
35 real :: shelf_speed !< The ice shelf speed of translation [km day-1]
36 logical :: first_call = .true. !< If true, this module has not been called before.
37end type user_ice_shelf_cs
38
39contains
40
41!> This subroutine sets up the initial mass and area covered by the ice shelf, based on user-provided code.
42subroutine user_initialize_shelf_mass(mass_shelf, area_shelf_h, h_shelf, hmask, G, US, CS, param_file, new_sim)
43
44 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
45 real, dimension(SZDI_(G),SZDJ_(G)), &
46 intent(out) :: mass_shelf !< The ice shelf mass per unit area averaged
47 !! over the full ocean cell [R Z ~> kg m-2].
48 real, dimension(SZDI_(G),SZDJ_(G)), &
49 intent(out) :: h_shelf !< The ice shelf thickness [Z ~> m].
50 real, dimension(SZDI_(G),SZDJ_(G)), &
51 intent(out) :: area_shelf_h !< The area per cell covered by the ice shelf [L2 ~> m2].
52 real, dimension(SZDI_(G),SZDJ_(G)), &
53 intent(out) :: hmask !< A mask indicating which tracer points are
54 !! partly or fully covered by an ice-shelf
55 type(unit_scale_type), intent(in) :: us !< A structure containing unit conversion factors
56 type(user_ice_shelf_cs), pointer :: cs !< A pointer to the user ice shelf control structure
57 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
58 logical, intent(in) :: new_sim !< If true, this is a new run; otherwise it is
59 !! being started from a restart file.
60
61! This subroutine sets up the initial mass and area covered by the ice shelf.
62 character(len=40) :: mdl = "USER_initialize_shelf_mass" ! This subroutine's name.
63
64 ! call MOM_error(FATAL, "USER_shelf_init.F90, USER_set_shelf_mass: " // &
65 ! "Unmodified user routine called - you must edit the routine to use it")
66
67 if (.not.associated(cs)) allocate(cs)
68
69 ! Read all relevant parameters and write them to the model log.
70 if (cs%first_call) call write_user_log(param_file)
71 cs%first_call = .false.
72 call get_param(param_file, mdl, "RHO_0", cs%Rho_ocean, &
73 "The mean ocean density used with BOUSSINESQ true to "//&
74 "calculate accelerations and the mass for conservation "//&
75 "properties, or with BOUSSINESQ false to convert some "//&
76 "parameters from vertical units of m to kg m-2.", &
77 units="kg m-3", default=1035.0, scale=us%kg_m3_to_R)
78 call get_param(param_file, mdl, "SHELF_MAX_DRAFT", cs%max_draft, &
79 units="m", default=1.0, scale=us%m_to_Z)
80 call get_param(param_file, mdl, "SHELF_MIN_DRAFT", cs%min_draft, &
81 units="m", default=1.0, scale=us%m_to_Z)
82 call get_param(param_file, mdl, "FLAT_SHELF_WIDTH", cs%flat_shelf_width, &
83 units="axis_units", default=0.0)
84 call get_param(param_file, mdl, "SHELF_SLOPE_SCALE", cs%shelf_slope_scale, &
85 units="axis_units", default=0.0)
86 call get_param(param_file, mdl, "SHELF_EDGE_POS_0", cs%pos_shelf_edge_0, &
87 units="axis_units", default=0.0)
88 call get_param(param_file, mdl, "SHELF_SPEED", cs%shelf_speed, &
89 units="axis_units day-1", default=0.0)
90
91 call user_update_shelf_mass(mass_shelf, area_shelf_h, h_shelf, hmask, g, cs, set_time(0,0), new_sim)
92
93end subroutine user_initialize_shelf_mass
94
95!> This subroutine updates the ice shelf thickness, as specified by user-provided code.
96subroutine user_init_ice_thickness(h_shelf, area_shelf_h, hmask, G, US, param_file)
97 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
98 real, dimension(SZDI_(G),SZDJ_(G)), &
99 intent(out) :: h_shelf !< The ice shelf thickness [Z ~> m].
100 real, dimension(SZDI_(G),SZDJ_(G)), &
101 intent(out) :: area_shelf_h !< The area per cell covered by the ice shelf [L2 ~> m2].
102 real, dimension(SZDI_(G),SZDJ_(G)), &
103 intent(out) :: hmask !< A mask indicating which tracer points are
104 !! partly or fully covered by an ice-shelf [nondim]
105 type(unit_scale_type), intent(in) :: us !< A structure containing unit conversion factors
106 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
107
108 ! This subroutine initializes the ice shelf thickness. Currently it does so
109 ! calling USER_initialize_shelf_mass, but this can be revised as needed.
110 real, dimension(SZI_(G),SZJ_(G)) :: mass_shelf ! The ice shelf mass per unit area averaged
111 ! over the full ocean cell [R Z ~> kg m-2].
112 type(user_ice_shelf_cs), pointer :: cs => null()
113
114 call user_initialize_shelf_mass(mass_shelf, area_shelf_h, h_shelf, hmask, g, us, cs, param_file, .true.)
115
116end subroutine user_init_ice_thickness
117
118!> This subroutine updates the ice shelf mass, as specified by user-provided code.
119subroutine user_update_shelf_mass(mass_shelf, area_shelf_h, h_shelf, hmask, G, CS, Time, new_sim)
120 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
121 real, dimension(SZDI_(G),SZDJ_(G)), &
122 intent(inout) :: mass_shelf !< The ice shelf mass per unit area averaged
123 !! over the full ocean cell [R Z ~> kg m-2].
124 real, dimension(SZDI_(G),SZDJ_(G)), &
125 intent(inout) :: area_shelf_h !< The area per cell covered by the ice shelf [L2 ~> m2].
126 real, dimension(SZDI_(G),SZDJ_(G)), &
127 intent(inout) :: h_shelf !< The ice shelf thickness [Z ~> m].
128 real, dimension(SZDI_(G),SZDJ_(G)), &
129 intent(inout) :: hmask !< A mask indicating which tracer points are
130 !! partly or fully covered by an ice-shelf [nondim]
131 type(user_ice_shelf_cs), pointer :: cs !< A pointer to the user ice shelf control structure
132 type(time_type), intent(in) :: time !< The current model time
133 logical, intent(in) :: new_sim !< If true, this the start of a new run.
134
135
136 real :: c1 ! The inverse of the range over which the shelf slopes [km-1]
137 real :: edge_pos ! The time-evolving position the ice shelf edge [km]
138 real :: slope_pos ! The time-evolving position of the start of the ice shelf slope [km]
139 integer :: i, j
140
141 edge_pos = cs%pos_shelf_edge_0 + cs%shelf_speed*(time_type_to_real(time) / 86400.0)
142
143 slope_pos = edge_pos - cs%flat_shelf_width
144 c1 = 0.0 ; if (cs%shelf_slope_scale > 0.0) c1 = 1.0 / cs%shelf_slope_scale
145
146
147 do j=g%jsd,g%jed
148
149 if (((j+g%jdg_offset) <= g%domain%njglobal+g%domain%njhalo) .AND. &
150 ((j+g%jdg_offset) >= g%domain%njhalo+1)) then
151
152 do i=g%isc,g%iec
153
154! if (((i+G%idg_offset) <= G%domain%niglobal+G%domain%nihalo) .AND. &
155! ((i+G%idg_offset) >= G%domain%nihalo+1)) then
156
157 if ((j >= g%jsc) .and. (j <= g%jec)) then
158 if (new_sim) then ; if (g%geoLonCu(i-1,j) >= edge_pos) then
159 ! Everything past the edge is open ocean.
160 mass_shelf(i,j) = 0.0
161 area_shelf_h(i,j) = 0.0
162 hmask(i,j) = 0.0
163 h_shelf(i,j) = 0.0
164 else
165 if (g%geoLonCu(i,j) > edge_pos) then
166 area_shelf_h(i,j) = g%areaT(i,j) * (edge_pos - g%geoLonCu(i-1,j)) / &
167 (g%geoLonCu(i,j) - g%geoLonCu(i-1,j))
168 hmask(i,j) = 2.0
169 else
170 area_shelf_h(i,j) = g%areaT(i,j)
171 hmask(i,j) = 1.0
172 endif
173
174 if (g%geoLonT(i,j) > slope_pos) then
175 h_shelf(i,j) = cs%min_draft
176 mass_shelf(i,j) = cs%Rho_ocean * cs%min_draft
177 else
178 mass_shelf(i,j) = cs%Rho_ocean * (cs%min_draft + &
179 (cs%max_draft - cs%min_draft) * &
180 min(1.0, (c1*(slope_pos - g%geoLonT(i,j)))**2) )
181 h_shelf(i,j) = (cs%min_draft + &
182 (cs%max_draft - cs%min_draft) * &
183 min(1.0, (c1*(slope_pos - g%geoLonT(i,j)))**2) )
184 endif
185 endif ; endif
186 endif
187
188 if ((i+g%idg_offset) == g%domain%nihalo+1) then
189 hmask(i-1,j) = 3.0
190 endif
191
192 enddo
193 endif
194 enddo
195
196end subroutine user_update_shelf_mass
197
198!> This subroutine writes out the user ice shelf code version number to the model log.
199subroutine write_user_log(param_file)
200 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
201
202 character(len=128) :: version = '$Id: user_shelf_init.F90,v 1.1.2.7 2012/06/19 22:15:52 Robert.Hallberg Exp $'
203 character(len=128) :: tagname = '$Name: MOM_ogrp $'
204 character(len=40) :: mdl = "user_shelf_init" ! This module's name.
205
206 call log_version(param_file, mdl, version, tagname)
207
208end subroutine write_user_log
209
210end module user_shelf_init