interpol_npts.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_npts (UG, U, &
                                hprogram,kluout,knpts,kcode,px,py,pfield,knear_nbr)
!     #########################################################
!
!!**** *INTERPOL_NPTS* interpolates with ###ine 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    03/2004
!!    Modification
!!    B. Decharme  2014  scan all point case if gaussien grid or NHALO = 0
!----------------------------------------------------------------------------
!
!*    0.     DECLARATION
!            -----------
!
!
USE modd_surf_atm_grid_n, ONLY : surf_atm_grid_t
USE modd_surf_atm_n, ONLY : surf_atm_t
!
USE modd_surf_par,       ONLY : xundef
!
USE modi_get_interp_halo
USE modi_get_near_meshes
USE modi_sum_on_all_procs
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*    0.1    Declaration of arguments
!            ------------------------
!
!
TYPE(surf_atm_grid_t), INTENT(INOUT) :: ug
TYPE(surf_atm_t), INTENT(INOUT) :: u
!
 CHARACTER(LEN=6),      INTENT(IN)     :: hprogram ! host program
INTEGER,               INTENT(IN)     :: kluout   ! output listing
INTEGER,               INTENT(IN)     :: knpts    ! number of points to interpolate with
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(IN) :: knear_nbr
!
!*    0.2    Declaration of local variables
!            ------------------------------
!
INTEGER                                     :: il ! number of points
INTEGER                                     :: jd ! data point index
INTEGER                                     :: js ! loop counter on data points
INTEGER                                     :: jl ! loop counter on points to initialize
INTEGER                                     :: jp, jpp ! loops counter on KNPTS
REAL :: zdist ! square distance between two interpolating and interpolated points
REAL, DIMENSION(0:KNPTS)                :: zndist ! 3 nearest square distances
REAL, DIMENSION(0:KNPTS,SIZE(PFIELD,2)) :: znval  ! 3 corresponding field values
REAL, DIMENSION(SIZE(PFIELD,2))         :: zsum
!
INTEGER :: ihalo
INTEGER                            :: jlist          ! loop counter on points to interpolate
INTEGER                            :: icount         ! counter
INTEGER                            :: inpts
INTEGER                            :: iscan          ! number of points to scan
INTEGER                            :: iscan_all      ! number of data points
INTEGER, DIMENSION(SIZE(PFIELD,1)) :: iindex       ! list of index to scan
INTEGER, DIMENSION(SIZE(PFIELD,1)) :: iindex_all   ! list of all data points
INTEGER, DIMENSION(SIZE(KCODE))    :: isize
INTEGER                            :: isize0
!
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('INTERPOL_NPTS',0,zhook_handle)
!
il = SIZE(pfield,1)
!
iindex(:) = 0
!
 CALL get_interp_halo(hprogram,ug%CGRID,ihalo)
!
IF(ug%CGRID=='GAUSS'.OR.ihalo==0)THEN
!
  isize(:) = 1.
  isize0 = sum_on_all_procs(hprogram,ug%CGRID,isize(:)==1)
!
  iindex_all(:) = 0
!        
  iscan_all = count(kcode(:)>0)
!
  js = 0
  DO jd=1,il
    IF (kcode(jd)>0) THEN
      js = js+1
      iindex_all(js) = jd
    END IF
  END DO 
!
ELSEIF (.NOT.ASSOCIATED(ug%NNEAR)) THEN
  ALLOCATE(ug%NNEAR(il,knear_nbr))
  ug%NNEAR(:,:) = 0
  CALL get_near_meshes(ug%CGRID,ug%NGRID_PAR,u%NSIZE_FULL,ug%XGRID_PAR,knear_nbr,ug%NNEAR)
ENDIF
!
DO jl=1,il

  IF (kcode(jl)/=0) cycle

  zndist(1:knpts) = 1.e20
  zndist(0) = 0.
  znval(0:knpts,:) = 0.
  !
  IF(ug%CGRID=='GAUSS'.OR.ihalo==0)THEN
    !
    IF (u%NDIM_FULL/=isize0) THEN
      ! point can not be interpolated further than halo in multiprocessor run
      kcode(jl) = -4
      cycle
    END IF
    inpts     = knpts
    iscan     = iscan_all
    iindex(:) = iindex_all(:)
    !
  ELSE
    !
    icount = 0
    DO jd=1,knear_nbr
      IF (ug%NNEAR(jl,jd)>0) THEN
        IF (kcode(ug%NNEAR(jl,jd))>0) THEN  
          icount = icount+1
          iindex(icount) = ug%NNEAR(jl,jd)
        END IF
      END IF
    END DO
    !
    !IF (ICOUNT>=1) THEN
    IF (icount>=knpts) THEN
      iscan = icount
      !INPTS = MIN(ICOUNT,KNPTS)
      inpts = knpts
    ELSEIF (knear_nbr>=u%NDIM_FULL .AND. icount>=1) THEN
      iscan = icount
      inpts = icount      
    ELSE
      kcode(jl) = -4
      cycle
    END IF
    !
  END IF
  !
  DO js=1,iscan
    !
    jd = iindex(js)
    !
    zdist=  ( ( px(jd)-px(jl) ) ** 2 ) + ( ( py(jd)-py(jl) ) ** 2 )
    !
    IF ( zdist>zndist(inpts) ) cycle
    !
    DO jp = inpts,1,-1
      !
      IF ( zdist>zndist(jp-1) ) THEN
        !
        IF ( jp<inpts ) THEN
          DO jpp = inpts,jp+1,-1
            zndist(jpp)  = zndist(jpp-1)
            znval(jpp,:) = znval(jpp-1,:)
          ENDDO
        ENDIF
        !
        zndist(jp)  = zdist
        znval(jp,:) = pfield(jd,:)
        !
        EXIT
        !
      ENDIF
      !
    ENDDO
    !
  ENDDO
  !
  zndist(:) = sqrt(zndist(:))
  !
  pfield(jl,:) = 0.
  zsum(:) = 0.
  DO jp = 1, inpts
    pfield(jl,:) = pfield(jl,:) + znval(jp,:)/zndist(jp)
    zsum(:) = zsum(:) + 1./zndist(jp)
  ENDDO
  pfield(jl,:) = pfield(jl,:) / zsum(:)
  !
END DO
!
IF (lhook) CALL dr_hook('INTERPOL_NPTS',1,zhook_handle)
!-------------------------------------------------------------------------------
!
END SUBROUTINE interpol_npts