hor_extrapol_surf.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 hor_extrapol_surf(KLUOUT,HCOORTYPE,KILEN,PILA1,PILA2,PILO1,PILO2,&
                                   KINLA,KINLO,KP,PFIELD_IN,PLAT,PLON,PFIELD,OINTERP,&
                                   PILATARRAY)
!     ###################################################################
!
!!**** *HOR_EXTRAPOL_SURF* extrapolate a surface field
!!
!!    PURPOSE
!!    -------
!!
!!    METHOD
!!    ------
!!       For each point to interpolate, the nearest valid point value is set.
!!
!!    EXTERNAL
!!    --------
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    REFERENCE
!!    ---------
!!
!!    AUTHOR
!!    ------
!!
!!    V. Masson          Meteo-France
!!
!!    MODIFICATION
!!    ------------
!!
!!    Original     01/12/98
!!     V. Masson    01/2004 extrapolation in latitude and longitude
!!     M. Jidane    11/2013 add OpenMP directives
!----------------------------------------------------------------------------
!
!*    0.     DECLARATION
!            -----------
!
USE modd_surfex_mpi, ONLY : nrank, nproc, npio, ncomm, idx_i
USE modd_surf_par,   ONLY : xundef
USE modd_csts,       ONLY : xpi
USE modn_prep_surf_atm, ONLY : nhalo_prep
!
USE modi_abor1_sfx
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
#ifdef SFX_MPI
include "mpif.h"
#endif
!
!*    0.1    Declaration of arguments
!            ------------------------
!
INTEGER,               INTENT(IN)     :: KLUOUT   ! output listing logical unit
 CHARACTER(LEN=4),      INTENT(IN)     :: HCOORTYPE! type of coordinate
 INTEGER, INTENT(IN) :: KILEN
REAL, INTENT(IN) :: PILA1
REAL, INTENT(IN) :: PILA2
REAL, INTENT(IN) :: PILO1
REAL, INTENT(IN) :: PILO2
INTEGER, INTENT(IN) :: KINLA
INTEGER, DIMENSION(:), INTENT(IN) :: KINLO
INTEGER, DIMENSION(:,:), INTENT(IN) :: KP
REAL,   DIMENSION(:,:),  INTENT(IN)     :: PFIELD_IN! input field on grid mesh
REAL,   DIMENSION(:),  INTENT(IN)     :: PLAT     ! latitude of each grid mesh.
REAL,   DIMENSION(:),  INTENT(IN)     :: PLON     ! longitude of each grid mesh.
REAL,   DIMENSION(:,:),  INTENT(INOUT)  :: PFIELD   ! field on grid mesh
LOGICAL,DIMENSION(:),  INTENT(IN)     :: OINTERP  ! .true. where physical value is needed
REAL, DIMENSION(:), INTENT(IN), OPTIONAL :: PILATARRAY
!
!*    0.2    Declaration of local variables
!            ------------------------------
!
REAL, DIMENSION(:), ALLOCATABLE :: ZTLONMIN, ZTLONMAX, ZTLATMIN, ZTLATMAX
REAL, DIMENSION(:,:), ALLOCATABLE :: ZFIELD
REAL :: ZLAT  ! latitude of point to define
REAL :: ZLON  ! longitude of point to define
REAL :: ZDIST ! current distance to valid point (in lat/lon grid)
REAL, DIMENSION(:), ALLOCATABLE :: ZNDIST! smallest distance to valid point
REAL :: ZCOSLA! cosine of latitude
REAL :: ZLONSC! longitude of valid point
REAL :: ZIDLO, ZIDLOMAX, ZIDLOMIN, ZIDLAMAX, ZIDLAMIN
REAL, DIMENSION(:,:), ALLOCATABLE :: ZCOOR
REAL, DIMENSION(:), ALLOCATABLE :: ZIDLA
REAL,    DIMENSION(:), ALLOCATABLE :: ZLA       ! input "latitude"  coordinate
REAL,    DIMENSION(:), ALLOCATABLE :: ZLO       ! input "longitude" coordinate
REAL(KIND=JPRB) :: ZRAD ! conversion degrees to radians
!
INTEGER, DIMENSION(:), ALLOCATABLE :: IMASK, IMASKR
INTEGER, DIMENSION(:,:), ALLOCATABLE :: IVAL_EXT
INTEGER, DIMENSION(NPROC) :: INO_TAB
INTEGER  :: INO     ! output array size
INTEGER, DIMENSION(2) :: ITSIZE, ITDIM
INTEGER, DIMENSION(2,0:NPROC-1) :: IBOR
INTEGER :: ISIZE, ISIZE_MAX, J, ID0, ICOMPT, ICPT
INTEGER :: INFOMPI, IDX, INL
INTEGER  :: JI, JL, JLAT, JLON, JIPOS, JP   ! loop index on points
INTEGER  :: JISC  ! loop index on valid points
#ifdef SFX_MPI
INTEGER, DIMENSION(MPI_STATUS_SIZE) :: ISTATUS
#endif
LOGICAL  :: GLALO ! flag true is second coordinate is a longitude or pseudo-lon.
                  !      false if metric coordinates
