21 implicit none ;
private
23 #include <MOM_memory.h>
36 character(len=40) ::
mdl =
"DOME2D_initialization"
43 real,
dimension(G%isd:G%ied,G%jsd:G%jed), &
46 real,
intent(in) :: max_depth
50 real :: x, bay_depth, l1, l2
51 real :: dome2d_width_bay, dome2d_width_bottom, dome2d_depth_bay
53 # include "version_variable.h"
56 call get_param(param_file,
mdl,
"DOME2D_SHELF_WIDTH", dome2d_width_bay, &
57 'Width of shelf, as fraction of domain, in 2d DOME configuration.', &
58 units=
'nondim',default=0.1)
59 call get_param(param_file,
mdl,
"DOME2D_BASIN_WIDTH", dome2d_width_bottom, &
60 'Width of deep ocean basin, as fraction of domain, in 2d DOME configuration.', &
61 units=
'nondim',default=0.3)
62 call get_param(param_file,
mdl,
"DOME2D_SHELF_DEPTH", dome2d_depth_bay, &
63 'Depth of shelf, as fraction of basin depth, in 2d DOME configuration.', &
64 units=
'nondim',default=0.2)
70 l2 = 1.0 - dome2d_width_bottom
72 bay_depth = dome2d_depth_bay
74 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
77 x = ( g%geoLonT(i,j) - g%west_lon ) / g%len_lon
80 d(i,j) = bay_depth * max_depth
81 elseif (( x > l1 ) .and. ( x < l2 ))
then
82 d(i,j) = bay_depth * max_depth + (1.0-bay_depth) * max_depth * &
83 ( x - l1 ) / (l2 - l1)
97 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
101 logical,
optional,
intent(in) :: just_read_params
107 real :: eta1d(szk_(gv)+1)
109 integer :: i, j, k, is, ie, js, je, nz
112 real :: min_thickness
113 real :: dome2d_width_bay, dome2d_width_bottom, dome2d_depth_bay
115 character(len=40) :: verticalcoordinate
117 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
119 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
121 if (.not.just_read) &
122 call mom_mesg(
"MOM_initialization.F90, DOME2d_initialize_thickness: setting thickness")
124 call get_param(param_file,
mdl,
"MIN_THICKNESS",min_thickness, &
125 default=1.e-3, units=
"m", do_not_log=.true., scale=us%m_to_Z)
126 call get_param(param_file,
mdl,
"REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
127 default=default_coordinate_mode, do_not_log=.true.)
128 call get_param(param_file,
mdl,
"DOME2D_SHELF_WIDTH", dome2d_width_bay, &
129 default=0.1, do_not_log=.true.)
130 call get_param(param_file,
mdl,
"DOME2D_BASIN_WIDTH", dome2d_width_bottom, &
131 default=0.3, do_not_log=.true.)
132 call get_param(param_file,
mdl,
"DOME2D_SHELF_DEPTH", dome2d_depth_bay, &
133 default=0.2, do_not_log=.true.)
135 if (just_read)
return
145 e0(k) = -g%max_depth * real(k-1) / real(nz)
148 select case ( coordinatemode(verticalcoordinate) )
152 do j=js,je ;
do i=is,ie
153 eta1d(nz+1) = -g%bathyT(i,j)
156 if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z))
then
157 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
158 h(i,j,k) = gv%Angstrom_H
160 h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
164 x = ( g%geoLonT(i,j) - g%west_lon ) / g%len_lon
165 if ( x <= dome2d_width_bay )
then
166 h(i,j,1:nz-1) = gv%Angstrom_H
167 h(i,j,nz) = gv%Z_to_H * dome2d_depth_bay * g%max_depth - (nz-1) * gv%Angstrom_H
194 case ( regridding_zstar )
196 do j=js,je ;
do i=is,ie
197 eta1d(nz+1) = -g%bathyT(i,j)
200 if (eta1d(k) < (eta1d(k+1) + min_thickness))
then
201 eta1d(k) = eta1d(k+1) + min_thickness
202 h(i,j,k) = gv%Z_to_H * min_thickness
204 h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
210 do j=js,je ;
do i=is,ie
211 h(i,j,:) = gv%Z_to_H*g%bathyT(i,j) / nz
215 call mom_error(fatal,
"dome2d_initialize: "// &
216 "Unrecognized i.c. setup - set REGRIDDING_COORDINATE_MODE")
225 eqn_of_state, just_read_params)
228 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(out) :: t
229 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(out) :: s
230 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(in) :: h
232 type(
eos_type),
pointer :: eqn_of_state
233 logical,
optional,
intent(in) :: just_read_params
237 integer :: i, j, k, is, ie, js, je, nz
239 integer :: index_bay_z
240 real :: delta_s, delta_t
241 real :: s_ref, t_ref;
242 real :: s_range, t_range;
245 character(len=40) :: verticalcoordinate
246 real :: dome2d_width_bay, dome2d_width_bottom, dome2d_depth_bay
248 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
250 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
252 call get_param(param_file,
mdl,
"REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
253 default=default_coordinate_mode, do_not_log=.true.)
254 call get_param(param_file,
mdl,
"DOME2D_SHELF_WIDTH", dome2d_width_bay, &
255 default=0.1, do_not_log=.true.)
256 call get_param(param_file,
mdl,
"DOME2D_BASIN_WIDTH", dome2d_width_bottom, &
257 default=0.3, do_not_log=.true.)
258 call get_param(param_file,
mdl,
"DOME2D_SHELF_DEPTH", dome2d_depth_bay, &
259 default=0.2, do_not_log=.true.)
260 call get_param(param_file,
mdl,
"S_REF",s_ref,
'Reference salinity',units=
'1e-3', &
261 fail_if_missing=.not.just_read, do_not_log=just_read)
262 call get_param(param_file,
mdl,
"T_REF",t_ref,
'Refernce temperature',units=
'C', &
263 fail_if_missing=.not.just_read, do_not_log=just_read)
264 call get_param(param_file,
mdl,
"S_RANGE",s_range,
'Initial salinity range', &
265 units=
'1e-3', default=2.0, do_not_log=just_read)
266 call get_param(param_file,
mdl,
"T_RANGE",t_range,
'Initial temperature range', &
267 units=
'1e-3', default=0.0, do_not_log=just_read)
269 if (just_read)
return
276 select case ( coordinatemode(verticalcoordinate) )
280 do j=js,je ;
do i=is,ie
283 xi1 = xi0 + (gv%H_to_Z * h(i,j,k)) / g%max_depth
284 s(i,j,k) = 34.0 + 0.5 * s_range * (xi0 + xi1)
291 do j=js,je ;
do i=is,ie
294 xi1 = xi0 + (gv%H_to_Z * h(i,j,k)) / g%max_depth
295 s(i,j,k) = 34.0 + 0.5 * s_range * (xi0 + xi1)
298 x = ( g%geoLonT(i,j) - g%west_lon ) / g%len_lon
299 if ( x <= dome2d_width_bay )
then
300 s(i,j,nz) = 34.0 + s_range
304 case ( regridding_layer )
306 delta_s = s_range / ( g%ke - 1.0 )
309 s(:,:,k) = s(:,:,k-1) + delta_s
313 call mom_error(fatal,
"dome2d_initialize: "// &
314 "Unrecognized i.c. setup - set REGRIDDING_COORDINATE_MODE")
319 if ( coordinatemode(verticalcoordinate) == regridding_zstar )
then
320 index_bay_z = nint( dome2d_depth_bay * g%ke )
321 do j = g%jsc,g%jec ;
do i = g%isc,g%iec
322 x = ( g%geoLonT(i,j) - g%west_lon ) / g%len_lon
323 if ( x <= dome2d_width_bay )
then
324 s(i,j,1:index_bay_z) = s_ref + s_range;
325 t(i,j,1:index_bay_z) = 1.0;
332 do i = g%isc,g%iec ;
do j = g%jsc,g%jec
333 x = ( g%geoLonT(i,j) - g%west_lon ) / g%len_lon
334 if ( x <= dome2d_width_bay )
then
335 s(i,j,1:g%ke) = s_ref + s_range;
342 t(g%isc:g%iec,g%jsc:g%jec,1:g%ke) = 0.0
343 if (( coordinatemode(verticalcoordinate) ==
regridding_rho ) .or. &
344 ( coordinatemode(verticalcoordinate) == regridding_layer ))
then
345 do i = g%isc,g%iec ;
do j = g%jsc,g%jec
346 x = ( g%geoLonT(i,j) - g%west_lon ) / g%len_lon
347 if ( x <= dome2d_width_bay )
then
361 logical,
intent(in) :: use_ale
365 real :: t(szi_(g),szj_(g),szk_(g))
366 real :: s(szi_(g),szj_(g),szk_(g))
367 real :: rho(szi_(g),szj_(g),szk_(g))
368 real :: h(szi_(g),szj_(g),szk_(g))
369 real :: eta(szi_(g),szj_(g),szk_(g)+1)
370 real :: idamp(szi_(g),szj_(g))
372 real :: s_range, t_range
373 real :: e0(szk_(g)+1)
375 real :: eta1d(szk_(g)+1)
377 real :: d_eta(szk_(g))
378 real :: dome2d_width_bay, dome2d_width_bottom, dome2d_depth_bay
379 real :: dome2d_west_sponge_time_scale, dome2d_east_sponge_time_scale
380 real :: dome2d_west_sponge_width, dome2d_east_sponge_width
382 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
384 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
385 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
387 call get_param(param_file,
mdl,
"DOME2D_WEST_SPONGE_TIME_SCALE", dome2d_west_sponge_time_scale, &
388 'The time-scale on the west edge of the domain for restoring T/S '//&
389 'in the sponge. If zero, the western sponge is disabled', &
390 units=
's', default=0.)
391 call get_param(param_file,
mdl,
"DOME2D_EAST_SPONGE_TIME_SCALE", dome2d_east_sponge_time_scale, &
392 'The time-scale on the east edge of the domain for restoring T/S '//&
393 'in the sponge. If zero, the eastern sponge is disabled', &
394 units=
's', default=0.)
395 call get_param(param_file,
mdl,
"DOME2D_WEST_SPONGE_WIDTH", dome2d_west_sponge_width, &
396 'The fraction of the domain in which the western sponge for restoring T/S '//&
398 units=
'nondim', default=0.1)
399 call get_param(param_file,
mdl,
"DOME2D_EAST_SPONGE_WIDTH", dome2d_east_sponge_width, &
400 'The fraction of the domain in which the eastern sponge for restoring T/S '//&
402 units=
'nondim', default=0.1)
405 if (dome2d_west_sponge_time_scale <= 0. .and. dome2d_east_sponge_time_scale <= 0.)
return
407 if (
associated(csp))
call mom_error(fatal, &
408 "DOME2d_initialize_sponges called with an associated control structure.")
409 if (
associated(acsp))
call mom_error(fatal, &
410 "DOME2d_initialize_sponges called with an associated ALE-sponge control structure.")
412 call get_param(param_file,
mdl,
"DOME2D_SHELF_WIDTH", dome2d_width_bay, &
413 default=0.1, do_not_log=.true.)
414 call get_param(param_file,
mdl,
"DOME2D_BASIN_WIDTH", dome2d_width_bottom, &
415 default=0.3, do_not_log=.true.)
416 call get_param(param_file,
mdl,
"DOME2D_SHELF_DEPTH", dome2d_depth_bay, &
417 default=0.2, do_not_log=.true.)
420 call get_param(param_file,
mdl,
"S_RANGE",s_range,default=2.0)
421 call get_param(param_file,
mdl,
"T_RANGE",t_range,default=0.0)
426 do j=js,je ;
do i=is,ie
427 if (g%mask2dT(i,j) > 0.)
then
428 x = ( g%geoLonT(i,j) - g%west_lon ) / g%len_lon
429 if ( dome2d_west_sponge_time_scale > 0. .and. x < dome2d_west_sponge_width )
then
431 dummy1 = 1. - x / dome2d_west_sponge_width
432 idamp(i,j) = 1./dome2d_west_sponge_time_scale * max(0., min(1., dummy1))
433 elseif ( dome2d_east_sponge_time_scale > 0. .and. x > ( 1. - dome2d_east_sponge_width ) )
then
435 dummy1 = 1. - ( 1. - x ) / dome2d_east_sponge_width
436 idamp(i,j) = 1./dome2d_east_sponge_time_scale * max(0., min(1., dummy1))
450 e0(k) = -g%max_depth * ( real(k-1) / real(nz) )
452 e0(nz+1) = -g%max_depth
453 do j=js,je ;
do i=is,ie
454 eta1d(nz+1) = -g%bathyT(i,j)
457 if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z))
then
458 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
459 h(i,j,k) = gv%Angstrom_H
461 h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
469 t(:,:,:) = 0.0 ; s(:,:,:) = 0.0
470 do j=js,je ;
do i=is,ie
473 z = z + 0.5 * gv%H_to_Z * h(i,j,k)
474 s(i,j,k) = 34.0 - 1.0 * (z / (g%max_depth))
475 if ( ( g%geoLonT(i,j) - g%west_lon ) / g%len_lon < dome2d_west_sponge_width ) &
476 s(i,j,k) = s_ref + s_range
477 z = z + 0.5 * gv%H_to_Z * h(i,j,k)
481 if (
associated(tv%T) )
then
484 if (
associated(tv%S) )
then
491 do j=js,je ;
do i=is,ie
492 eta1d(nz+1) = -g%bathyT(i,j)
494 eta1d(k) = -g%max_depth * real(k-1) / real(nz)
495 if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z))
then
496 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
497 d_eta(k) = gv%Angstrom_Z
499 d_eta(k) = (eta1d(k) - eta1d(k+1))
503 x = ( g%geoLonT(i,j) - g%west_lon ) / g%len_lon
504 if ( x <= dome2d_width_bay )
then
505 do k=1,nz-1 ; d_eta(k) = gv%Angstrom_Z ;
enddo
506 d_eta(nz) = dome2d_depth_bay * g%max_depth - (nz-1) * gv%Angstrom_Z
509 eta(i,j,nz+1) = -g%bathyT(i,j)
511 eta(i,j,k) = eta(i,j,k+1) + d_eta(k)