24 implicit none ;
private
26 #include <MOM_memory.h>
28 character(len=40) ::
mdl =
"ISOMIP_initialization"
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
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
153 real :: min_thickness, s_sur, s_bot, t_sur, t_bot
154 real :: rho_sur, rho_bot
157 character(len=256) :: mesg
158 character(len=40) :: verticalcoordinate
160 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
162 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
164 if (.not.just_read) &
165 call mom_mesg(
"MOM_initialization.F90, initialize_thickness_uniform: setting thickness")
167 call get_param(param_file,
mdl,
"MIN_THICKNESS", min_thickness, &
168 'Minimum layer thickness', units=
'm', default=1.e-3, do_not_log=just_read, scale=us%m_to_Z)
169 call get_param(param_file,
mdl,
"REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
170 default=default_coordinate_mode, do_not_log=just_read)
172 select case ( coordinatemode(verticalcoordinate) )
174 case ( regridding_layer, regridding_rho )
176 'Temperature at the surface (interface)', default=-1.9, do_not_log=just_read)
177 call get_param(param_file,
mdl,
"ISOMIP_S_SUR", s_sur, &
178 'Salinity at the surface (interface)', default=33.8, do_not_log=just_read)
179 call get_param(param_file,
mdl,
"ISOMIP_T_BOT", t_bot, &
180 'Temperature at the bottom (interface)', default=-1.9, do_not_log=just_read)
182 'Salinity at the bottom (interface)', default=34.55, do_not_log=just_read)
184 if (just_read)
return
187 call calculate_density(t_sur, s_sur, 0.0, rho_sur, tv%eqn_of_state, scale=us%kg_m3_to_R)
190 call calculate_density(t_bot, s_bot, 0.0, rho_bot, tv%eqn_of_state, scale=us%kg_m3_to_R)
193 rho_range = rho_bot - rho_sur
200 e0(k) = -g%max_depth * ( 0.5 * ( gv%Rlay(k-1) + gv%Rlay(k) ) - rho_sur ) / rho_range
201 e0(k) = min( 0., e0(k) )
202 e0(k) = max( -g%max_depth, e0(k) )
207 e0(nz+1) = -g%max_depth
210 do j=js,je ;
do i=is,ie
211 eta1d(nz+1) = -g%bathyT(i,j)
214 if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z))
then
215 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
216 h(i,j,k) = gv%Angstrom_H
218 h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
224 if (just_read)
return
225 do j=js,je ;
do i=is,ie
226 eta1d(nz+1) = -g%bathyT(i,j)
228 eta1d(k) = -g%max_depth * real(k-1) / real(nz)
229 if (eta1d(k) < (eta1d(k+1) + min_thickness))
then
230 eta1d(k) = eta1d(k+1) + min_thickness
231 h(i,j,k) = gv%Z_to_H * min_thickness
233 h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
238 case ( regridding_sigma )
239 if (just_read)
return
240 do j=js,je ;
do i=is,ie
241 h(i,j,:) = gv%Z_to_H * g%bathyT(i,j) / dfloat(nz)
245 call mom_error(fatal,
"isomip_initialize: "// &
246 "Unrecognized i.c. setup - set REGRIDDING_COORDINATE_MODE")
254 eqn_of_state, just_read_params)
258 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(out) :: t
259 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(out) :: s
260 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(in) :: h
262 type(
eos_type),
pointer :: eqn_of_state
263 logical,
optional,
intent(in) :: just_read_params
266 integer :: i, j, k, is, ie, js, je, nz, itt
268 real :: rho_sur, rho_bot
275 character(len=256) :: mesg
276 character(len=40) :: verticalcoordinate, density_profile
280 real :: t0(szk_(g)), s0(szk_(g))
281 real :: drho_dt(szk_(g))
282 real :: drho_ds(szk_(g))
283 real :: rho_guess(szk_(g))
284 real :: pres(szk_(g))
288 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
291 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
293 call get_param(param_file,
mdl,
"REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
294 default=default_coordinate_mode, do_not_log=just_read)
296 'Temperature at the surface (interface)', default=-1.9, do_not_log=just_read)
297 call get_param(param_file,
mdl,
"ISOMIP_S_SUR", s_sur, &
298 'Salinity at the surface (interface)', default=33.8, do_not_log=just_read)
299 call get_param(param_file,
mdl,
"ISOMIP_T_BOT", t_bot, &
300 'Temperature at the bottom (interface)', default=-1.9, do_not_log=just_read)
301 call get_param(param_file,
mdl,
"ISOMIP_S_BOT", s_bot, &
302 'Salinity at the bottom (interface)', default=34.55, do_not_log=just_read)
304 call calculate_density(t_sur,s_sur,0.0,rho_sur,eqn_of_state, scale=us%kg_m3_to_R)
307 call calculate_density(t_bot,s_bot,0.0,rho_bot,eqn_of_state, scale=us%kg_m3_to_R)
311 select case ( coordinatemode(verticalcoordinate) )
314 if (just_read)
return
316 ds_dz = (s_sur - s_bot) / g%max_depth
317 dt_dz = (t_sur - t_bot) / g%max_depth
318 do j=js,je ;
do i=is,ie
321 xi0 = xi0 + 0.5 * h(i,j,k) * gv%H_to_Z
322 s(i,j,k) = s_sur + ds_dz * xi0
323 t(i,j,k) = t_sur + dt_dz * xi0
324 xi0 = xi0 + 0.5 * h(i,j,k) * gv%H_to_Z
328 case ( regridding_layer )
329 call get_param(param_file,
mdl,
"FIT_SALINITY", fit_salin, &
330 "If true, accept the prescribed temperature and fit the "//&
331 "salinity; otherwise take salinity and fit temperature.", &
332 default=.false., do_not_log=just_read)
334 "Partial derivative of density with salinity.", &
335 units=
"kg m-3 PSU-1", scale=us%kg_m3_to_R, fail_if_missing=.not.just_read, do_not_log=just_read)
337 "Partial derivative of density with temperature.", &
338 units=
"kg m-3 K-1", scale=us%kg_m3_to_R, fail_if_missing=.not.just_read, do_not_log=just_read)
340 "A reference temperature used in initialization.", &
341 units=
"degC", fail_if_missing=.not.just_read, do_not_log=just_read)
343 "A reference salinity used in initialization.", units=
"PSU", &
344 default=35.0, do_not_log=just_read)
345 if (just_read)
return
350 ds_dz = (s_sur - s_bot) / g%max_depth
351 dt_dz = (t_sur - t_bot) / g%max_depth
353 do j=js,je ;
do i=is,ie
357 xi1 = xi0 + 0.5 * h(i,j,k) * gv%H_to_Z
358 s0(k) = s_sur - ds_dz * xi1
359 t0(k) = t_sur - dt_dz * xi1
360 xi0 = xi0 + h(i,j,k) * gv%H_to_Z
368 call calculate_density(t0(1),s0(1),0.,rho_guess(1),eqn_of_state, scale=us%kg_m3_to_R)
373 s0(k) = max(0.0, s0(1) + (gv%Rlay(k) - rho_guess(1)) / drho_ds1)
377 call calculate_density(t0,s0,pres,rho_guess,1,nz,eqn_of_state, scale=us%kg_m3_to_R)
380 s0(k) = max(0.0, s0(k) + (gv%Rlay(k) - rho_guess(k)) / drho_ds1)
387 t0(k) = t0(1) + (gv%Rlay(k) - rho_guess(1)) / drho_dt1
391 call calculate_density(t0,s0,pres,rho_guess,1,nz,eqn_of_state, scale=us%kg_m3_to_R)
394 t0(k) = t0(k) + (gv%Rlay(k) - rho_guess(k)) / drho_dt(k)
400 t(i,j,k) = t0(k) ; s(i,j,k) = s0(k)
406 call mom_error(fatal,
"isomip_initialize: "// &
407 "Unrecognized i.c. setup - set REGRIDDING_COORDINATE_MODE")
436 logical,
intent(in) :: use_ale
440 real :: t(szi_(g),szj_(g),szk_(g))
441 real :: s(szi_(g),szj_(g),szk_(g))
442 real :: rho(szi_(g),szj_(g),szk_(g))
443 real :: h(szi_(g),szj_(g),szk_(g))
444 real :: idamp(szi_(g),szj_(g))
449 real :: rho_sur, rho_bot
453 real :: e0(szk_(g)+1)
455 real :: eta1d(szk_(g)+1)
456 real :: eta(szi_(g),szj_(g),szk_(g)+1)
457 real :: min_depth, dummy1, z
458 real :: damp, rho_dummy, min_thickness, rho_tmp, xi0
459 character(len=40) :: verticalcoordinate, filename, state_file
460 character(len=40) :: temp_var, salt_var, eta_var, inputdir
462 character(len=40) ::
mdl =
"ISOMIP_initialize_sponges"
463 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
465 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
466 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
468 call get_param(pf,
mdl,
"MIN_THICKNESS", min_thickness,
"Minimum layer thickness", &
469 units=
"m", default=1.e-3, scale=us%m_to_Z)
471 call get_param(pf,
mdl,
"REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
472 default=default_coordinate_mode)
474 call get_param(pf,
mdl,
"ISOMIP_TNUDG", tnudg,
"Nudging time scale for sponge layers (days)", default=0.0)
476 call get_param(pf,
mdl,
"T_REF", t_ref,
"Reference temperature", default=10.0,&
479 call get_param(pf,
mdl,
"S_REF", s_ref,
"Reference salinity", default=35.0,&
483 'Surface salinity in sponge layer.', default=s_ref)
486 'Bottom salinity in sponge layer.', default=s_ref)
489 'Surface temperature in sponge layer.', default=t_ref)
492 'Bottom temperature in sponge layer.', default=t_ref)
494 t(:,:,:) = 0.0 ; s(:,:,:) = 0.0 ; idamp(:,:) = 0.0; rho(:,:,:) = 0.0
498 "The minimum depth of the ocean.", units=
"m", default=0.0, scale=us%m_to_Z)
500 if (
associated(csp))
call mom_error(fatal, &
501 "ISOMIP_initialize_sponges called with an associated control structure.")
502 if (
associated(acsp))
call mom_error(fatal, &
503 "ISOMIP_initialize_sponges called with an associated ALE-sponge control structure.")
510 do i=is,ie;
do j=js,je
511 if (g%geoLonT(i,j) >= 790.0 .AND. g%geoLonT(i,j) <= 800.0)
then
514 dummy1=(g%geoLonT(i,j)-790.0)/(800.0-790.0)
515 damp = 1.0/tnudg * max(0.0,dummy1)
521 if (g%bathyT(i,j) > min_depth)
then
522 idamp(i,j) = damp/86400.0
523 else ; idamp(i,j) = 0.0 ;
endif
528 call calculate_density(t_sur, s_sur, 0.0, rho_sur, tv%eqn_of_state, scale=us%kg_m3_to_R)
531 call calculate_density(t_bot, s_bot, 0.0, rho_bot, tv%eqn_of_state, scale=us%kg_m3_to_R)
534 rho_range = rho_bot - rho_sur
540 select case ( coordinatemode(verticalcoordinate) )
542 case ( regridding_rho )
546 e0(k) = -g%max_depth * ( 0.5 * ( gv%Rlay(k-1) + gv%Rlay(k) ) - rho_sur ) / rho_range
547 e0(k) = min( 0., e0(k) )
548 e0(k) = max( -g%max_depth, e0(k) )
553 e0(nz+1) = -g%max_depth
556 do j=js,je ;
do i=is,ie
557 eta1d(nz+1) = -g%bathyT(i,j)
560 if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z))
then
561 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
562 h(i,j,k) = gv%Angstrom_H
564 h(i,j,k) = gv%Z_to_H*(eta1d(k) - eta1d(k+1))
570 do j=js,je ;
do i=is,ie
571 eta1d(nz+1) = -g%bathyT(i,j)
573 eta1d(k) = -g%max_depth * real(k-1) / real(nz)
574 if (eta1d(k) < (eta1d(k+1) + min_thickness))
then
575 eta1d(k) = eta1d(k+1) + min_thickness
576 h(i,j,k) = min_thickness * gv%Z_to_H
578 h(i,j,k) = gv%Z_to_H*(eta1d(k) - eta1d(k+1))
583 case ( regridding_sigma )
584 do j=js,je ;
do i=is,ie
585 h(i,j,:) = gv%Z_to_H * (g%bathyT(i,j) / dfloat(nz))
589 call mom_error(fatal,
"ISOMIP_initialize_sponges: "// &
590 "Unrecognized i.c. setup - set REGRIDDING_COORDINATE_MODE")
598 ds_dz = (s_sur - s_bot) / g%max_depth
599 dt_dz = (t_sur - t_bot) / g%max_depth
600 do j=js,je ;
do i=is,ie
603 xi0 = xi0 + 0.5 * h(i,j,k) * gv%H_to_Z
604 s(i,j,k) = s_sur + ds_dz * xi0
605 t(i,j,k) = t_sur + dt_dz * xi0
606 xi0 = xi0 + 0.5 * h(i,j,k) * gv%H_to_Z
622 if (
associated(tv%T) )
then
625 if (
associated(tv%S) )
then
631 call get_param(pf,
mdl,
"INPUTDIR", inputdir, default=
".")
632 inputdir = slasher(inputdir)
638 call get_param(pf,
mdl,
"ISOMIP_SPONGE_FILE", state_file, &
639 "The name of the file with temps., salts. and interfaces to "//&
640 "damp toward.", fail_if_missing=.true.)
642 "The name of the potential temperature variable in "//&
643 "SPONGE_STATE_FILE.", default=
"Temp")
645 "The name of the salinity variable in "//&
646 "SPONGE_STATE_FILE.", default=
"Salt")
648 "The name of the interface height variable in "//&
649 "SPONGE_STATE_FILE.", default=
"eta")
652 filename = trim(inputdir)//trim(state_file)
654 "ISOMIP_initialize_sponges: Unable to open "//trim(filename))
655 call mom_read_data(filename, eta_var, eta(:,:,:), g%Domain, scale=us%m_to_Z)