interpol_splines.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
!     #########
      SUBROUTINE interpol_splines(KLUOUT,KCODE,PX,PY,PFIELD,KREP)
!     #########################################################
!
!!**** *INTERPOL_SPLINES* interpolates with spline f77 programs a 2D field
!!                           from all grid points valid values
!!
!!    PURPOSE
!!    -------
!!
!!    The points are all on only one grid (defined with the coordinates
!!    of all the points). The code to apply for each point is:
!!
!!    KCODE>0 : data point (with field valid for interpolation)
!!    KCODE=-1: point to ignore
!!    KCODE=0 : point to interpolate
!!
!!
!!
!!    METHOD
!!    ------
!!
!!    EXTERNAL
!!    --------
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!
!!
!!    REFERENCE
!!    ---------
!!
!!    AUTHOR
!!    ------
!!
!!    V. Masson          Meteo-France
!!
!!    MODIFICATION
!!    ------------
!!
!!    Original    19/03/95
!!    Modification
!!     25/05/96   (V Masson)  W and IW defined in MODD_SPLINESWORK  
!!     04/07/96   (V Masson)  call with 1d dummy arguments
!!     30/10/96   (V Masson)  add deallocations
!!     05/08/97   (V Masson)  output listing as dummy argument
!!     17/09/97   (V Masson)  routine performs ONLY the splines interpolation
!!     19/12/97   (V Masson)  possibility to use the same splines for collocated
!!                            fields
!----------------------------------------------------------------------------
!
!*    0.     DECLARATION
!            -----------
!
USE modd_surf_par,    ONLY : xundef
!
USE mode_splines
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE modi_abor1_sfx
!
IMPLICIT NONE
!
!*    0.1    Declaration of arguments
!            ------------------------
!
INTEGER,               INTENT(IN)     :: kluout   ! output listing
INTEGER,DIMENSION(:),  INTENT(INOUT)  :: kcode    ! code for each point
                                                  ! >0 point used for interpolation
                                                  !  0 point to interpolate
                                                  ! -1 point not used
                                                  ! -2 point not used
!                                                 ! -3 if spline is no computed
!                                                 ! for this point
REAL,   DIMENSION(:),  INTENT(IN)     :: px       ! x of each grid mesh.
REAL,   DIMENSION(:),  INTENT(IN)     :: py       ! y of each grid mesh.
REAL,   DIMENSION(:,:),INTENT(INOUT)  :: pfield   ! pgd field on grid mesh.
INTEGER,               INTENT(OUT)    :: krep     ! error flag
!
!*    0.2    Declaration of local variables
!            ------------------------------
!
INTEGER, PARAMETER      :: nd = 2                !
INTEGER                 :: iopt     ! 0: splines recomputed
!                                   ! 1: splines deduced from W and IW
INTEGER                 :: jk                   ! loop control
INTEGER                 :: jin                  ! loop control
INTEGER                 :: jout                 ! loop control
INTEGER                 :: iverb = 10           ! verbosity level
INTEGER                 :: ifield               ! number of colocated fields
INTEGER, PARAMETER      :: idatamax = 2000      ! maximum number of data per blok
!                                               ! WARNING, if IDATAMAX is modified
!                                               ! MODIFY nmax in all splines.f routines
INTEGER, PARAMETER      :: isdmax = 20          ! maximum number of subdomain in one direction
!                                               ! WARNING, if ISDMAX is modified
!                                               ! MODIFY ndmax in all splines.f routines
INTEGER, PARAMETER      :: immax = 3            ! maximum number of monomes (order 2)
INTEGER                 :: idata                ! total number of data
REAL, DIMENSION(2,ND)             :: zxdom     ! coordinates of the limits of the domain
REAL, DIMENSION(:,:), ALLOCATABLE  :: zxdata    ! coordinates of the available data
REAL, DIMENSION(:),   ALLOCATABLE  :: zdata     ! available data
INTEGER,DIMENSION(ISDMAX,ISDMAX)   :: ndomdata  ! number of data in a subdomain

