FMS  2025.04
Flexible Modeling System
horiz_interp_bilinear.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_bilinear_mod
19 !> @{
20  subroutine horiz_interp_bilinear_new_1d_ ( Interp, lon_in, lat_in, lon_out, lat_out, &
21  verbose, src_modulo )
22 
23  !-----------------------------------------------------------------------
24  type(horiz_interp_type), intent(inout) :: Interp
25  real(FMS_HI_KIND_), intent(in), dimension(:) :: lon_in , lat_in
26  real(FMS_HI_KIND_), intent(in), dimension(:,:) :: lon_out, lat_out
27  integer, intent(in), optional :: verbose
28  logical, intent(in), optional :: src_modulo
29 
30  logical :: src_is_modulo
31  integer :: nlon_in, nlat_in, nlon_out, nlat_out, n, m
32  integer :: ie, is, je, js, ln_err, lt_err, warns, iunit
33  real(FMS_HI_KIND_) :: wtw, wte, wts, wtn, lon, lat, tpi, hpi
34  real(FMS_HI_KIND_) :: glt_min, glt_max, gln_min, gln_max, min_lon, max_lon
35  integer,parameter :: kindl = fms_hi_kind_
36  logical :: decreasing_lat !< .True. if latitude is monotically decreasing
37  logical :: decreasing_lon !< .True. if longitude is monotically decreasing
38 
39  warns = 0
40  if(present(verbose)) warns = verbose
41  src_is_modulo = .true.
42  if (present(src_modulo)) src_is_modulo = src_modulo
43 
44  decreasing_lat = .false.
45  if (lat_in(1) > lat_in(2)) decreasing_lat = .true.
46 
47  decreasing_lon = .false.
48  if (lon_in(1) > lon_in(2)) decreasing_lon = .true.
49 
50  hpi = 0.5_kindl * real(pi, fms_hi_kind_)
51  tpi = 4.0_kindl * hpi
52  glt_min = hpi
53  glt_max = -hpi
54  gln_min = tpi
55  gln_max = -tpi
56  min_lon = 0.0_kindl
57  max_lon = tpi
58  ln_err = 0
59  lt_err = 0
60  !-----------------------------------------------------------------------
61 
62  allocate ( interp % HI_KIND_TYPE_ % wti (size(lon_out,1),size(lon_out,2),2), &
63  interp % HI_KIND_TYPE_ % wtj (size(lon_out,1),size(lon_out,2),2), &
64  interp % i_lon (size(lon_out,1),size(lon_out,2),2), &
65  interp % j_lat (size(lon_out,1),size(lon_out,2),2))
66  !-----------------------------------------------------------------------
67 
68  nlon_in = size(lon_in(:)) ; nlat_in = size(lat_in(:))
69  nlon_out = size(lon_out, 1); nlat_out = size(lon_out, 2)
70  interp%nlon_src = nlon_in; interp%nlat_src = nlat_in
71  interp%nlon_dst = nlon_out; interp%nlat_dst = nlat_out
72 
73  if(src_is_modulo) then
74  if(lon_in(nlon_in) - lon_in(1) .gt. tpi + real(epsln, fms_hi_kind_)) &
75  call mpp_error(fatal,'horiz_interp_bilinear_mod: '// &
76  'The range of source grid longitude should be no larger than tpi')
77 
78  if(lon_in(1) .lt. 0.0_kindl .OR. lon_in(nlon_in) > tpi ) then
79  min_lon = lon_in(1)
80  max_lon = lon_in(nlon_in)
81  endif
82  endif
83 
84  do n = 1, nlat_out
85  do m = 1, nlon_out
86  lon = lon_out(m,n)
87  lat = lat_out(m,n)
88 
89  if(src_is_modulo) then
90  if(lon .lt. min_lon) then
91  lon = lon + tpi
92  else if(lon .gt. max_lon) then
93  lon = lon - tpi
94  endif
95  else ! when the input grid is in not cyclic, the output grid should located inside
96  ! the input grid
97  if((lon .lt. lon_in(1)) .or. (lon .gt. lon_in(nlon_in))) &
98  call mpp_error(fatal,'horiz_interp_bilinear_mod: ' //&
99  'when input grid is not modulo, output grid should locate inside input grid')
100  endif
101 
102  glt_min = min(lat,glt_min); glt_max = max(lat,glt_max)
103  gln_min = min(lon,gln_min); gln_max = max(lon,gln_max)
104 
105  is = nearest_index(lon, lon_in )
106  if (decreasing_lon) then
107  ! Lon_in is increasing
108  ! This is so that is is the lower bound.
109  ! For example, if the array is [50 40 30 20 10] and lon is 11, `is` is going to be 5 from `nearest_index`
110  ! but it needs to be 4 so that it can use the data at lon_in(4) and lon_in(5)
111  if( lon_in(is) .lt. lon ) is = max(is-1,1)
112  else
113  ! Lon_in is increasing
114  ! This is so that is is the lower bound.
115  ! For example, if the array is [10 20 30 40 50] and lon is 49, `is` is going to be 5 from `nearest_index`
116  ! but it needs to be 4 so that it can use the data at lon_in(4) and lon_in(5)
117  if( lon_in(is) .gt. lon ) is = max(is-1,1)
118  endif
119  if( lon_in(is) .eq. lon .and. is .eq. nlon_in) is = max(is - 1,1)
120  ie = min(is+1,nlon_in)
121  if(lon_in(is) .ne. lon_in(ie) .and. (decreasing_lon .or. lon_in(is) .le. lon)) then
122  wtw = ( lon_in(ie) - lon) / (lon_in(ie) - lon_in(is) )
123  else
124  ! east or west of the last data value. this could be because a
125  ! cyclic condition is needed or the dataset is too small.
126  ln_err = 1
127  ie = 1
128  is = nlon_in
129  if (lon_in(ie) .ge. lon ) then
130  wtw = (lon_in(ie) -lon)/(lon_in(ie)-lon_in(is)+tpi+real(epsln,fms_hi_kind_))
131  else
132  wtw = (lon_in(ie)-lon+tpi+real(epsln,fms_hi_kind_))/(lon_in(ie)-lon_in(is)+tpi+real(epsln,fms_hi_kind_))
133  endif
134  endif
135  wte = 1.0_kindl - wtw
136 
137  js = nearest_index(lat, lat_in )
138  if (decreasing_lat) then
139  ! Lat_in is decreasing
140  if( lat_in(js) .lt. lat ) js = max(js - 1, 1)
141  else
142  ! Lat_in is increasing
143  if( lat_in(js) .gt. lat ) js = max(js - 1, 1)
144  endif
145  if( lat_in(js) .eq. lat .and. js .eq. nlat_in) js = max(js - 1, 1)
146  je = min(js + 1, nlat_in)
147 
148  if ( lat_in(js) .ne. lat_in(je) .and. (decreasing_lat .or. lat_in(js) .le. lat)) then
149  wts = ( lat_in(je) - lat )/(lat_in(je)-lat_in(js))
150  else
151  ! north or south of the last data value. this could be because a
152  ! pole is not included in the data set or the dataset is too small.
153  ! in either case extrapolate north or south
154  lt_err = 1
155  wts = 1.0_kindl
156  endif
157 
158  wtn = 1.0_kindl - wts
159 
160  interp % i_lon (m,n,1) = is; interp % i_lon (m,n,2) = ie
161  interp % j_lat (m,n,1) = js; interp % j_lat (m,n,2) = je
162  interp % HI_KIND_TYPE_ % wti (m,n,1) = wtw
163  interp % HI_KIND_TYPE_ % wti (m,n,2) = wte
164  interp % HI_KIND_TYPE_ % wtj (m,n,1) = wts
165  interp % HI_KIND_TYPE_ % wtj (m,n,2) = wtn
166 
167  enddo
168  enddo
169 
170  iunit = stdout()
171 
172  if (ln_err .eq. 1 .and. warns > 0) then
173  write (iunit,'(/,(1x,a))') &
174  '==> Warning: the geographic data set does not extend far ', &
175  ' enough east or west - a cyclic boundary ', &
176  ' condition was applied. check if appropriate '
177  write (iunit,'(/,(1x,a,2f8.4))') &
178  ' data required between longitudes:', gln_min, gln_max, &
179  ' data set is between longitudes:', lon_in(1), lon_in(nlon_in)
180  warns = warns - 1
181  endif
182 
183  if (lt_err .eq. 1 .and. warns > 0) then
184  write (iunit,'(/,(1x,a))') &
185  '==> Warning: the geographic data set does not extend far ',&
186  ' enough north or south - extrapolation from ',&
187  ' the nearest data was applied. this may create ',&
188  ' artificial gradients near a geographic pole '
189  write (iunit,'(/,(1x,a,2f8.4))') &
190  ' data required between latitudes:', glt_min, glt_max, &
191  ' data set is between latitudes:', lat_in(1), lat_in(nlat_in)
192  endif
193  interp% HI_KIND_TYPE_ % is_allocated = .true.
194  interp% interp_method = bilinear
195 
196  return
197 
198  end subroutine horiz_interp_bilinear_new_1d_
199 
200  !#######################################################################
201 
202  !> Initialization routine.
203  !!
204  !> Allocates space and initializes a derived-type variable
205  !! that contains pre-computed interpolation indices and weights.
206  subroutine horiz_interp_bilinear_new_2d_ ( Interp, lon_in, lat_in, lon_out, lat_out, &
207  verbose, src_modulo, new_search, no_crash_when_not_found )
208 
209  !-----------------------------------------------------------------------
210  type(horiz_interp_type), intent(inout) :: Interp !< A derived type variable containing indices
211  !! and weights for subsequent interpolations. To
212  !! reinitialize for different grid-to-grid interpolation
213  !! @ref horiz_interp_del must be used first.
214  real(FMS_HI_KIND_), intent(in), dimension(:,:) :: lon_in !< Latitude (radians) for source data grid
215  real(FMS_HI_KIND_), intent(in), dimension(:,:) :: lat_in !< Longitude (radians) for source data grid
216  real(FMS_HI_KIND_), intent(in), dimension(:,:) :: lon_out !< Longitude (radians) for output data grid
217  real(FMS_HI_KIND_), intent(in), dimension(:,:) :: lat_out !< Latitude (radians) for output data grid
218  integer, intent(in), optional :: verbose !< flag for amount of print output
219  logical, intent(in), optional :: src_modulo !< indicates if the boundary condition
220  !! along zonal boundary is cyclic or not. Cyclic when true
221  logical, intent(in), optional :: new_search
222  logical, intent(in), optional :: no_crash_when_not_found
223  integer :: warns
224  logical :: src_is_modulo
225  integer :: nlon_in, nlat_in, nlon_out, nlat_out
226  integer :: m, n, is, ie, js, je, num_solution
227  real(FMS_HI_KIND_) :: lon, lat, quadra, x, y, y1, y2
228  real(FMS_HI_KIND_) :: a1, b1, c1, d1, a2, b2, c2, d2, a, b, c
229  real(FMS_HI_KIND_) :: lon1, lat1, lon2, lat2, lon3, lat3, lon4, lat4
230  real(FMS_HI_KIND_) :: tpi, lon_min, lon_max
231  real(FMS_HI_KIND_) :: epsln2
232  logical :: use_new_search, no_crash
233 
234  integer, parameter :: kindl=fms_hi_kind_
235 
236  tpi = 2.0_kindl * real(pi, fms_hi_kind_)
237 
238  warns = 0
239  if(present(verbose)) warns = verbose
240  src_is_modulo = .true.
241  if (present(src_modulo)) src_is_modulo = src_modulo
242  use_new_search = .false.
243  if (present(new_search)) use_new_search = new_search
244  no_crash = .false.
245  if(present(no_crash_when_not_found)) no_crash = no_crash_when_not_found
246 
247  ! make sure lon and lat has the same dimension
248  if(size(lon_out,1) /= size(lat_out,1) .or. size(lon_out,2) /= size(lat_out,2) ) &
249  call mpp_error(fatal,'horiz_interp_bilinear_mod: when using bilinear ' // &
250  'interplation, the output grids should be geographical grids')
251 
252  if(size(lon_in,1) /= size(lat_in,1) .or. size(lon_in,2) /= size(lat_in,2) ) &
253  call mpp_error(fatal,'horiz_interp_bilinear_mod: when using bilinear '// &
254  'interplation, the input grids should be geographical grids')
255 
256  !--- get the grid size
257  nlon_in = size(lon_in,1) ; nlat_in = size(lat_in,2)
258  nlon_out = size(lon_out,1); nlat_out = size(lon_out,2)
259  interp%nlon_src = nlon_in; interp%nlat_src = nlat_in
260  interp%nlon_dst = nlon_out; interp%nlat_dst = nlat_out
261 
262  allocate ( interp % HI_KIND_TYPE_ % wti (size(lon_out,1),size(lon_out,2),2), &
263  interp % HI_KIND_TYPE_ % wtj (size(lon_out,1),size(lon_out,2),2), &
264  interp % i_lon (size(lon_out,1),size(lon_out,2),2), &
265  interp % j_lat (size(lon_out,1),size(lon_out,2),2))
266 
267  !--- first fine the neighbor points for the destination points.
268  if(use_new_search) then
269  epsln2 = real(epsln,fms_hi_kind_)* 1.0e5_kindl
270  call find_neighbor_new_(interp, lon_in, lat_in, lon_out, lat_out, src_is_modulo, no_crash)
271  else
272  if(kindl == r8_kind) then
273  epsln2 = epsln
274  else
275  epsln2 = epsln_r4
276  endif
277  call find_neighbor_(interp, lon_in, lat_in, lon_out, lat_out, src_is_modulo)
278  endif
279 
280  !***************************************************************************
281  ! Algorithm explanation (from disscussion with Steve Garner ) *
282  ! *
283  ! lon(x,y) = a1*x + b1*y + c1*x*y + d1 (1) *
284  ! lat(x,y) = a2*x + b2*y + c2*x*y + d2 (2) *
285  ! f (x,y) = a3*x + b3*y + c3*x*y + d3 (3) *
286  ! with x and y is between 0 and 1. *
287  ! lon1 = lon(0,0) = d1, lat1 = lat(0,0) = d2 *
288  ! lon2 = lon(1,0) = a1+d1, lat2 = lat(1,0) = a2+d2 *
289  ! lon3 = lon(1,1) = a1+b1+c1+d1, lat3 = lat(1,1) = a2+b2+c2+d2 *
290  ! lon4 = lon(0,1) = b1+d1, lat4 = lat(0,1) = b2+d2 *
291  ! where (lon1,lat1),(lon2,lat2),(lon3,lat3),(lon4,lat4) represents *
292  ! the four corners starting from the left lower corner of grid box *
293  ! that encloses a destination grid ( the rotation direction is *
294  ! counterclockwise ). With these conditions, we get *
295  ! a1 = lon2-lon1, a2 = lat2-lat1 *
296  ! b1 = lon4-lon1, b2 = lat4-lat1 *
297  ! c1 = lon3-lon2-lon4+lon1, c2 = lat3-lat2-lat4+lat1 *
298  ! d1 = lon1 d2 = lat1 *
299  ! So given any point (lon,lat), from equation (1) and (2) we can *
300  ! solve (x,y). *
301  ! From equation (3) *
302  ! f1 = f(0,0) = d3, f2 = f(1,0) = a3+d3 *
303  ! f3 = f(1,1) = a3+b3+c3+d3, f4 = f(0,1) = b3+d3 *
304  ! we obtain *
305  ! a3 = f2-f1, b3 = f4-f1 *
306  ! c3 = f3-f2-f4+f1, d3 = f1 *
307  ! at point (lon,lat) ---> (x,y) *
308  ! f(x,y) = (f2-f1)x + (f4-f1)y + (f3-f2-f4+f1)xy + f1 *
309  ! = f1*(1-x)*(1-y) + f2*x*(1-y) + f3*x*y + f4*y*(1-x) *
310  ! wtw=1-x; wte=x; wts=1-y; xtn=y *
311  ! *
312  !***************************************************************************
313 
314  lon_min = minval(lon_in);
315  lon_max = maxval(lon_in);
316  !--- calculate the weight
317  do n = 1, nlat_out
318  do m = 1, nlon_out
319  lon = lon_out(m,n)
320  lat = lat_out(m,n)
321  if(lon .lt. lon_min) then
322  lon = lon + tpi
323  else if(lon .gt. lon_max) then
324  lon = lon - tpi
325  endif
326  is = interp%i_lon(m,n,1); ie = interp%i_lon(m,n,2)
327  js = interp%j_lat(m,n,1); je = interp%j_lat(m,n,2)
328  if( is == dummy) cycle
329  lon1 = lon_in(is,js); lat1 = lat_in(is,js);
330  lon2 = lon_in(ie,js); lat2 = lat_in(ie,js);
331  lon3 = lon_in(ie,je); lat3 = lat_in(ie,je);
332  lon4 = lon_in(is,je); lat4 = lat_in(is,je);
333  if(lon .lt. lon_min) then
334  lon1 = lon1 -tpi; lon4 = lon4 - tpi
335  else if(lon .gt. lon_max) then
336  lon2 = lon2 +tpi; lon3 = lon3 + tpi
337  endif
338  a1 = lon2-lon1
339  b1 = lon4-lon1
340  c1 = lon1+lon3-lon4-lon2
341  d1 = lon1
342  a2 = lat2-lat1
343  b2 = lat4-lat1
344  c2 = lat1+lat3-lat4-lat2
345  d2 = lat1
346  !--- the coefficient of the quadratic equation
347  a = b2*c1-b1*c2
348  b = a1*b2-a2*b1+c1*d2-c2*d1+c2*lon-c1*lat
349  c = a2*lon-a1*lat+a1*d2-a2*d1
350  quadra = b*b-4._kindl*a*c
351  if(abs(quadra) < real(epsln, fms_hi_kind_)) quadra = 0.0_kindl
352  if(quadra < 0.0_kindl) call mpp_error(fatal, &
353  "horiz_interp_bilinear_mod: No solution existed for this quadratic equation")
354  if ( abs(a) .lt. epsln2) then ! a = 0 is a linear equation
355  if( abs(b) .lt. real(epsln,fms_hi_kind_)) call mpp_error(fatal, &
356  "horiz_interp_bilinear_mod: no unique solution existed for this linear equation")
357  y = -c/b
358  else
359  y1 = 0.5_kindl*(-b+sqrt(quadra))/a
360  y2 = 0.5_kindl*(-b-sqrt(quadra))/a
361  if(abs(y1) < epsln2) y1 = 0.0_kindl
362  if(abs(y2) < epsln2) y2 = 0.0_kindl
363  if(abs(1.0_kindl-y1) < epsln2) y1 = 1.0_kindl
364  if(abs(1.0_kindl-y2) < epsln2) y2 = 1.0_kindl
365  num_solution = 0
366  if(y1 >= 0.0_kindl .and. y1 <= 1.0_kindl) then
367  y = y1
368  num_solution = num_solution +1
369  endif
370  if(y2 >= 0.0_kindl .and. y2 <= 1.0_kindl) then
371  y = y2
372  num_solution = num_solution + 1
373  endif
374  if(num_solution == 0) then
375  call mpp_error(fatal, "horiz_interp_bilinear_mod: No solution found")
376  else if(num_solution == 2) then
377  call mpp_error(fatal, "horiz_interp_bilinear_mod: Two solutions found")
378  endif
379  endif
380  if(abs(a1+c1*y) < real(epsln,fms_hi_kind_)) call mpp_error(fatal, &
381  "horiz_interp_bilinear_mod: the denomenator is 0")
382  if(abs(y) < epsln2) y = 0.0_kindl
383  if(abs(1.0_kindl-y) < epsln2) y = 1.0_kindl
384  x = (lon-b1*y-d1)/(a1+c1*y)
385  if(abs(x) < epsln2) x = 0.0_kindl
386  if(abs(1.0_kindl-x) < epsln2) x = 1.0_kindl
387  ! x and y should be between 0 and 1.
388  !! Added for ECDA
389  if(use_new_search) then
390  if (x < 0.0_kindl) x = 0.0_kindl ! snz
391  if (y < 0.0_kindl) y = 0.0_kindl ! snz
392  if (x > 1.0_kindl) x = 1.0_kindl
393  if (y > 1.0_kindl) y = 1.0_kindl
394  endif
395  !! if using 4 byte reals, allow a bit more tolerance (1e-4) on weight values
396  if( kindl == r8_kind ) then
397  if( x>1.0_kindl .or. x<0.0_kindl .or. y>1.0_kindl .or. y < 0.0_kindl) &
398  call mpp_error(fatal, "horiz_interp_bilinear_mod: weight should be between 0 and 1, x=" &
399  //string(x)//" y="//string(y))
400  else
401  if( x>1.0_kindl + epsln2 .or. x<0.0_kindl - epsln2 .or. &
402  y>1.0_kindl + epsln2 .or. y < 0.0_kindl - epsln2) &
403  call mpp_error(fatal, "horiz_interp_bilinear_mod: weight should be between 0 and 1, x=" &
404  //string(x)//" y="//string(y))
405  endif
406  interp % HI_KIND_TYPE_ % wti(m,n,1)=1.0_kindl-x
407  interp % HI_KIND_TYPE_ % wti(m,n,2)=x
408  interp % HI_KIND_TYPE_ % wtj(m,n,1)=1.0_kindl-y
409  interp % HI_KIND_TYPE_ % wtj(m,n,2)=y
410  enddo
411  enddo
412 
413  interp% HI_KIND_TYPE_ % is_allocated = .true.
414  interp% interp_method = bilinear
415  end subroutine
416 
417  !#######################################################################
418  !> this routine will search the source grid to fine the grid box that encloses
419  !! each destination grid.
420  subroutine find_neighbor_ ( Interp, lon_in, lat_in, lon_out, lat_out, src_modulo )
421  type(horiz_interp_type), intent(inout) :: Interp
422  real(FMS_HI_KIND_), intent(in), dimension(:,:) :: lon_in , lat_in
423  real(FMS_HI_KIND_), intent(in), dimension(:,:) :: lon_out, lat_out
424  logical, intent(in) :: src_modulo
425  integer :: nlon_in, nlat_in, nlon_out, nlat_out
426  integer :: max_step, n, m, l, i, j, ip1, jp1, step
427  integer :: is, js, jstart, jend, istart, iend, npts
428  integer, allocatable, dimension(:) :: ilon, jlat
429  real(FMS_HI_KIND_) :: lon_min, lon_max, lon, lat, tpi
430  logical :: found
431  real(FMS_HI_KIND_) :: lon1, lat1, lon2, lat2, lon3, lat3, lon4, lat4
432 
433  integer, parameter :: kindl=fms_hi_kind_
434 
435  tpi = 2.0_kindl*real(pi, fms_hi_kind_)
436  nlon_in = size(lon_in,1) ; nlat_in = size(lat_in,2)
437  nlon_out = size(lon_out,1); nlat_out = size(lon_out,2)
438 
439  lon_min = minval(lon_in);
440  lon_max = maxval(lon_in);
441 
442  max_step = max(nlon_in,nlat_in) ! can be adjusted if needed
443  allocate(ilon(5*max_step), jlat(5*max_step) )
444 
445  do n = 1, nlat_out
446  do m = 1, nlon_out
447  found = .false.
448  lon = lon_out(m,n)
449  lat = lat_out(m,n)
450 
451  if(src_modulo) then
452  if(lon .lt. lon_min) then
453  lon = lon + tpi
454  else if(lon .gt. lon_max) then
455  lon = lon - tpi
456  endif
457  else
458  if(lon .lt. lon_min .or. lon .gt. lon_max ) &
459  call mpp_error(fatal,'horiz_interp_bilinear_mod: ' //&
460  'when input grid is not modulo, output grid should locate inside input grid')
461  endif
462  !--- search for the surrounding four points locatioon.
463  if(m==1 .and. n==1) then
464  j_loop: do j = 1, nlat_in-1
465  do i = 1, nlon_in
466  ip1 = i+1
467  jp1 = j+1
468  if(i==nlon_in) then
469  if(src_modulo)then
470  ip1 = 1
471  else
472  cycle
473  endif
474  endif
475  lon1 = lon_in(i, j); lat1 = lat_in(i,j)
476  lon2 = lon_in(ip1,j); lat2 = lat_in(ip1,j)
477  lon3 = lon_in(ip1,jp1); lat3 = lat_in(ip1,jp1)
478  lon4 = lon_in(i, jp1); lat4 = lat_in(i, jp1)
479 
480  if(lon .lt. lon_min .or. lon .gt. lon_max) then
481  if(i .ne. nlon_in) then
482  cycle
483  else
484  if(lon .lt. lon_min) then
485  lon1 = lon1 -tpi; lon4 = lon4 - tpi
486  else if(lon .gt. lon_max) then
487  lon2 = lon2 +tpi; lon3 = lon3 + tpi
488  endif
489  endif
490  endif
491 
492  if(lat .ge. intersect(lon1,lat1,lon2,lat2,lon))then ! south
493  if(lon .le. intersect(lat2,lon2,lat3,lon3,lat))then ! east
494  if(lat .le. intersect(lon3,lat3,lon4,lat4,lon))then ! north
495  if(lon .ge. intersect(lat4,lon4,lat1,lon1,lat))then ! west
496  found = .true.
497  interp % i_lon (m,n,1) = i; interp % i_lon (m,n,2) = ip1
498  interp % j_lat (m,n,1) = j; interp % j_lat (m,n,2) = jp1
499  exit j_loop
500  endif
501  endif
502  endif
503  endif
504  enddo
505  enddo j_loop
506  else
507  step = 0
508  do while ( .not. found .and. step .lt. max_step )
509  !--- take the adajcent point as the starting point
510  if(m == 1) then
511  is = interp % i_lon (m,n-1,1)
512  js = interp % j_lat (m,n-1,1)
513  else
514  is = interp % i_lon (m-1,n,1)
515  js = interp % j_lat (m-1,n,1)
516  endif
517  if(step==0) then
518  npts = 1
519  ilon(1) = is
520  jlat(1) = js
521  else
522  npts = 0
523  !--- bottom boundary
524  jstart = max(js-step,1)
525  jend = min(js+step,nlat_in)
526 
527  do l = -step, step
528  i = is+l
529  if(src_modulo)then
530  if( i < 1) then
531  i = i + nlon_in
532  else if (i > nlon_in) then
533  i = i - nlon_in
534  endif
535  if( i < 1 .or. i > nlon_in) call mpp_error(fatal, &
536  'horiz_interp_bilinear_mod: max_step is too big, decrease max_step' )
537  else
538  if( i < 1 .or. i > nlon_in) cycle
539  endif
540 
541  npts = npts + 1
542  ilon(npts) = i
543  jlat(npts) = jstart
544  enddo
545 
546  !--- right and left boundary -----------------------------------------------
547  istart = is - step
548  iend = is + step
549  if(src_modulo) then
550  if( istart < 1) istart = istart + nlon_in
551  if( iend > nlon_in) iend = iend - nlon_in
552  else
553  istart = max(istart,1)
554  iend = min(iend, nlon_in)
555  endif
556  do l = -step, step
557  j = js+l
558  if( j < 1 .or. j > nlat_in .or. j==jstart .or. j==jend) cycle
559  npts = npts+1
560  ilon(npts) = istart
561  jlat(npts) = j
562  npts = npts+1
563  ilon(npts) = iend
564  jlat(npts) = j
565  end do
566 
567  !--- top boundary
568 
569  do l = -step, step
570  i = is+l
571  if(src_modulo)then
572  if( i < 1) then
573  i = i + nlon_in
574  else if (i > nlon_in) then
575  i = i - nlon_in
576  endif
577  if( i < 1 .or. i > nlon_in) call mpp_error(fatal, &
578  'horiz_interp_bilinear_mod: max_step is too big, decrease max_step' )
579  else
580  if( i < 1 .or. i > nlon_in) cycle
581  endif
582 
583  npts = npts + 1
584  ilon(npts) = i
585  jlat(npts) = jend
586  enddo
587 
588 
589  end if
590 
591  !--- find the surrouding points
592  do l = 1, npts
593  i = ilon(l)
594  j = jlat(l)
595  ip1 = i+1
596  if(ip1>nlon_in) then
597  if(src_modulo) then
598  ip1 = 1
599  else
600  cycle
601  endif
602  endif
603  jp1 = j+1
604  if(jp1>nlat_in) cycle
605  lon1 = lon_in(i, j); lat1 = lat_in(i,j)
606  lon2 = lon_in(ip1,j); lat2 = lat_in(ip1,j)
607  lon3 = lon_in(ip1,jp1); lat3 = lat_in(ip1,jp1)
608  lon4 = lon_in(i, jp1); lat4 = lat_in(i, jp1)
609 
610  if(lon .lt. lon_min .or. lon .gt. lon_max) then
611  if(i .ne. nlon_in) then
612  cycle
613  else
614  if(lon .lt. lon_min) then
615  lon1 = lon1 -tpi; lon4 = lon4 - tpi
616  else if(lon .gt. lon_max) then
617  lon2 = lon2 +tpi; lon3 = lon3 + tpi
618  endif
619  endif
620  endif
621 
622  if(lat .ge. intersect(lon1,lat1,lon2,lat2,lon))then ! south
623  if(lon .le. intersect(lat2,lon2,lat3,lon3,lat))then ! east
624  if(lat .le. intersect(lon3,lat3,lon4,lat4,lon))then !north
625  if(lon .ge. intersect(lat4,lon4,lat1,lon1,lat))then ! west
626  found = .true.
627  is=i; js=j
628  interp % i_lon (m,n,1) = i; interp % i_lon (m,n,2) = ip1
629  interp % j_lat (m,n,1) = j; interp % j_lat (m,n,2) = jp1
630  exit
631  endif
632  endif
633  endif
634  endif
635  enddo
636  step = step + 1
637  enddo
638  endif
639  if(.not.found) then
640  print *,'lon,lat=',lon*180.0_kindl/real(pi,fms_hi_kind_),lat*180.0_kindl/real(pi,fms_hi_kind_)
641  print *,'npts=',npts
642  print *,'is,ie= ',istart,iend
643  print *,'js,je= ',jstart,jend
644  print *,'lon_in(is,js)=',lon_in(istart,jstart)*180.0_kindl/real(pi,fms_hi_kind_)
645  print *,'lon_in(ie,js)=',lon_in(iend,jstart)*180.0_kindl/real(pi,fms_hi_kind_)
646  print *,'lat_in(is,js)=',lat_in(istart,jstart)*180.0_kindl/real(pi,fms_hi_kind_)
647  print *,'lat_in(ie,js)=',lat_in(iend,jstart)*180.0_kindl/real(pi,fms_hi_kind_)
648  print *,'lon_in(is,je)=',lon_in(istart,jend)*180.0_kindl/real(pi,fms_hi_kind_)
649  print *,'lon_in(ie,je)=',lon_in(iend,jend)*180.0_kindl/real(pi,fms_hi_kind_)
650  print *,'lat_in(is,je)=',lat_in(istart,jend)*180.0_kindl/real(pi,fms_hi_kind_)
651  print *,'lat_in(ie,je)=',lat_in(iend,jend)*180.0_kindl/real(pi,fms_hi_kind_)
652 
653  call mpp_error(fatal, &
654  'FIND_NEIGHBOR_: the destination point is not inside the source grid' )
655  endif
656  enddo
657  enddo
658 
659  end subroutine
660 
661  !#######################################################################
662 
663  !> The function will return true if the point x,y is inside a polygon, or
664  !! false if it is not. If the point is exactly on the edge of a polygon,
665  !! the function will return .true.
666  function inside_polygon_(polyx, polyy, x, y)
667  real(fms_hi_kind_), dimension(:), intent(in) :: polyx !< longitude coordinates of corners
668  real(fms_hi_kind_), dimension(:), intent(in) :: polyy !< latitude coordinates of corners
669  real(fms_hi_kind_), intent(in) :: x !< x coordinate of point to be tested
670  real(fms_hi_kind_), intent(in) :: y !< y coordinate of point to be tested
671  logical :: inside_polygon_
672  integer :: i, j, nedges
673  real(fms_hi_kind_) :: xx
674 
675  inside_polygon_ = .false.
676  nedges = size(polyx(:))
677  j = nedges
678  do i = 1, nedges
679  if( (polyy(i) < y .AND. polyy(j) >= y) .OR. (polyy(j) < y .AND. polyy(i) >= y) ) then
680  xx = polyx(i)+(y-polyy(i))/(polyy(j)-polyy(i))*(polyx(j)-polyx(i))
681  if( xx == x ) then
682  inside_polygon_ = .true.
683  return
684  else if( xx < x ) then
685  inside_polygon_ = .not. inside_polygon_
686  endif
687  endif
688  j = i
689  enddo
690 
691  return
692 
693  end function
694 
695  !#######################################################################
696  !> this routine will search the source grid to fine the grid box that encloses
697  !! each destination grid.
698  subroutine find_neighbor_new_( Interp, lon_in, lat_in, lon_out, lat_out, src_modulo, no_crash )
699  type(horiz_interp_type), intent(inout) :: Interp
700  real(FMS_HI_KIND_), intent(in), dimension(:,:) :: lon_in , lat_in
701  real(FMS_HI_KIND_), intent(in), dimension(:,:) :: lon_out, lat_out
702  logical, intent(in) :: src_modulo, no_crash
703  integer :: nlon_in, nlat_in, nlon_out, nlat_out
704  integer :: max_step, n, m, l, i, j, ip1, jp1, step
705  integer :: is, js, jstart, jend, istart, iend, npts
706  integer, allocatable, dimension(:) :: ilon, jlat
707  real(FMS_HI_KIND_) :: lon_min, lon_max, lon, lat, tpi
708  logical :: found
709  real(FMS_HI_KIND_) :: polyx(4), polyy(4)
710  real(FMS_HI_KIND_) :: min_lon, min_lat, max_lon, max_lat
711 
712  integer, parameter :: step_div=8, kindl = fms_hi_kind_
713 
714  tpi = 2.0_kindl * real(pi, fms_hi_kind_)
715  nlon_in = size(lon_in,1) ; nlat_in = size(lat_in,2)
716  nlon_out = size(lon_out,1); nlat_out = size(lon_out,2)
717 
718  lon_min = minval(lon_in);
719  lon_max = maxval(lon_in);
720 
721  max_step = min(nlon_in,nlat_in)/step_div ! can be adjusted if needed
722  allocate(ilon(step_div*max_step), jlat(step_div*max_step) )
723 
724  do n = 1, nlat_out
725  do m = 1, nlon_out
726  found = .false.
727  lon = lon_out(m,n)
728  lat = lat_out(m,n)
729 
730  if(src_modulo) then
731  if(lon .lt. lon_min) then
732  lon = lon + tpi
733  else if(lon .gt. lon_max) then
734  lon = lon - tpi
735  endif
736  else
737  if(lon .lt. lon_min .or. lon .gt. lon_max ) &
738  call mpp_error(fatal,'horiz_interp_bilinear_mod: ' //&
739  'when input grid is not modulo, output grid should locate inside input grid')
740  endif
741  !--- search for the surrounding four points locatioon.
742  if(m==1 .and. n==1) then
743  j_loop: do j = 1, nlat_in-1
744  do i = 1, nlon_in
745  ip1 = i+1
746  jp1 = j+1
747  if(i==nlon_in) then
748  if(src_modulo)then
749  ip1 = 1
750  else
751  cycle
752  endif
753  endif
754 
755  polyx(1) = lon_in(i, j); polyy(1) = lat_in(i,j)
756  polyx(2) = lon_in(ip1,j); polyy(2) = lat_in(ip1,j)
757  polyx(3) = lon_in(ip1,jp1); polyy(3) = lat_in(ip1,jp1)
758  polyx(4) = lon_in(i, jp1); polyy(4) = lat_in(i, jp1)
759  if(lon .lt. lon_min .or. lon .gt. lon_max) then
760  if(i .ne. nlon_in) then
761  cycle
762  else
763  if(lon .lt. lon_min) then
764  polyx(1) = polyx(1) -tpi; polyx(4) = polyx(4) - tpi
765  else if(lon .gt. lon_max) then
766  polyx(2) = polyx(2) +tpi; polyx(3) = polyx(3) + tpi
767  endif
768  endif
769  endif
770 
771  min_lon = minval(polyx)
772  max_lon = maxval(polyx)
773  min_lat = minval(polyy)
774  max_lat = maxval(polyy)
775 ! if( lon .GE. min_lon .AND. lon .LE. max_lon .AND. &
776 ! lat .GE. min_lat .AND. lat .LE. max_lat ) then
777 ! print*, 'i =', i, 'j = ', j
778 ! print '(5f15.11)', lon, polyx
779 ! print '(5f15.11)', lat, polyy
780 ! endif
781 
782  if(inside_polygon_(polyx, polyy, lon, lat)) then
783  found = .true.
784 ! print*, " found ", i, j
785  interp % i_lon (m,n,1) = i; interp % i_lon (m,n,2) = ip1
786  interp % j_lat (m,n,1) = j; interp % j_lat (m,n,2) = jp1
787  exit j_loop
788  endif
789  enddo
790  enddo j_loop
791  else
792  step = 0
793  do while ( .not. found .and. step .lt. max_step )
794  !--- take the adajcent point as the starting point
795  if(m == 1) then
796  is = interp % i_lon (m,n-1,1)
797  js = interp % j_lat (m,n-1,1)
798  else
799  is = interp % i_lon (m-1,n,1)
800  js = interp % j_lat (m-1,n,1)
801  endif
802  if(step==0) then
803  npts = 1
804  ilon(1) = is
805  jlat(1) = js
806  else
807  npts = 0
808  !--- bottom and top boundary
809  jstart = max(js-step,1)
810  jend = min(js+step,nlat_in)
811 
812  do l = -step, step
813  i = is+l
814  if(src_modulo)then
815  if( i < 1) then
816  i = i + nlon_in
817  else if (i > nlon_in) then
818  i = i - nlon_in
819  endif
820  if( i < 1 .or. i > nlon_in) call mpp_error(fatal, &
821  'horiz_interp_bilinear_mod: max_step is too big, decrease max_step' )
822  else
823  if( i < 1 .or. i > nlon_in) cycle
824  endif
825 
826  npts = npts + 1
827  ilon(npts) = i
828  jlat(npts) = jstart
829  npts = npts + 1
830  ilon(npts) = i
831  jlat(npts) = jend
832  enddo
833 
834  !--- right and left boundary -----------------------------------------------
835  istart = is - step
836  iend = is + step
837  if(src_modulo) then
838  if( istart < 1) istart = istart + nlon_in
839  if( iend > nlon_in) iend = iend - nlon_in
840  else
841  istart = max(istart,1)
842  iend = min(iend, nlon_in)
843  endif
844  do l = -step, step
845  j = js+l
846  if( j < 1 .or. j > nlat_in) cycle
847  npts = npts+1
848  ilon(npts) = istart
849  jlat(npts) = j
850  npts = npts+1
851  ilon(npts) = iend
852  jlat(npts) = j
853  end do
854  end if
855 
856  !--- find the surrouding points
857  do l = 1, npts
858  i = ilon(l)
859  j = jlat(l)
860  ip1 = i+1
861  if(ip1>nlon_in) then
862  if(src_modulo) then
863  ip1 = 1
864  else
865  cycle
866  endif
867  endif
868  jp1 = j+1
869  if(jp1>nlat_in) cycle
870  polyx(1) = lon_in(i, j); polyy(1) = lat_in(i,j)
871  polyx(2) = lon_in(ip1,j); polyy(2) = lat_in(ip1,j)
872  polyx(3) = lon_in(ip1,jp1); polyy(3) = lat_in(ip1,jp1)
873  polyx(4) = lon_in(i, jp1); polyy(4) = lat_in(i, jp1)
874  if(inside_polygon_(polyx, polyy, lon, lat)) then
875  found = .true.
876  interp % i_lon (m,n,1) = i; interp % i_lon (m,n,2) = ip1
877  interp % j_lat (m,n,1) = j; interp % j_lat (m,n,2) = jp1
878  exit
879  endif
880  enddo
881  step = step + 1
882  enddo
883  endif
884  if(.not.found) then
885  if(no_crash) then
886  interp % i_lon (m,n,1:2) = dummy
887  interp % j_lat (m,n,1:2) = dummy
888  print*,'lon,lat=',lon,lat ! snz
889  else
890  call mpp_error(fatal, &
891  'horiz_interp_bilinear_mod: the destination point is not inside the source grid' )
892  endif
893  endif
894  enddo
895  enddo
896 
897  end subroutine
898 
899  !#######################################################################
900  function intersect_(x1, y1, x2, y2, x)
901  real(FMS_HI_KIND_), intent(in) :: x1, y1, x2, y2, x
902  real(FMS_HI_KIND_) :: INTERSECT_
903 
904  intersect_ = (y2-y1)*(x-x1)/(x2-x1) + y1
905 
906  return
907 
908  end function intersect_
909 
910  !#######################################################################
911 
912  !> Subroutine for performing the horizontal interpolation between two grids
913  !!
914  !! @ref horiz_interp_bilinear_new must be called before calling this routine.
915  subroutine horiz_interp_bilinear_ ( Interp, data_in, data_out, verbose, mask_in,mask_out, &
916  missing_value, missing_permit, new_handle_missing )
917  !-----------------------------------------------------------------------
918  type (horiz_interp_type), intent(in) :: Interp !< Derived type variable containing
919  !! interpolation indices and weights. Returned by a
920  !! previous call to horiz_interp_bilinear_new
921  real(FMS_HI_KIND_), intent(in), dimension(:,:) :: data_in !< input data on source grid
922  real(FMS_HI_KIND_), intent(out), dimension(:,:) :: data_out !< output data on source grid
923  integer, intent(in), optional :: verbose !< 0 = no output; 1 = min,max,means; 2 =
924  !! all output
925  real(FMS_HI_KIND_), intent(in), dimension(:,:), optional :: mask_in !< Input mask, must be the same size as
926  !! the input data. The real(FMS_HI_KIND_) value of mask_in must be in the
927  !! range (0.,1.). Set mask_in=0.0 for data points
928  !! that should not be used or have missing data
929  real(FMS_HI_KIND_), intent(out), dimension(:,:), optional :: mask_out !< output mask that specifies whether
930  !! data was computed
931  real(FMS_HI_KIND_), intent(in), optional :: missing_value
932  integer, intent(in), optional :: missing_permit
933  logical, intent(in), optional :: new_handle_missing
934  !-----------------------------------------------------------------------
935  integer :: nlon_in, nlat_in, nlon_out, nlat_out, n, m, &
936  is, ie, js, je, iverbose, max_missing, num_missing, &
937  miss_in, miss_out, iunit
938  real(FMS_HI_KIND_) :: dwtsum, wtsum, min_in, max_in, avg_in, &
939  min_out, max_out, avg_out, wtw, wte, wts, wtn
940  real(FMS_HI_KIND_) :: mask(size(data_in,1), size(data_in,2) )
941  logical :: set_to_missing, is_missing(4), new_handler
942  real(FMS_HI_KIND_) :: f1, f2, f3, f4, middle, w, s
943  integer, parameter :: kindl = fms_hi_kind_
944 
945  num_missing = 0
946 
947  nlon_in = interp%nlon_src; nlat_in = interp%nlat_src
948  nlon_out = interp%nlon_dst; nlat_out = interp%nlat_dst
949 
950  if(present(mask_in)) then
951  mask = mask_in
952  else
953  mask = 1.0_kindl
954  endif
955 
956  if (present(verbose)) then
957  iverbose = verbose
958  else
959  iverbose = 0
960  endif
961 
962  if(present(missing_permit)) then
963  max_missing = missing_permit
964  else
965  max_missing = 0
966  endif
967 
968  if(present(new_handle_missing)) then
969  new_handler = new_handle_missing
970  else
971  new_handler = .false.
972  endif
973 
974  if(max_missing .gt. 3 .or. max_missing .lt. 0) call mpp_error(fatal, &
975  'horiz_interp_bilinear_mod: missing_permit should be between 0 and 3')
976 
977  if (size(data_in,1) /= nlon_in .or. size(data_in,2) /= nlat_in) &
978  call mpp_error(fatal,'horiz_interp_bilinear_mod: size of input array incorrect')
979 
980  if (size(data_out,1) /= nlon_out .or. size(data_out,2) /= nlat_out) &
981  call mpp_error(fatal,'horiz_interp_bilinear_mod: size of output array incorrect')
982 
983  if(new_handler) then
984  if( .not. present(missing_value) ) call mpp_error(fatal, &
985  "horiz_interp_bilinear_mod: misisng_value must be present when new_handle_missing is .true.")
986  if( present(mask_in) ) call mpp_error(fatal, &
987  "horiz_interp_bilinear_mod: mask_in should not be present when new_handle_missing is .true.")
988  do n = 1, nlat_out
989  do m = 1, nlon_out
990  is = interp % i_lon (m,n,1); ie = interp % i_lon (m,n,2)
991  js = interp % j_lat (m,n,1); je = interp % j_lat (m,n,2)
992  wtw = interp % HI_KIND_TYPE_ % wti (m,n,1)
993  wte = interp % HI_KIND_TYPE_ % wti (m,n,2)
994  wts = interp % HI_KIND_TYPE_ % wtj (m,n,1)
995  wtn = interp % HI_KIND_TYPE_ % wtj (m,n,2)
996 
997  is_missing = .false.
998  num_missing = 0
999  set_to_missing = .false.
1000  if(data_in(is,js) == missing_value) then
1001  num_missing = num_missing+1
1002  is_missing(1) = .true.
1003  if(wtw .GE. 0.5_kindl .AND. wts .GE. 0.5_kindl) set_to_missing = .true.
1004  endif
1005 
1006  if(data_in(ie,js) == missing_value) then
1007  num_missing = num_missing+1
1008  is_missing(2) = .true.
1009  if(wte .GE. 0.5_kindl .AND. wts .GE. 0.5_kindl ) set_to_missing = .true.
1010  endif
1011  if(data_in(ie,je) == missing_value) then
1012  num_missing = num_missing+1
1013  is_missing(3) = .true.
1014  if(wte .GE. 0.5_kindl .AND. wtn .GE. 0.5_kindl ) set_to_missing = .true.
1015  endif
1016  if(data_in(is,je) == missing_value) then
1017  num_missing = num_missing+1
1018  is_missing(4) = .true.
1019  if(wtw .GE. 0.5_kindl .AND. wtn .GE. 0.5_kindl) set_to_missing = .true.
1020  endif
1021 
1022  if( num_missing == 4 .OR. set_to_missing ) then
1023  data_out(m,n) = missing_value
1024  if(present(mask_out)) mask_out(m,n) = 0.0_kindl
1025  cycle
1026  else if(num_missing == 0) then
1027  f1 = data_in(is,js)
1028  f2 = data_in(ie,js)
1029  f3 = data_in(ie,je)
1030  f4 = data_in(is,je)
1031  w = wtw
1032  s = wts
1033  else if(num_missing == 3) then !--- three missing value
1034  if(.not. is_missing(1) ) then
1035  data_out(m,n) = data_in(is,js)
1036  else if(.not. is_missing(2) ) then
1037  data_out(m,n) = data_in(ie,js)
1038  else if(.not. is_missing(3) ) then
1039  data_out(m,n) = data_in(ie,je)
1040  else if(.not. is_missing(4) ) then
1041  data_out(m,n) = data_in(is,je)
1042  endif
1043  if(present(mask_out) ) mask_out(m,n) = 1.0_kindl
1044  cycle
1045  else !--- one or two missing value
1046  if( num_missing == 1) then
1047  if( is_missing(1) .OR. is_missing(3) ) then
1048  middle = 0.5_kindl *(data_in(ie,js)+data_in(is,je))
1049  else
1050  middle = 0.5_kindl *(data_in(is,js)+data_in(ie,je))
1051  endif
1052  else ! num_missing = 2
1053  if( is_missing(1) .AND. is_missing(2) ) then
1054  middle = 0.5_kindl *(data_in(ie,je)+data_in(is,je))
1055  else if( is_missing(1) .AND. is_missing(3) ) then
1056  middle = 0.5_kindl *(data_in(ie,js)+data_in(is,je))
1057  else if( is_missing(1) .AND. is_missing(4) ) then
1058  middle = 0.5_kindl *(data_in(ie,js)+data_in(ie,je))
1059  else if( is_missing(2) .AND. is_missing(3) ) then
1060  middle = 0.5_kindl *(data_in(is,js)+data_in(is,je))
1061  else if( is_missing(2) .AND. is_missing(4) ) then
1062  middle = 0.5_kindl*(data_in(is,js)+data_in(ie,je))
1063  else if( is_missing(3) .AND. is_missing(4) ) then
1064  middle = 0.5_kindl*(data_in(is,js)+data_in(ie,js))
1065  endif
1066  endif
1067 
1068  if( wtw .GE. 0.5_kindl .AND. wts .GE. 0.5_kindl ) then ! zone 1
1069  w = 2.0_kindl*(wtw-0.5_kindl)
1070  s = 2.0_kindl*(wts-0.5_kindl)
1071  f1 = data_in(is,js)
1072  if(is_missing(2)) then
1073  f2 = f1
1074  else
1075  f2 = 0.5_kindl*(data_in(is,js)+data_in(ie,js))
1076  endif
1077  f3 = middle
1078  if(is_missing(4)) then
1079  f4 = f1
1080  else
1081  f4 = 0.5_kindl*(data_in(is,js)+data_in(is,je))
1082  endif
1083  else if( wte .GE. 0.5_kindl .AND. wts .GE. 0.5_kindl ) then ! zone 2
1084  w = 2.0_kindl*(1.0_kindl-wte)
1085  s = 2.0_kindl*(wts-0.5_kindl)
1086  f2 = data_in(ie,js)
1087  if(is_missing(1)) then
1088  f1 = f2
1089  else
1090  f1 = 0.5_kindl*(data_in(is,js)+data_in(ie,js))
1091  endif
1092  f4 = middle
1093  if(is_missing(3)) then
1094  f3 = f2
1095  else
1096  f3 = 0.5_kindl*(data_in(ie,js)+data_in(ie,je))
1097  endif
1098  else if( wte .GE. 0.5_kindl .AND. wtn .GE. 0.5_kindl ) then ! zone 3
1099  w = 2.0_kindl*(1.0_kindl-wte)
1100  s = 2.0_kindl*(1.0_kindl-wtn)
1101  f3 = data_in(ie,je)
1102  if(is_missing(2)) then
1103  f2 = f3
1104  else
1105  f2 = 0.5_kindl*(data_in(ie,js)+data_in(ie,je))
1106  endif
1107  f1 = middle
1108  if(is_missing(4)) then
1109  f4 = f3
1110  else
1111  f4 = 0.5_kindl*(data_in(ie,je)+data_in(is,je))
1112  endif
1113  else if( wtw .GE. 0.5_kindl .AND. wtn .GE. 0.5_kindl ) then ! zone 4
1114  w = 2.0_kindl*(wtw-0.5_kindl)
1115  s = 2.0_kindl*(1.0_kindl-wtn)
1116  f4 = data_in(is,je)
1117  if(is_missing(1)) then
1118  f1 = f4
1119  else
1120  f1 = 0.5_kindl*(data_in(is,js)+data_in(is,je))
1121  endif
1122  f2 = middle
1123  if(is_missing(3)) then
1124  f3 = f4
1125  else
1126  f3 = 0.5_kindl*(data_in(ie,je)+data_in(is,je))
1127  endif
1128  else
1129  call mpp_error(fatal, &
1130  "horiz_interp_bilinear_mod: the point should be in one of the four zone")
1131  endif
1132  endif
1133 
1134  data_out(m,n) = f3 + (f4-f3)*w + (f2-f3)*s + ((f1-f2)+(f3-f4))*w*s
1135  if(present(mask_out)) mask_out(m,n) = 1.0_kindl
1136  enddo
1137  enddo
1138  else
1139  do n = 1, nlat_out
1140  do m = 1, nlon_out
1141  is = interp % i_lon (m,n,1); ie = interp % i_lon (m,n,2)
1142  js = interp % j_lat (m,n,1); je = interp % j_lat (m,n,2)
1143  wtw = interp % HI_KIND_TYPE_ % wti (m,n,1)
1144  wte = interp % HI_KIND_TYPE_ % wti (m,n,2)
1145  wts = interp % HI_KIND_TYPE_ % wtj (m,n,1)
1146  wtn = interp % HI_KIND_TYPE_ % wtj (m,n,2)
1147 
1148  if(present(missing_value) ) then
1149  num_missing = 0
1150  if(data_in(is,js) == missing_value) then
1151  num_missing = num_missing+1
1152  mask(is,js) = 0.0_kindl
1153  endif
1154  if(data_in(ie,js) == missing_value) then
1155  num_missing = num_missing+1
1156  mask(ie,js) = 0.0_kindl
1157  endif
1158  if(data_in(ie,je) == missing_value) then
1159  num_missing = num_missing+1
1160  mask(ie,je) = 0.0_kindl
1161  endif
1162  if(data_in(is,je) == missing_value) then
1163  num_missing = num_missing+1
1164  mask(is,je) = 0.0_kindl
1165  endif
1166  endif
1167 
1168  dwtsum = data_in(is,js)*mask(is,js)*wtw*wts &
1169  + data_in(ie,js)*mask(ie,js)*wte*wts &
1170  + data_in(ie,je)*mask(ie,je)*wte*wtn &
1171  + data_in(is,je)*mask(is,je)*wtw*wtn
1172  wtsum = mask(is,js)*wtw*wts + mask(ie,js)*wte*wts &
1173  + mask(ie,je)*wte*wtn + mask(is,je)*wtw*wtn
1174 
1175  if(.not. present(mask_in) .and. .not. present(missing_value)) wtsum = 1.0_kindl
1176 
1177  if(num_missing .gt. max_missing ) then
1178  data_out(m,n) = missing_value
1179  if(present(mask_out)) mask_out(m,n) = 0.0_kindl
1180  else if(wtsum .lt. real(epsln, fms_hi_kind_)) then
1181  if(present(missing_value)) then
1182  data_out(m,n) = missing_value
1183  else
1184  data_out(m,n) = 0.0_kindl
1185  endif
1186  if(present(mask_out)) mask_out(m,n) = 0.0_kindl
1187  else
1188  data_out(m,n) = dwtsum/wtsum
1189  if(present(mask_out)) mask_out(m,n) = wtsum
1190  endif
1191  enddo
1192  enddo
1193  endif
1194  !***********************************************************************
1195  ! compute statistics: minimum, maximum, and mean
1196  !-----------------------------------------------------------------------
1197  if (iverbose > 0) then
1198 
1199  ! compute statistics of input data
1200 
1201  call stats (data_in, min_in, max_in, avg_in, miss_in, missing_value, mask_in)
1202 
1203  ! compute statistics of output data
1204  call stats (data_out, min_out, max_out, avg_out, miss_out, missing_value, mask_out)
1205 
1206  !---- output statistics ----
1207  iunit = stdout()
1208  write (iunit,900)
1209  write (iunit,901) min_in ,max_in, avg_in
1210  if (present(mask_in)) write (iunit,903) miss_in
1211  write (iunit,902) min_out,max_out,avg_out
1212  if (present(mask_out)) write (iunit,903) miss_out
1213 
1214 900 format (/,1x,10('-'),' output from horiz_interp ',10('-'))
1215 901 format (' input: min=',f16.9,' max=',f16.9,' avg=',f22.15)
1216 902 format (' output: min=',f16.9,' max=',f16.9,' avg=',f22.15)
1217 903 format (' number of missing points = ',i6)
1218 
1219  endif
1220 
1221  return
1222 
1223  end subroutine
1224 
1225  !> Subroutine for reading a weight file and use it to fill in the horiz interp type
1226  !! for the bilinear interpolation method.
1227  subroutine horiz_interp_read_weights_bilinear_(Interp, weight_filename, lon_out, lat_out, lon_in, lat_in, &
1228  weight_file_source, interp_method, isw, iew, jsw, jew, nglon, nglat)
1229  type(horiz_interp_type), intent(inout) :: Interp !< Horiz interp time to fill
1230  character(len=*), intent(in) :: weight_filename !< Name of the weight file
1231  real(FMS_HI_KIND_), target, intent(in) :: lat_out(:,:) !< Output (model) latitude
1232  real(FMS_HI_KIND_), target, intent(in) :: lon_out(:,:) !< Output (model) longitude
1233  real(FMS_HI_KIND_), intent(in) :: lat_in(:) !< Input (data) latitude
1234  real(FMS_HI_KIND_), intent(in) :: lon_in(:) !< Input (data) longitude
1235  character(len=*), intent(in) :: weight_file_source !< Source of the weight file
1236  character(len=*), intent(in) :: interp_method !< The interp method to use
1237  integer, intent(in) :: isw, iew, jsw, jew !< Starting and ending indices of the compute domain
1238  integer, intent(in) :: nglon !< Number of longitudes in the global domain
1239  integer, intent(in) :: nglat !< Number of latitudes in the globl domain
1240 
1241 
1242  real(FMS_HI_KIND_), allocatable :: var(:,:,:) !< Dummy variable to read the indices and weight into
1243  type(fmsnetcdffile_t) :: weight_fileobj !< FMS2io fileob for the weight file
1244  integer :: nlon !< Number of longitudes in the model grid as read
1245  !! from the weight file
1246  integer :: nlat !< Number of latitude in the model grid as read
1247  !! from the weight file
1248 
1249  if (.not. open_file(weight_fileobj, weight_filename, "read" )) &
1250  call mpp_error(fatal, "Error opening the weight file:"//&
1251  &trim(weight_filename))
1252 
1253  !< Check that weight file has the correct dimensions
1254  select case (trim(weight_file_source))
1255  case ("fregrid")
1256  call get_dimension_size(weight_fileobj, "nlon", nlon)
1257  if (nlon .ne. nglon) &
1258  call mpp_error(fatal, "The nlon from the weight file is not the same as in the input grid."//&
1259  &" From weight file:"//string(nlon)//" from input grid:"//string(size(lon_out,1)))
1260  call get_dimension_size(weight_fileobj, "nlat", nlat)
1261  if (nlat .ne. nglat) &
1262  call mpp_error(fatal, "The nlat from the weight file is not the same as in the input grid."//&
1263  &" From weight file:"//string(nlat)//" from input grid:"//string(size(lon_out,2)))
1264  case default
1265  call mpp_error(fatal, trim(weight_file_source)//&
1266  &" is not a supported weight file source. fregrid is the only supported weight file source." )
1267  end select
1268 
1269  interp%nlon_src = size(lon_in(:)) ; interp%nlat_src = size(lat_in(:))
1270  interp%nlon_dst = size(lon_out,1); interp%nlat_dst = size(lon_out,2)
1271 
1272  allocate ( interp % HI_KIND_TYPE_ % wti (interp%nlon_dst,interp%nlat_dst,2), &
1273  interp % HI_KIND_TYPE_ % wtj (interp%nlon_dst,interp%nlat_dst,2), &
1274  interp % i_lon (interp%nlon_dst,interp%nlat_dst,2), &
1275  interp % j_lat (interp%nlon_dst,interp%nlat_dst,2))
1276 
1277 
1278  !! Three is for lon, lat, tile
1279  !! Currently, interpolation is only supported from lat,lon input data
1280  allocate(var(interp%nlon_dst,interp%nlat_dst, 3))
1281  call read_data(weight_fileobj, "index", var, corner=(/isw, jsw, 1/), edge_lengths=(/iew-isw+1, jew-jsw+1, 3/))
1282 
1283  !! Each point has a lon (i), and lat(j) index
1284  !! From there the four corners are (i,j), (i,j+1) (i+1) (i+1,j+1)
1285  interp % i_lon (:,:,1) = var(:,:,1)
1286  interp % i_lon (:,:,2) = interp % i_lon (:,:,1) + 1
1287  where (interp % i_lon (:,:,2) > size(lon_in(:))) interp % i_lon (:,:,2) = 1
1288 
1289  interp % j_lat (:,:,1) = var(:,:,2)
1290  interp % j_lat (:,:,2) = interp % j_lat (:,:,1) + 1
1291  where (interp % j_lat (:,:,2) > size(lat_in(:))) interp % j_lat (:,:,2) = 1
1292 
1293  deallocate(var)
1294 
1295  allocate(var(interp%nlon_dst,interp%nlat_dst, 4))
1296  call read_data(weight_fileobj, "weight", var, corner=(/isw, jsw, 1/), edge_lengths=(/iew-isw+1, jew-jsw+1, 4/))
1297 
1298  !! The weights for the four corners
1299  !! var(:,:,1) -> (i,j)
1300  !! var(:,:,2) -> (i,j+1)
1301  !! var(:,:,3) -> (i+1,j)
1302  !! var(:,:,4) -> (i+1,j+1)
1303  interp % HI_KIND_TYPE_ % wti = var(:,:,1:2)
1304  interp % HI_KIND_TYPE_ % wtj = var(:,:,3:4)
1305  deallocate(var)
1306 
1307  interp% HI_KIND_TYPE_ % is_allocated = .true.
1308  interp% interp_method = bilinear
1309  interp% I_am_initialized = .true.
1310  call close_file(weight_fileobj)
1311  end subroutine horiz_interp_read_weights_bilinear_
1312 
1313 !> @}
find_neighbor_
subroutine find_neighbor_(Interp, lon_in, lat_in, lon_out, lat_out, src_modulo)
this routine will search the source grid to fine the grid box that encloses each destination grid.
Definition: horiz_interp_bilinear.inc:421
horiz_interp_bilinear_new_2d_
subroutine horiz_interp_bilinear_new_2d_(Interp, lon_in, lat_in, lon_out, lat_out, verbose, src_modulo, new_search, no_crash_when_not_found)
Initialization routine.
Definition: horiz_interp_bilinear.inc:208
find_neighbor_new_
subroutine find_neighbor_new_(Interp, lon_in, lat_in, lon_out, lat_out, src_modulo, no_crash)
this routine will search the source grid to fine the grid box that encloses each destination grid.
Definition: horiz_interp_bilinear.inc:699
horiz_interp_read_weights_bilinear_
subroutine horiz_interp_read_weights_bilinear_(Interp, weight_filename, lon_out, lat_out, lon_in, lat_in, weight_file_source, interp_method, isw, iew, jsw, jew, nglon, nglat)
Subroutine for reading a weight file and use it to fill in the horiz interp type for the bilinear int...
Definition: horiz_interp_bilinear.inc:1229
horiz_interp_bilinear_
subroutine horiz_interp_bilinear_(Interp, data_in, data_out, verbose, mask_in, mask_out, missing_value, missing_permit, new_handle_missing)
Subroutine for performing the horizontal interpolation between two grids.
Definition: horiz_interp_bilinear.inc:917
inside_polygon_
logical function inside_polygon_(polyx, polyy, x, y)
The function will return true if the point x,y is inside a polygon, or false if it is not....
Definition: horiz_interp_bilinear.inc:667
stdout
integer function stdout()
This function returns the current standard fortran unit numbers for output.
Definition: mpp_util.inc:42