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