20 use regrid_consts,
only : coordinatemode, default_coordinate_mode
24 implicit none ;
private
26 #include <MOM_memory.h>
28 character(len=40) :: mdl =
"ISOMIP_initialization"
31 public isomip_initialize_topography
32 public isomip_initialize_thickness
33 public isomip_initialize_temperature_salinity
34 public isomip_initialize_sponges
44 subroutine isomip_initialize_topography(D, G, param_file, max_depth, US)
46 real,
dimension(G%isd:G%ied,G%jsd:G%jed), &
49 real,
intent(in) :: max_depth
67 #include "version_variable.h"
68 character(len=40) :: mdl =
"ISOMIP_initialize_topography"
69 integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
70 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
71 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
73 call mom_mesg(
" ISOMIP_initialization.F90, ISOMIP_initialize_topography: setting topography", 5)
75 m_to_z = 1.0 ;
if (
present(us)) m_to_z = us%m_to_Z
78 call get_param(param_file, mdl,
"MINIMUM_DEPTH", min_depth, &
79 "The minimum depth of the ocean.", units=
"m", default=0.0, scale=m_to_z)
80 call get_param(param_file, mdl,
"ISOMIP_2D",is_2d,
'If true, use a 2D setup.', default=.false.)
83 bmax = 720.0*m_to_z ; dc = 500.0*m_to_z
84 b0 = -150.0*m_to_z ; b2 = -728.8*m_to_z ; b4 = 343.91*m_to_z ; b6 = -50.57*m_to_z
85 xbar = 300.0e3 ; fc = 4.0e3 ; wc = 24.0e3 ; ly = 80.0e3
86 bx = 0.0 ; by = 0.0 ; xtil = 0.0
90 do j=js,je ;
do i=is,ie
92 xtil = g%geoLonT(i,j)*1.0e3/xbar
94 bx = b0 + b2*xtil**2 + b4*xtil**4 + b6*xtil**6
99 by = (dc / (1.+exp(-2.*(40.0*1.0e3- ly/2. - wc)/fc))) + &
100 (dc / (1.+exp(2.*(40.0*1.0e3- ly/2. + wc)/fc)))
102 d(i,j) = -max(bx+by, -bmax)
103 if (d(i,j) > max_depth) d(i,j) = max_depth
104 if (d(i,j) < min_depth) d(i,j) = 0.5*min_depth
108 do j=js,je ;
do i=is,ie
118 xtil = g%geoLonT(i,j)*1.0e3/xbar
120 bx = b0 + b2*xtil**2 + b4*xtil**4 + b6*xtil**6
121 by = (dc / (1.+exp(-2.*(g%geoLatT(i,j)*1.0e3- ly/2. - wc)/fc))) + &
122 (dc / (1.+exp(2.*(g%geoLatT(i,j)*1.0e3- ly/2. + wc)/fc)))
124 d(i,j) = -max(bx+by, -bmax)
125 if (d(i,j) > max_depth) d(i,j) = max_depth
126 if (d(i,j) < min_depth) d(i,j) = 0.5*min_depth
130 end subroutine isomip_initialize_topography
133 subroutine isomip_initialize_thickness ( h, G, GV, US, param_file, tv, just_read_params)
137 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
144 logical,
optional,
intent(in) :: just_read_params
147 real :: e0(szk_(g)+1)
149 real :: eta1d(szk_(g)+1)
151 integer :: i, j, k, is, ie, js, je, nz, tmp1
154 real :: min_thickness, s_sur, s_bot, t_sur, t_bot, rho_sur, rho_bot
156 character(len=256) :: mesg
157 character(len=40) :: verticalcoordinate
159 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
161 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
163 if (.not.just_read) &
164 call mom_mesg(
"MOM_initialization.F90, initialize_thickness_uniform: setting thickness")
166 call get_param(param_file, mdl,
"MIN_THICKNESS", min_thickness, &
167 'Minimum layer thickness', units=
'm', default=1.e-3, do_not_log=just_read, scale=us%m_to_Z)
168 call get_param(param_file, mdl,
"REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
169 default=default_coordinate_mode, do_not_log=just_read)
171 select case ( coordinatemode(verticalcoordinate) )
173 case ( regridding_layer, regridding_rho )
174 call get_param(param_file, mdl,
"ISOMIP_T_SUR",t_sur, &
175 'Temperature at the surface (interface)', default=-1.9, do_not_log=just_read)
176 call get_param(param_file, mdl,
"ISOMIP_S_SUR", s_sur, &
177 'Salinity at the surface (interface)', default=33.8, do_not_log=just_read)
178 call get_param(param_file, mdl,
"ISOMIP_T_BOT", t_bot, &
179 'Temperature at the bottom (interface)', default=-1.9, do_not_log=just_read)
180 call get_param(param_file, mdl,
"ISOMIP_S_BOT", s_bot,&
181 'Salinity at the bottom (interface)', default=34.55, do_not_log=just_read)
183 if (just_read)
return
192 rho_range = rho_bot - rho_sur
199 e0(k) = -g%max_depth * ( 0.5 * ( gv%Rlay(k-1) + gv%Rlay(k) ) - rho_sur ) / rho_range
200 e0(k) = min( 0., e0(k) )
201 e0(k) = max( -g%max_depth, e0(k) )
205 e0(nz+1) = -g%max_depth
208 do j=js,je ;
do i=is,ie
209 eta1d(nz+1) = -g%bathyT(i,j)
212 if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z))
then
213 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
214 h(i,j,k) = gv%Angstrom_H
216 h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
221 case ( regridding_zstar, regridding_sigma_shelf_zstar )
222 if (just_read)
return
223 do j=js,je ;
do i=is,ie
224 eta1d(nz+1) = -g%bathyT(i,j)
226 eta1d(k) = -g%max_depth * real(k-1) / real(nz)
227 if (eta1d(k) < (eta1d(k+1) + min_thickness))
then
228 eta1d(k) = eta1d(k+1) + min_thickness
229 h(i,j,k) = gv%Z_to_H * min_thickness
231 h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
236 case ( regridding_sigma )
237 if (just_read)
return
238 do j=js,je ;
do i=is,ie
239 h(i,j,:) = gv%Z_to_H * g%bathyT(i,j) / dfloat(nz)
243 call mom_error(fatal,
"isomip_initialize: "// &
244 "Unrecognized i.c. setup - set REGRIDDING_COORDINATE_MODE")
248 end subroutine isomip_initialize_thickness
251 subroutine isomip_initialize_temperature_salinity ( T, S, h, G, GV, param_file, &
252 eqn_of_state, just_read_params)
255 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(out) :: t
256 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(out) :: s
257 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(in) :: h
259 type(
eos_type),
pointer :: eqn_of_state
260 logical,
optional,
intent(in) :: just_read_params
263 integer :: i, j, k, is, ie, js, je, nz, itt
264 real :: x, ds, dt, rho_sur, rho_bot
271 character(len=256) :: mesg
272 character(len=40) :: verticalcoordinate, density_profile
276 real :: t0(szk_(g)), s0(szk_(g))
277 real :: drho_dt(szk_(g))
278 real :: drho_ds(szk_(g))
279 real :: rho_guess(szk_(g))
280 real :: pres(szk_(g))
281 real :: drho_dt1, drho_ds1, t_ref, s_ref
282 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
285 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
287 call get_param(param_file, mdl,
"REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
288 default=default_coordinate_mode, do_not_log=just_read)
289 call get_param(param_file, mdl,
"ISOMIP_T_SUR",t_sur, &
290 'Temperature at the surface (interface)', default=-1.9, do_not_log=just_read)
291 call get_param(param_file, mdl,
"ISOMIP_S_SUR", s_sur, &
292 'Salinity at the surface (interface)', default=33.8, do_not_log=just_read)
293 call get_param(param_file, mdl,
"ISOMIP_T_BOT", t_bot, &
294 'Temperature at the bottom (interface)', default=-1.9, do_not_log=just_read)
295 call get_param(param_file, mdl,
"ISOMIP_S_BOT", s_bot, &
296 'Salinity at the bottom (interface)', default=34.55, do_not_log=just_read)
305 select case ( coordinatemode(verticalcoordinate) )
307 case ( regridding_rho, regridding_zstar, regridding_sigma_shelf_zstar, regridding_sigma )
308 if (just_read)
return
310 ds_dz = (s_sur - s_bot) / g%max_depth
311 dt_dz = (t_sur - t_bot) / g%max_depth
312 do j=js,je ;
do i=is,ie
315 xi0 = xi0 + 0.5 * h(i,j,k) * gv%H_to_Z
316 s(i,j,k) = s_sur + ds_dz * xi0
317 t(i,j,k) = t_sur + dt_dz * xi0
318 xi0 = xi0 + 0.5 * h(i,j,k) * gv%H_to_Z
322 case ( regridding_layer )
323 call get_param(param_file, mdl,
"FIT_SALINITY", fit_salin, &
324 "If true, accept the prescribed temperature and fit the "//&
325 "salinity; otherwise take salinity and fit temperature.", &
326 default=.false., do_not_log=just_read)
327 call get_param(param_file, mdl,
"DRHO_DS", drho_ds1, &
328 "Partial derivative of density with salinity.", &
329 units=
"kg m-3 PSU-1", fail_if_missing=.not.just_read, do_not_log=just_read)
330 call get_param(param_file, mdl,
"DRHO_DT", drho_dt1, &
331 "Partial derivative of density with temperature.", &
332 units=
"kg m-3 K-1", fail_if_missing=.not.just_read, do_not_log=just_read)
333 call get_param(param_file, mdl,
"T_REF", t_ref, &
334 "A reference temperature used in initialization.", &
335 units=
"degC", fail_if_missing=.not.just_read, do_not_log=just_read)
336 call get_param(param_file, mdl,
"S_REF", s_ref, &
337 "A reference salinity used in initialization.", units=
"PSU", &
338 default=35.0, do_not_log=just_read)
339 if (just_read)
return
344 ds_dz = (s_sur - s_bot) / g%max_depth
345 dt_dz = (t_sur - t_bot) / g%max_depth
347 do j=js,je ;
do i=is,ie
351 xi1 = xi0 + 0.5 * h(i,j,k) * gv%H_to_Z
352 s0(k) = s_sur - ds_dz * xi1
353 t0(k) = t_sur - dt_dz * xi1
354 xi0 = xi0 + h(i,j,k) * gv%H_to_Z
367 s0(k) = max(0.0, s0(1) + (gv%Rlay(k) - rho_guess(1)) / drho_ds1)
374 s0(k) = max(0.0, s0(k) + (gv%Rlay(k) - rho_guess(k)) / drho_ds1)
381 t0(k) = t0(1) + (gv%Rlay(k) - rho_guess(1)) / drho_dt1
388 t0(k) = t0(k) + (gv%Rlay(k) - rho_guess(k)) / drho_dt(k)
394 t(i,j,k) = t0(k) ; s(i,j,k) = s0(k)
400 call mom_error(fatal,
"isomip_initialize: "// &
401 "Unrecognized i.c. setup - set REGRIDDING_COORDINATE_MODE")
413 end subroutine isomip_initialize_temperature_salinity
418 subroutine isomip_initialize_sponges(G, GV, US, tv, PF, use_ALE, CSp, ACSp)
430 logical,
intent(in) :: use_ale
434 real :: t(szi_(g),szj_(g),szk_(g))
435 real :: s(szi_(g),szj_(g),szk_(g))
436 real :: rho(szi_(g),szj_(g),szk_(g))
437 real :: h(szi_(g),szj_(g),szk_(g))
438 real :: idamp(szi_(g),szj_(g))
443 real :: rho_sur, rho_bot, rho_range
446 real :: e0(szk_(g)+1)
448 real :: eta1d(szk_(g)+1)
449 real :: eta(szi_(g),szj_(g),szk_(g)+1)
450 real :: min_depth, dummy1, z
451 real :: damp, rho_dummy, min_thickness, rho_tmp, xi0
452 character(len=40) :: verticalcoordinate, filename, state_file
453 character(len=40) :: temp_var, salt_var, eta_var, inputdir
455 character(len=40) :: mdl =
"ISOMIP_initialize_sponges"
456 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
458 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
459 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
461 call get_param(pf, mdl,
"MIN_THICKNESS", min_thickness,
"Minimum layer thickness", &
462 units=
"m", default=1.e-3, scale=us%m_to_Z)
464 call get_param(pf, mdl,
"REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
465 default=default_coordinate_mode)
467 call get_param(pf, mdl,
"ISOMIP_TNUDG", tnudg,
"Nudging time scale for sponge layers (days)", default=0.0)
469 call get_param(pf, mdl,
"T_REF", t_ref,
"Reference temperature", default=10.0,&
472 call get_param(pf, mdl,
"S_REF", s_ref,
"Reference salinity", default=35.0,&
475 call get_param(pf, mdl,
"ISOMIP_S_SUR_SPONGE", s_sur, &
476 'Surface salinity in sponge layer.', default=s_ref)
478 call get_param(pf, mdl,
"ISOMIP_S_BOT_SPONGE", s_bot, &
479 'Bottom salinity in sponge layer.', default=s_ref)
481 call get_param(pf, mdl,
"ISOMIP_T_SUR_SPONGE", t_sur, &
482 'Surface temperature in sponge layer.', default=t_ref)
484 call get_param(pf, mdl,
"ISOMIP_T_BOT_SPONGE", t_bot, &
485 'Bottom temperature in sponge layer.', default=t_ref)
487 t(:,:,:) = 0.0 ; s(:,:,:) = 0.0 ; idamp(:,:) = 0.0; rho(:,:,:) = 0.0
490 call get_param(pf, mdl,
"MINIMUM_DEPTH", min_depth, &
491 "The minimum depth of the ocean.", units=
"m", default=0.0, scale=us%m_to_Z)
493 if (
associated(csp))
call mom_error(fatal, &
494 "ISOMIP_initialize_sponges called with an associated control structure.")
495 if (
associated(acsp))
call mom_error(fatal, &
496 "ISOMIP_initialize_sponges called with an associated ALE-sponge control structure.")
503 do i=is,ie;
do j=js,je
504 if (g%geoLonT(i,j) >= 790.0 .AND. g%geoLonT(i,j) <= 800.0)
then
507 dummy1=(g%geoLonT(i,j)-790.0)/(800.0-790.0)
508 damp = 1.0/tnudg * max(0.0,dummy1)
514 if (g%bathyT(i,j) > min_depth)
then
515 idamp(i,j) = damp/86400.0
516 else ; idamp(i,j) = 0.0 ;
endif
527 rho_range = rho_bot - rho_sur
533 select case ( coordinatemode(verticalcoordinate) )
535 case ( regridding_rho )
539 e0(k) = -g%max_depth * ( 0.5 * ( gv%Rlay(k-1) + gv%Rlay(k) ) - rho_sur ) / rho_range
540 e0(k) = min( 0., e0(k) )
541 e0(k) = max( -g%max_depth, e0(k) )
545 e0(nz+1) = -g%max_depth
548 do j=js,je ;
do i=is,ie
549 eta1d(nz+1) = -g%bathyT(i,j)
552 if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z))
then
553 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
554 h(i,j,k) = gv%Angstrom_H
556 h(i,j,k) = gv%Z_to_H*(eta1d(k) - eta1d(k+1))
561 case ( regridding_zstar, regridding_sigma_shelf_zstar )
562 do j=js,je ;
do i=is,ie
563 eta1d(nz+1) = -g%bathyT(i,j)
565 eta1d(k) = -g%max_depth * real(k-1) / real(nz)
566 if (eta1d(k) < (eta1d(k+1) + min_thickness))
then
567 eta1d(k) = eta1d(k+1) + min_thickness
568 h(i,j,k) = min_thickness * gv%Z_to_H
570 h(i,j,k) = gv%Z_to_H*(eta1d(k) - eta1d(k+1))
575 case ( regridding_sigma )
576 do j=js,je ;
do i=is,ie
577 h(i,j,:) = gv%Z_to_H * (g%bathyT(i,j) / dfloat(nz))
581 call mom_error(fatal,
"ISOMIP_initialize_sponges: "// &
582 "Unrecognized i.c. setup - set REGRIDDING_COORDINATE_MODE")
590 ds_dz = (s_sur - s_bot) / g%max_depth
591 dt_dz = (t_sur - t_bot) / g%max_depth
592 do j=js,je ;
do i=is,ie
595 xi0 = xi0 + 0.5 * h(i,j,k) * gv%H_to_Z
596 s(i,j,k) = s_sur + ds_dz * xi0
597 t(i,j,k) = t_sur + dt_dz * xi0
598 xi0 = xi0 + 0.5 * h(i,j,k) * gv%H_to_Z
614 if (
associated(tv%T) )
then
617 if (
associated(tv%S) )
then
623 call get_param(pf, mdl,
"INPUTDIR", inputdir, default=
".")
624 inputdir = slasher(inputdir)
630 call get_param(pf, mdl,
"ISOMIP_SPONGE_FILE", state_file, &
631 "The name of the file with temps., salts. and interfaces to "//&
632 "damp toward.", fail_if_missing=.true.)
633 call get_param(pf, mdl,
"SPONGE_PTEMP_VAR", temp_var, &
634 "The name of the potential temperature variable in "//&
635 "SPONGE_STATE_FILE.", default=
"Temp")
636 call get_param(pf, mdl,
"SPONGE_SALT_VAR", salt_var, &
637 "The name of the salinity variable in "//&
638 "SPONGE_STATE_FILE.", default=
"Salt")
639 call get_param(pf, mdl,
"SPONGE_ETA_VAR", eta_var, &
640 "The name of the interface height variable in "//&
641 "SPONGE_STATE_FILE.", default=
"eta")
644 filename = trim(inputdir)//trim(state_file)
645 if (.not.
file_exists(filename, g%Domain))
call mom_error(fatal, &
646 "ISOMIP_initialize_sponges: Unable to open "//trim(filename))
647 call mom_read_data(filename, eta_var, eta(:,:,:), g%Domain, scale=us%m_to_Z)
662 call initialize_sponge(idamp, eta, g, pf, csp, gv)
664 call set_up_sponge_field(t, tv%T, g, nz, csp)
665 call set_up_sponge_field(s, tv%S, g, nz, csp)
669 end subroutine isomip_initialize_sponges