REAL(KIND=JPRB) :: ZHOOK_HANDLE, ZHOOK_HANDLE_OMP
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('HOR_EXTRAPOL_SURF_1',0,zhook_handle)
!
ino = SIZE(pfield,1)
inl = SIZE(pfield,2)
!
ALLOCATE(zndist(inl))
!
!-------------------------------------------------------------------------------
!
glalo = hcoortype=='LALO'
!
!-------------------------------------------------------------------------------
!
ALLOCATE(zidla(kinla))
!
ALLOCATE(zla(kilen))
ALLOCATE(zlo(kilen))
!
zidlomax = 0.
zidlomin = xundef
jipos = 0
IF (PRESENT(pilatarray)) THEN
  zidla(1) = 0.
  DO jlat=2,kinla
    zidla(jlat) = pilatarray(jlat) - pilatarray(jlat-1)
  ENDDO
ELSE
  zidla(:) = (pila2 - pila1) / (kinla - 1)
ENDIF
!
zidlamax = maxval(abs(zidla))
zidlamin = minval(abs(zidla(2:kinla)))
!
DO jlat=1,kinla
  IF (glalo) THEN
    zidlo = (pilo2-pilo1) / kinlo(jlat)
  ELSE
    zidlo = (pilo2-pilo1) / (kinlo(jlat)-1)
  ENDIF
  DO jlon=1,kinlo(jlat)
    jipos = jipos + 1
    zla(jipos) = pila1 + sum(zidla(2:jlat))
    zlo(jipos) = pilo1 + (jlon-1) * zidlo 
  END DO
  zidlo = abs(zidlo)
  IF (zidlo>zidlomax) zidlomax = zidlo
  IF (zidlo<zidlomin) zidlomin = zidlo
END DO
!
!-------------------------------------------------------------------------------
!
!*      4.   Loop on points to define
!            ------------------------
!
ALLOCATE(ztlonmin(ino),ztlonmax(ino),ztlatmin(ino),ztlatmax(ino))
ztlonmin(:) = 0.
ztlonmax(:) = 0.
ztlatmin(:) = 0.
ztlatmax(:) = 0.
!
zrad=xpi/180.0_jprb
!
!1: ZTLONMIN, ZTLONMAX, ZTLATMIN, ZTLATMAX contain for each point to extrapol 
! the limits of the domain where to search for the valid points, according
! to NHALO_PREP
icpt = 0
isize_max = 0
isize = 0
DO ji=1,ino
  !
  IF (all(pfield(ji,:)/=xundef)) cycle
  IF (.NOT. ointerp(ji))  cycle
  icpt = icpt + 1
  !
  IF (nhalo_prep>0) THEN
    ztlonmin(ji) = minval(zlo(kp(ji,:)))-zidlomax*nhalo_prep
    ztlonmax(ji) = maxval(zlo(kp(ji,:)))+zidlomax*nhalo_prep
    ztlatmin(ji) = minval(zla(kp(ji,:)))-zidlamax*nhalo_prep
    ztlatmax(ji) = maxval(zla(kp(ji,:)))+zidlamax*nhalo_prep
  ELSE
    ztlonmin(ji) = minval(zlo(:))
    ztlonmax(ji) = maxval(zlo(:))
    ztlatmin(ji) = minval(zla(:))
    ztlatmax(ji) = maxval(zla(:))
  ENDIF
  isize = ceiling((ztlonmax(ji)-ztlonmin(ji)+1)/zidlomin)*&
          ceiling((ztlatmax(ji)-ztlatmin(ji)+1)/zidlamin)
  IF (isize>isize_max) isize_max = isize
  !
