| blob | f998b823f73ab3f3d5ec662a49c893fc2f67006b |
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 )
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
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
40 warns = 0
41 if(present(verbose)) warns = verbose
42 src_is_modulo = .true.
43 if (present(src_modulo)) src_is_modulo = src_modulo
45 decreasing_lat = .false.
46 if (lat_in(1) > lat_in(2)) decreasing_lat = .true.
48 decreasing_lon = .false.
49 if (lon_in(1) > lon_in(2)) decreasing_lon = .true.
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 !-----------------------------------------------------------------------
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 !-----------------------------------------------------------------------
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
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')
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
85 do n = 1, nlat_out
86 do m = 1, nlon_out
87 lon = lon_out(m,n)
88 lat = lat_out(m,n)
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
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)
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
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)
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
159 wtn = 1.0_kindl - wts
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
168 enddo
169 enddo
171 iunit = stdout()
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
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
197 return
199 end subroutine HORIZ_INTERP_BILINEAR_NEW_1D_
201 !#######################################################################
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 )
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
235 integer, parameter :: kindl=FMS_HI_KIND_
237 tpi = 2.0_kindl * real(pi, FMS_HI_KIND_)
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
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')
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')
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
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))
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
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 !***************************************************************************
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
401 Interp% HI_KIND_TYPE_ % is_allocated = .true.
402 Interp% interp_method = BILINEAR
403 end subroutine
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
421 integer, parameter :: kindl=FMS_HI_KIND_
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)
427 lon_min = minval(lon_in);
428 lon_max = maxval(lon_in);
430 max_step = max(nlon_in,nlat_in) ! can be adjusted if needed
431 allocate(ilon(5*max_step), jlat(5*max_step) )
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)
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)
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
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)
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
529 npts = npts + 1
530 ilon(npts) = i
531 jlat(npts) = jstart
532 enddo
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
555 !--- top boundary
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
571 npts = npts + 1
572 ilon(npts) = i
573 jlat(npts) = jend
574 enddo
577 end if
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)
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
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_)
641 call mpp_error(FATAL, &
642 'FIND_NEIGHBOR_: the destination point is not inside the source grid' )
643 endif
644 enddo
645 enddo
647 end subroutine
649 !#######################################################################
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
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
679 return
681 end function
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
700 integer, parameter :: step_div=8, kindl = FMS_HI_KIND_
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)
706 lon_min = minval(lon_in);
707 lon_max = maxval(lon_in);
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) )
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)
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
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
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
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)
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
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
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
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
885 end subroutine
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_
892 INTERSECT_ = (y2-y1)*(x-x1)/(x2-x1) + y1
894 return
896 end function INTERSECT_
898 !#######################################################################
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_
933 num_missing = 0
935 nlon_in = Interp%nlon_src; nlat_in = Interp%nlat_src
936 nlon_out = Interp%nlon_dst; nlat_out = Interp%nlat_dst
938 if(present(mask_in)) then
939 mask = mask_in
940 else
941 mask = 1.0_kindl
942 endif
944 if (present(verbose)) then
945 iverbose = verbose
946 else
947 iverbose = 0
948 endif
950 if(present(missing_permit)) then
951 max_missing = missing_permit
952 else
953 max_missing = 0
954 endif
956 if(present(new_handle_missing)) then
957 new_handler = new_handle_missing
958 else
959 new_handler = .false.
960 endif
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')
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')
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')
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)
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
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
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
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
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)
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
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
1163 if(.not. present(mask_in) .and. .not. present(missing_value)) wtsum = 1.0_kindl
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
1187 ! compute statistics of input data
1189 call stats (data_in, min_in, max_in, avg_in, miss_in, missing_value, mask_in)
1191 ! compute statistics of output data
1192 call stats (data_out, min_out, max_out, avg_out, miss_out, missing_value, mask_out)
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
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)
1207 endif
1209 return
1211 end subroutine
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
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
1237 if (.not. open_file(weight_fileobj, weight_filename, "read" )) &
1238 call mpp_error(FATAL, "Error opening the weight file:"//&
1239 &trim(weight_filename))
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
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)
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))
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/))
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
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
1281 deallocate(var)
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/))
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)
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_
1301 !> @}