FMS  2024.03
Flexible Modeling System
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 !> @{
21 subroutine 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)
146 710 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)
211 801 format (/,2x,'i',4x,'is',5x,'ie',4x,'facis',4x,'facie', &
212  /,(i4,2i7,2f10.5))
213 802 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 
798 900 format (/,1x,10('-'),' output from horiz_interp ',10('-'))
799 901 format (' input: min=',f16.9,' max=',f16.9,' avg=',f22.15)
800 902 format (' output: min=',f16.9,' max=',f16.9,' avg=',f22.15)
801 903 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:172
mpp_npes
integer function mpp_npes()
Returns processor count for current pelist.
Definition: mpp_util.inc:421
mpp_pe
integer function mpp_pe()
Returns processor ID.
Definition: mpp_util.inc:407