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