horibl_surf_gridin.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 horibl_surf_gridin(KINLA,KINLO,KILEN,PARIN,KOLEN,&
                                 ODVECT,KLUOUT,OGLOBS,OGLOBN,OGLOBLON,KP,&
                                 PARIN0_OUT,PARIN_OUT,KLSMIN_OUT,KLSMIN,KLSMOUT,KMASK )  
!   ###########################################################################
!
!!****  *HORIBL_SURF_GRIDIN* - horitontal bilinear interpolation
!!
!!    PURPOSE
!!    -------
!!
!!    Interpolates a field, supports masks.
!!
!!    METHOD
!!    ------
!!
!!    This routine performs a bilinear interpolation based on the 12 surrounding
!!    points. It begins with an interpolation along the latitudes (with third order
!!    polynoms interpolation with 4 points and linear interpolation for 2 points)
!!    and then a second along the longitude (third order polynoms interpolation).
!!    Two interpolations are performed : first along the parallels then between the
!!    four resulting points.
!!
!!    The disposition of the points is the following :
!!
!!
!!            N         1   2
!!
!!            ^     3   4   5   6
!!            |           x
!!            |     7   8   9  10
!!            |
!!                     11  12
!!            S
!!              W ---------------> E
!!
!!   Note : the name 'south', 'north', may not be exact if the last data point is
!!     to the south of first (delta latitude < 0). This does not affect computations.
!!
!!   The formula used to compute the weight is :
!!        (Lon   - Lon.i) . (Lon   - Lon.i) . (Lon   - Lon.i)
!!   Wi = ---------------------------------------------------
!!        (Lon.i - Lon.j) . (Lon.i - Lon.k) . (Lon.i - Lon.l)
!!   Where j,k,l are the other points of the line.
!!
!!   When masks are used, points with different types than the output points are
!!   not taken in account (in the formula, the corresponding coefficient is set
!!   to 1). If no points of the same nature are available, the interpolation is
!!   performed anyway with the 12 points. It is the task of the calling program
!!   to react to this situation.
!!
!!   When the inputs parameters define a circular map (or global), the inputs data
!!   are extended. The value of the parameter ODVECT is used to know if the datas
!!   are vectorial or scalar (this affects the sign of extended values).
!!
!!   EXTERNAL
!!   --------
!!
!!   subroutine FMLOOK_ll : to retrieve the logical unit number of the listing file
!!
!!   IMPLICIT ARGUMENTS
!!   ------------------
!!
!!   REFERENCE
!!   ---------
!!
!!   This routine is based on the one used by the software FULL-POS from Meteo France.
!!   More informations may be found in 'Book 1'
!!
!!   AUTHOR
!!   ------
!!
!!   J.Pettre & V.Bousquet
!!
!!   MODIFICATIONS
!!   -------------
!!
!!   Original       07/01/1999
!!                  21/04/1999 (V. Masson) set correct prefixes and bug in
!!                             a logical definition
!!                  21/04/1999 (V. Masson) bug in north and south poles
!!                             extension for input map land-sea mask
!!                  27/05/1999 (V. Masson) bug in 'grib south pole'
!!                             extrapolation (number of point per parallel)
!!                  27/05/1999 (V. Masson) bug in 'grib pole' extrapolation
!!                             extra latitudes are now computed symetrically
!!                             to the poles.
!!                  17/03/2010 (P. LeMoigne) bug in weights computations
!!                  16/06/2010 (G. Tanguy) bug in 'grib north pole"
!!                              extrapolation (tabular ZARIN not totaly filled)
!!
!------------------------------------------------------------------------------
!
!*      0. DECLARATIONS
!       ---------------
!
USE modd_surfex_mpi, ONLY : nrank, npio, ncomm, nproc, idx_i
USE modd_surf_par,  ONLY : xundef
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE modi_abor1_sfx
!
IMPLICIT NONE
!
#ifdef SFX_MPI
include "mpif.h"
#endif
!
!*      0.1. Declaration of arguments
!         
INTEGER,                   INTENT(IN)  :: KINLA   ! Number of parallels
INTEGER, DIMENSION(:), INTENT(IN)  :: KINLO   ! Number of point along a parallel
INTEGER,                   INTENT(IN)  :: KILEN   ! size of input arrays
REAL,    DIMENSION(:,:), INTENT(IN)  :: PARIN   ! input array
INTEGER,                   INTENT(IN)  :: KOLEN   ! size of output array
LOGICAL,                   INTENT(IN)  :: ODVECT  ! data is vectorial (True/False)
INTEGER,                   INTENT(IN)  :: KLUOUT  ! output listing logical unit
LOGICAL,                   INTENT(IN)  :: OGLOBS
LOGICAL,                   INTENT(IN)  :: OGLOBN
LOGICAL,                   INTENT(IN)  :: OGLOBLON
INTEGER, DIMENSION(:,:), INTENT(IN) :: KP
REAL,    DIMENSION(:,:), POINTER  :: PARIN0_OUT   ! input array
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PARIN_OUT
INTEGER, DIMENSION(:,:,:), INTENT(OUT) :: KLSMIN_OUT
INTEGER, DIMENSION(:), POINTER, OPTIONAL :: KMASK
INTEGER, DIMENSION(:,:), INTENT(IN), OPTIONAL  :: KLSMIN  ! input land/sea mask
INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL  :: KLSMOUT ! output land/sea mask
!
!*      0.2. Declaration of local variables
!            
REAL,   DIMENSION(:,:), ALLOCATABLE :: ZARIN     ! Extended input datas
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZOUT
REAL :: ZVECT     ! -1 if input is vectorial 
 ! Variables implied in the extension procedure
