Initialize ice shelf thickness for a channel configuration.
148 type(ocean_grid_type),
intent(in) :: G
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
156 type(unit_scale_type),
intent(in) :: US
157 type(param_file_type),
intent(in) :: PF
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) = 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) = 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