ISOMIP_initialization.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!> Configures the ISOMIP test case.
7
8use mom_ale_sponge, only : ale_sponge_cs, set_up_ale_sponge_field, initialize_ale_sponge
11use mom_error_handler, only : mom_mesg, mom_error, fatal, is_root_pe, warning
12use mom_file_parser, only : get_param, log_version, param_file_type
13use mom_get_input, only : directories
14use mom_grid, only : ocean_grid_type
15use mom_io, only : file_exists, mom_read_data, slasher
19use mom_eos, only : calculate_density, calculate_density_derivs, eos_type
20use regrid_consts, only : coordinatemode, default_coordinate_mode
21use regrid_consts, only : regridding_layer, regridding_zstar
22use regrid_consts, only : regridding_rho, regridding_sigma
23use regrid_consts, only : regridding_sigma_shelf_zstar
24implicit none ; private
25
26#include <MOM_memory.h>
27
28character(len=40) :: mdl = "ISOMIP_initialization" !< This module's name.
29
30! The following routines are visible to the outside world
35
36! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
37! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
38! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
39! vary with the Boussinesq approximation, the Boussinesq variant is given first.
40
41contains
42
43!> Initialization of topography for the ISOMIP configuration
44subroutine isomip_initialize_topography(D, G, param_file, max_depth, US)
45 type(dyn_horgrid_type), intent(in) :: g !< The dynamic horizontal grid type
46 real, dimension(G%isd:G%ied,G%jsd:G%jed), &
47 intent(out) :: d !< Ocean bottom depth [Z ~> m]
48 type(param_file_type), intent(in) :: param_file !< Parameter file structure
49 real, intent(in) :: max_depth !< Maximum model depth [Z ~> m]
50 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
51
52 ! Local variables
53 real :: min_depth ! The minimum depth of the ocean [Z ~> m].
54 ! The following variables are used to set up the bathymetry in the ISOMIP example.
55 real :: bmax ! maximum depth of bedrock topography [Z ~> m]
56 real :: b0, b2, b4, b6 ! first, second, third and fourth bedrock topography coeffs [Z ~> m]
57 real :: xbar ! characteristic along-flow length scale of the bedrock [L ~> m]
58 real :: dc ! depth of the trough compared with side walls [Z ~> m].
59 real :: fc ! characteristic width of the side walls of the channel [L ~> m]
60 real :: wc ! half-width of the trough [L ~> m]
61 real :: ly ! domain width (across ice flow) [L ~> m]
62 real :: bx, by ! The x- and y- contributions to the bathymetric profiles at a point [Z ~> m]
63 real :: xtil ! x-positon normalized by the characteristic along-flow length scale [nondim]
64 logical :: is_2d ! If true, use a 2D setup
65 ! This include declares and sets the variable "version".
66# include "version_variable.h"
67 character(len=40) :: mdl = "ISOMIP_initialize_topography" ! This subroutine's name.
68 integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
69 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
70 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
71
72 call mom_mesg(" ISOMIP_initialization.F90, ISOMIP_initialize_topography: setting topography", 5)
73
74 call log_version(param_file, mdl, version, "")
75 call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, &
76 "The minimum depth of the ocean.", units="m", default=0.0, scale=us%m_to_Z)
77 call get_param(param_file, mdl, "ISOMIP_2D", is_2d, 'If true, use a 2D setup.', default=.false.)
78 call get_param(param_file, mdl, "ISOMIP_MAX_BEDROCK", bmax, &
79 "Maximum depth of bedrock topography in the ISOMIP configuration.", &
80 units="m", default=720.0, scale=us%m_to_Z)
81 call get_param(param_file, mdl, "ISOMIP_TROUGH_DEPTH", dc, &
82 "Depth of the trough compared with side walls in the ISOMIP configuration.", &
83 units="m", default=500.0, scale=us%m_to_Z)
84 call get_param(param_file, mdl, "ISOMIP_BEDROCK_LENGTH", xbar, &
85 "Characteristic along-flow length scale of the bedrock in the ISOMIP configuration.", &
86 units="m", default=300.0e3, scale=us%m_to_L)
87 call get_param(param_file, mdl, "ISOMIP_TROUGH_WIDTH", wc, &
88 "Half-width of the trough in the ISOMIP configuration.", &
89 units="m", default=24.0e3, scale=us%m_to_L)
90 call get_param(param_file, mdl, "ISOMIP_DOMAIN_WIDTH", ly, &
91 "Domain width (across ice flow) in the ISOMIP configuration.", &
92 units="m", default=80.0e3, scale=us%m_to_L)
93 call get_param(param_file, mdl, "ISOMIP_SIDE_WIDTH", fc, &
94 "Characteristic width of the side walls of the channel in the ISOMIP configuration.", &
95 units="m", default=4.0e3, scale=us%m_to_L)
96
97 if (g%grid_unit_to_L <= 0.) call mom_error(fatal, "ISOMIP_initialization.F90: " //&
98 "ISOMIP_initialize_topography is only set to work with Cartesian axis units.")
99
100 ! The following variables should be transformed into runtime parameters.
101 b0 = -150.0*us%m_to_Z ; b2 = -728.8*us%m_to_Z ; b4 = 343.91*us%m_to_Z ; b6 = -50.57*us%m_to_Z
102
103 if (is_2d) then
104 do j=js,je ; do i=is,ie
105 ! For the 2D setup take a slice through the middle of the domain
106 xtil = g%geoLonT(i,j)*g%grid_unit_to_L / xbar
107 ! xtil = 450.0e3*US%m_to_L / xbar
108 bx = b0 + b2*xtil**2 + b4*xtil**4 + b6*xtil**6
109
110 by = 2.0 * dc / (1.0 + exp(2.0*wc / fc))
111
112 d(i,j) = -max(bx+by, -bmax)
113 if (d(i,j) > max_depth) d(i,j) = max_depth
114 if (d(i,j) < min_depth) d(i,j) = 0.5*min_depth
115 enddo ; enddo
116
117 else
118 do j=js,je ; do i=is,ie
119 ! 3D setup
120 ! ===== TEST =====
121 !if (G%geoLonT(i,j)<500.) then
122 ! xtil = 500.0e3*US%m_to_L / xbar
123 !else
124 ! xtil = G%geoLonT(i,j)*G%grid_unit_to_L / xbar
125 !endif
126 ! ===== TEST =====
127
128 xtil = g%geoLonT(i,j)*g%grid_unit_to_L / xbar
129
130 bx = b0 + b2*xtil**2 + b4*xtil**4 + b6*xtil**6
131 by = (dc / (1.0 + exp(-2.*(g%geoLatT(i,j)*g%grid_unit_to_L - 0.5*ly - wc) / fc))) + &
132 (dc / (1.0 + exp(2.*(g%geoLatT(i,j)*g%grid_unit_to_L - 0.5*ly + wc) / fc)))
133
134 d(i,j) = -max(bx+by, -bmax)
135 if (d(i,j) > max_depth) d(i,j) = max_depth
136 if (d(i,j) < min_depth) d(i,j) = 0.5*min_depth
137 enddo ; enddo
138 endif
139
140end subroutine isomip_initialize_topography
141
142!> Initialization of thicknesses
143subroutine isomip_initialize_thickness ( h, depth_tot, G, GV, US, param_file, tv, just_read)
144 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
145 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
146 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
147 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
148 intent(out) :: h !< The thickness that is being initialized [Z ~> m]
149 real, dimension(SZI_(G),SZJ_(G)), &
150 intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
151 type(param_file_type), intent(in) :: param_file !< A structure to parse for model parameter values
152 type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any
153 !! available thermodynamic fields, including
154 !! the eqn. of state.
155 logical, intent(in) :: just_read !< If true, this call will only read
156 !! parameters without changing h.
157 ! Local variables
158 real :: e0(szk_(gv)+1) ! The resting interface heights, in depth units [Z ~> m],
159 ! usually negative because it is positive upward.
160 real :: eta1d(szk_(gv)+1) ! Interface height relative to the sea surface
161 ! positive upward, in depth units [Z ~> m].
162 integer :: i, j, k, is, ie, js, je, nz
163 real :: min_thickness ! Minimum layer thicknesses [Z ~> m]
164 real :: s_sur, s_bot ! Surface and bottom salinities [S ~> ppt]
165 real :: t_sur, t_bot ! Surface and bottom temperatures [C ~> degC]
166 real :: rho_sur, rho_bot ! Surface and bottom densities [R ~> kg m-3]
167 real :: rho_range ! The range of densities [R ~> kg m-3]
168 !character(len=256) :: mesg ! The text of an error message
169 character(len=40) :: verticalcoordinate
170
171 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
172
173 if (.not.just_read) &
174 call mom_mesg("ISOMIP_initialization.F90, ISOMIP_initialize_thickness: setting thickness")
175
176 call get_param(param_file, mdl,"MIN_THICKNESS", min_thickness, &
177 'Minimum layer thickness', units='m', default=1.e-3, do_not_log=just_read, scale=us%m_to_Z)
178 call get_param(param_file, mdl,"REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
179 default=default_coordinate_mode, do_not_log=just_read)
180
181 select case ( coordinatemode(verticalcoordinate) )
182
183 case ( regridding_layer, regridding_rho ) ! Initial thicknesses for isopycnal coordinates
184 call get_param(param_file, mdl, "ISOMIP_T_SUR", t_sur, &
185 "Temperature at the surface (interface)", &
186 units="degC", default=-1.9, scale=us%degC_to_C, do_not_log=just_read)
187 call get_param(param_file, mdl, "ISOMIP_S_SUR", s_sur, &
188 "Salinity at the surface (interface)", &
189 units="ppt", default=33.8, scale=us%ppt_to_S, do_not_log=just_read)
190 call get_param(param_file, mdl, "ISOMIP_T_BOT", t_bot, &
191 "Temperature at the bottom (interface)", &
192 units="degC", default=-1.9, scale=us%degC_to_C, do_not_log=just_read)
193 call get_param(param_file, mdl, "ISOMIP_S_BOT", s_bot, &
194 "Salinity at the bottom (interface)", &
195 units="ppt", default=34.55, scale=us%ppt_to_S, do_not_log=just_read)
196
197 if (just_read) return ! All run-time parameters have been read, so return.
198
199 ! Compute min/max density using T_SUR/S_SUR and T_BOT/S_BOT
200 call calculate_density(t_sur, s_sur, 0.0, rho_sur, tv%eqn_of_state)
201 ! write(mesg,*) 'Surface density is:', rho_sur
202 ! call MOM_mesg(mesg,5)
203 call calculate_density(t_bot, s_bot, 0.0, rho_bot, tv%eqn_of_state)
204 ! write(mesg,*) 'Bottom density is:', rho_bot
205 ! call MOM_mesg(mesg,5)
206 rho_range = rho_bot - rho_sur
207 ! write(mesg,*) 'Density range is:', rho_range
208 ! call MOM_mesg(mesg,5)
209
210 ! Construct notional interface positions
211 e0(1) = 0.
212 do k=2,nz
213 e0(k) = -g%max_depth * ( 0.5 * ( gv%Rlay(k-1) + gv%Rlay(k) ) - rho_sur ) / rho_range
214 e0(k) = min( 0., e0(k) ) ! Bound by surface
215 e0(k) = max( -g%max_depth, e0(k) ) ! Bound by possible deepest point in model
216 ! write(mesg,*) 'G%max_depth,GV%Rlay(k-1),GV%Rlay(k),e0(k)', &
217 ! G%max_depth,GV%Rlay(k-1),GV%Rlay(k),e0(k)
218 ! call MOM_mesg(mesg,5)
219 enddo
220 e0(nz+1) = -g%max_depth
221
222 ! Calculate thicknesses
223 do j=js,je ; do i=is,ie
224 eta1d(nz+1) = -depth_tot(i,j)
225 do k=nz,1,-1
226 eta1d(k) = e0(k)
227 if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z)) then
228 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
229 h(i,j,k) = gv%Angstrom_Z
230 else
231 h(i,j,k) = eta1d(k) - eta1d(k+1)
232 endif
233 enddo
234 enddo ; enddo
235
236 case ( regridding_zstar, regridding_sigma_shelf_zstar ) ! Initial thicknesses for z coordinates
237 if (just_read) return ! All run-time parameters have been read, so return.
238 do j=js,je ; do i=is,ie
239 eta1d(nz+1) = -depth_tot(i,j)
240 do k=nz,1,-1
241 eta1d(k) = -g%max_depth * real(k-1) / real(nz)
242 if (eta1d(k) < (eta1d(k+1) + min_thickness)) then
243 eta1d(k) = eta1d(k+1) + min_thickness
244 h(i,j,k) = min_thickness
245 else
246 h(i,j,k) = eta1d(k) - eta1d(k+1)
247 endif
248 enddo
249 enddo ; enddo
250
251 case ( regridding_sigma ) ! Initial thicknesses for sigma coordinates
252 if (just_read) return ! All run-time parameters have been read, so return.
253 do j=js,je ; do i=is,ie
254 h(i,j,:) = depth_tot(i,j) / real(nz)
255 enddo ; enddo
256
257 case default
258 call mom_error(fatal,"isomip_initialize: "// &
259 "Unrecognized i.c. setup - set REGRIDDING_COORDINATE_MODE")
260
261 end select
262
263end subroutine isomip_initialize_thickness
264
265!> Initial values for temperature and salinity
266subroutine isomip_initialize_temperature_salinity ( T, S, h, depth_tot, G, GV, US, param_file, &
267 eqn_of_state, just_read)
268 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
269 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
270 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
271 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: t !< Potential temperature [C ~> degC]
272 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: s !< Salinity [S ~> ppt]
273 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [Z ~> m]
274 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: depth_tot !< The nominal total bottom-to-top
275 !! depth of the ocean [Z ~> m]
276 type(param_file_type), intent(in) :: param_file !< Parameter file structure
277 type(eos_type), intent(in) :: eqn_of_state !< Equation of state structure
278 logical, intent(in) :: just_read !< If true, this call will
279 !! only read parameters without changing T & S.
280 ! Local variables
281 real :: rho_sur, rho_bot ! Surface and bottom densities [R ~> kg m-3]
282 real :: xi0, xi1 ! Heights in depth units [Z ~> m].
283 real :: s_sur, s_bot ! Salinity at the surface and bottom [S ~> ppt]
284 real :: t_sur, t_bot ! Temperature at the surface and bottom [C ~> degC]
285 real :: dt_dz ! Vertical gradient of temperature [C Z-1 ~> degC m-1].
286 real :: ds_dz ! Vertical gradient of salinity [S Z-1 ~> ppt m-1].
287 real :: t0(szk_(gv)) ! A profile of temperatures [C ~> degC]
288 real :: s0(szk_(gv)) ! A profile of salinities [S ~> ppt]
289 real :: drho_dt(szk_(gv)) ! Derivative of density with temperature [R C-1 ~> kg m-3 degC-1].
290 real :: drho_ds(szk_(gv)) ! Derivative of density with salinity [R S-1 ~> kg m-3 ppt-1].
291 real :: rho_guess(szk_(gv)) ! Potential density at T0 & S0 [R ~> kg m-3].
292 real :: pres(szk_(gv)) ! An array of the reference pressure [R L2 T-2 ~> Pa]. (zero here)
293 real :: drho_dt1 ! A prescribed derivative of density with temperature [R C-1 ~> kg m-3 degC-1]
294 real :: drho_ds1 ! A prescribed derivative of density with salinity [R S-1 ~> kg m-3 ppt-1].
295 real :: t_ref ! Default value for other temperatures [C ~> degC]
296 real :: s_ref ! Default value for other salinities [S ~> ppt]
297 logical :: fit_salin ! If true, accept the prescribed temperature and fit the salinity.
298 !real :: rho_tmp ! A temporary density used for debugging [R ~> kg m-3]
299 !character(len=256) :: mesg ! The text of an error message
300 character(len=40) :: verticalcoordinate
301 integer :: i, j, k, is, ie, js, je, nz, itt
302
303 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
304 pres(:) = 0.0
305
306 call get_param(param_file, mdl, "REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
307 default=default_coordinate_mode, do_not_log=just_read)
308 call get_param(param_file, mdl, "ISOMIP_T_SUR",t_sur, &
309 "Temperature at the surface (interface)", &
310 units="degC", default=-1.9, scale=us%degC_to_C, do_not_log=just_read)
311 call get_param(param_file, mdl, "ISOMIP_S_SUR", s_sur, &
312 "Salinity at the surface (interface)", &
313 units="ppt", default=33.8, scale=us%ppt_to_S, do_not_log=just_read)
314 call get_param(param_file, mdl, "ISOMIP_T_BOT", t_bot, &
315 "Temperature at the bottom (interface)", &
316 units="degC", default=-1.9, scale=us%degC_to_C, do_not_log=just_read)
317 call get_param(param_file, mdl, "ISOMIP_S_BOT", s_bot, &
318 "Salinity at the bottom (interface)", &
319 units="ppt", default=34.55, scale=us%ppt_to_S, do_not_log=just_read)
320
321 call calculate_density(t_sur, s_sur, 0.0, rho_sur, eqn_of_state)
322 ! write(mesg,*) 'Density in the surface layer:', rho_sur
323 ! call MOM_mesg(mesg,5)
324 call calculate_density(t_bot, s_bot, 0.0, rho_bot, eqn_of_state)
325 ! write(mesg,*) 'Density in the bottom layer::', rho_bot
326 ! call MOM_mesg(mesg,5)
327
328 select case ( coordinatemode(verticalcoordinate) )
329
330 case ( regridding_rho, regridding_zstar, regridding_sigma_shelf_zstar, regridding_sigma )
331 if (just_read) return ! All run-time parameters have been read, so return.
332
333 ds_dz = (s_sur - s_bot) / g%max_depth
334 dt_dz = (t_sur - t_bot) / g%max_depth
335 do j=js,je ; do i=is,ie
336 xi0 = -depth_tot(i,j)
337 do k = nz,1,-1
338 xi0 = xi0 + 0.5 * h(i,j,k) ! Depth in middle of layer
339 s(i,j,k) = s_sur + ds_dz * xi0
340 t(i,j,k) = t_sur + dt_dz * xi0
341 xi0 = xi0 + 0.5 * h(i,j,k) ! Depth at top of layer
342 enddo
343 enddo ; enddo
344
345 case ( regridding_layer )
346 call get_param(param_file, mdl, "FIT_SALINITY", fit_salin, &
347 "If true, accept the prescribed temperature and fit the "//&
348 "salinity; otherwise take salinity and fit temperature.", &
349 default=.false., do_not_log=just_read)
350 call get_param(param_file, mdl, "DRHO_DS", drho_ds1, &
351 "Partial derivative of density with salinity.", &
352 units="kg m-3 ppt-1", scale=us%kg_m3_to_R*us%S_to_ppt, &
353 fail_if_missing=.not.just_read, do_not_log=just_read)
354 call get_param(param_file, mdl, "DRHO_DT", drho_dt1, &
355 "Partial derivative of density with temperature.", &
356 units="kg m-3 K-1", scale=us%kg_m3_to_R*us%C_to_degC, &
357 fail_if_missing=.not.just_read, do_not_log=just_read)
358 call get_param(param_file, mdl, "T_REF", t_ref, &
359 "A reference temperature used in initialization.", &
360 units="degC", scale=us%degC_to_C, fail_if_missing=.not.just_read, do_not_log=just_read)
361 call get_param(param_file, mdl, "S_REF", s_ref, &
362 "A reference salinity used in initialization.", &
363 units="ppt", default=35.0, scale=us%ppt_to_S, do_not_log=just_read)
364 if (just_read) return ! All run-time parameters have been read, so return.
365
366 ! write(mesg,*) 'read drho_dS, drho_dT', drho_dS1, drho_dT1
367 ! call MOM_mesg(mesg,5)
368
369 ds_dz = (s_sur - s_bot) / g%max_depth
370 dt_dz = (t_sur - t_bot) / g%max_depth
371
372 do j=js,je ; do i=is,ie
373 xi0 = 0.0
374 do k = 1,nz
375 !T0(k) = T_Ref ; S0(k) = S_Ref
376 xi1 = xi0 + 0.5 * h(i,j,k)
377 s0(k) = s_sur - ds_dz * xi1
378 t0(k) = t_sur - dt_dz * xi1
379 xi0 = xi0 + h(i,j,k)
380 ! write(mesg,*) 'S,T,xi0,xi1,k',S0(k),T0(k),xi0,xi1,k
381 ! call MOM_mesg(mesg,5)
382 enddo
383
384 call calculate_density_derivs(t0, s0, pres, drho_dt, drho_ds, eqn_of_state, (/1,1/) )
385 ! write(mesg,*) 'computed drho_dS, drho_dT', drho_dS(1), drho_dT(1)
386 ! call MOM_mesg(mesg,5)
387 call calculate_density(t0(1), s0(1), pres(1), rho_guess(1), eqn_of_state)
388
389 if (fit_salin) then
390 ! A first guess of the layers' salinity.
391 do k=nz,1,-1
392 s0(k) = max(0.0, s0(1) + (gv%Rlay(k) - rho_guess(1)) / drho_ds1)
393 enddo
394 ! Refine the guesses for each layer.
395 do itt=1,6
396 call calculate_density(t0, s0, pres, rho_guess, eqn_of_state)
397 call calculate_density_derivs(t0, s0, pres, drho_dt, drho_ds, eqn_of_state)
398 do k=1,nz
399 s0(k) = max(0.0, s0(k) + (gv%Rlay(k) - rho_guess(k)) / drho_ds1)
400 enddo
401 enddo
402
403 else
404 ! A first guess of the layers' temperatures.
405 do k=nz,1,-1
406 t0(k) = t0(1) + (gv%Rlay(k) - rho_guess(1)) / drho_dt1
407 enddo
408
409 do itt=1,6
410 call calculate_density(t0, s0, pres, rho_guess, eqn_of_state)
411 call calculate_density_derivs(t0, s0, pres, drho_dt, drho_ds, eqn_of_state)
412 do k=1,nz
413 t0(k) = t0(k) + (gv%Rlay(k) - rho_guess(k)) / drho_dt(k)
414 enddo
415 enddo
416 endif
417
418 do k=1,nz
419 t(i,j,k) = t0(k) ; s(i,j,k) = s0(k)
420 enddo
421
422 enddo ; enddo
423
424 case default
425 call mom_error(fatal,"isomip_initialize: "// &
426 "Unrecognized i.c. setup - set REGRIDDING_COORDINATE_MODE")
427
428 end select
429
430 ! for debugging
431 !i = G%iec ; j = G%jec
432 !do k = 1,nz
433 ! call calculate_density(T(i,j,k), S(i,j,k),0.0,rho_tmp,eqn_of_state, scale=US%kg_m3_to_R)
434 ! write(mesg,*) 'k,h,T,S,rho,Rlay',k,US%Z_to_m*h(i,j,k),US%C_to_degC*T(i,j,k),US%S_to_ppt*S(i,j,k),rho_tmp,GV%Rlay(k)
435 ! call MOM_mesg(mesg,5)
436 !enddo
437
439
440!> Sets up the inverse restoration time (Idamp), and
441!! the values towards which the interface heights and an arbitrary
442!! number of tracers should be restored within each sponge.
443subroutine isomip_initialize_sponges(G, GV, US, tv, depth_tot, PF, use_ALE, CSp, ACSp)
444 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
445 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
446 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
447 type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any available
448 !! thermodynamic fields, potential temperature and
449 !! salinity or mixed layer density.
450 !! Absent fields have NULL ptrs.
451 real, dimension(SZI_(G),SZJ_(G)), &
452 intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
453 type(param_file_type), intent(in) :: pf !< A structure to parse for model parameter values
454 logical, intent(in) :: use_ale !< If true, indicates model is in ALE mode
455 type(sponge_cs), pointer :: csp !< Layer-mode sponge structure
456 type(ale_sponge_cs), pointer :: acsp !< ALE-mode sponge structure
457 ! Local variables
458 real :: t(szi_(g),szj_(g),szk_(gv)) ! A temporary array for temp [C ~> degC]
459 real :: s(szi_(g),szj_(g),szk_(gv)) ! A temporary array for salt [S ~> ppt]
460 ! real :: RHO(SZI_(G),SZJ_(G),SZK_(GV)) ! A temporary array for RHO [R ~> kg m-3]
461 real :: dz(szi_(g),szj_(g),szk_(gv)) ! Sponge layer thicknesses in height units [Z ~> m]
462 real :: idamp(szi_(g),szj_(g)) ! The sponge damping rate [T-1 ~> s-1]
463 real :: tnudg ! Nudging time scale [T ~> s]
464 real :: s_sur, s_bot ! Surface and bottom salinities in the sponge region [S ~> ppt]
465 real :: t_sur, t_bot ! Surface and bottom temperatures in the sponge region [C ~> degC]
466 real :: t_ref ! Default value for other temperatures [C ~> degC]
467 real :: s_ref ! Default value for other salinities [S ~> ppt]
468 real :: rho_sur, rho_bot ! Surface and bottom densities [R ~> kg m-3]
469 real :: rho_range ! The range of densities [R ~> kg m-3]
470 real :: dt_dz ! Vertical gradient of temperature [C Z-1 ~> degC m-1]
471 real :: ds_dz ! Vertical gradient of salinity [S Z-1 ~> ppt m-1]
472
473 real :: e0(szk_(gv)+1) ! The resting interface heights [Z ~> m], usually
474 ! negative because it is positive upward.
475 real :: eta1d(szk_(gv)+1) ! Interface height relative to the sea surface, positive upward [Z ~> m].
476 real :: eta(szi_(g),szj_(g),szk_(gv)+1) ! A temporary array for interface heights [Z ~> m].
477 real :: min_depth ! The minimum depth of the ocean [Z ~> m]
478 real :: min_thickness ! The minimum layer thickness [Z ~> m]
479 real :: xi0 ! Interface heights in depth units [Z ~> m], usually negative.
480 !real :: rho_tmp ! A temporary density used for debugging [R ~> kg m-3]
481 character(len=40) :: verticalcoordinate, filename, state_file
482 character(len=40) :: temp_var, salt_var, eta_var, inputdir
483
484 character(len=40) :: mdl = "ISOMIP_initialize_sponges" ! This subroutine's name.
485 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
486
487 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
488 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
489
490 call get_param(pf, mdl, "MIN_THICKNESS", min_thickness, "Minimum layer thickness", &
491 units="m", default=1.e-3, scale=us%m_to_Z)
492
493 call get_param(pf, mdl, "REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
494 default=default_coordinate_mode)
495
496 call get_param(pf, mdl, "ISOMIP_TNUDG", tnudg, "Nudging time scale for sponge layers", &
497 units="days", default=0.0, scale=86400.0*us%s_to_T)
498
499 call get_param(pf, mdl, "T_REF", t_ref, "Reference temperature", &
500 units="degC", default=10.0, scale=us%degC_to_C, do_not_log=.true.)
501
502 call get_param(pf, mdl, "S_REF", s_ref, "Reference salinity", &
503 units="ppt", default=35.0, scale=us%ppt_to_S, do_not_log=.true.)
504
505 call get_param(pf, mdl, "ISOMIP_S_SUR_SPONGE", s_sur, &
506 "Surface salinity in sponge layer.", &
507 units="ppt", default=us%S_to_ppt*s_ref, scale=us%ppt_to_S)
508
509 call get_param(pf, mdl, "ISOMIP_S_BOT_SPONGE", s_bot, &
510 "Bottom salinity in sponge layer.", &
511 units="ppt", default=us%S_to_ppt*s_ref, scale=us%ppt_to_S)
512
513 call get_param(pf, mdl, "ISOMIP_T_SUR_SPONGE", t_sur, &
514 "Surface temperature in sponge layer.", &
515 units="degC", default=us%C_to_degC*t_ref, scale=us%degC_to_C)
516
517 call get_param(pf, mdl, "ISOMIP_T_BOT_SPONGE", t_bot, &
518 "Bottom temperature in sponge layer.", &
519 units="degC", default=us%C_to_degC*t_ref, scale=us%degC_to_C)
520
521 t(:,:,:) = 0.0 ; s(:,:,:) = 0.0 ; idamp(:,:) = 0.0 !; RHO(:,:,:) = 0.0
522
523! Set up sponges for ISOMIP configuration
524 call get_param(pf, mdl, "MINIMUM_DEPTH", min_depth, &
525 "The minimum depth of the ocean.", units="m", default=0.0, scale=us%m_to_Z)
526
527 if (associated(csp)) call mom_error(fatal, &
528 "ISOMIP_initialize_sponges called with an associated control structure.")
529 if (associated(acsp)) call mom_error(fatal, &
530 "ISOMIP_initialize_sponges called with an associated ALE-sponge control structure.")
531
532 ! Here the inverse damping time [T-1 ~> s-1], is set. Set Idamp to 0
533 ! wherever there is no sponge, and the subroutines that are called
534 ! will automatically set up the sponges only where Idamp is positive
535 ! and mask2dT is 1.
536
537 do j=js,je ; do i=is,ie
538 if (depth_tot(i,j) <= min_depth) then
539 idamp(i,j) = 0.0
540 elseif (g%geoLonT(i,j) >= 790.0 .AND. g%geoLonT(i,j) <= 800.0) then
541 idamp(i,j) = (1.0/tnudg) * max(0.0, (g%geoLonT(i,j)-790.0) / (800.0-790.0))
542 else
543 idamp(i,j) = 0.0
544 endif
545
546 enddo ; enddo
547
548 ! Compute min/max density using T_SUR/S_SUR and T_BOT/S_BOT
549 call calculate_density(t_sur, s_sur, 0.0, rho_sur, tv%eqn_of_state)
550 !write (mesg,*) 'Surface density in sponge:', rho_sur
551 ! call MOM_mesg(mesg,5)
552 call calculate_density(t_bot, s_bot, 0.0, rho_bot, tv%eqn_of_state)
553 !write (mesg,*) 'Bottom density in sponge:', rho_bot
554 ! call MOM_mesg(mesg,5)
555 rho_range = rho_bot - rho_sur
556 !write (mesg,*) 'Density range in sponge:', rho_range
557 ! call MOM_mesg(mesg,5)
558
559 if (use_ale) then
560
561 select case ( coordinatemode(verticalcoordinate) )
562
563 case ( regridding_rho )
564 ! Construct notional interface positions
565 e0(1) = 0.
566 do k=2,nz
567 e0(k) = -g%max_depth * ( 0.5 * ( gv%Rlay(k-1) + gv%Rlay(k) ) - rho_sur ) / rho_range
568 e0(k) = min( 0., e0(k) ) ! Bound by surface
569 e0(k) = max( -g%max_depth, e0(k) ) ! Bound by possible deepest point in model
570 ! write(mesg,*) 'G%max_depth,GV%Rlay(k-1),GV%Rlay(k),e0(k)',&
571 ! G%max_depth,GV%Rlay(k-1),GV%Rlay(k),e0(k)
572 ! call MOM_mesg(mesg,5)
573 enddo
574 e0(nz+1) = -g%max_depth
575
576 ! Calculate thicknesses
577 do j=js,je ; do i=is,ie
578 eta1d(nz+1) = -depth_tot(i,j)
579 do k=nz,1,-1
580 eta1d(k) = e0(k)
581 if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z)) then
582 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
583 dz(i,j,k) = gv%Angstrom_Z
584 else
585 dz(i,j,k) = eta1d(k) - eta1d(k+1)
586 endif
587 enddo
588 enddo ; enddo
589
590 case ( regridding_zstar, regridding_sigma_shelf_zstar ) ! Initial thicknesses for z coordinates
591 do j=js,je ; do i=is,ie
592 eta1d(nz+1) = -depth_tot(i,j)
593 do k=nz,1,-1
594 eta1d(k) = -g%max_depth * real(k-1) / real(nz)
595 if (eta1d(k) < (eta1d(k+1) + min_thickness)) then
596 eta1d(k) = eta1d(k+1) + min_thickness
597 dz(i,j,k) = min_thickness
598 else
599 dz(i,j,k) = eta1d(k) - eta1d(k+1)
600 endif
601 enddo
602 enddo ; enddo
603
604 case ( regridding_sigma ) ! Initial thicknesses for sigma coordinates
605 do j=js,je ; do i=is,ie
606 dz(i,j,:) = depth_tot(i,j) / real(nz)
607 enddo ; enddo
608
609 case default
610 call mom_error(fatal,"ISOMIP_initialize_sponges: "// &
611 "Unrecognized i.c. setup - set REGRIDDING_COORDINATE_MODE")
612
613 end select
614
615 ds_dz = (s_sur - s_bot) / g%max_depth
616 dt_dz = (t_sur - t_bot) / g%max_depth
617 do j=js,je ; do i=is,ie
618 xi0 = -depth_tot(i,j)
619 do k = nz,1,-1
620 xi0 = xi0 + 0.5 * dz(i,j,k) ! Depth in middle of layer
621 s(i,j,k) = s_sur + ds_dz * xi0
622 t(i,j,k) = t_sur + dt_dz * xi0
623 xi0 = xi0 + 0.5 * dz(i,j,k) ! Depth at top of layer
624 enddo
625 enddo ; enddo
626
627 ! for debugging
628 !i = G%iec ; j = G%jec
629 !do k = 1,nz
630 ! call calculate_density(T(i,j,k), S(i,j,k), 0.0, rho_tmp, tv%eqn_of_state, scale=US%kg_m3_to_R)
631 ! write(mesg,*) 'Sponge - k,h,T,S,rho,Rlay',k,h(i,j,k),T(i,j,k),S(i,j,k),rho_tmp,GV%Rlay(k)
632 ! call MOM_mesg(mesg,5)
633 !enddo
634
635 ! This call sets up the damping rates and interface heights in the sponges.
636 call initialize_ale_sponge(idamp, g, gv, pf, acsp, dz, nz, data_h_is_z=.true.)
637
638 ! Now register all of the fields which are damped in the sponge. !
639 ! By default, momentum is advected vertically within the sponge, but !
640 ! momentum is typically not damped within the sponge. !
641
642 ! The remaining calls to set_up_sponge_field can be in any order. !
643 if ( associated(tv%T) ) call set_up_ale_sponge_field(t, g, gv, tv%T, acsp, 'temp', &
644 sp_long_name='temperature', sp_unit='degC s-1')
645 if ( associated(tv%S) ) call set_up_ale_sponge_field(s, g, gv, tv%S, acsp, 'salt', &
646 sp_long_name='salinity', sp_unit='g kg-1 s-1')
647
648
649 else ! layer mode
650 ! 1) Read eta, salt and temp from IC file
651 call get_param(pf, mdl, "INPUTDIR", inputdir, default=".")
652 inputdir = slasher(inputdir)
653 ! GM: get two different files, one with temp and one with salt values
654 ! this is work around to avoid having wrong values near the surface
655 ! because of the FIT_SALINITY option. To get salt values right in the
656 ! sponge, FIT_SALINITY=False. The oposite is true for temp. One can
657 ! combined the *correct* temp and salt values in one file instead.
658 call get_param(pf, mdl, "ISOMIP_SPONGE_FILE", state_file, &
659 "The name of the file with temps., salts. and interfaces to "//&
660 "damp toward.", fail_if_missing=.true.)
661 call get_param(pf, mdl, "SPONGE_PTEMP_VAR", temp_var, &
662 "The name of the potential temperature variable in "//&
663 "SPONGE_STATE_FILE.", default="Temp")
664 call get_param(pf, mdl, "SPONGE_SALT_VAR", salt_var, &
665 "The name of the salinity variable in "//&
666 "SPONGE_STATE_FILE.", default="Salt")
667 call get_param(pf, mdl, "SPONGE_ETA_VAR", eta_var, &
668 "The name of the interface height variable in "//&
669 "SPONGE_STATE_FILE.", default="eta")
670
671 !read temp and eta
672 filename = trim(inputdir)//trim(state_file)
673 if (.not.file_exists(filename, g%Domain)) call mom_error(fatal, &
674 "ISOMIP_initialize_sponges: Unable to open "//trim(filename))
675 call mom_read_data(filename, eta_var, eta(:,:,:), g%Domain, scale=us%m_to_Z)
676 call mom_read_data(filename, temp_var, t(:,:,:), g%Domain, scale=us%degC_to_C)
677 call mom_read_data(filename, salt_var, s(:,:,:), g%Domain, scale=us%ppt_to_S)
678
679 ! for debugging
680 !i = G%iec ; j = G%jec
681 !do k = 1,nz
682 ! call calculate_density(T(i,j,k), S(i,j,k), 0.0, rho_tmp, tv%eqn_of_state, scale=US%kg_m3_to_R)
683 ! write(mesg,*) 'Sponge - k,eta,T,S,rho,Rlay',k,eta(i,j,k),T(i,j,k),&
684 ! S(i,j,k),rho_tmp,GV%Rlay(k)
685 ! call MOM_mesg(mesg,5)
686 !enddo
687
688 ! Set the sponge damping rates so that the model will know where to
689 ! apply the sponges, along with the interface heights.
690 call initialize_sponge(idamp, eta, g, pf, csp, gv)
691 ! Apply sponge in tracer fields
692 call set_up_sponge_field(t, tv%T, g, gv, nz, csp)
693 call set_up_sponge_field(s, tv%S, g, gv, nz, csp)
694
695 endif
696
697end subroutine isomip_initialize_sponges
698
699!> \namespace isomip_initialization
700!!
701!! See this paper for details: http://www.geosci-model-dev-discuss.net/8/9859/2015/gmdd-8-9859-2015.pdf
702end module isomip_initialization