7 implicit none ;
private
16 subroutine interpolate_column(nsrc, h_src, u_src, ndest, h_dest, missing_value, u_dest)
17 integer,
intent(in) :: nsrc
18 real,
dimension(nsrc),
intent(in) :: h_src
19 real,
dimension(nsrc+1),
intent(in) :: u_src
20 integer,
intent(in) :: ndest
21 real,
dimension(ndest),
intent(in) :: h_dest
22 real,
intent(in) :: missing_value
23 real,
dimension(ndest+1),
intent(inout) :: u_dest
28 real :: weight_a, weight_b
29 integer :: k_src, k_dest
30 logical :: still_vanished
33 still_vanished = .true.
41 do while (dh<=x_dest .and. k_src<nsrc)
53 weight_a = max(0., ( dh - x_dest ) / dh)
54 weight_b = min(1., x_dest / dh)
55 u_dest(k_dest) = weight_a * u1 + weight_b * u2
57 u_dest(k_dest) = 0.5 * ( u1 + u2 )
63 if (k_dest<=ndest)
then
64 x_dest = x_dest + h_dest(k_dest)
65 if (still_vanished .and. h_dest(k_dest)==0.)
then
68 u_dest(k_dest) = missing_value
70 still_vanished = .false.
77 still_vanished = .true.
78 do k_dest=ndest, 1, -1
79 if (still_vanished .and. h_dest(k_dest)==0.)
then
82 u_dest(k_dest+1) = missing_value
91 subroutine reintegrate_column(nsrc, h_src, uh_src, ndest, h_dest, missing_value, uh_dest)
92 integer,
intent(in) :: nsrc
93 real,
dimension(nsrc),
intent(in) :: h_src
94 real,
dimension(nsrc),
intent(in) :: uh_src
95 integer,
intent(in) :: ndest
96 real,
dimension(ndest),
intent(in) :: h_dest
97 real,
intent(in) :: missing_value
98 real,
dimension(ndest),
intent(inout) :: uh_dest
101 real :: h_src_rem, h_dest_rem, dh
102 real :: uh_src_rem, uh_dest_rem, duh
103 integer :: k_src, k_dest
105 logical :: src_ran_out, src_exists
107 uh_dest(:) = missing_value
113 src_ran_out = .false.
117 if (h_src_rem==0. .and. k_src<nsrc)
then
120 h_src_rem = h_src(k_src)
121 uh_src_rem = uh_src(k_src)
122 if (h_src_rem==0.) cycle
125 if (h_dest_rem==0. .and. k_dest<ndest)
then
128 h_dest_rem = h_dest(k_dest)
130 if (h_dest_rem==0.) cycle
132 if (k_src==nsrc .and. h_src_rem==0.)
then
133 if (src_ran_out)
exit
138 if (h_src_rem<h_dest_rem)
then
141 if (dh>0.) duh = uh_src_rem
144 h_dest_rem = max(0., h_dest_rem - dh)
145 elseif (h_src_rem>h_dest_rem)
then
148 duh = (dh / h_src_rem) * uh_src_rem
149 h_src_rem = max(0., h_src_rem - dh)
150 uh_src_rem = uh_src_rem - duh
159 uh_dest(k_dest) = uh_dest(k_dest) + duh
160 if (k_dest==ndest .and. (k_src==nsrc .or. h_dest_rem==0.))
exit
163 if (.not. src_exists) uh_dest(1:ndest) = missing_value
169 logical,
intent(in) :: verbose
171 real,
parameter :: mv=-9.999999999e9
176 write(0,*)
'==== MOM_diag_kernels: diag_vkernels_unit_tests =========='
177 if (v)
write(0,*)
'- - - - - - - - - - interpolation tests - - - - - - - - -'
180 3, (/1.,2.,3./), (/1.,2.,3.,4./), &
181 3, (/1.,2.,3./), (/1.,2.,3.,4./) )
185 3, (/1.,1.,1./), (/1.,2.,3.,4./), &
186 2, (/1.5,1.5/), (/1.,2.5,4./) )
190 2, (/1.5,1.5/), (/1.,4.,7./), &
191 3, (/1.,1.,1./), (/1.,3.,5.,7./) )
194 fail =
test_interp(v,mv,
'C: 3 layer (vanished middle) to 2', &
195 3, (/1.,0.,2./), (/1.,2.,2.,3./), &
196 2, (/1.,2./), (/1.,2.,3./) )
199 fail =
test_interp(v,mv,
'D: 3 layer (deep) to 3', &
200 3, (/1.,2.,3./), (/1.,2.,4.,7./), &
201 2, (/2.,2./), (/1.,3.,5./) )
204 fail =
test_interp(v,mv,
'E: 3 layer to 3 (deep)', &
205 3, (/1.,2.,4./), (/1.,2.,4.,8./), &
206 3, (/2.,3.,4./), (/1.,3.,6.,8./) )
209 fail =
test_interp(v,mv,
'F: 3 layer to 4 with vanished top/botton', &
210 3, (/1.,2.,4./), (/1.,2.,4.,8./), &
211 4, (/0.,2.,5.,0./), (/mv,1.,3.,8.,mv/) )
214 fail =
test_interp(v,mv,
'Fs: 3 layer to 4 with vanished top/botton (shallow)', &
215 3, (/1.,2.,4./), (/1.,2.,4.,8./), &
216 4, (/0.,2.,4.,0./), (/mv,1.,3.,7.,mv/) )
219 fail =
test_interp(v,mv,
'Fd: 3 layer to 4 with vanished top/botton (deep)', &
220 3, (/1.,2.,4./), (/1.,2.,4.,8./), &
221 4, (/0.,2.,6.,0./), (/mv,1.,3.,8.,mv/) )
224 if (v)
write(0,*)
'- - - - - - - - - - reintegration tests - - - - - - - - -'
227 3, (/1.,2.,3./), (/-5.,2.,1./), &
228 3, (/1.,2.,3./), (/-5.,2.,1./) )
232 3, (/2.,2.,2./), (/-5.,2.,1./), &
233 2, (/3.,3./), (/-4.,2./) )
237 3, (/2.,2.,2./), (/-5.,2.,1./), &
238 2, (/3.,4./), (/-4.,2./) )
242 3, (/2.,2.,2./), (/-5.,2.,1./), &
243 2, (/3.,2./), (/-4.,1.5/) )
247 3, (/2.,2.,2./), (/-5.,2.,1./), &
248 4, (/0.,3.,3.,0./), (/0.,-4.,2.,0./) )
251 fail =
test_reintegrate(v,mv,
'C: 3 layer to 4 with vanished top//middle/bottom', &
252 3, (/2.,2.,2./), (/-5.,2.,1./), &
253 5, (/0.,3.,0.,3.,0./), (/0.,-4.,0.,2.,0./) )
257 3, (/2.,2.,2./), (/-5.,2.,1./), &
258 3, (/0.,0.,0./), (/0.,0.,0./) )
262 3, (/0.,0.,0./), (/-5.,2.,1./), &
263 3, (/2.,2.,2./), (/mv, mv, mv/) )
267 3, (/0.,0.,0./), (/-5.,2.,1./), &
268 3, (/0.,0.,0./), (/mv, mv, mv/) )
272 3, (/0.,0.,0./), (/0.,0.,0./), &
273 3, (/0.,0.,0./), (/mv, mv, mv/) )
276 if (.not. fail)
write(*,*)
'Pass'
281 logical function test_interp(verbose, missing_value, msg, nsrc, h_src, u_src, ndest, h_dest, u_true)
282 logical,
intent(in) :: verbose
283 real,
intent(in) :: missing_value
284 character(len=*),
intent(in) :: msg
285 integer,
intent(in) :: nsrc
286 real,
dimension(nsrc),
intent(in) :: h_src
287 real,
dimension(nsrc+1),
intent(in) :: u_src
288 integer,
intent(in) :: ndest
289 real,
dimension(ndest),
intent(in) :: h_dest
290 real,
dimension(ndest+1),
intent(in) :: u_true
292 real,
dimension(ndest+1) :: u_dest
295 logical :: print_results
305 write(0,
'(2a)')
' Test: ',msg
306 write(0,
'(a3,3(a24))')
'k',
'u_result',
'u_true',
'error'
308 error = u_dest(k)-u_true(k)
310 write(0,
'(i3,3(1pe24.16))') k,u_dest(k),u_true(k),u_dest(k)-u_true(k)
312 write(0,
'(i3,3(1pe24.16),x,a)') k,u_dest(k),u_true(k),u_dest(k)-u_true(k),
'<--- WRONG!'
319 logical function test_reintegrate(verbose, missing_value, msg, nsrc, h_src, uh_src, ndest, h_dest, uh_true)
320 logical,
intent(in) :: verbose
321 real,
intent(in) :: missing_value
322 character(len=*),
intent(in) :: msg
323 integer,
intent(in) :: nsrc
324 real,
dimension(nsrc),
intent(in) :: h_src
325 real,
dimension(nsrc),
intent(in) :: uh_src
326 integer,
intent(in) :: ndest
327 real,
dimension(ndest),
intent(in) :: h_dest
328 real,
dimension(ndest),
intent(in) :: uh_true
330 real,
dimension(ndest) :: uh_dest
333 logical :: print_results
343 write(0,
'(2a)')
' Test: ',msg
344 write(0,
'(a3,3(a24))')
'k',
'uh_result',
'uh_true',
'error'
346 error = uh_dest(k)-uh_true(k)
348 write(0,
'(i3,3(1pe24.16))') k,uh_dest(k),uh_true(k),uh_dest(k)-uh_true(k)
350 write(0,
'(i3,3(1pe24.16),x,a)') k,uh_dest(k),uh_true(k),uh_dest(k)-uh_true(k),
'<--- WRONG!'