24 implicit none ;
private
33 character(len=40) ::
mdl =
"MOM_coord_initialization"
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")
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) )
105 case default ;
call mom_error(fatal,
"MOM_initialize_coord: "// &
106 "Unrecognized coordinate setup"//trim(config))
108 if (debug)
call chksum(us%R_to_kg_m3*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)
113 gv%max_depth = max_depth
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
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=us%R_to_kg_m3*gv%Rho0, scale=us%kg_m3_to_R)
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, scale=us%kg_m3_to_R)
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))
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
243 call calculate_density(t_ref, s_ref, p_ref, rlay(1), eqn_of_state, scale=us%kg_m3_to_R)
246 do k=2,nz ; rlay(k) = rlay(k-1) + g_prime(k)*(gv%Rho0/gv%g_Earth) ;
enddo
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
293 call calculate_density(t0, s0, pref, rlay, 1, nz, eqn_of_state, scale=us%kg_m3_to_R)
294 do k=2,nz; g_prime(k) = (gv%g_Earth/(gv%Rho0)) * (rlay(k) - rlay(k-1)) ;
enddo
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
373 call calculate_density(t0, s0, pref, rlay, k_light, nz-k_light+1, eqn_of_state, scale=us%kg_m3_to_R)
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
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)
419 do k=1,nz ; rlay(k) = us%kg_m3_to_R*rlay(k) ;
enddo
421 do k=2,nz ; g_prime(k) = (gv%g_Earth/(gv%Rho0)) * (rlay(k) - rlay(k-1)) ;
enddo
422 do k=1,nz ;
if (g_prime(k) <= 0.0)
then
423 call mom_error(fatal,
"MOM_initialization set_coord_from_file: "//&
424 "Zero or negative g_primes read from variable "//
"Layer"//
" in file "//&
437 real,
dimension(:),
intent(out) :: Rlay
439 real,
dimension(:),
intent(out) :: g_prime
445 character(len=40) :: mdl =
"set_coord_linear"
446 real :: Rlay_ref, Rlay_range, g_fs
450 call calltree_enter(trim(mdl)//
"(), MOM_coord_initialization.F90")
452 call get_param(param_file, mdl,
"LIGHTEST_DENSITY", rlay_ref, &
453 "The reference potential density used for the surface interface.", &
454 units=
"kg m-3", default=us%R_to_kg_m3*gv%Rho0, scale=us%kg_m3_to_R)
455 call get_param(param_file, mdl,
"DENSITY_RANGE", rlay_range, &
456 "The range of reference potential densities across all interfaces.", &
457 units=
"kg m-3", default=2.0, scale=us%kg_m3_to_R)
458 call get_param(param_file, mdl,
"GFS", g_fs, &
459 "The reduced gravity at the free surface.", units=
"m s-2", &
460 default=gv%mks_g_Earth, scale=us%m_s_to_L_T**2*us%Z_to_m)
466 rlay(k) = rlay_ref + rlay_range*((real(k)-0.5)/real(nz))
471 g_prime(k) = (gv%g_Earth/(gv%Rho0)) * (rlay(k) - rlay(k-1))
481 real,
dimension(:),
intent(out) :: Rlay
483 real,
dimension(:),
intent(out) :: g_prime
490 character(len=40) :: mdl =
"set_coord_to_none"
494 call calltree_enter(trim(mdl)//
"(), MOM_coord_initialization.F90")
496 call get_param(param_file, mdl,
"GFS" , g_fs, &
497 "The reduced gravity at the free surface.", units=
"m s-2", &
498 default=gv%mks_g_Earth, scale=us%m_s_to_L_T**2*us%Z_to_m)
501 do k=2,nz ; g_prime(k) = 0. ;
enddo
503 do k=2,nz ; rlay(k) = rlay(k-1) + g_prime(k)*(gv%Rho0/gv%g_Earth) ;
enddo
515 character(len=*),
intent(in) :: directory
517 character(len=240) :: filepath
519 type(fieldtype) :: fields(2)
522 filepath = trim(directory) // trim(
"Vertical_coordinate")
524 vars(1) =
var_desc(
"R",
"kilogram meter-3",
"Target Potential Density",
'1',
'L',
'1')
525 vars(2) =
var_desc(
"g",
"meter second-2",
"Reduced gravity",
'1',
'L',
'1')
527 call create_file(unit, trim(filepath), vars, 2, fields, single_file, gv=gv)
529 call write_field(unit, fields(1), us%R_to_kg_m3*gv%Rlay(:))
530 call write_field(unit, fields(2), us%L_T_to_m_s**2*us%m_to_Z*gv%g_prime(:))
532 call close_file(unit)