MOM_horizontal_regridding.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!> Horizontal interpolation
7
8use mom_debugging, only : hchksum
9use mom_coms, only : max_across_pes, min_across_pes, sum_across_pes, broadcast
10use mom_coms, only : reproducing_sum
11use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, clock_loop
12use mom_domains, only : pass_var
13use mom_error_handler, only : mom_mesg, mom_error, fatal, warning, is_root_pe
16use mom_file_parser, only : get_param, log_param, log_version, param_file_type
17use mom_grid, only : ocean_grid_type
18use mom_interpolate, only : time_interp_external
19use mom_interp_infra, only : run_horiz_interp, build_horiz_interp_weights
20use mom_interp_infra, only : horiz_interp_type, horizontal_interp_init
22use mom_interp_infra, only : external_field
23use mom_time_manager, only : time_type
24use mom_io, only : axis_info, get_axis_info, get_var_axes_info, mom_read_data
25use mom_io, only : read_attribute, read_variable
26
27implicit none ; private
28
29#include <MOM_memory.h>
30
31public :: horiz_interp_and_extrap_tracer, mystats, homogenize_field
32
33!> Extrapolate and interpolate data
34interface horiz_interp_and_extrap_tracer
37end interface
38
39! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
40! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
41! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
42! vary with the Boussinesq approximation, the Boussinesq variant is given first.
43! The functions in this module work with variables with arbitrary units, in which case the
44! arbitrary rescaled units are indicated with [A ~> a], while the unscaled units are just [a].
45
46contains
47
48!> Write to the terminal some basic statistics about the k-th level of an array
49subroutine mystats(array, missing, G, k, mesg, unscale, full_halo)
50 type(ocean_grid_type), intent(in) :: g !< Ocean grid type
51 real, dimension(SZI_(G),SZJ_(G)), &
52 intent(in) :: array !< input array in arbitrary units [A ~> a]
53 real, intent(in) :: missing !< missing value in arbitrary units [A ~> a]
54 integer, intent(in) :: k !< Level to calculate statistics for
55 character(len=*), intent(in) :: mesg !< Label to use in message
56 real, optional, intent(in) :: unscale !< A scaling factor for output that countacts
57 !! any internal dimesional scaling [a A-1 ~> 1]
58 logical, optional, intent(in) :: full_halo !< If present and true, test values on the whole
59 !! array rather than just the computational domain.
60 ! Local variables
61 real :: mina ! Minimum value in the array in the arbitrary units of the input array [A ~> a]
62 real :: maxa ! Maximum value in the array in the arbitrary units of the input array [A ~> a]
63 real :: scl ! A factor for undoing any scaling of the array statistics for output [a A-1 ~> 1]
64 integer :: i, j, is, ie, js, je
65 logical :: found
66 character(len=120) :: lmesg
67
68 scl = 1.0 ; if (present(unscale)) scl = unscale
69 mina = 9.e24 / scl ; maxa = -9.e24 / scl ; found = .false.
70
71 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
72 if (present(full_halo)) then ; if (full_halo) then
73 is = g%isd ; ie = g%ied ; js = g%jsd ; je = g%jed
74 endif ; endif
75
76 do j=js,je ; do i=is,ie
77 if (array(i,j) /= array(i,j)) stop 'Nan!'
78 if (abs(array(i,j)-missing) > 1.e-6*abs(missing)) then
79 if (found) then
80 mina = min(mina, array(i,j))
81 maxa = max(maxa, array(i,j))
82 else
83 found = .true.
84 mina = array(i,j)
85 maxa = array(i,j)
86 endif
87 endif
88 enddo ; enddo
89 call min_across_pes(mina)
90 call max_across_pes(maxa)
91 if (is_root_pe()) then
92 write(lmesg(1:120),'(2(a,es12.4),a,I0,1x,a)') &
93 'init_from_Z: min=',mina*scl,' max=',maxa*scl,' Level=',k,trim(mesg)
94 call mom_mesg(lmesg,2)
95 endif
96
97end subroutine mystats
98
99!> Use ICE-9 algorithm to populate points (fill=1) with valid data (good=1). If no information
100!! is available, use a previous guess (prev). Optionally (smooth) blend the filled points to
101!! achieve a more desirable result.
102subroutine fill_miss_2d(aout, good, fill, prev, G, acrit, num_pass, relc, debug, answer_date)
103 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
104 real, dimension(SZI_(G),SZJ_(G)), &
105 intent(inout) :: aout !< The array with missing values to fill [arbitrary]
106 real, dimension(SZI_(G),SZJ_(G)), &
107 intent(in) :: good !< Valid data mask for incoming array
108 !! (1==good data; 0==missing data) [nondim].
109 real, dimension(SZI_(G),SZJ_(G)), &
110 intent(in) :: fill !< Same shape array of points which need
111 !! filling (1==fill;0==dont fill) [nondim]
112 real, dimension(SZI_(G),SZJ_(G)), &
113 intent(in) :: prev !< First guess where isolated holes exist [arbitrary]
114 real, intent(in) :: acrit !< A minimal value for deltas between iterations that
115 !! determines when the smoothing has converged [arbitrary].
116 integer, optional, intent(in) :: num_pass !< The maximum number of iterations
117 real, optional, intent(in) :: relc !< A relaxation coefficient for Laplacian [nondim]
118 logical, optional, intent(in) :: debug !< If true, write verbose debugging messages.
119 integer, optional, intent(in) :: answer_date !< The vintage of the expressions in the code.
120 !! Dates before 20190101 give the same answers
121 !! as the code did in late 2018, while later versions
122 !! add parentheses for rotational symmetry.
123
124 real, dimension(SZI_(G),SZJ_(G)) :: a_filled ! The aout with missing values filled in [arbitrary]
125 real, dimension(SZI_(G),SZJ_(G)) :: a_chg ! The change in aout due to an iteration of smoothing [arbitrary]
126 real, dimension(SZI_(G),SZJ_(G)) :: fill_pts ! 1 for points that still need to be filled [nondim]
127 real, dimension(SZI_(G),SZJ_(G)) :: good_ ! The values that are valid for the current iteration [nondim]
128 real, dimension(SZI_(G),SZJ_(G)) :: good_new ! The values of good_ to use for the next iteration [nondim]
129
130 real :: east, west, north, south ! Valid neighboring values or 0 for invalid values [arbitrary]
131 real :: ge, gw, gn, gs ! Flags set to 0 or 1 indicating which neighbors have valid values [nondim]
132 real :: ngood ! The number of valid values in neighboring points [nondim]
133 real :: nfill ! The remaining number of points to fill [nondim]
134 real :: nfill_prev ! The previous value of nfill [nondim]
135 character(len=256) :: mesg ! The text of an error message
136 integer :: i, j, k
137 integer, parameter :: num_pass_default = 10000
138 real, parameter :: relc_default = 0.25 ! The default relaxation coefficient [nondim]
139
140 integer :: npass ! The maximum number of passes of the Laplacian smoother
141 integer :: is, ie, js, je
142 real :: relax_coeff ! The grid-scale Laplacian relaxation coefficient per timestep [nondim]
143 real :: ares ! The maximum magnitude change in aout [A]
144 logical :: debug_it, ans_2018
145
146 debug_it=.false.
147 if (PRESENT(debug)) debug_it=debug
148
149 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
150
151 npass = num_pass_default
152 if (PRESENT(num_pass)) npass = num_pass
153
154 relax_coeff = relc_default
155 if (PRESENT(relc)) relax_coeff = relc
156
157 ans_2018 = .true. ; if (PRESENT(answer_date)) ans_2018 = (answer_date < 20190101)
158
159 fill_pts(:,:) = fill(:,:)
160
161 nfill = sum(fill(is:ie,js:je))
162 call sum_across_pes(nfill)
163
164 nfill_prev = nfill
165 good_(:,:) = good(:,:)
166 a_chg(:,:) = 0.0
167
168 do while (nfill > 0.0)
169
170 call pass_var(good_,g%Domain)
171 call pass_var(aout,g%Domain)
172
173 a_filled(:,:) = aout(:,:)
174 good_new(:,:) = good_(:,:)
175
176 do j=js,je ; do i=is,ie
177
178 if (good_(i,j) == 1.0 .or. fill(i,j) == 0.) cycle
179
180 ge=good_(i+1,j) ; gw=good_(i-1,j)
181 gn=good_(i,j+1) ; gs=good_(i,j-1)
182 east=0.0 ; west=0.0 ; north=0.0 ; south=0.0
183 if (ge == 1.0) east = aout(i+1,j)*ge
184 if (gw == 1.0) west = aout(i-1,j)*gw
185 if (gn == 1.0) north = aout(i,j+1)*gn
186 if (gs == 1.0) south = aout(i,j-1)*gs
187
188 if (ans_2018) then
189 ngood = ge+gw+gn+gs
190 else
191 ngood = (ge+gw) + (gn+gs)
192 endif
193 if (ngood > 0.) then
194 if (ans_2018) then
195 a_filled(i,j) = (east+west+north+south)/ngood
196 else
197 a_filled(i,j) = ((east+west) + (north+south))/ngood
198 endif
199 fill_pts(i,j) = 0.0
200 good_new(i,j) = 1.0
201 endif
202 enddo ; enddo
203
204 aout(is:ie,js:je) = a_filled(is:ie,js:je)
205 good_(is:ie,js:je) = good_new(is:ie,js:je)
206 nfill_prev = nfill
207 nfill = sum(fill_pts(is:ie,js:je))
208 call sum_across_pes(nfill)
209
210 if (nfill == nfill_prev) then
211 do j=js,je ; do i=is,ie ; if (fill_pts(i,j) == 1.0) then
212 aout(i,j) = prev(i,j)
213 fill_pts(i,j) = 0.0
214 endif ; enddo ; enddo
215 endif
216
217 ! Determine the number of remaining points to fill globally.
218 nfill = sum(fill_pts(is:ie,js:je))
219 call sum_across_pes(nfill)
220
221 enddo ! while block for remaining points to fill.
222
223 ! Do Laplacian smoothing for the points that have been filled in.
224 do k=1,npass
225 call pass_var(aout,g%Domain)
226
227 a_chg(:,:) = 0.0
228 if (ans_2018) then
229 do j=js,je ; do i=is,ie
230 if (fill(i,j) == 1) then
231 east = max(good(i+1,j),fill(i+1,j)) ; west = max(good(i-1,j),fill(i-1,j))
232 north = max(good(i,j+1),fill(i,j+1)) ; south = max(good(i,j-1),fill(i,j-1))
233 a_chg(i,j) = relax_coeff*(south*aout(i,j-1)+north*aout(i,j+1) + &
234 west*aout(i-1,j)+east*aout(i+1,j) - &
235 (south+north+west+east)*aout(i,j))
236 endif
237 enddo ; enddo
238 else
239 do j=js,je ; do i=is,ie
240 if (fill(i,j) == 1) then
241 ge = max(good(i+1,j),fill(i+1,j)) ; gw = max(good(i-1,j),fill(i-1,j))
242 gn = max(good(i,j+1),fill(i,j+1)) ; gs = max(good(i,j-1),fill(i,j-1))
243 a_chg(i,j) = relax_coeff*( ((gs*aout(i,j-1) + gn*aout(i,j+1)) + &
244 (gw*aout(i-1,j) + ge*aout(i+1,j))) - &
245 ((gs + gn) + (gw + ge))*aout(i,j) )
246 endif
247 enddo ; enddo
248 endif
249
250 ares = 0.0
251 do j=js,je ; do i=is,ie
252 aout(i,j) = a_chg(i,j) + aout(i,j)
253 ares = max(ares, abs(a_chg(i,j)))
254 enddo ; enddo
255 call max_across_pes(ares)
256 if (ares <= acrit) exit
257 enddo
258
259 do j=js,je ; do i=is,ie
260 if (good_(i,j) == 0.0 .and. fill_pts(i,j) == 1.0) then
261 write(mesg,*) 'In fill_miss, fill, good,i,j= ',fill_pts(i,j),good_(i,j),i,j
262 call mom_error(warning, mesg, .true.)
263 call mom_error(fatal,"MOM_initialize: "// &
264 "fill is true and good is false after fill_miss, how did this happen? ")
265 endif
266 enddo ; enddo
267
268end subroutine fill_miss_2d
269
270!> Extrapolate and interpolate from a file record
271subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, recnum, G, tr_z, mask_z, &
272 z_in, z_edges_in, missing_value, scale, &
273 homogenize, m_to_Z, answers_2018, ongrid, tr_iter_tol, answer_date)
274
275 character(len=*), intent(in) :: filename !< Path to file containing tracer to be
276 !! interpolated.
277 character(len=*), intent(in) :: varnam !< Name of tracer in file.
278 integer, intent(in) :: recnum !< Record number of tracer to be read.
279 type(ocean_grid_type), intent(inout) :: G !< Grid object
280 real, allocatable, dimension(:,:,:), intent(out) :: tr_z
281 !< Allocatable tracer array on the horizontal
282 !! model grid and input-file vertical levels
283 !! in arbitrary units [A ~> a]
284 real, allocatable, dimension(:,:,:), intent(out) :: mask_z
285 !< Allocatable tracer mask array on the horizontal
286 !! model grid and input-file vertical levels [nondim]
287 real, allocatable, dimension(:), intent(out) :: z_in
288 !< Cell grid values for input data [Z ~> m]
289 real, allocatable, dimension(:), intent(out) :: z_edges_in
290 !< Cell grid edge values for input data [Z ~> m]
291 real, intent(out) :: missing_value !< The missing value in the returned array, scaled
292 !! to avoid accidentally having valid values match
293 !! missing values in the same units as tr_z [A ~> a]
294 real, intent(in) :: scale !< Scaling factor for tracer into the internal
295 !! units of the model for the units in the file [A a-1 ~> 1]
296 logical, optional, intent(in) :: homogenize !< If present and true, horizontally homogenize data
297 !! to produce perfectly "flat" initial conditions
298 real, optional, intent(in) :: m_to_Z !< A conversion factor from meters to the units
299 !! of depth [Z m-1 ~> 1]. If missing, G%bathyT must be in m.
300 logical, optional, intent(in) :: answers_2018 !< If true, use expressions that give the same
301 !! answers as the code did in late 2018. Otherwise
302 !! add parentheses for rotational symmetry.
303 logical, optional, intent(in) :: ongrid !< If true, then data are assumed to have been interpolated
304 !! to the model horizontal grid. In this case, only
305 !! extrapolation is performed by this routine
306 real, optional, intent(in) :: tr_iter_tol !< The tolerance for changes in tracer concentrations
307 !! between smoothing iterations that determines when to
308 !! stop iterating in the same units as tr_z [A ~> a]
309 integer, optional, intent(in) :: answer_date !< The vintage of the expressions in the code.
310 !! Dates before 20190101 give the same answers
311 !! as the code did in late 2018, while later versions
312 !! add parentheses for rotational symmetry.
313
314 ! Local variables
315 ! In the following comments, [A] is used to indicate the arbitrary, possibly rescaled units of the
316 ! input array while [a] indicates the unscaled (e.g., mks) units that can be used with the reproducing sums
317 real, dimension(:,:), allocatable :: tr_in !< A 2-d array for holding input data on its
318 !! native horizontal grid, with units that change
319 !! as the input data is interpreted [a] then [A ~> a]
320 real, dimension(:,:,:), allocatable :: tr_in_full !< A 3-d array for holding input data on the
321 !! model horizontal grid, with units that change
322 !! as the input data is interpreted [a] then [A ~> a]
323 real, dimension(:,:), allocatable :: tr_inp !< Native horizontal grid data extended to the poles
324 !! with units that change as the input data is
325 !! interpreted [a] then [A ~> a]
326 real, dimension(:,:), allocatable :: mask_in ! A 2-d mask for extended input grid [nondim]
327
328 real :: PI_180 ! A conversion factor from degrees to radians [radians degree-1]
329 integer :: id, jd, kd, jdp ! Input dataset data sizes
330 integer :: i, j, k
331 integer, dimension(4) :: start, count
332 real, dimension(:,:), allocatable :: x_in ! Input file longitudes [radians]
333 real, dimension(:,:), allocatable :: y_in ! Input file latitudes [radians]
334 real, dimension(:), allocatable :: lon_in ! The longitudes in the input file [degreesE] then [radians]
335 real, dimension(:), allocatable :: lat_in ! The latitudes in the input file [degreesN] then [radians]
336 real, dimension(:), allocatable :: lat_inp ! The input file latitudes expanded to the pole [degreesN] then [radians]
337 real :: max_lat ! The maximum latitude on the input grid [degreesN]
338 real :: pole ! The sum of tracer values at the pole [a]
339 real :: max_depth ! The maximum depth of the ocean [Z ~> m]
340 real :: npole ! The number of points contributing to the pole value [nondim]
341 real :: missing_val_in ! The missing value in the input field [a]
342 real :: roundoff ! The magnitude of roundoff, usually ~2e-16 [nondim]
343 real :: add_offset, scale_factor ! File-specific conversion factors [a] or [nondim]
344 integer :: ans_date ! The vintage of the expressions and order of arithmetic to use
345 logical :: found_attr
346 logical :: add_np
347 logical :: is_ongrid
348 type(horiz_interp_type) :: Interp
349 type(axis_info), dimension(4) :: axes_info ! Axis information used for regridding
350 integer :: is, ie, js, je ! compute domain indices
351 integer :: isg, ieg, jsg, jeg ! global extent
352 integer :: isd, ied, jsd, jed ! data domain indices
353 integer :: id_clock_read
354 logical :: debug=.false.
355 real :: I_scale ! The inverse of the scale factor for diagnostic output [a A-1 ~> 1]
356 real :: dtr_iter_stop ! The tolerance for changes in tracer concentrations between smoothing
357 ! iterations that determines when to stop iterating [A ~> a]
358 real, dimension(SZI_(G),SZJ_(G)) :: lon_out ! The longitude of points on the model grid [radians]
359 real, dimension(SZI_(G),SZJ_(G)) :: lat_out ! The latitude of points on the model grid [radians]
360 real, dimension(SZI_(G),SZJ_(G)) :: tr_out ! The tracer on the model grid [A ~> a]
361 real, dimension(SZI_(G),SZJ_(G)) :: mask_out ! The mask on the model grid [nondim]
362 real, dimension(SZI_(G),SZJ_(G)) :: good ! Where the data is valid, this is 1 [nondim]
363 real, dimension(SZI_(G),SZJ_(G)) :: fill ! 1 where the data needs to be filled in [nondim]
364 real, dimension(SZI_(G),SZJ_(G)) :: tr_outf ! The tracer concentrations after Ice-9 [A ~> a]
365 real, dimension(SZI_(G),SZJ_(G)) :: tr_prev ! The tracer concentrations in the layer above [A ~> a]
366 real, dimension(SZI_(G),SZJ_(G)) :: good2 ! 1 where the data is valid after Ice-9 [nondim]
367 real, dimension(SZI_(G),SZJ_(G)) :: fill2 ! 1 for points that still need to be filled after Ice-9 [nondim]
368
369 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
370 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
371 isg = g%isg ; ieg = g%ieg ; jsg = g%jsg ; jeg = g%jeg
372
373 id_clock_read = cpu_clock_id('(Initialize tracer from Z) read', grain=clock_loop)
374
375 is_ongrid = .false.
376 if (present(ongrid)) is_ongrid = ongrid
377
378 dtr_iter_stop = 1.0e-3*scale
379 if (present(tr_iter_tol)) dtr_iter_stop = tr_iter_tol
380
381 i_scale = 1.0 / scale
382
383 pi_180 = atan(1.0)/45.
384
385 ans_date = 20181231
386 if (present(answers_2018)) then ; if (.not.answers_2018) ans_date = 20190101 ; endif
387 if (present(answer_date)) ans_date = answer_date
388
389 ! Open NetCDF file and if present, extract data and spatial coordinate information
390 ! The convention adopted here requires that the data be written in (i,j,k) ordering.
391
392 call cpu_clock_begin(id_clock_read)
393
394 ! A note by MJH copied from elsewhere suggests that this code may be using the model connectivity
395 ! (e.g., reentrant or tripolar) but should use the dataset's connectivity instead.
396
397 call get_var_axes_info(trim(filename), trim(varnam), axes_info)
398
399 if (allocated(z_in)) deallocate(z_in)
400 if (allocated(z_edges_in)) deallocate(z_edges_in)
401 if (allocated(tr_z)) deallocate(tr_z)
402 if (allocated(mask_z)) deallocate(mask_z)
403
404 call get_axis_info(axes_info(1),ax_size=id)
405 call get_axis_info(axes_info(2),ax_size=jd)
406 call get_axis_info(axes_info(3),ax_size=kd)
407
408 allocate(lon_in(id), lat_in(jd), z_in(kd), z_edges_in(kd+1))
409 allocate(tr_z(isd:ied,jsd:jed,kd), source=0.0)
410 allocate(mask_z(isd:ied,jsd:jed,kd), source=0.0)
411
412 call get_axis_info(axes_info(1),ax_data=lon_in)
413 call get_axis_info(axes_info(2),ax_data=lat_in)
414 call get_axis_info(axes_info(3),ax_data=z_in)
415
416 call cpu_clock_end(id_clock_read)
417
418 if (present(m_to_z)) then ; do k=1,kd ; z_in(k) = m_to_z * z_in(k) ; enddo ; endif
419
420 add_np = .false.
421 jdp = jd
422 if (.not. is_ongrid) then
423 max_lat = maxval(lat_in)
424 if (max_lat < 90.0) then
425 ! Extrapolate the input data to the north pole using the northern-most latitude.
426 add_np = .true.
427 jdp = jd+1
428 allocate(lat_inp(jdp))
429 lat_inp(1:jd) = lat_in(:)
430 lat_inp(jd+1) = 90.0
431 deallocate(lat_in)
432 allocate(lat_in(1:jdp))
433 lat_in(:) = lat_inp(:)
434 endif
435 endif
436 ! construct level cell boundaries as the mid-point between adjacent centers
437
438 ! Set the I/O attributes
439 call read_attribute(trim(filename), "_FillValue", missing_val_in, &
440 varname=trim(varnam), found=found_attr)
441 if (.not. found_attr) call mom_error(fatal, &
442 "error finding missing value for " // trim(varnam) // &
443 " in file " // trim(filename) // " in hinterp_extrap")
444 missing_value = scale * missing_val_in
445
446 call read_attribute(trim(filename), "scale_factor", scale_factor, &
447 varname=trim(varnam), found=found_attr)
448 if (.not. found_attr) scale_factor = 1.
449
450 call read_attribute(trim(filename), "add_offset", add_offset, &
451 varname=trim(varnam), found=found_attr)
452 if (.not. found_attr) add_offset = 0.
453
454 z_edges_in(1) = 0.0
455 do k=2,kd
456 z_edges_in(k) = 0.5*(z_in(k-1)+z_in(k))
457 enddo
458 z_edges_in(kd+1) = 2.0*z_in(kd) - z_in(kd-1)
459
460 if (is_ongrid) then
461 allocate(tr_in(is:ie,js:je), source=0.0)
462 allocate(tr_in_full(is:ie,js:je,kd), source=0.0)
463 allocate(mask_in(is:ie,js:je), source=0.0)
464 else
465 call horizontal_interp_init()
466 lon_in = lon_in*pi_180
467 lat_in = lat_in*pi_180
468 allocate(x_in(id,jdp), y_in(id,jdp))
469 call meshgrid(lon_in, lat_in, x_in, y_in)
470 lon_out(:,:) = g%geoLonT(:,:)*pi_180
471 lat_out(:,:) = g%geoLatT(:,:)*pi_180
472 allocate(tr_in(id,jd), source=0.0)
473 allocate(tr_inp(id,jdp), source=0.0)
474 allocate(mask_in(id,jdp), source=0.0)
475 endif
476
477 max_depth = maxval(g%bathyT(:,:)) + g%Z_ref
478 call max_across_pes(max_depth)
479
480 if (z_edges_in(kd+1) < max_depth) z_edges_in(kd+1) = max_depth
481 roundoff = 3.0*epsilon(missing_val_in)
482
483 ! Loop through each data level and interpolate to model grid.
484 ! After interpolating, fill in points which will be needed to define the layers.
485
486 if (is_ongrid) then
487 start(1) = is+g%HI%idg_offset ; start(2) = js+g%HI%jdg_offset ; start(3) = 1
488 count(1) = ie-is+1 ; count(2) = je-js+1 ; count(3) = kd ; start(4) = 1 ; count(4) = 1
489 call mom_read_data(trim(filename), trim(varnam), tr_in_full, start, count, g%Domain)
490 endif
491
492 do k=1,kd
493 mask_in(:,:) = 0.0
494 tr_out(:,:) = 0.0
495
496 if (is_ongrid) then
497 tr_in(is:ie,js:je) = tr_in_full(is:ie,js:je,k)
498 do j=js,je
499 do i=is,ie
500 if (abs(tr_in(i,j)-missing_val_in) > abs(roundoff*missing_val_in)) then
501 mask_in(i,j) = 1.0
502 tr_in(i,j) = (tr_in(i,j)*scale_factor+add_offset) * scale
503 else
504 tr_in(i,j) = missing_value
505 endif
506 enddo
507 enddo
508
509 tr_out(is:ie,js:je) = tr_in(is:ie,js:je)
510
511 else ! .not.is_ongrid
512
513 start(:) = 1 ; start(3) = k
514 count(:) = 1 ; count(1) = id ; count(2) = jd
515 call read_variable(trim(filename), trim(varnam), tr_in, start=start, nread=count)
516
517 if (is_root_pe()) then
518 if (add_np) then
519 pole = 0.0 ; npole = 0.0
520 do i=1,id
521 if (abs(tr_in(i,jd)-missing_val_in) > abs(roundoff*missing_val_in)) then
522 pole = pole + tr_in(i,jd)
523 npole = npole + 1.0
524 endif
525 enddo
526 if (npole > 0) then
527 pole = pole / npole
528 else
529 pole = missing_val_in
530 endif
531 tr_inp(:,1:jd) = tr_in(:,:)
532 tr_inp(:,jdp) = pole
533 else
534 tr_inp(:,:) = tr_in(:,:)
535 endif
536 endif
537
538 call broadcast(tr_inp, id*jdp, blocking=.true.)
539
540 do j=1,jdp ; do i=1,id
541 if (abs(tr_inp(i,j)-missing_val_in) > abs(roundoff*missing_val_in)) then
542 mask_in(i,j) = 1.0
543 tr_inp(i,j) = (tr_inp(i,j)*scale_factor+add_offset) * scale
544 else
545 tr_inp(i,j) = missing_value
546 endif
547 enddo ; enddo
548
549 ! call fms routine horiz_interp to interpolate input level data to model horizontal grid
550 if (k == 1) then
551 call build_horiz_interp_weights(interp, x_in, y_in, lon_out(is:ie,js:je), lat_out(is:ie,js:je), &
552 interp_method='bilinear', src_modulo=.true.)
553 endif
554
555 if (debug) then
556 call mystats(tr_inp, missing_value, g, k, 'Tracer from file', unscale=i_scale, full_halo=.true.)
557 endif
558
559 call run_horiz_interp(interp, tr_inp, tr_out(is:ie,js:je), missing_value=missing_value)
560 endif ! End of .not.is_ongrid
561
562 mask_out(:,:) = 1.0
563 do j=js,je ; do i=is,ie
564 if (abs(tr_out(i,j)-missing_value) < abs(roundoff*missing_value)) mask_out(i,j) = 0.
565 enddo ; enddo
566
567 fill(:,:) = 0.0 ; good(:,:) = 0.0
568
569 do j=js,je ; do i=is,ie
570 if (mask_out(i,j) < 1.0) then
571 tr_out(i,j) = missing_value
572 else
573 good(i,j) = 1.0
574 endif
575 if ((g%mask2dT(i,j) == 1.0) .and. (z_edges_in(k) <= g%bathyT(i,j) + g%Z_ref) .and. &
576 (mask_out(i,j) < 1.0)) &
577 fill(i,j) = 1.0
578 enddo ; enddo
579
580 call pass_var(fill, g%Domain)
581 call pass_var(good, g%Domain)
582
583 if (debug) then
584 call mystats(tr_out, missing_value, g, k, 'variable from horiz_interp()', unscale=i_scale)
585 endif
586
587 ! Horizontally homogenize data to produce perfectly "flat" initial conditions
588 if (PRESENT(homogenize)) then ; if (homogenize) then
589 call homogenize_field(tr_out, g, tmp_scale=i_scale, weights=mask_out, answer_date=answer_date)
590 endif ; endif
591
592 ! tr_out contains input z-space data on the model grid with missing values
593 ! now fill in missing values using "ICE-nine" algorithm.
594 tr_outf(:,:) = tr_out(:,:)
595 if (k==1) tr_prev(:,:) = tr_outf(:,:)
596 good2(:,:) = good(:,:)
597 fill2(:,:) = fill(:,:)
598
599 call fill_miss_2d(tr_outf, good2, fill2, tr_prev, g, dtr_iter_stop, answer_date=ans_date)
600 if (debug) then
601 call mystats(tr_outf, missing_value, g, k, 'field from fill_miss_2d()', unscale=i_scale)
602 endif
603
604 tr_z(:,:,k) = tr_outf(:,:) * g%mask2dT(:,:)
605 mask_z(:,:,k) = good2(:,:) + fill2(:,:)
606
607 tr_prev(:,:) = tr_z(:,:,k)
608
609 if (debug) then
610 call hchksum(tr_prev, 'field after fill ', g%HI, unscale=i_scale)
611 endif
612
613 enddo ! kd
614
615 if (allocated(lat_inp)) deallocate(lat_inp)
616 deallocate(tr_in)
617 if (allocated(tr_inp)) deallocate(tr_inp)
618 if (allocated(tr_in_full)) deallocate(tr_in_full)
619
621
622!> Extrapolate and interpolate using a FMS time interpolation handle
623subroutine horiz_interp_and_extrap_tracer_fms_id(field, Time, G, tr_z, mask_z, &
624 z_in, z_edges_in, missing_value, scale, &
625 homogenize, spongeOngrid, m_to_Z, &
626 answers_2018, tr_iter_tol, answer_date, &
627 axes)
628
629 type(external_field), intent(in) :: field !< Handle for the time interpolated field
630 type(time_type), intent(in) :: Time !< A FMS time type
631 type(ocean_grid_type), intent(inout) :: G !< Grid object
632 real, allocatable, dimension(:,:,:), intent(out) :: tr_z
633 !< Allocatable tracer array on the horizontal
634 !! model grid and input-file vertical levels
635 !! in arbitrary units [A ~> a]
636 real, allocatable, dimension(:,:,:), intent(out) :: mask_z
637 !< Allocatable tracer mask array on the horizontal
638 !! model grid and input-file vertical levels [nondim]
639 real, allocatable, dimension(:), intent(out) :: z_in
640 !< Cell grid values for input data [Z ~> m]
641 real, allocatable, dimension(:), intent(out) :: z_edges_in
642 !< Cell grid edge values for input data [Z ~> m]
643 real, intent(out) :: missing_value !< The missing value in the returned array, scaled
644 !! to avoid accidentally having valid values match
645 !! missing values, in the same arbitrary units as tr_z [A ~> a]
646 real, intent(in) :: scale !< Scaling factor for tracer into the internal
647 !! units of the model [A a-1 ~> 1]
648 logical, optional, intent(in) :: homogenize !< If present and true, horizontally homogenize data
649 !! to produce perfectly "flat" initial conditions
650 logical, optional, intent(in) :: spongeOngrid !< If present and true, the sponge data are on the model grid
651 real, optional, intent(in) :: m_to_Z !< A conversion factor from meters to the units
652 !! of depth [Z m-1 ~> 1]. If missing, G%bathyT must be in m.
653 logical, optional, intent(in) :: answers_2018 !< If true, use expressions that give the same
654 !! answers as the code did in late 2018. Otherwise
655 !! add parentheses for rotational symmetry.
656 real, optional, intent(in) :: tr_iter_tol !< The tolerance for changes in tracer concentrations
657 !! between smoothing iterations that determines when to
658 !! stop iterating, in the same arbitrary units as tr_z [A ~> a]
659 integer, optional, intent(in) :: answer_date !< The vintage of the expressions in the code.
660 !! Dates before 20190101 give the same answers
661 !! as the code did in late 2018, while later versions
662 !! add parentheses for rotational symmetry.
663 type(axis_info), allocatable, dimension(:), optional, intent(inout) :: axes !< Axis types for the input data
664
665 ! Local variables
666 ! In the following comments, [A] is used to indicate the arbitrary, possibly rescaled units of the
667 ! input array while [a] indicates the unscaled (e.g., mks) units that can be used with the reproducing sums
668 real, dimension(:,:), allocatable :: tr_in !< A 2-d array for holding input data on its
669 !! native horizontal grid, with units that change
670 !! as the input data is interpreted [a] then [A ~> a]
671 real, dimension(:,:), allocatable :: tr_inp !< Native horizontal grid data extended to the poles
672 !! with units that change as the input data is
673 !! interpreted [a] then [A ~> a]
674 real, dimension(:,:,:), allocatable :: data_in !< A buffer for storing the full 3-d time-interpolated array
675 !! on the original grid [a]
676 real, dimension(:,:), allocatable :: mask_in !< A 2-d mask for extended input grid [nondim]
677
678 real :: PI_180 ! A conversion factor from degrees to radians [radians degree-1]
679 integer :: id, jd, kd, jdp ! Input dataset data sizes
680 integer :: i, j, k
681 real, dimension(:,:), allocatable :: x_in ! Input file longitudes [radians]
682 real, dimension(:,:), allocatable :: y_in ! Input file latitudes [radians]
683 real, dimension(:), allocatable :: lon_in ! The longitudes in the input file [degreesE] then [radians]
684 real, dimension(:), allocatable :: lat_in ! The latitudes in the input file [degreesN] then [radians]
685 real, dimension(:), allocatable :: lat_inp ! The input file latitudes expanded to the pole [degreesN] then [radians]
686 real :: max_lat ! The maximum latitude on the input grid [degreesN]
687 real :: pole ! The sum of tracer values at the pole [a]
688 real :: max_depth ! The maximum depth of the ocean [Z ~> m]
689 real :: npole ! The number of points contributing to the pole value [nondim]
690 real :: missing_val_in ! The missing value in the input field [a]
691 real :: roundoff ! The magnitude of roundoff, usually ~2e-16 [nondim]
692 logical :: add_np
693 type(horiz_interp_type) :: Interp
694 type(axis_info), dimension(4) :: axes_data
695 integer :: is, ie, js, je ! compute domain indices
696 integer :: isg, ieg, jsg, jeg ! global extent
697 integer :: isd, ied, jsd, jed ! data domain indices
698 integer :: id_clock_read
699 integer, dimension(4) :: fld_sz
700 logical :: debug=.false.
701 logical :: is_ongrid
702 integer :: ans_date ! The vintage of the expressions and order of arithmetic to use
703 real :: I_scale ! The inverse of the scale factor for diagnostic output [a A-1 ~> 1]
704 real :: dtr_iter_stop ! The tolerance for changes in tracer concentrations between smoothing
705 ! iterations that determines when to stop iterating [A ~> a]
706 real, dimension(SZI_(G),SZJ_(G)) :: lon_out ! The longitude of points on the model grid [radians]
707 real, dimension(SZI_(G),SZJ_(G)) :: lat_out ! The latitude of points on the model grid [radians]
708 real, dimension(SZI_(G),SZJ_(G)) :: tr_out ! The tracer on the model grid [A ~> a]
709 real, dimension(SZI_(G),SZJ_(G)) :: mask_out ! The mask on the model grid [nondim]
710 real, dimension(SZI_(G),SZJ_(G)) :: good ! Where the data is valid, this is 1 [nondim]
711 real, dimension(SZI_(G),SZJ_(G)) :: fill ! 1 where the data needs to be filled in [nondim]
712 real, dimension(SZI_(G),SZJ_(G)) :: tr_outf ! The tracer concentrations after Ice-9 [A ~> a]
713 real, dimension(SZI_(G),SZJ_(G)) :: tr_prev ! The tracer concentrations in the layer above [A ~> a]
714 real, dimension(SZI_(G),SZJ_(G)) :: good2 ! 1 where the data is valid after Ice-9 [nondim]
715 real, dimension(SZI_(G),SZJ_(G)) :: fill2 ! 1 for points that still need to be filled after Ice-9 [nondim]
716 integer :: turns
717 integer :: verbosity
718
719 turns = g%HI%turns
720
721 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
722 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
723 isg = g%isg ; ieg = g%ieg ; jsg = g%jsg ; jeg = g%jeg
724
725 id_clock_read = cpu_clock_id('(Initialize tracer from Z) read', grain=clock_loop)
726
727 dtr_iter_stop = 1.0e-3*scale
728 if (present(tr_iter_tol)) dtr_iter_stop = tr_iter_tol
729
730 i_scale = 1.0 / scale
731
732 pi_180 = atan(1.0)/45.
733
734 ans_date = 20181231
735 if (present(answers_2018)) then ; if (.not.answers_2018) ans_date = 20190101 ; endif
736 if (present(answer_date)) ans_date = answer_date
737
738 ! Open NetCDF file and if present, extract data and spatial coordinate information
739 ! The convention adopted here requires that the data be written in (i,j,k) ordering.
740
741 call cpu_clock_begin(id_clock_read)
742
743 if (present(axes) .and. allocated(axes)) then
744 call get_external_field_info(field, size=fld_sz, missing=missing_val_in)
745 axes_data = axes
746 else
747 call get_external_field_info(field, size=fld_sz, axes=axes_data, missing=missing_val_in)
748 if (present(axes)) then
749 allocate(axes(4))
750 axes = axes_data
751 endif
752 endif
753 missing_value = scale*missing_val_in
754
755 verbosity = mom_get_verbosity()
756
757 id = fld_sz(1) ; jd = fld_sz(2) ; kd = fld_sz(3)
758
759 is_ongrid = .false.
760 if (PRESENT(spongeongrid)) is_ongrid = spongeongrid
761 if (.not. is_ongrid) then
762 allocate(lon_in(id), lat_in(jd))
763 call get_axis_info(axes_data(1), ax_data=lon_in)
764 call get_axis_info(axes_data(2), ax_data=lat_in)
765 endif
766
767 allocate(z_in(kd), z_edges_in(kd+1))
768
769 allocate(tr_z(isd:ied,jsd:jed,kd), source=0.0)
770 allocate(mask_z(isd:ied,jsd:jed,kd), source=0.0)
771
772 call get_axis_info(axes_data(3), ax_data=z_in)
773
774 if (present(m_to_z)) then ; do k=1,kd ; z_in(k) = m_to_z * z_in(k) ; enddo ; endif
775
776 call cpu_clock_end(id_clock_read)
777
778 if (.not. is_ongrid) then
779 max_lat = maxval(lat_in)
780 add_np = .false.
781 if (max_lat < 90.0) then
782 ! Extrapolate the input data to the north pole using the northern-most latitude.
783 add_np = .true.
784 jdp = jd+1
785 allocate(lat_inp(jdp))
786 lat_inp(1:jd) = lat_in(:)
787 lat_inp(jd+1) = 90.0
788 deallocate(lat_in)
789 allocate(lat_in(1:jdp))
790 lat_in(:) = lat_inp(:)
791 else
792 jdp = jd
793 endif
794 call horizontal_interp_init()
795 lon_in = lon_in*pi_180
796 lat_in = lat_in*pi_180
797 allocate(x_in(id,jdp), y_in(id,jdp))
798 call meshgrid(lon_in, lat_in, x_in, y_in)
799 lon_out(:,:) = g%geoLonT(:,:)*pi_180
800 lat_out(:,:) = g%geoLatT(:,:)*pi_180
801 allocate(data_in(id,jd,kd), source=0.0)
802 allocate(tr_in(id,jd), source=0.0)
803 allocate(tr_inp(id,jdp), source=0.0)
804 allocate(mask_in(id,jdp), source=0.0)
805 else
806 allocate(data_in(isd:ied,jsd:jed,kd))
807 endif
808
809 ! Construct level cell boundaries as the mid-point between adjacent centers.
810 z_edges_in(1) = 0.0
811 do k=2,kd
812 z_edges_in(k) = 0.5*(z_in(k-1)+z_in(k))
813 enddo
814 z_edges_in(kd+1) = 2.0*z_in(kd) - z_in(kd-1)
815
816 max_depth = maxval(g%bathyT(:,:)) + g%Z_ref
817 call max_across_pes(max_depth)
818
819 if (z_edges_in(kd+1) < max_depth) z_edges_in(kd+1) = max_depth
820
821 ! roundoff = 3.0*EPSILON(missing_value)
822 roundoff = 1.e-4
823
824 if (.not.is_ongrid) then
825 if (is_root_pe()) &
826 call time_interp_external(field, time, data_in, verbose=(verbosity>5), turns=turns)
827
828 ! Loop through each data level and interpolate to model grid.
829 ! After interpolating, fill in points which will be needed to define the layers.
830 do k=1,kd
831 if (is_root_pe()) then
832 tr_in(1:id,1:jd) = data_in(1:id,1:jd,k)
833 if (add_np) then
834 pole = 0.0 ; npole = 0.0
835 do i=1,id
836 if (abs(tr_in(i,jd)-missing_val_in) > abs(roundoff*missing_val_in)) then
837 pole = pole + tr_in(i,jd)
838 npole = npole + 1.0
839 endif
840 enddo
841 if (npole > 0) then
842 pole = pole / npole
843 else
844 pole = missing_val_in
845 endif
846 tr_inp(:,1:jd) = tr_in(:,:)
847 tr_inp(:,jdp) = pole
848 else
849 tr_inp(:,:) = tr_in(:,:)
850 endif
851 endif
852
853 call broadcast(tr_inp, id*jdp, blocking=.true.)
854
855 mask_in(:,:) = 0.0
856
857 do j=1,jdp ; do i=1,id
858 if (abs(tr_inp(i,j)-missing_val_in) > abs(roundoff*missing_val_in)) then
859 mask_in(i,j) = 1.0
860 tr_inp(i,j) = tr_inp(i,j) * scale
861 else
862 tr_inp(i,j) = missing_value
863 endif
864 enddo ; enddo
865
866 ! call fms routine horiz_interp to interpolate input level data to model horizontal grid
867 if (k == 1) then
868 call build_horiz_interp_weights(interp, x_in, y_in, lon_out(is:ie,js:je), lat_out(is:ie,js:je), &
869 interp_method='bilinear', src_modulo=.true.)
870 endif
871
872 if (debug) then
873 call mystats(tr_inp, missing_value, g, k, 'Tracer from file', unscale=i_scale, full_halo=.true.)
874 endif
875
876 tr_out(:,:) = 0.0
877
878 call run_horiz_interp(interp, tr_inp, tr_out(is:ie,js:je), missing_value=missing_value)
879
880 mask_out(:,:) = 1.0
881 do j=js,je ; do i=is,ie
882 if (abs(tr_out(i,j)-missing_value) < abs(roundoff*missing_value)) mask_out(i,j) = 0.
883 enddo ; enddo
884
885 fill(:,:) = 0.0 ; good(:,:) = 0.0
886
887 do j=js,je ; do i=is,ie
888 if (mask_out(i,j) < 1.0) then
889 tr_out(i,j) = missing_value
890 else
891 good(i,j) = 1.0
892 endif
893 if ((g%mask2dT(i,j) == 1.0) .and. (z_edges_in(k) <= g%bathyT(i,j) + g%Z_ref) .and. &
894 (mask_out(i,j) < 1.0)) &
895 fill(i,j) = 1.0
896 enddo ; enddo
897 call pass_var(fill, g%Domain)
898 call pass_var(good, g%Domain)
899
900 if (debug) then
901 call mystats(tr_out, missing_value, g, k, 'variable from horiz_interp()', unscale=i_scale)
902 endif
903
904 ! Horizontally homogenize data to produce perfectly "flat" initial conditions
905 if (PRESENT(homogenize)) then ; if (homogenize) then
906 call homogenize_field(tr_out, g, tmp_scale=i_scale, weights=mask_out, answer_date=answer_date)
907 endif ; endif
908
909 ! tr_out contains input z-space data on the model grid with missing values
910 ! now fill in missing values using "ICE-nine" algorithm.
911 tr_outf(:,:) = tr_out(:,:)
912 if (k==1) tr_prev(:,:) = tr_outf(:,:)
913 good2(:,:) = good(:,:)
914 fill2(:,:) = fill(:,:)
915
916 call fill_miss_2d(tr_outf, good2, fill2, tr_prev, g, dtr_iter_stop, answer_date=ans_date)
917
918! if (debug) then
919! call hchksum(tr_outf, 'field from fill_miss_2d ', G%HI, unscale=I_scale)
920! call myStats(tr_outf, missing_value, G, k, 'field from fill_miss_2d()', unscale=I_scale)
921! endif
922
923 tr_z(:,:,k) = tr_outf(:,:) * g%mask2dT(:,:)
924 mask_z(:,:,k) = good2(:,:) + fill2(:,:)
925 tr_prev(:,:) = tr_z(:,:,k)
926
927 if (debug) then
928 call hchksum(tr_prev, 'field after fill ', g%HI, unscale=i_scale)
929 endif
930
931 enddo ! kd
932 else
933 call time_interp_external(field, time, data_in, verbose=(verbosity>5), turns=turns)
934 do k=1,kd
935 do j=js,je
936 do i=is,ie
937 tr_z(i,j,k) = data_in(i,j,k) * scale
938 if (ans_date >= 20190101) mask_z(i,j,k) = 1.
939 if (abs(tr_z(i,j,k)-missing_value) < abs(roundoff*missing_value)) mask_z(i,j,k) = 0.
940 enddo
941 enddo
942 enddo
943 endif
944
946
947!> Replace all values of a 2-d field with the weighted average over the valid points.
948subroutine homogenize_field(field, G, tmp_scale, weights, answer_date, wt_unscale)
949 type(ocean_grid_type), intent(inout) :: g !< Ocean grid type
950 real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: field !< The tracer on the model grid in arbitrary units [A ~> a]
951 real, optional, intent(in) :: tmp_scale !< A temporary rescaling factor for the
952 !! variable that is reversed in the
953 !! return value [a A-1 ~> 1]
954 real, dimension(SZI_(G),SZJ_(G)), &
955 optional, intent(in) :: weights !< The weights for the tracer in arbitrary units that
956 !! typically differ from those used by field [B ~> b]
957 integer, optional, intent(in) :: answer_date !< The vintage of the expressions in the code.
958 !! Dates before 20230101 use non-reproducing sums
959 !! in their averages, while later versions use
960 !! reproducing sums for rotational symmetry and
961 !! consistency across PE layouts.
962 real, optional, intent(in) :: wt_unscale !< A factor that undoes any dimensional scaling
963 !! of the weights so that they can be used with
964 !! reproducing sums [b B-1 ~> 1]
965
966 ! Local variables
967 ! In the following comments, [A] and [B] are used to indicate the arbitrary, possibly rescaled
968 ! units of the input field and the weighting array, while [a] and [b] indicate the corresponding
969 ! unscaled (e.g., mks) units that can be used with the reproducing sums
970 real, dimension(G%isc:G%iec, G%jsc:G%jec) :: field_for_sums ! The field times the weights [A B ~> a b]
971 real, dimension(G%isc:G%iec, G%jsc:G%jec) :: weight ! A copy of weights, if it is present, or the
972 ! tracer-point grid mask if it weights is absent [B ~> b]
973 real :: var_unscale ! The reciprocal of the scaling factor for the field and weights [a b A-1 B-1 ~> 1]
974 real :: wt_sum ! The sum of the weights, in [B ~> b]
975 real :: varsum ! The weighted sum of field being averaged [A B ~> a b]
976 real :: varavg ! The average of the field [A ~> a]
977 logical :: use_repro_sums ! If true, use reproducing sums.
978 integer :: i, j, is, ie, js, je
979
980 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
981
982 varavg = 0.0 ! This value will be used if wt_sum is 0.
983
984 use_repro_sums = .false. ; if (present(answer_date)) use_repro_sums = (answer_date >= 20230101)
985
986 if (present(weights)) then
987 do j=js,je ; do i=is,ie
988 weight(i,j) = weights(i,j)
989 enddo ; enddo
990 else
991 do j=js,je ; do i=is,ie
992 weight(i,j) = g%mask2dT(i,j)
993 enddo ; enddo
994 endif
995
996 if (use_repro_sums) then
997 var_unscale = 1.0 ; if (present(tmp_scale)) var_unscale = tmp_scale
998 if (present(wt_unscale)) var_unscale = wt_unscale * var_unscale
999
1000 do j=js,je ; do i=is,ie
1001 field_for_sums(i,j) = field(i,j) * weight(i,j)
1002 enddo ; enddo
1003
1004 wt_sum = reproducing_sum(weight, unscale=wt_unscale)
1005 if (abs(wt_sum) > 0.0) &
1006 varavg = reproducing_sum(field_for_sums, unscale=var_unscale) * (1.0 / wt_sum)
1007
1008 else ! Do the averages with order-dependent sums to reproduce older answers.
1009 wt_sum = 0 ; varsum = 0.
1010 do j=js,je ; do i=is,ie
1011 if (weight(i,j) > 0.0) then
1012 wt_sum = wt_sum + weight(i,j)
1013 varsum = varsum + field(i,j) * weight(i,j)
1014 endif
1015 enddo ; enddo
1016
1017 ! Note that these averages will not reproduce across PE layouts or grid rotation.
1018 call sum_across_pes(wt_sum)
1019 if (wt_sum > 0.0) then
1020 call sum_across_pes(varsum)
1021 varavg = varsum / wt_sum
1022 endif
1023
1024 endif
1025
1026 ! This seems like an unlikely case to ever be used, but it is needed to recreate previous behavior.
1027 if (present(tmp_scale)) then ; if (tmp_scale == 0.0) varavg = 0.0 ; endif
1028
1029 field(:,:) = varavg
1030
1031end subroutine homogenize_field
1032
1033
1034!> Create a 2d-mesh of grid coordinates from 1-d arrays.
1035subroutine meshgrid(x, y, x_T, y_T)
1036 real, dimension(:), intent(in) :: x !< input 1-dimensional vector [arbitrary]
1037 real, dimension(:), intent(in) :: y !< input 1-dimensional vector [arbitrary]
1038 real, dimension(size(x,1),size(y,1)), intent(inout) :: x_T !< output 2-dimensional array [arbitrary]
1039 real, dimension(size(x,1),size(y,1)), intent(inout) :: y_T !< output 2-dimensional array [arbitrary]
1040
1041 integer :: ni, nj, i, j
1042
1043 ni = size(x,1) ; nj = size(y,1)
1044
1045 do j=1,nj ; do i=1,ni
1046 x_t(i,j) = x(i)
1047 enddo ; enddo
1048
1049 do j=1,nj ; do i=1,ni
1050 y_t(i,j) = y(j)
1051 enddo ; enddo
1052
1053end subroutine meshgrid
1054