24 implicit none ;
private
26 public mom_initialize_coord
33 character(len=40) :: mdl =
"MOM_coord_initialization"
39 subroutine mom_initialize_coord(GV, US, PF, write_geom, output_dir, tv, max_depth)
44 logical,
intent(in) :: write_geom
45 character(len=*),
intent(in) :: output_dir
47 real,
intent(in) :: max_depth
49 character(len=200) :: config
52 #include "version_variable.h"
54 type(
eos_type),
pointer :: eos => null()
56 if (
associated(tv%eqn_of_state)) eos => tv%eqn_of_state
60 call calltree_enter(
"MOM_initialize_coord(), MOM_coord_initialization.F90")
62 call get_param(pf, mdl,
"DEBUG", debug, default=.false.)
65 call get_param(pf, mdl,
"COORD_CONFIG", config, &
66 "This specifies how layers are to be defined: \n"//&
67 " \t ALE or none - used to avoid defining layers in ALE mode \n"//&
68 " \t file - read coordinate information from the file \n"//&
69 " \t\t specified by (COORD_FILE).\n"//&
70 " \t BFB - Custom coords for buoyancy-forced basin case \n"//&
71 " \t\t based on SST_S, T_BOT and DRHO_DT.\n"//&
72 " \t linear - linear based on interfaces not layers \n"//&
73 " \t layer_ref - linear based on layer densities \n"//&
74 " \t ts_ref - use reference temperature and salinity \n"//&
75 " \t ts_range - use range of temperature and salinity \n"//&
76 " \t\t (T_REF and S_REF) to determine surface density \n"//&
77 " \t\t and GINT calculate internal densities. \n"//&
78 " \t gprime - use reference density (RHO_0) for surface \n"//&
79 " \t\t density and GINT calculate internal densities. \n"//&
80 " \t ts_profile - use temperature and salinity profiles \n"//&
81 " \t\t (read from COORD_FILE) to set layer densities. \n"//&
82 " \t USER - call a user modified routine.", &
83 fail_if_missing=.true.)
84 select case ( trim(config) )
86 call set_coord_from_gprime(gv%Rlay, gv%g_prime, gv, us, pf)
88 call set_coord_from_layer_density(gv%Rlay, gv%g_prime, gv, us, pf)
90 call set_coord_linear(gv%Rlay, gv%g_prime, gv, us, pf)
92 call set_coord_from_ts_ref(gv%Rlay, gv%g_prime, gv, us, pf, eos, tv%P_Ref)
94 call set_coord_from_ts_profile(gv%Rlay, gv%g_prime, gv, us, pf, eos, tv%P_Ref)
96 call set_coord_from_ts_range(gv%Rlay, gv%g_prime, gv, us, pf, eos, tv%P_Ref)
98 call set_coord_from_file(gv%Rlay, gv%g_prime, gv, us, pf)
100 call user_set_coord(gv%Rlay, gv%g_prime, gv, pf, eos)
102 call bfb_set_coord(gv%Rlay, gv%g_prime, gv, pf, eos)
104 call set_coord_to_none(gv%Rlay, gv%g_prime, gv, us, pf)
105 case default ;
call mom_error(fatal,
"MOM_initialize_coord: "// &
106 "Unrecognized coordinate setup"//trim(config))
108 if (debug)
call chksum(gv%Rlay,
"MOM_initialize_coord: Rlay ", 1, nz)
109 if (debug)
call chksum(us%m_to_Z*us%L_to_m**2*us%s_to_T**2*gv%g_prime(:),
"MOM_initialize_coord: g_prime ", 1, nz)
110 call setverticalgridaxes( gv%Rlay, gv )
113 gv%max_depth = max_depth
116 if (write_geom)
call write_vertgrid_file(gv, us, pf, output_dir)
118 call calltree_leave(
'MOM_initialize_coord()')
120 end subroutine mom_initialize_coord
125 subroutine set_coord_from_gprime(Rlay, g_prime, GV, US, param_file)
126 real,
dimension(:),
intent(out) :: Rlay
128 real,
dimension(:),
intent(out) :: g_prime
136 character(len=40) :: mdl =
"set_coord_from_gprime"
140 call calltree_enter(trim(mdl)//
"(), MOM_coord_initialization.F90")
142 call get_param(param_file, mdl,
"GFS" , g_fs, &
143 "The reduced gravity at the free surface.", units=
"m s-2", &
144 default=gv%mks_g_Earth, scale=us%m_s_to_L_T**2*us%Z_to_m)
145 call get_param(param_file, mdl,
"GINT", g_int, &
146 "The reduced gravity across internal interfaces.", &
147 units=
"m s-2", fail_if_missing=.true., scale=us%m_s_to_L_T**2*us%Z_to_m)
150 do k=2,nz ; g_prime(k) = g_int ;
enddo
152 do k=2,nz ; rlay(k) = rlay(k-1) + g_prime(k)*(gv%Rho0/gv%g_Earth) ;
enddo
154 call calltree_leave(trim(mdl)//
'()')
156 end subroutine set_coord_from_gprime
159 subroutine set_coord_from_layer_density(Rlay, g_prime, GV, US, param_file)
160 real,
dimension(:),
intent(out) :: Rlay
162 real,
dimension(:),
intent(out) :: g_prime
171 character(len=40) :: mdl =
"set_coord_from_layer_density"
175 call calltree_enter(trim(mdl)//
"(), MOM_coord_initialization.F90")
177 call get_param(param_file, mdl,
"GFS", g_fs, &
178 "The reduced gravity at the free surface.", units=
"m s-2", &
179 default=gv%mks_g_Earth, scale=us%m_s_to_L_T**2*us%Z_to_m)
180 call get_param(param_file, mdl,
"LIGHTEST_DENSITY", rlay_ref, &
181 "The reference potential density used for layer 1.", &
182 units=
"kg m-3", default=gv%Rho0)
183 call get_param(param_file, mdl,
"DENSITY_RANGE", rlay_range, &
184 "The range of reference potential densities in the layers.", &
185 units=
"kg m-3", default=2.0)
190 rlay(k) = rlay(k-1) + rlay_range/(real(nz-1))
194 g_prime(k) = (gv%g_Earth/gv%Rho0) * (rlay(k) - rlay(k-1))
197 call calltree_leave(trim(mdl)//
'()')
198 end subroutine set_coord_from_layer_density
201 subroutine set_coord_from_ts_ref(Rlay, g_prime, GV, US, param_file, eqn_of_state, &
203 real,
dimension(:),
intent(out) :: Rlay
205 real,
dimension(:),
intent(out) :: g_prime
211 type(
eos_type),
pointer :: eqn_of_state
212 real,
intent(in) :: P_Ref
218 character(len=40) :: mdl =
"set_coord_from_TS_ref"
222 call calltree_enter(trim(mdl)//
"(), MOM_coord_initialization.F90")
224 call get_param(param_file, mdl,
"T_REF", t_ref, &
225 "The initial temperature of the lightest layer.", units=
"degC", &
226 fail_if_missing=.true.)
227 call get_param(param_file, mdl,
"S_REF", s_ref, &
228 "The initial salinities.", units=
"PSU", default=35.0)
229 call get_param(param_file, mdl,
"GFS", g_fs, &
230 "The reduced gravity at the free surface.", units=
"m s-2", &
231 default=gv%mks_g_Earth, scale=us%m_s_to_L_T**2*us%Z_to_m)
232 call get_param(param_file, mdl,
"GINT", g_int, &
233 "The reduced gravity across internal interfaces.", &
234 units=
"m s-2", fail_if_missing=.true., scale=us%m_s_to_L_T**2*us%Z_to_m)
238 do k=2,nz ; g_prime(k) = g_int ;
enddo
246 do k=2,nz ; rlay(k) = rlay(k-1) + g_prime(k)*(gv%Rho0/gv%g_Earth) ;
enddo
248 call calltree_leave(trim(mdl)//
'()')
249 end subroutine set_coord_from_ts_ref
252 subroutine set_coord_from_ts_profile(Rlay, g_prime, GV, US, param_file, &
254 real,
dimension(:),
intent(out) :: Rlay
256 real,
dimension(:),
intent(out) :: g_prime
262 type(
eos_type),
pointer :: eqn_of_state
263 real,
intent(in) :: P_Ref
265 real,
dimension(GV%ke) :: T0, S0, Pref
268 character(len=40) :: mdl =
"set_coord_from_TS_profile"
269 character(len=200) :: filename, coord_file, inputdir
272 call calltree_enter(trim(mdl)//
"(), MOM_coord_initialization.F90")
274 call get_param(param_file, mdl,
"GFS", g_fs, &
275 "The reduced gravity at the free surface.", units=
"m s-2", &
276 default=gv%mks_g_Earth, scale=us%m_s_to_L_T**2*us%Z_to_m)
277 call get_param(param_file, mdl,
"COORD_FILE", coord_file, &
278 "The file from which the coordinate temperatures and "//&
279 "salinities are read.", fail_if_missing=.true.)
281 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
282 filename = trim(slasher(inputdir))//trim(coord_file)
283 call log_param(param_file, mdl,
"INPUTDIR/COORD_FILE", filename)
288 if (.not.
file_exists(filename))
call mom_error(fatal, &
289 " set_coord_from_TS_profile: Unable to open " //trim(filename))
292 do k=1,nz ; pref(k) = p_ref ;
enddo
294 do k=2,nz; g_prime(k) = (gv%g_Earth/gv%Rho0) * (rlay(k) - rlay(k-1)) ;
enddo
296 call calltree_leave(trim(mdl)//
'()')
297 end subroutine set_coord_from_ts_profile
300 subroutine set_coord_from_ts_range(Rlay, g_prime, GV, US, param_file, &
302 real,
dimension(:),
intent(out) :: Rlay
304 real,
dimension(:),
intent(out) :: g_prime
310 type(
eos_type),
pointer :: eqn_of_state
311 real,
intent(in) :: P_Ref
314 real,
dimension(GV%ke) :: T0, S0, Pref
315 real :: S_Ref, S_Light, S_Dense
316 real :: T_Ref, T_Light, T_Dense
322 real :: a1, frac_dense, k_frac
323 integer :: k, nz, k_light
324 character(len=40) :: mdl =
"set_coord_from_TS_range"
325 character(len=200) :: filename, coord_file, inputdir
328 call calltree_enter(trim(mdl)//
"(), MOM_coord_initialization.F90")
330 call get_param(param_file, mdl,
"T_REF", t_ref, &
331 "The default initial temperatures.", units=
"degC", default=10.0)
332 call get_param(param_file, mdl,
"TS_RANGE_T_LIGHT", t_light, &
333 "The initial temperature of the lightest layer when "//&
334 "COORD_CONFIG is set to ts_range.", units=
"degC", default=t_ref)
335 call get_param(param_file, mdl,
"TS_RANGE_T_DENSE", t_dense, &
336 "The initial temperature of the densest layer when "//&
337 "COORD_CONFIG is set to ts_range.", units=
"degC", default=t_ref)
339 call get_param(param_file, mdl,
"S_REF", s_ref, &
340 "The default initial salinities.", units=
"PSU", default=35.0)
341 call get_param(param_file, mdl,
"TS_RANGE_S_LIGHT", s_light, &
342 "The initial lightest salinities when COORD_CONFIG "//&
343 "is set to ts_range.", default = s_ref, units=
"PSU")
344 call get_param(param_file, mdl,
"TS_RANGE_S_DENSE", s_dense, &
345 "The initial densest salinities when COORD_CONFIG "//&
346 "is set to ts_range.", default = s_ref, units=
"PSU")
348 call get_param(param_file, mdl,
"TS_RANGE_RESOLN_RATIO", res_rat, &
349 "The ratio of density space resolution in the densest "//&
350 "part of the range to that in the lightest part of the "//&
351 "range when COORD_CONFIG is set to ts_range. Values "//&
352 "greater than 1 increase the resolution of the denser water.",&
353 default=1.0, units=
"nondim")
355 call get_param(param_file, mdl,
"GFS", g_fs, &
356 "The reduced gravity at the free surface.", units=
"m s-2", &
357 default=gv%mks_g_Earth, scale=us%m_s_to_L_T**2*us%Z_to_m)
359 k_light = gv%nk_rho_varies + 1
362 t0(k_light) = t_light ; s0(k_light) = s_light
363 a1 = 2.0 * res_rat / (1.0 + res_rat)
365 k_frac = real(k-k_light)/real(nz-k_light)
366 frac_dense = a1 * k_frac + (1.0 - a1) * k_frac**2
367 t0(k) = frac_dense * (t_dense - t_light) + t_light
368 s0(k) = frac_dense * (s_dense - s_light) + s_light
372 do k=1,nz ; pref(k) = p_ref ;
enddo
376 rlay(k) = 2.0*rlay(k+1) - rlay(k+2)
378 do k=2,nz; g_prime(k) = (gv%g_Earth/gv%Rho0) * (rlay(k) - rlay(k-1));
enddo
380 call calltree_leave(trim(mdl)//
'()')
381 end subroutine set_coord_from_ts_range
384 subroutine set_coord_from_file(Rlay, g_prime, GV, US, param_file)
385 real,
dimension(:),
intent(out) :: Rlay
387 real,
dimension(:),
intent(out) :: g_prime
395 character(len=40) :: mdl =
"set_coord_from_file"
396 character(len=40) :: coord_var
397 character(len=200) :: filename,coord_file,inputdir
400 call calltree_enter(trim(mdl)//
"(), MOM_coord_initialization.F90")
402 call get_param(param_file, mdl,
"GFS", g_fs, &
403 "The reduced gravity at the free surface.", units=
"m s-2", &
404 default=gv%mks_g_Earth, scale=us%m_s_to_L_T**2*us%Z_to_m)
405 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
406 inputdir = slasher(inputdir)
407 call get_param(param_file, mdl,
"COORD_FILE", coord_file, &
408 "The file from which the coordinate densities are read.", &
409 fail_if_missing=.true.)
410 call get_param(param_file, mdl,
"COORD_VAR", coord_var, &
411 "The variable in COORD_FILE that is to be used for the "//&
412 "coordinate densities.", default=
"Layer")
413 filename = trim(inputdir)//trim(coord_file)
414 call log_param(param_file, mdl,
"INPUTDIR/COORD_FILE", filename)
415 if (.not.
file_exists(filename))
call mom_error(fatal, &
416 " set_coord_from_file: Unable to open "//trim(filename))
418 call read_axis_data(filename, coord_var, rlay)
420 do k=2,nz ; g_prime(k) = (gv%g_Earth/gv%Rho0) * (rlay(k) - rlay(k-1)) ;
enddo
421 do k=1,nz ;
if (g_prime(k) <= 0.0)
then
422 call mom_error(fatal,
"MOM_initialization set_coord_from_file: "//&
423 "Zero or negative g_primes read from variable "//
"Layer"//
" in file "//&
427 call calltree_leave(trim(mdl)//
'()')
428 end subroutine set_coord_from_file
435 subroutine set_coord_linear(Rlay, g_prime, GV, US, param_file)
436 real,
dimension(:),
intent(out) :: Rlay
438 real,
dimension(:),
intent(out) :: g_prime
444 character(len=40) :: mdl =
"set_coord_linear"
445 real :: Rlay_ref, Rlay_range, g_fs
449 call calltree_enter(trim(mdl)//
"(), MOM_coord_initialization.F90")
451 call get_param(param_file, mdl,
"LIGHTEST_DENSITY", rlay_ref, &
452 "The reference potential density used for the surface interface.", &
453 units=
"kg m-3", default=gv%Rho0)
454 call get_param(param_file, mdl,
"DENSITY_RANGE", rlay_range, &
455 "The range of reference potential densities across all interfaces.", &
456 units=
"kg m-3", default=2.0)
457 call get_param(param_file, mdl,
"GFS", g_fs, &
458 "The reduced gravity at the free surface.", units=
"m s-2", &
459 default=gv%mks_g_Earth, scale=us%m_s_to_L_T**2*us%Z_to_m)
465 rlay(k) = rlay_ref + rlay_range*((real(k)-0.5)/real(nz))
470 g_prime(k) = (gv%g_Earth/gv%Rho0) * (rlay(k) - rlay(k-1))
473 call calltree_leave(trim(mdl)//
'()')
474 end subroutine set_coord_linear
479 subroutine set_coord_to_none(Rlay, g_prime, GV, US, param_file)
480 real,
dimension(:),
intent(out) :: Rlay
482 real,
dimension(:),
intent(out) :: g_prime
489 character(len=40) :: mdl =
"set_coord_to_none"
493 call calltree_enter(trim(mdl)//
"(), MOM_coord_initialization.F90")
495 call get_param(param_file, mdl,
"GFS" , g_fs, &
496 "The reduced gravity at the free surface.", units=
"m s-2", &
497 default=gv%mks_g_Earth, scale=us%m_s_to_L_T**2*us%Z_to_m)
500 do k=2,nz ; g_prime(k) = 0. ;
enddo
502 do k=2,nz ; rlay(k) = rlay(k-1) + g_prime(k)*(gv%Rho0/gv%g_Earth) ;
enddo
504 call calltree_leave(trim(mdl)//
'()')
506 end subroutine set_coord_to_none
510 subroutine write_vertgrid_file(GV, US, param_file, directory)
514 character(len=*),
intent(in) :: directory
516 character(len=240) :: filepath
518 type(fieldtype) :: fields(2)
521 filepath = trim(directory) // trim(
"Vertical_coordinate")
523 vars(1) = var_desc(
"R",
"kilogram meter-3",
"Target Potential Density",
'1',
'L',
'1')
524 vars(2) = var_desc(
"g",
"meter second-2",
"Reduced gravity",
'1',
'L',
'1')
526 call create_file(unit, trim(filepath), vars, 2, fields, single_file, gv=gv)
528 call write_field(unit, fields(1), gv%Rlay)
529 call write_field(unit, fields(2), us%L_T_to_m_s**2*us%m_to_Z*gv%g_prime(:))
531 call close_file(unit)
533 end subroutine write_vertgrid_file