MOM_ice_shelf_initialize.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!> Initialize ice shelf variables
7
9use mom_array_transform, only : rotate_array
11use mom_file_parser, only : get_param, read_param, log_param, param_file_type
12use mom_io, only: mom_read_data, file_exists, field_exists, slasher, corner
13use mom_error_handler, only : mom_error, mom_mesg, fatal, warning, is_root_pe
16
17implicit none ; private
18
19#include <MOM_memory.h>
20
28! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
29! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
30! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
31! vary with the Boussinesq approximation, the Boussinesq variant is given first.
32
33contains
34
35!> Initialize ice shelf thickness
36subroutine initialize_ice_thickness(h_shelf, area_shelf_h, hmask, melt_mask, G, G_in, US, PF, rotate_index, turns)
37 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
38 type(ocean_grid_type), intent(in) :: g_in !< The ocean's unrotated grid structure
39 real, dimension(SZDI_(G),SZDJ_(G)), &
40 intent(inout) :: h_shelf !< The ice shelf thickness [Z ~> m].
41 real, dimension(SZDI_(G),SZDJ_(G)), &
42 intent(inout) :: area_shelf_h !< The area per cell covered by the ice shelf [L2 ~> m2].
43 real, dimension(SZDI_(G),SZDJ_(G)), &
44 intent(inout) :: hmask !< A mask indicating which tracer points are
45 !! partly or fully covered by an ice-shelf [nondim]
46 real, dimension(SZDI_(G),SZDJ_(G)), &
47 intent(inout) :: melt_mask !< A mask indicating where to allow ice-shelf melting [nondim]
48 type(unit_scale_type), intent(in) :: us !< A structure containing unit conversion factors
49 type(param_file_type), intent(in) :: pf !< A structure to parse for run-time parameters
50 logical, intent(in), optional :: rotate_index !< If true, this is a rotation test
51 integer, intent(in), optional :: turns !< Number of turns for rotation test
52
53 character(len=40) :: mdl = "initialize_ice_thickness" ! This subroutine's name.
54 character(len=200) :: config
55 logical :: rotate = .false.
56 real, allocatable, dimension(:,:) :: tmp1_2d ! Temporary array for storing ice shelf input data [Z~>m]
57 real, allocatable, dimension(:,:) :: tmp2_2d ! Temporary array for storing ice shelf input data [L2~>m2]
58 real, allocatable, dimension(:,:) :: tmp3_2d ! Temporary array for storing ice shelf input data [nondim]
59 real, allocatable, dimension(:,:) :: tmp4_2d ! Temporary array for storing ice shelf input data [nondim]
60
61 call get_param(pf, mdl, "ICE_PROFILE_CONFIG", config, &
62 "This specifies how the initial ice profile is specified. "//&
63 "Valid values are: CHANNEL, FILE, and USER.", &
64 fail_if_missing=.true.)
65
66 if (PRESENT(rotate_index)) rotate=rotate_index
67
68 if (rotate) then
69 allocate(tmp1_2d(g_in%isd:g_in%ied,g_in%jsd:g_in%jed), source=0.0)
70 allocate(tmp2_2d(g_in%isd:g_in%ied,g_in%jsd:g_in%jed), source=0.0)
71 allocate(tmp3_2d(g_in%isd:g_in%ied,g_in%jsd:g_in%jed), source=0.0)
72 allocate(tmp4_2d(g_in%isd:g_in%ied,g_in%jsd:g_in%jed), source=1.0)
73 select case ( trim(config) )
74 case ("CHANNEL") ; call initialize_ice_thickness_channel (tmp1_2d, tmp2_2d, tmp3_2d, g_in, us, pf)
75 case ("FILE") ; call initialize_ice_thickness_from_file (tmp1_2d, tmp2_2d, tmp3_2d, tmp4_2d, g_in, us, pf)
76 case ("USER") ; call user_init_ice_thickness (tmp1_2d, tmp2_2d, tmp3_2d, g_in, us, pf)
77 case default ; call mom_error(fatal,"MOM_initialize: Unrecognized ice profile setup "//trim(config))
78 end select
79 call rotate_array(tmp1_2d,turns, h_shelf)
80 call rotate_array(tmp2_2d,turns, area_shelf_h)
81 call rotate_array(tmp3_2d,turns, hmask)
82 call rotate_array(tmp4_2d,turns, melt_mask)
83 deallocate(tmp1_2d,tmp2_2d,tmp3_2d)
84 else
85 select case ( trim(config) )
86 case ("CHANNEL") ; call initialize_ice_thickness_channel (h_shelf, area_shelf_h, hmask, g, us, pf)
87 case ("FILE") ; call initialize_ice_thickness_from_file (h_shelf, area_shelf_h, hmask, melt_mask, g, us, pf)
88 case ("USER") ; call user_init_ice_thickness (h_shelf, area_shelf_h, hmask, g, us, pf)
89 case default ; call mom_error(fatal,"MOM_initialize: Unrecognized ice profile setup "//trim(config))
90 end select
91 endif
92
93end subroutine initialize_ice_thickness
94
95!> Initialize ice shelf thickness from file
96subroutine initialize_ice_thickness_from_file(h_shelf, area_shelf_h, hmask, melt_mask, G, US, PF)
97 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
98 real, dimension(SZDI_(G),SZDJ_(G)), &
99 intent(inout) :: h_shelf !< The ice shelf thickness [Z ~> m].
100 real, dimension(SZDI_(G),SZDJ_(G)), &
101 intent(inout) :: area_shelf_h !< The area per cell covered by the ice shelf [L2 ~> m2].
102 real, dimension(SZDI_(G),SZDJ_(G)), &
103 intent(inout) :: hmask !< A mask indicating which tracer points are
104 !! partly or fully covered by an ice-shelf [nondim]
105 real, dimension(SZDI_(G),SZDJ_(G)), &
106 intent(inout) :: melt_mask !< A mask indicating where to allow ice-shelf melting [nondim]
107 type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors
108 type(param_file_type), intent(in) :: PF !< A structure to parse for run-time parameters
109
110 ! This subroutine reads ice thickness and area from a file and puts it into
111 ! h_shelf [Z ~> m] and area_shelf_h [L2 ~> m2] (and dimensionless) and updates hmask
112 character(len=200) :: filename,thickness_file,inputdir ! Strings for file/path
113 character(len=200) :: thickness_varname, area_varname, hmask_varname, melt_mask_varname ! Variable name in file
114 character(len=40) :: mdl = "initialize_ice_thickness_from_file" ! This subroutine's name.
115 integer :: i, j, isc, jsc, iec, jec
116 logical :: hmask_set
117 real :: len_sidestress, udh
118
119 call mom_mesg("Initialize_ice_thickness_from_file: reading thickness")
120
121 call get_param(pf, mdl, "INPUTDIR", inputdir, default=".")
122 inputdir = slasher(inputdir)
123 call get_param(pf, mdl, "ICE_THICKNESS_FILE", thickness_file, &
124 "The file from which the bathymetry is read.", &
125 default="ice_shelf_h.nc")
126 call get_param(pf, mdl, "LEN_SIDE_STRESS", len_sidestress, &
127 "position past which shelf sides are stress free.", &
128 default=0.0, units="axis_units")
129
130 filename = trim(inputdir)//trim(thickness_file)
131 call log_param(pf, mdl, "INPUTDIR/THICKNESS_FILE", filename)
132 call get_param(pf, mdl, "ICE_THICKNESS_VARNAME", thickness_varname, &
133 "The name of the thickness variable in ICE_THICKNESS_FILE.", &
134 default="h_shelf")
135 call get_param(pf, mdl, "ICE_AREA_VARNAME", area_varname, &
136 "The name of the area variable in ICE_THICKNESS_FILE.", &
137 default="area_shelf_h")
138 hmask_varname="h_mask"
139 call get_param(pf, mdl, "MELT_MASK_VARNAME", melt_mask_varname, &
140 "The name of the melt mask variable in ICE_THICKNESS_FILE.", &
141 default="melt_mask")
142 if (.not.file_exists(filename, g%Domain)) call mom_error(fatal, &
143 " initialize_topography_from_file: Unable to open "//trim(filename))
144 call mom_read_data(filename, trim(thickness_varname), h_shelf, g%Domain, scale=us%m_to_Z)
145 call mom_read_data(filename,trim(area_varname), area_shelf_h, g%Domain, scale=us%m_to_L**2)
146 if (field_exists(filename, trim(hmask_varname), mom_domain=g%Domain)) then
147 call mom_read_data(filename, trim(hmask_varname), hmask, g%Domain)
148 hmask_set = .true.
149 else
150 call mom_error(warning, "Ice shelf thickness initialized without setting the shelf mask "//&
151 "from variable "//trim(hmask_varname)//", which does not exist in "//trim(filename))
152 hmask_set = .false.
153 endif
154 if (field_exists(filename, trim(melt_mask_varname), mom_domain=g%Domain)) then
155 call mom_read_data(filename, trim(melt_mask_varname), melt_mask, g%Domain)
156 else
157 melt_mask(:,:)=1.0
158 endif
159
160 isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
161
162 if (.not.hmask_set) then
163 ! Set hmask based on the values in h_shelf.
164 do j=jsc,jec ; do i=isc,iec
165 hmask(i,j) = 0.0
166 if (h_shelf(i,j) > 0.0) hmask(i,j) = 1.0
167 enddo ; enddo
168 endif
169
170 do j=jsc,jec
171 do i=isc,iec
172
173 ! taper ice shelf in area where there is no sidestress -
174 ! but do not interfere with hmask
175
176 if ((len_sidestress > 0.) .and. (g%geoLonCv(i,j) > len_sidestress)) then
177 udh = exp(-(g%geoLonCv(i,j)-len_sidestress)/5.0) * h_shelf(i,j)
178 if (udh <= 25.0) then
179 h_shelf(i,j) = 0.0
180 area_shelf_h(i,j) = 0.0
181 else
182 h_shelf(i,j) = udh
183 endif
184 endif
185
186 ! update thickness mask
187
188 if (area_shelf_h(i,j) >= g%areaT(i,j)) then
189 hmask(i,j) = 1.
190 area_shelf_h(i,j)=g%areaT(i,j)
191 elseif (area_shelf_h(i,j) == 0.0) then
192 hmask(i,j) = 0.
193 elseif ((area_shelf_h(i,j) > 0) .and. (area_shelf_h(i,j) <= g%areaT(i,j))) then
194 hmask(i,j) = 2.
195 else
196 call mom_error(fatal,mdl// " AREA IN CELL OUT OF RANGE")
197 endif
198 enddo
199 enddo
201
202!> Initialize ice shelf thickness for a channel configuration
203subroutine initialize_ice_thickness_channel(h_shelf, area_shelf_h, hmask, G, US, PF)
204 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
205 real, dimension(SZDI_(G),SZDJ_(G)), &
206 intent(inout) :: h_shelf !< The ice shelf thickness [Z ~> m].
207 real, dimension(SZDI_(G),SZDJ_(G)), &
208 intent(inout) :: area_shelf_h !< The area per cell covered by the ice shelf [L2 ~> m2].
209 real, dimension(SZDI_(G),SZDJ_(G)), &
210 intent(inout) :: hmask !< A mask indicating which tracer points are
211 !! partly or fully covered by an ice-shelf
212 type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors
213 type(param_file_type), intent(in) :: PF !< A structure to parse for run-time parameters
214
215 character(len=40) :: mdl = "initialize_ice_shelf_thickness_channel" ! This subroutine's name.
216 real :: max_draft, min_draft, flat_shelf_width, c1, slope_pos
217 real :: edge_pos, shelf_slope_scale
218 integer :: i, j, jsc, jec, jsd, jed, jedg, nyh, isc, iec, isd, ied
219 integer :: j_off
220
221 jsc = g%jsc ; jec = g%jec ; isc = g%isc ; iec = g%iec
222 jsd = g%jsd ; jed = g%jed ; isd = g%isd ; ied = g%ied
223 nyh = g%domain%njhalo ; jedg = g%domain%njglobal+nyh
224 j_off = g%jdg_offset
225
226 call mom_mesg(mdl//": setting thickness")
227
228 call get_param(pf, mdl, "SHELF_MAX_DRAFT", max_draft, &
229 units="m", default=1.0, scale=us%m_to_Z)
230 call get_param(pf, mdl, "SHELF_MIN_DRAFT", min_draft, &
231 units="m", default=1.0, scale=us%m_to_Z)
232 call get_param(pf, mdl, "FLAT_SHELF_WIDTH", flat_shelf_width, &
233 units="axis_units", default=0.0)
234 call get_param(pf, mdl, "SHELF_SLOPE_SCALE", shelf_slope_scale, &
235 units="axis_units", default=0.0)
236 call get_param(pf, mdl, "SHELF_EDGE_POS_0", edge_pos, &
237 units="axis_units", default=0.0)
238! call get_param(param_file, mdl, "RHO_0", Rho_ocean, &
239! "The mean ocean density used with BOUSSINESQ true to "//&
240! "calculate accelerations and the mass for conservation "//&
241! "properties, or with BOUSSINESQ false to convert some "//&
242! "parameters from vertical units of m to kg m-2.", &
243! units="kg m-3", default=1035.0, scale=US%Z_to_m)
244
245 slope_pos = edge_pos - flat_shelf_width
246 c1 = 0.0 ; if (shelf_slope_scale > 0.0) c1 = 1.0 / shelf_slope_scale
247
248
249 do j=g%jsd,g%jed
250
251 if (((j+j_off) <= jedg) .AND. ((j+j_off) >= nyh+1)) then
252
253 do i=g%isc,g%iec
254
255 if ((j >= jsc) .and. (j <= jec)) then
256
257 if (g%geoLonCu(i-1,j) >= edge_pos) then
258 ! Everything past the edge is open ocean.
259 area_shelf_h(i,j) = 0.0
260 hmask(i,j) = 0.0
261 h_shelf(i,j) = 0.0
262 else
263 if (g%geoLonCu(i,j) > edge_pos) then
264 area_shelf_h(i,j) = g%areaT(i,j) * (edge_pos - g%geoLonCu(i-1,j)) / &
265 (g%geoLonCu(i,j) - g%geoLonCu(i-1,j))
266 hmask(i,j) = 2.0
267 else
268 area_shelf_h(i,j) = g%areaT(i,j)
269 hmask(i,j) = 1.0
270 endif
271
272 if (g%geoLonT(i,j) > slope_pos) then
273 h_shelf(i,j) = min_draft
274 else
275 h_shelf(i,j) = (min_draft + &
276 (max_draft - min_draft) * &
277 min(1.0, (c1*(slope_pos - g%geoLonT(i,j)))**2) )
278 endif
279
280 endif
281 endif
282
283 if ((i+g%idg_offset) == g%domain%nihalo+1) then
284 hmask(i-1,j) = 3.0
285 endif
286
287 enddo
288 endif ; enddo
289
291
292!> Initialize ice shelf boundary conditions for a channel configuration
293subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_bdry, &
294 u_flux_bdry_val, v_flux_bdry_val, u_bdry_val, v_bdry_val, u_shelf, v_shelf, h_bdry_val, &
295 hmask, h_shelf, G, US, PF )
296
297 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
298 real, dimension(SZIB_(G),SZJB_(G)), &
299 intent(inout) :: u_face_mask_bdry !< A boundary-type mask at C-grid u faces
300
301 real, dimension(SZIB_(G),SZJ_(G)), &
302 intent(inout) :: u_flux_bdry_val !< The boundary thickness flux through
303 !! C-grid u faces [L Z T-1 ~> m2 s-1].
304 real, dimension(SZIB_(G),SZJB_(G)), &
305 intent(inout) :: v_face_mask_bdry !< A boundary-type mask at C-grid v faces
306
307 real, dimension(SZI_(G),SZJB_(G)), &
308 intent(inout) :: v_flux_bdry_val !< The boundary thickness flux through
309 !! C-grid v faces [L Z T-1 ~> m2 s-1].
310 real, dimension(SZIB_(G),SZJB_(G)), &
311 intent(inout) :: u_bdry_val !< The zonal ice shelf velocity at open
312 !! boundary vertices [L T-1 ~> m s-1].
313 real, dimension(SZIB_(G),SZJB_(G)), &
314 intent(inout) :: v_bdry_val !< The meridional ice shelf velocity at open
315 real, dimension(SZIB_(G),SZJB_(G)), &
316 intent(inout) :: u_shelf !< The zonal ice shelf velocity [L T-1 ~> m s-1].
317 real, dimension(SZIB_(G),SZJB_(G)), &
318 intent(inout) :: v_shelf !< The meridional ice shelf velocity [L T-1 ~> m s-1].
319 real, dimension(SZDI_(G),SZDJ_(G)), &
320 intent(inout) :: h_bdry_val !< The ice shelf thickness at open boundaries [Z ~> m]
321 real, dimension(SZDI_(G),SZDJ_(G)), &
322 intent(inout) :: hmask !< A mask indicating which tracer points are
323 !! partly or fully covered by an ice-shelf
324 real, dimension(SZDI_(G),SZDJ_(G)), &
325 intent(inout) :: h_shelf !< Ice-shelf thickness [Z ~> m]
326 type(unit_scale_type), intent(in) :: us !< A structure containing unit conversion factors
327 type(param_file_type), intent(in) :: pf !< A structure to parse for run-time parameters
328
329 character(len=40) :: mdl = "initialize_ice_shelf_boundary_channel" ! This subroutine's name.
330 integer :: i, j, isd, jsd, giec, gjec, gisc, gjsc,gisd,gjsd, isc, jsc, iec, jec, ied, jed
331 real :: input_thick ! The input ice shelf thickness [Z ~> m]
332 real :: input_vel ! The input ice velocity at the upstream boundary [L T-1 ~> m s-1]
333 real :: lenlat, len_stress, westlon, lenlon, southlat ! The input positions of the channel boundarises
334
335 lenlat = g%len_lat
336 lenlon = g%len_lon
337 westlon = g%west_lon
338 southlat = g%south_lat
339
340 call get_param(pf, mdl, "INPUT_VEL_ICE_SHELF", input_vel, &
341 "inflow ice velocity at upstream boundary", &
342 units="m s-1", default=0., scale=us%m_s_to_L_T)
343 call get_param(pf, mdl, "INPUT_THICK_ICE_SHELF", input_thick, &
344 "flux thickness at upstream boundary", &
345 units="m", default=1000., scale=us%m_to_Z)
346 call get_param(pf, mdl, "LEN_SIDE_STRESS", len_stress, &
347 "maximum position of no-flow condition in along-flow direction", &
348 units="km", default=0.)
349
350 call mom_mesg(mdl//": setting boundary")
351
352 isd = g%isd ; ied = g%ied
353 jsd = g%jsd ; jed = g%jed
354 isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
355 gjsd = g%Domain%njglobal ; gisd = g%Domain%niglobal
356 gisc = g%Domain%nihalo ; gjsc = g%Domain%njhalo
357 giec = g%Domain%niglobal+gisc ; gjec = g%Domain%njglobal+gjsc
358
359 !---------b.c.s based on geopositions -----------------
360 do j=jsc,jec+1
361 do i=isc-1,iec+1
362 ! upstream boundary - set either dirichlet or flux condition
363
364 if (g%geoLonBu(i,j) == westlon) then
365 hmask(i+1,j) = 3.0
366 !---
367 !OLD: thickness_bdry_val was used for ice dynamics, and h_bdry_val was not used anywhere except here:
368 !h_bdry_val(i+1,j) = h_shelf(i+1,j) ; thickness_bdry_val(i+1,j) = h_bdry_val(i+0*1,j)
369 !---
370 !NEW: h_bdry_val is used for ice dynamics instead of thickness_bdry_val, which was removed
371 h_bdry_val(i+1,j) = h_shelf(i+0*1,j) !why 0*1
372 !---
373 u_face_mask_bdry(i+1,j) = 5.0
374 u_bdry_val(i+1,j) = input_vel*(1-16.0*((g%geoLatBu(i-1,j)/lenlat-0.5))**4) !velocity distribution
375 endif
376
377
378 ! side boundaries: no flow
379 if (g%geoLatBu(i,j-1) == southlat) then !bot boundary
380 if (len_stress == 0. .OR. g%geoLonCv(i,j) <= len_stress) then
381 v_face_mask_bdry(i,j+1) = 0.
382 u_face_mask_bdry(i,j) = 3.
383 u_bdry_val(i,j) = 0.
384 v_bdry_val(i,j) = 0.
385 else
386 v_face_mask_bdry(i,j+1) = 1.
387 u_face_mask_bdry(i,j) = 3.
388 u_bdry_val(i,j) = 0.
389 v_bdry_val(i,j) = 0.
390 endif
391 elseif (g%geoLatBu(i,j-1) == southlat+lenlat) then !top boundary
392 if (len_stress == 0. .OR. g%geoLonCv(i,j) <= len_stress) then
393 v_face_mask_bdry(i,j-1) = 0.
394 u_face_mask_bdry(i,j-1) = 3.
395 else
396 v_face_mask_bdry(i,j-1) = 3.
397 u_face_mask_bdry(i,j-1) = 3.
398 endif
399 endif
400
401 ! downstream boundary - CFBC
402 if (g%geoLonBu(i,j) == westlon+lenlon) then
403 u_face_mask_bdry(i-1,j) = 2.0
404 endif
405
406 enddo
407 enddo
409
410
411!> Initialize ice shelf flow from file
412subroutine initialize_ice_flow_from_file(bed_elev,u_shelf, v_shelf,float_cond,&
413 G, US, PF)
414 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
415 real, dimension(SZDI_(G),SZDJ_(G)), &
416 intent(inout) :: bed_elev !< The bed elevation [Z ~> m].
417 real, dimension(SZIB_(G),SZJB_(G)), &
418 intent(inout) :: u_shelf !< The zonal ice shelf velocity [L T-1 ~> m s-1].
419 real, dimension(SZIB_(G),SZJB_(G)), &
420 intent(inout) :: v_shelf !< The meridional ice shelf velocity [L T-1 ~> m s-1].
421 real, dimension(SZDI_(G),SZDJ_(G)), &
422 intent(inout) :: float_cond !< An array indicating where the ice
423 !! shelf is floating: 0 if floating, 1 if not. [nondim]
424 type(unit_scale_type), intent(in) :: us !< A structure containing unit conversion factors
425 type(param_file_type), intent(in) :: pf !< A structure to parse for run-time parameters
426
427 ! This subroutine reads ice thickness and area from a file and puts it into
428 ! h_shelf [Z ~> m] and area_shelf_h [L2 ~> m2] (and dimensionless) and updates hmask
429 character(len=200) :: filename,vel_file,inputdir,bed_topo_file ! Strings for file/path
430 character(len=200) :: ushelf_varname, vshelf_varname, &
431 floatfr_varname, bed_varname ! Variable name in file
432 character(len=40) :: mdl = "initialize_ice_velocity_from_file" ! This subroutine's name.
433
434 call mom_mesg(" MOM_ice_shelf_init_profile.F90, initialize_velocity_from_file: reading velocity")
435
436 call get_param(pf, mdl, "INPUTDIR", inputdir, default=".")
437 inputdir = slasher(inputdir)
438 call get_param(pf, mdl, "ICE_VELOCITY_FILE", vel_file, &
439 "The file from which the velocity is read.", &
440 default="ice_shelf_vel.nc")
441
442 filename = trim(inputdir)//trim(vel_file)
443 call log_param(pf, mdl, "INPUTDIR/THICKNESS_FILE", filename)
444 call get_param(pf, mdl, "ICE_U_VEL_VARNAME", ushelf_varname, &
445 "The name of the u velocity variable in ICE_VELOCITY_FILE.", &
446 default="u_shelf")
447 call get_param(pf, mdl, "ICE_V_VEL_VARNAME", vshelf_varname, &
448 "The name of the v velocity variable in ICE_VELOCITY_FILE.", &
449 default="v_shelf")
450 call get_param(pf, mdl, "ICE_FLOAT_FRAC_VARNAME", floatfr_varname, &
451 "The name of the ice float fraction (grounding fraction) variable in ICE_VELOCITY_FILE.", &
452 default="float_frac")
453 call get_param(pf, mdl, "BED_TOPO_FILE", bed_topo_file, &
454 "The file from which the bed elevation is read.", &
455 default="ice_shelf_vel.nc")
456 call get_param(pf, mdl, "BED_TOPO_VARNAME", bed_varname, &
457 "The name of the bed elevation variable in ICE_INPUT_FILE.", &
458 default="depth")
459 if (.not.file_exists(filename, g%Domain)) call mom_error(fatal, &
460 " initialize_ice_shelf_velocity_from_file: Unable to open "//trim(filename))
461
462 call mom_read_data(filename, trim(ushelf_varname), u_shelf, g%Domain, position=corner, scale=us%m_s_to_L_T)
463 call mom_read_data(filename, trim(vshelf_varname), v_shelf, g%Domain, position=corner, scale=us%m_s_to_L_T)
464 call mom_read_data(filename, trim(floatfr_varname), float_cond, g%Domain, scale=1.)
465
466 filename = trim(inputdir)//trim(bed_topo_file)
467 call mom_read_data(filename, trim(bed_varname), bed_elev, g%Domain, scale=us%m_to_Z)
468
469
471
472!> Initialize ice shelf b.c.s from file
473subroutine initialize_ice_shelf_boundary_from_file(u_face_mask_bdry, v_face_mask_bdry, &
474 u_bdry_val, v_bdry_val, umask, vmask, h_bdry_val, &
475 hmask, h_shelf, G, US, PF )
476
477 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
478 real, dimension(SZIB_(G),SZJB_(G)), &
479 intent(inout) :: u_face_mask_bdry !< A boundary-type mask at B-grid u faces [nondim]
480 real, dimension(SZIB_(G),SZJB_(G)), &
481 intent(inout) :: v_face_mask_bdry !< A boundary-type mask at B-grid v faces [nondim]
482 real, dimension(SZIB_(G),SZJB_(G)), &
483 intent(inout) :: u_bdry_val !< The zonal ice shelf velocity at open
484 !! boundary vertices [L T-1 ~> m s-1].
485 real, dimension(SZIB_(G),SZJB_(G)), &
486 intent(inout) :: v_bdry_val !< The meridional ice shelf velocity at open
487 !! boundary vertices [L T-1 ~> m s-1].
488 real, dimension(SZDIB_(G),SZDJB_(G)), &
489 intent(inout) :: umask !< A mask for ice shelf velocity [nondim]
490 real, dimension(SZDIB_(G),SZDJB_(G)), &
491 intent(inout) :: vmask !< A mask for ice shelf velocity [nondim]
492 real, dimension(SZDI_(G),SZDJ_(G)), &
493 intent(inout) :: h_bdry_val !< The ice shelf thickness at open boundaries [Z ~> m]
494 real, dimension(SZDI_(G),SZDJ_(G)), &
495 intent(inout) :: hmask !< A mask indicating which tracer points are
496 !! partly or fully covered by an ice-shelf [nondim]
497 real, dimension(SZDI_(G),SZDJ_(G)), &
498 intent(in) :: h_shelf !< Ice-shelf thickness [Z ~> m]
499 type(unit_scale_type), intent(in) :: us !< A structure containing unit conversion factors
500 type(param_file_type), intent(in) :: pf !< A structure to parse for run-time parameters
501
502 character(len=200) :: filename, bc_file, inputdir, icethick_file ! Strings for file/path
503 character(len=200) :: ufcmskbdry_varname, vfcmskbdry_varname, &
504 ubdryv_varname, vbdryv_varname, umask_varname, vmask_varname, &
505 hmsk_varname ! Variable name in file
506 character(len=40) :: mdl = "initialize_ice_shelf_boundary_from_file" ! This subroutine's name.
507
508 integer :: i, j, isc, jsc, iec, jec
509
510 h_bdry_val(:,:) = 0.
511
512 call mom_mesg(" MOM_ice_shelf_init_profile.F90, initialize_b_c_s_from_file: reading b.c.s")
513
514 call get_param(pf, mdl, "INPUTDIR", inputdir, default=".")
515 inputdir = slasher(inputdir)
516 call get_param(pf, mdl, "ICE_SHELF_BC_FILE", bc_file, &
517 "The file from which the boundary conditions are read.", &
518 default="ice_shelf_bc.nc")
519 call get_param(pf, mdl, "ICE_THICKNESS_FILE", icethick_file, &
520 "The file from which the ice-shelf thickness is read.", &
521 default="ice_shelf_thick.nc")
522 call get_param(pf, mdl, "ICE_THICKNESS_MASK_VARNAME", hmsk_varname, &
523 "The name of the icethickness mask variable in ICE_THICKNESS_FILE.", &
524 default="h_mask")
525
526 filename = trim(inputdir)//trim(bc_file)
527 call log_param(pf, mdl, "INPUTDIR/ICE_SHELF_BC_FILE", filename)
528 call get_param(pf, mdl, "ICE_UBDRYMSK_VARNAME", ufcmskbdry_varname, &
529 "The name of the ice-shelf ubdrymask variable in ICE_SHELF_BC_FILE.", &
530 default="ufacemask")
531 call get_param(pf, mdl, "ICE_VBDRYMSK_VARNAME", vfcmskbdry_varname, &
532 "The name of the ice-shelf vbdrymask variable in ICE_SHELF_BC_FILE.", &
533 default="vfacemask")
534 call get_param(pf, mdl, "ICE_UMASK_VARNAME", umask_varname, &
535 "The name of the ice-shelf ubdrymask variable in ICE_SHELF_BC_FILE.", &
536 default="umask")
537 call get_param(pf, mdl, "ICE_VMASK_VARNAME", vmask_varname, &
538 "The name of the ice-shelf vbdrymask variable in ICE_SHELF_BC_FILE.", &
539 default="vmask")
540 call get_param(pf, mdl, "ICE_UBDRYVAL_VARNAME", ubdryv_varname, &
541 "The name of the ice-shelf ice_shelf ubdry variable in ICE_SHELF_BC_FILE.", &
542 default="ubdry_val")
543 call get_param(pf, mdl, "ICE_VBDRYVAL_VARNAME", vbdryv_varname, &
544 "The name of the ice-shelf ice_shelf vbdry variable in ICE_SHELF_BC_FILE.", &
545 default="vbdry_val")
546 if (.not.file_exists(filename, g%Domain)) call mom_error(fatal, &
547 " initialize_ice_shelf_velocity_from_file: Unable to open "//trim(filename))
548
549
550 call mom_read_data(filename, trim(ufcmskbdry_varname), u_face_mask_bdry, g%Domain, position=corner, &
551 scale=1.)
552 call mom_read_data(filename, trim(vfcmskbdry_varname), v_face_mask_bdry, g%Domain, position=corner, &
553 scale=1.)
554 call mom_read_data(filename, trim(ubdryv_varname), u_bdry_val, g%Domain, position=corner, scale=us%m_s_to_L_T)
555 call mom_read_data(filename, trim(vbdryv_varname), v_bdry_val, g%Domain, position=corner, scale=us%m_s_to_L_T)
556 call mom_read_data(filename, trim(umask_varname), umask, g%Domain, position=corner, scale=1.)
557 call mom_read_data(filename, trim(vmask_varname), vmask, g%Domain, position=corner, scale=1.)
558 filename = trim(inputdir)//trim(icethick_file)
559
560 call mom_read_data(filename,trim(hmsk_varname), hmask, g%Domain, scale=1.)
561 isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
562
563 do j=jsc,jec
564 do i=isc,iec
565 if (hmask(i,j) == 3.) then
566 h_bdry_val(i,j) = h_shelf(i,j)
567 endif
568 enddo
569 enddo
570
572
573!> Initialize ice basal friction
574subroutine initialize_ice_c_basal_friction(C_basal_friction, G, US, PF)
575 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
576 real, dimension(SZDI_(G),SZDJ_(G)), &
577 intent(inout) :: c_basal_friction !< Ice-stream basal friction
578 !! in units of [R L Z T-2 (s m-1)^n_basal_fric ~> Pa (s m-1)^n_basal_fric]
579 type(unit_scale_type), intent(in) :: us !< A structure containing unit conversion factors
580 type(param_file_type), intent(in) :: pf !< A structure to parse for run-time parameters
581
582! integer :: i, j
583 real :: c_friction ! Constant ice-stream basal friction in units of
584 ! [R L Z T-2 (s m-1)^n_basal_fric ~> Pa (s m-1)^n_basal_fric]
585 character(len=40) :: mdl = "initialize_ice_basal_friction" ! This subroutine's name.
586 character(len=200) :: config
587 character(len=200) :: varname
588 character(len=200) :: inputdir, filename, c_friction_file
589
590 call get_param(pf, mdl, "ICE_BASAL_FRICTION_CONFIG", config, &
591 "This specifies how the initial basal friction profile is specified. "//&
592 "Valid values are: CONSTANT and FILE.", &
593 fail_if_missing=.true.)
594
595 if (trim(config)=="CONSTANT") then
596 call get_param(pf, mdl, "BASAL_FRICTION_COEFF", c_friction, &
597 "Coefficient in sliding law.", units="Pa (s m-1)^(n_basal_fric)", default=5.e10, scale=us%Pa_to_RLZ_T2)
598
599 c_basal_friction(:,:) = c_friction
600 elseif (trim(config)=="FILE") then
601 call mom_mesg(" MOM_ice_shelf.F90, initialize_ice_shelf: reading friction coefficients")
602 call get_param(pf, mdl, "INPUTDIR", inputdir, default=".")
603 inputdir = slasher(inputdir)
604
605 call get_param(pf, mdl, "BASAL_FRICTION_FILE", c_friction_file, &
606 "The file from which basal friction coefficients are read.", &
607 default="ice_basal_friction.nc")
608 filename = trim(inputdir)//trim(c_friction_file)
609 call log_param(pf, mdl, "INPUTDIR/BASAL_FRICTION_FILE", filename)
610
611 call get_param(pf, mdl, "BASAL_FRICTION_VARNAME", varname, &
612 "The variable to use in basal traction.", &
613 default="tau_b_beta")
614
615 if (.not.file_exists(filename, g%Domain)) call mom_error(fatal, &
616 " initialize_ice_basal_friction_from_file: Unable to open "//trim(filename))
617
618 call mom_read_data(filename, trim(varname), c_basal_friction, g%Domain, scale=us%Pa_to_RLZ_T2)
619
620 endif
621end subroutine
622
623
624!> Initialize ice-stiffness parameter
625subroutine initialize_ice_aglen(AGlen, ice_viscosity_compute, G, US, PF)
626 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
627 real, dimension(SZDI_(G),SZDJ_(G)), &
628 intent(inout) :: aglen !< The ice-stiffness parameter A_Glen, often in [Pa-3 s-1]
629 character(len=40) :: ice_viscosity_compute !< Specifies whether the ice viscosity is computed internally
630 !! according to Glen's flow law; is constant (for debugging purposes)
631 !! or using observed strain rates and read from a file
632 type(unit_scale_type), intent(in) :: us !< A structure containing unit conversion factors
633 type(param_file_type), intent(in) :: pf !< A structure to parse for run-time parameters
634
635 real :: a_glen ! Ice-stiffness parameter, often in [Pa-3 s-1]
636 character(len=40) :: mdl = "initialize_ice_stiffness" ! This subroutine's name.
637 character(len=200) :: config
638 character(len=200) :: varname
639 character(len=200) :: inputdir, filename, aglen_file
640
641 call get_param(pf, mdl, "ICE_A_GLEN_CONFIG", config, &
642 "This specifies how the initial ice-stiffness parameter is specified. "//&
643 "Valid values are: CONSTANT and FILE.", &
644 fail_if_missing=.true.)
645
646 if (trim(config)=="CONSTANT") then
647 call get_param(pf, mdl, "A_GLEN", a_glen, &
648 "Ice-stiffness parameter.", units="Pa-n_g s-1", default=2.261e-25)
649
650 aglen(:,:) = a_glen
651
652 elseif (trim(config)=="FILE") then
653 call mom_mesg(" MOM_ice_shelf.F90, initialize_ice_shelf: reading ice-stiffness parameter")
654 call get_param(pf, mdl, "INPUTDIR", inputdir, default=".")
655 inputdir = slasher(inputdir)
656
657 call get_param(pf, mdl, "ICE_STIFFNESS_FILE", aglen_file, &
658 "The file from which the ice-stiffness is read.", &
659 default="ice_AGlen.nc")
660 filename = trim(inputdir)//trim(aglen_file)
661 call log_param(pf, mdl, "INPUTDIR/ICE_STIFFNESS_FILE", filename)
662 call get_param(pf, mdl, "A_GLEN_VARNAME", varname, &
663 "The variable to use as ice-stiffness.", &
664 default="A_GLEN")
665
666 if (.not.file_exists(filename, g%Domain)) call mom_error(fatal, &
667 " initialize_ice_stiffness_from_file: Unable to open "//trim(filename))
668
669 if (trim(ice_viscosity_compute) == "OBS") then
670 ! AGlen is the ice viscosity [R L2 T-1 ~> Pa s] computed from obs and read from a file
671 call mom_read_data(filename, trim(varname), aglen, g%Domain, scale=us%Pa_to_RL2_T2*us%s_to_T)
672 else
673 ! AGlen is the ice stiffness parameter [Pa-n_g s-1]
674 call mom_read_data(filename, trim(varname), aglen, g%Domain)
675 endif
676 endif
677end subroutine initialize_ice_aglen
678
679!> Initialize ice surface mass balance field that is held constant over time
680subroutine initialize_ice_smb(SMB, G, US, PF)
681 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
682 real, dimension(SZDI_(G),SZDJ_(G)), &
683 intent(inout) :: smb !< Ice surface mass balance parameter, often in [R Z T-1 ~> kg m-2 s-1]
684 type(unit_scale_type), intent(in) :: us !< A structure containing unit conversion factors
685 type(param_file_type), intent(in) :: pf !< A structure to parse for run-time parameters
686
687 real :: smb_val ! Constant ice surface mass balance parameter, often in [R Z T-1 ~> kg m-2 s-1]
688 character(len=40) :: mdl = "initialize_ice_SMB" ! This subroutine's name.
689 character(len=200) :: config
690 character(len=200) :: varname
691 character(len=200) :: inputdir, filename, smb_file
692
693 call get_param(pf, mdl, "ICE_SMB_CONFIG", config, &
694 "This specifies how the initial ice surface mass balance parameter is specified. "//&
695 "Valid values are: CONSTANT and FILE.", &
696 default="CONSTANT")
697
698 if (trim(config)=="CONSTANT") then
699 call get_param(pf, mdl, "SMB", smb_val, &
700 "Surface mass balance.", units="kg m-2 s-1", default=0.0, scale=us%kg_m2s_to_RZ_T)
701
702 smb(:,:) = smb_val
703
704 elseif (trim(config)=="FILE") then
705 call mom_mesg(" MOM_ice_shelf.F90, initialize_ice_shelf: reading SMB parameter")
706 call get_param(pf, mdl, "INPUTDIR", inputdir, default=".")
707 inputdir = slasher(inputdir)
708
709 call get_param(pf, mdl, "ICE_SMB_FILE", smb_file, &
710 "The file from which the ice surface mass balance is read.", &
711 default="ice_SMB.nc")
712 filename = trim(inputdir)//trim(smb_file)
713 call log_param(pf, mdl, "INPUTDIR/ICE_SMB_FILE", filename)
714 call get_param(pf, mdl, "ICE_SMB_VARNAME", varname, &
715 "The variable to use as surface mass balance.", &
716 default="SMB")
717
718 if (.not.file_exists(filename, g%Domain)) call mom_error(fatal, &
719 " initialize_ice_SMV_from_file: Unable to open "//trim(filename))
720 call mom_read_data(filename,trim(varname), smb, g%Domain, scale=us%kg_m2s_to_RZ_T)
721
722 endif
723end subroutine initialize_ice_smb