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 elseif (nfill == nfill_prev) then
216 call mom_error(warning, &
217 'Unable to fill missing points using either data at the same vertical level from a connected basin '//&
218 'or using a point from a previous vertical level. Make sure that the original data has some valid '//&
219 'data in all basins.', .true.)
220 write(mesg,*) 'nfill=',nfill
221 call mom_error(warning, mesg, .true.)
222 endif
223
224 ! Determine the number of remaining points to fill globally.
225 nfill = sum(fill_pts(is:ie,js:je))
226 call sum_across_pes(nfill)
227
228 enddo ! while block for remaining points to fill.
229
230 ! Do Laplacian smoothing for the points that have been filled in.
231 do k=1,npass
232 call pass_var(aout,g%Domain)
233
234 a_chg(:,:) = 0.0
235 if (ans_2018) then
236 do j=js,je ; do i=is,ie
237 if (fill(i,j) == 1) then
238 east = max(good(i+1,j),fill(i+1,j)) ; west = max(good(i-1,j),fill(i-1,j))
239 north = max(good(i,j+1),fill(i,j+1)) ; south = max(good(i,j-1),fill(i,j-1))
240 a_chg(i,j) = relax_coeff*(south*aout(i,j-1)+north*aout(i,j+1) + &
241 west*aout(i-1,j)+east*aout(i+1,j) - &
242 (south+north+west+east)*aout(i,j))
243 endif
244 enddo ; enddo
245 else
246 do j=js,je ; do i=is,ie
247 if (fill(i,j) == 1) then
248 ge = max(good(i+1,j),fill(i+1,j)) ; gw = max(good(i-1,j),fill(i-1,j))
249 gn = max(good(i,j+1),fill(i,j+1)) ; gs = max(good(i,j-1),fill(i,j-1))
250 a_chg(i,j) = relax_coeff*( ((gs*aout(i,j-1) + gn*aout(i,j+1)) + &
251 (gw*aout(i-1,j) + ge*aout(i+1,j))) - &
252 ((gs + gn) + (gw + ge))*aout(i,j) )
253 endif
254 enddo ; enddo
255 endif
256
257 ares = 0.0
258 do j=js,je ; do i=is,ie
259 aout(i,j) = a_chg(i,j) + aout(i,j)
260 ares = max(ares, abs(a_chg(i,j)))
261 enddo ; enddo
262 call max_across_pes(ares)
263 if (ares <= acrit) exit
264 enddo
265
266 do j=js,je ; do i=is,ie
267 if (good_(i,j) == 0.0 .and. fill_pts(i,j) == 1.0) then
268 write(mesg,*) 'In fill_miss, fill, good,i,j= ',fill_pts(i,j),good_(i,j),i,j
269 call mom_error(warning, mesg, .true.)
270 call mom_error(fatal,"MOM_initialize: "// &
271 "fill is true and good is false after fill_miss, how did this happen? ")
272 endif
273 enddo ; enddo
274
275end subroutine fill_miss_2d
276
277!> Extrapolate and interpolate from a file record
278subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, recnum, G, tr_z, mask_z, &
279 z_in, z_edges_in, missing_value, scale, &
280 homogenize, m_to_Z, answers_2018, ongrid, tr_iter_tol, answer_date)
281
282 character(len=*), intent(in) :: filename !< Path to file containing tracer to be
283 !! interpolated.
284 character(len=*), intent(in) :: varnam !< Name of tracer in file.
285 integer, intent(in) :: recnum !< Record number of tracer to be read.
286 type(ocean_grid_type), intent(inout) :: G !< Grid object
287 real, allocatable, dimension(:,:,:), intent(out) :: tr_z
288 !< Allocatable tracer array on the horizontal
289 !! model grid and input-file vertical levels
290 !! in arbitrary units [A ~> a]
291 real, allocatable, dimension(:,:,:), intent(out) :: mask_z
292 !< Allocatable tracer mask array on the horizontal
293 !! model grid and input-file vertical levels [nondim]
294 real, allocatable, dimension(:), intent(out) :: z_in
295 !< Cell grid values for input data [Z ~> m]
296 real, allocatable, dimension(:), intent(out) :: z_edges_in
297 !< Cell grid edge values for input data [Z ~> m]
298 real, intent(out) :: missing_value !< The missing value in the returned array, scaled
299 !! to avoid accidentally having valid values match
300 !! missing values in the same units as tr_z [A ~> a]
301 real, intent(in) :: scale !< Scaling factor for tracer into the internal
302 !! units of the model for the units in the file [A a-1 ~> 1]
303 logical, optional, intent(in) :: homogenize !< If present and true, horizontally homogenize data
304 !! to produce perfectly "flat" initial conditions
305 real, optional, intent(in) :: m_to_Z !< A conversion factor from meters to the units
306 !! of depth [Z m-1 ~> 1]. If missing, G%bathyT must be in m.
307 logical, optional, intent(in) :: answers_2018 !< If true, use expressions that give the same
308 !! answers as the code did in late 2018. Otherwise
309 !! add parentheses for rotational symmetry.
310 logical, optional, intent(in) :: ongrid !< If true, then data are assumed to have been interpolated
311 !! to the model horizontal grid. In this case, only
312 !! extrapolation is performed by this routine
313 real, optional, intent(in) :: tr_iter_tol !< The tolerance for changes in tracer concentrations
314 !! between smoothing iterations that determines when to
315 !! stop iterating in the same units as tr_z [A ~> a]
316 integer, optional, intent(in) :: answer_date !< The vintage of the expressions in the code.
317 !! Dates before 20190101 give the same answers
318 !! as the code did in late 2018, while later versions
319 !! add parentheses for rotational symmetry.
320
321 ! Local variables
322 ! In the following comments, [A] is used to indicate the arbitrary, possibly rescaled units of the
323 ! input array while [a] indicates the unscaled (e.g., mks) units that can be used with the reproducing sums
324 real, dimension(:,:), allocatable :: tr_in !< A 2-d array for holding input data on its
325 !! native horizontal grid, with units that change
326 !! as the input data is interpreted [a] then [A ~> a]
327 real, dimension(:,:,:), allocatable :: tr_in_full !< A 3-d array for holding input data on the
328 !! model horizontal grid, with units that change
329 !! as the input data is interpreted [a] then [A ~> a]
330 real, dimension(:,:), allocatable :: tr_inp !< Native horizontal grid data extended to the poles
331 !! with units that change as the input data is
332 !! interpreted [a] then [A ~> a]
333 real, dimension(:,:), allocatable :: mask_in ! A 2-d mask for extended input grid [nondim]
334
335 real :: PI_180 ! A conversion factor from degrees to radians [radians degree-1]
336 integer :: id, jd, kd, jdp ! Input dataset data sizes
337 integer :: i, j, k
338 integer, dimension(4) :: start, count
339 real, dimension(:,:), allocatable :: x_in ! Input file longitudes [radians]
340 real, dimension(:,:), allocatable :: y_in ! Input file latitudes [radians]
341 real, dimension(:), allocatable :: lon_in ! The longitudes in the input file [degreesE] then [radians]
342 real, dimension(:), allocatable :: lat_in ! The latitudes in the input file [degreesN] then [radians]
343 real, dimension(:), allocatable :: lat_inp ! The input file latitudes expanded to the pole [degreesN] then [radians]
344 real :: max_lat ! The maximum latitude on the input grid [degreesN]
345 real :: pole ! The sum of tracer values at the pole [a]
346 real :: max_depth ! The maximum depth of the ocean [Z ~> m]
347 real :: npole ! The number of points contributing to the pole value [nondim]
348 real :: missing_val_in ! The missing value in the input field [a]
349 real :: roundoff ! The magnitude of roundoff, usually ~2e-16 [nondim]
350 real :: add_offset, scale_factor ! File-specific conversion factors [a] or [nondim]
351 integer :: ans_date ! The vintage of the expressions and order of arithmetic to use
352 logical :: found_attr
353 logical :: add_np
354 logical :: is_ongrid
355 type(horiz_interp_type) :: Interp
356 type(axis_info), dimension(4) :: axes_info ! Axis information used for regridding
357 integer :: is, ie, js, je ! compute domain indices
358 integer :: isg, ieg, jsg, jeg ! global extent
359 integer :: isd, ied, jsd, jed ! data domain indices
360 integer :: id_clock_read
361 logical :: debug=.false.
362 real :: I_scale ! The inverse of the scale factor for diagnostic output [a A-1 ~> 1]
363 real :: dtr_iter_stop ! The tolerance for changes in tracer concentrations between smoothing
364 ! iterations that determines when to stop iterating [A ~> a]
365 real, dimension(SZI_(G),SZJ_(G)) :: lon_out ! The longitude of points on the model grid [radians]
366 real, dimension(SZI_(G),SZJ_(G)) :: lat_out ! The latitude of points on the model grid [radians]
367 real, dimension(SZI_(G),SZJ_(G)) :: tr_out ! The tracer on the model grid [A ~> a]
368 real, dimension(SZI_(G),SZJ_(G)) :: mask_out ! The mask on the model grid [nondim]
369 real, dimension(SZI_(G),SZJ_(G)) :: good ! Where the data is valid, this is 1 [nondim]
370 real, dimension(SZI_(G),SZJ_(G)) :: fill ! 1 where the data needs to be filled in [nondim]
371 real, dimension(SZI_(G),SZJ_(G)) :: tr_outf ! The tracer concentrations after Ice-9 [A ~> a]
372 real, dimension(SZI_(G),SZJ_(G)) :: tr_prev ! The tracer concentrations in the layer above [A ~> a]
373 real, dimension(SZI_(G),SZJ_(G)) :: good2 ! 1 where the data is valid after Ice-9 [nondim]
374 real, dimension(SZI_(G),SZJ_(G)) :: fill2 ! 1 for points that still need to be filled after Ice-9 [nondim]
375
376 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
377 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
378 isg = g%isg ; ieg = g%ieg ; jsg = g%jsg ; jeg = g%jeg
379
380 id_clock_read = cpu_clock_id('(Initialize tracer from Z) read', grain=clock_loop)
381
382 is_ongrid = .false.
383 if (present(ongrid)) is_ongrid = ongrid
384
385 dtr_iter_stop = 1.0e-3*scale
386 if (present(tr_iter_tol)) dtr_iter_stop = tr_iter_tol
387
388 i_scale = 1.0 / scale
389
390 pi_180 = atan(1.0)/45.
391
392 ans_date = 20181231
393 if (present(answers_2018)) then ; if (.not.answers_2018) ans_date = 20190101 ; endif
394 if (present(answer_date)) ans_date = answer_date
395
396 ! Open NetCDF file and if present, extract data and spatial coordinate information
397 ! The convention adopted here requires that the data be written in (i,j,k) ordering.
398
399 call cpu_clock_begin(id_clock_read)
400
401 ! A note by MJH copied from elsewhere suggests that this code may be using the model connectivity
402 ! (e.g., reentrant or tripolar) but should use the dataset's connectivity instead.
403
404 call get_var_axes_info(trim(filename), trim(varnam), axes_info)
405
406 if (allocated(z_in)) deallocate(z_in)
407 if (allocated(z_edges_in)) deallocate(z_edges_in)
408 if (allocated(tr_z)) deallocate(tr_z)
409 if (allocated(mask_z)) deallocate(mask_z)
410
411 call get_axis_info(axes_info(1),ax_size=id)
412 call get_axis_info(axes_info(2),ax_size=jd)
413 call get_axis_info(axes_info(3),ax_size=kd)
414
415 allocate(lon_in(id), lat_in(jd), z_in(kd), z_edges_in(kd+1))
416 allocate(tr_z(isd:ied,jsd:jed,kd), source=0.0)
417 allocate(mask_z(isd:ied,jsd:jed,kd), source=0.0)
418
419 call get_axis_info(axes_info(1),ax_data=lon_in)
420 call get_axis_info(axes_info(2),ax_data=lat_in)
421 call get_axis_info(axes_info(3),ax_data=z_in)
422
423 call cpu_clock_end(id_clock_read)
424
425 if (present(m_to_z)) then ; do k=1,kd ; z_in(k) = m_to_z * z_in(k) ; enddo ; endif
426
427 add_np = .false.
428 jdp = jd
429 if (.not. is_ongrid) then
430 max_lat = maxval(lat_in)
431 if (max_lat < 90.0) then
432 ! Extrapolate the input data to the north pole using the northern-most latitude.
433 add_np = .true.
434 jdp = jd+1
435 allocate(lat_inp(jdp))
436 lat_inp(1:jd) = lat_in(:)
437 lat_inp(jd+1) = 90.0
438 deallocate(lat_in)
439 allocate(lat_in(1:jdp))
440 lat_in(:) = lat_inp(:)
441 endif
442 endif
443 ! construct level cell boundaries as the mid-point between adjacent centers
444
445 ! Set the I/O attributes
446 call read_attribute(trim(filename), "_FillValue", missing_val_in, &
447 varname=trim(varnam), found=found_attr)
448 if (.not. found_attr) call mom_error(fatal, &
449 "error finding missing value for " // trim(varnam) // &
450 " in file " // trim(filename) // " in hinterp_extrap")
451 missing_value = scale * missing_val_in
452
453 call read_attribute(trim(filename), "scale_factor", scale_factor, &
454 varname=trim(varnam), found=found_attr)
455 if (.not. found_attr) scale_factor = 1.
456
457 call read_attribute(trim(filename), "add_offset", add_offset, &
458 varname=trim(varnam), found=found_attr)
459 if (.not. found_attr) add_offset = 0.
460
461 z_edges_in(1) = 0.0
462 do k=2,kd
463 z_edges_in(k) = 0.5*(z_in(k-1)+z_in(k))
464 enddo
465 z_edges_in(kd+1) = 2.0*z_in(kd) - z_in(kd-1)
466
467 if (is_ongrid) then
468 allocate(tr_in(is:ie,js:je), source=0.0)
469 allocate(tr_in_full(is:ie,js:je,kd), source=0.0)
470 allocate(mask_in(is:ie,js:je), source=0.0)
471 else
472 call horizontal_interp_init()
473 lon_in = lon_in*pi_180
474 lat_in = lat_in*pi_180
475 allocate(x_in(id,jdp), y_in(id,jdp))
476 call meshgrid(lon_in, lat_in, x_in, y_in)
477 lon_out(:,:) = g%geoLonT(:,:)*pi_180
478 lat_out(:,:) = g%geoLatT(:,:)*pi_180
479 allocate(tr_in(id,jd), source=0.0)
480 allocate(tr_inp(id,jdp), source=0.0)
481 allocate(mask_in(id,jdp), source=0.0)
482 endif
483
484 max_depth = maxval(g%bathyT(:,:)) + g%Z_ref
485 call max_across_pes(max_depth)
486
487 if (z_edges_in(kd+1) < max_depth) z_edges_in(kd+1) = max_depth
488 roundoff = 3.0*epsilon(missing_val_in)
489
490 ! Loop through each data level and interpolate to model grid.
491 ! After interpolating, fill in points which will be needed to define the layers.
492
493 if (is_ongrid) then
494 start(1) = is+g%HI%idg_offset ; start(2) = js+g%HI%jdg_offset ; start(3) = 1
495 count(1) = ie-is+1 ; count(2) = je-js+1 ; count(3) = kd ; start(4) = 1 ; count(4) = 1
496 call mom_read_data(trim(filename), trim(varnam), tr_in_full, start, count, g%Domain)
497 endif
498
499 do k=1,kd
500 mask_in(:,:) = 0.0
501 tr_out(:,:) = 0.0
502
503 if (is_ongrid) then
504 tr_in(is:ie,js:je) = tr_in_full(is:ie,js:je,k)
505 do j=js,je
506 do i=is,ie
507 if (abs(tr_in(i,j)-missing_val_in) > abs(roundoff*missing_val_in)) then
508 mask_in(i,j) = 1.0
509 tr_in(i,j) = (tr_in(i,j)*scale_factor+add_offset) * scale
510 else
511 tr_in(i,j) = missing_value
512 endif
513 enddo
514 enddo
515
516 tr_out(is:ie,js:je) = tr_in(is:ie,js:je)
517
518 else ! .not.is_ongrid
519
520 start(:) = 1 ; start(3) = k
521 count(:) = 1 ; count(1) = id ; count(2) = jd
522 call read_variable(trim(filename), trim(varnam), tr_in, start=start, nread=count)
523
524 if (is_root_pe()) then
525 if (add_np) then
526 pole = 0.0 ; npole = 0.0
527 do i=1,id
528 if (abs(tr_in(i,jd)-missing_val_in) > abs(roundoff*missing_val_in)) then
529 pole = pole + tr_in(i,jd)
530 npole = npole + 1.0
531 endif
532 enddo
533 if (npole > 0) then
534 pole = pole / npole
535 else
536 pole = missing_val_in
537 endif
538 tr_inp(:,1:jd) = tr_in(:,:)
539 tr_inp(:,jdp) = pole
540 else
541 tr_inp(:,:) = tr_in(:,:)
542 endif
543 endif
544
545 call broadcast(tr_inp, id*jdp, blocking=.true.)
546
547 do j=1,jdp ; do i=1,id
548 if (abs(tr_inp(i,j)-missing_val_in) > abs(roundoff*missing_val_in)) then
549 mask_in(i,j) = 1.0
550 tr_inp(i,j) = (tr_inp(i,j)*scale_factor+add_offset) * scale
551 else
552 tr_inp(i,j) = missing_value
553 endif
554 enddo ; enddo
555
556 ! call fms routine horiz_interp to interpolate input level data to model horizontal grid
557 if (k == 1) then
558 call build_horiz_interp_weights(interp, x_in, y_in, lon_out(is:ie,js:je), lat_out(is:ie,js:je), &
559 interp_method='bilinear', src_modulo=.true.)
560 endif
561
562 if (debug) then
563 call mystats(tr_inp, missing_value, g, k, 'Tracer from file', unscale=i_scale, full_halo=.true.)
564 endif
565
566 call run_horiz_interp(interp, tr_inp, tr_out(is:ie,js:je), missing_value=missing_value)
567 endif ! End of .not.is_ongrid
568
569 mask_out(:,:) = 1.0
570 do j=js,je ; do i=is,ie
571 if (abs(tr_out(i,j)-missing_value) < abs(roundoff*missing_value)) mask_out(i,j) = 0.
572 enddo ; enddo
573
574 fill(:,:) = 0.0 ; good(:,:) = 0.0
575
576 do j=js,je ; do i=is,ie
577 if (mask_out(i,j) < 1.0) then
578 tr_out(i,j) = missing_value
579 else
580 good(i,j) = 1.0
581 endif
582 if ((g%mask2dT(i,j) == 1.0) .and. (z_edges_in(k) <= g%bathyT(i,j) + g%Z_ref) .and. &
583 (mask_out(i,j) < 1.0)) &
584 fill(i,j) = 1.0
585 enddo ; enddo
586
587 call pass_var(fill, g%Domain)
588 call pass_var(good, g%Domain)
589
590 if (debug) then
591 call mystats(tr_out, missing_value, g, k, 'variable from horiz_interp()', unscale=i_scale)
592 endif
593
594 ! Horizontally homogenize data to produce perfectly "flat" initial conditions
595 if (PRESENT(homogenize)) then ; if (homogenize) then
596 call homogenize_field(tr_out, g, tmp_scale=i_scale, weights=mask_out, answer_date=answer_date)
597 endif ; endif
598
599 ! tr_out contains input z-space data on the model grid with missing values
600 ! now fill in missing values using "ICE-nine" algorithm.
601 tr_outf(:,:) = tr_out(:,:)
602 if (k==1) tr_prev(:,:) = tr_outf(:,:)
603 good2(:,:) = good(:,:)
604 fill2(:,:) = fill(:,:)
605
606 call fill_miss_2d(tr_outf, good2, fill2, tr_prev, g, dtr_iter_stop, answer_date=ans_date)
607 if (debug) then
608 call mystats(tr_outf, missing_value, g, k, 'field from fill_miss_2d()', unscale=i_scale)
609 endif
610
611 tr_z(:,:,k) = tr_outf(:,:) * g%mask2dT(:,:)
612 mask_z(:,:,k) = good2(:,:) + fill2(:,:)
613
614 tr_prev(:,:) = tr_z(:,:,k)
615
616 if (debug) then
617 call hchksum(tr_prev, 'field after fill ', g%HI, unscale=i_scale)
618 endif
619
620 enddo ! kd
621
622 if (allocated(lat_inp)) deallocate(lat_inp)
623 deallocate(tr_in)
624 if (allocated(tr_inp)) deallocate(tr_inp)
625 if (allocated(tr_in_full)) deallocate(tr_in_full)
626
628
629!> Extrapolate and interpolate using a FMS time interpolation handle
630subroutine horiz_interp_and_extrap_tracer_fms_id(field, Time, G, tr_z, mask_z, &
631 z_in, z_edges_in, missing_value, scale, &
632 homogenize, spongeOngrid, m_to_Z, &
633 answers_2018, tr_iter_tol, answer_date, &
634 axes)
635
636 type(external_field), intent(in) :: field !< Handle for the time interpolated field
637 type(time_type), intent(in) :: Time !< A FMS time type
638 type(ocean_grid_type), intent(inout) :: G !< Grid object
639 real, allocatable, dimension(:,:,:), intent(out) :: tr_z
640 !< Allocatable tracer array on the horizontal
641 !! model grid and input-file vertical levels
642 !! in arbitrary units [A ~> a]
643 real, allocatable, dimension(:,:,:), intent(out) :: mask_z
644 !< Allocatable tracer mask array on the horizontal
645 !! model grid and input-file vertical levels [nondim]
646 real, allocatable, dimension(:), intent(out) :: z_in
647 !< Cell grid values for input data [Z ~> m]
648 real, allocatable, dimension(:), intent(out) :: z_edges_in
649 !< Cell grid edge values for input data [Z ~> m]
650 real, intent(out) :: missing_value !< The missing value in the returned array, scaled
651 !! to avoid accidentally having valid values match
652 !! missing values, in the same arbitrary units as tr_z [A ~> a]
653 real, intent(in) :: scale !< Scaling factor for tracer into the internal
654 !! units of the model [A a-1 ~> 1]
655 logical, optional, intent(in) :: homogenize !< If present and true, horizontally homogenize data
656 !! to produce perfectly "flat" initial conditions
657 logical, optional, intent(in) :: spongeOngrid !< If present and true, the sponge data are on the model grid
658 real, optional, intent(in) :: m_to_Z !< A conversion factor from meters to the units
659 !! of depth [Z m-1 ~> 1]. If missing, G%bathyT must be in m.
660 logical, optional, intent(in) :: answers_2018 !< If true, use expressions that give the same
661 !! answers as the code did in late 2018. Otherwise
662 !! add parentheses for rotational symmetry.
663 real, optional, intent(in) :: tr_iter_tol !< The tolerance for changes in tracer concentrations
664 !! between smoothing iterations that determines when to
665 !! stop iterating, in the same arbitrary units as tr_z [A ~> a]
666 integer, optional, intent(in) :: answer_date !< The vintage of the expressions in the code.
667 !! Dates before 20190101 give the same answers
668 !! as the code did in late 2018, while later versions
669 !! add parentheses for rotational symmetry.
670 type(axis_info), allocatable, dimension(:), optional, intent(inout) :: axes !< Axis types for the input data
671
672 ! Local variables
673 ! In the following comments, [A] is used to indicate the arbitrary, possibly rescaled units of the
674 ! input array while [a] indicates the unscaled (e.g., mks) units that can be used with the reproducing sums
675 real, dimension(:,:), allocatable :: tr_in !< A 2-d array for holding input data on its
676 !! native horizontal grid, with units that change
677 !! as the input data is interpreted [a] then [A ~> a]
678 real, dimension(:,:), allocatable :: tr_inp !< Native horizontal grid data extended to the poles
679 !! with units that change as the input data is
680 !! interpreted [a] then [A ~> a]
681 real, dimension(:,:,:), allocatable :: data_in !< A buffer for storing the full 3-d time-interpolated array
682 !! on the original grid [a]
683 real, dimension(:,:), allocatable :: mask_in !< A 2-d mask for extended input grid [nondim]
684
685 real :: PI_180 ! A conversion factor from degrees to radians [radians degree-1]
686 integer :: id, jd, kd, jdp ! Input dataset data sizes
687 integer :: i, j, k
688 real, dimension(:,:), allocatable :: x_in ! Input file longitudes [radians]
689 real, dimension(:,:), allocatable :: y_in ! Input file latitudes [radians]
690 real, dimension(:), allocatable :: lon_in ! The longitudes in the input file [degreesE] then [radians]
691 real, dimension(:), allocatable :: lat_in ! The latitudes in the input file [degreesN] then [radians]
692 real, dimension(:), allocatable :: lat_inp ! The input file latitudes expanded to the pole [degreesN] then [radians]
693 real :: max_lat ! The maximum latitude on the input grid [degreesN]
694 real :: pole ! The sum of tracer values at the pole [a]
695 real :: max_depth ! The maximum depth of the ocean [Z ~> m]
696 real :: npole ! The number of points contributing to the pole value [nondim]
697 real :: missing_val_in ! The missing value in the input field [a]
698 real :: roundoff ! The magnitude of roundoff, usually ~2e-16 [nondim]
699 logical :: add_np
700 type(horiz_interp_type) :: Interp
701 type(axis_info), dimension(4) :: axes_data
702 integer :: is, ie, js, je ! compute domain indices
703 integer :: isg, ieg, jsg, jeg ! global extent
704 integer :: isd, ied, jsd, jed ! data domain indices
705 integer :: id_clock_read
706 integer, dimension(4) :: fld_sz
707 logical :: debug=.false.
708 logical :: is_ongrid
709 integer :: ans_date ! The vintage of the expressions and order of arithmetic to use
710 real :: I_scale ! The inverse of the scale factor for diagnostic output [a A-1 ~> 1]
711 real :: dtr_iter_stop ! The tolerance for changes in tracer concentrations between smoothing
712 ! iterations that determines when to stop iterating [A ~> a]
713 real, dimension(SZI_(G),SZJ_(G)) :: lon_out ! The longitude of points on the model grid [radians]
714 real, dimension(SZI_(G),SZJ_(G)) :: lat_out ! The latitude of points on the model grid [radians]
715 real, dimension(SZI_(G),SZJ_(G)) :: tr_out ! The tracer on the model grid [A ~> a]
716 real, dimension(SZI_(G),SZJ_(G)) :: mask_out ! The mask on the model grid [nondim]
717 real, dimension(SZI_(G),SZJ_(G)) :: good ! Where the data is valid, this is 1 [nondim]
718 real, dimension(SZI_(G),SZJ_(G)) :: fill ! 1 where the data needs to be filled in [nondim]
719 real, dimension(SZI_(G),SZJ_(G)) :: tr_outf ! The tracer concentrations after Ice-9 [A ~> a]
720 real, dimension(SZI_(G),SZJ_(G)) :: tr_prev ! The tracer concentrations in the layer above [A ~> a]
721 real, dimension(SZI_(G),SZJ_(G)) :: good2 ! 1 where the data is valid after Ice-9 [nondim]
722 real, dimension(SZI_(G),SZJ_(G)) :: fill2 ! 1 for points that still need to be filled after Ice-9 [nondim]
723 integer :: turns
724 integer :: verbosity
725
726 turns = g%HI%turns
727
728 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
729 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
730 isg = g%isg ; ieg = g%ieg ; jsg = g%jsg ; jeg = g%jeg
731
732 id_clock_read = cpu_clock_id('(Initialize tracer from Z) read', grain=clock_loop)
733
734 dtr_iter_stop = 1.0e-3*scale
735 if (present(tr_iter_tol)) dtr_iter_stop = tr_iter_tol
736
737 i_scale = 1.0 / scale
738
739 pi_180 = atan(1.0)/45.
740
741 ans_date = 20181231
742 if (present(answers_2018)) then ; if (.not.answers_2018) ans_date = 20190101 ; endif
743 if (present(answer_date)) ans_date = answer_date
744
745 ! Open NetCDF file and if present, extract data and spatial coordinate information
746 ! The convention adopted here requires that the data be written in (i,j,k) ordering.
747
748 call cpu_clock_begin(id_clock_read)
749
750 if (present(axes) .and. allocated(axes)) then
751 call get_external_field_info(field, size=fld_sz, missing=missing_val_in)
752 axes_data = axes
753 else
754 call get_external_field_info(field, size=fld_sz, axes=axes_data, missing=missing_val_in)
755 if (present(axes)) then
756 allocate(axes(4))
757 axes = axes_data
758 endif
759 endif
760 missing_value = scale*missing_val_in
761
762 verbosity = mom_get_verbosity()
763
764 id = fld_sz(1) ; jd = fld_sz(2) ; kd = fld_sz(3)
765
766 is_ongrid = .false.
767 if (PRESENT(spongeongrid)) is_ongrid = spongeongrid
768 if (.not. is_ongrid) then
769 allocate(lon_in(id), lat_in(jd))
770 call get_axis_info(axes_data(1), ax_data=lon_in)
771 call get_axis_info(axes_data(2), ax_data=lat_in)
772 endif
773
774 allocate(z_in(kd), z_edges_in(kd+1))
775
776 allocate(tr_z(isd:ied,jsd:jed,kd), source=0.0)
777 allocate(mask_z(isd:ied,jsd:jed,kd), source=0.0)
778
779 call get_axis_info(axes_data(3), ax_data=z_in)
780
781 if (present(m_to_z)) then ; do k=1,kd ; z_in(k) = m_to_z * z_in(k) ; enddo ; endif
782
783 call cpu_clock_end(id_clock_read)
784
785 if (.not. is_ongrid) then
786 max_lat = maxval(lat_in)
787 add_np = .false.
788 if (max_lat < 90.0) then
789 ! Extrapolate the input data to the north pole using the northern-most latitude.
790 add_np = .true.
791 jdp = jd+1
792 allocate(lat_inp(jdp))
793 lat_inp(1:jd) = lat_in(:)
794 lat_inp(jd+1) = 90.0
795 deallocate(lat_in)
796 allocate(lat_in(1:jdp))
797 lat_in(:) = lat_inp(:)
798 else
799 jdp = jd
800 endif
801 call horizontal_interp_init()
802 lon_in = lon_in*pi_180
803 lat_in = lat_in*pi_180
804 allocate(x_in(id,jdp), y_in(id,jdp))
805 call meshgrid(lon_in, lat_in, x_in, y_in)
806 lon_out(:,:) = g%geoLonT(:,:)*pi_180
807 lat_out(:,:) = g%geoLatT(:,:)*pi_180
808 allocate(data_in(id,jd,kd), source=0.0)
809 allocate(tr_in(id,jd), source=0.0)
810 allocate(tr_inp(id,jdp), source=0.0)
811 allocate(mask_in(id,jdp), source=0.0)
812 else
813 allocate(data_in(isd:ied,jsd:jed,kd))
814 endif
815
816 ! Construct level cell boundaries as the mid-point between adjacent centers.
817 z_edges_in(1) = 0.0
818 do k=2,kd
819 z_edges_in(k) = 0.5*(z_in(k-1)+z_in(k))
820 enddo
821 z_edges_in(kd+1) = 2.0*z_in(kd) - z_in(kd-1)
822
823 max_depth = maxval(g%bathyT(:,:)) + g%Z_ref
824 call max_across_pes(max_depth)
825
826 if (z_edges_in(kd+1) < max_depth) z_edges_in(kd+1) = max_depth
827
828 ! roundoff = 3.0*EPSILON(missing_value)
829 roundoff = 1.e-4
830
831 if (.not.is_ongrid) then
832 if (is_root_pe()) &
833 call time_interp_external(field, time, data_in, verbose=(verbosity>5), turns=turns)
834
835 ! Loop through each data level and interpolate to model grid.
836 ! After interpolating, fill in points which will be needed to define the layers.
837 do k=1,kd
838 if (is_root_pe()) then
839 tr_in(1:id,1:jd) = data_in(1:id,1:jd,k)
840 if (add_np) then
841 pole = 0.0 ; npole = 0.0
842 do i=1,id
843 if (abs(tr_in(i,jd)-missing_val_in) > abs(roundoff*missing_val_in)) then
844 pole = pole + tr_in(i,jd)
845 npole = npole + 1.0
846 endif
847 enddo
848 if (npole > 0) then
849 pole = pole / npole
850 else
851 pole = missing_val_in
852 endif
853 tr_inp(:,1:jd) = tr_in(:,:)
854 tr_inp(:,jdp) = pole
855 else
856 tr_inp(:,:) = tr_in(:,:)
857 endif
858 endif
859
860 call broadcast(tr_inp, id*jdp, blocking=.true.)
861
862 mask_in(:,:) = 0.0
863
864 do j=1,jdp ; do i=1,id
865 if (abs(tr_inp(i,j)-missing_val_in) > abs(roundoff*missing_val_in)) then
866 mask_in(i,j) = 1.0
867 tr_inp(i,j) = tr_inp(i,j) * scale
868 else
869 tr_inp(i,j) = missing_value
870 endif
871 enddo ; enddo
872
873 ! call fms routine horiz_interp to interpolate input level data to model horizontal grid
874 if (k == 1) then
875 call build_horiz_interp_weights(interp, x_in, y_in, lon_out(is:ie,js:je), lat_out(is:ie,js:je), &
876 interp_method='bilinear', src_modulo=.true.)
877 endif
878
879 if (debug) then
880 call mystats(tr_inp, missing_value, g, k, 'Tracer from file', unscale=i_scale, full_halo=.true.)
881 endif
882
883 tr_out(:,:) = 0.0
884
885 call run_horiz_interp(interp, tr_inp, tr_out(is:ie,js:je), missing_value=missing_value)
886
887 mask_out(:,:) = 1.0
888 do j=js,je ; do i=is,ie
889 if (abs(tr_out(i,j)-missing_value) < abs(roundoff*missing_value)) mask_out(i,j) = 0.
890 enddo ; enddo
891
892 fill(:,:) = 0.0 ; good(:,:) = 0.0
893
894 do j=js,je ; do i=is,ie
895 if (mask_out(i,j) < 1.0) then
896 tr_out(i,j) = missing_value
897 else
898 good(i,j) = 1.0
899 endif
900 if ((g%mask2dT(i,j) == 1.0) .and. (z_edges_in(k) <= g%bathyT(i,j) + g%Z_ref) .and. &
901 (mask_out(i,j) < 1.0)) &
902 fill(i,j) = 1.0
903 enddo ; enddo
904 call pass_var(fill, g%Domain)
905 call pass_var(good, g%Domain)
906
907 if (debug) then
908 call mystats(tr_out, missing_value, g, k, 'variable from horiz_interp()', unscale=i_scale)
909 endif
910
911 ! Horizontally homogenize data to produce perfectly "flat" initial conditions
912 if (PRESENT(homogenize)) then ; if (homogenize) then
913 call homogenize_field(tr_out, g, tmp_scale=i_scale, weights=mask_out, answer_date=answer_date)
914 endif ; endif
915
916 ! tr_out contains input z-space data on the model grid with missing values
917 ! now fill in missing values using "ICE-nine" algorithm.
918 tr_outf(:,:) = tr_out(:,:)
919 if (k==1) tr_prev(:,:) = tr_outf(:,:)
920 good2(:,:) = good(:,:)
921 fill2(:,:) = fill(:,:)
922
923 call fill_miss_2d(tr_outf, good2, fill2, tr_prev, g, dtr_iter_stop, answer_date=ans_date)
924
925! if (debug) then
926! call hchksum(tr_outf, 'field from fill_miss_2d ', G%HI, unscale=I_scale)
927! call myStats(tr_outf, missing_value, G, k, 'field from fill_miss_2d()', unscale=I_scale)
928! endif
929
930 tr_z(:,:,k) = tr_outf(:,:) * g%mask2dT(:,:)
931 mask_z(:,:,k) = good2(:,:) + fill2(:,:)
932 tr_prev(:,:) = tr_z(:,:,k)
933
934 if (debug) then
935 call hchksum(tr_prev, 'field after fill ', g%HI, unscale=i_scale)
936 endif
937
938 enddo ! kd
939 else
940 call time_interp_external(field, time, data_in, verbose=(verbosity>5), turns=turns)
941 do k=1,kd
942 do j=js,je
943 do i=is,ie
944 tr_z(i,j,k) = data_in(i,j,k) * scale
945 if (ans_date >= 20190101) mask_z(i,j,k) = 1.
946 if (abs(tr_z(i,j,k)-missing_value) < abs(roundoff*missing_value)) mask_z(i,j,k) = 0.
947 enddo
948 enddo
949 enddo
950 endif
951
953
954!> Replace all values of a 2-d field with the weighted average over the valid points.
955subroutine homogenize_field(field, G, tmp_scale, weights, answer_date, wt_unscale)
956 type(ocean_grid_type), intent(inout) :: g !< Ocean grid type
957 real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: field !< The tracer on the model grid in arbitrary units [A ~> a]
958 real, optional, intent(in) :: tmp_scale !< A temporary rescaling factor for the
959 !! variable that is reversed in the
960 !! return value [a A-1 ~> 1]
961 real, dimension(SZI_(G),SZJ_(G)), &
962 optional, intent(in) :: weights !< The weights for the tracer in arbitrary units that
963 !! typically differ from those used by field [B ~> b]
964 integer, optional, intent(in) :: answer_date !< The vintage of the expressions in the code.
965 !! Dates before 20230101 use non-reproducing sums
966 !! in their averages, while later versions use
967 !! reproducing sums for rotational symmetry and
968 !! consistency across PE layouts.
969 real, optional, intent(in) :: wt_unscale !< A factor that undoes any dimensional scaling
970 !! of the weights so that they can be used with
971 !! reproducing sums [b B-1 ~> 1]
972
973 ! Local variables
974 ! In the following comments, [A] and [B] are used to indicate the arbitrary, possibly rescaled
975 ! units of the input field and the weighting array, while [a] and [b] indicate the corresponding
976 ! unscaled (e.g., mks) units that can be used with the reproducing sums
977 real, dimension(G%isc:G%iec, G%jsc:G%jec) :: field_for_sums ! The field times the weights [A B ~> a b]
978 real, dimension(G%isc:G%iec, G%jsc:G%jec) :: weight ! A copy of weights, if it is present, or the
979 ! tracer-point grid mask if it weights is absent [B ~> b]
980 real :: var_unscale ! The reciprocal of the scaling factor for the field and weights [a b A-1 B-1 ~> 1]
981 real :: wt_sum ! The sum of the weights, in [B ~> b]
982 real :: varsum ! The weighted sum of field being averaged [A B ~> a b]
983 real :: varavg ! The average of the field [A ~> a]
984 logical :: use_repro_sums ! If true, use reproducing sums.
985 integer :: i, j, is, ie, js, je
986
987 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
988
989 varavg = 0.0 ! This value will be used if wt_sum is 0.
990
991 use_repro_sums = .false. ; if (present(answer_date)) use_repro_sums = (answer_date >= 20230101)
992
993 if (present(weights)) then
994 do j=js,je ; do i=is,ie
995 weight(i,j) = weights(i,j)
996 enddo ; enddo
997 else
998 do j=js,je ; do i=is,ie
999 weight(i,j) = g%mask2dT(i,j)
1000 enddo ; enddo
1001 endif
1002
1003 if (use_repro_sums) then
1004 var_unscale = 1.0 ; if (present(tmp_scale)) var_unscale = tmp_scale
1005 if (present(wt_unscale)) var_unscale = wt_unscale * var_unscale
1006
1007 do j=js,je ; do i=is,ie
1008 field_for_sums(i,j) = field(i,j) * weight(i,j)
1009 enddo ; enddo
1010
1011 wt_sum = reproducing_sum(weight, unscale=wt_unscale)
1012 if (abs(wt_sum) > 0.0) &
1013 varavg = reproducing_sum(field_for_sums, unscale=var_unscale) * (1.0 / wt_sum)
1014
1015 else ! Do the averages with order-dependent sums to reproduce older answers.
1016 wt_sum = 0 ; varsum = 0.
1017 do j=js,je ; do i=is,ie
1018 if (weight(i,j) > 0.0) then
1019 wt_sum = wt_sum + weight(i,j)
1020 varsum = varsum + field(i,j) * weight(i,j)
1021 endif
1022 enddo ; enddo
1023
1024 ! Note that these averages will not reproduce across PE layouts or grid rotation.
1025 call sum_across_pes(wt_sum)
1026 if (wt_sum > 0.0) then
1027 call sum_across_pes(varsum)
1028 varavg = varsum / wt_sum
1029 endif
1030
1031 endif
1032
1033 ! This seems like an unlikely case to ever be used, but it is needed to recreate previous behavior.
1034 if (present(tmp_scale)) then ; if (tmp_scale == 0.0) varavg = 0.0 ; endif
1035
1036 field(:,:) = varavg
1037
1038end subroutine homogenize_field
1039
1040
1041!> Create a 2d-mesh of grid coordinates from 1-d arrays.
1042subroutine meshgrid(x, y, x_T, y_T)
1043 real, dimension(:), intent(in) :: x !< input 1-dimensional vector [arbitrary]
1044 real, dimension(:), intent(in) :: y !< input 1-dimensional vector [arbitrary]
1045 real, dimension(size(x,1),size(y,1)), intent(inout) :: x_T !< output 2-dimensional array [arbitrary]
1046 real, dimension(size(x,1),size(y,1)), intent(inout) :: y_T !< output 2-dimensional array [arbitrary]
1047
1048 integer :: ni, nj, i, j
1049
1050 ni = size(x,1) ; nj = size(y,1)
1051
1052 do j=1,nj ; do i=1,ni
1053 x_t(i,j) = x(i)
1054 enddo ; enddo
1055
1056 do j=1,nj ; do i=1,ni
1057 y_t(i,j) = y(j)
1058 enddo ; enddo
1059
1060end subroutine meshgrid
1061