DOME_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!> Configures the model for the "DOME" experiment.
6!! DOME = Dynamics of Overflows and Mixing Experiment
8
11use mom_error_handler, only : mom_mesg, mom_error, fatal, warning, is_root_pe
12use mom_file_parser, only : get_param, log_version, param_file_type
13use mom_get_input, only : directories
14use mom_grid, only : ocean_grid_type
15use mom_open_boundary, only : ocean_obc_type, obc_none
17use mom_tracer_registry, only : tracer_registry_type, tracer_type
22use mom_eos, only : calculate_density, calculate_density_derivs
23
24implicit none ; private
25
26#include <MOM_memory.h>
27
32
33! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
34! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
35! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
36! vary with the Boussinesq approximation, the Boussinesq variant is given first.
37
38contains
39
40! -----------------------------------------------------------------------------
41!> This subroutine sets up the DOME topography
42subroutine dome_initialize_topography(D, G, param_file, max_depth, US)
43 type(dyn_horgrid_type), intent(in) :: g !< The dynamic horizontal grid type
44 real, dimension(G%isd:G%ied,G%jsd:G%jed), &
45 intent(out) :: d !< Ocean bottom depth [Z ~> m]
46 type(param_file_type), intent(in) :: param_file !< Parameter file structure
47 real, intent(in) :: max_depth !< Maximum model depth [Z ~> m]
48 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
49
50 ! Local variables
51 real :: min_depth ! The minimum ocean depth [Z ~> m]
52 real :: shelf_depth ! The ocean depth on the shelf in the DOME configuration [Z ~> m]
53 real :: slope ! The bottom slope in the DOME configuration [Z L-1 ~> nondim]
54 real :: shelf_edge_lat ! The latitude of the edge of the topographic shelf in the same units as geolat, often [km]
55 real :: inflow_lon ! The edge longitude of the DOME inflow in the same units as geolon, often [km]
56 real :: inflow_width ! The longitudinal width of the DOME inflow channel in the same units as geolat, often [km]
57 real :: km_to_grid_unit ! The conversion factor from km to the units of latitude often 1 [nondim],
58 ! but this could be 1000 [m km-1]
59 ! This include declares and sets the variable "version".
60# include "version_variable.h"
61 character(len=40) :: mdl = "DOME_initialize_topography" ! This subroutine's name.
62 integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
63 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
64 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
65
66 if (g%grid_unit_to_L <= 0.) call mom_error(fatal, "DOME_initialization: "//&
67 "DOME_initialize_topography is only set to work with Cartesian axis units.")
68 if (abs(g%grid_unit_to_L*us%L_to_m - 1000.0) < 1.0e-3) then ! The grid latitudes are in km.
69 km_to_grid_unit = 1.0
70 elseif (abs(g%grid_unit_to_L*us%L_to_m - 1.0) < 1.0e-6) then ! The grid latitudes are in m.
71 km_to_grid_unit = 1000.0
72 else
73 call mom_error(fatal, "DOME_initialization: "//&
74 "DOME_initialize_topography is not recognizing the value of G%grid_unit_to_L.")
75 endif
76
77 call mom_mesg(" DOME_initialization.F90, DOME_initialize_topography: setting topography", 5)
78
79 call log_version(param_file, mdl, version, "")
80 call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, &
81 "The minimum depth of the ocean.", default=0.0, units="m", scale=us%m_to_Z)
82 call get_param(param_file, mdl, "DOME_TOPOG_SLOPE", slope, &
83 "The slope of the bottom topography in the DOME configuration.", &
84 default=0.01, units="nondim", scale=us%L_to_Z)
85 call get_param(param_file, mdl, "DOME_SHELF_DEPTH", shelf_depth, &
86 "The bottom depth in the shelf inflow region in the DOME configuration.", &
87 default=600.0, units="m", scale=us%m_to_Z)
88 call get_param(param_file, mdl, "DOME_SHELF_EDGE_LAT", shelf_edge_lat, &
89 "The latitude of the shelf edge in the DOME configuration.", &
90 default=600.0, units="km", scale=km_to_grid_unit)
91 call get_param(param_file, mdl, "DOME_INFLOW_LON", inflow_lon, &
92 "The edge longitude of the DOME inflow.", units="km", default=1000.0, scale=km_to_grid_unit)
93 call get_param(param_file, mdl, "DOME_INFLOW_WIDTH", inflow_width, &
94 "The longitudinal width of the DOME inflow channel.", &
95 units="km", default=100.0, scale=km_to_grid_unit)
96
97 do j=js,je ; do i=is,ie
98 if (g%geoLatT(i,j) < shelf_edge_lat) then
99 d(i,j) = min(shelf_depth - slope * (g%geoLatT(i,j)-shelf_edge_lat)*g%grid_unit_to_L, max_depth)
100 else
101 if ((g%geoLonT(i,j) > inflow_lon) .AND. (g%geoLonT(i,j) < inflow_lon+inflow_width)) then
102 d(i,j) = shelf_depth
103 else
104 d(i,j) = 0.5*min_depth
105 endif
106 endif
107
108 if (d(i,j) > max_depth) d(i,j) = max_depth
109 if (d(i,j) < min_depth) d(i,j) = 0.5*min_depth
110 enddo ; enddo
111
112end subroutine dome_initialize_topography
113! -----------------------------------------------------------------------------
114
115! -----------------------------------------------------------------------------
116!> This subroutine initializes layer thicknesses for the DOME experiment
117subroutine dome_initialize_thickness(h, depth_tot, G, GV, param_file, just_read)
118 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
119 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
120 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
121 intent(out) :: h !< The thickness that is being initialized [Z ~> m]
122 real, dimension(SZI_(G),SZJ_(G)), &
123 intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
124 type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
125 !! to parse for model parameter values.
126 logical, intent(in) :: just_read !< If true, this call will only read
127 !! parameters without changing h.
128
129 real :: e0(szk_(gv)+1) ! The resting interface heights [Z ~> m], usually
130 ! negative because it is positive upward [Z ~> m].
131 real :: eta1d(szk_(gv)+1) ! Interface height relative to the sea surface
132 ! positive upward [Z ~> m].
133 integer :: i, j, k, is, ie, js, je, nz
134
135 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
136
137 if (just_read) return ! This subroutine has no run-time parameters.
138
139 call mom_mesg(" DOME_initialization.F90, DOME_initialize_thickness: setting thickness", 5)
140
141 e0(1)=0.0
142 do k=2,nz
143 e0(k) = -g%max_depth * (real(k-1)-0.5)/real(nz-1)
144 enddo
145
146 do j=g%jsc,g%jec ; do i=g%isc,g%iec
147! This sets the initial thickness (in m) of the layers. The !
148! thicknesses are set to insure that: 1. each layer is at least an !
149! Angstrom thick, and 2. the interfaces are where they should be !
150! based on the resting depths and interface height perturbations, !
151! as long at this doesn't interfere with 1. !
152 eta1d(nz+1) = -depth_tot(i,j)
153 do k=nz,1,-1
154 eta1d(k) = e0(k)
155 if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z)) then
156 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
157 h(i,j,k) = gv%Angstrom_Z
158 else
159 h(i,j,k) = eta1d(k) - eta1d(k+1)
160 endif
161 enddo
162 enddo ; enddo
163
164end subroutine dome_initialize_thickness
165! -----------------------------------------------------------------------------
166
167! -----------------------------------------------------------------------------
168!> This subroutine sets the inverse restoration time (Idamp), and the values
169!! toward which the interface heights and an arbitrary number of tracers will be
170!! restored within the sponges for the DOME configuration. !
171subroutine dome_initialize_sponges(G, GV, US, tv, depth_tot, PF, CSp)
172 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
173 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
174 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
175 type(thermo_var_ptrs), intent(in) :: tv !< A structure containing any available
176 !! thermodynamic fields, including potential
177 !! temperature and salinity or mixed layer density.
178 real, dimension(SZI_(G),SZJ_(G)), &
179 intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
180 type(param_file_type), intent(in) :: pf !< A structure indicating the open file to
181 !! parse for model parameter values.
182 type(sponge_cs), pointer :: csp !< A pointer that is set to point to the control
183 !! structure for this module.
184
185 real :: eta(szi_(g),szj_(g),szk_(gv)+1) ! A temporary array for interface heights [Z ~> m].
186 real :: temp(szi_(g),szj_(g),szk_(gv)) ! A temporary array for other variables [various]
187 real :: idamp(szi_(g),szj_(g)) ! The sponge damping rate [T-1 ~> s-1]
188
189 real :: e_tgt(szk_(gv)+1) ! Target interface heights [Z ~> m].
190 real :: min_depth ! The minimum depth at which to apply damping [Z ~> m]
191 real :: damp_w, damp_e ! Damping rates in the western and eastern sponges [T-1 ~> s-1]
192 real :: peak_damping ! The maximum sponge damping rates as the edges [T-1 ~> s-1]
193 real :: edge_dist ! The distance to an edge [L ~> m]
194 real :: sponge_width ! The width of the sponges [L ~> m]
195 character(len=40) :: mdl = "DOME_initialize_sponges" ! This subroutine's name.
196 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
197
198 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
199 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
200
201 if (g%grid_unit_to_L <= 0.) call mom_error(fatal, "DOME_initialization: "//&
202 "DOME_initialize_sponges is only set to work with Cartesian axis units.")
203
204 ! Set up sponges for the DOME configuration
205 call get_param(pf, mdl, "MINIMUM_DEPTH", min_depth, &
206 "The minimum depth of the ocean.", units="m", default=0.0, scale=us%m_to_Z)
207 call get_param(pf, mdl, "DOME_SPONGE_DAMP_RATE", peak_damping, &
208 "The largest damping rate in the DOME sponges.", &
209 default=10.0, units="day-1", scale=1.0/(86400.0*us%s_to_T))
210 call get_param(pf, mdl, "DOME_SPONGE_WIDTH", sponge_width, &
211 "The width of the DOME sponges.", &
212 default=200.0, units="km", scale=1.0e3*us%m_to_L)
213
214 ! Here the inverse damping time [T-1 ~> s-1], is set. Set Idamp to 0 wherever
215 ! there is no sponge, and the subroutines that are called will automatically
216 ! set up the sponges only where Idamp is positive and mask2dT is 1.
217
218 idamp(:,:) = 0.0
219 do j=js,je ; do i=is,ie ; if (depth_tot(i,j) > min_depth) then
220 edge_dist = (g%geoLonT(i,j) - g%west_lon) * g%grid_unit_to_L
221 if (edge_dist < 0.5*sponge_width) then
222 damp_w = peak_damping
223 elseif (edge_dist < sponge_width) then
224 damp_w = peak_damping * (sponge_width - edge_dist) / (0.5*sponge_width)
225 else
226 damp_w = 0.0
227 endif
228
229 edge_dist = ((g%len_lon + g%west_lon) - g%geoLonT(i,j)) * g%grid_unit_to_L
230 if (edge_dist < 0.5*sponge_width) then
231 damp_e = peak_damping
232 elseif (edge_dist < sponge_width) then
233 damp_e = peak_damping * (sponge_width - edge_dist) / (0.5*sponge_width)
234 else
235 damp_e = 0.0
236 endif
237
238 idamp(i,j) = max(damp_w, damp_e)
239 endif ; enddo ; enddo
240
241 e_tgt(1) = 0.0
242 do k=2,nz ; e_tgt(k) = -(real(k-1)-0.5)*g%max_depth / real(nz-1) ; enddo
243 e_tgt(nz+1) = -g%max_depth
244 eta(:,:,:) = 0.0
245 do k=1,nz+1 ; do j=js,je ; do i=is,ie
246 ! These target interface heights will be rescaled inside of apply_sponge, so
247 ! they can be in depth space for Boussinesq or non-Boussinesq models.
248 eta(i,j,k) = max(e_tgt(k), gv%Angstrom_Z*(nz+1-k) - depth_tot(i,j))
249 enddo ; enddo ; enddo
250
251 ! This call stores the sponge damping rates and target interface heights.
252 call initialize_sponge(idamp, eta, g, pf, csp, gv)
253
254 ! Now register all of the fields which are damped in the sponge.
255 ! By default, momentum is advected vertically within the sponge, but
256 ! momentum is typically not damped within the layer-mode sponge.
257
258! At this point, the layer-mode DOME configuration is done. The following are here as a
259! template for other configurations.
260
261 ! The remaining calls to set_up_sponge_field can be in any order.
262 if ( associated(tv%T) ) then
263 temp(:,:,:) = 0.0
264 call mom_error(fatal, "DOME_initialize_sponges is not set up for use with "//&
265 "temperatures defined.")
266 ! This should use the target values of T in temp.
267 call set_up_sponge_field(temp, tv%T, g, gv, nz, csp)
268 ! This should use the target values of S in temp.
269 call set_up_sponge_field(temp, tv%S, g, gv, nz, csp)
270 endif
271
272end subroutine dome_initialize_sponges
273
274!> Add DOME to the OBC registry and set up some variables that will be used to guide
275!! code setting up the restart fields related to the OBCs.
276subroutine register_dome_obc(param_file, US, OBC, tr_Reg)
277 type(param_file_type), intent(in) :: param_file !< parameter file.
278 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
279 type(ocean_obc_type), pointer :: obc !< OBC registry.
280 type(tracer_registry_type), pointer :: tr_reg !< Tracer registry.
281
282 if (obc%number_of_segments /= 1) then
283 call mom_error(fatal, 'Error in register_DOME_OBC - DOME should have 1 OBC segment', .true.)
284 endif
285
286 ! Store this information for use in setting up the OBC restarts for tracer reservoirs.
287 obc%ntr = tr_reg%ntr
288 if (.not. allocated(obc%tracer_x_reservoirs_used)) then
289 allocate(obc%tracer_x_reservoirs_used(obc%ntr))
290 allocate(obc%tracer_y_reservoirs_used(obc%ntr))
291 obc%tracer_x_reservoirs_used(:) = .false.
292 obc%tracer_y_reservoirs_used(:) = .false.
293 obc%tracer_y_reservoirs_used(1) = .true.
294 endif
295
296end subroutine register_dome_obc
297
298!> This subroutine sets the properties of flow at open boundary conditions.
299!! This particular example is for the DOME inflow describe in Legg et al. 2006.
300subroutine dome_set_obc_data(OBC, tv, G, GV, US, PF, tr_Reg)
301 type(ocean_obc_type), pointer :: obc !< This open boundary condition type specifies
302 !! whether, where, and what open boundary
303 !! conditions are used.
304 type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any
305 !! available thermodynamic fields, including potential
306 !! temperature and salinity or mixed layer density. Absent
307 !! fields have NULL ptrs.
308 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
309 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
310 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
311 type(param_file_type), intent(in) :: pf !< A structure indicating the open file
312 !! to parse for model parameter values.
313 type(tracer_registry_type), pointer :: tr_reg !< Tracer registry.
314
315 ! Local variables
316 real :: t0(szk_(gv)) ! A profile of target temperatures [C ~> degC]
317 real :: s0(szk_(gv)) ! A profile of target salinities [S ~> ppt]
318 real :: pres(szk_(gv)) ! An array of the reference pressure [R L2 T-2 ~> Pa].
319 real :: drho_dt(szk_(gv)) ! Derivative of density with temperature [R C-1 ~> kg m-3 degC-1].
320 real :: drho_ds(szk_(gv)) ! Derivative of density with salinity [R S-1 ~> kg m-3 ppt-1].
321 real :: rho_guess(szk_(gv)) ! Potential density at T0 & S0 [R ~> kg m-3].
322 real :: s_ref ! A default value for salinities [S ~> ppt]
323 real :: t_light ! A first guess at the temperature of the lightest layer [C ~> degC]
324 ! The following variables are used to set up the transport in the DOME example.
325 real :: tr_0 ! The total integrated inflow transport [H L2 T-1 ~> m3 s-1 or kg s-1]
326 real :: tr_k ! The integrated inflow transport of a layer [H L2 T-1 ~> m3 s-1 or kg s-1]
327 real :: v_k ! The velocity of a layer at the edge [L T-1 ~> m s-1]
328 real :: yt, yb ! The log of these variables gives the fractional velocities at the
329 ! top and bottom of a layer [nondim]
330 real :: rst, rsb ! The relative position of the top and bottom of a layer [nondim],
331 ! with a range from 0 for the densest water to -1 for the lightest
332 real :: rc ! The relative position of the center of a layer [nondim]
333 real :: lon_im1 ! An extrapolated value for the longitude of the western edge of a
334 ! v-velocity face, in the same units as G%geoLon [km]
335 real :: d_edge ! The thickness [Z ~> m] of the dense fluid at the
336 ! inner edge of the inflow
337 real :: rlay_range ! The range of densities [R ~> kg m-3].
338 real :: rlay_ref ! The surface layer's target density [R ~> kg m-3].
339 real :: f_0 ! The reference value of the Coriolis parameter [T-1 ~> s-1]
340 real :: f_inflow ! The value of the Coriolis parameter used to determine DOME inflow
341 ! properties [T-1 ~> s-1]
342 real :: g_prime_tot ! The reduced gravity across all layers [L2 Z-1 T-2 ~> m s-2]
343 real :: def_rad ! The deformation radius, based on fluid of thickness D_edge [L ~> m]
344 real :: inflow_lon ! The edge longitude of the DOME inflow in the same units as geolon, often [km]
345 real :: i_def_rad ! The inverse of the deformation radius in the same units as G%geoLon [km-1]
346 real :: ri_trans ! The shear Richardson number in the transition
347 ! region of the specified shear profile [nondim]
348 real :: km_to_grid_unit ! The conversion factor from km to the units of latitude often 1 [nondim],
349 ! but this could be 1000 [m km-1]
350 character(len=32) :: name ! The name of a tracer field.
351 character(len=40) :: mdl = "DOME_set_OBC_data" ! This subroutine's name.
352 integer :: i, j, k, itt, is, ie, js, je, isd, ied, jsd, jed, m, nz, ntherm, ntr_id
353 integer :: isdb, iedb, jsdb, jedb
354 type(obc_segment_type), pointer :: segment => null()
355 type(tracer_type), pointer :: tr_ptr => null()
356
357 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
358 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
359 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
360
361 if (g%grid_unit_to_L <= 0.) call mom_error(fatal, "DOME_initialization: "//&
362 "DOME_initialize_topography is only set to work with Cartesian axis units.")
363 if (abs(g%grid_unit_to_L*us%L_to_m - 1000.0) < 1.0e-3) then ! The grid latitudes are in km.
364 km_to_grid_unit = 1.0
365 elseif (abs(g%grid_unit_to_L*us%L_to_m - 1.0) < 1.0e-6) then ! The grid latitudes are in m.
366 km_to_grid_unit = 1000.0
367 else
368 call mom_error(fatal, "DOME_initialization: "//&
369 "DOME_initialize_topography is not recognizing the value of G%grid_unit_to_L.")
370 endif
371
372 call get_param(pf, mdl, "DOME_INFLOW_THICKNESS", d_edge, &
373 "The thickness of the dense DOME inflow at the inner edge.", &
374 default=300.0, units="m", scale=us%m_to_Z)
375 call get_param(pf, mdl, "DOME_INFLOW_RI_TRANS", ri_trans, &
376 "The shear Richardson number in the transition region of the specified "//&
377 "DOME inflow shear profile.", default=(1.0/3.0), units="nondim")
378 call get_param(pf, mdl, "DENSITY_RANGE", rlay_range, &
379 "The range of reference potential densities in the layers.", &
380 units="kg m-3", default=2.0, scale=us%kg_m3_to_R)
381 call get_param(pf, mdl, "LIGHTEST_DENSITY", rlay_ref, &
382 "The reference potential density used for layer 1.", &
383 units="kg m-3", default=us%R_to_kg_m3*gv%Rho0, scale=us%kg_m3_to_R)
384 call get_param(pf, mdl, "F_0", f_0, &
385 "The reference value of the Coriolis parameter with the betaplane option.", &
386 units="s-1", default=0.0, scale=us%T_to_s)
387 call get_param(pf, mdl, "DOME_INFLOW_F", f_inflow, &
388 "The value of the Coriolis parameter that is used to determine the DOME "//&
389 "inflow properties.", units="s-1", default=f_0*us%s_to_T, scale=us%T_to_s)
390 call get_param(pf, mdl, "DOME_INFLOW_LON", inflow_lon, &
391 "The edge longitude of the DOME inflow.", units="km", default=1000.0, scale=km_to_grid_unit)
392 if (associated(tv%S) .or. associated(tv%T)) then
393 call get_param(pf, mdl, "S_REF", s_ref, &
394 units="ppt", default=35.0, scale=us%ppt_to_S, do_not_log=.true.)
395 call get_param(pf, mdl, "DOME_T_LIGHT", t_light, &
396 "A first guess at the temperature of the lightest layer in the DOME test case.", &
397 units="degC", default=25.0, scale=us%degC_to_C)
398 endif
399
400 if (.not.associated(obc)) return
401
402 if (gv%Boussinesq) then
403 g_prime_tot = (gv%g_Earth / gv%Rho0) * rlay_range
404 def_rad = sqrt(d_edge*g_prime_tot) / abs(f_inflow)
405 tr_0 = (-d_edge*sqrt(d_edge*g_prime_tot)*0.5*def_rad) * gv%Z_to_H
406 else
407 g_prime_tot = (gv%g_Earth / (rlay_ref + 0.5*rlay_range)) * rlay_range
408 def_rad = sqrt(d_edge*g_prime_tot) / abs(f_inflow)
409 tr_0 = (-d_edge*sqrt(d_edge*g_prime_tot)*0.5*def_rad) * (rlay_ref + 0.5*rlay_range) * gv%RZ_to_H
410 endif
411
412 i_def_rad = 1.0 / ((1.0e-3*us%L_to_m*km_to_grid_unit) * def_rad)
413 ! This is mathematically equivalent to
414 ! I_Def_Rad = G%grid_unit_to_L / Def_Rad
415
416 if (obc%number_of_segments /= 1) then
417 call mom_error(warning, 'Error in DOME OBC segment setup', .true.)
418 return !!! Need a better error message here
419 endif
420
421 segment => obc%segment(1)
422 if (.not. segment%on_pe) return
423
424 ! Set up space for the OBCs to use for all the tracers.
425 ntherm = 0
426 if (associated(tv%S)) ntherm = ntherm + 1
427 if (associated(tv%T)) ntherm = ntherm + 1
428 allocate(segment%field(ntherm+tr_reg%ntr))
429
430 do k=1,nz
431 rst = -1.0
432 if (k>1) rst = -1.0 + (real(k-1)-0.5)/real(nz-1)
433
434 rsb = 0.0
435 if (k<nz) rsb = -1.0 + (real(k-1)+0.5)/real(nz-1)
436 rc = -1.0 + real(k-1)/real(nz-1)
437
438 ! These come from assuming geostrophy and a constant Ri profile.
439 yt = (2.0*ri_trans*rst + ri_trans + 2.0)/(2.0 - ri_trans)
440 yb = (2.0*ri_trans*rsb + ri_trans + 2.0)/(2.0 - ri_trans)
441 tr_k = tr_0 * (2.0/(ri_trans*(2.0-ri_trans))) * &
442 ((log(yt)+1.0)/yt - (log(yb)+1.0)/yb)
443 v_k = -sqrt(d_edge*g_prime_tot)*log((2.0 + ri_trans*(1.0 + 2.0*rc)) / &
444 (2.0 - ri_trans))
445 if (k == nz) tr_k = tr_k + tr_0 * (2.0/(ri_trans*(2.0+ri_trans))) * &
446 log((2.0+ri_trans)/(2.0-ri_trans))
447 ! New way
448 isd = segment%HI%isd ; ied = segment%HI%ied
449 jsdb = segment%HI%JsdB ; jedb = segment%HI%JedB
450 do j=jsdb,jedb ; do i=isd,ied
451 ! Here lon_im1 estimates G%geoLonBu(I-1,J), which may not have been set if
452 ! the symmetric memory mode is not being used.
453 lon_im1 = 2.0*g%geoLonCv(i,j) - g%geoLonBu(i,j)
454 segment%normal_trans(i,j,k) = tr_k * (exp(-2.0*(lon_im1 - inflow_lon) * i_def_rad) - &
455 exp(-2.0*(g%geoLonBu(i,j) - inflow_lon) * i_def_rad))
456 segment%normal_vel(i,j,k) = v_k * exp(-2.0*(g%geoLonCv(i,j) - inflow_lon) * i_def_rad)
457 enddo ; enddo
458 enddo
459
460 ! The inflow values of temperature and salinity also need to be set here if
461 ! these variables are used. The following code is just a naive example.
462 if (associated(tv%S)) then
463 ! In this example, all S inflows have values given by S_ref.
464 name = 'salt'
465 call tracer_name_lookup(tr_reg, ntr_id, tr_ptr, name)
466 call register_segment_tracer(tr_ptr, ntr_id, pf, gv, segment, obc_scalar=s_ref, scale=us%ppt_to_S)
467 endif
468 if (associated(tv%T)) then
469 ! In this example, the T values are set to be consistent with the layer
470 ! target density and a salinity of S_ref. This code is taken from
471 ! USER_initialize_temp_sal.
472 pres(:) = tv%P_Ref ; s0(:) = s_ref ; t0(1) = t_light
473 call calculate_density(t0(1), s0(1), pres(1), rho_guess(1), tv%eqn_of_state)
474 call calculate_density_derivs(t0, s0, pres, drho_dt, drho_ds, tv%eqn_of_state, (/1,1/) )
475
476 do k=1,nz ; t0(k) = t0(1) + (gv%Rlay(k)-rho_guess(1)) / drho_dt(1) ; enddo
477 do itt=1,6
478 call calculate_density(t0, s0, pres, rho_guess, tv%eqn_of_state)
479 call calculate_density_derivs(t0, s0, pres, drho_dt, drho_ds, tv%eqn_of_state)
480 do k=1,nz ; t0(k) = t0(k) + (gv%Rlay(k)-rho_guess(k)) / drho_dt(k) ; enddo
481 enddo
482
483 ! Temperature is tracer 1 for the OBCs.
484 allocate(segment%field(1)%buffer_src(segment%HI%isd:segment%HI%ied,segment%HI%JsdB:segment%HI%JedB,nz))
485 do k=1,nz ; do j=jsdb,jedb ; do i=isd,ied
486 ! With the revised OBC code, buffer_src uses the same rescaled units as for tracers.
487 segment%field(1)%buffer_src(i,j,k) = t0(k)
488 enddo ; enddo ; enddo
489 name = 'temp'
490 call tracer_name_lookup(tr_reg, ntr_id, tr_ptr, name)
491 call register_segment_tracer(tr_ptr, ntr_id, pf, gv, segment, obc_array=.true., scale=us%degC_to_C)
492 endif
493
494 ! Set up dye tracers
495 ! First dye - only one with OBC values
496 ! This field(ntherm+1) requires tr_D1 to be the first tracer after temperature and salinity.
497 allocate(segment%field(ntherm+1)%buffer_src(segment%HI%isd:segment%HI%ied,segment%HI%JsdB:segment%HI%JedB,nz))
498 do k=1,nz ; do j=segment%HI%jsd,segment%HI%jed ; do i=segment%HI%isd,segment%HI%ied
499 if (k < nz/2) then ; segment%field(ntherm+1)%buffer_src(i,j,k) = 0.0
500 else ; segment%field(ntherm+1)%buffer_src(i,j,k) = 1.0 ; endif
501 enddo ; enddo ; enddo
502 name = 'tr_D1'
503 call tracer_name_lookup(tr_reg, ntr_id, tr_ptr, name)
504 call register_segment_tracer(tr_ptr, ntr_id, pf, gv, obc%segment(1), obc_array=.true.)
505
506 ! All tracers but the first have 0 concentration in their inflows. As 0 is the
507 ! default value for the inflow concentrations, the following calls are unnecessary.
508 do m=2,tr_reg%ntr
509 write(name,'("tr_D",I0)') m
510 call tracer_name_lookup(tr_reg, ntr_id, tr_ptr, name)
511 call register_segment_tracer(tr_ptr, ntr_id, pf, gv, obc%segment(1), obc_scalar=0.0)
512 enddo
513
514end subroutine dome_set_obc_data
515
516end module dome_initialization