FMS 2025.01-dev
Flexible Modeling System
Loading...
Searching...
No Matches
horiz_interp_bicubic.inc
1!***********************************************************************
2!* GNU Lesser General Public License
3!*
4!* This file is part of the GFDL Flexible Modeling System (FMS).
5!*
6!* FMS is free software: you can redistribute it and/or modify it under
7!* the terms of the GNU Lesser General Public License as published by
8!* the Free Software Foundation, either version 3 of the License, or (at
9!* your option) any later version.
10!*
11!* FMS is distributed in the hope that it will be useful, but WITHOUT
12!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13!* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
14!* for more details.
15!*
16!* You should have received a copy of the GNU Lesser General Public
17!* License along with FMS. If not, see <http://www.gnu.org/licenses/>.
18!***********************************************************************
19!> @addtogroup horiz_interp_bicubic_mod
20!> @{
21
22 !> @brief Creates a new @ref horiz_interp_type
23 !!
24 !> Allocates space and initializes a derived-type variable
25 !! that contains pre-computed interpolation indices and weights.
26 subroutine horiz_interp_bicubic_new_1d_s_ ( Interp, lon_in, lat_in, lon_out, lat_out, &
27 verbose, src_modulo )
28
29 !-----------------------------------------------------------------------
30 type(horiz_interp_type), intent(inout) :: Interp !< A derived-type variable containing indices
31 !! and weights used for subsequent interpolations. To
32 !! reinitialize this variable for a different grid-to-grid
33 !! interpolation you must first use the
34 !! @ref HORIZ_INTERP_BICUBIC_NEW__del interface.
35 real(FMS_HI_KIND_), intent(in), dimension(:) :: lon_in !< Longitude (radians) for source data grid
36 real(FMS_HI_KIND_), intent(in), dimension(:) :: lat_in !< Latitude (radians) for source data grid
37 real(FMS_HI_KIND_), intent(in), dimension(:,:) :: lon_out !< Longitude (radians) for output data grid
38 real(FMS_HI_KIND_), intent(in), dimension(:,:) :: lat_out !< Latitude (radians) for output data grid
39 integer, intent(in), optional :: verbose !< flag for print output amount
40 logical, intent(in), optional :: src_modulo !< indicates if the boundary condition along
41 !! zonal boundary is cyclic or not. Zonal boundary condition
42 !!is cyclic when true
43 integer :: i, j, ip1, im1, jp1, jm1
44 logical :: src_is_modulo
45 integer :: nlon_in, nlat_in, nlon_out, nlat_out
46 integer :: jcl, jcu, icl, icu, jj
47 real(FMS_HI_KIND_) :: xz, yz
48 integer :: iunit
49 integer, parameter :: kindl = fms_hi_kind_ !< real size at compile time
50
51 if(present(verbose)) verbose_bicubic = verbose
52 src_is_modulo = .false.
53 if (present(src_modulo)) src_is_modulo = src_modulo
54
55 if(size(lon_out,1) /= size(lat_out,1) .or. size(lon_out,2) /= size(lat_out,2) ) &
56 call mpp_error(fatal,'horiz_interp_bilinear_mod: when using bilinear ' // &
57 'interplation, the output grids should be geographical grids')
58
59 !--- get the grid size
60 nlon_in = size(lon_in) ; nlat_in = size(lat_in)
61 nlon_out = size(lon_out,1); nlat_out = size(lat_out,2)
62 interp%nlon_src = nlon_in; interp%nlat_src = nlat_in
63 interp%nlon_dst = nlon_out; interp%nlat_dst = nlat_out
64! use wti(:,:,1) for x-derivative, wti(:,:,2) for y-derivative, wti(:,:,3) for xy-derivative
65 allocate ( interp%HI_KIND_TYPE_%wti (nlon_in, nlat_in, 3) )
66 allocate ( interp%HI_KIND_TYPE_%lon_in (nlon_in) )
67 allocate ( interp%HI_KIND_TYPE_%lat_in (nlat_in) )
68 allocate ( interp%HI_KIND_TYPE_%rat_x (nlon_out, nlat_out) )
69 allocate ( interp%HI_KIND_TYPE_%rat_y (nlon_out, nlat_out) )
70 allocate ( interp%i_lon (nlon_out, nlat_out, 2) )
71 allocate ( interp%j_lat (nlon_out, nlat_out, 2) )
72
73 interp%HI_KIND_TYPE_%lon_in = lon_in
74 interp%HI_KIND_TYPE_%lat_in = lat_in
75
76 if ( verbose_bicubic > 0 ) then
77 iunit = stdout()
78 write (iunit,'(/,"Initialising bicubic interpolation, interface horiz_interp_bicubic_new_1d_s")')
79 write (iunit,'(/," Longitude of coarse grid points (radian): xc(i) i=1, ",i4)') interp%nlon_src
80 write (iunit,'(1x,10f10.4)') (interp%HI_KIND_TYPE_%lon_in(jj),jj=1,interp%nlon_src)
81 write (iunit,'(/," Latitude of coarse grid points (radian): yc(j) j=1, ",i4)') interp%nlat_src
82 write (iunit,'(1x,10f10.4)') (interp%HI_KIND_TYPE_%lat_in(jj),jj=1,interp%nlat_src)
83 do i=1, interp%nlat_dst
84 write (iunit,*)
85 write (iunit,'(/," Longitude of fine grid points (radian): xf(i) i=1, ",i4)') interp%nlat_dst
86 write (iunit,'(1x,10f10.4)') (lon_out(jj,i),jj=1,interp%nlon_dst)
87 enddo
88 do i=1, interp%nlon_dst
89 write (iunit,*)
90 write (iunit,'(/," Latitude of fine grid points (radian): yf(j) j=1, ",i4)') interp%nlon_dst
91 write (iunit,'(1x,10f10.4)') (lat_out(i,jj),jj=1,interp%nlat_dst)
92 enddo
93 endif
94
95
96!---------------------------------------------------------------------------
97! Find the x-derivative. Use central differences and forward or
98! backward steps at the boundaries
99
100 do j=1,nlat_in
101 do i=1,nlon_in
102 ip1=min(i+1,nlon_in)
103 im1=max(i-1,1)
104 interp%HI_KIND_TYPE_%wti(i,j,1) = 1.0_kindl/(interp%HI_KIND_TYPE_%lon_in(ip1)-interp%HI_KIND_TYPE_%lon_in(im1))
105 enddo
106 enddo
107
108
109!---------------------------------------------------------------------------
110
111! Find the y-derivative. Use central differences and forward or
112! backward steps at the boundaries
113 do j=1,nlat_in
114 jp1=min(j+1,nlat_in)
115 jm1=max(j-1,1)
116 do i=1,nlon_in
117 interp%HI_KIND_TYPE_%wti(i,j,2) =1.0_kindl/(interp%HI_KIND_TYPE_%lat_in(jp1)-interp%HI_KIND_TYPE_%lat_in(jm1))
118 enddo
119 enddo
120
121!---------------------------------------------------------------------------
122
123! Find the xy-derivative. Use central differences and forward or
124! backward steps at the boundaries
125 do j=1,nlat_in
126 jp1=min(j+1,nlat_in)
127 jm1=max(j-1,1)
128 do i=1,nlon_in
129 ip1=min(i+1,nlon_in)
130 im1=max(i-1,1)
131 interp%HI_KIND_TYPE_%wti(i,j,3) = 1.0_kindl / &
132 ((interp%HI_KIND_TYPE_%lon_in(ip1)-interp%HI_KIND_TYPE_%lon_in(im1)) * &
133 (interp%HI_KIND_TYPE_%lat_in(jp1)-interp%HI_KIND_TYPE_%lat_in(jm1)))
134 enddo
135 enddo
136!---------------------------------------------------------------------------
137! Now for each point at the dest-grid find the boundary points of
138! the source grid
139 do j=1, nlat_out
140 do i=1,nlon_out
141 yz = lat_out(i,j)
142 xz = lon_out(i,j)
143
144 jcl = 0
145 jcu = 0
146 if( yz .le. interp%HI_KIND_TYPE_%lat_in(1) ) then
147 jcl = 1
148 jcu = 1
149 else if( yz .ge. interp%HI_KIND_TYPE_%lat_in(nlat_in) ) then
150 jcl = nlat_in
151 jcu = nlat_in
152 else
153 jcl = indl(interp%HI_KIND_TYPE_%lat_in, yz)
154 jcu = indu(interp%HI_KIND_TYPE_%lat_in, yz)
155 endif
156
157 icl = 0
158 icu = 0
159 !--- cyclic condition, do we need to use do while
160 if( xz .gt. interp%HI_KIND_TYPE_%lon_in(nlon_in) ) xz = xz - real(tpi,fms_hi_kind_)
161 if( xz .le. interp%HI_KIND_TYPE_%lon_in(1) ) xz = xz + real(tpi,fms_hi_kind_)
162 if( xz .ge. interp%HI_KIND_TYPE_%lon_in(nlon_in) ) then
163 icl = nlon_in
164 icu = 1
165 interp%HI_KIND_TYPE_%rat_x(i,j) = (xz - interp%HI_KIND_TYPE_%lon_in(icl))/(interp%HI_KIND_TYPE_%lon_in(icu)&
166 & - interp%HI_KIND_TYPE_%lon_in(icl) + real(tpi,fms_hi_kind_))
167 else
168 icl = indl(interp%HI_KIND_TYPE_%lon_in, xz)
169 icu = indu(interp%HI_KIND_TYPE_%lon_in, xz)
170 interp%HI_KIND_TYPE_%rat_x(i,j) = (xz - interp%HI_KIND_TYPE_%lon_in(icl))/(interp%HI_KIND_TYPE_%lon_in(icu)&
171 & - interp%HI_KIND_TYPE_%lon_in(icl))
172 endif
173 interp%j_lat(i,j,1) = jcl
174 interp%j_lat(i,j,2) = jcu
175 interp%i_lon(i,j,1) = icl
176 interp%i_lon(i,j,2) = icu
177 if(jcl == jcu) then
178 interp%HI_KIND_TYPE_%rat_y(i,j) = 0.0_kindl
179 else
180 interp%HI_KIND_TYPE_%rat_y(i,j) = (yz-interp%HI_KIND_TYPE_%lat_in(jcl))/(interp%HI_KIND_TYPE_%lat_in(jcu)&
181 & - interp%HI_KIND_TYPE_%lat_in(jcl))
182 endif
183! if(yz.gt.Interp%HI_KIND_TYPE_%lat_in(jcu)) call mpp_error(FATAL, ' HORIZ_INTERP_BICUBIC_NEW_1D_S_:
184! yf < ycl, no valid boundary point')
185! if(yz.lt.Interp%HI_KIND_TYPE_%lat_in(jcl)) call mpp_error(FATAL, ' HORIZ_INTERP_BICUBIC_NEW_1D_S_:
186! yf > ycu, no valid boundary point')
187! if(xz.gt.Interp%HI_KIND_TYPE_%lon_in(icu)) call mpp_error(FATAL, ' HORIZ_INTERP_BICUBIC_NEW_1D_S_:
188! xf < xcl, no valid boundary point')
189! if(xz.lt.Interp%HI_KIND_TYPE_%lon_in(icl)) call mpp_error(FATAL, ' HORIZ_INTERP_BICUBIC_NEW_1D_S_:
190! xf > xcu, no valid boundary point')
191 enddo
192 enddo
193 interp% HI_KIND_TYPE_ % is_allocated = .true.
194 interp%interp_method = bicubic
195 end subroutine horiz_interp_bicubic_new_1d_s_
196
197 !> @brief Creates a new @ref horiz_interp_type
198 !!
199 !> Allocates space and initializes a derived-type variable
200 !! that contains pre-computed interpolation indices and weights.
201 subroutine horiz_interp_bicubic_new_1d_ ( Interp, lon_in, lat_in, lon_out, lat_out, &
202 verbose, src_modulo )
203
204 !-----------------------------------------------------------------------
205 type(horiz_interp_type), intent(inout) :: Interp
206 real(FMS_HI_KIND_), intent(in), dimension(:) :: lon_in , lat_in
207 real(FMS_HI_KIND_), intent(in), dimension(:) :: lon_out, lat_out
208 integer, intent(in), optional :: verbose
209 logical, intent(in), optional :: src_modulo
210 integer :: i, j, ip1, im1, jp1, jm1
211 logical :: src_is_modulo
212 integer :: nlon_in, nlat_in, nlon_out, nlat_out
213 integer :: jcl, jcu, icl, icu, jj
214 real(FMS_HI_KIND_) :: xz, yz
215 integer :: iunit
216 integer, parameter :: kindl = fms_hi_kind_ !< real size at compile time
217
218 if(present(verbose)) verbose_bicubic = verbose
219 src_is_modulo = .false.
220 if (present(src_modulo)) src_is_modulo = src_modulo
221
222 !--- get the grid size
223 nlon_in = size(lon_in) ; nlat_in = size(lat_in)
224 nlon_out = size(lon_out); nlat_out = size(lat_out)
225 interp%nlon_src = nlon_in; interp%nlat_src = nlat_in
226 interp%nlon_dst = nlon_out; interp%nlat_dst = nlat_out
227 allocate ( interp%HI_KIND_TYPE_%wti (nlon_in, nlat_in, 3) )
228 allocate ( interp%HI_KIND_TYPE_%lon_in (nlon_in) )
229 allocate ( interp%HI_KIND_TYPE_%lat_in (nlat_in) )
230 allocate ( interp%HI_KIND_TYPE_%rat_x (nlon_out, nlat_out) )
231 allocate ( interp%HI_KIND_TYPE_%rat_y (nlon_out, nlat_out) )
232 allocate ( interp%i_lon (nlon_out, nlat_out, 2) )
233 allocate ( interp%j_lat (nlon_out, nlat_out, 2) )
234
235 interp%HI_KIND_TYPE_%lon_in = lon_in
236 interp%HI_KIND_TYPE_%lat_in = lat_in
237
238 if ( verbose_bicubic > 0 ) then
239 iunit = stdout()
240 write (iunit,'(/,"Initialising bicubic interpolation, interface HORIZ_INTERP_BICUBIC_NEW_1D_")')
241 write (iunit,'(/," Longitude of coarse grid points (radian): xc(i) i=1, ",i4)') interp%nlon_src
242 write (iunit,'(1x,10f10.4)') (interp%HI_KIND_TYPE_%lon_in(jj),jj=1,interp%nlon_src)
243 write (iunit,'(/," Latitude of coarse grid points (radian): yc(j) j=1, ",i4)') interp%nlat_src
244 write (iunit,'(1x,10f10.4)') (interp%HI_KIND_TYPE_%lat_in(jj),jj=1,interp%nlat_src)
245 write (iunit,*)
246 write (iunit,'(/," Longitude of fine grid points (radian): xf(i) i=1, ",i4)') interp%nlat_dst
247 write (iunit,'(1x,10f10.4)') (lon_out(jj),jj=1,interp%nlon_dst)
248 write (iunit,'(/," Latitude of fine grid points (radian): yf(j) j=1, ",i4)') interp%nlon_dst
249 write (iunit,'(1x,10f10.4)') (lat_out(jj),jj=1,interp%nlat_dst)
250 endif
251
252
253!---------------------------------------------------------------------------
254! Find the x-derivative. Use central differences and forward or
255! backward steps at the boundaries
256
257 do j=1,nlat_in
258 do i=1,nlon_in
259 ip1=min(i+1,nlon_in)
260 im1=max(i-1,1)
261 interp%HI_KIND_TYPE_%wti(i,j,1) = 1.0_kindl /(lon_in(ip1)-lon_in(im1))
262 enddo
263 enddo
264
265
266!---------------------------------------------------------------------------
267
268! Find the y-derivative. Use central differences and forward or
269! backward steps at the boundaries
270 do j=1,nlat_in
271 jp1=min(j+1,nlat_in)
272 jm1=max(j-1,1)
273 do i=1,nlon_in
274 interp%HI_KIND_TYPE_%wti(i,j,2) = 1.0_kindl /(lat_in(jp1)-lat_in(jm1))
275 enddo
276 enddo
277
278!---------------------------------------------------------------------------
279
280! Find the xy-derivative. Use central differences and forward or
281! backward steps at the boundaries
282 do j=1,nlat_in
283 jp1=min(j+1,nlat_in)
284 jm1=max(j-1,1)
285 do i=1,nlon_in
286 ip1=min(i+1,nlon_in)
287 im1=max(i-1,1)
288 interp%HI_KIND_TYPE_%wti(i,j,3) = 1.0_kindl /((lon_in(ip1)-lon_in(im1))*(lat_in(jp1)-lat_in(jm1)))
289 enddo
290 enddo
291!---------------------------------------------------------------------------
292! Now for each point at the dest-grid find the boundary points of
293! the source grid
294 do j=1, nlat_out
295 yz = lat_out(j)
296 jcl = 0
297 jcu = 0
298 if( yz .le. lat_in(1) ) then
299 jcl = 1
300 jcu = 1
301 else if( yz .ge. lat_in(nlat_in) ) then
302 jcl = nlat_in
303 jcu = nlat_in
304 else
305 jcl = indl(lat_in, yz)
306 jcu = indu(lat_in, yz)
307 endif
308 do i=1,nlon_out
309 xz = lon_out(i)
310 icl = 0
311 icu = 0
312 !--- cyclic condition, do we need to use do while
313 if( xz .gt. lon_in(nlon_in) ) xz = xz - real(tpi,fms_hi_kind_)
314 if( xz .le. lon_in(1) ) xz = xz + real(tpi, fms_hi_kind_)
315 if( xz .ge. lon_in(nlon_in) ) then
316 icl = nlon_in
317 icu = 1
318 interp%HI_KIND_TYPE_%rat_x(i,j) = (xz - interp%HI_KIND_TYPE_%lon_in(icl))/(interp%HI_KIND_TYPE_%lon_in(icu)&
319 & - interp%HI_KIND_TYPE_%lon_in(icl) + real(tpi,fms_hi_kind_))
320 else
321 icl = indl(lon_in, xz)
322 icu = indu(lon_in, xz)
323 interp%HI_KIND_TYPE_%rat_x(i,j) = (xz - interp%HI_KIND_TYPE_%lon_in(icl))/(interp%HI_KIND_TYPE_%lon_in(icu)&
324 & - interp%HI_KIND_TYPE_%lon_in(icl))
325 endif
326 icl = indl(lon_in, xz)
327 icu = indu(lon_in, xz)
328 interp%j_lat(i,j,1) = jcl
329 interp%j_lat(i,j,2) = jcu
330 interp%i_lon(i,j,1) = icl
331 interp%i_lon(i,j,2) = icu
332 if(jcl == jcu) then
333 interp%HI_KIND_TYPE_%rat_y(i,j) = 0.0_kindl
334 else
335 interp%HI_KIND_TYPE_%rat_y(i,j) = (yz- interp%HI_KIND_TYPE_%lat_in(jcl))/(interp%HI_KIND_TYPE_%lat_in(jcu)&
336 & - interp%HI_KIND_TYPE_%lat_in(jcl))
337 endif
338! if(yz.gt.lat_in(jcu)) call mpp_error(FATAL, ' HORIZ_INTERP_BICUBIC_NEW_1D_: yf <
339! ycl, no valid boundary point')
340! if(yz.lt.lat_in(jcl)) call mpp_error(FATAL, ' HORIZ_INTERP_BICUBIC_NEW_1D_: yf >
341! ycu, no valid boundary point')
342! if(xz.gt.lon_in(icu)) call mpp_error(FATAL, ' HORIZ_INTERP_BICUBIC_NEW_1D_: xf <
343! xcl, no valid boundary point')
344! if(xz.lt.lon_in(icl)) call mpp_error(FATAL, ' HORIZ_INTERP_BICUBIC_NEW_1D_: xf >
345! xcu, no valid boundary point')
346 enddo
347 enddo
348 interp% HI_KIND_TYPE_ % is_allocated = .true.
349 interp%interp_method = bicubic
350
351 end subroutine horiz_interp_bicubic_new_1d_
352
353 !> @brief Perform bicubic horizontal interpolation
354 subroutine horiz_interp_bicubic_( Interp, data_in, data_out, verbose, mask_in, mask_out, missing_value, &
355 & missing_permit)
356 type (horiz_interp_type), intent(in) :: Interp
357 real(FMS_HI_KIND_), intent(in), dimension(:,:) :: data_in
358 real(FMS_HI_KIND_), intent(out), dimension(:,:) :: data_out
359 integer, intent(in), optional :: verbose
360 real(FMS_HI_KIND_), intent(in), dimension(:,:), optional :: mask_in
361 real(FMS_HI_KIND_), intent(out), dimension(:,:), optional :: mask_out
362 real(FMS_HI_KIND_), intent(in), optional :: missing_value
363 integer, intent(in), optional :: missing_permit
364 real(FMS_HI_KIND_) :: yz, ycu, ycl
365 real(FMS_HI_KIND_) :: xz, xcu, xcl
366 real(FMS_HI_KIND_) :: val, val1, val2
367 real(FMS_HI_KIND_), dimension(4) :: y, y1, y2, y12
368 integer :: icl, icu, jcl, jcu
369 integer :: iclp1, icup1, jclp1, jcup1
370 integer :: iclm1, icum1, jclm1, jcum1
371 integer :: i,j
372 integer, parameter :: kindl = fms_hi_kind_ !< set kind size at compile time
373
374 if ( present(verbose) ) verbose_bicubic = verbose
375
376 do j=1, interp%nlat_dst
377 do i=1, interp%nlon_dst
378 yz = interp%HI_KIND_TYPE_%rat_y(i,j)
379 xz = interp%HI_KIND_TYPE_%rat_x(i,j)
380 jcl = interp%j_lat(i,j,1)
381 jcu = interp%j_lat(i,j,2)
382 icl = interp%i_lon(i,j,1)
383 icu = interp%i_lon(i,j,2)
384 if( icl > icu ) then
385 iclp1 = icu
386 icum1 = icl
387 xcl = interp%HI_KIND_TYPE_%lon_in(icl)
388 xcu = interp%HI_KIND_TYPE_%lon_in(icu)+real(tpi, fms_hi_kind_)
389 else
390 iclp1 = min(icl+1,interp%nlon_src)
391 icum1 = max(icu-1,1)
392 xcl = interp%HI_KIND_TYPE_%lon_in(icl)
393 xcu = interp%HI_KIND_TYPE_%lon_in(icu)
394 endif
395 iclm1 = max(icl-1,1)
396 icup1 = min(icu+1,interp%nlon_src)
397 jclp1 = min(jcl+1,interp%nlat_src)
398 jclm1 = max(jcl-1,1)
399 jcup1 = min(jcu+1,interp%nlat_src)
400 jcum1 = max(jcu-1,1)
401 ycl = interp%HI_KIND_TYPE_%lat_in(jcl)
402 ycu = interp%HI_KIND_TYPE_%lat_in(jcu)
403! xcl = Interp%HI_KIND_TYPE_%lon_in(icl)
404! xcu = Interp%HI_KIND_TYPE_%lon_in(icu)
405 y(1) = data_in(icl,jcl)
406 y(2) = data_in(icu,jcl)
407 y(3) = data_in(icu,jcu)
408 y(4) = data_in(icl,jcu)
409 y1(1) = ( data_in(iclp1,jcl) - data_in(iclm1,jcl) ) * interp%HI_KIND_TYPE_%wti(icl,jcl,1)
410 y1(2) = ( data_in(icup1,jcl) - data_in(icum1,jcl) ) * interp%HI_KIND_TYPE_%wti(icu,jcl,1)
411 y1(3) = ( data_in(icup1,jcu) - data_in(icum1,jcu) ) * interp%HI_KIND_TYPE_%wti(icu,jcu,1)
412 y1(4) = ( data_in(iclp1,jcu) - data_in(iclm1,jcu) ) * interp%HI_KIND_TYPE_%wti(icl,jcu,1)
413 y2(1) = ( data_in(icl,jclp1) - data_in(icl,jclm1) ) * interp%HI_KIND_TYPE_%wti(icl,jcl,2)
414 y2(2) = ( data_in(icu,jclp1) - data_in(icu,jclm1) ) * interp%HI_KIND_TYPE_%wti(icu,jcl,2)
415 y2(3) = ( data_in(icu,jcup1) - data_in(icu,jcum1) ) * interp%HI_KIND_TYPE_%wti(icu,jcu,2)
416 y2(4) = ( data_in(icl,jcup1) - data_in(icl,jcum1) ) * interp%HI_KIND_TYPE_%wti(icl,jcu,2)
417 y12(1)= ( data_in(iclp1,jclp1) + data_in(iclm1,jclm1) - data_in(iclm1,jclp1) &
418 - data_in(iclp1,jclm1) ) * interp%HI_KIND_TYPE_%wti(icl,jcl,3)
419 y12(2)= ( data_in(icup1,jclp1) + data_in(icum1,jclm1) - data_in(icum1,jclp1) &
420 - data_in(icup1,jclm1) ) * interp%HI_KIND_TYPE_%wti(icu,jcl,3)
421 y12(3)= ( data_in(icup1,jcup1) + data_in(icum1,jcum1) - data_in(icum1,jcup1) &
422 - data_in(icup1,jcum1) ) * interp%HI_KIND_TYPE_%wti(icu,jcu,3)
423 y12(4)= ( data_in(iclp1,jcup1) + data_in(iclm1,jcum1) - data_in(iclm1,jcup1) &
424 - data_in(iclp1,jcum1) ) * interp%HI_KIND_TYPE_%wti(icl,jcu,3)
425
426 call bcuint(y,y1,y2,y12,xcl,xcu,ycl,ycu,xz,yz,val,val1,val2)
427 data_out(i,j) = val
428 if(present(mask_out)) mask_out(i,j) = 1.0_kindl
429!! dff_x(i,j) = val1
430!! dff_y(i,j) = val2
431 enddo
432 enddo
433 return
434 end subroutine horiz_interp_bicubic_
435
436!---------------------------------------------------------------------------
437
438 subroutine bcuint_(y,y1,y2,y12,x1l,x1u,x2l,x2u,t,u,ansy,ansy1,ansy2)
439 real(FMS_HI_KIND_) ansy,ansy1,ansy2,x1l,x1u,x2l,x2u,y(4),y1(4),y12(4),y2(4)
440! uses BCUCOF_
441 integer i
442 integer, parameter :: kindl = fms_hi_kind_ !< compiled kind size
443 real(FMS_HI_KIND_) t,u,c(4,4)
444 call bcucof(y,y1,y2,y12,x1u-x1l,x2u-x2l,c)
445 ansy=0.0_kindl
446 ansy2=0.0_kindl
447 ansy1=0.0_kindl
448 do i=4,1,-1
449 ansy=t*ansy+((c(i,4)*u+c(i,3))*u+c(i,2))*u+c(i,1)
450! ansy2=t*ansy2+(3.*c(i,4)*u+2.*c(i,3))*u+c(i,2)
451! ansy1=u*ansy1+(3.*c(4,i)*t+2.*c(3,i))*t+c(2,i)
452 enddo
453! ansy1=ansy1/(x1u-x1l) ! could be used for accuracy checks
454! ansy2=ansy2/(x2u-x2l) ! could be used for accuracy checks
455 return
456! (c) copr. 1986-92 numerical recipes software -3#(-)f.
457 end subroutine bcuint_
458!---------------------------------------------------------------------------
459
460 subroutine bcucof_(y,y1,y2,y12,d1,d2,c)
461 real(FMS_HI_KIND_) d1,d2,c(4,4),y(4),y1(4),y12(4),y2(4)
462 integer i,j,k,l
463 real(FMS_HI_KIND_) d1d2,xx,cl(16),x(16)
464 integer, parameter :: kindl = fms_hi_kind_!< compiled kind type
465 !! n*0.0 represents n consecutive 0.0's
466 real(FMS_HI_KIND_), save, dimension(16,16) :: wt !< weights use
467 data wt/1.0_kindl, 0.0_kindl, -3.0_kindl, 2.0_kindl, 4*0.0_kindl, -3.0_kindl, 0.0_kindl, 9.0_kindl, -6.0_kindl, &
468 2.0_kindl, 0.0_kindl, -6.0_kindl, 4.0_kindl, 8*0.0_kindl, 3.0_kindl, 0.0_kindl, -9.0_kindl, 6.0_kindl, &
469 -2.0_kindl, 0.0_kindl, 6.0_kindl, -4.0_kindl, 10*0.0_kindl, 9.0_kindl, -6.0_kindl, 2*0.0_kindl, &
470 -6.0_kindl, 4.0_kindl, 2*0.0_kindl, 3.0_kindl, -2.0_kindl, 6*0.0_kindl, -9.0_kindl, 6.0_kindl, &
471 2*0.0_kindl, 6.0_kindl, -4.0_kindl, 4*0.0_kindl, 1.0_kindl, 0.0_kindl, -3.0_kindl, 2.0_kindl,-2.0_kindl,&
472 0.0_kindl, 6.0_kindl, -4.0_kindl, 1.0_kindl, 0.0_kindl, -3.0_kindl, 2.0_kindl, 8*0.0_kindl, -1.0_kindl, &
473 0.0_kindl, 3.0_kindl, -2.0_kindl, 1.0_kindl, 0.0_kindl, -3.0_kindl, 2.0_kindl, 10*0.0_kindl, -3.0_kindl,&
474 2.0_kindl, 2*0.0_kindl, 3.0_kindl, -2.0_kindl, 6*0.0_kindl, 3.0_kindl, -2.0_kindl, 2*0.0_kindl, &
475 -6.0_kindl, 4.0_kindl, 2*0.0_kindl, 3.0_kindl, -2.0_kindl, 0.0_kindl, 1.0_kindl, -2.0_kindl, 1.0_kindl, &
476 5*0.0_kindl, -3.0_kindl, 6.0_kindl, -3.0_kindl, 0.0_kindl, 2.0_kindl, -4.0_kindl, 2.0_kindl,9*0.0_kindl,&
477 3.0_kindl, -6.0_kindl, 3.0_kindl, 0.0_kindl, -2.0_kindl, 4.0_kindl, -2.0_kindl, 10*0.0_kindl,-3.0_kindl,&
478 3.0_kindl, 2*0.0_kindl, 2.0_kindl, -2.0_kindl, 2*0.0_kindl, -1.0_kindl, 1.0_kindl,6*0.0_kindl,3.0_kindl,&
479 -3.0_kindl, 2*0.0_kindl, -2.0_kindl, 2.0_kindl, 5*0.0_kindl, 1.0_kindl, -2.0_kindl, 1.0_kindl,0.0_kindl,&
480 -2.0_kindl, 4.0_kindl, -2.0_kindl, 0.0_kindl, 1.0_kindl, -2.0_kindl, 1.0_kindl, 9*0.0_kindl, -1.0_kindl,&
481 2.0_kindl, -1.0_kindl, 0.0_kindl, 1.0_kindl, -2.0_kindl, 1.0_kindl, 10*0.0_kindl, 1.0_kindl, -1.0_kindl,&
482 2*0.0_kindl, -1.0_kindl, 1.0_kindl, 6*0.0_kindl, -1.0_kindl, 1.0_kindl, 2*0.0_kindl, 2.0_kindl, &
483 -2.0_kindl, 2*0.0_kindl, -1.0_kindl, 1.0_kindl/
484
485
486
487 d1d2=d1*d2
488 do i=1,4
489 x(i)=y(i)
490 x(i+4)=y1(i)*d1
491 x(i+8)=y2(i)*d2
492 x(i+12)=y12(i)*d1d2
493 enddo
494 do i=1,16
495 xx=0.0_kindl
496 do k=1,16
497 xx=xx+wt(i,k)*x(k)
498 enddo
499 cl(i)=xx
500 enddo
501 l=0
502 do i=1,4
503 do j=1,4
504 l=l+1
505 c(i,j)=cl(l)
506 enddo
507 enddo
508 return
509! (c) copr. 1986-92 numerical recipes software -3#(-)f.
510 end subroutine bcucof_
511
512!-----------------------------------------------------------------------
513
514!! TODO These routines are redundant, we can find the lower neighbor and add 1
515!> find the lower neighbour of xf in field xc, return is the index
516 function indl_(xc, xf)
517 real(fms_hi_kind_), intent(in) :: xc(1:)
518 real(fms_hi_kind_), intent(in) :: xf
519 integer :: indl_
520 integer :: ii
521 indl_ = 1
522 do ii=1, size(xc)
523 if(xc(ii).gt.xf) return
524 indl_ = ii
525 enddo
526 call mpp_error(fatal,'Error in INDL_')
527 return
528 end function indl_
529
530!-----------------------------------------------------------------------
531
532!> find the upper neighbour of xf in field xc, return is the index
533 function indu_(xc, xf)
534 real(fms_hi_kind_), intent(in) :: xc(1:)
535 real(fms_hi_kind_), intent(in) :: xf
536 integer :: indu_
537 integer :: ii
538 do ii=1, size(xc)
539 indu_ = ii
540 if(xc(ii).gt.xf) return
541 enddo
542 call mpp_error(fatal,'Error in INDU_')
543 return
544 end function indu_
545
546!-----------------------------------------------------------------------
547
548 subroutine fill_xy_(fi, ics, ice, jcs, jce, mask, maxpass)
549 integer, intent(in) :: ics,ice,jcs,jce
550 real(fms_hi_kind_), intent(inout) :: fi(ics:ice,jcs:jce)
551 real(fms_hi_kind_), intent(in), optional :: mask(ics:ice,jcs:jce)
552 integer, intent(in) :: maxpass
553 real(fms_hi_kind_) :: work_old(ics:ice,jcs:jce)
554 real(fms_hi_kind_) :: work_new(ics:ice,jcs:jce)
555 logical :: ready
556 integer, parameter :: kindl = fms_hi_kind_
557 real(fms_hi_kind_), parameter :: blank = real(-1.e30, fms_hi_kind_)
558 real(fms_hi_kind_) :: tavr
559 integer :: ipass
560 integer :: inl, inr, jnl, jnu, i, j, is, js, iavr
561
562
563 ready = .false.
564
565 work_new(:,:) = fi(:,:)
566 work_old(:,:) = work_new(:,:)
567 ipass = 0
568 if ( present(mask) ) then
569 do while (.not.ready)
570 ipass = ipass+1
571 ready = .true.
572 do j=jcs, jce
573 do i=ics, ice
574 if (work_old(i,j).le.blank) then
575 tavr=0.0_kindl
576 iavr=0
577 inl = max(i-1,ics)
578 inr = min(i+1,ice)
579 jnl = max(j-1,jcs)
580 jnu = min(j+1,jce)
581 do js=jnl,jnu
582 do is=inl,inr
583 if (work_old(is,js) .ne. blank .and. mask(is,js).ne.0.0_kindl) then
584 tavr = tavr + work_old(is,js)
585 iavr = iavr+1
586 endif
587 enddo
588 enddo
589 if (iavr.gt.0) then
590 if (iavr.eq.1) then
591! spreading is not allowed if the only valid neighbor is a corner point
592! otherwise an ill posed cellular automaton is established leading to
593! a spreading of constant values in diagonal direction
594! if all corner points are blanked the valid neighbor must be a direct one
595! and spreading is allowed
596 if (work_old(inl,jnu).eq.blank.and.&
597 work_old(inr,jnu).eq.blank.and.&
598 work_old(inr,jnl).eq.blank.and.&
599 work_old(inl,jnl).eq.blank) then
600 work_new(i,j)=tavr/real(iavr,fms_hi_kind_)
601 ready = .false.
602 endif
603 else
604 work_new(i,j)=tavr/real(iavr,fms_hi_kind_)
605 ready = .false.
606 endif
607 endif
608 endif
609 enddo ! j
610 enddo ! i
611! save changes made during this pass to work_old
612 work_old(:,:)=work_new(:,:)
613 if(ipass.eq.maxpass) ready=.true.
614 enddo !while (.not.ready)
615 fi(:,:) = work_new(:,:)
616 else
617 do while (.not.ready)
618 ipass = ipass+1
619 ready = .true.
620 do j=jcs, jce
621 do i=ics, ice
622 if (work_old(i,j).le.blank) then
623 tavr=0.0_kindl
624 iavr=0
625 inl = max(i-1,ics)
626 inr = min(i+1,ice)
627 jnl = max(j-1,jcs)
628 jnu = min(j+1,jce)
629 do is=inl,inr
630 do js=jnl,jnu
631 if (work_old(is,js).gt.blank) then
632 tavr = tavr + work_old(is,js)
633 iavr = iavr+1
634 endif
635 enddo
636 enddo
637 if (iavr.gt.0) then
638 if (iavr.eq.1) then
639! spreading is not allowed if the only valid neighbor is a corner point
640! otherwise an ill posed cellular automaton is established leading to
641! a spreading of constant values in diagonal direction
642! if all corner points are blanked the valid neighbor must be a direct one
643! and spreading is allowed
644 if (work_old(inl,jnu).le.blank.and. &
645 work_old(inr,jnu).le.blank.and. &
646 work_old(inr,jnl).le.blank.and. &
647 work_old(inl,jnl).le.blank) then
648 work_new(i,j)=tavr/real(iavr,fms_hi_kind_)
649 ready = .false.
650 endif
651 else
652 work_new(i,j)=tavr/real(iavr,fms_hi_kind_)
653 ready = .false.
654 endif
655 endif
656 endif
657 enddo ! j
658 enddo ! i
659! save changes made during this pass to work_old
660 work_old(:,:)=work_new(:,:)
661 if(ipass.eq.maxpass) ready=.true.
662 enddo !while (.not.ready)
663 fi(:,:) = work_new(:,:)
664 endif
665 return
666 end subroutine fill_xy_
667!> @}
horiz_interp_bicubic_new_1d_
subroutine horiz_interp_bicubic_new_1d_(interp, lon_in, lat_in, lon_out, lat_out, verbose, src_modulo)
Creates a new horiz_interp_type.
Definition horiz_interp_bicubic.inc:203
horiz_interp_bicubic_
subroutine horiz_interp_bicubic_(interp, data_in, data_out, verbose, mask_in, mask_out, missing_value, missing_permit)
Perform bicubic horizontal interpolation.
Definition horiz_interp_bicubic.inc:356
indl_
integer function indl_(xc, xf)
find the lower neighbour of xf in field xc, return is the index
Definition horiz_interp_bicubic.inc:517
horiz_interp_bicubic_new_1d_s_
subroutine horiz_interp_bicubic_new_1d_s_(interp, lon_in, lat_in, lon_out, lat_out, verbose, src_modulo)
Creates a new horiz_interp_type.
Definition horiz_interp_bicubic.inc:28
indu_
integer function indu_(xc, xf)
find the upper neighbour of xf in field xc, return is the index
Definition horiz_interp_bicubic.inc:534