MOM_coord_initialization.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Initializes fixed aspects of the related to its vertical coordinate.
7
8use mom_debugging, only : chksum
9use mom_eos, only : calculate_density, eos_type
10use mom_error_handler, only : mom_mesg, mom_error, fatal, warning, is_root_pe
12use mom_file_parser, only : get_param, read_param, log_param, param_file_type, log_version
13use mom_io, only : create_mom_file, file_exists
14use mom_io, only : mom_netcdf_file, mom_field
15use mom_io, only : mom_read_data, mom_write_field, vardesc, var_desc, single_file
22
23implicit none ; private
24
26
27! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
28! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
29! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
30! vary with the Boussinesq approximation, the Boussinesq variant is given first.
31
32character(len=40) :: mdl = "MOM_coord_initialization" !< This module's name.
33
34contains
35
36!> MOM_initialize_coord sets up time-invariant quantities related to MOM6's
37!! vertical coordinate.
38subroutine mom_initialize_coord(GV, US, PF, tv, max_depth)
39 type(verticalgrid_type), intent(inout) :: gv !< Ocean vertical grid structure.
40 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
41 type(param_file_type), intent(in) :: pf !< A structure indicating the open file
42 !! to parse for model parameter values.
43 type(thermo_var_ptrs), intent(inout) :: tv !< The thermodynamic variable structure.
44 real, intent(in) :: max_depth !< The ocean's maximum depth [Z ~> m].
45 ! Local
46 character(len=200) :: config
47 logical :: debug
48 ! This include declares and sets the variable "version".
49# include "version_variable.h"
50 integer :: nz
51
52 nz = gv%ke
53
54 call calltree_enter("MOM_initialize_coord(), MOM_coord_initialization.F90")
55 call log_version(pf, mdl, version, "")
56 call get_param(pf, mdl, "DEBUG", debug, default=.false.)
57
58! Set-up the layer densities, GV%Rlay, and reduced gravities, GV%g_prime.
59 call get_param(pf, mdl, "COORD_CONFIG", config, &
60 "This specifies how layers are to be defined: \n"//&
61 " \t ALE or none - used to avoid defining layers in ALE mode \n"//&
62 " \t file - read coordinate information from the file \n"//&
63 " \t\t specified by (COORD_FILE).\n"//&
64 " \t BFB - Custom coords for buoyancy-forced basin case \n"//&
65 " \t\t based on SST_S, T_BOT and DRHO_DT.\n"//&
66 " \t linear - linear based on interfaces not layers \n"//&
67 " \t layer_ref - linear based on layer densities \n"//&
68 " \t ts_ref - use reference temperature and salinity \n"//&
69 " \t ts_range - use range of temperature and salinity \n"//&
70 " \t\t (T_REF and S_REF) to determine surface density \n"//&
71 " \t\t and GINT calculate internal densities. \n"//&
72 " \t gprime - use reference density (RHO_0) for surface \n"//&
73 " \t\t density and GINT calculate internal densities. \n"//&
74 " \t ts_profile - use temperature and salinity profiles \n"//&
75 " \t\t (read from COORD_FILE) to set layer densities. \n"//&
76 " \t USER - call a user modified routine.", &
77 default="none")
78 select case ( trim(config) )
79 case ("gprime")
80 call set_coord_from_gprime(gv%Rlay, gv%g_prime, gv, us, pf)
81 case ("layer_ref")
82 call set_coord_from_layer_density(gv%Rlay, gv%g_prime, gv, us, pf)
83 case ("linear")
84 call set_coord_linear(gv%Rlay, gv%g_prime, gv, us, pf)
85 case ("ts_ref")
86 call set_coord_from_ts_ref(gv%Rlay, gv%g_prime, gv, us, pf, tv%eqn_of_state, tv%P_Ref)
87 case ("ts_profile")
88 call set_coord_from_ts_profile(gv%Rlay, gv%g_prime, gv, us, pf, tv%eqn_of_state, tv%P_Ref)
89 case ("ts_range")
90 call set_coord_from_ts_range(gv%Rlay, gv%g_prime, gv, us, pf, tv%eqn_of_state, tv%P_Ref)
91 case ("file")
92 call set_coord_from_file(gv%Rlay, gv%g_prime, gv, us, pf)
93 case ("USER")
94 call user_set_coord(gv%Rlay, gv%g_prime, gv, us, pf)
95 case ("BFB")
96 call bfb_set_coord(gv%Rlay, gv%g_prime, gv, us, pf)
97 case ("none", "ALE")
98 call set_coord_to_none(gv%Rlay, gv%g_prime, gv, us, pf)
99 case default ; call mom_error(fatal,"MOM_initialize_coord: "// &
100 "Unrecognized coordinate setup"//trim(config))
101 end select
102 ! There are nz+1 values of g_prime because it is an interface field, but the value at the bottom
103 ! should not matter. This is here just to avoid having an uninitialized value in some output.
104 gv%g_prime(nz+1) = 10.0*gv%g_Earth
105
106 if (debug) call chksum(us%R_to_kg_m3*gv%Rlay(:), "MOM_initialize_coord: Rlay ", 1, nz)
107 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)
108 call setverticalgridaxes( gv%Rlay, gv, scale=us%R_to_kg_m3 )
109
110 ! Copy the maximum depth across from the input argument
111 gv%max_depth = max_depth
112
113 call calltree_leave('MOM_initialize_coord()')
114
115end subroutine mom_initialize_coord
116
117! The set_coord routines deal with initializing aspects of the vertical grid.
118
119!> Sets the layer densities (Rlay) and the interface reduced gravities (g).
120subroutine set_coord_from_gprime(Rlay, g_prime, GV, US, param_file)
121 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
122 real, dimension(GV%ke), intent(out) :: Rlay !< The layers' target coordinate values
123 !! (potential density) [R ~> kg m-3].
124 real, dimension(GV%ke+1), intent(out) :: g_prime !< The reduced gravity across the interfaces
125 !! [L2 Z-1 T-2 ~> m s-2].
126 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
127 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
128 ! Local variables
129 real :: g_int ! Reduced gravities across the internal interfaces [L2 Z-1 T-2 ~> m s-2].
130 real :: g_fs ! Reduced gravity across the free surface [L2 Z-1 T-2 ~> m s-2].
131 real :: Rlay_Ref ! The target density of the surface layer [R ~> kg m-3].
132 character(len=40) :: mdl = "set_coord_from_gprime" ! This subroutine's name.
133 integer :: k, nz
134 nz = gv%ke
135
136 call calltree_enter(trim(mdl)//"(), MOM_coord_initialization.F90")
137
138 call get_param(param_file, mdl, "GFS" , g_fs, &
139 "The reduced gravity at the free surface.", units="m s-2", &
140 default=gv%g_Earth*us%L_T_to_m_s**2*us%m_to_Z, scale=us%m_s_to_L_T**2*us%Z_to_m)
141 call get_param(param_file, mdl, "GINT", g_int, &
142 "The reduced gravity across internal interfaces.", &
143 units="m s-2", fail_if_missing=.true., scale=us%m_s_to_L_T**2*us%Z_to_m)
144 call get_param(param_file, mdl, "LIGHTEST_DENSITY", rlay_ref, &
145 "The reference potential density used for layer 1.", &
146 units="kg m-3", default=us%R_to_kg_m3*gv%Rho0, scale=us%kg_m3_to_R)
147
148 g_prime(1) = g_fs
149 do k=2,nz ; g_prime(k) = g_int ; enddo
150 rlay(1) = rlay_ref
151 if (gv%Boussinesq .or. gv%semi_Boussinesq) then
152 do k=2,nz ; rlay(k) = rlay(k-1) + g_prime(k)*(gv%Rho0/gv%g_Earth) ; enddo
153 else
154 do k=2,nz
155 rlay(k) = rlay(k-1) * ((gv%g_Earth + 0.5*g_prime(k)) / (gv%g_Earth - 0.5*g_prime(k)))
156 enddo
157 endif
158
159 call calltree_leave(trim(mdl)//'()')
160
161end subroutine set_coord_from_gprime
162
163!> Sets the layer densities (Rlay) and the interface reduced gravities (g).
164subroutine set_coord_from_layer_density(Rlay, g_prime, GV, US, param_file)
165 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
166 real, dimension(GV%ke), intent(out) :: Rlay !< The layers' target coordinate values
167 !! (potential density) [R ~> kg m-3].
168 real, dimension(GV%ke+1), intent(out) :: g_prime !< The reduced gravity across the interfaces
169 !! [L2 Z-1 T-2 ~> m s-2].
170 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
171 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
172
173 ! Local variables
174 real :: g_fs ! Reduced gravity across the free surface [L2 Z-1 T-2 ~> m s-2].
175 real :: Rlay_Ref! The surface layer's target density [R ~> kg m-3].
176 real :: RLay_range ! The range of densities [R ~> kg m-3].
177 character(len=40) :: mdl = "set_coord_from_layer_density" ! This subroutine's name.
178 integer :: k, nz
179 nz = gv%ke
180
181 call calltree_enter(trim(mdl)//"(), MOM_coord_initialization.F90")
182
183 call get_param(param_file, mdl, "GFS", g_fs, &
184 "The reduced gravity at the free surface.", units="m s-2", &
185 default=gv%g_Earth*us%L_T_to_m_s**2*us%m_to_Z, scale=us%m_s_to_L_T**2*us%Z_to_m)
186 call get_param(param_file, mdl, "LIGHTEST_DENSITY", rlay_ref, &
187 "The reference potential density used for layer 1.", &
188 units="kg m-3", default=us%R_to_kg_m3*gv%Rho0, scale=us%kg_m3_to_R)
189 call get_param(param_file, mdl, "DENSITY_RANGE", rlay_range, &
190 "The range of reference potential densities in the layers.", &
191 units="kg m-3", default=2.0, scale=us%kg_m3_to_R)
192
193 rlay(1) = rlay_ref
194 do k=2,nz
195 rlay(k) = rlay(k-1) + rlay_range/(real(nz-1))
196 enddo
197! These statements set the interface reduced gravities. !
198 g_prime(1) = g_fs
199 if (gv%Boussinesq .or. gv%semi_Boussinesq) then
200 do k=2,nz
201 g_prime(k) = (gv%g_Earth/gv%Rho0) * (rlay(k) - rlay(k-1))
202 enddo
203 else
204 do k=2,nz
205 g_prime(k) = 2.0*gv%g_Earth * (rlay(k) - rlay(k-1)) / (rlay(k) + rlay(k-1))
206 enddo
207 endif
208
209 call calltree_leave(trim(mdl)//'()')
210end subroutine set_coord_from_layer_density
211
212!> Sets the layer densities (Rlay) and the interface reduced gravities (g) from a profile of g'.
213subroutine set_coord_from_ts_ref(Rlay, g_prime, GV, US, param_file, eqn_of_state, P_Ref)
214 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
215 real, dimension(GV%ke), intent(out) :: Rlay !< The layers' target coordinate values
216 !! (potential density) [R ~> kg m-3].
217 real, dimension(GV%ke+1), intent(out) :: g_prime !< The reduced gravity across the interfaces
218 !! [L2 Z-1 T-2 ~> m s-2].
219 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
220 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
221 type(eos_type), intent(in) :: eqn_of_state !< Equation of state structure
222 real, intent(in) :: P_Ref !< The coordinate-density reference pressure
223 !! [R L2 T-2 ~> Pa].
224
225 ! Local variables
226 real :: T_ref ! Reference temperature [C ~> degC]
227 real :: S_ref ! Reference salinity [S ~> ppt]
228 real :: g_int ! Reduced gravities across the internal interfaces [L2 Z-1 T-2 ~> m s-2].
229 real :: g_fs ! Reduced gravity across the free surface [L2 Z-1 T-2 ~> m s-2].
230 character(len=40) :: mdl = "set_coord_from_TS_ref" ! This subroutine's name.
231 integer :: k, nz
232 nz = gv%ke
233
234 call calltree_enter(trim(mdl)//"(), MOM_coord_initialization.F90")
235
236 call get_param(param_file, mdl, "T_REF", t_ref, &
237 "The initial temperature of the lightest layer.", &
238 units="degC", scale=us%degC_to_C, fail_if_missing=.true.)
239 call get_param(param_file, mdl, "S_REF", s_ref, &
240 "The initial salinities.", units="ppt", default=35.0, scale=us%ppt_to_S)
241 call get_param(param_file, mdl, "GFS", g_fs, &
242 "The reduced gravity at the free surface.", units="m s-2", &
243 default=gv%g_Earth*us%L_T_to_m_s**2*us%m_to_Z, scale=us%m_s_to_L_T**2*us%Z_to_m)
244 call get_param(param_file, mdl, "GINT", g_int, &
245 "The reduced gravity across internal interfaces.", &
246 units="m s-2", fail_if_missing=.true., scale=us%m_s_to_L_T**2*us%Z_to_m)
247
248! These statements set the interface reduced gravities. !
249 g_prime(1) = g_fs
250 do k=2,nz ; g_prime(k) = g_int ; enddo
251
252! The uppermost layer's density is set here. Subsequent layers' !
253! densities are determined from this value and the g values. !
254! T0 = 28.228 ; S0 = 34.5848 ; Pref = P_Ref
255 call calculate_density(t_ref, s_ref, p_ref, rlay(1), eqn_of_state)
256
257! These statements set the layer densities. !
258 if (gv%Boussinesq .or. gv%semi_Boussinesq) then
259 do k=2,nz ; rlay(k) = rlay(k-1) + g_prime(k)*(gv%Rho0/gv%g_Earth) ; enddo
260 else
261 do k=2,nz
262 rlay(k) = rlay(k-1) * ((gv%g_Earth + 0.5*g_prime(k)) / (gv%g_Earth - 0.5*g_prime(k)))
263 enddo
264 endif
265
266 call calltree_leave(trim(mdl)//'()')
267end subroutine set_coord_from_ts_ref
268
269!> Sets the layer densities (Rlay) and the interface reduced gravities (g) from a T-S profile.
270subroutine set_coord_from_ts_profile(Rlay, g_prime, GV, US, param_file, eqn_of_state, P_Ref)
271 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
272 real, dimension(GV%ke), intent(out) :: Rlay !< Layer potential density [R ~> kg m-3].
273 real, dimension(GV%ke+1), intent(out) :: g_prime !< The reduced gravity at each
274 !! interface [L2 Z-1 T-2 ~> m s-2].
275 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
276 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
277 type(eos_type), intent(in) :: eqn_of_state !< Equation of state structure
278 real, intent(in) :: P_Ref !< The coordinate-density reference pressure
279 !! [R L2 T-2 ~> Pa].
280
281 ! Local variables
282 real, dimension(GV%ke) :: T0 ! A profile of temperatures [C ~> degC]
283 real, dimension(GV%ke) :: S0 ! A profile of salinities [S ~> ppt]
284 real, dimension(GV%ke) :: Pref ! A array of reference pressures [R L2 T-2 ~> Pa]
285 real :: g_fs ! Reduced gravity across the free surface [L2 Z-1 T-2 ~> m s-2].
286 integer :: k, nz
287 character(len=40) :: mdl = "set_coord_from_TS_profile" ! This subroutine's name.
288 character(len=200) :: filename, coord_file, inputdir ! Strings for file/path
289 character(len=64) :: temp_var, salt_var ! Temperature and salinity names in files
290
291 nz = gv%ke
292
293 call calltree_enter(trim(mdl)//"(), MOM_coord_initialization.F90")
294
295 call get_param(param_file, mdl, "GFS", g_fs, &
296 "The reduced gravity at the free surface.", units="m s-2", &
297 default=gv%g_Earth*us%L_T_to_m_s**2*us%m_to_Z, scale=us%m_s_to_L_T**2*us%Z_to_m)
298 call get_param(param_file, mdl, "COORD_FILE", coord_file, &
299 "The file from which the coordinate temperatures and salinities are read.", &
300 fail_if_missing=.true.)
301 call get_param(param_file, mdl, "TEMP_COORD_VAR", temp_var, &
302 "The coordinate reference profile variable name for potential temperature.", &
303 default="PTEMP")
304 call get_param(param_file, mdl, "SALT_COORD_VAR", salt_var, &
305 "The coordinate reference profile variable name for salinity.", &
306 default="SALT")
307
308 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
309 filename = trim(slasher(inputdir))//trim(coord_file)
310 call log_param(param_file, mdl, "INPUTDIR/COORD_FILE", filename)
311
312 call mom_read_data(filename, temp_var, t0(:), scale=us%degC_to_C)
313 call mom_read_data(filename, salt_var, s0(:), scale=us%ppt_to_S)
314
315 if (.not.file_exists(filename)) call mom_error(fatal, &
316 " set_coord_from_TS_profile: Unable to open " //trim(filename))
317! These statements set the interface reduced gravities. !
318 g_prime(1) = g_fs
319 do k=1,nz ; pref(k) = p_ref ; enddo
320 call calculate_density(t0, s0, pref, rlay, eqn_of_state, (/1,nz/) )
321 if (gv%Boussinesq .or. gv%semi_Boussinesq) then
322 do k=2,nz
323 g_prime(k) = (gv%g_Earth/gv%Rho0) * (rlay(k) - rlay(k-1))
324 enddo
325 else
326 do k=2,nz
327 g_prime(k) = 2.0*gv%g_Earth * (rlay(k) - rlay(k-1)) / (rlay(k) + rlay(k-1))
328 enddo
329 endif
330
331 call calltree_leave(trim(mdl)//'()')
332end subroutine set_coord_from_ts_profile
333
334!> Sets the layer densities (Rlay) and the interface reduced gravities (g) from a linear T-S profile.
335subroutine set_coord_from_ts_range(Rlay, g_prime, GV, US, param_file, eqn_of_state, P_Ref)
336 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
337 real, dimension(GV%ke), intent(out) :: Rlay !< Layer potential density [R ~> kg m-3].
338 real, dimension(GV%ke+1), intent(out) :: g_prime !< The reduced gravity at each
339 !! interface [L2 Z-1 T-2 ~> m s-2].
340 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
341 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
342 type(eos_type), intent(in) :: eqn_of_state !< Equation of state structure
343 real, intent(in) :: P_Ref !< The coordinate-density reference pressure
344 !! [R L2 T-2 ~> Pa].
345
346 ! Local variables
347 real, dimension(GV%ke) :: T0 ! A profile of temperatures [C ~> degC]
348 real, dimension(GV%ke) :: S0 ! A profile of salinities [S ~> ppt]
349 real, dimension(GV%ke) :: Pref ! A array of reference pressures [R L2 T-2 ~> Pa]
350 real :: S_Ref ! Default salinity range parameters [S ~> ppt].
351 real :: T_Ref ! Default temperature range parameters [C ~> degC].
352 real :: S_Light, S_Dense ! Salinity range parameters [S ~> ppt].
353 real :: T_Light, T_Dense ! Temperature range parameters [C ~> degC].
354 real :: res_rat ! The ratio of density space resolution in the denser part
355 ! of the range to that in the lighter part of the range.
356 ! Setting this greater than 1 increases the resolution for
357 ! the denser water [nondim].
358 real :: g_fs ! Reduced gravity across the free surface [L2 Z-1 T-2 ~> m s-2].
359 real :: a1, frac_dense, k_frac ! Nondimensional temporary variables [nondim]
360 integer :: k, nz, k_light
361 character(len=40) :: mdl = "set_coord_from_TS_range" ! This subroutine's name.
362
363 nz = gv%ke
364
365 call calltree_enter(trim(mdl)//"(), MOM_coord_initialization.F90")
366
367 call get_param(param_file, mdl, "T_REF", t_ref, &
368 "The default initial temperatures.", &
369 units="degC", default=10.0, scale=us%degC_to_C)
370 call get_param(param_file, mdl, "TS_RANGE_T_LIGHT", t_light, &
371 "The initial temperature of the lightest layer when "//&
372 "COORD_CONFIG is set to ts_range.", &
373 units="degC", default=us%C_to_degC*t_ref, scale=us%degC_to_C)
374 call get_param(param_file, mdl, "TS_RANGE_T_DENSE", t_dense, &
375 "The initial temperature of the densest layer when "//&
376 "COORD_CONFIG is set to ts_range.", &
377 units="degC", default=us%C_to_degC*t_ref, scale=us%degC_to_C)
378
379 call get_param(param_file, mdl, "S_REF", s_ref, &
380 "The default initial salinities.", &
381 units="ppt", default=35.0, scale=us%ppt_to_S)
382 call get_param(param_file, mdl, "TS_RANGE_S_LIGHT", s_light, &
383 "The initial lightest salinities when COORD_CONFIG is set to ts_range.", &
384 units="ppt", default=us%S_to_ppt*s_ref, scale=us%ppt_to_S)
385 call get_param(param_file, mdl, "TS_RANGE_S_DENSE", s_dense, &
386 "The initial densest salinities when COORD_CONFIG is set to ts_range.", &
387 units="ppt", default=us%S_to_ppt*s_ref, scale=us%ppt_to_S)
388
389 call get_param(param_file, mdl, "TS_RANGE_RESOLN_RATIO", res_rat, &
390 "The ratio of density space resolution in the densest "//&
391 "part of the range to that in the lightest part of the "//&
392 "range when COORD_CONFIG is set to ts_range. Values "//&
393 "greater than 1 increase the resolution of the denser water.",&
394 default=1.0, units="nondim")
395
396 call get_param(param_file, mdl, "GFS", g_fs, &
397 "The reduced gravity at the free surface.", units="m s-2", &
398 default=gv%g_Earth*us%L_T_to_m_s**2*us%m_to_Z, scale=us%m_s_to_L_T**2*us%Z_to_m)
399
400 if ((gv%nk_rho_varies > 0) .and. (nz < gv%nk_rho_varies+2)) &
401 call mom_error(fatal, "set_coord_from_TS_range requires that NZ >= NKML+NKBL+2.")
402
403 k_light = gv%nk_rho_varies + 1
404
405 ! Set T0(k) to range from T_LIGHT to T_DENSE, and similarly for S0(k).
406 t0(k_light) = t_light ; s0(k_light) = s_light
407 a1 = 2.0 * res_rat / (1.0 + res_rat)
408 do k=k_light+1,nz
409 k_frac = real(k-k_light)/real(nz-k_light)
410 frac_dense = a1 * k_frac + (1.0 - a1) * k_frac**2
411 t0(k) = frac_dense * (t_dense - t_light) + t_light
412 s0(k) = frac_dense * (s_dense - s_light) + s_light
413 enddo
414
415 g_prime(1) = g_fs
416 do k=1,nz ; pref(k) = p_ref ; enddo
417 call calculate_density(t0, s0, pref, rlay, eqn_of_state, (/k_light,nz/) )
418 ! Extrapolate target densities for the variable density mixed and buffer layers.
419 do k=k_light-1,1,-1
420 rlay(k) = 2.0*rlay(k+1) - rlay(k+2)
421 enddo
422 if (gv%Boussinesq .or. gv%semi_Boussinesq) then
423 do k=2,nz
424 g_prime(k) = (gv%g_Earth/gv%Rho0) * (rlay(k) - rlay(k-1))
425 enddo
426 else
427 do k=2,nz
428 g_prime(k) = 2.0*gv%g_Earth * (rlay(k) - rlay(k-1)) / (rlay(k) + rlay(k-1))
429 enddo
430 endif
431
432 call calltree_leave(trim(mdl)//'()')
433end subroutine set_coord_from_ts_range
434
435! Sets the layer densities (Rlay) and the interface reduced gravities (g) from data in file.
436subroutine set_coord_from_file(Rlay, g_prime, GV, US, param_file)
437 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
438 real, dimension(GV%ke), intent(out) :: Rlay !< Layer potential density [R ~> kg m-3].
439 real, dimension(GV%ke+1), intent(out) :: g_prime !< The reduced gravity at each
440 !! interface [L2 Z-1 T-2 ~> m s-2].
441 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
442 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
443
444 ! Local variables
445 real :: g_fs ! Reduced gravity across the free surface [L2 Z-1 T-2 ~> m s-2].
446 integer :: k, nz
447 character(len=40) :: mdl = "set_coord_from_file" ! This subroutine's name.
448 character(len=40) :: coord_var
449 character(len=200) :: filename,coord_file,inputdir ! Strings for file/path
450 nz = gv%ke
451
452 call calltree_enter(trim(mdl)//"(), MOM_coord_initialization.F90")
453
454 call get_param(param_file, mdl, "GFS", g_fs, &
455 "The reduced gravity at the free surface.", units="m s-2", &
456 default=gv%g_Earth*us%L_T_to_m_s**2*us%m_to_Z, scale=us%m_s_to_L_T**2*us%Z_to_m)
457 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
458 inputdir = slasher(inputdir)
459 call get_param(param_file, mdl, "COORD_FILE", coord_file, &
460 "The file from which the coordinate densities are read.", &
461 fail_if_missing=.true.)
462 call get_param(param_file, mdl, "COORD_VAR", coord_var, &
463 "The variable in COORD_FILE that is to be used for the "//&
464 "coordinate densities.", default="Layer")
465 filename = trim(inputdir)//trim(coord_file)
466 call log_param(param_file, mdl, "INPUTDIR/COORD_FILE", filename)
467 if (.not.file_exists(filename)) call mom_error(fatal, &
468 " set_coord_from_file: Unable to open "//trim(filename))
469
470 call mom_read_data(filename, coord_var, rlay, scale=us%kg_m3_to_R)
471 g_prime(1) = g_fs
472 if (gv%Boussinesq .or. gv%semi_Boussinesq) then
473 do k=2,nz
474 g_prime(k) = (gv%g_Earth/gv%Rho0) * (rlay(k) - rlay(k-1))
475 enddo
476 else
477 do k=2,nz
478 g_prime(k) = 2.0*gv%g_Earth * (rlay(k) - rlay(k-1)) / (rlay(k) + rlay(k-1))
479 enddo
480 endif
481 do k=1,nz ; if (g_prime(k) <= 0.0) then
482 call mom_error(fatal, "MOM_initialization set_coord_from_file: "//&
483 "Zero or negative g_primes read from variable "//"Layer"//" in file "//&
484 trim(filename))
485 endif ; enddo
486
487 call calltree_leave(trim(mdl)//'()')
488end subroutine set_coord_from_file
489
490!> Sets the layer densities (Rlay) and the interface
491!! reduced gravities (g) according to a linear profile starting at a
492!! reference surface layer density and spanning a range of densities
493!! to the bottom defined by the parameter RLAY_RANGE
494!! (defaulting to 2.0 if not defined)
495subroutine set_coord_linear(Rlay, g_prime, GV, US, param_file)
496 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
497 real, dimension(GV%ke), intent(out) :: Rlay !< Layer potential density [R ~> kg m-3].
498 real, dimension(GV%ke+1), intent(out) :: g_prime !< The reduced gravity at each
499 !! interface [L2 Z-1 T-2 ~> m s-2].
500 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
501 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
502
503 ! Local variables
504 character(len=40) :: mdl = "set_coord_linear" ! This subroutine
505 real :: Rlay_ref, Rlay_range ! A reference density and its range [R ~> kg m-3]
506 real :: g_fs ! The reduced gravity across the free surface [L2 Z-1 T-2 ~> m s-2]
507 integer :: k, nz
508 nz = gv%ke
509
510 call calltree_enter(trim(mdl)//"(), MOM_coord_initialization.F90")
511
512 call get_param(param_file, mdl, "LIGHTEST_DENSITY", rlay_ref, &
513 "The reference potential density used for the surface interface.", &
514 units="kg m-3", default=us%R_to_kg_m3*gv%Rho0, scale=us%kg_m3_to_R)
515 call get_param(param_file, mdl, "DENSITY_RANGE", rlay_range, &
516 "The range of reference potential densities across all interfaces.", &
517 units="kg m-3", default=2.0, scale=us%kg_m3_to_R)
518 call get_param(param_file, mdl, "GFS", g_fs, &
519 "The reduced gravity at the free surface.", units="m s-2", &
520 default=gv%g_Earth*us%L_T_to_m_s**2*us%m_to_Z, scale=us%m_s_to_L_T**2*us%Z_to_m)
521
522 ! This following sets the target layer densities such that the
523 ! surface interface has density Rlay_ref and the bottom
524 ! is Rlay_range larger
525 do k=1,nz
526 rlay(k) = rlay_ref + rlay_range*((real(k)-0.5)/real(nz))
527 enddo
528 ! These statements set the interface reduced gravities.
529 g_prime(1) = g_fs
530 if (gv%Boussinesq .or. gv%semi_Boussinesq) then
531 do k=2,nz
532 g_prime(k) = (gv%g_Earth/gv%Rho0) * (rlay(k) - rlay(k-1))
533 enddo
534 else
535 do k=2,nz
536 g_prime(k) = 2.0*gv%g_Earth * (rlay(k) - rlay(k-1)) / (rlay(k) + rlay(k-1))
537 enddo
538 endif
539
540 call calltree_leave(trim(mdl)//'()')
541end subroutine set_coord_linear
542
543!> Sets Rlay to Rho0 and g_prime to zero except for the free surface.
544!! This is for use only in ALE mode where Rlay should not be used and g_prime(1) alone
545!! might be used.
546subroutine set_coord_to_none(Rlay, g_prime, GV, US, param_file)
547 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
548 real, dimension(GV%ke), intent(out) :: Rlay !< Layer potential density [R ~> kg m-3].
549 real, dimension(GV%ke+1), intent(out) :: g_prime !< The reduced gravity at each
550 !! interface [L2 Z-1 T-2 ~> m s-2].
551 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
552 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
553 ! Local variables
554 real :: g_fs ! Reduced gravity across the free surface [L2 Z-1 T-2 ~> m s-2].
555 real :: Rlay_Ref ! The target density of the surface layer [R ~> kg m-3].
556 character(len=40) :: mdl = "set_coord_to_none" ! This subroutine's name.
557 integer :: k, nz
558 nz = gv%ke
559
560 call calltree_enter(trim(mdl)//"(), MOM_coord_initialization.F90")
561
562 call get_param(param_file, mdl, "GFS" , g_fs, &
563 "The reduced gravity at the free surface.", units="m s-2", &
564 default=gv%g_Earth*us%L_T_to_m_s**2*us%m_to_Z, scale=us%m_s_to_L_T**2*us%Z_to_m)
565 call get_param(param_file, mdl, "LIGHTEST_DENSITY", rlay_ref, &
566 "The reference potential density used for layer 1.", &
567 units="kg m-3", default=us%R_to_kg_m3*gv%Rho0, scale=us%kg_m3_to_R)
568
569 g_prime(1) = g_fs
570 do k=2,nz ; g_prime(k) = 0. ; enddo
571 rlay(1) = rlay_ref
572 if (gv%Boussinesq .or. gv%semi_Boussinesq) then
573 do k=2,nz ; rlay(k) = rlay(k-1) + g_prime(k)*(gv%Rho0/gv%g_Earth) ; enddo
574 else
575 do k=2,nz
576 rlay(k) = rlay(k-1) * ((gv%g_Earth + 0.5*g_prime(k)) / (gv%g_Earth - 0.5*g_prime(k)))
577 enddo
578 endif
579
580 call calltree_leave(trim(mdl)//'()')
581
582end subroutine set_coord_to_none
583
584!> Writes out a file containing any available data related
585!! to the vertical grid used by the MOM ocean model.
586subroutine write_vertgrid_file(GV, US, param_file, directory)
587 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
588 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
589 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
590 character(len=*), intent(in) :: directory !< The directory into which to place the file.
591 ! Local variables
592 character(len=240) :: filepath
593 type(vardesc) :: vars(2)
594 type(mom_field) :: fields(2)
595 type(mom_netcdf_file) :: io_handle ! The I/O handle of the fileset
596
597 filepath = trim(directory) // trim("Vertical_coordinate.nc")
598
599 vars(1) = var_desc("R","kilogram meter-3","Target Potential Density",'1','L','1')
600 vars(2) = var_desc("g","meter second-2","Reduced gravity",'1','i','1')
601
602 call create_mom_file(io_handle, trim(filepath), vars, 2, fields, &
603 single_file, gv=gv)
604
605 call mom_write_field(io_handle, fields(1), gv%Rlay, unscale=us%R_to_kg_m3)
606 call mom_write_field(io_handle, fields(2), gv%g_prime, unscale=us%L_T_to_m_s**2*us%m_to_Z)
607
608 call io_handle%close()
609
610end subroutine write_vertgrid_file
611