INTEGER                            :: iout      ! total number of points to interpolate
REAL, DIMENSION(:,:), ALLOCATABLE  :: zxout     ! coordinates of these points
REAL, DIMENSION(:),   ALLOCATABLE  :: zvalout   ! values of these points
REAL                               :: zminval   ! minimum value of the field
REAL                               :: zmaxval   ! maximum value of the field
!
!*    0.3    Declaration for f77 routines
!            ----------------------------
!
INTEGER, PARAMETER    :: iorder = 2                 ! order of the spline
INTEGER, PARAMETER    :: im = (iorder*(iorder+1))/2 ! number of monomes
INTEGER               :: isdi,isdj                  ! number of subdomains in each
!                                                   ! direction for spline computation
INTEGER               :: iiorder
INTEGER               :: iim
INTEGER               :: jsdi,jsdj                  ! loop controls on subdomains
INTEGER               :: iinter                     ! inverse ratio of the intersection of
!                                                   ! 2 subdomains over one of the subdomain
REAL, DIMENSION(ND,ND)  :: zg                         ! matrix of correspondance between
!                                                   ! x and y units
REAL                  :: zp, zs2
!
REAL,    DIMENSION(:,:,:), ALLOCATABLE :: zc        ! coefficients of the spline
REAL,    DIMENSION(IDATAMAX+3*IMMAX)   :: zwe       ! work array
!
!INTEGER, DIMENSION(2,IM)   :: IW
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
!
!*    1.     Miscellaneous Initializations
!            -----------------------------
!
IF (lhook) CALL dr_hook('INTERPOL_SPLINES',0,zhook_handle)
ifield=SIZE(pfield,2)
!
!IW = 0
iiorder=iorder
iim=im
!-------------------------------------------------------------------------------
!
!*    2.     Splines initializations
!            -----------------------
!
!
!*    2.1    Coordinates of input data points
!            --------------------------------
!
idata=count(kcode>0)
ALLOCATE(zxdata(nd,idata))
ALLOCATE(zdata(idata))
!
zxdata(1,:)=pack(px(:),mask=(kcode(:)>0))
zxdata(2,:)=pack(py(:),mask=(kcode(:)>0))
!
!*    2.2    Coordinates of domain limits
!            ----------------------------
!
zxdom(1,1)=minval(px(:),mask=kcode(:)>=0)
zxdom(2,1)=maxval(px(:),mask=kcode(:)>=0)
zxdom(1,2)=minval(py(:),mask=kcode(:)>=0)
zxdom(2,2)=maxval(py(:),mask=kcode(:)>=0)
!
!*    2.3    Coordinates of output interpolated points
!            -----------------------------------------
!
iout=count(kcode==0)
ALLOCATE(zxout(nd,iout))
ALLOCATE(zvalout(iout))
!
zxout(1,:)=pack(px(:),mask=(kcode(:)==0))
zxout(2,:)=pack(py(:),mask=(kcode(:)==0))
!
!*    2.4    Matrix of unit correspondance
!            -----------------------------
!
zg(1,1)=1.
zg(2,2)=1.
zg(1,2)=0.
zg(2,1)=0.
!
!-------------------------------------------------------------------------------
!
!*    3.     Subdomains generation
!            ---------------------
!
iinter=4
!
isdi=max(1,min(int( sqrt( idata/100+1 -1.e-6 )),isdmax))
isdj=isdi
!
!
subdomains: DO
  DO jsdi=1,isdi
    DO jsdj=1,isdj
      ndomdata(jsdi,jsdj)= count( (kcode(:)>0) .AND.                                          &
         (      ( px(:) >= zxdom(1,1) + (jsdi-1-0.5/(iinter-1))*(zxdom(2,1)-zxdom(1,1))/isdi)   &
         .AND.  ( px(:) <= zxdom(1,1) + (jsdi  +0.5/(iinter-1))*(zxdom(2,1)-zxdom(1,1))/isdi)   &
         .AND.  ( py(:) >= zxdom(1,2) + (jsdj-1-0.5/(iinter-1))*(zxdom(2,2)-zxdom(1,2))/isdj)   &
         .AND.  ( py(:) <= zxdom(1,2) + (jsdj  +0.5/(iinter-1))*(zxdom(2,2)-zxdom(1,2))/isdj) )            &
                                                           )  
    END DO
  END DO
  IF ( (isdi<isdmax) .AND. (isdj<isdmax) .AND. any(ndomdata(1:isdi,1:isdj) > idatamax/2) ) THEN
    isdi=isdi+1
    isdj=isdj+1
    cycle subdomains
  END IF
  EXIT subdomains
END DO subdomains
!
IF ( any(ndomdata(1:isdi,1:isdj) > idatamax) ) THEN
  WRITE(kluout,*) '**********************************************************'
  WRITE(kluout,*) '* ERROR during data interpolation                        *'
  WRITE(kluout,*) '* interpolation by splines routines is impossible        *'
  WRITE(kluout,*) '* please use an input data file with a better resolution *'
  WRITE(kluout,*) '**********************************************************'
  CALL abor1_sfx('INTERPOL_SPLINES: ERROR DURING DATA INTERPOLATION BY SPLINES')
