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