INTEGER, DIMENSION(:,:), ALLOCATABLE :: ILSMIN    ! Extended land/sea mask
INTEGER, DIMENSION(:,:), ALLOCATABLE :: IP
INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: IOUT
#ifdef SFX_MPI
INTEGER, DIMENSION(MPI_STATUS_SIZE) :: ISTATUS
#endif
REAL, DIMENSION(:,:), ALLOCATABLE :: ZDARIN
INTEGER, DIMENSION(:,:), ALLOCATABLE :: IDLS
INTEGER, DIMENSION(2) :: IEXT
INTEGER :: IBIGSIZE  ! Size of the extended map
INTEGER :: IMID   ! Used for extensions around the poles
INTEGER :: IOF1  ! Offset in map
INTEGER :: IOF2  ! Offset in map
 ! Loop counters
INTEGER :: JOP     ! Output position
INTEGER :: JIP     ! Input position
INTEGER :: JL, JI, JI2, JL2, JT, INL, JC   ! Dummy counter
INTEGER :: ID1, ID2, ISIZE
INTEGER :: INFOMPI, I, J
INTEGER :: IKINLO
LOGICAL :: LDLSM     ! Specify if land/sea mask is present or not
!
!------------------------------------------------------------------------------
REAL(KIND=JPRB) :: ZHOOK_HANDLE, ZHOOK_HANDLE_OMP
!------------------------------------------------------------------------------
!
!*     1. DETERMINATION  of the latitude of the poles (depending of the latitude
!         -------------                                 of the first data point)
!
IF (lhook) CALL dr_hook('HORIBL_SURF_GRIDIN_1',0,zhook_handle)
!
IF (PRESENT(kmask)) NULLIFY(kmask)
!
inl = SIZE(parin_out,2)
!
ldlsm = .false.
IF (PRESENT(klsmin) .AND. PRESENT(klsmout)) ldlsm = .true.
!
!* 2.6 Compute the resulting map
!
zvect = 1.
IF (odvect) zvect=-1.
!
IF (nrank==npio) THEN
       
  ibigsize = kilen
  IF (oglobs  ) ibigsize=ibigsize+(4+kinlo(    1))+(4+kinlo(      2))
  IF (oglobn  ) ibigsize=ibigsize+(4+kinlo(kinla))+(4+kinlo(kinla-1))
  IF (ogloblon) ibigsize=ibigsize+ 4*kinla
  !
  ALLOCATE (zarin(ibigsize,inl))
  ALLOCATE (ilsmin(ibigsize,inl))
  ALLOCATE (kmask(ibigsize))
  kmask(:) = -1
!
! 2.6.1 Compute the longitude extension
!
! This is a basic copy of the data. If extension is possible, the first and last
! two lines are copied twice this way :
!
!    /---------------\
!    |               |
!   [.] [.] [....   ...] [.] [.] 
!        |            |
!        \------------/
!
! A point represent a data.
!

  DO jt = 1,inl

    jip = 1
    jop = 1
    IF (oglobs) jop = jop + (4+kinlo(1)) + (4+kinlo(2))

    IF (ogloblon) THEN

      DO jl = 1,kinla
        id1 = jip+kinlo(jl)
        id2 = jop+2+kinlo(jl)

        kmask(jop  ) = id1-2
        kmask(jop+1) = id1-1
        DO ji = jop+2,id2-1
          kmask(ji) = jip + ji - (jop+2)
        ENDDO
        kmask(id2  ) = jip
        kmask(id2+1) = jip+1
        
        zarin(jop  ,jt) = parin(id1-2,jt)
        zarin(jop+1,jt) = parin(id1-1,jt)
        zarin(jop+2:id2-1,jt) = parin(jip:id1-1,jt)
        zarin(id2  ,jt) = parin(jip  ,jt)
        zarin(id2+1,jt) = parin(jip+1,jt)

        IF (ldlsm) THEN
          ilsmin(jop  ,jt) = klsmin(id1-2,jt)
          ilsmin(jop+1,jt) = klsmin(id1-1,jt)
          ilsmin(jop+2:id2-1,jt) = klsmin(jip:id1-1,jt)
          ilsmin(id2  ,jt) = klsmin(jip  ,jt)
          ilsmin(id2+1,jt) = klsmin(jip+1,jt)
        END IF

        jip = jip + kinlo(jl)
        jop = jop + kinlo(jl) + 4
        
      END DO

    ELSE

      DO ji = jop,jop+kilen-1
        kmask(ji) = jip + ji - jop
      ENDDO

      zarin(jop:jop+kilen-1,jt) = parin(jip:jip+kilen-1,jt)
      IF (ldlsm) THEN
        ilsmin(jop:jop+kilen-1,jt) = klsmin(jip:jip+kilen-1,jt)
      END IF
    END IF
