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_surfex_omp, ONLY : nblock
USE modd_surfex_mpi, ONLY : nrank, npio, ncomm, nproc, idx_i, nindex,nnum
USE modd_surf_par,       ONLY : xundef
!
USE modi_gather_and_write_mpi
USE modi_read_and_send_mpi
USE modi_get_near_meshes
USE modi_sum_on_all_procs
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
#ifdef SFX_MNH
USE modd_io_ll, ONLY : isp, isnproc
USE modd_var_ll, ONLY : nmnh_comm_world
USE mode_gather_ll
USE mode_fd_ll, ONLY : getfd, fd_ll
USE mode_tools_ll,     ONLY : get_globaldims_ll
USE modd_parameters_ll,ONLY : jphext
!
USE modd_io_surf_mnh, ONLY : niu, nju
#endif
!
IMPLICIT NONE
!
#if defined(SFX_MPI) || defined(SFX_MNH)
include "mpif.h"
#endif
!
!*    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
!            ------------------------------
!
REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZFIELD, ZFIELD2
REAL, DIMENSION(:,:,:), ALLOCATABLE   :: ZFIELD3
REAL, DIMENSION(:,:), ALLOCATABLE     :: ZNDIST  ! 3 nearest square distances
REAL, DIMENSION(:,:), ALLOCATABLE     :: ZNVAL  ! 3 corresponding field values
REAL, DIMENSION(:), ALLOCATABLE       :: ZX, ZY
REAL, DIMENSION(:), ALLOCATABLE       :: ZDIST ! square distance between two interpolating and interpolated points
!
INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: ININD0, ININD_ALL  ! 3 corresponding field values
INTEGER, DIMENSION(:,:), ALLOCATABLE   :: ININD
INTEGER, DIMENSION(:), ALLOCATABLE     :: IINDEX       ! list of index to scan
INTEGER, DIMENSION(:), ALLOCATABLE     :: ISIZE, ISIZE_TOT
!
REAL    :: ZSUM
INTEGER :: IP, ICPT, IL1, JL, JP, JK, JKK 
!
INTEGER :: ICOUNT, IL2   ! counter
INTEGER :: INPTS
INTEGER :: INFOMPI, IDIM_FULL, INEAR_NBR, IOLD
!
#if defined(SFX_MPI) || defined(SFX_MNH)
INTEGER, DIMENSION(MPI_STATUS_SIZE) :: ISTATUS
#endif
!
#ifdef SFX_MNH
TYPE(fd_ll), POINTER                     :: TZFD
!
REAL, DIMENSION(:,:), ALLOCATABLE :: ZCOORD_2D, ZCOORD_2D_ALL
INTEGER, DIMENSION(:,:), ALLOCATABLE :: ISIZE_2D, ISIZE_2D_ALL, INUM_2D, INUM_2D_ALL, IINDEX_2D, IINDEX_2D_ALL
INTEGER, DIMENSION(:), ALLOCATABLE :: INUM_1D, IINDEX_1D, INUM_TOT, IINDEX_TOT
INTEGER :: IIMAX, IJMAX, IIU, IJU, JI
INTEGER :: IRANK_SAVE, IPROC_SAVE, IPIO_SAVE, ICOMM_SAVE
#endif
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('INTERPOL_NPTS_1',0,zhook_handle)
!
iold = 0
!NHALO /= 0
IF ( knear_nbr/=u%NDIM_FULL ) THEN
  iold = 1
ELSE
! case mesonh and NHALO = 0
#ifdef SFX_MNH
  iold = 2
! case NHALO = 0
#else
  iold = 1
#endif
ENDIF
!
IF (iold==2) THEN
  !
