MOM6
MOM_tracer_Z_init.F90
Go to the documentation of this file.
1 !> Used to initialize tracers from a depth- (or z*-) space file.
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
6 !use MOM_diag_to_Z, only : find_overlap, find_limited_slope
7 use mom_error_handler, only : mom_error, fatal, warning, mom_mesg, is_root_pe
8 ! use MOM_file_parser, only : get_param, log_version, param_file_type
9 use mom_grid, only : ocean_grid_type
10 use mom_io, only : mom_read_data
12 
13 use netcdf
14 
15 implicit none ; private
16 
17 #include <MOM_memory.h>
18 
19 public tracer_z_init
20 
21 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
22 ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
23 ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
24 ! vary with the Boussinesq approximation, the Boussinesq variant is given first.
25 
26 contains
27 
28 !> This function initializes a tracer by reading a Z-space file, returning
29 !! .true. if this appears to have been successful, and false otherwise.
30 function tracer_z_init(tr, h, filename, tr_name, G, US, missing_val, land_val)
31  logical :: tracer_z_init !< A return code indicating if the initialization has been successful
32  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
33  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
34  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
35  intent(out) :: tr !< The tracer to initialize
36  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
37  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
38  character(len=*), intent(in) :: filename !< The name of the file to read from
39  character(len=*), intent(in) :: tr_name !< The name of the tracer in the file
40 ! type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
41  real, optional, intent(in) :: missing_val !< The missing value for the tracer
42  real, optional, intent(in) :: land_val !< A value to use to fill in land points
43 
44  ! This function initializes a tracer by reading a Z-space file, returning true if this
45  ! appears to have been successful, and false otherwise.
46 !
47  integer, save :: init_calls = 0
48 ! This include declares and sets the variable "version".
49 #include "version_variable.h"
50  character(len=40) :: mdl = "MOM_tracer_Z_init" ! This module's name.
51  character(len=256) :: mesg ! Message for error messages.
52 
53  real, allocatable, dimension(:,:,:) :: &
54  tr_in ! The z-space array of tracer concentrations that is read in.
55  real, allocatable, dimension(:) :: &
56  z_edges, & ! The depths of the cell edges or cell centers (depending on
57  ! the value of has_edges) in the input z* data [Z ~> m].
58  tr_1d, & ! A copy of the input tracer concentrations in a column.
59  wt, & ! The fractional weight for each layer in the range between
60  ! k_top and k_bot, nondim.
61  z1, & ! z1 and z2 are the depths of the top and bottom limits of the part
62  z2 ! of a z-cell that contributes to a layer, relative to the cell
63  ! center and normalized by the cell thickness, nondim.
64  ! Note that -1/2 <= z1 <= z2 <= 1/2.
65  real :: e(szk_(g)+1) ! The z-star interface heights [Z ~> m].
66  real :: landval ! The tracer value to use in land points.
67  real :: sl_tr ! The normalized slope of the tracer
68  ! within the cell, in tracer units.
69  real :: htot(szi_(g)) ! The vertical sum of h [H ~> m or kg m-2].
70  real :: dilate ! The amount by which the thicknesses are dilated to
71  ! create a z-star coordinate, nondim or in m3 kg-1.
72  real :: missing ! The missing value for the tracer.
73 
74  logical :: has_edges, use_missing, zero_surface
75  character(len=80) :: loc_msg
76  integer :: k_top, k_bot, k_bot_prev
77  integer :: i, j, k, kz, is, ie, js, je, nz, nz_in
78  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
79 
80  landval = 0.0 ; if (present(land_val)) landval = land_val
81 
82  zero_surface = .false. ! Make this false for errors to be fatal.
83 
84  use_missing = .false.
85  if (present(missing_val)) then
86  use_missing = .true. ; missing = missing_val
87  endif
88 
89  ! Find out the number of input levels and read the depth of the edges,
90  ! also modifying their sign convention to be monotonically decreasing.
91  call read_z_edges(filename, tr_name, z_edges, nz_in, has_edges, use_missing, &
92  missing, scale=us%m_to_Z)
93  if (nz_in < 1) then
94  tracer_z_init = .false.
95  return
96  endif
97 
98  allocate(tr_in(g%isd:g%ied,g%jsd:g%jed,nz_in)) ; tr_in(:,:,:) = 0.0
99  allocate(tr_1d(nz_in)) ; tr_1d(:) = 0.0
100  call mom_read_data(filename, tr_name, tr_in(:,:,:), g%Domain)
101 
102  ! Fill missing values from above? Use a "close" test to avoid problems
103  ! from type-conversion rounoff.
104  if (present(missing_val)) then
105  do j=js,je ; do i=is,ie
106  if (g%mask2dT(i,j) == 0.0) then
107  tr_in(i,j,1) = landval
108  elseif (abs(tr_in(i,j,1) - missing_val) <= 1e-6*abs(missing_val)) then
109  write(loc_msg,'(f7.2," N ",f7.2," E")') g%geoLatT(i,j), g%geoLonT(i,j)
110  if (zero_surface) then
111  call mom_error(warning, "tracer_Z_init: Missing value of "// &
112  trim(tr_name)//" found in an ocean point at "//trim(loc_msg)// &
113  " in "//trim(filename) )
114  tr_in(i,j,1) = 0.0
115  else
116  call mom_error(fatal, "tracer_Z_init: Missing value of "// &
117  trim(tr_name)//" found in an ocean point at "//trim(loc_msg)// &
118  " in "//trim(filename) )
119  endif
120  endif
121  enddo ; enddo
122  do k=2,nz_in ; do j=js,je ; do i=is,ie
123  if (abs(tr_in(i,j,k) - missing_val) <= 1e-6*abs(missing_val)) &
124  tr_in(i,j,k) = tr_in(i,j,k-1)
125  enddo ; enddo ; enddo
126  endif
127 
128  allocate(wt(nz_in+1)) ; allocate(z1(nz_in+1)) ; allocate(z2(nz_in+1))
129 
130  ! This is a placeholder, and will be replaced with our full vertical
131  ! interpolation machinery when it is in place.
132  if (has_edges) then
133  do j=js,je
134  do i=is,ie ; htot(i) = 0.0 ; enddo
135  do k=1,nz ; do i=is,ie ; htot(i) = htot(i) + h(i,j,k) ; enddo ; enddo
136 
137  do i=is,ie ; if (g%mask2dT(i,j)*htot(i) > 0.0) then
138  ! Determine the z* heights of the model interfaces.
139  dilate = (g%bathyT(i,j) - 0.0) / htot(i)
140  e(nz+1) = -g%bathyT(i,j)
141  do k=nz,1,-1 ; e(k) = e(k+1) + dilate * h(i,j,k) ; enddo
142 
143  ! Create a single-column copy of tr_in. ### CHANGE THIS LATER?
144  do k=1,nz_in ; tr_1d(k) = tr_in(i,j,k) ; enddo
145  k_bot = 1 ; k_bot_prev = -1
146  do k=1,nz
147  if (e(k+1) > z_edges(1)) then
148  tr(i,j,k) = tr_1d(1)
149  elseif (e(k) < z_edges(nz_in+1)) then
150  tr(i,j,k) = tr_1d(nz_in)
151  else
152  call find_overlap(z_edges, e(k), e(k+1), nz_in, &
153  k_bot, k_top, k_bot, wt, z1, z2)
154  kz = k_top
155  if (kz /= k_bot_prev) then
156  ! Calculate the intra-cell profile.
157  sl_tr = 0.0 ! ; cur_tr = 0.0
158  if ((kz < nz_in) .and. (kz > 1)) call &
159  find_limited_slope(tr_1d, z_edges, sl_tr, kz)
160  endif
161  ! This is the piecewise linear form.
162  tr(i,j,k) = wt(kz) * &
163  (tr_1d(kz) + 0.5*sl_tr*(z2(kz) + z1(kz)))
164  ! For the piecewise parabolic form add the following...
165  ! + C1_3*cur_tr*(z2(kz)**2 + z2(kz)*z1(kz) + z1(kz)**2))
166  do kz=k_top+1,k_bot-1
167  tr(i,j,k) = tr(i,j,k) + wt(kz)*tr_1d(kz)
168  enddo
169  if (k_bot > k_top) then
170  kz = k_bot
171  ! Calculate the intra-cell profile.
172  sl_tr = 0.0 ! ; cur_tr = 0.0
173  if ((kz < nz_in) .and. (kz > 1)) call &
174  find_limited_slope(tr_1d, z_edges, sl_tr, kz)
175  ! This is the piecewise linear form.
176  tr(i,j,k) = tr(i,j,k) + wt(kz) * &
177  (tr_1d(kz) + 0.5*sl_tr*(z2(kz) + z1(kz)))
178  ! For the piecewise parabolic form add the following...
179  ! + C1_3*cur_tr*(z2(kz)**2 + z2(kz)*z1(kz) + z1(kz)**2))
180  endif
181  k_bot_prev = k_bot
182 
183  ! Now handle the unlikely case where the layer partially extends
184  ! past the valid range of the input data by extrapolating using
185  ! the top or bottom value.
186  if ((e(k) > z_edges(1)) .and. (z_edges(nz_in+1) > e(k+1))) then
187  tr(i,j,k) = (((e(k) - z_edges(1)) * tr_1d(1) + &
188  (z_edges(1) - z_edges(nz_in)) * tr(i,j,k)) + &
189  (z_edges(nz_in+1) - e(k+1)) * tr_1d(nz_in)) / &
190  (e(k) - e(k+1))
191  elseif (e(k) > z_edges(1)) then
192  tr(i,j,k) = ((e(k) - z_edges(1)) * tr_1d(1) + &
193  (z_edges(1) - e(k+1)) * tr(i,j,k)) / &
194  (e(k) - e(k+1))
195  elseif (z_edges(nz_in) > e(k+1)) then
196  tr(i,j,k) = ((e(k) - z_edges(nz_in+1)) * tr(i,j,k) + &
197  (z_edges(nz_in+1) - e(k+1)) * tr_1d(nz_in)) / &
198  (e(k) - e(k+1))
199  endif
200  endif
201  enddo ! k-loop
202  else
203  do k=1,nz ; tr(i,j,k) = landval ; enddo
204  endif ; enddo ! i-loop
205  enddo ! j-loop
206  else
207  ! Without edge values, integrate a linear interpolation between cell centers.
208  do j=js,je
209  do i=is,ie ; htot(i) = 0.0 ; enddo
210  do k=1,nz ; do i=is,ie ; htot(i) = htot(i) + h(i,j,k) ; enddo ; enddo
211 
212  do i=is,ie ; if (g%mask2dT(i,j)*htot(i) > 0.0) then
213  ! Determine the z* heights of the model interfaces.
214  dilate = (g%bathyT(i,j) - 0.0) / htot(i)
215  e(nz+1) = -g%bathyT(i,j)
216  do k=nz,1,-1 ; e(k) = e(k+1) + dilate * h(i,j,k) ; enddo
217 
218  ! Create a single-column copy of tr_in. ### CHANGE THIS LATER?
219  do k=1,nz_in ; tr_1d(k) = tr_in(i,j,k) ; enddo
220  k_bot = 1
221  do k=1,nz
222  if (e(k+1) > z_edges(1)) then
223  tr(i,j,k) = tr_1d(1)
224  elseif (z_edges(nz_in) > e(k)) then
225  tr(i,j,k) = tr_1d(nz_in)
226  else
227  call find_overlap(z_edges, e(k), e(k+1), nz_in-1, &
228  k_bot, k_top, k_bot, wt, z1, z2)
229 
230  kz = k_top
231  if (k_top < nz_in) then
232  tr(i,j,k) = wt(kz)*0.5*((tr_1d(kz) + tr_1d(kz+1)) + &
233  (tr_1d(kz+1) - tr_1d(kz))*(z2(kz)+z1(kz)))
234  else
235  tr(i,j,k) = wt(kz)*tr_1d(nz_in)
236  endif
237  do kz=k_top+1,k_bot-1
238  tr(i,j,k) = tr(i,j,k) + wt(kz)*0.5*(tr_1d(kz) + tr_1d(kz+1))
239  enddo
240  if (k_bot > k_top) then
241  kz = k_bot
242  tr(i,j,k) = tr(i,j,k) + wt(kz)*0.5*((tr_1d(kz) + tr_1d(kz+1)) + &
243  (tr_1d(kz+1) - tr_1d(kz))*(z2(kz)+z1(kz)))
244  endif
245 
246  ! Now handle the case where the layer partially extends past
247  ! the valid range of the input data.
248  if ((e(k) > z_edges(1)) .and. (z_edges(nz_in) > e(k+1))) then
249  tr(i,j,k) = (((e(k) - z_edges(1)) * tr_1d(1) + &
250  (z_edges(1) - z_edges(nz_in)) * tr(i,j,k)) + &
251  (z_edges(nz_in) - e(k+1)) * tr_1d(nz_in)) / &
252  (e(k) - e(k+1))
253  elseif (e(k) > z_edges(1)) then
254  tr(i,j,k) = ((e(k) - z_edges(1)) * tr_1d(1) + &
255  (z_edges(1) - e(k+1)) * tr(i,j,k)) / &
256  (e(k) - e(k+1))
257  elseif (z_edges(nz_in) > e(k+1)) then
258  tr(i,j,k) = ((e(k) - z_edges(nz_in)) * tr(i,j,k) + &
259  (z_edges(nz_in) - e(k+1)) * tr_1d(nz_in)) / &
260  (e(k) - e(k+1))
261  endif
262  endif
263  enddo
264  else
265  do k=1,nz ; tr(i,j,k) = landval ; enddo
266  endif ; enddo ! i-loop
267  enddo ! j-loop
268  endif
269 
270  deallocate(tr_in) ; deallocate(tr_1d) ; deallocate(z_edges)
271  deallocate(wt) ; deallocate(z1) ; deallocate(z2)
272 
273  tracer_z_init = .true.
274 
275 end function tracer_z_init
276 
277 !> This subroutine reads the vertical coordinate data for a field from a NetCDF file.
278 !! It also might read the missing value attribute for that same field.
279 subroutine read_z_edges(filename, tr_name, z_edges, nz_out, has_edges, &
280  use_missing, missing, scale)
281  character(len=*), intent(in) :: filename !< The name of the file to read from.
282  character(len=*), intent(in) :: tr_name !< The name of the tracer in the file.
283  real, dimension(:), allocatable, &
284  intent(out) :: z_edges !< The depths of the vertical edges of the tracer array
285  integer, intent(out) :: nz_out !< The number of vertical layers in the tracer array
286  logical, intent(out) :: has_edges !< If true the values in z_edges are the edges of the
287  !! tracer cells, otherwise they are the cell centers
288  logical, intent(inout) :: use_missing !< If false on input, see whether the tracer has a
289  !! missing value, and if so return true
290  real, intent(inout) :: missing !< The missing value, if one has been found
291  real, intent(in) :: scale !< A scaling factor for z_edges into new units.
292 
293  ! This subroutine reads the vertical coordinate data for a field from a
294  ! NetCDF file. It also might read the missing value attribute for that same field.
295  character(len=32) :: mdl
296  character(len=120) :: dim_name, edge_name, tr_msg, dim_msg
297  logical :: monotonic
298  integer :: ncid, status, intid, tr_id, layid, k
299  integer :: nz_edge, ndim, tr_dim_ids(NF90_MAX_VAR_DIMS)
300 
301  mdl = "MOM_tracer_Z_init read_Z_edges: "
302  tr_msg = trim(tr_name)//" in "//trim(filename)
303 
304  status = nf90_open(filename, nf90_nowrite, ncid)
305  if (status /= nf90_noerr) then
306  call mom_error(warning,mdl//" Difficulties opening "//trim(filename)//&
307  " - "//trim(nf90_strerror(status)))
308  nz_out = -1 ; return
309  endif
310 
311  status = nf90_inq_varid(ncid, tr_name, tr_id)
312  if (status /= nf90_noerr) then
313  call mom_error(warning,mdl//" Difficulties finding variable "//&
314  trim(tr_msg)//" - "//trim(nf90_strerror(status)))
315  nz_out = -1 ; status = nf90_close(ncid) ; return
316  endif
317  status = nf90_inquire_variable(ncid, tr_id, ndims=ndim, dimids=tr_dim_ids)
318  if (status /= nf90_noerr) then
319  call mom_error(warning,mdl//" cannot inquire about "//trim(tr_msg))
320  elseif ((ndim < 3) .or. (ndim > 4)) then
321  call mom_error(warning,mdl//" "//trim(tr_msg)//&
322  " has too many or too few dimensions.")
323  nz_out = -1 ; status = nf90_close(ncid) ; return
324  endif
325 
326  if (.not.use_missing) then
327  ! Try to find the missing value from the dataset.
328  status = nf90_get_att(ncid, tr_id, "missing_value", missing)
329  if (status /= nf90_noerr) use_missing = .true.
330  endif
331 
332  ! Get the axis name and length.
333  status = nf90_inquire_dimension(ncid, tr_dim_ids(3), dim_name, len=nz_out)
334  if (status /= nf90_noerr) then
335  call mom_error(warning,mdl//" cannot inquire about dimension(3) of "//&
336  trim(tr_msg))
337  endif
338 
339  dim_msg = trim(dim_name)//" in "//trim(filename)
340  status = nf90_inq_varid(ncid, dim_name, layid)
341  if (status /= nf90_noerr) then
342  call mom_error(warning,mdl//" Difficulties finding variable "//&
343  trim(dim_msg)//" - "//trim(nf90_strerror(status)))
344  nz_out = -1 ; status = nf90_close(ncid) ; return
345  endif
346  ! Find out if the Z-axis has an edges attribute
347  status = nf90_get_att(ncid, layid, "edges", edge_name)
348  if (status /= nf90_noerr) then
349  call mom_mesg(mdl//" "//trim(dim_msg)//&
350  " has no readable edges attribute - "//trim(nf90_strerror(status)))
351  has_edges = .false.
352  else
353  has_edges = .true.
354  status = nf90_inq_varid(ncid, edge_name, intid)
355  if (status /= nf90_noerr) then
356  call mom_error(warning,mdl//" Difficulties finding edge variable "//&
357  trim(edge_name)//" in "//trim(filename)//" - "//trim(nf90_strerror(status)))
358  has_edges = .false.
359  endif
360  endif
361 
362  nz_edge = nz_out ; if (has_edges) nz_edge = nz_out+1
363  allocate(z_edges(nz_edge)) ; z_edges(:) = 0.0
364 
365  if (nz_out < 1) return
366 
367  ! Read the right variable.
368  if (has_edges) then
369  dim_msg = trim(edge_name)//" in "//trim(filename)
370  status = nf90_get_var(ncid, intid, z_edges)
371  if (status /= nf90_noerr) then
372  call mom_error(warning,mdl//" Difficulties reading variable "//&
373  trim(dim_msg)//" - "//trim(nf90_strerror(status)))
374  nz_out = -1 ; status = nf90_close(ncid) ; return
375  endif
376  else
377  status = nf90_get_var(ncid, layid, z_edges)
378  if (status /= nf90_noerr) then
379  call mom_error(warning,mdl//" Difficulties reading variable "//&
380  trim(dim_msg)//" - "//trim(nf90_strerror(status)))
381  nz_out = -1 ; status = nf90_close(ncid) ; return
382  endif
383  endif
384 
385  status = nf90_close(ncid)
386  if (status /= nf90_noerr) call mom_error(warning, mdl// &
387  " Difficulties closing "//trim(filename)//" - "//trim(nf90_strerror(status)))
388 
389  ! z_edges should be montonically decreasing with our sign convention.
390  ! Change the sign sign convention if it looks like z_edges is increasing.
391  if (z_edges(1) < z_edges(2)) then
392  do k=1,nz_edge ; z_edges(k) = -z_edges(k) ; enddo
393  endif
394  ! Check that z_edges is now monotonically decreasing.
395  monotonic = .true.
396  do k=2,nz_edge ; if (z_edges(k) >= z_edges(k-1)) monotonic = .false. ; enddo
397  if (.not.monotonic) &
398  call mom_error(warning,mdl//" "//trim(dim_msg)//" is not monotonic.")
399 
400  if (scale /= 1.0) then ; do k=1,nz_edge ; z_edges(k) = scale*z_edges(k) ; enddo ; endif
401 
402 end subroutine read_z_edges
403 
404 !### `find_overlap` and `find_limited_slope` were previously part of
405 ! MOM_diag_to_Z.F90, and are nearly identical to `find_overlap` in
406 ! `midas_vertmap.F90` with some slight differences. We keep it here for
407 ! reproducibility, but the two should be merged at some point
408 
409 !> Determines the layers bounded by interfaces e that overlap
410 !! with the depth range between Z_top and Z_bot, and the fractional weights
411 !! of each layer. It also calculates the normalized relative depths of the range
412 !! of each layer that overlaps that depth range.
413 
414 ! ### TODO: Merge with midas_vertmap.F90:find_overlap()
415 subroutine find_overlap(e, Z_top, Z_bot, k_max, k_start, k_top, k_bot, wt, z1, z2)
416  real, dimension(:), intent(in) :: e !< Column interface heights, in arbitrary units.
417  real, intent(in) :: Z_top !< Top of range being mapped to, in the units of e.
418  real, intent(in) :: Z_bot !< Bottom of range being mapped to, in the units of e.
419  integer, intent(in) :: k_max !< Number of valid layers.
420  integer, intent(in) :: k_start !< Layer at which to start searching.
421  integer, intent(inout) :: k_top !< Indices of top layers that overlap with the depth
422  !! range.
423  integer, intent(inout) :: k_bot !< Indices of bottom layers that overlap with the
424  !! depth range.
425  real, dimension(:), intent(out) :: wt !< Relative weights of each layer from k_top to k_bot.
426  real, dimension(:), intent(out) :: z1 !< Depth of the top limits of the part of
427  !! a layer that contributes to a depth level, relative to the cell center and normalized
428  !! by the cell thickness [nondim]. Note that -1/2 <= z1 < z2 <= 1/2.
429  real, dimension(:), intent(out) :: z2 !< Depths of the bottom limit of the part of
430  !! a layer that contributes to a depth level, relative to the cell center and normalized
431  !! by the cell thickness [nondim]. Note that -1/2 <= z1 < z2 <= 1/2.
432  ! Local variables
433  real :: Ih, e_c, tot_wt, I_totwt
434  integer :: k
435 
436  do k=k_start,k_max ; if (e(k+1)<z_top) exit ; enddo
437  k_top = k
438  if (k>k_max) return
439 
440  ! Determine the fractional weights of each layer.
441  ! Note that by convention, e and Z_int decrease with increasing k.
442  if (e(k+1)<=z_bot) then
443  wt(k) = 1.0 ; k_bot = k
444  ih = 0.0 ; if (e(k) /= e(k+1)) ih = 1.0 / (e(k)-e(k+1))
445  e_c = 0.5*(e(k)+e(k+1))
446  z1(k) = (e_c - min(e(k),z_top)) * ih
447  z2(k) = (e_c - z_bot) * ih
448  else
449  wt(k) = min(e(k),z_top) - e(k+1) ; tot_wt = wt(k) ! These are always > 0.
450  if (e(k) /= e(k+1)) then
451  z1(k) = (0.5*(e(k)+e(k+1)) - min(e(k), z_top)) / (e(k)-e(k+1))
452  else ; z1(k) = -0.5 ; endif
453  z2(k) = 0.5
454  k_bot = k_max
455  do k=k_top+1,k_max
456  if (e(k+1)<=z_bot) then
457  k_bot = k
458  wt(k) = e(k) - z_bot ; z1(k) = -0.5
459  if (e(k) /= e(k+1)) then
460  z2(k) = (0.5*(e(k)+e(k+1)) - z_bot) / (e(k)-e(k+1))
461  else ; z2(k) = 0.5 ; endif
462  else
463  wt(k) = e(k) - e(k+1) ; z1(k) = -0.5 ; z2(k) = 0.5
464  endif
465  tot_wt = tot_wt + wt(k) ! wt(k) is always > 0.
466  if (k>=k_bot) exit
467  enddo
468 
469  i_totwt = 1.0 / tot_wt
470  do k=k_top,k_bot ; wt(k) = i_totwt*wt(k) ; enddo
471  endif
472 
473 end subroutine find_overlap
474 
475 !> This subroutine determines a limited slope for val to be advected with
476 !! a piecewise limited scheme.
477 ! ### TODO: Merge with midas_vertmap.F90:find_limited_slope()
478 subroutine find_limited_slope(val, e, slope, k)
479  real, dimension(:), intent(in) :: val !< A column of values that are being interpolated.
480  real, dimension(:), intent(in) :: e !< Column interface heights in arbitrary units
481  real, intent(out) :: slope !< Normalized slope in the intracell distribution of val.
482  integer, intent(in) :: k !< Layer whose slope is being determined.
483  ! Local variables
484  real :: d1, d2 ! Thicknesses in the units of e.
485 
486  d1 = 0.5*(e(k-1)-e(k+1)) ; d2 = 0.5*(e(k)-e(k+2))
487  if (((val(k)-val(k-1)) * (val(k)-val(k+1)) >= 0.0) .or. (d1*d2 <= 0.0)) then
488  slope = 0.0 ! ; curvature = 0.0
489  else
490  slope = (d1**2*(val(k+1) - val(k)) + d2**2*(val(k) - val(k-1))) * &
491  ((e(k) - e(k+1)) / (d1*d2*(d1+d2)))
492  ! slope = 0.5*(val(k+1) - val(k-1))
493  ! This is S.J. Lin's form of the PLM limiter.
494  slope = sign(1.0,slope) * min(abs(slope), &
495  2.0*(max(val(k-1),val(k),val(k+1)) - val(k)), &
496  2.0*(val(k) - min(val(k-1),val(k),val(k+1))))
497  ! curvature = 0.0
498  endif
499 
500 end subroutine find_limited_slope
501 
502 
503 end module mom_tracer_z_init
mom_tracer_z_init::find_limited_slope
subroutine find_limited_slope(val, e, slope, k)
This subroutine determines a limited slope for val to be advected with a piecewise limited scheme.
Definition: MOM_tracer_Z_init.F90:479
mom_tracer_z_init::read_z_edges
subroutine read_z_edges(filename, tr_name, z_edges, nz_out, has_edges, use_missing, missing, scale)
This subroutine reads the vertical coordinate data for a field from a NetCDF file....
Definition: MOM_tracer_Z_init.F90:281
mom_error_handler::mom_mesg
subroutine, public mom_mesg(message, verb, all_print)
This provides a convenient interface for writing an informative comment.
Definition: MOM_error_handler.F90:53
mom_io
This module contains I/O framework code.
Definition: MOM_io.F90:2
mom_tracer_z_init
Used to initialize tracers from a depth- (or z*-) space file.
Definition: MOM_tracer_Z_init.F90:2
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
mom_io::mom_read_data
Read a data field from a file.
Definition: MOM_io.F90:74
mom_tracer_z_init::tracer_z_init
logical function, public tracer_z_init(tr, h, filename, tr_name, G, US, missing_val, land_val)
This function initializes a tracer by reading a Z-space file, returning .true. if this appears to hav...
Definition: MOM_tracer_Z_init.F90:31
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
mom_error_handler::is_root_pe
logical function, public is_root_pe()
This returns .true. if the current PE is the root PE.
Definition: MOM_error_handler.F90:44
mom_error_handler::mom_error
subroutine, public mom_error(level, message, all_print)
This provides a convenient interface for writing an mpp_error message with run-time filter based on a...
Definition: MOM_error_handler.F90:72
mom_tracer_z_init::find_overlap
subroutine find_overlap(e, Z_top, Z_bot, k_max, k_start, k_top, k_bot, wt, z1, z2)
Determines the layers bounded by interfaces e that overlap with the depth range between Z_top and Z_b...
Definition: MOM_tracer_Z_init.F90:416
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:26