!
! 2.6.2 Compute the south pole extension
!
! Pole extension is performed by copying the first half datas to the last half 
! datas of the extension parallel :
!
!  [.] [.] [....] [....] [.] [.]
!                  ||||
!            /-------/
!           ||||
!  [.] [.] [....] [....] [.] [.]
!
    DO jc = 1,4
      !
      IF (jc<3.AND.oglobs.OR.jc>2.AND.oglobn) THEN
        !
        IF (jc==1) THEN
          iof1 = 4 + kinlo(2)
          iof2 = iof1 + 4 + kinlo(1)
          imid = (kinlo(1)+4) / 2
          ikinlo = kinlo(1) 
        ELSEIF (jc==2) THEN
          iof1 = 0
          iof2 = iof2 + 4 + kinlo(2)
          imid = (kinlo(2)+4) / 2
          ikinlo = kinlo(2)
        ELSEIF (jc==3) THEN
          iof1 = ibigsize - (4+kinlo(kinla-1)) - (4+kinlo(kinla))
          iof2 = iof1 - (4+kinlo(kinla))
          imid = (kinlo(kinla)+4) / 2
          ikinlo =  kinlo(kinla)           
        ELSEIF (jc==4) THEN
          iof1 = iof1 + (4+kinlo(kinla))
          iof2 = iof2 - (4+kinlo(kinla-1))
          imid = (kinlo(kinla-1)+4) / 2
          ikinlo = kinlo(kinla-1)
        ENDIF
        !
        DO ji = 1,imid
          kmask(iof1+ji) = iof2 + imid + ji - 2
        ENDDO
        DO ji = imid+1,ikinlo+4
          kmask(iof1+ji) = iof2 + 1 + 2 + ji - (imid + 1)
        ENDDO
        !
        zarin(iof1+1:iof1+imid,jt) = zvect*zarin(iof2+1+imid-2:iof2+2*imid-2,jt)  
        zarin(iof1+imid+1:iof1+ikinlo+4,jt) = zvect*zarin(iof2+1+2:iof2+ikinlo+4-imid+2,jt)  
        IF (ldlsm) THEN
          ilsmin(iof1+1:iof1+imid,jt) = ilsmin(iof2+1+imid-2:iof2+2*imid-2,jt)  
          ilsmin(iof1+imid+1:iof1+ikinlo+4,jt) = ilsmin(iof2+1+2:iof2+ikinlo+4-imid+2,jt)  
        END IF
        !
      END IF
      !
    ENDDO
    !
  ENDDO
  !
ENDIF
!
IF (nproc>1) THEN
#ifdef SFX_MPI
  CALL mpi_bcast(ibigsize,kind(ibigsize)/4,mpi_integer,npio,ncomm,infompi)
  IF (nrank/=npio) ALLOCATE(kmask(ibigsize))
  CALL mpi_bcast(kmask,SIZE(kmask)*kind(kmask)/4,mpi_integer,npio,ncomm,infompi)
#endif
ENDIF
!
IF (nrank==npio .AND. (ogloblon.OR.oglobn.OR.oglobs)) THEN
  !
  ALLOCATE(parin0_out(SIZE(zarin,1),SIZE(zarin,2)))
  parin0_out = zarin
  !
ELSE
  !
  ALLOCATE(parin0_out(0,0))
  !
ENDIF
!
IF (lhook) CALL dr_hook('HORIBL_SURF_GRIDIN_1',1,zhook_handle)

IF (nrank/=npio) THEN
IF (lhook) CALL dr_hook('HORIBL_SURF_GRIDIN_3',0,zhook_handle)
  iext(1) = minval(kp)
  iext(2) = maxval(kp)
  idx_i = idx_i + 1
