interpol_field2d.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_field2d (UG, U, &
                                   HPROGRAM,KLUOUT,KCODE,PFIELD,HFIELD,PDEF,KNPTS)
!     ################################################
!
!!**** *INTERPOL_FIELD* initializes coordinate system for spline interpolation
!!
!!    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
!!
!!
!!
!!    AUTHOR
!!    ------
!!
!!    V. Masson          Meteo-France
!!
!!    MODIFICATION
!!    ------------
!!
!!    Original    01/2004
!!    Modification
!!      A. Alias        07/2013 add MODI_ABOR1_SFX
!----------------------------------------------------------------------------
!
!*    0.     DECLARATION
!            -----------
!
USE modd_surfex_mpi, ONLY : nrank, npio
!
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_grid_coord
USE modi_interpol_npts
USE modi_sum_on_all_procs
USE modi_abor1_sfx
USE modi_get_interp_halo
!
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,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(INOUT)  :: PFIELD   ! pgd field on grid mesh.
 CHARACTER(LEN=*),        INTENT(IN)   :: HFIELD   ! name of the field for prints
REAL,DIMENSION(:),OPTIONAL, INTENT(IN):: PDEF     ! default value if not enough data
INTEGER, OPTIONAL,       INTENT(IN)   :: KNPTS    ! number of points to interpolate with

!
!*    0.2    Declaration of local variables
!            ------------------------------
!
REAL, DIMENSION(U%NDIM_FULL)   :: ZX             ! coordinate used for
REAL, DIMENSION(U%NDIM_FULL)   :: ZY             ! splines interpolation
REAL, DIMENSION(SIZE(PFIELD,2)):: ZDEF           ! default value for field
INTEGER                        :: INPTS         ! number of points to interpolate with
INTEGER :: IHALO, INEAR_NBR
!
INTEGER                        :: JLOOP          ! loop counter
!
INTEGER :: IERR0
INTEGER                        :: IERR1          ! number of points interpolated
INTEGER                        :: IERR2          ! number of points not interpolated in the end
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('INTERPOL_FIELD:INTERPOL_FIELD2D',0,zhook_handle)
!
inpts = 3
IF (PRESENT(knpts)) inpts = knpts
!
zdef = xundef
IF (PRESENT(pdef)) zdef = pdef
!
!*    2.     Miscellaneous Initializations
!            -----------------------------
!
 CALL get_grid_coord(ug%G%CGRID, ug%G%NGRID_PAR, ug%G%XGRID_PAR, u%NSIZE_FULL, &
                     kluout,kl=u%NDIM_FULL,hgrid=ug%G%CGRID,pgrid_par=ug%XGRID_FULL_PAR,&
                        px=zx,py=zy)
!
!-------------------------------------------------------------------------------
!
!*    5.     Interpolation with 3 nearest points
!            -----------------------------------
!
ierr0 = sum_on_all_procs(hprogram,ug%G%CGRID,kcode(:)==0)
!
 CALL get_interp_halo(hprogram,ug%G%CGRID,ihalo)
!
IF (ihalo/=0) THEN
  inear_nbr = (2*ihalo+1)**2
ELSE
  inear_nbr = u%NDIM_FULL
ENDIF
!
 CALL interpol_npts(ug, u, &
                    hprogram,kluout,inpts,kcode,zx,zy,pfield,inear_nbr)
!
!-------------------------------------------------------------------------------
!
!*    6.     Final check
!            -----------
!
ierr1 = sum_on_all_procs(hprogram,ug%G%CGRID,kcode(:)==0)
ierr2 = sum_on_all_procs(hprogram,ug%G%CGRID,kcode(:)==-4)
!
IF (nrank==npio) THEN
  !
  IF (ierr1>0 .OR. ierr2>0) THEN
    !
    WRITE(kluout,*) ' '
    WRITE(kluout,*) ' Interpolation of field : ',hfield
    WRITE(kluout,*) ' ----------------------'
    WRITE(kluout,*) ' '
    WRITE(kluout,*) ' Number of points interpolated with ',inpts,' nearest points: ', &
                      ierr1
    !
    !
    IF (ierr2>0) THEN
      WRITE(kluout,*) ' Number of points that could not be interpolated : ', &
                        ierr2
                !if all points were scanned or if no point could be interpolated 
      IF (PRESENT(pdef) .AND. (inear_nbr>=u%NDIM_FULL .OR. ierr2==ierr0)) THEN          
        DO jloop=1,SIZE(pfield,2)
          WRITE(kluout,*) ' For these points, the default value (',zdef(jloop),') is set.'
        ENDDO
      ELSE
        WRITE(kluout,*) ' Please provide data with better resolution'
        WRITE(kluout,*) ' Or define a higher halo value             '
      END IF
    END IF
    !
  END IF
  !
END IF
!
IF (ierr2>0) THEN
  !
  IF (.NOT.PRESENT(pdef) .OR. (inear_nbr<u%NDIM_FULL .AND. ierr2/=ierr0)) &
          CALL abor1_sfx('Some points lack data and are too far away from other points. &
                Please define a higher halo value in NAM_IO_OFFLINE.')
  !
ENDIF
!
IF (count(kcode(:)==-4)>0) THEN
  !
  DO jloop=1,SIZE(pfield,2)
    IF (zdef(jloop)/=xundef) THEN
      WHERE(kcode(:)==-4)
        pfield(:,jloop)=zdef(jloop)
      END WHERE
    ENDIF
  END DO
  !
END IF
!
IF (lhook) CALL dr_hook('INTERPOL_FIELD:INTERPOL_FIELD2D',1,zhook_handle)
!-------------------------------------------------------------------------------
!
END SUBROUTINE interpol_field2d