END IF
!
WRITE(kluout,*) '----------------------------------'
!
!-------------------------------------------------------------------------------
!
!*    4.     Computation of the spline in the subdomains (f77 routine)
!            -------------------------------------------
!
!*    4.1    Output printing
!            ---------------
!
  IF (iverb>=7) THEN
    DO jsdi=1,isdi
      DO jsdj=1,isdj
        WRITE(kluout,*) '  domain sdi=',jsdi,' sdj=',jsdj,' : ',     &
                           ndomdata(jsdi,jsdj),' data points available'  
      END DO
    END DO
  END IF
!
!*    4.2    Choice between splines computation or deduction
!            -----------------------------------------------
!
  IF (iverb>=7) WRITE(kluout,*) ' splines computations:'
!
!*    4.3    Allocations
!            -----------
!
  ALLOCATE(zc( idatamax+immax , isdi , isdj ))
  iopt=0
!
!*    4.4    Loop on fields
!            --------------
!
  DO jk=1,ifield
!
    zdata(:)=pack(pfield(:,jk),mask=(kcode(:)>0))
!
!*    4.5    special case of uniform input data
!            -----------------------------
!
    zminval = minval(zdata)
    zmaxval = maxval(zdata)
!-------------------------------------------------------------------------------
    IF ( abs(zmaxval - zminval) <= 1.e-10 * max(abs(zmaxval),abs(zminval)) ) THEN
!-------------------------------------------------------------------------------
!
      zvalout(:)=maxval(zdata)
!
!-------------------------------------------------------------------------------
    ELSE
!-------------------------------------------------------------------------------
!
!*    4.6    Call for splines coefficients
!            -----------------------------
!
!
      zp = 0.
      zs2 = 0.
!
      CALL splb2c(iiorder,iim,zxdata,zg,zdata,zs2,zp,0,iopt,isdi,isdj,iinter,zxdom,zc,krep)
!      CALL SPLB2C(IDATA,ZXDATA,ZG,ZDATA,0,0.,ZP,IORDER,IM,IOPT,ISDI,ISDJ,IINTER,ZXDOM, &
!                   ZC,KREP,IW,W) 
!
      iopt=1
!
      IF (krep/=0) THEN
        WRITE(kluout,*) 'Routine INTERPOL_SPLINES: error in SPLB2C, KREP=',krep
        WRITE(kluout,*) 'The field is interpolated from 3 nearest points'
        IF (ALLOCATED(zxdata))  DEALLOCATE(zxdata)
        IF (ALLOCATED(zdata))   DEALLOCATE(zdata)
        IF (ALLOCATED(zc))      DEALLOCATE(zc)
        IF (ALLOCATED(zvalout)) DEALLOCATE(zvalout)
        IF (ALLOCATED(zxout))   DEALLOCATE(zxout)
        IF (lhook) CALL dr_hook('INTERPOL_SPLINES',1,zhook_handle)
        RETURN
      END IF
!
!-------------------------------------------------------------------------------
!
!*    5.     Evaluation of the spline (f77 routine)
!            ------------------------
!
      zvalout(:)=xundef
    !
      IF (iverb>=7) WRITE(kluout,*) ' splines evaluations'
      !CALL SPLB2E(IORDER,IM,ZG,ISDI,ISDJ,ZC,IOUT,ZXOUT,ZVALOUT,ZWE,IW,W)
      CALL splb2e(iiorder,iim,isdi,isdj,zg,zc,zxout,zvalout)
!
!-------------------------------------------------------------------------------
    END IF
!-------------------------------------------------------------------------------
!
!*    6.     filling grid points
!            -------------------
!
    jout=0
    DO jin=1,SIZE(pfield,1)
      IF (kcode(jin)/=0) cycle
      jout=jout+1
      pfield(jin,jk) = zvalout(jout)
      pfield(jin,jk) = max(zminval, min(zmaxval,pfield(jin,jk))) ! control of extrema
    END DO
!
!-------------------------------------------------------------------------------
  END DO
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
!
!*    8.     Deallocations
!            -------------
!
IF (ALLOCATED(zxdata))  DEALLOCATE(zxdata)
IF (ALLOCATED(zdata))   DEALLOCATE(zdata)
IF (ALLOCATED(zc))      DEALLOCATE(zc)
IF (ALLOCATED(zvalout)) DEALLOCATE(zvalout)
IF (ALLOCATED(zxout))   DEALLOCATE(zxout)
IF (lhook) CALL dr_hook('INTERPOL_SPLINES',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE interpol_splines