#ifdef SFX_MPI  
  CALL mpi_send(iext,2*kind(iext)/4,mpi_integer,npio,idx_i,ncomm,infompi)
#endif  
    IF (lhook) CALL dr_hook('HORIBL_SURF_GRIDIN_3',1,zhook_handle)
    IF (lhook) CALL dr_hook('HORIBL_SURF_GRIDIN_4',0,zhook_handle)
    isize = iext(2)-iext(1)+1
    ALLOCATE(zdarin(isize,inl))
    ALLOCATE(idls(isize,inl))
    idx_i = idx_i + 1
#ifdef SFX_MPI    
  CALL mpi_recv(zdarin,SIZE(zdarin)*kind(zdarin)/4,mpi_real,npio,idx_i,ncomm,istatus,infompi)
#endif  
  idx_i = idx_i + 1
#ifdef SFX_MPI  
  CALL mpi_recv(idls,SIZE(idls)*kind(idls)/4,mpi_integer,npio,idx_i,ncomm,istatus,infompi)
#endif  
  IF (lhook) CALL dr_hook('HORIBL_SURF_GRIDIN_4',1,zhook_handle)
ELSE
#ifdef SFX_MPI
!$OMP PARALLEL DO SCHEDULE(DYNAMIC,1) PRIVATE(J,IEXT,ISIZE,ZDARIN,IDLS,JT,JL,ISTATUS,INFOMPI,ZHOOK_HANDLE_OMP)
#endif
  DO j=0,nproc-1
    IF (lhook) CALL dr_hook('HORIBL_SURF_GRIDIN_2',0,zhook_handle_omp)
    IF (j/=npio) THEN
#ifdef SFX_MPI            
      CALL mpi_recv(iext,2*kind(iext)/4,mpi_integer,j,idx_i+1,ncomm,istatus,infompi)
#endif      
    ELSE
      iext(1) = minval(kp)
      iext(2) = maxval(kp)
    ENDIF
    IF (lhook) CALL dr_hook('HORIBL_SURF_GRIDIN_2',1,zhook_handle_omp)
    IF (lhook) CALL dr_hook('HORIBL_SURF_GRIDIN_30',0,zhook_handle_omp)
    isize = iext(2)-iext(1)+1
    ALLOCATE(zdarin(isize,inl))
    ALLOCATE(idls(isize,inl))
    DO jt = 1,inl
      DO jl = iext(1),iext(2)
        zdarin(jl-iext(1)+1,jt) = zarin(jl,jt)
        idls(jl-iext(1)+1,jt) = ilsmin(jl,jt)
      ENDDO
    ENDDO
    IF (lhook) CALL dr_hook('HORIBL_SURF_GRIDIN_30',1,zhook_handle_omp)
    IF (lhook) CALL dr_hook('HORIBL_SURF_GRIDIN_40',0,zhook_handle_omp)
    IF (j/=npio) THEN
#ifdef SFX_MPI            
      CALL mpi_send(zdarin,SIZE(zdarin)*kind(zdarin)/4,mpi_real,j,idx_i+2,ncomm,infompi)
      CALL mpi_send(idls,SIZE(idls)*kind(idls)/4,mpi_integer,j,idx_i+3,ncomm,infompi)
#endif
    ELSE
      DO jt = 1,inl
        DO jl2 = 1,12
          DO jl = 1,SIZE(kp,1)
            id1 = kp(jl,jl2) - iext(1) + 1
            parin_out(jl,jt,jl2) = zdarin(id1,jt)
            klsmin_out(jl,jt,jl2) = idls(id1,jt)
          ENDDO
        ENDDO
      ENDDO
    ENDIF
    DEALLOCATE(zdarin,idls)
    IF (lhook) CALL dr_hook('HORIBL_SURF_GRIDIN_40',1,zhook_handle_omp)
  ENDDO
#ifdef SFX_MPI
!$OMP END PARALLEL DO
#endif
  idx_i = idx_i+3
ENDIF
!
IF (lhook) CALL dr_hook('HORIBL_SURF_GRIDIN_5',0,zhook_handle)
IF (nrank/=npio) THEN
DO jt = 1,inl
  DO jl2 = 1,12
    DO jl = 1,SIZE(kp,1)
      id1 = kp(jl,jl2) - iext(1) + 1
      parin_out(jl,jt,jl2) = zdarin(id1,jt)
      klsmin_out(jl,jt,jl2) = idls(id1,jt)
    ENDDO
  ENDDO
ENDDO
ENDIF
!
IF (lhook) CALL dr_hook('HORIBL_SURF_GRIDIN_5',1,zhook_handle)
!
END SUBROUTINE horibl_surf_gridin