FV3 Bundle
horiz_interp_bilinear.F90
Go to the documentation of this file.
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 module horiz_interp_bilinear_mod
20 
21  ! <CONTACT EMAIL="Zhi.Liang@noaa.gov"> Zhi Liang </CONTACT>
22 
23  ! <HISTORY SRC="http://www.gfdl.noaa.gov/fms-cgi-bin/cvsweb.cgi/FMS/"/>
24 
25  ! <OVERVIEW>
26  ! Performs spatial interpolation between grids using bilinear interpolation
27  ! </OVERVIEW>
28 
29  ! <DESCRIPTION>
30  ! This module can interpolate data from regular rectangular grid
31  ! to rectangular/tripolar grid. The interpolation scheme is bilinear interpolation.
32  ! There is an optional mask field for missing input data.
33  ! An optional output mask field may be used in conjunction with
34  ! the input mask to show where output data exists.
35  ! </DESCRIPTION>
36 
37  use mpp_mod, only: mpp_error, fatal, stdout, mpp_pe, mpp_root_pe
38  use fms_mod, only: write_version_number
39  use constants_mod, only: pi
40  use horiz_interp_type_mod, only: horiz_interp_type, stats
41 
42  implicit none
43  private
44 
45 
46  public :: horiz_interp_bilinear_new, horiz_interp_bilinear, horiz_interp_bilinear_del
47  public :: horiz_interp_bilinear_init
48 
49  !--- public interface
50  interface horiz_interp_bilinear_new
51  module procedure horiz_interp_bilinear_new_1d
52  module procedure horiz_interp_bilinear_new_2d
53  end interface
54 
55 
56  real, parameter :: epsln=1.e-10
57  integer, parameter :: dummy = -999
58 
59  !-----------------------------------------------------------------------
60 ! Include variable "version" to be written to log file.
61 #include<file_version.h>
62  logical :: module_is_initialized = .false.
63 
64 contains
65 
66  !#######################################################################
67  ! <SUBROUTINE NAME="horiz_interp_bilinear_init">
68  ! <OVERVIEW>
69  ! writes version number to logfile.out
70  ! </OVERVIEW>
71  ! <DESCRIPTION>
72  ! writes version number to logfile.out
73  ! </DESCRIPTION>
74 
75  subroutine horiz_interp_bilinear_init
76 
77  if(module_is_initialized) return
78  call write_version_number("HORIZ_INTERP_BILINEAR_MOD", version)
79  module_is_initialized = .true.
80 
81  end subroutine horiz_interp_bilinear_init
82 
83  ! </SUBROUTINE>
84 
85  !########################################################################
86 
87  subroutine horiz_interp_bilinear_new_1d ( Interp, lon_in, lat_in, lon_out, lat_out, &
88  verbose, src_modulo )
89 
90  !-----------------------------------------------------------------------
91  type(horiz_interp_type), intent(inout) :: Interp
92  real, intent(in), dimension(:) :: lon_in , lat_in
93  real, intent(in), dimension(:,:) :: lon_out, lat_out
94  integer, intent(in), optional :: verbose
95  logical, intent(in), optional :: src_modulo
96 
97  logical :: src_is_modulo
98  integer :: nlon_in, nlat_in, nlon_out, nlat_out, n, m
99  integer :: ie, is, je, js, ln_err, lt_err, warns, unit
100  real :: wtw, wte, wts, wtn, lon, lat, tpi, hpi
101  real :: glt_min, glt_max, gln_min, gln_max, min_lon, max_lon
102 
103  warns = 0
104  if(present(verbose)) warns = verbose
105  src_is_modulo = .true.
106  if (present(src_modulo)) src_is_modulo = src_modulo
107 
108  hpi = 0.5*pi
109  tpi = 4.0*hpi
110  glt_min = hpi
111  glt_max = -hpi
112  gln_min = tpi
113  gln_max = -tpi
114  min_lon = 0.0
115  max_lon = tpi
116  ln_err = 0
117  lt_err = 0
118  !-----------------------------------------------------------------------
119 
120  allocate ( interp % wti (size(lon_out,1),size(lon_out,2),2), &
121  interp % wtj (size(lon_out,1),size(lon_out,2),2), &
122  interp % i_lon (size(lon_out,1),size(lon_out,2),2), &
123  interp % j_lat (size(lon_out,1),size(lon_out,2),2))
124  !-----------------------------------------------------------------------
125 
126  nlon_in = size(lon_in(:)) ; nlat_in = size(lat_in(:))
127  nlon_out = size(lon_out, 1); nlat_out = size(lon_out, 2)
128  interp%nlon_src = nlon_in; interp%nlat_src = nlat_in
129  interp%nlon_dst = nlon_out; interp%nlat_dst = nlat_out
130 
131  if(src_is_modulo) then
132  if(lon_in(nlon_in) - lon_in(1) .gt. tpi + epsln) &
133  call mpp_error(fatal,'horiz_interp_bilinear_mod: '// &
134  'The range of source grid longitude should be no larger than tpi')
135 
136  if(lon_in(1) .lt. 0.0 .OR. lon_in(nlon_in) > tpi ) then
137  min_lon = lon_in(1)
138  max_lon = lon_in(nlon_in)
139  endif
140  endif
141 
142  do n = 1, nlat_out
143  do m = 1, nlon_out
144  lon = lon_out(m,n)
145  lat = lat_out(m,n)
146 
147  if(src_is_modulo) then
148  if(lon .lt. min_lon) then
149  lon = lon + tpi
150  else if(lon .gt. max_lon) then
151  lon = lon - tpi
152  endif
153  else ! when the input grid is in not cyclic, the output grid should located inside
154  ! the input grid
155  if((lon .lt. lon_in(1)) .or. (lon .gt. lon_in(nlon_in))) &
156  call mpp_error(fatal,'horiz_interp_bilinear_mod: ' //&
157  'when input grid is not modulo, output grid should locate inside input grid')
158  endif
159 
160  glt_min = min(lat,glt_min); glt_max = max(lat,glt_max)
161  gln_min = min(lon,gln_min); gln_max = max(lon,gln_max)
162 
163  is = indp(lon, lon_in )
164  if( lon_in(is) .gt. lon ) is = max(is-1,1)
165  if( lon_in(is) .eq. lon .and. is .eq. nlon_in) is = max(is - 1,1)
166  ie = min(is+1,nlon_in)
167  if(lon_in(is) .ne. lon_in(ie) .and. lon_in(is) .le. lon) then
168  wtw = ( lon_in(ie) - lon) / (lon_in(ie) - lon_in(is) )
169  else
170  ! east or west of the last data value. this could be because a
171  ! cyclic condition is needed or the dataset is too small.
172  ln_err = 1
173  ie = 1
174  is = nlon_in
175  if (lon_in(ie) .ge. lon ) then
176  wtw = (lon_in(ie) -lon)/(lon_in(ie)-lon_in(is)+tpi+epsln)
177  else
178  wtw = (lon_in(ie) -lon+tpi+epsln)/(lon_in(ie)-lon_in(is)+tpi+epsln)
179  endif
180  endif
181  wte = 1. - wtw
182 
183  js = indp(lat, lat_in )
184 
185  if( lat_in(js) .gt. lat ) js = max(js - 1, 1)
186  if( lat_in(js) .eq. lat .and. js .eq. nlat_in) js = max(js - 1, 1)
187  je = min(js + 1, nlat_in)
188 
189  if ( lat_in(js) .ne. lat_in(je) .and. lat_in(js) .le. lat) then
190  wts = ( lat_in(je) - lat )/(lat_in(je)-lat_in(js))
191  else
192  ! north or south of the last data value. this could be because a
193  ! pole is not included in the data set or the dataset is too small.
194  ! in either case extrapolate north or south
195  lt_err = 1
196  wts = 1.
197  endif
198 
199  wtn = 1. - wts
200 
201  interp % i_lon (m,n,1) = is; interp % i_lon (m,n,2) = ie
202  interp % j_lat (m,n,1) = js; interp % j_lat (m,n,2) = je
203  interp % wti (m,n,1) = wtw
204  interp % wti (m,n,2) = wte
205  interp % wtj (m,n,1) = wts
206  interp % wtj (m,n,2) = wtn
207 
208  enddo
209  enddo
210 
211  unit = stdout()
212 
213  if (ln_err .eq. 1 .and. warns > 0) then
214  write (unit,'(/,(1x,a))') &
215  '==> Warning: the geographic data set does not extend far ', &
216  ' enough east or west - a cyclic boundary ', &
217  ' condition was applied. check if appropriate '
218  write (unit,'(/,(1x,a,2f8.4))') &
219  ' data required between longitudes:', gln_min, gln_max, &
220  ' data set is between longitudes:', lon_in(1), lon_in(nlon_in)
221  warns = warns - 1
222  endif
223 
224  if (lt_err .eq. 1 .and. warns > 0) then
225  write (unit,'(/,(1x,a))') &
226  '==> Warning: the geographic data set does not extend far ',&
227  ' enough north or south - extrapolation from ',&
228  ' the nearest data was applied. this may create ',&
229  ' artificial gradients near a geographic pole '
230  write (unit,'(/,(1x,a,2f8.4))') &
231  ' data required between latitudes:', glt_min, glt_max, &
232  ' data set is between latitudes:', lat_in(1), lat_in(nlat_in)
233  endif
234 
235  return
236 
237  end subroutine horiz_interp_bilinear_new_1d
238 
239  !#######################################################################
240  ! <SUBROUTINE NAME="horiz_interp_bilinear_new">
241 
242  ! <OVERVIEW>
243  ! Initialization routine.
244  ! </OVERVIEW>
245  ! <DESCRIPTION>
246  ! Allocates space and initializes a derived-type variable
247  ! that contains pre-computed interpolation indices and weights.
248  ! </DESCRIPTION>
249  ! <TEMPLATE>
250  ! call horiz_interp_bilinear_new ( Interp, lon_in, lat_in, lon_out, lat_out, verbose, src_modulo )
251 
252  ! </TEMPLATE>
253  !
254  ! <IN NAME="lon_in" TYPE="real, dimension(:,:)" UNITS="radians">
255  ! Longitude (in radians) for source data grid.
256  ! </IN>
257 
258  ! <IN NAME="lat_in" TYPE="real, dimension(:,:)" UNITS="radians">
259  ! Latitude (in radians) for source data grid.
260  ! </IN>
261 
262  ! <IN NAME="lon_out" TYPE="real, dimension(:,:)" UNITS="radians" >
263  ! Longitude (in radians) for source data grid.
264  ! </IN>
265 
266  ! <IN NAME="lat_out" TYPE="real, dimension(:,:)" UNITS="radians" >
267  ! Latitude (in radians) for source data grid.
268  ! </IN>
269 
270  ! <IN NAME="src_modulo" TYPE="logical, optional">
271  ! logical variable to indicate if the boundary condition along zonal boundary
272  ! is cyclic or not. When true, the zonal boundary condition is cyclic.
273  ! </IN>
274 
275  ! <IN NAME="verbose" TYPE="integer, optional" >
276  ! flag for the amount of print output.
277  ! </IN>
278 
279  ! <INOUT NAME="Interp" TYPE="type(horiz_interp_type)" >
280  ! A derived-type variable containing indices and weights used for subsequent
281  ! interpolations. To reinitialize this variable for a different grid-to-grid
282  ! interpolation you must first use the "horiz_interp_del" interface.
283  ! </INOUT>
284 
285  subroutine horiz_interp_bilinear_new_2d ( Interp, lon_in, lat_in, lon_out, lat_out, &
286  verbose, src_modulo, new_search, no_crash_when_not_found )
287 
288  !-----------------------------------------------------------------------
289  type(horiz_interp_type), intent(inout) :: Interp
290  real, intent(in), dimension(:,:) :: lon_in , lat_in
291  real, intent(in), dimension(:,:) :: lon_out, lat_out
292  integer, intent(in), optional :: verbose
293  logical, intent(in), optional :: src_modulo
294  logical, intent(in), optional :: new_search
295  logical, intent(in), optional :: no_crash_when_not_found
296  integer :: warns
297  logical :: src_is_modulo
298  integer :: nlon_in, nlat_in, nlon_out, nlat_out
299  integer :: m, n, is, ie, js, je, num_solution
300  real :: lon, lat, quadra, x, y, y1, y2
301  real :: a1, b1, c1, d1, a2, b2, c2, d2, a, b, c
302  real :: lon1, lat1, lon2, lat2, lon3, lat3, lon4, lat4
303  real :: tpi, lon_min, lon_max
304  real :: epsln2
305  logical :: use_new_search, no_crash
306 
307  tpi = 2.0*pi
308 
309  warns = 0
310  if(present(verbose)) warns = verbose
311  src_is_modulo = .true.
312  if (present(src_modulo)) src_is_modulo = src_modulo
313  use_new_search = .false.
314  if (present(new_search)) use_new_search = new_search
315  no_crash = .false.
316  if(present(no_crash_when_not_found)) no_crash = no_crash_when_not_found
317 
318  ! make sure lon and lat has the same dimension
319  if(size(lon_out,1) /= size(lat_out,1) .or. size(lon_out,2) /= size(lat_out,2) ) &
320  call mpp_error(fatal,'horiz_interp_bilinear_mod: when using bilinear ' // &
321  'interplation, the output grids should be geographical grids')
322 
323  if(size(lon_in,1) /= size(lat_in,1) .or. size(lon_in,2) /= size(lat_in,2) ) &
324  call mpp_error(fatal,'horiz_interp_bilinear_mod: when using bilinear '// &
325  'interplation, the input grids should be geographical grids')
326 
327  !--- get the grid size
328  nlon_in = size(lon_in,1) ; nlat_in = size(lat_in,2)
329  nlon_out = size(lon_out,1); nlat_out = size(lon_out,2)
330  interp%nlon_src = nlon_in; interp%nlat_src = nlat_in
331  interp%nlon_dst = nlon_out; interp%nlat_dst = nlat_out
332 
333  allocate ( interp % wti (size(lon_out,1),size(lon_out,2),2), &
334  interp % wtj (size(lon_out,1),size(lon_out,2),2), &
335  interp % i_lon (size(lon_out,1),size(lon_out,2),2), &
336  interp % j_lat (size(lon_out,1),size(lon_out,2),2))
337 
338  !--- first fine the neighbor points for the destination points.
339  if(use_new_search) then
340  epsln2 = epsln*1e5
341  call find_neighbor_new(interp, lon_in, lat_in, lon_out, lat_out, src_is_modulo, no_crash)
342  else
343  epsln2 = epsln
344  call find_neighbor(interp, lon_in, lat_in, lon_out, lat_out, src_is_modulo)
345  endif
346 
347  !***************************************************************************
348  ! Algorithm explanation (from disscussion with Steve Garner ) *
349  ! *
350  ! lon(x,y) = a1*x + b1*y + c1*x*y + d1 (1) *
351  ! lat(x,y) = a2*x + b2*y + c2*x*y + d2 (2) *
352  ! f (x,y) = a3*x + b3*y + c3*x*y + d3 (3) *
353  ! with x and y is between 0 and 1. *
354  ! lon1 = lon(0,0) = d1, lat1 = lat(0,0) = d2 *
355  ! lon2 = lon(1,0) = a1+d1, lat2 = lat(1,0) = a2+d2 *
356  ! lon3 = lon(1,1) = a1+b1+c1+d1, lat3 = lat(1,1) = a2+b2+c2+d2 *
357  ! lon4 = lon(0,1) = b1+d1, lat4 = lat(0,1) = b2+d2 *
358  ! where (lon1,lat1),(lon2,lat2),(lon3,lat3),(lon4,lat4) represents *
359  ! the four corners starting from the left lower corner of grid box *
360  ! that encloses a destination grid ( the rotation direction is *
361  ! counterclockwise ). With these conditions, we get *
362  ! a1 = lon2-lon1, a2 = lat2-lat1 *
363  ! b1 = lon4-lon1, b2 = lat4-lat1 *
364  ! c1 = lon3-lon2-lon4+lon1, c2 = lat3-lat2-lat4+lat1 *
365  ! d1 = lon1 d2 = lat1 *
366  ! So given any point (lon,lat), from equation (1) and (2) we can *
367  ! solve (x,y). *
368  ! From equation (3) *
369  ! f1 = f(0,0) = d3, f2 = f(1,0) = a3+d3 *
370  ! f3 = f(1,1) = a3+b3+c3+d3, f4 = f(0,1) = b3+d3 *
371  ! we obtain *
372  ! a3 = f2-f1, b3 = f4-f1 *
373  ! c3 = f3-f2-f4+f1, d3 = f1 *
374  ! at point (lon,lat) ---> (x,y) *
375  ! f(x,y) = (f2-f1)x + (f4-f1)y + (f3-f2-f4+f1)xy + f1 *
376  ! = f1*(1-x)*(1-y) + f2*x*(1-y) + f3*x*y + f4*y*(1-x) *
377  ! wtw=1-x; wte=x; wts=1-y; xtn=y *
378  ! *
379  !***************************************************************************
380 
381  lon_min = minval(lon_in);
382  lon_max = maxval(lon_in);
383  !--- calculate the weight
384  do n = 1, nlat_out
385  do m = 1, nlon_out
386  lon = lon_out(m,n)
387  lat = lat_out(m,n)
388  if(lon .lt. lon_min) then
389  lon = lon + tpi
390  else if(lon .gt. lon_max) then
391  lon = lon - tpi
392  endif
393  is = interp%i_lon(m,n,1); ie = interp%i_lon(m,n,2)
394  js = interp%j_lat(m,n,1); je = interp%j_lat(m,n,2)
395  if( is == dummy) cycle
396  lon1 = lon_in(is,js); lat1 = lat_in(is,js);
397  lon2 = lon_in(ie,js); lat2 = lat_in(ie,js);
398  lon3 = lon_in(ie,je); lat3 = lat_in(ie,je);
399  lon4 = lon_in(is,je); lat4 = lat_in(is,je);
400  if(lon .lt. lon_min) then
401  lon1 = lon1 -tpi; lon4 = lon4 - tpi
402  else if(lon .gt. lon_max) then
403  lon2 = lon2 +tpi; lon3 = lon3 + tpi
404  endif
405  a1 = lon2-lon1
406  b1 = lon4-lon1
407  c1 = lon1+lon3-lon4-lon2
408  d1 = lon1
409  a2 = lat2-lat1
410  b2 = lat4-lat1
411  c2 = lat1+lat3-lat4-lat2
412  d2 = lat1
413  !--- the coefficient of the quadratic equation
414  a = b2*c1-b1*c2
415  b = a1*b2-a2*b1+c1*d2-c2*d1+c2*lon-c1*lat
416  c = a2*lon-a1*lat+a1*d2-a2*d1
417  quadra = b*b-4.*a*c
418  if(abs(quadra) < epsln) quadra = 0.0
419  if(quadra < 0.0) call mpp_error(fatal, &
420  "horiz_interp_bilinear_mod: No solution existed for this quadratic equation")
421  if ( abs(a) .lt. epsln2) then ! a = 0 is a linear equation
422  if( abs(b) .lt. epsln) call mpp_error(fatal, &
423  "horiz_interp_bilinear_mod: no unique solution existed for this linear equation")
424  y = -c/b
425  else
426  y1 = 0.5*(-b+sqrt(quadra))/a
427  y2 = 0.5*(-b-sqrt(quadra))/a
428  if(abs(y1) < epsln2) y1 = 0.0
429  if(abs(y2) < epsln2) y2 = 0.0
430  if(abs(1.0-y1) < epsln2) y1 = 1.0
431  if(abs(1.0-y2) < epsln2) y2 = 1.0
432  num_solution = 0
433  if(y1 >= 0.0 .and. y1 <= 1.0) then
434  y = y1
435  num_solution = num_solution +1
436  endif
437  if(y2 >= 0.0 .and. y2 <= 1.0) then
438  y = y2
439  num_solution = num_solution + 1
440  endif
441  if(num_solution == 0) then
442  call mpp_error(fatal, "horiz_interp_bilinear_mod: No solution found")
443  else if(num_solution == 2) then
444  call mpp_error(fatal, "horiz_interp_bilinear_mod: Two solutions found")
445  endif
446  endif
447  if(abs(a1+c1*y) < epsln) call mpp_error(fatal, &
448  "horiz_interp_bilinear_mod: the denomenator is 0")
449  if(abs(y) < epsln2) y = 0.0
450  if(abs(1.0-y) < epsln2) y = 1.0
451  x = (lon-b1*y-d1)/(a1+c1*y)
452  if(abs(x) < epsln2) x = 0.0
453  if(abs(1.0-x) < epsln2) x = 1.0
454  ! x and y should be between 0 and 1.
455  !! Added for ECDA
456  if(use_new_search) then
457  if (x < 0.0) x = 0.0 ! snz
458  if (y < 0.0) y = 0.0 ! snz
459  if (x > 1.0) x = 1.0
460  if (y > 1.0) y = 1.0
461  endif
462  if( x>1. .or. x<0. .or. y>1. .or. y < 0.) call mpp_error(fatal, &
463  "horiz_interp_bilinear_mod: weight should be between 0 and 1")
464  interp % wti(m,n,1)=1.0-x; interp % wti(m,n,2)=x
465  interp % wtj(m,n,1)=1.0-y; interp % wtj(m,n,2)=y
466  enddo
467  enddo
468 
469  end subroutine horiz_interp_bilinear_new_2d
470  ! </SUBROUTINE>
471 
472  !#######################################################################
473  ! this routine will search the source grid to fine the grid box that encloses
474  ! each destination grid.
475  subroutine find_neighbor( Interp, lon_in, lat_in, lon_out, lat_out, src_modulo )
476  type(horiz_interp_type), intent(inout) :: Interp
477  real, intent(in), dimension(:,:) :: lon_in , lat_in
478  real, intent(in), dimension(:,:) :: lon_out, lat_out
479  logical, intent(in) :: src_modulo
480  integer :: nlon_in, nlat_in, nlon_out, nlat_out
481  integer :: max_step, n, m, l, i, j, ip1, jp1, step
482  integer :: is, js, jstart, jend, istart, iend, npts
483  integer, allocatable, dimension(:) :: ilon, jlat
484  real :: lon_min, lon_max, lon, lat, tpi
485  logical :: found
486  real :: lon1, lat1, lon2, lat2, lon3, lat3, lon4, lat4
487 
488  tpi = 2.0*pi
489  nlon_in = size(lon_in,1) ; nlat_in = size(lat_in,2)
490  nlon_out = size(lon_out,1); nlat_out = size(lon_out,2)
491 
492  lon_min = minval(lon_in);
493  lon_max = maxval(lon_in);
494 
495  max_step = max(nlon_in,nlat_in) ! can be adjusted if needed
496  allocate(ilon(5*max_step), jlat(5*max_step) )
497 
498  do n = 1, nlat_out
499  do m = 1, nlon_out
500  found = .false.
501  lon = lon_out(m,n)
502  lat = lat_out(m,n)
503 
504  if(src_modulo) then
505  if(lon .lt. lon_min) then
506  lon = lon + tpi
507  else if(lon .gt. lon_max) then
508  lon = lon - tpi
509  endif
510  else
511  if(lon .lt. lon_min .or. lon .gt. lon_max ) &
512  call mpp_error(fatal,'horiz_interp_bilinear_mod: ' //&
513  'when input grid is not modulo, output grid should locate inside input grid')
514  endif
515  !--- search for the surrounding four points locatioon.
516  if(m==1 .and. n==1) then
517  j_loop: do j = 1, nlat_in-1
518  do i = 1, nlon_in
519  ip1 = i+1
520  jp1 = j+1
521  if(i==nlon_in) then
522  if(src_modulo)then
523  ip1 = 1
524  else
525  cycle
526  endif
527  endif
528  lon1 = lon_in(i, j); lat1 = lat_in(i,j)
529  lon2 = lon_in(ip1,j); lat2 = lat_in(ip1,j)
530  lon3 = lon_in(ip1,jp1); lat3 = lat_in(ip1,jp1)
531  lon4 = lon_in(i, jp1); lat4 = lat_in(i, jp1)
532 
533  if(lon .lt. lon_min .or. lon .gt. lon_max) then
534  if(i .ne. nlon_in) then
535  cycle
536  else
537  if(lon .lt. lon_min) then
538  lon1 = lon1 -tpi; lon4 = lon4 - tpi
539  else if(lon .gt. lon_max) then
540  lon2 = lon2 +tpi; lon3 = lon3 + tpi
541  endif
542  endif
543  endif
544 
545  if(lat .ge. intersect(lon1,lat1,lon2,lat2,lon))then ! south
546  if(lon .le. intersect(lat2,lon2,lat3,lon3,lat))then ! east
547  if(lat .le. intersect(lon3,lat3,lon4,lat4,lon))then ! north
548  if(lon .ge. intersect(lat4,lon4,lat1,lon1,lat))then ! west
549  found = .true.
550  interp % i_lon (m,n,1) = i; interp % i_lon (m,n,2) = ip1
551  interp % j_lat (m,n,1) = j; interp % j_lat (m,n,2) = jp1
552  exit j_loop
553  endif
554  endif
555  endif
556  endif
557  enddo
558  enddo j_loop
559  else
560  step = 0
561  do while ( .not. found .and. step .lt. max_step )
562  !--- take the adajcent point as the starting point
563  if(m == 1) then
564  is = interp % i_lon (m,n-1,1)
565  js = interp % j_lat (m,n-1,1)
566  else
567  is = interp % i_lon (m-1,n,1)
568  js = interp % j_lat (m-1,n,1)
569  endif
570  if(step==0) then
571  npts = 1
572  ilon(1) = is
573  jlat(1) = js
574  else
575  npts = 0
576  !--- bottom boundary
577  jstart = max(js-step,1)
578  jend = min(js+step,nlat_in)
579 
580  do l = -step, step
581  i = is+l
582  if(src_modulo)then
583  if( i < 1) then
584  i = i + nlon_in
585  else if (i > nlon_in) then
586  i = i - nlon_in
587  endif
588  if( i < 1 .or. i > nlon_in) call mpp_error(fatal, &
589  'horiz_interp_bilinear_mod: max_step is too big, decrease max_step' )
590  else
591  if( i < 1 .or. i > nlon_in) cycle
592  endif
593 
594  npts = npts + 1
595  ilon(npts) = i
596  jlat(npts) = jstart
597  enddo
598 
599  !--- right and left boundary -----------------------------------------------
600  istart = is - step
601  iend = is + step
602  if(src_modulo) then
603  if( istart < 1) istart = istart + nlon_in
604  if( iend > nlon_in) iend = iend - nlon_in
605  else
606  istart = max(istart,1)
607  iend = min(iend, nlon_in)
608  endif
609  do l = -step, step
610  j = js+l
611  if( j < 1 .or. j > nlat_in .or. j==jstart .or. j==jend) cycle
612  npts = npts+1
613  ilon(npts) = istart
614  jlat(npts) = j
615  npts = npts+1
616  ilon(npts) = iend
617  jlat(npts) = j
618  end do
619 
620  !--- top boundary
621 
622  do l = -step, step
623  i = is+l
624  if(src_modulo)then
625  if( i < 1) then
626  i = i + nlon_in
627  else if (i > nlon_in) then
628  i = i - nlon_in
629  endif
630  if( i < 1 .or. i > nlon_in) call mpp_error(fatal, &
631  'horiz_interp_bilinear_mod: max_step is too big, decrease max_step' )
632  else
633  if( i < 1 .or. i > nlon_in) cycle
634  endif
635 
636  npts = npts + 1
637  ilon(npts) = i
638  jlat(npts) = jend
639  enddo
640 
641 
642  end if
643 
644  !--- find the surrouding points
645  do l = 1, npts
646  i = ilon(l)
647  j = jlat(l)
648  ip1 = i+1
649  if(ip1>nlon_in) then
650  if(src_modulo) then
651  ip1 = 1
652  else
653  cycle
654  endif
655  endif
656  jp1 = j+1
657  if(jp1>nlat_in) cycle
658  lon1 = lon_in(i, j); lat1 = lat_in(i,j)
659  lon2 = lon_in(ip1,j); lat2 = lat_in(ip1,j)
660  lon3 = lon_in(ip1,jp1); lat3 = lat_in(ip1,jp1)
661  lon4 = lon_in(i, jp1); lat4 = lat_in(i, jp1)
662 
663  if(lon .lt. lon_min .or. lon .gt. lon_max) then
664  if(i .ne. nlon_in) then
665  cycle
666  else
667  if(lon .lt. lon_min) then
668  lon1 = lon1 -tpi; lon4 = lon4 - tpi
669  else if(lon .gt. lon_max) then
670  lon2 = lon2 +tpi; lon3 = lon3 + tpi
671  endif
672  endif
673  endif
674 
675  if(lat .ge. intersect(lon1,lat1,lon2,lat2,lon))then ! south
676  if(lon .le. intersect(lat2,lon2,lat3,lon3,lat))then ! east
677  if(lat .le. intersect(lon3,lat3,lon4,lat4,lon))then !north
678  if(lon .ge. intersect(lat4,lon4,lat1,lon1,lat))then ! west
679  found = .true.
680  is=i; js=j
681  interp % i_lon (m,n,1) = i; interp % i_lon (m,n,2) = ip1
682  interp % j_lat (m,n,1) = j; interp % j_lat (m,n,2) = jp1
683  exit
684  endif
685  endif
686  endif
687  endif
688  enddo
689  step = step + 1
690  enddo
691  endif
692  if(.not.found) then
693  print *,'lon,lat=',lon*180./pi,lat*180./pi
694  print *,'npts=',npts
695  print *,'is,ie= ',istart,iend
696  print *,'js,je= ',jstart,jend
697  print *,'lon_in(is,js)=',lon_in(istart,jstart)*180./pi
698  print *,'lon_in(ie,js)=',lon_in(iend,jstart)*180./pi
699  print *,'lat_in(is,js)=',lat_in(istart,jstart)*180./pi
700  print *,'lat_in(ie,js)=',lat_in(iend,jstart)*180./pi
701  print *,'lon_in(is,je)=',lon_in(istart,jend)*180./pi
702  print *,'lon_in(ie,je)=',lon_in(iend,jend)*180./pi
703  print *,'lat_in(is,je)=',lat_in(istart,jend)*180./pi
704  print *,'lat_in(ie,je)=',lat_in(iend,jend)*180./pi
705 
706  call mpp_error(fatal, &
707  'find_neighbor: the destination point is not inside the source grid' )
708  endif
709  enddo
710  enddo
711 
712  end subroutine find_neighbor
713 
714  !#######################################################################
715  !
716  ! The function will return true if the point x,y is inside a polygon, or
717  ! NO if it is not. If the point is exactly on the edge of a polygon,
718  ! the function will return .true.
719  !
720  ! real polyx(:) : longitude coordinates of corners
721  ! real polyx(:) : latitude coordinates of corners
722  ! real x,y : point to be tested
723  ! ??? How to deal with truncation error.
724  !
725  function inside_polygon(polyx, polyy, x, y)
726  real, dimension(:), intent(in) :: polyx, polyy
727  real, intent(in) :: x, y
728  logical :: inside_polygon
729  integer :: i, j, nedges
730  real :: xx
731 
732  inside_polygon = .false.
733  nedges = size(polyx(:))
734  j = nedges
735  do i = 1, nedges
736  if( (polyy(i) < y .AND. polyy(j) >= y) .OR. (polyy(j) < y .AND. polyy(i) >= y) ) then
737  xx = polyx(i)+(y-polyy(i))/(polyy(j)-polyy(i))*(polyx(j)-polyx(i))
738  if( xx == x ) then
739  inside_polygon = .true.
740  return
741  else if( xx < x ) then
742  inside_polygon = .not. inside_polygon
743  endif
744  endif
745  j = i
746  enddo
747 
748  return
749 
750  end function inside_polygon
751 
752  !#######################################################################
753  ! this routine will search the source grid to fine the grid box that encloses
754  ! each destination grid.
755  subroutine find_neighbor_new( Interp, lon_in, lat_in, lon_out, lat_out, src_modulo, no_crash )
756  type(horiz_interp_type), intent(inout) :: Interp
757  real, intent(in), dimension(:,:) :: lon_in , lat_in
758  real, intent(in), dimension(:,:) :: lon_out, lat_out
759  logical, intent(in) :: src_modulo, no_crash
760  integer :: nlon_in, nlat_in, nlon_out, nlat_out
761  integer :: max_step, n, m, l, i, j, ip1, jp1, step
762  integer :: is, js, jstart, jend, istart, iend, npts
763  integer, allocatable, dimension(:) :: ilon, jlat
764  real :: lon_min, lon_max, lon, lat, tpi
765  logical :: found
766  real :: polyx(4), polyy(4)
767  real :: min_lon, min_lat, max_lon, max_lat
768 
769  integer, parameter :: step_div=8
770 
771  tpi = 2.0*pi
772  nlon_in = size(lon_in,1) ; nlat_in = size(lat_in,2)
773  nlon_out = size(lon_out,1); nlat_out = size(lon_out,2)
774 
775  lon_min = minval(lon_in);
776  lon_max = maxval(lon_in);
777 
778  max_step = min(nlon_in,nlat_in)/step_div ! can be adjusted if needed
779  allocate(ilon(step_div*max_step), jlat(step_div*max_step) )
780 
781  do n = 1, nlat_out
782  do m = 1, nlon_out
783  found = .false.
784  lon = lon_out(m,n)
785  lat = lat_out(m,n)
786 
787  if(src_modulo) then
788  if(lon .lt. lon_min) then
789  lon = lon + tpi
790  else if(lon .gt. lon_max) then
791  lon = lon - tpi
792  endif
793  else
794  if(lon .lt. lon_min .or. lon .gt. lon_max ) &
795  call mpp_error(fatal,'horiz_interp_bilinear_mod: ' //&
796  'when input grid is not modulo, output grid should locate inside input grid')
797  endif
798  !--- search for the surrounding four points locatioon.
799  if(m==1 .and. n==1) then
800  j_loop: do j = 1, nlat_in-1
801  do i = 1, nlon_in
802  ip1 = i+1
803  jp1 = j+1
804  if(i==nlon_in) then
805  if(src_modulo)then
806  ip1 = 1
807  else
808  cycle
809  endif
810  endif
811 
812  polyx(1) = lon_in(i, j); polyy(1) = lat_in(i,j)
813  polyx(2) = lon_in(ip1,j); polyy(2) = lat_in(ip1,j)
814  polyx(3) = lon_in(ip1,jp1); polyy(3) = lat_in(ip1,jp1)
815  polyx(4) = lon_in(i, jp1); polyy(4) = lat_in(i, jp1)
816  if(lon .lt. lon_min .or. lon .gt. lon_max) then
817  if(i .ne. nlon_in) then
818  cycle
819  else
820  if(lon .lt. lon_min) then
821  polyx(1) = polyx(1) -tpi; polyx(4) = polyx(4) - tpi
822  else if(lon .gt. lon_max) then
823  polyx(2) = polyx(2) +tpi; polyx(3) = polyx(3) + tpi
824  endif
825  endif
826  endif
827 
828  min_lon = minval(polyx)
829  max_lon = maxval(polyx)
830  min_lat = minval(polyy)
831  max_lat = maxval(polyy)
832 ! if( lon .GE. min_lon .AND. lon .LE. max_lon .AND. &
833 ! lat .GE. min_lat .AND. lat .LE. max_lat ) then
834 ! print*, 'i =', i, 'j = ', j
835 ! print '(5f15.11)', lon, polyx
836 ! print '(5f15.11)', lat, polyy
837 ! endif
838 
839  if(inside_polygon(polyx, polyy, lon, lat)) then
840  found = .true.
841 ! print*, " found ", i, j
842  interp % i_lon (m,n,1) = i; interp % i_lon (m,n,2) = ip1
843  interp % j_lat (m,n,1) = j; interp % j_lat (m,n,2) = jp1
844  exit j_loop
845  endif
846  enddo
847  enddo j_loop
848  else
849  step = 0
850  do while ( .not. found .and. step .lt. max_step )
851  !--- take the adajcent point as the starting point
852  if(m == 1) then
853  is = interp % i_lon (m,n-1,1)
854  js = interp % j_lat (m,n-1,1)
855  else
856  is = interp % i_lon (m-1,n,1)
857  js = interp % j_lat (m-1,n,1)
858  endif
859  if(step==0) then
860  npts = 1
861  ilon(1) = is
862  jlat(1) = js
863  else
864  npts = 0
865  !--- bottom and top boundary
866  jstart = max(js-step,1)
867  jend = min(js+step,nlat_in)
868 
869  do l = -step, step
870  i = is+l
871  if(src_modulo)then
872  if( i < 1) then
873  i = i + nlon_in
874  else if (i > nlon_in) then
875  i = i - nlon_in
876  endif
877  if( i < 1 .or. i > nlon_in) call mpp_error(fatal, &
878  'horiz_interp_bilinear_mod: max_step is too big, decrease max_step' )
879  else
880  if( i < 1 .or. i > nlon_in) cycle
881  endif
882 
883  npts = npts + 1
884  ilon(npts) = i
885  jlat(npts) = jstart
886  npts = npts + 1
887  ilon(npts) = i
888  jlat(npts) = jend
889  enddo
890 
891  !--- right and left boundary -----------------------------------------------
892  istart = is - step
893  iend = is + step
894  if(src_modulo) then
895  if( istart < 1) istart = istart + nlon_in
896  if( iend > nlon_in) iend = iend - nlon_in
897  else
898  istart = max(istart,1)
899  iend = min(iend, nlon_in)
900  endif
901  do l = -step, step
902  j = js+l
903  if( j < 1 .or. j > nlat_in) cycle
904  npts = npts+1
905  ilon(npts) = istart
906  jlat(npts) = j
907  npts = npts+1
908  ilon(npts) = iend
909  jlat(npts) = j
910  end do
911  end if
912 
913  !--- find the surrouding points
914  do l = 1, npts
915  i = ilon(l)
916  j = jlat(l)
917  ip1 = i+1
918  if(ip1>nlon_in) then
919  if(src_modulo) then
920  ip1 = 1
921  else
922  cycle
923  endif
924  endif
925  jp1 = j+1
926  if(jp1>nlat_in) cycle
927  polyx(1) = lon_in(i, j); polyy(1) = lat_in(i,j)
928  polyx(2) = lon_in(ip1,j); polyy(2) = lat_in(ip1,j)
929  polyx(3) = lon_in(ip1,jp1); polyy(3) = lat_in(ip1,jp1)
930  polyx(4) = lon_in(i, jp1); polyy(4) = lat_in(i, jp1)
931  if(inside_polygon(polyx, polyy, lon, lat)) then
932  found = .true.
933  interp % i_lon (m,n,1) = i; interp % i_lon (m,n,2) = ip1
934  interp % j_lat (m,n,1) = j; interp % j_lat (m,n,2) = jp1
935  exit
936  endif
937  enddo
938  step = step + 1
939  enddo
940  endif
941  if(.not.found) then
942  if(no_crash) then
943  interp % i_lon (m,n,1:2) = dummy
944  interp % j_lat (m,n,1:2) = dummy
945  print*,'lon,lat=',lon,lat ! snz
946  else
947  call mpp_error(fatal, &
948  'horiz_interp_bilinear_mod: the destination point is not inside the source grid' )
949  endif
950  endif
951  enddo
952  enddo
953 
954  end subroutine find_neighbor_new
955 
956  !#######################################################################
957  function intersect(x1, y1, x2, y2, x)
958  real, intent(in) :: x1, y1, x2, y2, x
959  real :: intersect
960 
961  intersect = (y2-y1)*(x-x1)/(x2-x1) + y1
962 
963  return
964 
965  end function intersect
966 
967  !#######################################################################
968  ! <SUBROUTINE NAME="horiz_interp_bilinear">
969 
970  ! <OVERVIEW>
971  ! Subroutine for performing the horizontal interpolation between two grids.
972  ! </OVERVIEW>
973  ! <DESCRIPTION>
974  ! Subroutine for performing the horizontal interpolation between two grids.
975  ! horiz_interp_bilinear_new must be called before calling this routine.
976  ! </DESCRIPTION>
977  ! <TEMPLATE>
978  ! call horiz_interp_bilinear ( Interp, data_in, data_out, verbose, mask_in,mask_out, missing_value, missing_permit)
979  ! </TEMPLATE>
980  !
981  ! <IN NAME="Interp" TYPE="type(horiz_interp_type)">
982  ! Derived-type variable containing interpolation indices and weights.
983  ! Returned by a previous call to horiz_interp_bilinear_new.
984  ! </IN>
985  ! <IN NAME="data_in" TYPE="real, dimension(:,:)">
986  ! Input data on source grid.
987  ! </IN>
988  ! <IN NAME="verbose" TYPE="integer, optional">
989  ! flag for the amount of print output.
990  ! verbose = 0, no output; = 1, min,max,means; = 2, still more
991  ! </IN>
992  ! <IN NAME="mask_in" TYPE="real, dimension(:,:),optional">
993  ! Input mask, must be the same size as the input data. The real value of
994  ! mask_in must be in the range (0.,1.). Set mask_in=0.0 for data points
995  ! that should not be used or have missing data.
996  ! </IN>
997  ! <IN NAME="missing_value" TYPE="real, optional">
998  ! Use the missing_value to indicate missing data.
999  ! </IN>
1000 
1001  ! <IN NAME="missing_permit" TUPE="integer, optional">
1002  ! numbers of points allowed to miss for the bilinear interpolation. The value
1003  ! should be between 0 and 3.
1004  ! </IN>
1005 
1006  ! <OUT NAME="data_out" TYPE="real, dimension(:,:)">
1007  ! Output data on destination grid.
1008  ! </OUT>
1009  ! <OUT NAME="mask_out" TYPE="real, dimension(:,:),optional">
1010  ! Output mask that specifies whether data was computed.
1011  ! </OUT>
1012 
1013  subroutine horiz_interp_bilinear ( Interp, data_in, data_out, verbose, mask_in,mask_out, &
1014  missing_value, missing_permit, new_handle_missing )
1015  !-----------------------------------------------------------------------
1016  type(horiz_interp_type), intent(in) :: interp
1017  real, intent(in), dimension(:,:) :: data_in
1018  real, intent(out), dimension(:,:) :: data_out
1019  integer, intent(in), optional :: verbose
1020  real, intent(in), dimension(:,:), optional :: mask_in
1021  real, intent(out), dimension(:,:), optional :: mask_out
1022  real, intent(in), optional :: missing_value
1023  integer, intent(in), optional :: missing_permit
1024  logical, intent(in), optional :: new_handle_missing
1025  !-----------------------------------------------------------------------
1026  integer :: nlon_in, nlat_in, nlon_out, nlat_out, n, m, &
1027  is, ie, js, je, iverbose, max_missing, num_missing, &
1028  miss_in, miss_out, unit
1029  real :: dwtsum, wtsum, min_in, max_in, avg_in, &
1030  min_out, max_out, avg_out, wtw, wte, wts, wtn
1031  real :: mask(size(data_in,1), size(data_in,2) )
1032  logical :: set_to_missing, is_missing(4), new_handler
1033  real :: f1, f2, f3, f4, middle, w, s
1034 
1035  num_missing = 0
1036 
1037  nlon_in = interp%nlon_src; nlat_in = interp%nlat_src
1038  nlon_out = interp%nlon_dst; nlat_out = interp%nlat_dst
1039 
1040  if(present(mask_in)) then
1041  mask = mask_in
1042  else
1043  mask = 1.0
1044  endif
1045 
1046  if (present(verbose)) then
1047  iverbose = verbose
1048  else
1049  iverbose = 0
1050  endif
1051 
1052  if(present(missing_permit)) then
1053  max_missing = missing_permit
1054  else
1055  max_missing = 0
1056  endif
1057 
1058  if(present(new_handle_missing)) then
1059  new_handler = new_handle_missing
1060  else
1061  new_handler = .false.
1062  endif
1063 
1064  if(max_missing .gt. 3 .or. max_missing .lt. 0) call mpp_error(fatal, &
1065  'horiz_interp_bilinear_mod: missing_permit should be between 0 and 3')
1066 
1067  if (size(data_in,1) /= nlon_in .or. size(data_in,2) /= nlat_in) &
1068  call mpp_error(fatal,'horiz_interp_bilinear_mod: size of input array incorrect')
1069 
1070  if (size(data_out,1) /= nlon_out .or. size(data_out,2) /= nlat_out) &
1071  call mpp_error(fatal,'horiz_interp_bilinear_mod: size of output array incorrect')
1072 
1073  if(new_handler) then
1074  if( .not. present(missing_value) ) call mpp_error(fatal, &
1075  "horiz_interp_bilinear_mod: misisng_value must be present when new_handle_missing is .true.")
1076  if( present(mask_in) ) call mpp_error(fatal, &
1077  "horiz_interp_bilinear_mod: mask_in should not be present when new_handle_missing is .true.")
1078  do n = 1, nlat_out
1079  do m = 1, nlon_out
1080  is = interp % i_lon (m,n,1); ie = interp % i_lon (m,n,2)
1081  js = interp % j_lat (m,n,1); je = interp % j_lat (m,n,2)
1082  wtw = interp % wti (m,n,1)
1083  wte = interp % wti (m,n,2)
1084  wts = interp % wtj (m,n,1)
1085  wtn = interp % wtj (m,n,2)
1086 
1087  is_missing = .false.
1088  num_missing = 0
1089  set_to_missing = .false.
1090  if(data_in(is,js) == missing_value) then
1091  num_missing = num_missing+1
1092  is_missing(1) = .true.
1093  if(wtw .GE. 0.5 .AND. wts .GE. 0.5) set_to_missing = .true.
1094  endif
1095 
1096  if(data_in(ie,js) == missing_value) then
1097  num_missing = num_missing+1
1098  is_missing(2) = .true.
1099  if(wte .GE. 0.5 .AND. wts .GE. 0.5) set_to_missing = .true.
1100  endif
1101  if(data_in(ie,je) == missing_value) then
1102  num_missing = num_missing+1
1103  is_missing(3) = .true.
1104  if(wte .GE. 0.5 .AND. wtn .GE. 0.5) set_to_missing = .true.
1105  endif
1106  if(data_in(is,je) == missing_value) then
1107  num_missing = num_missing+1
1108  is_missing(4) = .true.
1109  if(wtw .GE. 0.5 .AND. wtn .GE. 0.5) set_to_missing = .true.
1110  endif
1111 
1112  if( num_missing == 4 .OR. set_to_missing ) then
1113  data_out(m,n) = missing_value
1114  if(present(mask_out)) mask_out(m,n) = 0.0
1115  cycle
1116  else if(num_missing == 0) then
1117  f1 = data_in(is,js)
1118  f2 = data_in(ie,js)
1119  f3 = data_in(ie,je)
1120  f4 = data_in(is,je)
1121  w = wtw
1122  s = wts
1123  else if(num_missing == 3) then !--- three missing value
1124  if(.not. is_missing(1) ) then
1125  data_out(m,n) = data_in(is,js)
1126  else if(.not. is_missing(2) ) then
1127  data_out(m,n) = data_in(ie,js)
1128  else if(.not. is_missing(3) ) then
1129  data_out(m,n) = data_in(ie,je)
1130  else if(.not. is_missing(4) ) then
1131  data_out(m,n) = data_in(is,je)
1132  endif
1133  if(present(mask_out) ) mask_out(m,n) = 1.0
1134  cycle
1135  else !--- one or two missing value
1136  if( num_missing == 1) then
1137  if( is_missing(1) .OR. is_missing(3) ) then
1138  middle = 0.5*(data_in(ie,js)+data_in(is,je))
1139  else
1140  middle = 0.5*(data_in(is,js)+data_in(ie,je))
1141  endif
1142  else ! num_missing = 2
1143  if( is_missing(1) .AND. is_missing(2) ) then
1144  middle = 0.5*(data_in(ie,je)+data_in(is,je))
1145  else if( is_missing(1) .AND. is_missing(3) ) then
1146  middle = 0.5*(data_in(ie,js)+data_in(is,je))
1147  else if( is_missing(1) .AND. is_missing(4) ) then
1148  middle = 0.5*(data_in(ie,js)+data_in(ie,je))
1149  else if( is_missing(2) .AND. is_missing(3) ) then
1150  middle = 0.5*(data_in(is,js)+data_in(is,je))
1151  else if( is_missing(2) .AND. is_missing(4) ) then
1152  middle = 0.5*(data_in(is,js)+data_in(ie,je))
1153  else if( is_missing(3) .AND. is_missing(4) ) then
1154  middle = 0.5*(data_in(is,js)+data_in(ie,js))
1155  endif
1156  endif
1157 
1158  if( wtw .GE. 0.5 .AND. wts .GE. 0.5 ) then ! zone 1
1159  w = 2.0*(wtw-0.5)
1160  s = 2.0*(wts-0.5)
1161  f1 = data_in(is,js)
1162  if(is_missing(2)) then
1163  f2 = f1
1164  else
1165  f2 = 0.5*(data_in(is,js)+data_in(ie,js))
1166  endif
1167  f3 = middle
1168  if(is_missing(4)) then
1169  f4 = f1
1170  else
1171  f4 = 0.5*(data_in(is,js)+data_in(is,je))
1172  endif
1173  else if( wte .GE. 0.5 .AND. wts .GE. 0.5 ) then ! zone 2
1174  w = 2.0*(1.0-wte)
1175  s = 2.0*(wts-0.5)
1176  f2 = data_in(ie,js)
1177  if(is_missing(1)) then
1178  f1 = f2
1179  else
1180  f1 = 0.5*(data_in(is,js)+data_in(ie,js))
1181  endif
1182  f4 = middle
1183  if(is_missing(3)) then
1184  f3 = f2
1185  else
1186  f3 = 0.5*(data_in(ie,js)+data_in(ie,je))
1187  endif
1188  else if( wte .GE. 0.5 .AND. wtn .GE. 0.5 ) then ! zone 3
1189  w = 2.0*(1.0-wte)
1190  s = 2.0*(1.0-wtn)
1191  f3 = data_in(ie,je)
1192  if(is_missing(2)) then
1193  f2 = f3
1194  else
1195  f2 = 0.5*(data_in(ie,js)+data_in(ie,je))
1196  endif
1197  f1 = middle
1198  if(is_missing(4)) then
1199  f4 = f3
1200  else
1201  f4 = 0.5*(data_in(ie,je)+data_in(is,je))
1202  endif
1203  else if( wtw .GE. 0.5 .AND. wtn .GE. 0.5 ) then ! zone 4
1204  w = 2.0*(wtw-0.5)
1205  s = 2.0*(1.0-wtn)
1206  f4 = data_in(is,je)
1207  if(is_missing(1)) then
1208  f1 = f4
1209  else
1210  f1 = 0.5*(data_in(is,js)+data_in(is,je))
1211  endif
1212  f2 = middle
1213  if(is_missing(3)) then
1214  f3 = f4
1215  else
1216  f3 = 0.5*(data_in(ie,je)+data_in(is,je))
1217  endif
1218  else
1219  call mpp_error(fatal, &
1220  "horiz_interp_bilinear_mod: the point should be in one of the four zone")
1221  endif
1222  endif
1223 
1224  data_out(m,n) = f3 + (f4-f3)*w + (f2-f3)*s + ((f1-f2)+(f3-f4))*w*s
1225  if(present(mask_out)) mask_out(m,n) = 1.0
1226  enddo
1227  enddo
1228  else
1229  do n = 1, nlat_out
1230  do m = 1, nlon_out
1231  is = interp % i_lon (m,n,1); ie = interp % i_lon (m,n,2)
1232  js = interp % j_lat (m,n,1); je = interp % j_lat (m,n,2)
1233  wtw = interp % wti (m,n,1)
1234  wte = interp % wti (m,n,2)
1235  wts = interp % wtj (m,n,1)
1236  wtn = interp % wtj (m,n,2)
1237 
1238  if(present(missing_value) ) then
1239  num_missing = 0
1240  if(data_in(is,js) == missing_value) then
1241  num_missing = num_missing+1
1242  mask(is,js) = 0.0
1243  endif
1244  if(data_in(ie,js) == missing_value) then
1245  num_missing = num_missing+1
1246  mask(ie,js) = 0.0
1247  endif
1248  if(data_in(ie,je) == missing_value) then
1249  num_missing = num_missing+1
1250  mask(ie,je) = 0.0
1251  endif
1252  if(data_in(is,je) == missing_value) then
1253  num_missing = num_missing+1
1254  mask(is,je) = 0.0
1255  endif
1256  endif
1257 
1258  dwtsum = data_in(is,js)*mask(is,js)*wtw*wts &
1259  + data_in(ie,js)*mask(ie,js)*wte*wts &
1260  + data_in(ie,je)*mask(ie,je)*wte*wtn &
1261  + data_in(is,je)*mask(is,je)*wtw*wtn
1262  wtsum = mask(is,js)*wtw*wts + mask(ie,js)*wte*wts &
1263  + mask(ie,je)*wte*wtn + mask(is,je)*wtw*wtn
1264 
1265  if(.not. present(mask_in) .and. .not. present(missing_value)) wtsum = 1.0
1266 
1267  if(num_missing .gt. max_missing ) then
1268  data_out(m,n) = missing_value
1269  if(present(mask_out)) mask_out(m,n) = 0.0
1270  else if(wtsum .lt. epsln) then
1271  if(present(missing_value)) then
1272  data_out(m,n) = missing_value
1273  else
1274  data_out(m,n) = 0.0
1275  endif
1276  if(present(mask_out)) mask_out(m,n) = 0.0
1277  else
1278  data_out(m,n) = dwtsum/wtsum
1279  if(present(mask_out)) mask_out(m,n) = wtsum
1280  endif
1281  enddo
1282  enddo
1283  endif
1284  !***********************************************************************
1285  ! compute statistics: minimum, maximum, and mean
1286  !-----------------------------------------------------------------------
1287  if (iverbose > 0) then
1288 
1289  ! compute statistics of input data
1290 
1291  call stats (data_in, min_in, max_in, avg_in, miss_in, missing_value, mask_in)
1292 
1293  ! compute statistics of output data
1294  call stats (data_out, min_out, max_out, avg_out, miss_out, missing_value, mask_out)
1295 
1296  !---- output statistics ----
1297  unit = stdout()
1298  write (unit,900)
1299  write (unit,901) min_in ,max_in, avg_in
1300  if (present(mask_in)) write (unit,903) miss_in
1301  write (unit,902) min_out,max_out,avg_out
1302  if (present(mask_out)) write (unit,903) miss_out
1303 
1304 900 format (/,1x,10('-'),' output from horiz_interp ',10('-'))
1305 901 format (' input: min=',f16.9,' max=',f16.9,' avg=',f22.15)
1306 902 format (' output: min=',f16.9,' max=',f16.9,' avg=',f22.15)
1307 903 format (' number of missing points = ',i6)
1308 
1309  endif
1310 
1311  return
1312 
1313  end subroutine horiz_interp_bilinear
1314  ! </SUBROUTINE>
1315 
1316  !#######################################################################
1317  ! <SUBROUTINE NAME="horiz_interp_bilinear_del">
1318 
1319  ! <OVERVIEW>
1320  ! Deallocates memory used by "horiz_interp_type" variables.
1321  ! Must be called before reinitializing with horiz_interp_bilinear_new.
1322  ! </OVERVIEW>
1323  ! <DESCRIPTION>
1324  ! Deallocates memory used by "horiz_interp_type" variables.
1325  ! Must be called before reinitializing with horiz_interp_bilinear_new.
1326  ! </DESCRIPTION>
1327  ! <TEMPLATE>
1328  ! call horiz_interp_bilinear_del ( Interp )
1329  ! </TEMPLATE>
1330 
1331  ! <INOUT NAME="Interp" TYPE="horiz_interp_type">
1332  ! A derived-type variable returned by previous call
1333  ! to horiz_interp_bilinear_new. The input variable must have
1334  ! allocated arrays. The returned variable will contain
1335  ! deallocated arrays.
1336  ! </INOUT>
1337 
1338  subroutine horiz_interp_bilinear_del( Interp )
1339 
1340  type(horiz_interp_type), intent(inout) :: interp
1341 
1342  if(associated(interp%wti)) deallocate(interp%wti)
1343  if(associated(interp%wtj)) deallocate(interp%wtj)
1344  if(associated(interp%i_lon)) deallocate(interp%i_lon)
1345  if(associated(interp%j_lat)) deallocate(interp%j_lat)
1346 
1347  end subroutine horiz_interp_bilinear_del
1348  ! </SUBROUTINE>
1349 
1350  !#######################################################################
1351 
1352  function indp (value, array)
1353  integer :: indp
1354  real, dimension(:), intent(in) :: array
1355  real, intent(in) :: value
1356  !
1357  !=======================================================================
1358  !
1359  ! indp = index of nearest data point within "array" corresponding to
1360  ! "value".
1361 
1362  ! inputs:
1363  ! value = arbitrary data...same units as elements in "array"
1364  ! array = array of data points (must be monotonically increasing)
1365 
1366  ! output:
1367  ! indp = index of nearest data point to "value"
1368  ! if "value" is outside the domain of "array" then indp = 1
1369  ! or "ia" depending on whether array(1) or array(ia) is
1370  ! closest to "value"
1371  !=======================================================================
1372  !
1373  integer i, ia, unit
1374  logical keep_going
1375  !
1376  ia = size(array(:))
1377  do i=2,ia
1378  if (array(i) .lt. array(i-1)) then
1379  unit = stdout()
1380  write (unit,*) &
1381  ' => Error: array must be monotonically increasing in "indp"' , &
1382  ' when searching for nearest element to value=',value
1383  write (unit,*) ' array(i) < array(i-1) for i=',i
1384  write (unit,*) ' array(i) for i=1..ia follows:'
1385  call abort()
1386  endif
1387  enddo
1388  if (value .lt. array(1) .or. value .gt. array(ia)) then
1389  if (value .lt. array(1)) indp = 1
1390  if (value .gt. array(ia)) indp = ia
1391  else
1392  i=1
1393  keep_going = .true.
1394  do while (i .le. ia .and. keep_going)
1395  i = i+1
1396  if (value .le. array(i)) then
1397  indp = i
1398  if (array(i)-value .gt. value-array(i-1)) indp = i-1
1399  keep_going = .false.
1400  endif
1401  enddo
1402  endif
1403  return
1404  end function indp
1405 
1406  !######################################################################
1407 
1408 end module horiz_interp_bilinear_mod
fms_mod
Definition: fms.F90:20
constants_mod
Definition: constants.F90:24
horiz_interp_bilinear_mod::find_neighbor_new
subroutine find_neighbor_new(Interp, lon_in, lat_in, lon_out, lat_out, src_modulo, no_crash)
Definition: horiz_interp_bilinear.F90:756
horiz_interp_bilinear_mod::epsln
real, parameter epsln
Definition: horiz_interp_bilinear.F90:56
horiz_interp_type_mod
Definition: horiz_interp_type.F90:19
horiz_interp_bilinear_mod::intersect
real function intersect(x1, y1, x2, y2, x)
Definition: horiz_interp_bilinear.F90:958
horiz_interp_bilinear_mod::find_neighbor
subroutine find_neighbor(Interp, lon_in, lat_in, lon_out, lat_out, src_modulo)
Definition: horiz_interp_bilinear.F90:476
horiz_interp_bilinear_mod::horiz_interp_bilinear
subroutine, public horiz_interp_bilinear(Interp, data_in, data_out, verbose, mask_in, mask_out, missing_value, missing_permit, new_handle_missing)
Definition: horiz_interp_bilinear.F90:1015
horiz_interp_type_mod::stats
subroutine, public stats(dat, low, high, avg, miss, missing_value, mask)
Definition: horiz_interp_type.F90:98
horiz_interp_bilinear_mod::horiz_interp_bilinear_del
subroutine, public horiz_interp_bilinear_del(Interp)
Definition: horiz_interp_bilinear.F90:1339
mpp_mod
Definition: mpp.F90:39
horiz_interp_bilinear_mod::module_is_initialized
logical module_is_initialized
Definition: horiz_interp_bilinear.F90:62
horiz_interp_bilinear_mod::inside_polygon
logical function inside_polygon(polyx, polyy, x, y)
Definition: horiz_interp_bilinear.F90:726
horiz_interp_bilinear_mod::indp
integer function indp(value, array)
Definition: horiz_interp_bilinear.F90:1353
horiz_interp_bilinear_mod::horiz_interp_bilinear_init
subroutine, public horiz_interp_bilinear_init
Definition: horiz_interp_bilinear.F90:76
horiz_interp_bilinear_mod
Definition: horiz_interp_bilinear.F90:19
horiz_interp_bilinear_mod::horiz_interp_bilinear_new
Definition: horiz_interp_bilinear.F90:50
horiz_interp_bilinear_mod::dummy
integer, parameter dummy
Definition: horiz_interp_bilinear.F90:57
horiz_interp_type_mod::horiz_interp_type
Definition: horiz_interp_type.F90:57
mpp_mod::mpp_error
Definition: mpp.F90:402
horiz_interp_bilinear_mod::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)
Definition: horiz_interp_bilinear.F90:287
size
************************************************************************GNU Lesser General Public License **This file is part of the GFDL Flexible Modeling System(FMS). ! *! *FMS is free software without even the implied warranty of MERCHANTABILITY or *FITNESS FOR A PARTICULAR PURPOSE See the GNU General Public License *for more details **You should have received a copy of the GNU Lesser General Public *License along with FMS If see< http:! ***********************************************************************subroutine READ_RECORD_CORE_(unit, field, nwords, data, start, axsiz) integer, intent(in) ::unit type(fieldtype), intent(in) ::field integer, intent(in) ::nwords MPP_TYPE_, intent(inout) ::data(nwords) integer, intent(in) ::start(:), axsiz(:) integer(SHORT_KIND) ::i2vals(nwords)!rab used in conjunction with transfer intrinsic to determine size of a variable integer(KIND=1) ::one_byte(8) integer ::word_sz!#ifdef __sgi integer(INT_KIND) ::ivals(nwords) real(FLOAT_KIND) ::rvals(nwords)!#else! integer ::ivals(nwords)! real ::rvals(nwords)!#endif real(DOUBLE_KIND) ::r8vals(nwords) pointer(ptr1, i2vals) pointer(ptr2, ivals) pointer(ptr3, rvals) pointer(ptr4, r8vals) if(mpp_io_stack_size< nwords) call mpp_io_set_stack_size(nwords) call mpp_error(FATAL, 'MPP_READ currently requires use_netCDF option') end subroutine READ_RECORD_CORE_ subroutine READ_RECORD_(unit, field, nwords, data, time_level, domain, position, tile_count, start_in, axsiz_in)!routine that is finally called by all mpp_read routines to perform the read!a non-netCDF record contains:! field ID! a set of 4 coordinates(is:ie, js:je) giving the data subdomain! a timelevel and a timestamp(=NULLTIME if field is static)! 3D real data(stored as 1D)!if you are using direct access I/O, the RECL argument to OPEN must be large enough for the above!in a global direct access file, record position on PE is given by %record.!Treatment of timestamp:! We assume that static fields have been passed without a timestamp.! Here that is converted into a timestamp of NULLTIME.! For non-netCDF fields, field is treated no differently, but is written! with a timestamp of NULLTIME. There is no check in the code to prevent! the user from repeatedly writing a static field. integer, intent(in) ::unit, nwords type(fieldtype), intent(in) ::field MPP_TYPE_, intent(inout) ::data(nwords) integer, intent(in), optional ::time_level type(domain2D), intent(in), optional ::domain integer, intent(in), optional ::position, tile_count integer, intent(in), optional ::start_in(:), axsiz_in(:) integer, dimension(size(field%axes(:))) ::start, axsiz integer ::tlevel !, subdomain(4) integer ::i, error, is, ie, js, je, isg, ieg, jsg, jeg type(domain2d), pointer ::io_domain=> tlevel if(PRESENT(start_in) .AND. PRESENT(axsiz_in)) then if(size(start(! the data domain and compute domain must refer to the subdomain being passed ! In this ! since that attempts to gather all data on PE size(field%axes(:)) axsiz(i)
max
#define max(a, b)
Definition: mosaic_util.h:33
horiz_interp_bilinear_mod::horiz_interp_bilinear_new_1d
subroutine horiz_interp_bilinear_new_1d(Interp, lon_in, lat_in, lon_out, lat_out, verbose, src_modulo)
Definition: horiz_interp_bilinear.F90:89
min
#define min(a, b)
Definition: mosaic_util.h:32
crtm_parameters::pi
real(fp), parameter, public pi
Definition: CRTM_Parameters.f90:86