ENDDO
!
!number of points to extrapol and maximum number of points to sound
itsize(1) = icpt
itsize(2) = isize_max
!
!NPIO knows the numbers of points to extrapolate for all tasks
IF (nproc>1) THEN
#ifdef SFX_MPI
  CALL mpi_gather(itsize,2*kind(itsize)/4,mpi_integer,&
                  ibor,2*kind(ibor)/4,mpi_integer,& 
                  npio,ncomm,infompi)
#endif
ELSE
  ibor(:,0) = itsize(:)
ENDIF
!
!imask associated the number of the point to extrapolate to its real index in
!the field
ALLOCATE(imask(itsize(1)))
imask(:) = 0
!
IF (nrank==npio) THEN
  ALLOCATE(ival_ext(maxval(ibor(1,:)),maxval(ibor(2,:))))
  ALLOCATE(zcoor(maxval(ibor(1,:)),2))
ELSE
  ALLOCATE(ival_ext(itsize(1),itsize(2)))
  ALLOCATE(zcoor(itsize(1),2))
ENDIF
ival_ext(:,:) = 0
zcoor(:,:) = 0.
!
!
icpt = 0
!2: loop on the points 
DO ji=1,ino
  !
  IF (all(pfield(ji,:)/=xundef)) cycle
  IF (.NOT. ointerp(ji))  cycle
  !
  icpt = icpt + 1
  !imask associated the number of the point to extrapolate to its real index in
  !the field
  imask(icpt) = ji
  !
  icompt = 0
  jisc = 0
  !
  !coordinates of the point in the grid 
  zcoor(icpt,1) = plat(ji)
  zcoor(icpt,2) = plon(ji)
  !
  !loop on the grid
  DO jlat = 1,kinla
    zlat = pila1 + sum(zidla(2:jlat))
    IF (zlat>=ztlatmin(ji) .AND. zlat<=ztlatmax(ji)) THEN
      IF (glalo) THEN
        zidlo = (pilo2-pilo1) / kinlo(jlat)
      ELSE
        zidlo = (pilo2-pilo1) / (kinlo(jlat)-1)
      ENDIF
      IF (jlat>1) jisc = sum(kinlo(1:jlat-1))
      DO jlon = 1,kinlo(jlat)
        zlon = pilo1 + (jlon-1) * zidlo
        IF (zlon>=ztlonmin(ji) .AND. zlon<=ztlonmax(ji)) THEN
          icompt = icompt + 1
          !ival_ext contains the indexes of the points needed to interpolate
          !in the complete grid
          ival_ext(icpt,icompt) = jisc + jlon
        ENDIF
      ENDDO
    ENDIF
  ENDDO
  !
ENDDO
!
DEALLOCATE(zidla)
!
!
IF (lhook) CALL dr_hook('HOR_EXTRAPOL_SURF_1',1,zhook_handle)
!
IF (nrank/=npio) THEN

  IF (lhook) CALL dr_hook('HOR_EXTRAPOL_SURF_2',0,zhook_handle)

  !if some points of this task need to be extrapolated
  IF (sum(itsize)/=0) THEN

    !zfield will contain the values of the field
    ALLOCATE(zfield(itsize(1),inl))

    idx_i = idx_i + 1
    !sends indexes to npio
#ifdef SFX_MPI    
    CALL mpi_send(ival_ext,SIZE(ival_ext)*kind(ival_ext)/4,mpi_integer,npio,idx_i,ncomm,infompi)
#endif

    idx_i = idx_i + 1
    !send coords of the points to extrapolate
#ifdef SFX_MPI    
    CALL mpi_send(zcoor,SIZE(zcoor)*kind(zcoor)/4,mpi_real,npio,idx_i,ncomm,infompi)
#endif

    idx_i = idx_i + 1
    !receives values of the field from NPIO
#ifdef SFX_MPI    
    CALL mpi_recv(zfield,SIZE(zfield)*kind(zfield)/4,mpi_real,npio,idx_i,ncomm,istatus,infompi)
#endif

    DO ji=1,itsize(1)
      pfield(imask(ji),:) = zfield(ji,:)
    ENDDO
    DEALLOCATE(zfield)

  ELSE
    idx_i = idx_i + 3
  ENDIF

  IF (lhook) CALL dr_hook('HOR_EXTRAPOL_SURF_2',1,zhook_handle)

ELSE

IF (lhook) CALL dr_hook('HOR_EXTRAPOL_SURF_31',0,zhook_handle_omp)
  DO jp = npio,nproc-1+npio

    j = jp
    IF (jp>nproc-1) j = jp-nproc

    IF (sum(ibor(:,j))/=0) THEN

      ALLOCATE(zfield(ibor(1,j),inl))
      zfield=xundef

      IF (j/=npio) THEN

        !receives indexes and coordinaites