#ifdef SFX_MNH
  !physical dimensions of the task
  CALL get_dim_phys_ll('B',iiu,iju)
  ! total dimensions
  CALL get_globaldims_ll(iimax,ijmax)
  idim_full = (iimax) * (ijmax)
  inear_nbr = idim_full
  !
  tzfd=>getfd(nmnh_comm_world)
  !
  ! on sauve les infos de bases
  irank_save = nrank
  iproc_save = nproc
  ipio_save  = npio
  icomm_save = ncomm
  !
  ! on met les infos de mésonh
  nrank = isp-1
  nproc = isnproc
  npio  = tzfd%OWNER-1
  ncomm = tzfd%COMM
  ! 
  !KCODE to 2D
  ALLOCATE(isize_2d(iiu+2*jphext,iju+2*jphext))
  isize_2d(:,:) = -1
  isize_2d(1+jphext:iiu+jphext,1+jphext:iju+jphext) = reshape(kcode, (/ iiu,iju /) )
  ! tasks to whole domaine
  ALLOCATE(isize_2d_all(iimax+2*jphext,ijmax+2*jphext))
  CALL gather_xyfield(isize_2d,isize_2d_all,tzfd%OWNER,tzfd%COMM)
  DEALLOCATE(isize_2d)
  ALLOCATE(isize_tot(idim_full))
  isize_tot = reshape(isize_2d_all(1+jphext:iimax+jphext,1+jphext:ijmax+jphext),(/idim_full/))
  DEALLOCATE(isize_2d_all)
#endif
  !
ELSEIF (iold==1) THEN
  !
  idim_full = u%NDIM_FULL
  inear_nbr = knear_nbr
  !
  ALLOCATE(isize_tot(idim_full))
  CALL gather_and_write_mpi(kcode,isize_tot)
  !
ENDIF
!
!...known by all tasks
IF (nproc>1) THEN
#if defined(SFX_MPI) || defined(SFX_MNH)
  CALL mpi_bcast(isize_tot,idim_full*kind(isize_tot)/4,mpi_integer,npio,ncomm,infompi)
#endif
ENDIF
!
IF (all(isize_tot/=0)) THEN
  IF (iold==2) THEN
#ifdef SFX_MNH
    nrank = irank_save
    nproc = iproc_save
    npio = ipio_save
    ncomm = icomm_save
#endif
  ENDIF
  DEALLOCATE(isize_tot)
  CALL dr_hook('INTERPOL_NPTS_1',1,zhook_handle)
  RETURN
ENDIF
!
!
ip = count(kcode(:)==0)
!
il1 = SIZE(pfield,1)
il2 = SIZE(pfield,2)
!
IF (iold==2) THEN
  !
#ifdef SFX_MNH
!
  ! NNUM
  ALLOCATE(inum_1d(u%NDIM_FULL))
  DO ji = 1,u%NDIM_FULL
    inum_1d(ji) = ji
  ENDDO
  ! INUM_1D to 2D 
  ALLOCATE(inum_2d(iiu+2*jphext,iju+2*jphext))
  inum_2d = 0
  inum_2d(1+jphext:iiu+jphext,1+jphext:iju+jphext) = reshape(inum_1d, (/ iiu,iju /) )
  DEALLOCATE(inum_1d)
  ! tasks to whole domaine
  ALLOCATE(inum_2d_all(iimax+2*jphext,ijmax+2*jphext))
  CALL gather_xyfield(inum_2d,inum_2d_all,tzfd%OWNER,tzfd%COMM)
  DEALLOCATE(inum_2d)
  ALLOCATE(inum_tot(idim_full))
  inum_tot = reshape(inum_2d_all(1+jphext:iimax+jphext,1+jphext:ijmax+jphext),(/idim_full/))
  DEALLOCATE(inum_2d_all)
  !
  ! NINDEX
  ALLOCATE(iindex_1d(u%NDIM_FULL))
  iindex_1d(:) = isp - 1
  ! INDEX_1D to 2D
  ALLOCATE(iindex_2d(iiu+2*jphext,iju+2*jphext))
  iindex_2d(:,:) = 0
  iindex_2d(1+jphext:iiu+jphext,1+jphext:iju+jphext) = reshape(iindex_1d, (/ iiu,iju /) )
  DEALLOCATE(iindex_1d)
  ! tasks to whole domaine
  ALLOCATE(iindex_2d_all(iimax+2*jphext,ijmax+2*jphext))
  CALL gather_xyfield(iindex_2d,iindex_2d_all,tzfd%OWNER,tzfd%COMM)
  DEALLOCATE(iindex_2d)
  ALLOCATE(iindex_tot(idim_full))
  iindex_tot = reshape(iindex_2d_all(1+jphext:iimax+jphext,1+jphext:ijmax+jphext),(/idim_full/))
  DEALLOCATE(iindex_2d_all)
  !
  !PX, PY: coordinates of all the points
  !ZX, ZY: coordinates of the points in this task
  ALLOCATE(zcoord_2d(iiu+2*jphext,iju+2*jphext))
  ALLOCATE(zcoord_2d_all(iimax+2*jphext,ijmax+2*jphext))
  !
  zcoord_2d = 0.
  zcoord_2d(1+jphext:iiu+jphext,1+jphext:iju+jphext) = reshape(px, (/ iiu,iju /) )
  CALL gather_xyfield(zcoord_2d,zcoord_2d_all,tzfd%OWNER,tzfd%COMM)
  ALLOCATE(zx(idim_full))
  zx = reshape(zcoord_2d_all(1+jphext:iimax+jphext,1+jphext:ijmax+jphext),(/idim_full/))
  !
  zcoord_2d = 0.
  zcoord_2d(1+jphext:iiu+jphext,1+jphext:iju+jphext) = reshape(py, (/ iiu,iju  /) )
  CALL gather_xyfield(zcoord_2d,zcoord_2d_all,tzfd%OWNER,tzfd%COMM)
  ALLOCATE(zy(idim_full))
  zy = reshape(zcoord_2d_all(1+jphext:iimax+jphext,1+jphext:ijmax+jphext),(/idim_full/))
  !
  DEALLOCATE(zcoord_2d,zcoord_2d_all)
  !
  IF (nproc>1) THEN
    CALL mpi_bcast(inum_tot,idim_full*kind(inum_tot)/4,mpi_integer,npio,ncomm,infompi)
    CALL mpi_bcast(iindex_tot,idim_full*kind(iindex_tot)/4,mpi_integer,npio,ncomm,infompi)
    CALL mpi_bcast(zx,idim_full*kind(zx)/4,mpi_float,npio,ncomm,infompi)
    CALL mpi_bcast(zy,idim_full*kind(zy)/4,mpi_float,npio,ncomm,infompi)
  ENDIF
!
#endif
  !
ELSEIF (iold==1) THEN
  !
  !PX, PY: coordinates of all the points
  !ZX, ZY: coordinates of the points in this task
  ALLOCATE(zx(il1),zy(il1))
  CALL read_and_send_mpi(px,zx)
  CALL read_and_send_mpi(py,zy)
  !
ENDIF
!
IF (inear_nbr/=idim_full) THEN
  IF (.NOT.ASSOCIATED(ug%NNEAR)) THEN
    ALLOCATE(ug%NNEAR(il1,inear_nbr))
    !seach near meshes in the complete grid (xgrid_full_par) for this task
    CALL get_near_meshes(ug%G%CGRID,ug%NGRID_FULL_PAR,idim_full,ug%XGRID_FULL_PAR,inear_nbr,ug%NNEAR)
  ENDIF
ENDIF 
!
IF (lhook) CALL dr_hook('INTERPOL_NPTS_1',1,zhook_handle)
IF (lhook) CALL dr_hook('INTERPOL_NPTS_2',0,zhook_handle)
!
ALLOCATE(iindex(inear_nbr))
iindex(:) = 0
!
IF (inear_nbr==idim_full) THEN

  icount = 0
  DO jl=1,idim_full
    !is the neareast point available to interpolation
    IF (isize_tot(jl)>0) THEN  
      icount = icount+1
      iindex(icount) = jl
    END IF
  END DO

  !did we found enough points for interpolate
  IF (icount>=1) THEN
    inpts = min(knpts,icount)
  ELSE
    WHERE(kcode(:)==0) kcode(:) = -4
  END IF
  ALLOCATE(zdist(icount))
ELSE
  ALLOCATE(zdist(inear_nbr))
ENDIF
!
zdist(:) = 0.
!
IF (lhook) CALL dr_hook('INTERPOL_NPTS_2',1,zhook_handle)
IF (lhook) CALL dr_hook('INTERPOL_NPTS_3',0,zhook_handle)
!
!indexes of points used for interpolation, for each point to interpolate
ALLOCATE(inind(ip,knpts))
inind(:,:) = 0
!distances of the points used for interpolation
ALLOCATE(zndist(ip,0:knpts))
zndist(:,1:knpts) = 1.e20
zndist(:,0) = 0.
!
icpt=0
!loop on points for this task
DO jl=1,il1
  !
  !does this point need to be interpolated? 
  IF (kcode(jl)/=0) cycle
  !
  IF (inear_nbr/=idim_full) THEN

    icount = 0
    DO jk=1,inear_nbr
      IF (ug%NNEAR(jl,jk)>0) THEN
        !is the neareast point available to interpolation
        IF (isize_tot(ug%NNEAR(jl,jk))>0) THEN  
          icount = icount+1
          iindex(icount) = ug%NNEAR(jl,jk)
        END IF
      END IF
    END DO
    !
    !did we found enough points for interpolate
    IF (icount>=knpts) THEN
      inpts = knpts
    ELSE
      kcode(jl) = -4
      cycle
    END IF
    !
  ENDIF
  !
  !one point more to interpolate
  icpt = icpt + 1
  !
  IF (iold==2) THEN
#ifdef SFX_MNH
    zdist(1:icount) = (zx(iindex(1:icount))-px(jl))**2 + (zy(iindex(1:icount))-py(jl))**2
#endif
  ELSE
    zdist(1:icount) = (px(iindex(1:icount))-zx(jl))**2 + (py(iindex(1:icount))-zy(jl))**2
  ENDIF
  !
  DO jp = 1,icount
    !
    IF (zdist(jp)>zndist(icpt,inpts)) cycle
    !
    DO jk = inpts,1,-1
      !
      IF ( zdist(jp)>zndist(icpt,jk-1) ) THEN
        !
        IF ( jk<inpts ) THEN
          DO jkk = inpts,jk+1,-1
            zndist(icpt,jkk) = zndist(icpt,jkk-1)
            inind(icpt,jkk) = inind(icpt,jkk-1)
          ENDDO
        ENDIF
        !
        !distances and indexes of points used to interpolate are saved
        zndist(icpt,jk) = zdist(jp)
        inind(icpt,jk) = iindex(jp)
        !
        EXIT
        !
      ENDIF
      !
    ENDDO
    !
  ENDDO
  !
ENDDO
!
IF (lhook) CALL dr_hook('INTERPOL_NPTS_3',1,zhook_handle)
IF (lhook) CALL dr_hook('INTERPOL_NPTS_4',0,zhook_handle)
!
DEALLOCATE(iindex,zx,zy,isize_tot,zdist)
!
zndist(:,:) = sqrt(zndist(:,:))
!
ALLOCATE(isize(0:nproc-1))

!numbers of points to interpolated are gathered
IF (nproc>1) THEN
#if defined(SFX_MPI) || defined(SFX_MNH)
  CALL mpi_allgather(icpt,kind(icpt)/4,mpi_integer,&
                     isize,kind(isize)/4,mpi_integer,ncomm,infompi)
#endif
ELSE
  isize(:) = icpt
ENDIF
!
!this array contains, for each point to interpolate, 
!the correspondance between the task where is the point to use for interpolation
!(NINDEX) and its index in this task (NNUM)
ALLOCATE(inind0(maxval(isize),knpts,0:nproc-1))
inind0(:,:,:) = 0
!
!number of points needed to interpolated
DO jl=1,knpts
  !number of points to interpolated
  DO jp=1,icpt
    !index of the point needed in the whole grid
    jk = inind(jp,jl)
    !inind0 contains the task and the index in this task for this point
    IF (iold==2) THEN
#ifdef SFX_MNH
      IF (jk/=0) inind0(jp,jl,iindex_tot(jk)) = inum_tot(jk)
#endif
    ELSEIF (iold==1) THEN
      IF (jk/=0) inind0(jp,jl,nindex(jk)) = nnum(jk)
    ENDIF
  ENDDO
ENDDO
!
!for each task, points needed and located in this task are gathered
ALLOCATE(inind_all(maxval(isize),knpts,0:nproc-1))
!
IF (nproc>1) THEN
  !for each task
  DO jp=0,nproc-1
#if defined(SFX_MPI) || defined(SFX_MNH)  
    !inind_all receives from all tasks the points they need that are
    !located in it
    CALL mpi_gather(inind0(:,:,jp),maxval(isize)*knpts*kind(inind0)/4,mpi_integer,&
                    inind_all,maxval(isize)*knpts*kind(inind_all)/4,mpi_integer,&
                    jp,ncomm,infompi)
#endif
  ENDDO
  !
ELSE
  inind_all(:,:,:) = inind0(:,:,:)
ENDIF
!
DEALLOCATE(inind0)
!
IF (lhook) CALL dr_hook('INTERPOL_NPTS_4',1,zhook_handle)
IF (lhook) CALL dr_hook('INTERPOL_NPTS_5',0,zhook_handle)
!
!zfield contains the values of the points needed located in this task
!(ie values for indexes of ININD_ALL)
ALLOCATE(zfield(maxval(isize),knpts,SIZE(pfield,2),0:nproc-1))
zfield(:,:,:,:) = xundef
DO jp=0,nproc-1
  DO jk=1,maxval(isize)
    DO jl=1,knpts
      IF (inind_all(jk,jl,jp)/=0) THEN
        !pfield in only on this task
        zfield(jk,jl,:,jp) = pfield(inind_all(jk,jl,jp),:)
      ENDIF
    ENDDO
  ENDDO
ENDDO
!
DEALLOCATE(inind_all)
!
!ZFIELD2 gathers values needed for this task from all other tasks
!(inverse operation than before)
ALLOCATE(zfield2(icpt,knpts,SIZE(pfield,2),0:nproc-1))
IF (nproc>1) THEN
  DO jp=0,nproc-1
#if defined(SFX_MPI) || defined(SFX_MNH)
    CALL mpi_gather(zfield(1:isize(jp),:,:,jp),SIZE(zfield(1:isize(jp),:,:,jp))*kind(zfield)/4,mpi_real,&
                    zfield2,isize(jp)*knpts*SIZE(pfield,2)*kind(zfield2)/4,mpi_real,jp,ncomm,infompi)
#endif
  ENDDO
ELSE
  zfield2(:,:,:,:) = zfield(:,:,:,:)
ENDIF
!
DEALLOCATE(zfield)
!
!zfield3 contains the values of the points needed for interpolation, gathered from all tasks
ALLOCATE(zfield3(icpt,knpts,SIZE(pfield,2)))
DO jp=0,nproc-1
  WHERE (zfield2(:,:,:,jp)/=xundef) zfield3(:,:,:) = zfield2(:,:,:,jp)
ENDDO
DEALLOCATE(zfield2)
!
IF (lhook) CALL dr_hook('INTERPOL_NPTS_5',1,zhook_handle)
IF (lhook) CALL dr_hook('INTERPOL_NPTS_6',0,zhook_handle)
!
!values of the points used for interpolation
ALLOCATE(znval(ip,il2))
znval(:,:) = xundef
!
!znval contains the averaged values for the knpts points
znval(:,:) = 0.
DO jl=1,icpt
  zsum = 0.
  DO jp=1,knpts
    IF (inind(jl,jp)/=0) THEN
      znval(jl,:) = znval(jl,:) + zfield3(jl,jp,:)/zndist(jl,jp)
      zsum = zsum + 1./zndist(jl,jp)
    ENDIF
  ENDDO
  IF (zsum/=0.) znval(jl,:) = znval(jl,:) / zsum
ENDDO
!
DEALLOCATE(inind, zndist, zfield3)
!
!
!finally, pfield contains the interpolated values! 
icpt=0
DO jl=1,il1

  IF (kcode(jl)/=0) cycle

  icpt = icpt + 1
  pfield(jl,:) = znval(icpt,:)
  
ENDDO
!
DEALLOCATE(znval)
!
DEALLOCATE(isize)
!
IF (iold==2) THEN
#ifdef SFX_MNH
  nrank = irank_save
  nproc = iproc_save
  npio = ipio_save
  ncomm = icomm_save
  DEALLOCATE(iindex_tot,inum_tot)
#endif
ENDIF
!
IF (lhook) CALL dr_hook('INTERPOL_NPTS_6',1,zhook_handle)
!-------------------------------------------------------------------------------
!
END SUBROUTINE interpol_npts