MOM6
MOM_diag_vkernels.F90
1 !> Provides kernels for single-column interpolation, re-integration (re-mapping of integrated quantities)
2 !! and intensive-variable remapping in the vertical
3 module mom_diag_vkernels
4 
5 ! This file is part of MOM6. See LICENSE.md for the license.
6 
7 use iso_fortran_env, only : stdout=>output_unit, stderr=>error_unit
8 
9 implicit none ; private
10 
11 public diag_vkernels_unit_tests
12 public interpolate_column
13 public reintegrate_column
14 
15 contains
16 
17 !> Linearly interpolate interface data, u_src, from grid h_src to a grid h_dest
18 subroutine interpolate_column(nsrc, h_src, u_src, ndest, h_dest, missing_value, u_dest)
19  integer, intent(in) :: nsrc !< Number of source cells
20  real, dimension(nsrc), intent(in) :: h_src !< Thickness of source cells
21  real, dimension(nsrc+1), intent(in) :: u_src !< Values at source cell interfaces
22  integer, intent(in) :: ndest !< Number of destination cells
23  real, dimension(ndest), intent(in) :: h_dest !< Thickness of destination cells
24  real, intent(in) :: missing_value !< Value to assign in vanished cells
25  real, dimension(ndest+1), intent(inout) :: u_dest !< Interpolated value at destination cell interfaces
26  ! Local variables
27  real :: x_dest ! Relative position of target interface
28  real :: dh ! Source cell thickness
29  real :: u1, u2 ! Values to interpolate between
30  real :: weight_a, weight_b ! Weights for interpolation
31  integer :: k_src, k_dest ! Index of cell in src and dest columns
32  logical :: still_vanished ! Used for figuring out what to mask as missing
33 
34  ! Initial values for the loop
35  still_vanished = .true.
36 
37  ! The following forces the "do while" loop to do one cycle that will set u1, u2, dh.
38  k_src = 0
39  dh = 0.
40  x_dest = 0.
41 
42  do k_dest=1, ndest+1
43  do while (dh<=x_dest .and. k_src<nsrc)
44  ! Move positions pointers forward until the interval 0 .. dh spans x_dest.
45  x_dest = x_dest - dh
46  k_src = k_src + 1
47  dh = h_src(k_src) ! Thickness of layer k_src
48 
49  ! Values that will be used for the interpolation.
50  u1 = u_src(k_src) ! Value on left of source cell
51  u2 = u_src(k_src+1) ! Value on right of source cell
52  enddo
53 
54  if (dh>0.) then
55  weight_a = max(0., ( dh - x_dest ) / dh) ! Weight of u1
56  weight_b = min(1., x_dest / dh) ! Weight of u2
57  u_dest(k_dest) = weight_a * u1 + weight_b * u2 ! Linear interpolation between u1 and u2
58  else
59  u_dest(k_dest) = 0.5 * ( u1 + u2 ) ! For a vanished layer we need to do something reasonable...
60  endif
61 
62  ! Mask vanished layers at the surface which would be under an ice-shelf.
63  ! TODO: Need to figure out what to do for an isopycnal coordinate diagnostic that could
64  ! also have vanished layers at the surface.
65  if (k_dest<=ndest) then
66  x_dest = x_dest + h_dest(k_dest) ! Position of interface k_dest+1
67  if (still_vanished .and. h_dest(k_dest)==0.) then
68  ! When the layer k_dest is vanished and all layers above are also vanished, the k_dest
69  ! interface value should be missing.
70  u_dest(k_dest) = missing_value
71  else
72  still_vanished = .false.
73  endif
74  endif
75 
76  enddo
77 
78  ! Mask vanished layers on topography
79  still_vanished = .true.
80  do k_dest=ndest, 1, -1
81  if (still_vanished .and. h_dest(k_dest)==0.) then
82  ! When the layer k_dest is vanished and all layers below are also vanished, the k_dest+1
83  ! interface value should be missing.
84  u_dest(k_dest+1) = missing_value
85  else
86  exit
87  endif
88  enddo
89 
90 end subroutine interpolate_column
91 
92 !> Conservatively calculate integrated data, uh_dest, on grid h_dest, from layer-integrated data, uh_src, on grid h_src
93 subroutine reintegrate_column(nsrc, h_src, uh_src, ndest, h_dest, missing_value, uh_dest)
94  integer, intent(in) :: nsrc !< Number of source cells
95  real, dimension(nsrc), intent(in) :: h_src !< Thickness of source cells
96  real, dimension(nsrc), intent(in) :: uh_src !< Values at source cell interfaces
97  integer, intent(in) :: ndest !< Number of destination cells
98  real, dimension(ndest), intent(in) :: h_dest !< Thickness of destination cells
99  real, intent(in) :: missing_value !< Value to assign in vanished cells
100  real, dimension(ndest), intent(inout) :: uh_dest !< Interpolated value at destination cell interfaces
101  ! Local variables
102  real :: x_dest ! Relative position of target interface
103  real :: h_src_rem, h_dest_rem, dh ! Incremental thicknesses
104  real :: uh_src_rem, uh_dest_rem, duh ! Incremental amounts of stuff
105  integer :: k_src, k_dest ! Index of cell in src and dest columns
106  integer :: iter
107  logical :: src_ran_out, src_exists
108 
109  uh_dest(:) = missing_value
110 
111  k_src = 0
112  k_dest = 0
113  h_dest_rem = 0.
114  h_src_rem = 0.
115  src_ran_out = .false.
116  src_exists = .false.
117 
118  do while(.true.)
119  if (h_src_rem==0. .and. k_src<nsrc) then
120  ! Supply is empty so move to the next source cell
121  k_src = k_src + 1
122  h_src_rem = h_src(k_src)
123  uh_src_rem = uh_src(k_src)
124  if (h_src_rem==0.) cycle
125  src_exists = .true. ! This stops us masking out the entire column
126  endif
127  if (h_dest_rem==0. .and. k_dest<ndest) then
128  ! Sink has no capacity so move to the next destination cell
129  k_dest = k_dest + 1
130  h_dest_rem = h_dest(k_dest)
131  uh_dest(k_dest) = 0.
132  if (h_dest_rem==0.) cycle
133  endif
134  if (k_src==nsrc .and. h_src_rem==0.) then
135  if (src_ran_out) exit ! This is the second time implying there is no more src
136  src_ran_out = .true.
137  cycle
138  endif
139  duh = 0.
140  if (h_src_rem<h_dest_rem) then
141  ! The source cell is fully within the destination cell
142  dh = h_src_rem
143  if (dh>0.) duh = uh_src_rem
144  h_src_rem = 0.
145  uh_src_rem = 0.
146  h_dest_rem = max(0., h_dest_rem - dh)
147  elseif (h_src_rem>h_dest_rem) then
148  ! Only part of the source cell can be used up
149  dh = h_dest_rem
150  duh = (dh / h_src_rem) * uh_src_rem
151  h_src_rem = max(0., h_src_rem - dh)
152  uh_src_rem = uh_src_rem - duh
153  h_dest_rem = 0.
154  else ! h_src_rem==h_dest_rem
155  ! The source cell exactly fits the destination cell
156  duh = uh_src_rem
157  h_src_rem = 0.
158  uh_src_rem = 0.
159  h_dest_rem = 0.
160  endif
161  uh_dest(k_dest) = uh_dest(k_dest) + duh
162  if (k_dest==ndest .and. (k_src==nsrc .or. h_dest_rem==0.)) exit
163  enddo
164 
165  if (.not. src_exists) uh_dest(1:ndest) = missing_value
166 
167 end subroutine reintegrate_column
168 
169 !> Returns true if any unit tests for module MOM_diag_vkernels fail
170 logical function diag_vkernels_unit_tests(verbose)
171  logical, intent(in) :: verbose !< If true, write results to stdout
172  ! Local variables
173  real, parameter :: mv=-9.999999999e9 ! Value to use for vanished layers
174  logical :: fail, v
175 
176  v = verbose
177 
178  write(stdout,*) '==== MOM_diag_kernels: diag_vkernels_unit_tests =========='
179  if (v) write(stdout,*) '- - - - - - - - - - interpolation tests - - - - - - - - -'
180 
181  fail = test_interp(v,mv,'Identity: 3 layer', &
182  3, (/1.,2.,3./), (/1.,2.,3.,4./), &
183  3, (/1.,2.,3./), (/1.,2.,3.,4./) )
184  diag_vkernels_unit_tests = fail
185 
186  fail = test_interp(v,mv,'A: 3 layer to 2', &
187  3, (/1.,1.,1./), (/1.,2.,3.,4./), &
188  2, (/1.5,1.5/), (/1.,2.5,4./) )
189  diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail
190 
191  fail = test_interp(v,mv,'B: 2 layer to 3', &
192  2, (/1.5,1.5/), (/1.,4.,7./), &
193  3, (/1.,1.,1./), (/1.,3.,5.,7./) )
194  diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail
195 
196  fail = test_interp(v,mv,'C: 3 layer (vanished middle) to 2', &
197  3, (/1.,0.,2./), (/1.,2.,2.,3./), &
198  2, (/1.,2./), (/1.,2.,3./) )
199  diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail
200 
201  fail = test_interp(v,mv,'D: 3 layer (deep) to 3', &
202  3, (/1.,2.,3./), (/1.,2.,4.,7./), &
203  2, (/2.,2./), (/1.,3.,5./) )
204  diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail
205 
206  fail = test_interp(v,mv,'E: 3 layer to 3 (deep)', &
207  3, (/1.,2.,4./), (/1.,2.,4.,8./), &
208  3, (/2.,3.,4./), (/1.,3.,6.,8./) )
209  diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail
210 
211  fail = test_interp(v,mv,'F: 3 layer to 4 with vanished top/botton', &
212  3, (/1.,2.,4./), (/1.,2.,4.,8./), &
213  4, (/0.,2.,5.,0./), (/mv,1.,3.,8.,mv/) )
214  diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail
215 
216  fail = test_interp(v,mv,'Fs: 3 layer to 4 with vanished top/botton (shallow)', &
217  3, (/1.,2.,4./), (/1.,2.,4.,8./), &
218  4, (/0.,2.,4.,0./), (/mv,1.,3.,7.,mv/) )
219  diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail
220 
221  fail = test_interp(v,mv,'Fd: 3 layer to 4 with vanished top/botton (deep)', &
222  3, (/1.,2.,4./), (/1.,2.,4.,8./), &
223  4, (/0.,2.,6.,0./), (/mv,1.,3.,8.,mv/) )
224  diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail
225 
226  if (v) write(stdout,*) '- - - - - - - - - - reintegration tests - - - - - - - - -'
227 
228  fail = test_reintegrate(v,mv,'Identity: 3 layer', &
229  3, (/1.,2.,3./), (/-5.,2.,1./), &
230  3, (/1.,2.,3./), (/-5.,2.,1./) )
231  diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail
232 
233  fail = test_reintegrate(v,mv,'A: 3 layer to 2', &
234  3, (/2.,2.,2./), (/-5.,2.,1./), &
235  2, (/3.,3./), (/-4.,2./) )
236  diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail
237 
238  fail = test_reintegrate(v,mv,'A: 3 layer to 2 (deep)', &
239  3, (/2.,2.,2./), (/-5.,2.,1./), &
240  2, (/3.,4./), (/-4.,2./) )
241  diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail
242 
243  fail = test_reintegrate(v,mv,'A: 3 layer to 2 (shallow)', &
244  3, (/2.,2.,2./), (/-5.,2.,1./), &
245  2, (/3.,2./), (/-4.,1.5/) )
246  diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail
247 
248  fail = test_reintegrate(v,mv,'B: 3 layer to 4 with vanished top/bottom', &
249  3, (/2.,2.,2./), (/-5.,2.,1./), &
250  4, (/0.,3.,3.,0./), (/0.,-4.,2.,0./) )
251  diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail
252 
253  fail = test_reintegrate(v,mv,'C: 3 layer to 4 with vanished top//middle/bottom', &
254  3, (/2.,2.,2./), (/-5.,2.,1./), &
255  5, (/0.,3.,0.,3.,0./), (/0.,-4.,0.,2.,0./) )
256  diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail
257 
258  fail = test_reintegrate(v,mv,'D: 3 layer to 3 (vanished)', &
259  3, (/2.,2.,2./), (/-5.,2.,1./), &
260  3, (/0.,0.,0./), (/0.,0.,0./) )
261  diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail
262 
263  fail = test_reintegrate(v,mv,'D: 3 layer (vanished) to 3', &
264  3, (/0.,0.,0./), (/-5.,2.,1./), &
265  3, (/2.,2.,2./), (/mv, mv, mv/) )
266  diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail
267 
268  fail = test_reintegrate(v,mv,'D: 3 layer (vanished) to 3 (vanished)', &
269  3, (/0.,0.,0./), (/-5.,2.,1./), &
270  3, (/0.,0.,0./), (/mv, mv, mv/) )
271  diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail
272 
273  fail = test_reintegrate(v,mv,'D: 3 layer (vanished) to 3 (vanished)', &
274  3, (/0.,0.,0./), (/0.,0.,0./), &
275  3, (/0.,0.,0./), (/mv, mv, mv/) )
276  diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail
277 
278  if (.not. fail) write(stdout,*) 'Pass'
279 
280 end function diag_vkernels_unit_tests
281 
282 !> Returns true if a test of interpolate_column() produces the wrong answer
283 logical function test_interp(verbose, missing_value, msg, nsrc, h_src, u_src, ndest, h_dest, u_true)
284  logical, intent(in) :: verbose !< If true, write results to stdout
285  real, intent(in) :: missing_value !< Value to indicate missing data
286  character(len=*), intent(in) :: msg !< Message to label test
287  integer, intent(in) :: nsrc !< Number of source cells
288  real, dimension(nsrc), intent(in) :: h_src !< Thickness of source cells
289  real, dimension(nsrc+1), intent(in) :: u_src !< Values at source cell interfaces
290  integer, intent(in) :: ndest !< Number of destination cells
291  real, dimension(ndest), intent(in) :: h_dest !< Thickness of destination cells
292  real, dimension(ndest+1), intent(in) :: u_true !< Correct value at destination cell interfaces
293  ! Local variables
294  real, dimension(ndest+1) :: u_dest ! Interpolated value at destination cell interfaces
295  integer :: k
296  real :: error
297  logical :: print_results
298 
299  ! Interpolate from src to dest
300  call interpolate_column(nsrc, h_src, u_src, ndest, h_dest, missing_value, u_dest)
301 
302  test_interp = .false.
303  do k=1,ndest+1
304  if (u_dest(k)/=u_true(k)) test_interp = .true.
305  enddo
306  if (verbose .or. test_interp) then
307  write(stdout,'(2a)') ' Test: ',msg
308  write(stdout,'(a3,3(a24))') 'k','u_result','u_true','error'
309  do k=1,ndest+1
310  error = u_dest(k)-u_true(k)
311  if (error==0.) then
312  write(stdout,'(i3,3(1pe24.16))') k,u_dest(k),u_true(k),u_dest(k)-u_true(k)
313  else
314  write(stdout,'(i3,3(1pe24.16),x,a)') k,u_dest(k),u_true(k),u_dest(k)-u_true(k),'<--- WRONG!'
315  write(stderr,'(i3,3(1pe24.16),x,a)') k,u_dest(k),u_true(k),u_dest(k)-u_true(k),'<--- WRONG!'
316  endif
317  enddo
318  endif
319 end function test_interp
320 
321 !> Returns true if a test of reintegrate_column() produces the wrong answer
322 logical function test_reintegrate(verbose, missing_value, msg, nsrc, h_src, uh_src, ndest, h_dest, uh_true)
323  logical, intent(in) :: verbose !< If true, write results to stdout
324  real, intent(in) :: missing_value !< Value to indicate missing data
325  character(len=*), intent(in) :: msg !< Message to label test
326  integer, intent(in) :: nsrc !< Number of source cells
327  real, dimension(nsrc), intent(in) :: h_src !< Thickness of source cells
328  real, dimension(nsrc), intent(in) :: uh_src !< Values of source cell stuff
329  integer, intent(in) :: ndest !< Number of destination cells
330  real, dimension(ndest), intent(in) :: h_dest !< Thickness of destination cells
331  real, dimension(ndest), intent(in) :: uh_true !< Correct value of destination cell stuff
332  ! Local variables
333  real, dimension(ndest) :: uh_dest ! Reintegrated value on destination cells
334  integer :: k
335  real :: error
336  logical :: print_results
337 
338  ! Interpolate from src to dest
339  call reintegrate_column(nsrc, h_src, uh_src, ndest, h_dest, missing_value, uh_dest)
340 
341  test_reintegrate = .false.
342  do k=1,ndest
343  if (uh_dest(k)/=uh_true(k)) test_reintegrate = .true.
344  enddo
345  if (verbose .or. test_reintegrate) then
346  write(stdout,'(2a)') ' Test: ',msg
347  write(stdout,'(a3,3(a24))') 'k','uh_result','uh_true','error'
348  do k=1,ndest
349  error = uh_dest(k)-uh_true(k)
350  if (error==0.) then
351  write(stdout,'(i3,3(1pe24.16))') k,uh_dest(k),uh_true(k),uh_dest(k)-uh_true(k)
352  else
353  write(stdout,'(i3,3(1pe24.16),x,a)') k,uh_dest(k),uh_true(k),uh_dest(k)-uh_true(k),'<--- WRONG!'
354  write(stderr,'(i3,3(1pe24.16),x,a)') k,uh_dest(k),uh_true(k),uh_dest(k)-uh_true(k),'<--- WRONG!'
355  endif
356  enddo
357  endif
358 end function test_reintegrate
359 
360 end module mom_diag_vkernels
mom_diag_vkernels
Provides kernels for single-column interpolation, re-integration (re-mapping of integrated quantities...
Definition: MOM_diag_vkernels.F90:3