MOM6
Functions/Subroutines
mom_tracer_z_init Module Reference

Detailed Description

Used to initialize tracers from a depth- (or z*-) space file.

Functions/Subroutines

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 have been successful, and false otherwise. More...
 
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. It also might read the missing value attribute for that same field. More...
 
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_bot, and the fractional weights of each layer. It also calculates the normalized relative depths of the range of each layer that overlaps that depth range. More...
 
subroutine find_limited_slope (val, e, slope, k)
 This subroutine determines a limited slope for val to be advected with a piecewise limited scheme. More...
 

Function/Subroutine Documentation

◆ find_limited_slope()

subroutine mom_tracer_z_init::find_limited_slope ( real, dimension(:), intent(in)  val,
real, dimension(:), intent(in)  e,
real, intent(out)  slope,
integer, intent(in)  k 
)
private

This subroutine determines a limited slope for val to be advected with a piecewise limited scheme.

Parameters
[in]valA column of values that are being interpolated.
[in]eColumn interface heights in arbitrary units
[out]slopeNormalized slope in the intracell distribution of val.
[in]kLayer whose slope is being determined.

Definition at line 479 of file MOM_tracer_Z_init.F90.

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 

Referenced by tracer_z_init().

Here is the caller graph for this function:

◆ find_overlap()

subroutine mom_tracer_z_init::find_overlap ( real, dimension(:), intent(in)  e,
real, intent(in)  Z_top,
real, intent(in)  Z_bot,
integer, intent(in)  k_max,
integer, intent(in)  k_start,
integer, intent(inout)  k_top,
integer, intent(inout)  k_bot,
real, dimension(:), intent(out)  wt,
real, dimension(:), intent(out)  z1,
real, dimension(:), intent(out)  z2 
)
private

Determines the layers bounded by interfaces e that overlap with the depth range between Z_top and Z_bot, and the fractional weights of each layer. It also calculates the normalized relative depths of the range of each layer that overlaps that depth range.

Parameters
[in]eColumn interface heights, in arbitrary units.
[in]z_topTop of range being mapped to, in the units of e.
[in]z_botBottom of range being mapped to, in the units of e.
[in]k_maxNumber of valid layers.
[in]k_startLayer at which to start searching.
[in,out]k_topIndices of top layers that overlap with the depth range.
[in,out]k_botIndices of bottom layers that overlap with the depth range.
[out]wtRelative weights of each layer from k_top to k_bot.
[out]z1Depth of the top limits of the part of a layer that contributes to a depth level, relative to the cell center and normalized by the cell thickness [nondim]. Note that -1/2 <= z1 < z2 <= 1/2.
[out]z2Depths of the bottom limit of the part of a layer that contributes to a depth level, relative to the cell center and normalized by the cell thickness [nondim]. Note that -1/2 <= z1 < z2 <= 1/2.

Definition at line 416 of file MOM_tracer_Z_init.F90.

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 

Referenced by tracer_z_init().

Here is the caller graph for this function:

◆ read_z_edges()

subroutine mom_tracer_z_init::read_z_edges ( character(len=*), intent(in)  filename,
character(len=*), intent(in)  tr_name,
real, dimension(:), intent(out), allocatable  z_edges,
integer, intent(out)  nz_out,
logical, intent(out)  has_edges,
logical, intent(inout)  use_missing,
real, intent(inout)  missing,
real, intent(in)  scale 
)
private

This subroutine reads the vertical coordinate data for a field from a NetCDF file. It also might read the missing value attribute for that same field.

Parameters
[in]filenameThe name of the file to read from.
[in]tr_nameThe name of the tracer in the file.
[out]z_edgesThe depths of the vertical edges of the tracer array
[out]nz_outThe number of vertical layers in the tracer array
[out]has_edgesIf true the values in z_edges are the edges of the tracer cells, otherwise they are the cell centers
[in,out]use_missingIf false on input, see whether the tracer has a missing value, and if so return true
[in,out]missingThe missing value, if one has been found
[in]scaleA scaling factor for z_edges into new units.

Definition at line 281 of file MOM_tracer_Z_init.F90.

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 

References mom_error_handler::mom_error(), and mom_error_handler::mom_mesg().

Referenced by tracer_z_init().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ tracer_z_init()

logical function, public mom_tracer_z_init::tracer_z_init ( real, dimension(szi_(g),szj_(g),szk_(g)), intent(out)  tr,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h,
character(len=*), intent(in)  filename,
character(len=*), intent(in)  tr_name,
type(ocean_grid_type), intent(in)  G,
type(unit_scale_type), intent(in)  US,
real, intent(in), optional  missing_val,
real, intent(in), optional  land_val 
)

This function initializes a tracer by reading a Z-space file, returning .true. if this appears to have been successful, and false otherwise.

Returns
A return code indicating if the initialization has been successful
Parameters
[in]gThe ocean's grid structure
[in]usA dimensional unit scaling type
[out]trThe tracer to initialize
[in]hLayer thicknesses [H ~> m or kg m-2]
[in]filenameThe name of the file to read from
[in]tr_nameThe name of the tracer in the file
[in]missing_valThe missing value for the tracer
[in]land_valA value to use to fill in land points

Definition at line 31 of file MOM_tracer_Z_init.F90.

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 

References find_limited_slope(), find_overlap(), mom_error_handler::mom_error(), and read_z_edges().

Referenced by mom_ocmip2_cfc::init_tracer_cfc(), ideal_age_example::initialize_ideal_age_tracer(), and oil_tracer::initialize_oil_tracer().

Here is the call graph for this function:
Here is the caller graph for this function: