13 implicit none ;
private
15 #include <MOM_memory.h>
30 real,
dimension(SZDI_(G),SZDJ_(G)), &
31 intent(inout) :: h_shelf
32 real,
dimension(SZDI_(G),SZDJ_(G)), &
33 intent(inout) :: area_shelf_h
34 real,
dimension(SZDI_(G),SZDJ_(G)), &
35 intent(inout) :: hmask
41 character(len=40) :: mdl =
"initialize_ice_thickness"
42 character(len=200) :: config
44 call get_param(pf, mdl,
"ICE_PROFILE_CONFIG", config, &
45 "This specifies how the initial ice profile is specified. "//&
46 "Valid values are: CHANNEL, FILE, and USER.", &
47 fail_if_missing=.true.)
49 select case ( trim(config) )
53 case default ;
call mom_error(fatal,
"MOM_initialize: "// &
54 "Unrecognized ice profile setup "//trim(config))
62 real,
dimension(SZDI_(G),SZDJ_(G)), &
63 intent(inout) :: h_shelf
64 real,
dimension(SZDI_(G),SZDJ_(G)), &
65 intent(inout) :: area_shelf_h
66 real,
dimension(SZDI_(G),SZDJ_(G)), &
67 intent(inout) :: hmask
74 character(len=200) :: filename,thickness_file,inputdir
75 character(len=200) :: thickness_varname, area_varname
76 character(len=40) :: mdl =
"initialize_ice_thickness_from_file"
77 integer :: i, j, isc, jsc, iec, jec
78 real :: len_sidestress, mask, udh
80 call mom_mesg(
" MOM_ice_shelf_init_profile.F90, initialize_thickness_from_file: reading thickness")
82 call get_param(pf, mdl,
"INPUTDIR", inputdir, default=
".")
83 inputdir = slasher(inputdir)
84 call get_param(pf, mdl,
"ICE_THICKNESS_FILE", thickness_file, &
85 "The file from which the bathymetry is read.", &
86 default=
"ice_shelf_h.nc")
87 call get_param(pf, mdl,
"LEN_SIDE_STRESS", len_sidestress, &
88 "position past which shelf sides are stress free.", &
89 default=0.0, units=
"axis_units")
91 filename = trim(inputdir)//trim(thickness_file)
92 call log_param(pf, mdl,
"INPUTDIR/THICKNESS_FILE", filename)
93 call get_param(pf, mdl,
"ICE_THICKNESS_VARNAME", thickness_varname, &
94 "The name of the thickness variable in ICE_THICKNESS_FILE.", &
96 call get_param(pf, mdl,
"ICE_AREA_VARNAME", area_varname, &
97 "The name of the area variable in ICE_THICKNESS_FILE.", &
98 default=
"area_shelf_h")
101 " initialize_topography_from_file: Unable to open "//trim(filename))
103 call mom_read_data(filename, trim(thickness_varname), h_shelf, g%Domain, scale=us%m_to_Z)
104 call mom_read_data(filename,trim(area_varname),area_shelf_h,g%Domain)
110 isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
118 if ((g%geoLonCv(i,j) > len_sidestress).and. &
119 (len_sidestress > 0.))
then
120 udh = exp(-(g%geoLonCv(i,j)-len_sidestress)/5.0) * h_shelf(i,j)
121 if (udh <= 25.0)
then
123 area_shelf_h(i,j) = 0.0
131 if (area_shelf_h(i,j) >= us%L_to_m**2*g%areaT(i,j))
then
133 elseif (area_shelf_h(i,j) == 0.0)
then
135 elseif ((area_shelf_h(i,j) > 0) .and. (area_shelf_h(i,j) <= us%L_to_m**2*g%areaT(i,j)))
then
138 call mom_error(fatal,mdl//
" AREA IN CELL OUT OF RANGE")
149 real,
dimension(SZDI_(G),SZDJ_(G)), &
150 intent(inout) :: h_shelf
151 real,
dimension(SZDI_(G),SZDJ_(G)), &
152 intent(inout) :: area_shelf_h
153 real,
dimension(SZDI_(G),SZDJ_(G)), &
154 intent(inout) :: hmask
159 character(len=40) :: mdl =
"initialize_ice_shelf_thickness_channel"
160 real :: max_draft, min_draft, flat_shelf_width, c1, slope_pos
161 real :: edge_pos, shelf_slope_scale, Rho_ocean
162 integer :: i, j, jsc, jec, jsd, jed, jedg, nyh, isc, iec, isd, ied
165 jsc = g%jsc ; jec = g%jec ; isc = g%isc ; iec = g%iec
166 jsd = g%jsd ; jed = g%jed ; isd = g%isd ; ied = g%ied
167 nyh = g%domain%njhalo ; jedg = g%domain%njglobal+nyh
170 call mom_mesg(mdl//
": setting thickness")
172 call get_param(pf, mdl,
"SHELF_MAX_DRAFT", max_draft, &
173 units=
"m", default=1.0, scale=us%m_to_Z)
174 call get_param(pf, mdl,
"SHELF_MIN_DRAFT", min_draft, &
175 units=
"m", default=1.0, scale=us%m_to_Z)
176 call get_param(pf, mdl,
"FLAT_SHELF_WIDTH", flat_shelf_width, &
177 units=
"axis_units", default=0.0)
178 call get_param(pf, mdl,
"SHELF_SLOPE_SCALE", shelf_slope_scale, &
179 units=
"axis_units", default=0.0)
180 call get_param(pf, mdl,
"SHELF_EDGE_POS_0", edge_pos, &
181 units=
"axis_units", default=0.0)
189 slope_pos = edge_pos - flat_shelf_width
190 c1 = 0.0 ;
if (shelf_slope_scale > 0.0) c1 = 1.0 / shelf_slope_scale
195 if (((j+j_off) <= jedg) .AND. ((j+j_off) >= nyh+1))
then
199 if ((j >= jsc) .and. (j <= jec))
then
201 if (g%geoLonCu(i-1,j) >= edge_pos)
then
204 area_shelf_h(i,j) = 0.0
208 if (g%geoLonCu(i,j) > edge_pos)
then
209 area_shelf_h(i,j) = us%L_to_m**2*g%areaT(i,j) * (edge_pos - g%geoLonCu(i-1,j)) / &
210 (g%geoLonCu(i,j) - g%geoLonCu(i-1,j))
213 area_shelf_h(i,j) = us%L_to_m**2*g%areaT(i,j)
217 if (g%geoLonT(i,j) > slope_pos)
then
218 h_shelf(i,j) = min_draft
224 h_shelf(i,j) = (min_draft + &
225 (max_draft - min_draft) * &
226 min(1.0, (c1*(slope_pos - g%geoLonT(i,j)))**2) )
232 if ((i+g%idg_offset) == g%domain%nihalo+1)
then