#ifdef SFX_MPI        
        CALL mpi_recv(ival_ext(1:ibor(1,j),1:ibor(2,j)), ibor(1,j)*ibor(2,j)*kind(ival_ext)/4, &
                        mpi_integer, j, idx_i+1, ncomm, istatus, infompi)
        CALL mpi_recv(zcoor(1:ibor(1,j),:), ibor(1,j)*SIZE(zcoor,2)*kind(zcoor)/4,&
                        mpi_real, j, idx_i+2, ncomm, istatus, infompi)
#endif

      ENDIF
    
!$OMP PARALLEL DO PRIVATE(JI,ZNDIST,IDX,ZCOSLA,JISC,ID0,ZLONSC,ZDIST)
        DO ji=1,ibor(1,j)
          zndist(:) = xundef
          idx = ibor(2,j)+1
          zcosla=cos(zcoor(ji,1)*zrad)
          DO jisc = 1,ibor(2,j)
            !index in the whole grid of the point used to interpolate
            id0 = ival_ext(ji,jisc)
            IF (id0==0) EXIT
            IF (any(pfield_in(id0,:)/=xundef)) THEN
              zlonsc = zlo(id0)
              IF (glalo) THEN
                IF (zlonsc-zcoor(ji,2)> 180.) zlonsc = zlonsc - 360.
                IF (zlonsc-zcoor(ji,2)<-180.) zlonsc = zlonsc + 360.
                zdist= (zla(id0)-zcoor(ji,1)) ** 2 + ((zlonsc-zcoor(ji,2))*zcosla) ** 2
              ELSE
                zdist= (zla(id0)-zcoor(ji,1)) ** 2 + (zlonsc-zcoor(ji,2)) ** 2
              END IF
              DO jl=1,inl
                IF (zdist<=zndist(jl)) THEN
                  IF (pfield_in(id0,jl)/=xundef) THEN
                    zfield(ji,jl) = pfield_in(id0,jl)
                    zndist(jl) = zdist
                  ENDIF
                ENDIF
              ENDDO
            ENDIF
          END DO   
        ENDDO
!$OMP END PARALLEL DO
      !
      IF (j/=npio) THEN
        !send values found to extrapolate
#ifdef SFX_MPI        
        CALL mpi_send(zfield,SIZE(zfield)*kind(zfield)/4,mpi_real,j,idx_i+3,ncomm,infompi)
#endif
      ELSE
        DO ji = 1,ibor(1,j)
          pfield(imask(ji),:) = zfield(ji,:)
        ENDDO
      ENDIF
      !
      DEALLOCATE(zfield)
      !
    ENDIF
    !
  ENDDO
IF (lhook) CALL dr_hook('HOR_EXTRAPOL_SURF_31',1,zhook_handle_omp)
!
DEALLOCATE(zndist)
  !
IF (lhook) CALL dr_hook('HOR_EXTRAPOL_SURF_32',0,zhook_handle)  
  !
  idx_i = idx_i + 3
  !
IF (lhook) CALL dr_hook('HOR_EXTRAPOL_SURF_32',1,zhook_handle)  
  !
ENDIF
!
IF (lhook) CALL dr_hook('HOR_EXTRAPOL_SURF_4',0,zhook_handle)  
!
DEALLOCATE(zcoor)
DEALLOCATE(ival_ext)
DEALLOCATE(imask)
!
IF (ALLOCATED(zla)) DEALLOCATE(zla)
IF (ALLOCATED(zlo)) DEALLOCATE(zlo)
DEALLOCATE(ztlonmin,ztlonmax,ztlatmin,ztlatmax)
!
DO jl=1,inl
  IF (any(pfield(:,jl)==xundef .AND. ointerp(:))) THEN
    WRITE(*,*) 'LAYER ',jl,': NO EXTRAPOLATION : INCREASE YOUR HALO_PREP IN NAM_PREP_SURF_ATM'
    CALL abor1_sfx('NO EXTRAPOLATION : INCREASE YOUR HALO_PREP IN NAM_PREP_SURF_ATM')
  ENDIF
  WHERE (.NOT. ointerp(:)) pfield(:,jl) = xundef
ENDDO
!
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('HOR_EXTRAPOL_SURF_4',1,zhook_handle)
!
END SUBROUTINE hor_extrapol_surf