FMS  2024.03
Flexible Modeling System
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_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
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
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_
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
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
stdout
integer function stdout()
This function returns the current standard fortran unit numbers for output.
Definition: mpp_util.inc:43