FMS 2025.01-dev
Flexible Modeling System
Loading...
Searching...
No Matches
horiz_interp_conserve.inc
1!***********************************************************************
2!* GNU Lesser General Public License
3!*
4!* This file is part of the GFDL Flexible Modeling System (FMS).
5!*
6!* FMS is free software: you can redistribute it and/or modify it under
7!* the terms of the GNU Lesser General Public License as published by
8!* the Free Software Foundation, either version 3 of the License, or (at
9!* your option) any later version.
10!*
11!* FMS is distributed in the hope that it will be useful, but WITHOUT
12!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13!* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
14!* for more details.
15!*
16!* You should have received a copy of the GNU Lesser General Public
17!* License along with FMS. If not, see <http://www.gnu.org/licenses/>.
18!***********************************************************************
19!> @addtogroup horiz_interp_conserve_mod
20!> @{
21subroutine horiz_interp_conserve_new_1dx1d_ ( Interp, lon_in, lat_in, lon_out, lat_out, verbose)
22 type(horiz_interp_type), intent(inout) :: Interp
23 real(FMS_HI_KIND_), intent(in), dimension(:) :: lon_in , lat_in
24 real(FMS_HI_KIND_), intent(in), dimension(:) :: lon_out, lat_out
25 integer, intent(in), optional :: verbose
26
27 !-----------------------------------------------------------------------
28 real(FMS_HI_KIND_), dimension(size(lat_out(:))-1,2) :: sph
29 real(FMS_HI_KIND_), dimension(size(lon_out(:))-1,2) :: theta
30 real(FMS_HI_KIND_), dimension(size(lat_in(:))) :: slat_in
31 real(FMS_HI_KIND_), dimension(size(lon_in(:))-1) :: dlon_in
32 real(FMS_HI_KIND_), dimension(size(lat_in(:))-1) :: dsph_in
33 real(FMS_HI_KIND_), dimension(size(lon_out(:))-1) :: dlon_out
34 real(FMS_HI_KIND_), dimension(size(lat_out(:))-1) :: dsph_out
35 real(FMS_HI_KIND_) :: blon, fac, hpi, tpi, eps
36 integer, parameter :: num_iters = 4
37 integer :: i, j, m, n, nlon_in, nlat_in, nlon_out, nlat_out, &
38 iverbose, m2, n2, iter
39 logical :: s2n
40 character(len=64) :: mesg
41 integer, parameter :: kindl = fms_hi_kind_ !< compiled kind size
42
43 if(.not. module_is_initialized) call mpp_error(fatal, &
44 'HORIZ_INTERP_CONSERVE_NEW_1DX1D_: horiz_interp_conserve_init is not called')
45
46 if(great_circle_algorithm) call mpp_error(fatal, &
47 'HORIZ_INTERP_CONSERVE_NEW_1DX1D_: great_circle_algorithm is not implemented, contact developer')
48 !-----------------------------------------------------------------------
49 iverbose = 0; if (present(verbose)) iverbose = verbose
50
51 pe = mpp_pe()
52 root_pe = mpp_root_pe()
53 !-----------------------------------------------------------------------
54 hpi = 0.5_kindl * real(pi, fms_hi_kind_)
55 tpi = 4.0_kindl * real(hpi, fms_hi_kind_)
56 interp%version = 1
57 nlon_in = size(lon_in(:))-1; nlat_in = size(lat_in(:))-1
58 nlon_out = size(lon_out(:))-1; nlat_out = size(lat_out(:))-1
59
60 allocate ( interp % HI_KIND_TYPE_ % facj (nlat_out,2), interp % jlat (nlat_out,2), &
61 interp % HI_KIND_TYPE_ % faci (nlon_out,2), interp % ilon (nlon_out,2), &
62 interp % HI_KIND_TYPE_ % area_src (nlon_in, nlat_in), &
63 interp % HI_KIND_TYPE_ % area_dst (nlon_out, nlat_out) )
64
65 !-----------------------------------------------------------------------
66 ! --- set-up for input grid boxes ---
67
68 do j = 1, nlat_in+1
69 slat_in(j) = sin(lat_in(j))
70 enddo
71
72 do j = 1, nlat_in
73 dsph_in(j) = abs(slat_in(j+1)-slat_in(j))
74 enddo
75
76 do i = 1,nlon_in
77 dlon_in(i) = abs(lon_in(i+1)-lon_in(i))
78 enddo
79
80 ! set south to north flag
81 s2n = .true.
82 if (lat_in(1) > lat_in(nlat_in+1)) s2n = .false.
83
84 !-----------------------------------------------------------------------
85 ! --- set-up for output grid boxes ---
86
87 do n = 1, nlat_out
88 dsph_out(n) = abs(sin(lat_out(n+1))-sin(lat_out(n)))
89 enddo
90
91 do m = 1,nlon_out
92 theta(m,1) = lon_out(m)
93 theta(m,2) = lon_out(m+1)
94 dlon_out(m) = abs(lon_out(m+1)-lon_out(m))
95 enddo
96
97 interp%nlon_src = nlon_in; interp%nlat_src = nlat_in
98 interp%nlon_dst = nlon_out; interp%nlat_dst = nlat_out
99 !***********************************************************************
100
101 !------ set up latitudinal indexing ------
102 !------ make sure output grid goes south to north ------
103
104 do n = 1, nlat_out
105 if (lat_out(n) < lat_out(n+1)) then
106 sph(n,1) = sin(lat_out(n))
107 sph(n,2) = sin(lat_out(n+1))
108 else
109 sph(n,1) = sin(lat_out(n+1))
110 sph(n,2) = sin(lat_out(n))
111 endif
112 enddo
113
114 interp%jlat = 0
115 do n2 = 1, 2 ! looping on grid box edges
116 do n = 1, nlat_out ! looping on output latitudes
117 eps = 0.0_kindl
118 do iter=1,num_iters
119 ! find indices from input latitudes
120 do j = 1, nlat_in
121 if ( (s2n .and. (slat_in(j)-sph(n,n2)) <= eps .and. &
122 (sph(n,n2)-slat_in(j+1)) <= eps) .or. &
123 (.not.s2n .and. (slat_in(j+1)-sph(n,n2)) <= eps .and. &
124 (sph(n,n2)-slat_in(j)) <= eps) ) then
125 interp%jlat(n,n2) = j
126 ! weight with sin(lat) to exactly conserve area-integral
127 fac = (sph(n,n2)-slat_in(j))/(slat_in(j+1)-slat_in(j))
128 if (s2n) then
129 if (n2 == 1) interp%HI_KIND_TYPE_%facj(n,n2) = 1.0_kindl - fac
130 if (n2 == 2) interp%HI_KIND_TYPE_%facj(n,n2) = fac
131 else
132 if (n2 == 1) interp%HI_KIND_TYPE_%facj(n,n2) = fac
133 if (n2 == 2) interp%HI_KIND_TYPE_%facj(n,n2) = 1.0_kindl - fac
134 endif
135 exit
136 endif
137 enddo
138 if ( interp%jlat(n,n2) /= 0 ) exit
139 ! did not find this output grid edge in the input grid
140 ! increase tolerance for multiple passes
141 eps = epsilon(sph)*real(10.0_kindl**iter, kindl)
142 enddo
143 ! no match
144 if ( interp%jlat(n,n2) == 0 ) then
145 write (mesg,710) n,sph(n,n2)
146710 format (': n,sph=',i3,f14.7,40x)
147 call mpp_error(fatal, 'horiz_interp_conserve_mod:no latitude index found'//trim(mesg))
148 endif
149 enddo
150 enddo
151
152 !------ set up longitudinal indexing ------
153
154 interp%ilon = 0
155 do m2 = 1, 2 ! looping on grid box edges
156 do m = 1, nlon_out ! looping on output longitudes
157 blon = theta(m,m2)
158 if ( blon < lon_in(1) ) blon = blon + tpi
159 if ( blon > lon_in(nlon_in+1) ) blon = blon - tpi
160 eps = 0.0_kindl
161 do iter=1,num_iters
162 ! find indices from input longitudes
163 do i = 1, nlon_in
164 if ( (lon_in(i)-blon) <= eps .and. &
165 (blon-lon_in(i+1)) <= eps ) then
166 interp%ilon(m,m2) = i
167 fac = (blon-lon_in(i))/(lon_in(i+1)-lon_in(i))
168 if (m2 == 1) interp%HI_KIND_TYPE_%faci(m,m2) = 1.0_kindl - fac
169 if (m2 == 2) interp%HI_KIND_TYPE_%faci(m,m2) = fac
170 exit
171 endif
172 enddo
173 if ( interp%ilon(m,m2) /= 0 ) exit
174 ! did not find this output grid edge in the input grid
175 ! increase tolerance for multiple passes
176 eps = epsilon(blon)*real(10.0_kindl**iter, kindl)
177 enddo
178 ! no match
179 if ( interp%ilon(m,m2) == 0 ) then
180 print *, 'lon_out,blon,blon_in,eps=', &
181 theta(m,m2),blon,lon_in(1),lon_in(nlon_in+1),eps
182 call mpp_error(fatal, 'horiz_interp_conserve_mod: no longitude index found')
183 endif
184 enddo
185 enddo
186
187 ! --- area of input grid boxes ---
188
189 do j = 1,nlat_in
190 do i = 1,nlon_in
191 interp%HI_KIND_TYPE_%area_src(i,j) = dlon_in(i) * dsph_in(j)
192 enddo
193 enddo
194
195 ! --- area of output grid boxes ---
196
197 do n = 1, nlat_out
198 do m = 1, nlon_out
199 interp%HI_KIND_TYPE_%area_dst(m,n) = dlon_out(m) * dsph_out(n)
200 enddo
201 enddo
202
203 !-----------------------------------------------------------------------
204 ! this output may be quite lengthy and is not recommended
205 ! when using more than one processor
206 if (iverbose > 2) then
207 write (*,801) (i,interp%ilon(i,1),interp%ilon(i,2), &
208 interp%HI_KIND_TYPE_%faci(i,1),interp%HI_KIND_TYPE_%faci(i,2),i=1,nlon_out)
209 write (*,802) (j,interp%jlat(j,1),interp%jlat(j,2), &
210 interp%HI_KIND_TYPE_%facj(j,1),interp%HI_KIND_TYPE_%facj(j,2),j=1,nlat_out)
211801 format (/,2x,'i',4x,'is',5x,'ie',4x,'facis',4x,'facie', &
212 /,(i4,2i7,2f10.5))
213802 format (/,2x,'j',4x,'js',5x,'je',4x,'facjs',4x,'facje', &
214 /,(i4,2i7,2f10.5))
215 endif
216 !-----------------------------------------------------------------------
217
218 interp% HI_KIND_TYPE_ % is_allocated = .true.
219 interp% interp_method = conserve
220
221 end subroutine horiz_interp_conserve_new_1dx1d_
222
223 !#######################################################################
224
225 subroutine horiz_interp_conserve_new_1dx2d_ ( Interp, lon_in, lat_in, lon_out, lat_out, &
226 mask_in, mask_out, verbose)
227 type(horiz_interp_type), intent(inout) :: Interp
228 real(FMS_HI_KIND_), intent(in), dimension(:) :: lon_in , lat_in
229 real(FMS_HI_KIND_), intent(in), dimension(:,:) :: lon_out, lat_out
230 real(FMS_HI_KIND_), intent(in), optional, dimension(:,:) :: mask_in
231 real(FMS_HI_KIND_), intent(inout), optional, dimension(:,:) :: mask_out
232 integer, intent(in), optional :: verbose
233
234
235 integer :: create_xgrid_1DX2D_order1, get_maxxgrid, maxxgrid
236 integer :: create_xgrid_great_circle
237 integer :: nlon_in, nlat_in, nlon_out, nlat_out, nxgrid, i, j
238 real(r8_kind), dimension(size(lon_in(:))-1, size(lat_in(:))-1) :: mask_src
239 integer, allocatable, dimension(:) :: i_src, j_src, i_dst, j_dst
240 real(r8_kind), allocatable, dimension(:) :: xgrid_area, clon, clat
241 real(r8_kind), allocatable, dimension(:,:) :: dst_area, lon_src, lat_src
242 real(r8_kind), allocatable, dimension(:) :: lat_in_flip
243 real(r8_kind), allocatable, dimension(:,:) :: mask_src_flip
244 real(r8_kind), allocatable, dimension(:) :: lon_in_r8, lat_in_r8
245 real(r8_kind), allocatable, dimension(:,:) :: lon_out_r8, lat_out_r8
246
247 integer :: nincrease, ndecrease
248 logical :: flip_lat
249 integer :: wordsz
250 integer(kind=1) :: one_byte(8)
251 integer, parameter :: kindl = fms_hi_kind_ !< compiled kind size
252
253 if(.not. module_is_initialized) call mpp_error(fatal, &
254 'HORIZ_INTERP_CONSERVE_NEW_1DX2D_: horiz_interp_conserve_init is not called')
255
256 wordsz=size(transfer(lon_in(1), one_byte))
257 if(wordsz .NE. 4 .AND. wordsz .NE. 8) call mpp_error(fatal, &
258 'HORIZ_INTERP_CONSERVE_NEW_1DX2D_: wordsz should be 4 or 8')
259
260 if( (size(lon_out,1) .NE. size(lat_out,1)) .OR. (size(lon_out,2) .NE. size(lat_out,2)) ) &
261 call mpp_error(fatal, 'horiz_interp_conserve_mod: size mismatch between lon_out and lat_out')
262 nlon_in = size(lon_in(:)) - 1; nlat_in = size(lat_in(:)) - 1
263 nlon_out = size(lon_out,1) - 1; nlat_out = size(lon_out,2) - 1
264
265 mask_src = 1.0_r8_kind
266 if(present(mask_in)) then
267 if( (size(mask_in,1) .NE. nlon_in) .OR. (size(mask_in,2) .NE. nlat_in)) call mpp_error(fatal, &
268 'horiz_interp_conserve_mod: size mismatch between mask_in and lon_in/lat_in')
269 mask_src = real(mask_in, r8_kind)
270 end if
271
272 maxxgrid = get_maxxgrid()
273 allocate(i_src(maxxgrid), j_src(maxxgrid), i_dst(maxxgrid), j_dst(maxxgrid) )
274 allocate( xgrid_area(maxxgrid), dst_area(nlon_out, nlat_out) )
275
276 !--- check if source latitude is flipped
277 nincrease = 0
278 ndecrease = 0
279 do j = 1, nlat_in
280 if( lat_in(j+1) > lat_in(j) ) then
281 nincrease = nincrease + 1
282 else if ( lat_in(j+1) < lat_in(j) ) then
283 ndecrease = ndecrease + 1
284 endif
285 enddo
286
287 if(nincrease == nlat_in) then
288 flip_lat = .false.
289 else if(ndecrease == nlat_in) then
290 flip_lat = .true.
291 else
292 call mpp_error(fatal, 'horiz_interp_conserve_mod: nlat_in should be equal to nincreaase or ndecrease')
293 endif
294
295 allocate(lon_out_r8(size(lon_out,1),size(lon_out,2)))
296 allocate(lat_out_r8(size(lat_out,1),size(lat_out,2)))
297 lon_out_r8 = real(lon_out, r8_kind)
298 lat_out_r8 = real(lat_out, r8_kind)
299
300 if( .not. great_circle_algorithm ) then
301 if(flip_lat) then
302 allocate(lat_in_flip(nlat_in+1), mask_src_flip(nlon_in,nlat_in))
303 do j = 1, nlat_in+1
304 lat_in_flip(j) = real(lat_in(nlat_in+2-j), r8_kind)
305 enddo
306 do j = 1, nlat_in
307 mask_src_flip(:,j) = mask_src(:,nlat_in+1-j)
308 enddo
309 allocate(lon_in_r8(size(lon_in)))
310 lon_in_r8 = real(lon_in, r8_kind)
311 nxgrid = create_xgrid_1dx2d_order1(nlon_in, nlat_in, nlon_out, nlat_out, lon_in_r8, lat_in_flip, &
312 lon_out_r8, lat_out_r8, mask_src_flip, i_src, j_src, i_dst, j_dst, xgrid_area)
313 deallocate(lon_in_r8, lat_in_flip, mask_src_flip)
314 else
315 allocate(lon_in_r8(size(lon_in)))
316 allocate(lat_in_r8(size(lat_in)))
317 lon_in_r8 = real(lon_in, r8_kind)
318 lat_in_r8 = real(lat_in, r8_kind)
319 nxgrid = create_xgrid_1dx2d_order1(nlon_in, nlat_in, nlon_out, nlat_out, lon_in_r8, lat_in_r8, lon_out_r8, &
320 & lat_out_r8, mask_src, i_src, j_src, i_dst, j_dst, xgrid_area)
321 deallocate(lon_in_r8,lat_in_r8)
322 endif
323 else
324 allocate(lon_src(nlon_in+1,nlat_in+1), lat_src(nlon_in+1,nlat_in+1))
325 allocate(clon(maxxgrid), clat(maxxgrid))
326 if(flip_lat) then
327 allocate(mask_src_flip(nlon_in,nlat_in))
328 do j = 1, nlat_in+1
329 do i = 1, nlon_in+1
330 lon_src(i,j) = real(lon_in(i), r8_kind)
331 lat_src(i,j) = real(lat_in(nlat_in+2-j), r8_kind)
332 enddo
333 enddo
334 do j = 1, nlat_in
335 mask_src_flip(:,j) = mask_src(:,nlat_in+1-j)
336 enddo
337 nxgrid = create_xgrid_great_circle(nlon_in, nlat_in, nlon_out, nlat_out, lon_src, lat_src, lon_out_r8, &
338 & lat_out_r8, mask_src_flip, i_src, j_src, i_dst, j_dst, xgrid_area, clon, clat)
339 deallocate(mask_src_flip)
340 else
341 do j = 1, nlat_in+1
342 do i = 1, nlon_in+1
343 lon_src(i,j) = real(lon_in(i), r8_kind)
344 lat_src(i,j) = real(lat_in(j), r8_kind)
345 enddo
346 enddo
347 nxgrid = create_xgrid_great_circle(nlon_in, nlat_in, nlon_out, nlat_out, lon_src, lat_src, lon_out_r8, &
348 & lat_out_r8, mask_src, i_src, j_src, i_dst, j_dst, xgrid_area, clon, clat)
349 endif
350 deallocate(lon_src, lat_src, clon, clat)
351 endif
352
353 deallocate(lon_out_r8, lat_out_r8)
354
355 allocate(interp%i_src(nxgrid), interp%j_src(nxgrid) )
356 allocate(interp%i_dst(nxgrid), interp%j_dst(nxgrid) )
357 allocate(interp%HI_KIND_TYPE_%area_frac_dst(nxgrid) )
358 interp%version = 2
359 interp%nxgrid = nxgrid
360 interp%i_src = i_src(1:nxgrid)+1 ! in C, the starting index is 0
361 interp%j_src = j_src(1:nxgrid)+1
362 if(flip_lat) interp%j_src = nlat_in+1-interp%j_src
363 interp%i_dst = i_dst(1:nxgrid)+1
364 interp%j_dst = j_dst(1:nxgrid)+1
365
366 ! sum over exchange grid area to get destination grid area
367 dst_area = 0.0_r8_kind
368 do i = 1, nxgrid
369 dst_area(interp%i_dst(i), interp%j_dst(i)) = dst_area(interp%i_dst(i), interp%j_dst(i)) + xgrid_area(i)
370 end do
371
372 do i = 1, nxgrid
373 interp%HI_KIND_TYPE_%area_frac_dst(i) = real(xgrid_area(i)/dst_area(interp%i_dst(i), interp%j_dst(i) ), &
374 fms_hi_kind_)
375 end do
376 interp%nlon_src = nlon_in; interp%nlat_src = nlat_in
377 interp%nlon_dst = nlon_out; interp%nlat_dst = nlat_out
378 if(present(mask_out)) then
379 if( (size(mask_out,1) .NE. nlon_out) .OR. (size(mask_out,2) .NE. nlat_out) ) call mpp_error(fatal, &
380 'horiz_interp_conserve_mod: size mismatch between mask_out and lon_out/lat_out')
381 mask_out = 0.0_kindl
382 do i = 1, nxgrid
383 mask_out(interp%i_dst(i),interp%j_dst(i)) = mask_out(interp%i_dst(i), &
384 & interp%j_dst(i)) + interp%HI_KIND_TYPE_%area_frac_dst(i)
385 end do
386 end if
387
388 deallocate(i_src, j_src, i_dst, j_dst, xgrid_area, dst_area )
389
390 interp% HI_KIND_TYPE_ % is_allocated = .true.
391 interp% interp_method = conserve
392
393 end subroutine horiz_interp_conserve_new_1dx2d_
394
395 !#######################################################################
396
397 subroutine horiz_interp_conserve_new_2dx1d_ ( Interp, lon_in, lat_in, lon_out, lat_out, &
398 mask_in, mask_out, verbose)
399 type(horiz_interp_type), intent(inout) :: Interp
400 real(FMS_HI_KIND_), intent(in), dimension(:,:) :: lon_in , lat_in
401 real(FMS_HI_KIND_), intent(in), dimension(:) :: lon_out, lat_out
402 real(FMS_HI_KIND_), intent(in), optional, dimension(:,:) :: mask_in
403 real(FMS_HI_KIND_), intent(inout), optional, dimension(:,:) :: mask_out
404 integer, intent(in), optional :: verbose
405
406 integer :: create_xgrid_2DX1D_order1, get_maxxgrid, maxxgrid
407 integer :: create_xgrid_great_circle
408 integer :: nlon_in, nlat_in, nlon_out, nlat_out, nxgrid, i, j
409 integer, allocatable, dimension(:) :: i_src, j_src, i_dst, j_dst
410 real(r8_kind), allocatable, dimension(:,:) :: dst_area
411 real(r8_kind), dimension(size(lon_in,1)-1, size(lon_in,2)-1) :: mask_src
412 real(r8_kind), allocatable, dimension(:) :: xgrid_area, clon, clat
413 real(r8_kind), allocatable, dimension(:) :: lon_out_r8, lat_out_r8
414 real(r8_kind), allocatable, dimension(:,:) :: lon_in_r8, lat_in_r8
415 real(r8_kind), allocatable, dimension(:,:) :: lon_dst, lat_dst
416 integer :: wordsz
417 integer(kind=1) :: one_byte(8)
418 integer, parameter :: kindl = fms_hi_kind_ !< compiled kind size
419
420 if(.not. module_is_initialized) call mpp_error(fatal, &
421 'HORIZ_INTERP_CONSERVE_NEW_2DX1D_: horiz_interp_conserve_init is not called')
422
423 wordsz=size(transfer(lon_in(1,1), one_byte))
424 if(wordsz .NE. 8) call mpp_error(fatal, &
425 'HORIZ_INTERP_CONSERVE_NEW_2DX1D_: currently only support 64-bit real(FMS_HI_KIND_), contact developer')
426
427 if( (size(lon_in,1) .NE. size(lat_in,1)) .OR. (size(lon_in,2) .NE. size(lat_in,2)) ) &
428 call mpp_error(fatal, 'horiz_interp_conserve_mod: size mismatch between lon_in and lat_in')
429 nlon_in = size(lon_in,1) - 1; nlat_in = size(lon_in,2) - 1
430 nlon_out = size(lon_out(:)) - 1; nlat_out = size(lat_out(:)) - 1
431
432 mask_src = 1.0_r8_kind
433 if(present(mask_in)) then
434 if( (size(mask_in,1) .NE. nlon_in) .OR. (size(mask_in,2) .NE. nlat_in)) call mpp_error(fatal, &
435 'horiz_interp_conserve_mod: size mismatch between mask_in and lon_in/lat_in')
436 mask_src = real(mask_in, r8_kind)
437 end if
438
439 maxxgrid = get_maxxgrid()
440 allocate(i_src(maxxgrid), j_src(maxxgrid), i_dst(maxxgrid), j_dst(maxxgrid) )
441 allocate( xgrid_area(maxxgrid), dst_area(nlon_out, nlat_out) )
442
443 allocate(lon_in_r8(size(lon_in,1), size(lon_in, 2)))
444 allocate(lat_in_r8(size(lat_in,1), size(lat_in, 2)))
445 allocate(lon_out_r8(size(lon_out)))
446 allocate(lat_out_r8(size(lat_out)))
447 lon_out_r8 = real(lon_out, r8_kind)
448 lat_out_r8 = real(lat_out, r8_kind)
449 lon_in_r8 = real(lon_in, r8_kind)
450 lat_in_r8 = real(lat_in, r8_kind)
451
452 if( .not. great_circle_algorithm ) then
453 nxgrid = create_xgrid_2dx1d_order1(nlon_in, nlat_in, nlon_out, nlat_out, lon_in_r8, lat_in_r8, &
454 lon_out_r8, lat_out_r8, mask_src, i_src, j_src, i_dst, j_dst, xgrid_area)
455 else
456 allocate(lon_dst(nlon_out+1, nlat_out+1), lat_dst(nlon_out+1, nlat_out+1) )
457 allocate(clon(maxxgrid), clat(maxxgrid))
458 do j = 1, nlat_out+1
459 do i = 1, nlon_out+1
460 lon_dst(i,j) = real(lon_out(i), r8_kind)
461 lat_dst(i,j) = real(lat_out(j), r8_kind)
462 enddo
463 enddo
464 nxgrid = create_xgrid_great_circle(nlon_in, nlat_in, nlon_out, nlat_out, lon_in_r8, lat_in_r8, lon_dst, &
465 & lat_dst, mask_src, i_src, j_src, i_dst, j_dst, xgrid_area, clon, clat)
466 endif
467 deallocate(lon_out_r8,lat_out_r8, lon_in_r8, lat_in_r8)
468 allocate(interp%i_src(nxgrid), interp%j_src(nxgrid) )
469 allocate(interp%i_dst(nxgrid), interp%j_dst(nxgrid) )
470 allocate(interp%HI_KIND_TYPE_%area_frac_dst(nxgrid) )
471 interp%version = 2
472 interp%nxgrid = nxgrid
473 interp%i_src = i_src(1:nxgrid)+1 ! in C, the starting index is 0
474 interp%j_src = j_src(1:nxgrid)+1
475 interp%i_dst = i_dst(1:nxgrid)+1
476 interp%j_dst = j_dst(1:nxgrid)+1
477
478 ! sum over exchange grid area to get destination grid area
479 dst_area = 0.0_r8_kind
480 do i = 1, nxgrid
481 dst_area(interp%i_dst(i), interp%j_dst(i)) = dst_area(interp%i_dst(i), interp%j_dst(i)) + xgrid_area(i)
482 end do
483
484 do i = 1, nxgrid
485 interp%HI_KIND_TYPE_%area_frac_dst(i) = real(xgrid_area(i)/dst_area(interp%i_dst(i), interp%j_dst(i) ), &
486 fms_hi_kind_)
487 end do
488 interp%nlon_src = nlon_in; interp%nlat_src = nlat_in
489 interp%nlon_dst = nlon_out; interp%nlat_dst = nlat_out
490 if(present(mask_out)) then
491 if( (size(mask_out,1) .NE. nlon_out) .OR. (size(mask_out,2) .NE. nlat_out) ) call mpp_error(fatal, &
492 'horiz_interp_conserve_mod: size mismatch between mask_out and lon_out/lat_out')
493 mask_out = 0.0_kindl
494 do i = 1, nxgrid
495 mask_out(interp%i_dst(i),interp%j_dst(i)) = mask_out(interp%i_dst(i), &
496 & interp%j_dst(i)) + interp%HI_KIND_TYPE_%area_frac_dst(i)
497 end do
498 end if
499
500 deallocate(i_src, j_src, i_dst, j_dst, xgrid_area, dst_area)
501
502 interp% HI_KIND_TYPE_ % is_allocated = .true.
503 interp% interp_method = conserve
504
505 end subroutine horiz_interp_conserve_new_2dx1d_
506
507 !#######################################################################
508
509 subroutine horiz_interp_conserve_new_2dx2d_ ( Interp, lon_in, lat_in, lon_out, lat_out, &
510 mask_in, mask_out, verbose)
511 type(horiz_interp_type), intent(inout) :: Interp
512 real(FMS_HI_KIND_), intent(in), dimension(:,:) :: lon_in , lat_in
513 real(FMS_HI_KIND_), intent(in), dimension(:,:) :: lon_out, lat_out
514 real(FMS_HI_KIND_), intent(in), optional, dimension(:,:) :: mask_in
515 real(FMS_HI_KIND_), intent(inout), optional, dimension(:,:) :: mask_out
516 integer, intent(in), optional :: verbose
517
518 integer :: create_xgrid_2DX2D_order1, get_maxxgrid, maxxgrid
519 integer :: create_xgrid_great_circle
520 integer :: nlon_in, nlat_in, nlon_out, nlat_out, nxgrid, i
521 integer, allocatable, dimension(:) :: i_src, j_src, i_dst, j_dst
522 real(r8_kind), dimension(size(lon_in,1)-1, size(lon_in,2)-1) :: mask_src
523 real(r8_kind), allocatable, dimension(:) :: xgrid_area, clon, clat
524 real(r8_kind), allocatable, dimension(:,:) :: dst_area
525 real(r8_kind), allocatable, dimension(:,:) :: lon_in_r8, lat_in_r8
526 real(r8_kind), allocatable, dimension(:,:) :: lon_out_r8, lat_out_r8
527 integer :: wordsz
528 integer(kind=1) :: one_byte(8)
529 integer, parameter :: kindl = fms_hi_kind_ !< compiled kind size
530
531 if(.not. module_is_initialized) call mpp_error(fatal, &
532 'HORIZ_INTERP_CONSERVE_NEW_2DX2D_: horiz_interp_conserve_init is not called')
533
534 wordsz=size(transfer(lon_in(1,1), one_byte))
535 if(wordsz .NE. 4 .AND. wordsz .NE. 8) call mpp_error(fatal, &
536 'HORIZ_INTERP_CONSERVE_NEW_2DX2D_: wordsz should be 4 or 8')
537
538 if( (size(lon_in,1) .NE. size(lat_in,1)) .OR. (size(lon_in,2) .NE. size(lat_in,2)) ) &
539 call mpp_error(fatal, 'horiz_interp_conserve_mod: size mismatch between lon_in and lat_in')
540 if( (size(lon_out,1) .NE. size(lat_out,1)) .OR. (size(lon_out,2) .NE. size(lat_out,2)) ) &
541 call mpp_error(fatal, 'horiz_interp_conserve_mod: size mismatch between lon_out and lat_out')
542 nlon_in = size(lon_in,1) - 1; nlat_in = size(lon_in,2) - 1
543 nlon_out = size(lon_out,1) - 1; nlat_out = size(lon_out,2) - 1
544
545 mask_src = 1.0_r8_kind
546 if(present(mask_in)) then
547 if( (size(mask_in,1) .NE. nlon_in) .OR. (size(mask_in,2) .NE. nlat_in)) call mpp_error(fatal, &
548 'horiz_interp_conserve_mod: size mismatch between mask_in and lon_in/lat_in')
549 mask_src = real(mask_in, r8_kind)
550 end if
551
552 maxxgrid = get_maxxgrid()
553 allocate(i_src(maxxgrid), j_src(maxxgrid), i_dst(maxxgrid), j_dst(maxxgrid) )
554 allocate( xgrid_area(maxxgrid), dst_area(nlon_out, nlat_out) )
555
556 allocate(lon_in_r8(size(lon_in,1),size(lon_in,2)))
557 allocate(lat_in_r8(size(lat_in,1),size(lat_in,2)))
558 allocate(lon_out_r8(size(lon_out,1),size(lon_out,2)))
559 allocate(lat_out_r8(size(lat_out,1),size(lat_out,2)))
560 lon_in_r8 = real(lon_in,r8_kind)
561 lat_in_r8 = real(lat_in, r8_kind)
562 lon_out_r8 = real(lon_out, r8_kind)
563 lat_out_r8 = real(lat_out, r8_kind)
564
565 if( .not. great_circle_algorithm ) then
566 nxgrid = create_xgrid_2dx2d_order1(nlon_in, nlat_in, nlon_out, nlat_out, lon_in_r8, lat_in_r8, lon_out_r8, &
567 & lat_out_r8, mask_src, i_src, j_src, i_dst, j_dst, xgrid_area)
568 else
569 allocate(clon(maxxgrid), clat(maxxgrid))
570 nxgrid = create_xgrid_great_circle(nlon_in, nlat_in, nlon_out, nlat_out, lon_in_r8, lat_in_r8, lon_out_r8, &
571 & lat_out_r8, mask_src, i_src, j_src, i_dst, j_dst, xgrid_area, clon, clat)
572 deallocate(clon, clat)
573 endif
574
575 deallocate(lon_in_r8, lat_in_r8, lon_out_r8, lat_out_r8)
576
577 allocate(interp%i_src(nxgrid), interp%j_src(nxgrid) )
578 allocate(interp%i_dst(nxgrid), interp%j_dst(nxgrid) )
579 allocate(interp%HI_KIND_TYPE_%area_frac_dst(nxgrid) )
580 interp%version = 2
581 interp%nxgrid = nxgrid
582 interp%i_src = i_src(1:nxgrid)+1 ! in C, the starting index is 0
583 interp%j_src = j_src(1:nxgrid)+1
584 interp%i_dst = i_dst(1:nxgrid)+1
585 interp%j_dst = j_dst(1:nxgrid)+1
586
587 ! sum over exchange grid area to get destination grid area
588 dst_area = 0.0_r8_kind
589 do i = 1, nxgrid
590 dst_area(interp%i_dst(i), interp%j_dst(i)) = dst_area(interp%i_dst(i), interp%j_dst(i)) + xgrid_area(i)
591 end do
592
593 do i = 1, nxgrid
594 interp%HI_KIND_TYPE_%area_frac_dst(i) = real(xgrid_area(i)/dst_area(interp%i_dst(i), interp%j_dst(i)), &
595 fms_hi_kind_)
596 end do
597
598 interp%nlon_src = nlon_in; interp%nlat_src = nlat_in
599 interp%nlon_dst = nlon_out; interp%nlat_dst = nlat_out
600 if(present(mask_out)) then
601 if( (size(mask_out,1) .NE. nlon_out) .OR. (size(mask_out,2) .NE. nlat_out) ) call mpp_error(fatal, &
602 'horiz_interp_conserve_mod: size mismatch between mask_out and lon_out/lat_out')
603 mask_out = 0.0_kindl
604 do i = 1, nxgrid
605 mask_out(interp%i_dst(i),interp%j_dst(i)) = mask_out(interp%i_dst(i), &
606 & interp%j_dst(i)) + interp%HI_KIND_TYPE_%area_frac_dst(i)
607 end do
608 end if
609
610 deallocate(i_src, j_src, i_dst, j_dst, xgrid_area, dst_area )
611
612 interp% HI_KIND_TYPE_ % is_allocated = .true.
613 interp% interp_method = conserve
614
615 end subroutine horiz_interp_conserve_new_2dx2d_
616
617 !########################################################################
618
619 !> @brief Subroutine for performing the horizontal interpolation between two grids.
620 !!
621 !> Subroutine for performing the horizontal interpolation between two grids.
622 !! horiz_interp_conserve_new must be called before calling this routine.
623 subroutine horiz_interp_conserve_( Interp, data_in, data_out, verbose, &
624 mask_in, mask_out)
625 !-----------------------------------------------------------------------
626 type (horiz_interp_type), intent(in) :: Interp
627 real(FMS_HI_KIND_), intent(in), dimension(:,:) :: data_in !< Input data on source grid
628 real(FMS_HI_KIND_), intent(out), dimension(:,:) :: data_out !< Output data on destination grid
629 integer, intent(in), optional :: verbose !< 0 = no output; 1 = min,max,means;
630 !! 2 = max output
631 real(FMS_HI_KIND_), intent(in), dimension(:,:), optional :: mask_in !< Input mask, must be the same size as
632 !! the input data. The real value of mask_in must be in the range (0.,1.).
633 !! Set mask_in=0.0 for data points that should not be used or have missing
634 !! data. mask_in will be applied only when horiz_interp_conserve_new_1d is
635 !! called. mask_in will be passed into horiz_interp_conserve_new_2d
636 real(FMS_HI_KIND_), intent(out), dimension(:,:), optional :: mask_out !< Output mask that specifies whether
637 !! data was computed. mask_out will be computed only when
638 !! horiz_interp_conserve_new_1d is called. mask_out will be computed in
639 !! horiz_interp_conserve_new_2d
640
641 ! --- error checking ---
642 if (size(data_in,1) /= interp%nlon_src .or. size(data_in,2) /= interp%nlat_src) &
643 call mpp_error(fatal, 'horiz_interp_conserve_mod: size of input array incorrect')
644
645 if (size(data_out,1) /= interp%nlon_dst .or. size(data_out,2) /= interp%nlat_dst) &
646 call mpp_error(fatal, 'horiz_interp_conserve_mod: size of output array incorrect')
647
648 select case ( interp%version)
649 case (1)
650 call horiz_interp_conserve_version1(interp, data_in, data_out, verbose, mask_in, mask_out)
651 case (2)
652 if(present(mask_in) .OR. present(mask_out) ) call mpp_error(fatal, 'HORIZ_INTERP_CONSERVE_:'// &
653 & ' for version 2, mask_in and mask_out must be passed in horiz_interp_new, not in horiz_interp')
654 call horiz_interp_conserve_version2(interp, data_in, data_out, verbose)
655 end select
656
657 end subroutine horiz_interp_conserve_
658
659 !##############################################################################
660 subroutine horiz_interp_conserve_version1_ ( Interp, data_in, data_out, verbose, &
661 mask_in, mask_out)
662 !-----------------------------------------------------------------------
663 type (horiz_interp_type), intent(in) :: Interp
664 real(FMS_HI_KIND_), intent(in), dimension(:,:) :: data_in
665 real(FMS_HI_KIND_), intent(out), dimension(:,:) :: data_out
666 integer, intent(in), optional :: verbose
667 real(FMS_HI_KIND_), intent(in), dimension(:,:), optional :: mask_in
668 real(FMS_HI_KIND_), intent(out), dimension(:,:), optional :: mask_out
669 !----------local variables----------------------------------------------------
670 integer :: m, n, nlon_in, nlat_in, nlon_out, nlat_out, &
671 miss_in, miss_out, is, ie, js, je, &
672 np, npass, iverbose
673 real(FMS_HI_KIND_) :: dsum, wsum, avg_in, min_in, max_in, &
674 avg_out, min_out, max_out, eps, asum, &
675 dwtsum, wtsum, arsum, fis, fie, fjs, fje
676 integer, parameter :: kindl = fms_hi_kind_ !< compiled kind size
677 !-----------------------------------------------------------------------
678 iverbose = 0; if (present(verbose)) iverbose = verbose
679
680 eps = epsilon(wtsum)
681
682 nlon_in = interp%nlon_src; nlat_in = interp%nlat_src
683 nlon_out = interp%nlon_dst; nlat_out = interp%nlat_dst
684
685 if (present(mask_in)) then
686 if ( count(mask_in < -.0001_kindl .or. mask_in > 1.0001_kindl) > 0 ) &
687 call mpp_error(fatal, 'horiz_interp_conserve_mod: input mask not between 0,1')
688 endif
689
690 !-----------------------------------------------------------------------
691 !---- loop through output grid boxes ----
692
693 data_out = 0.0_kindl
694 do n = 1, nlat_out
695 ! latitude window
696 ! setup ascending latitude indices and weights
697 if (interp%jlat(n,1) <= interp%jlat(n,2)) then
698 js = interp%jlat(n,1); je = interp%jlat(n,2)
699 fjs = interp%HI_KIND_TYPE_%facj(n,1); fje = interp%HI_KIND_TYPE_%facj(n,2)
700 else
701 js = interp%jlat(n,2); je = interp%jlat(n,1)
702 fjs = interp%HI_KIND_TYPE_%facj(n,2); fje = interp%HI_KIND_TYPE_%facj(n,1)
703 endif
704
705 do m = 1, nlon_out
706 ! longitude window
707 is = interp%ilon(m,1); ie = interp%ilon(m,2)
708 fis = interp%HI_KIND_TYPE_%faci(m,1); fie = interp%HI_KIND_TYPE_%faci(m,2)
709 npass = 1
710 dwtsum = 0.0_kindl
711 wtsum = 0.0_kindl
712 arsum = 0.0_kindl
713
714 ! wrap-around on input grid
715 ! sum using 2 passes (pass 1: end of input grid)
716 if ( ie < is ) then
717 ie = nlon_in
718 fie = 1.0_kindl
719 npass = 2
720 endif
721
722 do np = 1, npass
723 ! pass 2: beginning of input grid
724 if ( np == 2 ) then
725 is = 1
726 fis = 1.0_kindl
727 ie = interp%ilon(m,2)
728 fie = interp%HI_KIND_TYPE_%faci(m,2)
729 endif
730
731 ! summing data*weight and weight for single grid point
732 if (present(mask_in)) then
733 call data_sum( data_in(is:ie,js:je), interp%HI_KIND_TYPE_%area_src(is:ie,js:je), &
734 fis, fie, fjs,fje, dwtsum, wtsum, arsum, mask_in(is:ie,js:je) )
735 else if( allocated(interp%HI_KIND_TYPE_%mask_in) ) then
736 call data_sum( data_in(is:ie,js:je), interp%HI_KIND_TYPE_%area_src(is:ie,js:je), &
737 fis, fie, fjs,fje, dwtsum, wtsum, arsum, interp%HI_KIND_TYPE_%mask_in(is:ie,js:je) )
738 else
739 call data_sum( data_in(is:ie,js:je), interp%HI_KIND_TYPE_%area_src(is:ie,js:je), &
740 fis, fie, fjs,fje, dwtsum, wtsum, arsum )
741 endif
742 enddo
743
744 if (wtsum > eps) then
745 data_out(m,n) = dwtsum/wtsum
746 if (present(mask_out)) mask_out(m,n) = wtsum/arsum
747 else
748 data_out(m,n) = 0.0_kindl
749 if (present(mask_out)) mask_out(m,n) = 0.0_kindl
750 endif
751
752 enddo
753 enddo
754
755 !***********************************************************************
756 ! compute statistics: minimum, maximum, and mean
757 !-----------------------------------------------------------------------
758
759 if (iverbose > 0) then
760
761 ! compute statistics of input data
762
763 call stats(data_in, interp%HI_KIND_TYPE_%area_src, asum, dsum, wsum, min_in, max_in, miss_in, mask_in)
764 ! diagnostic messages
765 ! on the root_pe, we can calculate the global mean, minimum and maximum.
766 if(pe == root_pe) then
767 if (wsum > 0.0_kindl) then
768 avg_in=dsum/wsum
769 else
770 print *, 'horiz_interp stats: input area equals zero '
771 avg_in=0.0_kindl
772 endif
773 if (iverbose > 1) print '(2f16.11)', 'global sum area_in = ', asum, wsum
774 endif
775
776 ! compute statistics of output data
777 call stats(data_out, interp%HI_KIND_TYPE_%area_dst, asum, dsum, wsum, min_out, max_out, miss_out, mask_out)
778 ! diagnostic messages
779 if(pe == root_pe) then
780 if (wsum > 0.0_kindl ) then
781 avg_out=dsum/wsum
782 else
783 print *, 'horiz_interp stats: output area equals zero '
784 avg_out=0.0_kindl
785 endif
786 if (iverbose > 1) print '(2f16.11)', 'global sum area_out = ', asum, wsum
787 endif
788 !---- output statistics ----
789 ! the global mean, min and max are calculated on the root pe.
790 if(pe == root_pe) then
791 write (*,900)
792 write (*,901) min_in ,max_in ,avg_in
793 if (present(mask_in)) write (*,903) miss_in
794 write (*,902) min_out,max_out,avg_out
795 if (present(mask_out)) write (*,903) miss_out
796 endif
797
798900 format (/,1x,10('-'),' output from horiz_interp ',10('-'))
799901 format (' input: min=',f16.9,' max=',f16.9,' avg=',f22.15)
800902 format (' output: min=',f16.9,' max=',f16.9,' avg=',f22.15)
801903 format (' number of missing points = ',i6)
802
803 endif
804
805 !-----------------------------------------------------------------------
806 end subroutine horiz_interp_conserve_version1_
807
808 !#############################################################################
809 subroutine horiz_interp_conserve_version2_ ( Interp, data_in, data_out, verbose )
810 !-----------------------------------------------------------------------
811 type (horiz_interp_type), intent(in) :: Interp
812 real(FMS_HI_KIND_), intent(in), dimension(:,:) :: data_in
813 real(FMS_HI_KIND_), intent(out), dimension(:,:) :: data_out
814 integer, intent(in), optional :: verbose
815 integer :: i, i_src, j_src, i_dst, j_dst
816 integer, parameter :: kindl = fms_hi_kind_ !< compiled kind size
817
818 data_out = 0.0_kindl
819 do i = 1, interp%nxgrid
820 i_src = interp%i_src(i); j_src = interp%j_src(i)
821 i_dst = interp%i_dst(i); j_dst = interp%j_dst(i)
822 data_out(i_dst, j_dst) = data_out(i_dst, j_dst) + data_in(i_src,j_src)*interp%HI_KIND_TYPE_%area_frac_dst(i)
823 end do
824
825 end subroutine horiz_interp_conserve_version2_
826
827
828 !#######################################################################
829 !> This statistics is for conservative scheme
830 subroutine stats_ ( dat, area, asum, dsum, wsum, low, high, miss, mask )
831 real(FMS_HI_KIND_), intent(in) :: dat(:,:), area(:,:)
832 real(FMS_HI_KIND_), intent(out) :: asum, dsum, wsum, low, high
833 integer, intent(out) :: miss
834 real(FMS_HI_KIND_), intent(in), optional :: mask(:,:)
835 integer, parameter :: kindl = fms_hi_kind_ !< compiled kind size
836
837 integer :: pe, root_pe, npes, p, buffer_int(1)
838 real(FMS_HI_KIND_) :: buffer_real(5)
839
840 pe = mpp_pe()
841 root_pe = mpp_root_pe()
842 npes = mpp_npes()
843
844 ! sum data, data*area; and find min,max on each pe.
845
846 if (present(mask)) then
847 asum = sum(area(:,:))
848 dsum = sum(area(:,:)*dat(:,:)*mask(:,:))
849 wsum = sum(area(:,:)*mask(:,:))
850 miss = count(mask(:,:) <= 0.5_kindl)
851 low = minval(dat(:,:),mask=mask(:,:) > 0.5_kindl )
852 high = maxval(dat(:,:),mask=mask(:,:) > 0.5_kindl )
853 else
854 asum = sum(area(:,:))
855 dsum = sum(area(:,:)*dat(:,:))
856 wsum = sum(area(:,:))
857 miss = 0
858 low = minval(dat(:,:))
859 high = maxval(dat(:,:))
860 endif
861
862 ! other pe send local min, max, avg to the root pe and
863 ! root pe receive these information
864
865 if(pe == root_pe) then
866 do p = 1, npes - 1
867 ! Force use of "scalar", integer pointer mpp interface
868 call mpp_recv(buffer_real(1),glen=5,from_pe=root_pe+p, tag=comm_tag_1)
869 asum = asum + buffer_real(1)
870 dsum = dsum + buffer_real(2)
871 wsum = wsum + buffer_real(3)
872 low = min(low, buffer_real(4))
873 high = max(high, buffer_real(5))
874 call mpp_recv(buffer_int(1),glen=1,from_pe=root_pe+p, tag=comm_tag_2)
875 miss = miss + buffer_int(1)
876 enddo
877 else
878 buffer_real(1) = asum
879 buffer_real(2) = dsum
880 buffer_real(3) = wsum
881 buffer_real(4) = low
882 buffer_real(5) = high
883 ! Force use of "scalar", integer pointer mpp interface
884 call mpp_send(buffer_real(1),plen=5,to_pe=root_pe, tag=comm_tag_1)
885 buffer_int(1) = miss
886 call mpp_send(buffer_int(1),plen=1,to_pe=root_pe, tag=comm_tag_2)
887 endif
888
889 call mpp_sync_self()
890
891 end subroutine stats_
892
893 !#######################################################################
894
895 !> sums up the data and weights for a single output grid box
896 subroutine data_sum_( grid_data, area, facis, facie, facjs, facje, &
897 dwtsum, wtsum, arsum, mask )
898
899 !-----------------------------------------------------------------------
900 real(FMS_HI_KIND_), intent(in), dimension(:,:) :: grid_data, area
901 real(FMS_HI_KIND_), intent(in) :: facis, facie, facjs, facje
902 real(FMS_HI_KIND_), intent(inout) :: dwtsum, wtsum, arsum
903 real(FMS_HI_KIND_), intent(in), optional :: mask(:,:)
904
905 ! fac__ = fractional portion of each boundary grid box included
906 ! in the integral
907 ! dwtsum = sum(grid_data*area*mask)
908 ! wtsum = sum(area*mask)
909 ! arsum = sum(area)
910 !-----------------------------------------------------------------------
911 real(FMS_HI_KIND_), dimension(size(area,1),size(area,2)) :: wt
912 real(FMS_HI_KIND_) :: asum
913 integer :: id, jd
914 !-----------------------------------------------------------------------
915
916 id=size(area,1); jd=size(area,2)
917
918 wt=area
919 wt( 1,:)=wt( 1,:)*facis
920 wt(id,:)=wt(id,:)*facie
921 wt(:, 1)=wt(:, 1)*facjs
922 wt(:,jd)=wt(:,jd)*facje
923
924 asum = sum(wt)
925 arsum = arsum + asum
926
927 if (present(mask)) then
928 wt = wt * mask
929 dwtsum = dwtsum + sum(wt*grid_data)
930 wtsum = wtsum + sum(wt)
931 else
932 dwtsum = dwtsum + sum(wt*grid_data)
933 wtsum = wtsum + asum
934 endif
935 !-----------------------------------------------------------------------
936
937 end subroutine data_sum_
938!> @}
horiz_interp_conserve_
subroutine horiz_interp_conserve_(interp, data_in, data_out, verbose, mask_in, mask_out)
Subroutine for performing the horizontal interpolation between two grids.
Definition horiz_interp_conserve.inc:625
stats_
subroutine stats_(dat, area, asum, dsum, wsum, low, high, miss, mask)
This statistics is for conservative scheme.
Definition horiz_interp_conserve.inc:831
data_sum_
subroutine data_sum_(grid_data, area, facis, facie, facjs, facje, dwtsum, wtsum, arsum, mask)
sums up the data and weights for a single output grid box
Definition horiz_interp_conserve.inc:898
mpp_sync_self
subroutine mpp_sync_self(pelist, check, request, msg_size, msg_type)
This is to check if current PE's outstanding puts are complete but we can't use shmem_fence because w...
Definition mpp_util_mpi.inc:175