horiz_interp/horiz_interp_bicubic.F90

   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 !> @defgroup horiz_interp_bicubic_mod horiz_interp_bicubic_mod
  20 !> @ingroup horiz_interp
  21 !> @brief Delivers methods for bicubic interpolation from a coarse regular grid
  22 !! on a fine regular grid
  23 !!
  24 !> This module delivers methods for bicubic interpolation from a
  25 !! coarse regular grid on a fine regular grid.
  26 !! Subroutines
  27 !!
  28 !! - @ref bcuint
  29 !! - @ref bcucof
  30 !!
  31 !! are methods taken from
  32 !!
  33 !!       W. H. Press, S. A. Teukolski, W. T. Vetterling and B. P. Flannery,
  34 !!       Numerical Recipies in FORTRAN, The Art of Scientific Computing.
  35 !!       Cambridge University Press, 1992
  36 !!
  37 !! written by
  38 !!       martin.schmidt@io-warnemuende.de (2004)
  39 !! revised by
  40 !!       martin.schmidt@io-warnemuende.de (2004)
  41 !!
  42 !! Version 1.0.0.2005-07-06
  43 !! The module is thought to interact with MOM-4.
  44 !! Alle benotigten Felder werden extern von MOM verwaltet, da sie
  45 !! nicht fur alle interpolierten Daten die gleiche Dimension haben mussen.
  46
  47 !> @file
  48 !> @brief File for @ref horiz_interp_bicubic_mod
  49
  50 module horiz_interp_bicubic_mod
  51
  52   use mpp_mod,               only: mpp_error, FATAL, stdout, mpp_pe, mpp_root_pe
  53   use fms_mod,               only: write_version_number
  54   use horiz_interp_type_mod, only: horiz_interp_type
  55   use constants_mod,         only: PI
  56
  57
  58  implicit none
  59
  60    private
  61
  62    public  :: horiz_interp_bicubic, horiz_interp_bicubic_new, horiz_interp_bicubic_del, fill_xy
  63    public  :: horiz_interp_bicubic_init
  64
  65   !> Creates a new @ref horiz_interp_type for bicubic interpolation.
  66   !> @ingroup horiz_interp_bicubic_mod
  67   interface horiz_interp_bicubic_new
  68     module procedure horiz_interp_bicubic_new_1d
  69     module procedure horiz_interp_bicubic_new_1d_s
  70   end interface
  71
  72 !> @addtogroup horiz_interp_bicubic_mod
  73 !> @{
  74
  75 ! Include variable "version" to be written to log file.
  76 #include<file_version.h>
  77    logical            :: module_is_initialized = .FALSE.
  78    integer            :: verbose_bicubic = 0
  79
  80 !     Grid variables
  81 !     xc, yc : co-ordinates of the coarse grid
  82 !     xf, yf : co-ordinates of the fine grid
  83 !     fc     : variable to be interpolated at the coarse grid
  84 !     dfc_x  : x-derivative of fc at the coarse grid
  85 !     dfc_y  : y-derivative of fc at the coarse grid
  86 !     dfc_xy : x-y-derivative of fc at the coarse grid
  87 !     ff     : variable to be interpolated at the fine grid
  88 !     dff_x  : x-derivative of fc at the fine grid
  89 !     dff_y  : y-derivative of fc at the fine grid
  90 !     dff_xy : x-y-derivative of fc at the fine grid
  91
  92
  93    real               :: tpi
  94
  95    interface fill_xy
  96       module procedure fill_xy
  97    end interface
  98
  99
 100    contains
 101
 102   !> @brief Initializes module and writes version number to logfile.out
 103   subroutine horiz_interp_bicubic_init
 104
 105      if(module_is_initialized) return
 106      call write_version_number("HORIZ_INTERP_BICUBIC_MOD", version)
 107      module_is_initialized = .true.
 108      tpi = 2.0*PI
 109
 110   end subroutine horiz_interp_bicubic_init
 111
 112   !#######################################################################
 113
 114   !> @brief Creates a new @ref horiz_interp_type
 115   !!
 116   !> Allocates space and initializes a derived-type variable
 117   !! that contains pre-computed interpolation indices and weights.
 118   subroutine horiz_interp_bicubic_new_1d_s ( Interp, lon_in, lat_in, lon_out, lat_out, &
 119        verbose, src_modulo )
 120
 121     !-----------------------------------------------------------------------
 122     type(horiz_interp_type), intent(inout) :: Interp !< A derived-type variable containing indices
 123                                        !! and weights used for subsequent interpolations. To
 124                                        !! reinitialize this variable for a different grid-to-grid
 125                                        !! interpolation you must first use the
 126                                        !! @ref horiz_interp_bicubic_del interface.
 127     real, intent(in),  dimension(:)        :: lon_in !< Longitude (radians) for source data grid
 128     real, intent(in),  dimension(:)        :: lat_in !< Latitude (radians) for source data grid
 129     real, intent(in),  dimension(:,:)      :: lon_out !< Longitude (radians) for output data grid
 130     real, intent(in),  dimension(:,:)      :: lat_out !< Latitude (radians) for output data grid
 131     integer, intent(in),          optional :: verbose !< flag for print output amount
 132     logical, intent(in),          optional :: src_modulo !< indicates if the boundary condition along
 133                                        !! zonal boundary is cyclic or not. Zonal boundary condition
 134                                        !!is cyclic when true
 135     integer                                :: i, j, ip1, im1, jp1, jm1
 136     logical                                :: src_is_modulo
 137     integer                                :: nlon_in, nlat_in, nlon_out, nlat_out
 138     integer                                :: jcl, jcu, icl, icu, jj
 139     real                                   :: xz, yz
 140     integer                                :: unit
 141
 142     if(present(verbose)) verbose_bicubic = verbose
 143     src_is_modulo = .false.
 144     if (present(src_modulo)) src_is_modulo = src_modulo
 145
 146     if(size(lon_out,1) /= size(lat_out,1) .or. size(lon_out,2) /= size(lat_out,2) ) &
 147          call mpp_error(FATAL,'horiz_interp_bilinear_mod: when using bilinear ' // &
 148          'interplation, the output grids should be geographical grids')
 149
 150     !--- get the grid size
 151     nlon_in  = size(lon_in)   ; nlat_in  = size(lat_in)
 152     nlon_out = size(lon_out,1); nlat_out = size(lat_out,2)
 153     Interp%nlon_src = nlon_in;  Interp%nlat_src = nlat_in
 154     Interp%nlon_dst = nlon_out; Interp%nlat_dst = nlat_out
 155 !   use wti(:,:,1) for x-derivative, wti(:,:,2) for y-derivative, wti(:,:,3) for xy-derivative
 156     allocate ( Interp%wti    (nlon_in, nlat_in, 3) )
 157     allocate ( Interp%lon_in (nlon_in) )
 158     allocate ( Interp%lat_in (nlat_in) )
 159     allocate ( Interp%rat_x  (nlon_out, nlat_out) )
 160     allocate ( Interp%rat_y  (nlon_out, nlat_out) )
 161     allocate ( Interp%i_lon  (nlon_out, nlat_out, 2) )
 162     allocate ( Interp%j_lat  (nlon_out, nlat_out, 2) )
 163
 164     Interp%lon_in = lon_in
 165     Interp%lat_in = lat_in
 166
 167     if ( verbose_bicubic > 0 ) then
 168        unit = stdout()
 169        write (unit,'(/,"Initialising bicubic interpolation, interface horiz_interp_bicubic_new_1d_s")')
 170        write (unit,'(/," Longitude of coarse grid points (radian): xc(i) i=1, ",i4)') Interp%nlon_src
 171        write (unit,'(1x,10f10.4)') (Interp%lon_in(jj),jj=1,Interp%nlon_src)
 172        write (unit,'(/," Latitude of coarse grid points (radian):  yc(j) j=1, ",i4)') Interp%nlat_src
 173        write (unit,'(1x,10f10.4)') (Interp%lat_in(jj),jj=1,Interp%nlat_src)
 174        do i=1, Interp%nlat_dst
 175          write (unit,*)
 176          write (unit,'(/," Longitude of fine grid points (radian): xf(i) i=1, ",i4)') Interp%nlat_dst
 177          write (unit,'(1x,10f10.4)') (lon_out(jj,i),jj=1,Interp%nlon_dst)
 178        enddo
 179        do i=1, Interp%nlon_dst
 180          write (unit,*)
 181          write (unit,'(/," Latitude of fine grid points (radian):  yf(j) j=1, ",i4)') Interp%nlon_dst
 182          write (unit,'(1x,10f10.4)') (lat_out(i,jj),jj=1,Interp%nlat_dst)
 183        enddo
 184     endif
 185
 186
 187 !---------------------------------------------------------------------------
 188 !     Find the x-derivative. Use central differences and forward or
 189 !     backward steps at the boundaries
 190
 191     do j=1,nlat_in
 192       do i=1,nlon_in
 193         ip1=min(i+1,nlon_in)
 194         im1=max(i-1,1)
 195         Interp%wti(i,j,1) = 1./(Interp%lon_in(ip1)-Interp%lon_in(im1))
 196       enddo
 197     enddo
 198
 199
 200 !---------------------------------------------------------------------------
 201
 202 !     Find the y-derivative. Use central differences and forward or
 203 !     backward steps at the boundaries
 204       do j=1,nlat_in
 205         jp1=min(j+1,nlat_in)
 206         jm1=max(j-1,1)
 207         do i=1,nlon_in
 208           Interp%wti(i,j,2) = 1./(Interp%lat_in(jp1)-Interp%lat_in(jm1))
 209         enddo
 210       enddo
 211
 212 !---------------------------------------------------------------------------
 213
 214 !     Find the xy-derivative. Use central differences and forward or
 215 !     backward steps at the boundaries
 216       do j=1,nlat_in
 217         jp1=min(j+1,nlat_in)
 218         jm1=max(j-1,1)
 219         do i=1,nlon_in
 220           ip1=min(i+1,nlon_in)
 221           im1=max(i-1,1)
 222           Interp%wti(i,j,3) = 1./((Interp%lon_in(ip1)-Interp%lon_in(im1))*(Interp%lat_in(jp1)-Interp%lat_in(jm1)))
 223         enddo
 224       enddo
 225 !---------------------------------------------------------------------------
 226 !     Now for each point at the dest-grid find the boundary points of
 227 !     the source grid
 228       do j=1, nlat_out
 229         do i=1,nlon_out
 230           yz  = lat_out(i,j)
 231           xz  = lon_out(i,j)
 232
 233           jcl = 0
 234           jcu = 0
 235           if( yz .le. Interp%lat_in(1) ) then
 236              jcl = 1
 237              jcu = 1
 238           else if( yz .ge. Interp%lat_in(nlat_in) ) then
 239              jcl = nlat_in
 240              jcu = nlat_in
 241           else
 242              jcl = indl(Interp%lat_in, yz)
 243              jcu = indu(Interp%lat_in, yz)
 244           endif
 245
 246           icl = 0
 247           icu = 0
 248           !--- cyclic condition, do we need to use do while
 249           if( xz .gt. Interp%lon_in(nlon_in) ) xz = xz - tpi
 250           if( xz .le. Interp%lon_in(1) ) xz = xz + tpi
 251           if( xz .ge. Interp%lon_in(nlon_in) ) then
 252             icl = nlon_in
 253             icu = 1
 254             Interp%rat_x(i,j) = (xz - Interp%lon_in(icl))/(Interp%lon_in(icu) - Interp%lon_in(icl) + tpi)
 255           else
 256             icl = indl(Interp%lon_in, xz)
 257             icu = indu(Interp%lon_in, xz)
 258             Interp%rat_x(i,j) = (xz - Interp%lon_in(icl))/(Interp%lon_in(icu) - Interp%lon_in(icl))
 259           endif
 260           Interp%j_lat(i,j,1) = jcl
 261           Interp%j_lat(i,j,2) = jcu
 262           Interp%i_lon(i,j,1) = icl
 263           Interp%i_lon(i,j,2) = icu
 264           if(jcl == jcu) then
 265              Interp%rat_y(i,j) = 0.0
 266           else
 267              Interp%rat_y(i,j) = (yz - Interp%lat_in(jcl))/(Interp%lat_in(jcu) - Interp%lat_in(jcl))
 268           endif
 269 !          if(yz.gt.Interp%lat_in(jcu)) call mpp_error(FATAL, ' horiz_interp_bicubic_new_1d_s:
 270 !          yf < ycl, no valid boundary point')
 271 !          if(yz.lt.Interp%lat_in(jcl)) call mpp_error(FATAL, ' horiz_interp_bicubic_new_1d_s:
 272 !          yf > ycu, no valid boundary point')
 273 !          if(xz.gt.Interp%lon_in(icu)) call mpp_error(FATAL, ' horiz_interp_bicubic_new_1d_s:
 274 !          xf < xcl, no valid boundary point')
 275 !          if(xz.lt.Interp%lon_in(icl)) call mpp_error(FATAL, ' horiz_interp_bicubic_new_1d_s:
 276 !          xf > xcu, no valid boundary point')
 277         enddo
 278       enddo
 279   end subroutine horiz_interp_bicubic_new_1d_s
 280
 281   !> @brief Creates a new @ref horiz_interp_type
 282   !!
 283   !> Allocates space and initializes a derived-type variable
 284   !! that contains pre-computed interpolation indices and weights.
 285   subroutine horiz_interp_bicubic_new_1d ( Interp, lon_in, lat_in, lon_out, lat_out, &
 286        verbose, src_modulo )
 287
 288     !-----------------------------------------------------------------------
 289     type(horiz_interp_type), intent(inout) :: Interp
 290     real, intent(in),  dimension(:)        :: lon_in , lat_in
 291     real, intent(in),  dimension(:)        :: lon_out, lat_out
 292     integer, intent(in),          optional :: verbose
 293     logical, intent(in),          optional :: src_modulo
 294     integer                                :: i, j, ip1, im1, jp1, jm1
 295     logical                                :: src_is_modulo
 296     integer                                :: nlon_in, nlat_in, nlon_out, nlat_out
 297     integer                                :: jcl, jcu, icl, icu, jj
 298     real                                   :: xz, yz
 299     integer                                :: unit
 300
 301     if(present(verbose)) verbose_bicubic = verbose
 302     src_is_modulo = .false.
 303     if (present(src_modulo)) src_is_modulo = src_modulo
 304
 305     !--- get the grid size
 306     nlon_in  = size(lon_in) ; nlat_in  = size(lat_in)
 307     nlon_out = size(lon_out); nlat_out = size(lat_out)
 308     Interp%nlon_src = nlon_in;  Interp%nlat_src = nlat_in
 309     Interp%nlon_dst = nlon_out; Interp%nlat_dst = nlat_out
 310     allocate ( Interp%wti     (nlon_in, nlat_in, 3) )
 311     allocate ( Interp%lon_in  (nlon_in) )
 312     allocate ( Interp%lat_in  (nlat_in) )
 313     allocate ( Interp%rat_x   (nlon_out, nlat_out) )
 314     allocate ( Interp%rat_y   (nlon_out, nlat_out) )
 315     allocate ( Interp%i_lon   (nlon_out, nlat_out, 2) )
 316     allocate ( Interp%j_lat   (nlon_out, nlat_out, 2) )
 317
 318     Interp%lon_in = lon_in
 319     Interp%lat_in = lat_in
 320
 321     if ( verbose_bicubic > 0 ) then
 322        unit = stdout()
 323        write (unit,'(/,"Initialising bicubic interpolation, interface horiz_interp_bicubic_new_1d")')
 324        write (unit,'(/," Longitude of coarse grid points (radian): xc(i) i=1, ",i4)') Interp%nlon_src
 325        write (unit,'(1x,10f10.4)') (Interp%lon_in(jj),jj=1,Interp%nlon_src)
 326        write (unit,'(/," Latitude of coarse grid points (radian):  yc(j) j=1, ",i4)') Interp%nlat_src
 327        write (unit,'(1x,10f10.4)') (Interp%lat_in(jj),jj=1,Interp%nlat_src)
 328        write (unit,*)
 329        write (unit,'(/," Longitude of fine grid points (radian): xf(i) i=1, ",i4)') Interp%nlat_dst
 330        write (unit,'(1x,10f10.4)') (lon_out(jj),jj=1,Interp%nlon_dst)
 331        write (unit,'(/," Latitude of fine grid points (radian):  yf(j) j=1, ",i4)') Interp%nlon_dst
 332        write (unit,'(1x,10f10.4)') (lat_out(jj),jj=1,Interp%nlat_dst)
 333     endif
 334
 335
 336 !---------------------------------------------------------------------------
 337 !     Find the x-derivative. Use central differences and forward or
 338 !     backward steps at the boundaries
 339
 340     do j=1,nlat_in
 341       do i=1,nlon_in
 342         ip1=min(i+1,nlon_in)
 343         im1=max(i-1,1)
 344         Interp%wti(i,j,1) = 1./(lon_in(ip1)-lon_in(im1))
 345       enddo
 346     enddo
 347
 348
 349 !---------------------------------------------------------------------------
 350
 351 !     Find the y-derivative. Use central differences and forward or
 352 !     backward steps at the boundaries
 353       do j=1,nlat_in
 354         jp1=min(j+1,nlat_in)
 355         jm1=max(j-1,1)
 356         do i=1,nlon_in
 357           Interp%wti(i,j,2) = 1./(lat_in(jp1)-lat_in(jm1))
 358         enddo
 359       enddo
 360
 361 !---------------------------------------------------------------------------
 362
 363 !     Find the xy-derivative. Use central differences and forward or
 364 !     backward steps at the boundaries
 365       do j=1,nlat_in
 366         jp1=min(j+1,nlat_in)
 367         jm1=max(j-1,1)
 368         do i=1,nlon_in
 369           ip1=min(i+1,nlon_in)
 370           im1=max(i-1,1)
 371           Interp%wti(i,j,3) = 1./((lon_in(ip1)-lon_in(im1))*(lat_in(jp1)-lat_in(jm1)))
 372         enddo
 373       enddo
 374 !---------------------------------------------------------------------------
 375 !     Now for each point at the dest-grid find the boundary points of
 376 !     the source grid
 377       do j=1, nlat_out
 378         yz  = lat_out(j)
 379         jcl = 0
 380         jcu = 0
 381         if( yz .le. lat_in(1) ) then
 382            jcl = 1
 383            jcu = 1
 384         else if( yz .ge. lat_in(nlat_in) ) then
 385            jcl = nlat_in
 386            jcu = nlat_in
 387         else
 388            jcl = indl(lat_in, yz)
 389            jcu = indu(lat_in, yz)
 390         endif
 391         do i=1,nlon_out
 392           xz = lon_out(i)
 393           icl = 0
 394           icu = 0
 395          !--- cyclic condition, do we need to use do while
 396           if( xz .gt. lon_in(nlon_in) ) xz = xz - tpi
 397           if( xz .le. lon_in(1) ) xz = xz + tpi
 398           if( xz .ge. lon_in(nlon_in) ) then
 399             icl = nlon_in
 400             icu = 1
 401             Interp%rat_x(i,j) = (xz - Interp%lon_in(icl))/(Interp%lon_in(icu) - Interp%lon_in(icl) + tpi)
 402           else
 403             icl = indl(lon_in, xz)
 404             icu = indu(lon_in, xz)
 405             Interp%rat_x(i,j) = (xz - Interp%lon_in(icl))/(Interp%lon_in(icu) - Interp%lon_in(icl))
 406           endif
 407           icl = indl(lon_in, xz)
 408           icu = indu(lon_in, xz)
 409           Interp%j_lat(i,j,1) = jcl
 410           Interp%j_lat(i,j,2) = jcu
 411           Interp%i_lon(i,j,1) = icl
 412           Interp%i_lon(i,j,2) = icu
 413           if(jcl == jcu) then
 414              Interp%rat_y(i,j) = 0.0
 415           else
 416              Interp%rat_y(i,j) = (yz - Interp%lat_in(jcl))/(Interp%lat_in(jcu) - Interp%lat_in(jcl))
 417           endif
 418 !          if(yz.gt.lat_in(jcu)) call mpp_error(FATAL, ' horiz_interp_bicubic_new_1d: yf <
 419 !          ycl, no valid boundary point')
 420 !          if(yz.lt.lat_in(jcl)) call mpp_error(FATAL, ' horiz_interp_bicubic_new_1d: yf >
 421 !          ycu, no valid boundary point')
 422 !          if(xz.gt.lon_in(icu)) call mpp_error(FATAL, ' horiz_interp_bicubic_new_1d: xf <
 423 !          xcl, no valid boundary point')
 424 !          if(xz.lt.lon_in(icl)) call mpp_error(FATAL, ' horiz_interp_bicubic_new_1d: xf >
 425 !          xcu, no valid boundary point')
 426         enddo
 427       enddo
 428
 429   end subroutine horiz_interp_bicubic_new_1d
 430
 431   !> @brief Perform bicubic horizontal interpolation
 432   subroutine horiz_interp_bicubic( Interp, data_in, data_out, verbose, mask_in, mask_out, missing_value, &
 433                                  &  missing_permit)
 434     type (horiz_interp_type), intent(in)        :: Interp
 435     real, intent(in),  dimension(:,:)           :: data_in
 436     real, intent(out), dimension(:,:)           :: data_out
 437     integer, intent(in),               optional :: verbose
 438     real, intent(in),  dimension(:,:), optional :: mask_in
 439     real, intent(out), dimension(:,:), optional :: mask_out
 440     real, intent(in),                  optional :: missing_value
 441     integer, intent(in),               optional :: missing_permit
 442     real :: yz, ycu, ycl
 443     real :: xz, xcu, xcl
 444     real :: val, val1, val2
 445     real, dimension(4) :: y, y1, y2, y12
 446     integer :: icl, icu, jcl, jcu
 447     integer :: iclp1, icup1, jclp1, jcup1
 448     integer :: iclm1, icum1, jclm1, jcum1
 449     integer :: i,j
 450
 451     if ( present(verbose) ) verbose_bicubic = verbose
 452 !    fill_in = .false.
 453 !    if ( present(fill) ) fill_in = fill
 454 !   use dfc_x and dfc_y as workspace
 455 !    if ( fill_in ) call fill_xy(fc(ics:ice,jcs:jce), ics, ice, jcs, jce, maxpass=2)
 456 !    where ( data_in .le. missing ) data_in(:,:) = 0.
 457 !!
 458     do j=1, Interp%nlat_dst
 459       do i=1, Interp%nlon_dst
 460         yz  = Interp%rat_y(i,j)
 461         xz  = Interp%rat_x(i,j)
 462         jcl = Interp%j_lat(i,j,1)
 463         jcu = Interp%j_lat(i,j,2)
 464         icl = Interp%i_lon(i,j,1)
 465         icu = Interp%i_lon(i,j,2)
 466         if( icl > icu ) then
 467           iclp1 = icu
 468           icum1 = icl
 469           xcl = Interp%lon_in(icl)
 470           xcu = Interp%lon_in(icu)+tpi
 471         else
 472           iclp1 = min(icl+1,Interp%nlon_src)
 473           icum1 = max(icu-1,1)
 474           xcl = Interp%lon_in(icl)
 475           xcu = Interp%lon_in(icu)
 476         endif
 477         iclm1 = max(icl-1,1)
 478         icup1 = min(icu+1,Interp%nlon_src)
 479         jclp1 = min(jcl+1,Interp%nlat_src)
 480         jclm1 = max(jcl-1,1)
 481         jcup1 = min(jcu+1,Interp%nlat_src)
 482         jcum1 = max(jcu-1,1)
 483         ycl = Interp%lat_in(jcl)
 484         ycu = Interp%lat_in(jcu)
 485 !        xcl = Interp%lon_in(icl)
 486 !        xcu = Interp%lon_in(icu)
 487         y(1)  =  data_in(icl,jcl)
 488         y(2)  =  data_in(icu,jcl)
 489         y(3)  =  data_in(icu,jcu)
 490         y(4)  =  data_in(icl,jcu)
 491         y1(1) = ( data_in(iclp1,jcl) - data_in(iclm1,jcl) ) * Interp%wti(icl,jcl,1)
 492         y1(2) = ( data_in(icup1,jcl) - data_in(icum1,jcl) ) * Interp%wti(icu,jcl,1)
 493         y1(3) = ( data_in(icup1,jcu) - data_in(icum1,jcu) ) * Interp%wti(icu,jcu,1)
 494         y1(4) = ( data_in(iclp1,jcu) - data_in(iclm1,jcu) ) * Interp%wti(icl,jcu,1)
 495         y2(1) = ( data_in(icl,jclp1) - data_in(icl,jclm1) ) * Interp%wti(icl,jcl,2)
 496         y2(2) = ( data_in(icu,jclp1) - data_in(icu,jclm1) ) * Interp%wti(icu,jcl,2)
 497         y2(3) = ( data_in(icu,jcup1) - data_in(icu,jcum1) ) * Interp%wti(icu,jcu,2)
 498         y2(4) = ( data_in(icl,jcup1) - data_in(icl,jcum1) ) * Interp%wti(icl,jcu,2)
 499         y12(1)= ( data_in(iclp1,jclp1) + data_in(iclm1,jclm1) - data_in(iclm1,jclp1) &
 500                 - data_in(iclp1,jclm1) ) * Interp%wti(icl,jcl,3)
 501         y12(2)= ( data_in(icup1,jclp1) + data_in(icum1,jclm1) - data_in(icum1,jclp1) &
 502                 - data_in(icup1,jclm1) ) * Interp%wti(icu,jcl,3)
 503         y12(3)= ( data_in(icup1,jcup1) + data_in(icum1,jcum1) - data_in(icum1,jcup1) &
 504                 - data_in(icup1,jcum1) ) * Interp%wti(icu,jcu,3)
 505         y12(4)= ( data_in(iclp1,jcup1) + data_in(iclm1,jcum1) - data_in(iclm1,jcup1) &
 506                 - data_in(iclp1,jcum1) ) * Interp%wti(icl,jcu,3)
 507
 508         call bcuint(y,y1,y2,y12,xcl,xcu,ycl,ycu,xz,yz,val,val1,val2)
 509         data_out   (i,j) = val
 510         if(present(mask_out)) mask_out(i,j) = 1.
 511 !!        dff_x(i,j) = val1
 512 !!        dff_y(i,j) = val2
 513       enddo
 514     enddo
 515   return
 516   end subroutine horiz_interp_bicubic
 517
 518
 519 !---------------------------------------------------------------------------
 520
 521    subroutine bcuint(y,y1,y2,y12,x1l,x1u,x2l,x2u,t,u,ansy,ansy1,ansy2)
 522       real ansy,ansy1,ansy2,x1l,x1u,x2l,x2u,y(4),y1(4),y12(4),y2(4)
 523 !     uses bcucof
 524       integer i
 525       real t,u,c(4,4)
 526       call bcucof(y,y1,y2,y12,x1u-x1l,x2u-x2l,c)
 527       ansy=0.
 528       ansy2=0.
 529       ansy1=0.
 530       do i=4,1,-1
 531         ansy=t*ansy+((c(i,4)*u+c(i,3))*u+c(i,2))*u+c(i,1)
 532 !        ansy2=t*ansy2+(3.*c(i,4)*u+2.*c(i,3))*u+c(i,2)
 533 !        ansy1=u*ansy1+(3.*c(4,i)*t+2.*c(3,i))*t+c(2,i)
 534       enddo
 535 !      ansy1=ansy1/(x1u-x1l) ! could be used for accuracy checks
 536 !      ansy2=ansy2/(x2u-x2l) ! could be used for accuracy checks
 537       return
 538 !  (c) copr. 1986-92 numerical recipes software -3#(-)f.
 539    end subroutine bcuint
 540 !---------------------------------------------------------------------------
 541
 542    subroutine bcucof(y,y1,y2,y12,d1,d2,c)
 543       real d1,d2,c(4,4),y(4),y1(4),y12(4),y2(4)
 544       integer i,j,k,l
 545       real d1d2,xx,cl(16),wt(16,16),x(16)
 546       save wt
 547       data wt/1., 0., -3., 2., 4*0., -3., 0., 9., -6., 2., 0., -6., 4., 8*0., &
 548               3., 0., -9., 6., -2., 0., 6., -4., 10*0., 9., -6., 2*0., -6., &
 549               4., 2*0., 3., -2., 6*0., -9., 6., 2*0., 6., -4., 4*0., 1., 0., &
 550               -3., 2., -2., 0., 6., -4., 1., 0., -3., 2., 8*0., -1., 0., 3., &
 551               -2., 1., 0., -3., 2., 10*0., -3., 2., 2*0., 3., -2., 6*0., 3., &
 552               -2., 2*0., -6., 4., 2*0., 3., -2., 0., 1., -2., 1., 5*0., -3., &
 553               6., -3., 0., 2., -4., 2., 9*0., 3., -6., 3., 0., -2., 4., -2., &
 554               10*0., -3., 3., 2*0., 2., -2., 2*0., -1., 1., 6*0., 3., -3., &
 555               2*0., -2., 2., 5*0., 1., -2., 1., 0., -2., 4., -2., 0., 1., -2., &
 556               1., 9*0., -1., 2., -1., 0., 1., -2., 1., 10*0., 1., -1., 2*0., &
 557               -1., 1., 6*0., -1., 1., 2*0., 2., -2., 2*0., -1., 1./
 558
 559       d1d2=d1*d2
 560       do i=1,4
 561         x(i)=y(i)
 562         x(i+4)=y1(i)*d1
 563         x(i+8)=y2(i)*d2
 564         x(i+12)=y12(i)*d1d2
 565       enddo
 566       do i=1,16
 567         xx=0.
 568         do k=1,16
 569           xx=xx+wt(i,k)*x(k)
 570         enddo
 571         cl(i)=xx
 572       enddo
 573       l=0
 574       do i=1,4
 575         do j=1,4
 576           l=l+1
 577           c(i,j)=cl(l)
 578         enddo
 579       enddo
 580       return
 581 !  (c) copr. 1986-92 numerical recipes software -3#(-)f.
 582    end subroutine bcucof
 583
 584 !-----------------------------------------------------------------------
 585
 586 !> find the lower neighbour of xf in field xc, return is the index
 587     function indl(xc, xf)
 588     real, intent(in) :: xc(1:)
 589     real, intent(in) :: xf
 590     integer             :: indl
 591     integer             :: ii
 592        indl = 1
 593        do ii=1, size(xc)
 594          if(xc(ii).gt.xf) return
 595          indl = ii
 596        enddo
 597        call mpp_error(FATAL,'Error in indl')
 598     return
 599     end function indl
 600
 601 !-----------------------------------------------------------------------
 602
 603 !> find the upper neighbour of xf in field xc, return is the index
 604     function indu(xc, xf)
 605     real, intent(in) :: xc(1:)
 606     real, intent(in) :: xf
 607     integer             :: indu
 608     integer             :: ii
 609        do ii=1, size(xc)
 610          indu = ii
 611          if(xc(ii).gt.xf) return
 612        enddo
 613        call mpp_error(FATAL,'Error in indu')
 614     return
 615     end function indu
 616
 617 !-----------------------------------------------------------------------
 618
 619     subroutine fill_xy(fi, ics, ice, jcs, jce, mask, maxpass)
 620       integer, intent(in)        :: ics,ice,jcs,jce
 621       real, intent(inout)        :: fi(ics:ice,jcs:jce)
 622       real, intent(in), optional :: mask(ics:ice,jcs:jce)
 623       integer, intent(in)        :: maxpass
 624       real                       :: work_old(ics:ice,jcs:jce)
 625       real                       :: work_new(ics:ice,jcs:jce)
 626       logical :: ready
 627       real    :: blank = -1.e30
 628       real    :: tavr
 629       integer :: ipass = 0
 630       integer :: inl, inr, jnl, jnu, i, j, is, js,  iavr
 631
 632
 633       ready = .false.
 634
 635       work_new(:,:) = fi(:,:)
 636       work_old(:,:) = work_new(:,:)
 637       ipass = 0
 638       if ( present(mask) ) then
 639          do while (.not.ready)
 640            ipass = ipass+1
 641            ready = .true.
 642            do j=jcs, jce
 643              do i=ics, ice
 644                if (work_old(i,j).le.blank) then
 645                  tavr=0.
 646                  iavr=0
 647                  inl = max(i-1,ics)
 648                  inr = min(i+1,ice)
 649                  jnl = max(j-1,jcs)
 650                  jnu = min(j+1,jce)
 651                  do js=jnl,jnu
 652                    do is=inl,inr
 653                      if (work_old(is,js) .ne. blank .and. mask(is,js).ne.0.) then
 654                        tavr = tavr + work_old(is,js)
 655                        iavr = iavr+1
 656                      endif
 657                    enddo
 658                  enddo
 659                  if (iavr.gt.0) then
 660                    if (iavr.eq.1) then
 661 ! spreading is not allowed if the only valid neighbor is a corner point
 662 ! otherwise an ill posed cellular automaton is established leading to
 663 ! a spreading of constant values in diagonal direction
 664 ! if all corner points are blanked the valid neighbor must be a direct one
 665 ! and spreading is allowed
 666                      if (work_old(inl,jnu).eq.blank.and.&
 667                          work_old(inr,jnu).eq.blank.and.&
 668                          work_old(inr,jnl).eq.blank.and.&
 669                          work_old(inl,jnl).eq.blank) then
 670                            work_new(i,j)=tavr/iavr
 671                            ready = .false.
 672                      endif
 673                   else
 674                     work_new(i,j)=tavr/iavr
 675                     ready = .false.
 676                   endif
 677                 endif
 678               endif
 679             enddo ! j
 680           enddo   ! i
 681 ! save changes made during this pass to work_old
 682           work_old(:,:)=work_new(:,:)
 683           if(ipass.eq.maxpass) ready=.true.
 684         enddo !while (.not.ready)
 685         fi(:,:) = work_new(:,:)
 686       else
 687          do while (.not.ready)
 688            ipass = ipass+1
 689            ready = .true.
 690            do j=jcs, jce
 691              do i=ics, ice
 692                if (work_old(i,j).le.blank) then
 693                  tavr=0.
 694                  iavr=0
 695                  inl = max(i-1,ics)
 696                  inr = min(i+1,ice)
 697                  jnl = max(j-1,jcs)
 698                  jnu = min(j+1,jce)
 699                  do is=inl,inr
 700                    do js=jnl,jnu
 701                      if (work_old(is,js).gt.blank) then
 702                        tavr = tavr + work_old(is,js)
 703                        iavr = iavr+1
 704                      endif
 705                    enddo
 706                  enddo
 707                  if (iavr.gt.0) then
 708                    if (iavr.eq.1) then
 709 ! spreading is not allowed if the only valid neighbor is a corner point
 710 ! otherwise an ill posed cellular automaton is established leading to
 711 ! a spreading of constant values in diagonal direction
 712 ! if all corner points are blanked the valid neighbor must be a direct one
 713 ! and spreading is allowed
 714                      if (work_old(inl,jnu).le.blank.and. &
 715                          work_old(inr,jnu).le.blank.and. &
 716                          work_old(inr,jnl).le.blank.and. &
 717                          work_old(inl,jnl).le.blank) then
 718                            work_new(i,j)=tavr/iavr
 719                            ready = .false.
 720                      endif
 721                   else
 722                     work_new(i,j)=tavr/iavr
 723                     ready = .false.
 724                   endif
 725                 endif
 726               endif
 727             enddo ! j
 728           enddo   ! i
 729 ! save changes made during this pass to work_old
 730           work_old(:,:)=work_new(:,:)
 731           if(ipass.eq.maxpass) ready=.true.
 732         enddo !while (.not.ready)
 733         fi(:,:) = work_new(:,:)
 734       endif
 735       return
 736     end subroutine fill_xy
 737
 738   subroutine horiz_interp_bicubic_del( Interp )
 739
 740     type (horiz_interp_type), intent(inout) :: Interp
 741
 742     if(associated(Interp%rat_x))  deallocate ( Interp%rat_x )
 743     if(associated(Interp%rat_y))  deallocate ( Interp%rat_y )
 744     if(associated(Interp%lon_in)) deallocate ( Interp%lon_in )
 745     if(associated(Interp%lat_in)) deallocate ( Interp%lat_in )
 746     if(associated(Interp%i_lon))  deallocate ( Interp%i_lon )
 747     if(associated(Interp%j_lat))  deallocate ( Interp%j_lat )
 748     if(associated(Interp%wti))    deallocate ( Interp%wti )
 749
 750   end subroutine horiz_interp_bicubic_del
 751
 752 end module horiz_interp_bicubic_mod
 753 !> @}
 754 ! close documentation