MOM_tracer_Z_init.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!> Used to initialize tracers from a depth- (or z*-) space file.
7
8use mom_error_handler, only : mom_error, fatal, warning, mom_mesg, is_root_pe
9use mom_file_parser, only : get_param, log_version, param_file_type
10use mom_grid, only : ocean_grid_type
11use mom_io, only : mom_read_data, get_var_sizes, read_attribute, read_variable
13use mom_eos, only : eos_type, calculate_density, calculate_density_derivs, eos_domain
16
17implicit none ; private
18
19#include <MOM_memory.h>
20
22
23! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
24! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
25! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
26! vary with the Boussinesq approximation, the Boussinesq variant is given first.
27
28contains
29
30!> This function initializes a tracer by reading a Z-space file, returning
31!! .true. if this appears to have been successful, and false otherwise.
32function tracer_z_init(tr, h, filename, tr_name, G, GV, US, missing_val, land_val, scale)
33 logical :: tracer_z_init !< A return code indicating if the initialization has been successful
34 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
35 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
36 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
37 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
38 intent(out) :: tr !< The tracer to initialize [CU ~> conc]
39 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
40 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] or other
41 !! arbitrary units such as [Z ~> m]
42 character(len=*), intent(in) :: filename !< The name of the file to read from
43 character(len=*), intent(in) :: tr_name !< The name of the tracer in the file
44 real, optional, intent(in) :: missing_val !< The missing value for the tracer [CU ~> conc]
45 real, optional, intent(in) :: land_val !< A value to use to fill in land points [CU ~> conc]
46 real, optional, intent(in) :: scale !< A factor by which to scale the output tracers from the
47 !! their units in the file [CU conc-1 ~> 1]
48
49 ! Local variables
50 real, allocatable, dimension(:,:,:) :: &
51 tr_in ! The z-space array of tracer concentrations that is read in [CU ~> conc]
52 real, allocatable, dimension(:) :: &
53 z_edges, & ! The depths of the cell edges or cell centers (depending on
54 ! the value of has_edges) in the input z* data [Z ~> m].
55 tr_1d, & ! A copy of the input tracer concentrations in a column [CU ~> conc]
56 wt, & ! The fractional weight for each layer in the range between
57 ! k_top and k_bot [nondim]
58 z1, z2 ! z1 and z2 are the depths of the top and bottom limits of the part
59 ! of a z-cell that contributes to a layer, relative to the cell
60 ! center and normalized by the cell thickness [nondim].
61 ! Note that -1/2 <= z1 <= z2 <= 1/2.
62 real :: e(szk_(gv)+1) ! The z-star interface heights [Z ~> m].
63 real :: landval ! The tracer value to use in land points [CU ~> conc]
64 real :: sl_tr ! The normalized slope of the tracer
65 ! within the cell, in tracer units [CU ~> conc]
66 real :: htot(szi_(g)) ! The vertical sum of h [H ~> m or kg m-2].
67 real :: dilate ! The amount by which the thicknesses are dilated to
68 ! create a z-star coordinate [Z H-1 ~> nondim or m3 kg-1]
69 ! or other units reflecting those of h
70 real :: missing ! The missing value for the tracer [CU ~> conc]
71 real :: scale_fac ! A factor by which to scale the output tracers from the units in the
72 ! input file [CU conc-1 ~> 1]
73 ! This include declares and sets the variable "version".
74# include "version_variable.h"
75 logical :: has_edges, use_missing, zero_surface
76 character(len=80) :: loc_msg
77 integer :: k_top, k_bot, k_bot_prev, k_start
78 integer :: i, j, k, kz, is, ie, js, je, nz, nz_in
79
80 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
81
82 scale_fac = 1.0 ; if (present(scale)) then ; scale_fac = scale ; endif
83
84 landval = 0.0 ; if (present(land_val)) landval = land_val
85
86 zero_surface = .false. ! Make this false for errors to be fatal.
87
88 use_missing = .false.
89 if (present(missing_val)) then
90 use_missing = .true. ; missing = missing_val
91 endif
92
93 ! Find out the number of input levels and read the depth of the edges,
94 ! also modifying their sign convention to be monotonically decreasing.
95 call read_z_edges(filename, tr_name, z_edges, nz_in, has_edges, use_missing, &
96 missing, scale=us%m_to_Z, missing_scale=scale_fac)
97 if (nz_in < 1) then
98 tracer_z_init = .false.
99 return
100 endif
101
102 allocate(tr_in(g%isd:g%ied,g%jsd:g%jed,nz_in), source=0.0)
103 allocate(tr_1d(nz_in), source=0.0)
104 call mom_read_data(filename, tr_name, tr_in(:,:,:), g%Domain, scale=scale_fac)
105
106 ! Fill missing values from above? Use a "close" test to avoid problems
107 ! from type-conversion rounoff.
108 if (present(missing_val)) then
109 do j=js,je ; do i=is,ie
110 if (g%mask2dT(i,j) == 0.0) then
111 tr_in(i,j,1) = landval
112 elseif (abs(tr_in(i,j,1) - missing_val) <= 1e-6*abs(missing_val)) then
113 write(loc_msg,'(f7.2," N ",f7.2," E")') g%geoLatT(i,j), g%geoLonT(i,j)
114 if (zero_surface) then
115 call mom_error(warning, "tracer_Z_init: Missing value of "// &
116 trim(tr_name)//" found in an ocean point at "//trim(loc_msg)// &
117 " in "//trim(filename) )
118 tr_in(i,j,1) = 0.0
119 else
120 call mom_error(fatal, "tracer_Z_init: Missing value of "// &
121 trim(tr_name)//" found in an ocean point at "//trim(loc_msg)// &
122 " in "//trim(filename) )
123 endif
124 endif
125 enddo ; enddo
126 do k=2,nz_in ; do j=js,je ; do i=is,ie
127 if (abs(tr_in(i,j,k) - missing_val) <= 1e-6*abs(missing_val)) &
128 tr_in(i,j,k) = tr_in(i,j,k-1)
129 enddo ; enddo ; enddo
130 endif
131
132 allocate(wt(nz_in+1)) ; allocate(z1(nz_in+1)) ; allocate(z2(nz_in+1))
133
134 ! This is a placeholder, and will be replaced with our full vertical
135 ! interpolation machinery when it is in place.
136 if (has_edges) then
137 do j=js,je
138 do i=is,ie ; htot(i) = 0.0 ; enddo
139 do k=1,nz ; do i=is,ie ; htot(i) = htot(i) + h(i,j,k) ; enddo ; enddo
140
141 do i=is,ie ; if (g%mask2dT(i,j)*htot(i) > 0.0) then
142 ! Determine the z* heights of the model interfaces.
143 dilate = (g%bathyT(i,j) + g%Z_ref) / htot(i)
144 e(nz+1) = -g%bathyT(i,j) - g%Z_ref
145 do k=nz,1,-1 ; e(k) = e(k+1) + dilate * h(i,j,k) ; enddo
146
147 ! Create a single-column copy of tr_in. Efficiency is not an issue here.
148 do k=1,nz_in ; tr_1d(k) = tr_in(i,j,k) ; enddo
149 k_bot = 1 ; k_bot_prev = -1
150 do k=1,nz
151 if (e(k+1) > z_edges(1)) then
152 tr(i,j,k) = tr_1d(1)
153 elseif (e(k) < z_edges(nz_in+1)) then
154 tr(i,j,k) = tr_1d(nz_in)
155 else
156 k_start = k_bot ! The starting point for this search
157 call find_overlap(z_edges, e(k), e(k+1), nz_in, &
158 k_start, k_top, k_bot, wt, z1, z2)
159 kz = k_top
160 if (kz /= k_bot_prev) then
161 ! Calculate the intra-cell profile.
162 sl_tr = 0.0 ! ; cur_tr = 0.0
163 if ((kz < nz_in) .and. (kz > 1)) &
164 sl_tr = find_limited_slope(tr_1d, z_edges, kz)
165 endif
166 ! This is the piecewise linear form.
167 tr(i,j,k) = wt(kz) * (tr_1d(kz) + 0.5*sl_tr*(z2(kz) + z1(kz)))
168 ! For the piecewise parabolic form add the following...
169 ! + C1_3*cur_tr*(z2(kz)**2 + z2(kz)*z1(kz) + z1(kz)**2))
170 do kz=k_top+1,k_bot-1
171 tr(i,j,k) = tr(i,j,k) + wt(kz)*tr_1d(kz)
172 enddo
173 if (k_bot > k_top) then
174 kz = k_bot
175 ! Calculate the intra-cell profile.
176 sl_tr = 0.0 ! ; cur_tr = 0.0
177 if ((kz < nz_in) .and. (kz > 1)) &
178 sl_tr = find_limited_slope(tr_1d, z_edges, kz)
179 ! This is the piecewise linear form.
180 tr(i,j,k) = tr(i,j,k) + wt(kz) * &
181 (tr_1d(kz) + 0.5*sl_tr*(z2(kz) + z1(kz)))
182 ! For the piecewise parabolic form add the following...
183 ! + C1_3*cur_tr*(z2(kz)**2 + z2(kz)*z1(kz) + z1(kz)**2))
184 endif
185 k_bot_prev = k_bot
186
187 ! Now handle the unlikely case where the layer partially extends
188 ! past the valid range of the input data by extrapolating using
189 ! the top or bottom value.
190 if ((e(k) > z_edges(1)) .and. (z_edges(nz_in+1) > e(k+1))) then
191 tr(i,j,k) = (((e(k) - z_edges(1)) * tr_1d(1) + &
192 (z_edges(1) - z_edges(nz_in)) * tr(i,j,k)) + &
193 (z_edges(nz_in+1) - e(k+1)) * tr_1d(nz_in)) / &
194 (e(k) - e(k+1))
195 elseif (e(k) > z_edges(1)) then
196 tr(i,j,k) = ((e(k) - z_edges(1)) * tr_1d(1) + &
197 (z_edges(1) - e(k+1)) * tr(i,j,k)) / &
198 (e(k) - e(k+1))
199 elseif (z_edges(nz_in) > e(k+1)) then
200 tr(i,j,k) = ((e(k) - z_edges(nz_in+1)) * tr(i,j,k) + &
201 (z_edges(nz_in+1) - e(k+1)) * tr_1d(nz_in)) / &
202 (e(k) - e(k+1))
203 endif
204 endif
205 enddo ! k-loop
206 else
207 do k=1,nz ; tr(i,j,k) = landval ; enddo
208 endif ; enddo ! i-loop
209 enddo ! j-loop
210 else
211 ! Without edge values, integrate a linear interpolation between cell centers.
212 do j=js,je
213 do i=is,ie ; htot(i) = 0.0 ; enddo
214 do k=1,nz ; do i=is,ie ; htot(i) = htot(i) + h(i,j,k) ; enddo ; enddo
215
216 do i=is,ie ; if (g%mask2dT(i,j)*htot(i) > 0.0) then
217 ! Determine the z* heights of the model interfaces.
218 dilate = (g%bathyT(i,j) + g%Z_ref) / htot(i)
219 e(nz+1) = -g%bathyT(i,j) - g%Z_ref
220 do k=nz,1,-1 ; e(k) = e(k+1) + dilate * h(i,j,k) ; enddo
221
222 ! Create a single-column copy of tr_in. Efficiency is not an issue here.
223 do k=1,nz_in ; tr_1d(k) = tr_in(i,j,k) ; enddo
224 k_bot = 1
225 do k=1,nz
226 if (e(k+1) > z_edges(1)) then
227 tr(i,j,k) = tr_1d(1)
228 elseif (z_edges(nz_in) > e(k)) then
229 tr(i,j,k) = tr_1d(nz_in)
230 else
231 k_start = k_bot ! The starting point for this search
232 call find_overlap(z_edges, e(k), e(k+1), nz_in-1, &
233 k_start, k_top, k_bot, wt, z1, z2)
234
235 kz = k_top
236 if (k_top < nz_in) then
237 tr(i,j,k) = wt(kz)*0.5*((tr_1d(kz) + tr_1d(kz+1)) + &
238 (tr_1d(kz+1) - tr_1d(kz))*(z2(kz)+z1(kz)))
239 else
240 tr(i,j,k) = wt(kz)*tr_1d(nz_in)
241 endif
242 do kz=k_top+1,k_bot-1
243 tr(i,j,k) = tr(i,j,k) + wt(kz)*0.5*(tr_1d(kz) + tr_1d(kz+1))
244 enddo
245 if (k_bot > k_top) then
246 kz = k_bot
247 tr(i,j,k) = tr(i,j,k) + wt(kz)*0.5*((tr_1d(kz) + tr_1d(kz+1)) + &
248 (tr_1d(kz+1) - tr_1d(kz))*(z2(kz)+z1(kz)))
249 endif
250
251 ! Now handle the case where the layer partially extends past
252 ! the valid range of the input data.
253 if ((e(k) > z_edges(1)) .and. (z_edges(nz_in) > e(k+1))) then
254 tr(i,j,k) = (((e(k) - z_edges(1)) * tr_1d(1) + &
255 (z_edges(1) - z_edges(nz_in)) * tr(i,j,k)) + &
256 (z_edges(nz_in) - e(k+1)) * tr_1d(nz_in)) / &
257 (e(k) - e(k+1))
258 elseif (e(k) > z_edges(1)) then
259 tr(i,j,k) = ((e(k) - z_edges(1)) * tr_1d(1) + &
260 (z_edges(1) - e(k+1)) * tr(i,j,k)) / &
261 (e(k) - e(k+1))
262 elseif (z_edges(nz_in) > e(k+1)) then
263 tr(i,j,k) = ((e(k) - z_edges(nz_in)) * tr(i,j,k) + &
264 (z_edges(nz_in) - e(k+1)) * tr_1d(nz_in)) / &
265 (e(k) - e(k+1))
266 endif
267 endif
268 enddo
269 else
270 do k=1,nz ; tr(i,j,k) = landval ; enddo
271 endif ; enddo ! i-loop
272 enddo ! j-loop
273 endif
274
275 deallocate(tr_in) ; deallocate(tr_1d) ; deallocate(z_edges)
276 deallocate(wt) ; deallocate(z1) ; deallocate(z2)
277
278 tracer_z_init = .true.
279
280end function tracer_z_init
281
282!> Layer model routine for remapping tracers from pseudo-z coordinates into layers defined
283!! by target interface positions.
284subroutine tracer_z_init_array(tr_in, z_edges, nk_data, e, land_fill, G, nlay, nlevs, &
285 eps_z, tr, scale)
286 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
287 integer, intent(in) :: nk_data !< The number of levels in the input data
288 real, dimension(SZI_(G),SZJ_(G),nk_data), &
289 intent(in) :: tr_in !< The z-space array of tracer concentrations
290 !! that is read in [A]
291 real, dimension(nk_data+1), intent(in) :: z_edges !< The depths of the cell edges in the input z* data
292 !! [Z ~> m] or [m]
293 integer, intent(in) :: nlay !< The number of vertical layers in the target grid
294 real, dimension(SZI_(G),SZJ_(G),nlay+1), &
295 intent(in) :: e !< The depths of the target layer interfaces [Z ~> m] or [m]
296 real, intent(in) :: land_fill !< fill in data over land [B]
297 integer, dimension(SZI_(G),SZJ_(G)), &
298 intent(in) :: nlevs !< The number of input levels with valid data
299 real, intent(in) :: eps_z !< A negligibly thin layer thickness [Z ~> m].
300 real, dimension(SZI_(G),SZJ_(G),nlay), &
301 intent(out) :: tr !< tracers in model space [B]
302 real, optional, intent(in) :: scale !< A factor by which to scale the output tracers from the
303 !! input tracers [B A-1 ~> 1]
304
305 ! Local variables
306 real :: tr_1d(nk_data) ! A copy of the input tracer concentrations in a column [B]
307 real :: e_1d(nlay+1) ! A 1-d column of interface heights, in the same units as e [Z ~> m] or [m]
308 real :: sl_tr ! The tracer concentration slope times the layer thickness, in tracer units [B]
309 real :: wt(nk_data) ! The fractional weight for each layer in the range between z1 and z2 [nondim]
310 real :: z1(nk_data) ! The fractional depth of the top limit of the part of a z-cell that contributes to
311 ! a layer, relative to the cell center and normalized by the cell thickness [nondim].
312 real :: z2(nk_data) ! The fractional depth of the bottom limit of the part of a z-cell that contributes to
313 ! a layer, relative to the cell center and normalized by the cell thickness [nondim].
314 ! Note that -1/2 <= z1 <= z2 <= 1/2.
315 real :: scale_fac ! A factor by which to scale the output tracers from the input tracers [B A-1 ~> 1]
316 integer :: k_top, k_bot, k_bot_prev, kstart
317 integer :: i, j, k, kz, is, ie, js, je
318
319 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
320
321 scale_fac = 1.0 ; if (present(scale)) then ; scale_fac = scale ; endif
322
323 do j=js,je
324 i_loop: do i=is,ie
325 if (nlevs(i,j) == 0 .or. g%mask2dT(i,j) == 0.) then
326 tr(i,j,:) = land_fill
327 cycle i_loop
328 endif
329
330 do k=1,nk_data
331 tr_1d(k) = scale_fac*tr_in(i,j,k)
332 enddo
333
334 do k=1,nlay+1
335 e_1d(k) = e(i,j,k)
336 enddo
337 k_bot = 1 ; k_bot_prev = -1
338 do k=1,nlay
339 if (e_1d(k+1) > z_edges(1)) then
340 tr(i,j,k) = tr_1d(1)
341 elseif (e_1d(k) < z_edges(nlevs(i,j)+1)) then
342 tr(i,j,k) = tr_1d(nlevs(i,j))
343
344 else
345 kstart = k_bot
346 call find_overlap(z_edges, e_1d(k), e_1d(k+1), nlevs(i,j), &
347 kstart, k_top, k_bot, wt, z1, z2)
348 kz = k_top
349 sl_tr = 0.0 ! ; cur_tr=0.0
350 if (kz /= k_bot_prev) then
351 ! Calculate the intra-cell profile.
352 if ((kz < nlevs(i,j)) .and. (kz > 1)) then
353 sl_tr = find_limited_slope(tr_1d, z_edges, kz)
354 endif
355 endif
356 if (kz > nlevs(i,j)) kz = nlevs(i,j)
357 ! This is the piecewise linear form.
358 tr(i,j,k) = wt(kz) * (tr_1d(kz) + 0.5*sl_tr*(z2(kz) + z1(kz)))
359 ! For the piecewise parabolic form add the following...
360 ! + C1_3*wt(kz) * cur_tr*(z2(kz)**2 + z2(kz)*z1(kz) + z1(kz)**2))
361 do kz=k_top+1,k_bot-1
362 tr(i,j,k) = tr(i,j,k) + wt(kz)*tr_1d(kz)
363 enddo
364
365 if (k_bot > k_top) then
366 kz = k_bot
367 ! Calculate the intra-cell profile.
368 sl_tr = 0.0 ! ; cur_tr = 0.0
369 if ((kz < nlevs(i,j)) .and. (kz > 1)) then
370 sl_tr = find_limited_slope(tr_1d, z_edges, kz)
371 endif
372 ! This is the piecewise linear form.
373 tr(i,j,k) = tr(i,j,k) + wt(kz) * (tr_1d(kz) + 0.5*sl_tr*(z2(kz) + z1(kz)))
374 ! For the piecewise parabolic form add the following...
375 ! + C1_3*cur_tr*(z2(kz)**2 + z2(kz)*z1(kz) + z1(kz)**2))
376 endif
377 k_bot_prev = k_bot
378
379 endif
380 enddo ! k-loop
381
382 do k=2,nlay ! simply fill vanished layers with adjacent value
383 if (e_1d(k)-e_1d(k+1) <= eps_z) tr(i,j,k) = tr(i,j,k-1)
384 enddo
385
386 enddo i_loop
387 enddo
388
389end subroutine tracer_z_init_array
390
391!> This subroutine reads the vertical coordinate data for a field from a NetCDF file.
392!! It also might read the missing value attribute for that same field.
393subroutine read_z_edges(filename, tr_name, z_edges, nz_out, has_edges, &
394 use_missing, missing, scale, missing_scale)
395 character(len=*), intent(in) :: filename !< The name of the file to read from.
396 character(len=*), intent(in) :: tr_name !< The name of the tracer in the file.
397 real, dimension(:), allocatable, &
398 intent(out) :: z_edges !< The depths of the vertical edges of the tracer array [Z ~> m]
399 integer, intent(out) :: nz_out !< The number of vertical layers in the tracer array
400 logical, intent(out) :: has_edges !< If true the values in z_edges are the edges of the
401 !! tracer cells, otherwise they are the cell centers
402 logical, intent(inout) :: use_missing !< If false on input, see whether the tracer has a
403 !! missing value, and if so return true
404 real, intent(inout) :: missing !< The missing value, if one has been found [CU ~> conc]
405 real, intent(in) :: scale !< A scaling factor for z_edges into new units [Z m-1 ~> 1]
406 real, intent(in) :: missing_scale !< A scaling factor to use to convert the
407 !! tracers and their missing value from the units in
408 !! the file into their internal units [CU conc-1 ~> 1]
409
410 ! This subroutine reads the vertical coordinate data for a field from a
411 ! NetCDF file. It also might read the missing value attribute for that same field.
412 character(len=32) :: mdl
413 character(len=120) :: tr_msg, dim_msg
414 character(:), allocatable :: edge_name
415 character(len=256) :: dim_names(4)
416 logical :: monotonic
417 integer :: ncid, k
418 integer :: nz_edge, ndim, sizes(4)
419
420 mdl = "MOM_tracer_Z_init read_Z_edges: "
421 tr_msg = trim(tr_name)//" in "//trim(filename)
422
423 if (is_root_pe()) then
424 call open_file_to_read(filename, ncid)
425 else
426 ncid = -1
427 endif
428
429 call get_var_sizes(filename, tr_name, ndim, sizes, dim_names=dim_names, ncid_in=ncid)
430 if ((ndim < 3) .or. (ndim > 4)) &
431 call mom_error(fatal, mdl//" "//trim(tr_msg)//" has too many or too few dimensions.")
432 nz_out = sizes(3)
433
434 if (.not.use_missing) then ! Try to find the missing value from the dataset.
435 call read_attribute(filename, "missing_value", missing, varname=tr_name, found=use_missing, ncid_in=ncid)
436 if (use_missing) missing = missing * missing_scale
437 endif
438 ! Find out if the Z-axis has an edges attribute
439 call read_attribute(filename, "edges", edge_name, varname=dim_names(3), found=has_edges, ncid_in=ncid)
440
441 nz_edge = sizes(3) ; if (has_edges) nz_edge = sizes(3)+1
442 allocate(z_edges(nz_edge), source=0.0)
443
444 if (nz_out < 1) return
445
446 ! Read the right variable.
447 if (has_edges) then
448 call read_variable(filename, edge_name, z_edges, ncid)
449 else
450 call read_variable(filename, dim_names(3), z_edges, ncid)
451 endif
452 call close_file_to_read(ncid, filename)
453 if (allocated(edge_name)) deallocate(edge_name)
454
455 ! z_edges should be montonically decreasing with our sign convention.
456 ! Change the sign sign convention if it looks like z_edges is increasing.
457 if (z_edges(1) < z_edges(2)) then
458 do k=1,nz_edge ; z_edges(k) = -z_edges(k) ; enddo
459 endif
460 ! Check that z_edges is now monotonically decreasing.
461 monotonic = .true.
462 do k=2,nz_edge ; if (z_edges(k) >= z_edges(k-1)) monotonic = .false. ; enddo
463 if (.not.monotonic) call mom_error(warning,mdl//" "//trim(dim_msg)//" is not monotonic.")
464
465 if (scale /= 1.0) then ; do k=1,nz_edge ; z_edges(k) = scale*z_edges(k) ; enddo ; endif
466
467end subroutine read_z_edges
468
469!> Determines the layers bounded by interfaces e that overlap
470!! with the depth range between Z_top and Z_bot, and the fractional weights
471!! of each layer. It also calculates the normalized relative depths of the range
472!! of each layer that overlaps that depth range.
473subroutine find_overlap(e, Z_top, Z_bot, k_max, k_start, k_top, k_bot, wt, z1, z2)
474 real, dimension(:), intent(in) :: e !< Column interface heights, [Z ~> m] or other units.
475 real, intent(in) :: Z_top !< Top of range being mapped to, in the units of e [Z ~> m].
476 real, intent(in) :: Z_bot !< Bottom of range being mapped to, in the units of e [Z ~> m].
477 integer, intent(in) :: k_max !< Number of valid layers.
478 integer, intent(in) :: k_start !< Layer at which to start searching.
479 integer, intent(out) :: k_top !< Indices of top layers that overlap with the depth range.
480 integer, intent(out) :: k_bot !< Indices of bottom layers that overlap with the depth range.
481 real, dimension(:), intent(out) :: wt !< Relative weights of each layer from k_top to k_bot [nondim].
482 real, dimension(:), intent(out) :: z1 !< Depth of the top limits of the part of
483 !! a layer that contributes to a depth level, relative to the cell center and normalized
484 !! by the cell thickness [nondim]. Note that -1/2 <= z1 < z2 <= 1/2.
485 real, dimension(:), intent(out) :: z2 !< Depths of the bottom limit of the part of
486 !! a layer that contributes to a depth level, relative to the cell center and normalized
487 !! by the cell thickness [nondim]. Note that -1/2 <= z1 < z2 <= 1/2.
488
489 ! Local variables
490 real :: Ih ! The inverse of the vertical distance across a layer, in the inverse of the units of e [Z-1 ~> m-1]
491 real :: e_c ! The height of the layer center, in the units of e [Z ~> m]
492 real :: tot_wt ! The sum of the thicknesses contributing to a layer [Z ~> m]
493 real :: I_totwt ! The Adcroft reciprocal of tot_wt [Z-1 ~> m-1]
494 integer :: k
495
496 wt(:) = 0.0 ; z1(:) = 0.0 ; z2(:) = 0.0 ; k_bot = k_max
497
498 do k=k_start,k_max ; if (e(k+1) < z_top) exit ; enddo
499 k_top = k
500 if (k_top > k_max) return
501
502 ! Determine the fractional weights of each layer.
503 ! Note that by convention, e and Z_int decrease with increasing k.
504 if (e(k+1) <= z_bot) then
505 wt(k) = 1.0 ; k_bot = k
506 ih = 0.0 ; if (e(k) /= e(k+1)) ih = 1.0 / (e(k)-e(k+1))
507 e_c = 0.5*(e(k)+e(k+1))
508 z1(k) = (e_c - min(e(k), z_top)) * ih
509 z2(k) = (e_c - z_bot) * ih
510 else
511 ! Note that in theis branch, wt temporarily has units of [Z ~> m]
512 wt(k) = min(e(k),z_top) - e(k+1) ; tot_wt = wt(k) ! These are always > 0.
513 if (e(k) /= e(k+1)) then
514 z1(k) = (0.5*(e(k)+e(k+1)) - min(e(k), z_top)) / (e(k)-e(k+1))
515 else ; z1(k) = -0.5 ; endif
516 z2(k) = 0.5
517 k_bot = k_max
518 do k=k_top+1,k_max
519 if (e(k+1) <= z_bot) then
520 k_bot = k
521 wt(k) = e(k) - z_bot ; z1(k) = -0.5
522 if (e(k) /= e(k+1)) then
523 z2(k) = (0.5*(e(k)+e(k+1)) - z_bot) / (e(k)-e(k+1))
524 else ; z2(k) = 0.5 ; endif
525 else
526 wt(k) = e(k) - e(k+1) ; z1(k) = -0.5 ; z2(k) = 0.5
527 endif
528 tot_wt = tot_wt + wt(k) ! wt(k) is always > 0.
529 if (k>=k_bot) exit
530 enddo
531
532 i_totwt = 0.0 ; if (tot_wt > 0.0) i_totwt = 1.0 / tot_wt
533 ! This loop changes the units of wt from [Z ~> m] to [nondim].
534 do k=k_top,k_bot ; wt(k) = i_totwt*wt(k) ; enddo
535 endif
536
537end subroutine find_overlap
538
539!> This subroutine determines a limited slope for val to be advected with
540!! a piecewise limited scheme.
541function find_limited_slope(val, e, k) result(slope)
542 real, dimension(:), intent(in) :: val !< A column of the values that are being interpolated, in arbitrary units [A]
543 real, dimension(:), intent(in) :: e !< A column's interface heights [Z ~> m] or other units.
544 integer, intent(in) :: k !< The layer whose slope is being determined.
545 real :: slope !< The normalized slope in the intracell distribution of val [A]
546 ! Local variables
547 real :: amn, cmn ! Limited differences and curvatures in the values [A]
548 real :: d1, d2 ! Layer thicknesses, in the units of e [Z ~> m]
549
550 if ((val(k)-val(k-1)) * (val(k)-val(k+1)) >= 0.0) then
551 slope = 0.0 ! ; curvature = 0.0
552 else
553 d1 = 0.5*(e(k-1)-e(k+1)) ; d2 = 0.5*(e(k)-e(k+2))
554 if (d1*d2 > 0.0) then
555 slope = ((d1**2)*(val(k+1) - val(k)) + (d2**2)*(val(k) - val(k-1))) * &
556 (e(k) - e(k+1)) / (d1*d2*(d1+d2))
557 ! slope = 0.5*(val(k+1) - val(k-1))
558 ! This is S.J. Lin's form of the PLM limiter.
559 amn = min(abs(slope), 2.0*(max(val(k-1), val(k), val(k+1)) - val(k)))
560 cmn = 2.0*(val(k) - min(val(k-1), val(k), val(k+1)))
561 slope = sign(1.0, slope) * min(amn, cmn)
562
563 ! min(abs(slope), 2.0*(max(val(k-1),val(k),val(k+1)) - val(k)), &
564 ! 2.0*(val(k) - min(val(k-1),val(k),val(k+1))))
565 ! curvature = 0.0
566 else
567 slope = 0.0 ! ; curvature = 0.0
568 endif
569 endif
570
571end function find_limited_slope
572
573!> This subroutine determines the potential temperature and salinity that
574!! is consistent with the target density using provided initial guess
575subroutine determine_temperature(temp, salt, R_tgt, EOS, p_ref, niter, k_start, G, GV, US, PF, &
576 just_read)
577 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
578 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
579 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
580 intent(inout) :: temp !< potential temperature [C ~> degC]
581 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
582 intent(inout) :: salt !< salinity [S ~> ppt]
583 real, dimension(SZK_(GV)), intent(in) :: r_tgt !< desired potential density [R ~> kg m-3].
584 type(eos_type), intent(in) :: eos !< seawater equation of state control structure
585 real, intent(in) :: p_ref !< reference pressure [R L2 T-2 ~> Pa].
586 integer, intent(in) :: niter !< maximum number of iterations
587 integer, intent(in) :: k_start !< starting index (i.e. below the buffer layer)
588 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
589 type(param_file_type), intent(in) :: pf !< A structure indicating the open file
590 !! to parse for model parameter values.
591 logical, intent(in) :: just_read !< If true, this call will only read
592 !! parameters without changing T or S.
593
594 ! Local variables (All of which need documentation!)
595 real, dimension(SZI_(G),SZK_(GV)) :: &
596 t, & ! A 2-d working copy of the layer temperatures [C ~> degC]
597 s, & ! A 2-d working copy of the layer salinities [S ~> ppt]
598 dt, & ! An estimated change in temperature before bounding [C ~> degC]
599 ds, & ! An estimated change in salinity before bounding [S ~> ppt]
600 rho, & ! Layer densities with the current estimate of temperature and salinity [R ~> kg m-3]
601 drho_dt, & ! Partial derivative of density with temperature [R C-1 ~> kg m-3 degC-1]
602 drho_ds ! Partial derivative of density with salinity [R S-1 ~> kg m-3 ppt-1]
603 real, dimension(SZI_(G)) :: press ! Reference pressures [R L2 T-2 ~> Pa]
604 real :: dt_ds_gauge ! The relative penalizing of temperature to salinity changes when
605 ! minimizing property changes while correcting density [C S-1 ~> degC ppt-1].
606 real :: i_denom ! The inverse of the magnitude squared of the density gradient in
607 ! T-S space when stretched with dT_dS_gauge [S2 R-2 ~> ppt2 m6 kg-2]
608 real :: t_min, t_max ! The minimum and maximum temperatures [C ~> degC]
609 real :: s_min, s_max ! Minimum and maximum salinities [S ~> ppt]
610 real :: tol_t ! The tolerance for temperature matches [C ~> degC]
611 real :: tol_s ! The tolerance for salinity matches [S ~> ppt]
612 real :: tol_rho ! The tolerance for density matches [R ~> kg m-3]
613 real :: max_t_adj ! The largest permitted temperature changes with each iteration
614 ! when old_fit is true [C ~> degC]
615 real :: max_s_adj ! The largest permitted salinity changes with each iteration
616 ! when old_fit is true [S ~> ppt]
617 ! This include declares and sets the variable "version".
618# include "version_variable.h"
619 character(len=40) :: mdl = "determine_temperature" ! This subroutine's name.
620 logical :: domore(szk_(gv)) ! Records which layers need additional iterations
621 logical :: adjust_salt, fit_together, convergence_bug, do_any
622 integer, dimension(2) :: eosdom ! The i-computational domain for the equation of state
623 integer :: i, j, k, is, ie, js, je, nz, itt
624
625 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
626
627 call log_version(pf, mdl, version, "")
628
629 ! We should switch the default to the newer method which simultaneously adjusts
630 ! temp and salt based on the ratio of the thermal and haline coefficients, once it is tested.
631 call get_param(pf, mdl, "DETERMINE_TEMP_ADJUST_T_AND_S", fit_together, &
632 "If true, simltaneously adjust the estimates of the temperature and salinity "//&
633 "based on the ratio of the thermal and haline coefficients. Otherwise try to "//&
634 "match the density by only adjusting temperatures within a maximum range before "//&
635 "revising estimates of the salinity.", default=.false., do_not_log=just_read)
636 call get_param(pf, mdl, "DETERMINE_TEMP_CONVERGENCE_BUG", convergence_bug, &
637 "If true, use layout-dependent tests on the changes in temperature and salinity "//&
638 "to determine when the iterations have converged when DETERMINE_TEMP_ADJUST_T_AND_S "//&
639 "is false. For realistic equations of state and the default values of the "//&
640 "various tolerances, this bug does not impact the solutions.", &
641 default=.false., do_not_log=just_read)
642
643 call get_param(pf, mdl, "DETERMINE_TEMP_T_MIN", t_min, &
644 "The minimum temperature that can be found by determine_temperature.", &
645 units="degC", default=-2.0, scale=us%degC_to_C, do_not_log=just_read)
646 call get_param(pf, mdl, "DETERMINE_TEMP_T_MAX", t_max, &
647 "The maximum temperature that can be found by determine_temperature.", &
648 units="degC", default=31.0, scale=us%degC_to_C, do_not_log=just_read)
649 call get_param(pf, mdl, "DETERMINE_TEMP_S_MIN", s_min, &
650 "The minimum salinity that can be found by determine_temperature.", &
651 units="ppt", default=0.5, scale=us%ppt_to_S, do_not_log=just_read)
652 call get_param(pf, mdl, "DETERMINE_TEMP_S_MAX", s_max, &
653 "The maximum salinity that can be found by determine_temperature.", &
654 units="ppt", default=65.0, scale=us%ppt_to_S, do_not_log=just_read)
655 call get_param(pf, mdl, "DETERMINE_TEMP_T_TOLERANCE", tol_t, &
656 "The convergence tolerance for temperature in determine_temperature.", &
657 units="degC", default=1.0e-4, scale=us%degC_to_C, &
658 do_not_log=just_read.or.(.not.convergence_bug))
659 call get_param(pf, mdl, "DETERMINE_TEMP_S_TOLERANCE", tol_s, &
660 "The convergence tolerance for temperature in determine_temperature.", &
661 units="ppt", default=1.0e-4, scale=us%ppt_to_S, &
662 do_not_log=just_read.or.(.not.convergence_bug))
663 call get_param(pf, mdl, "DETERMINE_TEMP_RHO_TOLERANCE", tol_rho, &
664 "The convergence tolerance for density in determine_temperature.", &
665 units="kg m-3", default=1.0e-4, scale=us%kg_m3_to_R, do_not_log=just_read)
666 if (fit_together) then
667 ! By default 10 degC is weighted equivalently to 1 ppt when minimizing changes.
668 call get_param(pf, mdl, "DETERMINE_TEMP_DT_DS_WEIGHT", dt_ds_gauge, &
669 "When extrapolating T & S to match the layer target densities, this "//&
670 "factor (in degC / ppt) is combined with the derivatives of density "//&
671 "with T & S to determine what direction is orthogonal to density contours. "//&
672 "It could be based on a typical value of (dR/dS) / (dR/dT) in oceanic profiles.", &
673 units="degC ppt-1", default=10.0, scale=us%degC_to_C*us%S_to_ppt)
674 else
675 call get_param(pf, mdl, "DETERMINE_TEMP_T_ADJ_RANGE", max_t_adj, &
676 "The maximum amount by which the initial layer temperatures can be "//&
677 "modified in determine_temperature.", &
678 units="degC", default=1.0, scale=us%degC_to_C, do_not_log=just_read)
679 call get_param(pf, mdl, "DETERMINE_TEMP_S_ADJ_RANGE", max_s_adj, &
680 "The maximum amount by which the initial layer salinities can be "//&
681 "modified in determine_temperature.", &
682 units="ppt", default=0.5, scale=us%ppt_to_S, do_not_log=just_read)
683 endif
684
685 if (just_read) return ! All run-time parameters have been read, so return.
686
687 press(:) = p_ref
688 eosdom(:) = eos_domain(g%HI)
689
690 do j=js,je
691 ds(:,:) = 0. ! Needs to be zero everywhere since there is a maxval(abs(dS)) later...
692 t(:,:) = temp(:,j,:)
693 s(:,:) = salt(:,j,:)
694 dt(:,:) = 0.0
695 domore(:) = .true.
696 adjust_salt = .true.
697 iter_loop: do itt = 1,niter
698 do k=k_start,nz ; if (domore(k)) then
699 domore(k) = .false.
700 call calculate_density(t(:,k), s(:,k), press, rho(:,k), eos, eosdom )
701 call calculate_density_derivs(t(:,k), s(:,k), press, drho_dt(:,k), drho_ds(:,k), &
702 eos, eosdom )
703 do i=is,ie
704! if (abs(rho(i,k)-R_tgt(k))>tol_rho .and. abs(T(i,k)-land_fill) < epsln) then
705 if (abs(rho(i,k)-r_tgt(k))>tol_rho) then
706 domore(k) = .true.
707 if (.not.fit_together) then
708 dt(i,k) = max(min((r_tgt(k)-rho(i,k)) / drho_dt(i,k), max_t_adj), -max_t_adj)
709 t(i,k) = max(min(t(i,k)+dt(i,k), t_max), t_min)
710 else
711 i_denom = 1.0 / (drho_ds(i,k)**2 + dt_ds_gauge**2*drho_dt(i,k)**2)
712 ds(i,k) = (r_tgt(k)-rho(i,k)) * drho_ds(i,k) * i_denom
713 dt(i,k) = (r_tgt(k)-rho(i,k)) * dt_ds_gauge**2*drho_dt(i,k) * i_denom
714
715 t(i,k) = max(min(t(i,k)+dt(i,k), t_max), t_min)
716 s(i,k) = max(min(s(i,k)+ds(i,k), s_max), s_min)
717 endif
718 endif
719 enddo
720 endif ; enddo
721 if (convergence_bug) then
722 ! If this test does anything, it is layout-dependent.
723 if (maxval(abs(dt)) < tol_t) then
724 adjust_salt = .false.
725 exit iter_loop
726 endif
727 endif
728
729 do_any = .false.
730 do k=k_start,nz ; if (domore(k)) do_any = .true. ; enddo
731 if (.not.do_any) exit iter_loop ! Further iterations will not change anything.
732 enddo iter_loop
733
734 if (adjust_salt .and. .not.fit_together) then ; do itt = 1,niter
735 do k=k_start,nz ; if (domore(k)) then
736 domore(k) = .false.
737 call calculate_density(t(:,k), s(:,k), press, rho(:,k), eos, eosdom )
738 call calculate_density_derivs(t(:,k), s(:,k), press, drho_dt(:,k), drho_ds(:,k), &
739 eos, eosdom )
740 do i=is,ie
741! if (abs(rho(i,k)-R_tgt(k))>tol_rho .and. abs(T(i,k)-land_fill) < epsln ) then
742 if (abs(rho(i,k)-r_tgt(k)) > tol_rho) then
743 ds(i,k) = max(min((r_tgt(k)-rho(i,k)) / drho_ds(i,k), max_s_adj), -max_s_adj)
744 s(i,k) = max(min(s(i,k)+ds(i,k), s_max), s_min)
745 domore(k) = .true.
746 endif
747 enddo
748 endif ; enddo
749
750 if (convergence_bug) then
751 ! If this test does anything, it is layout-dependent.
752 if (maxval(abs(ds)) < tol_s) exit
753 endif
754
755 do_any = .false.
756 do k=k_start,nz ; if (domore(k)) do_any = .true. ; enddo
757 if (.not.do_any) exit ! Further iterations will not change anything
758 enddo ; endif
759
760 temp(:,j,:) = t(:,:)
761 salt(:,j,:) = s(:,:)
762 enddo
763
764end subroutine determine_temperature
765
766end module mom_tracer_z_init