horibl_surf_coef.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_coef(KOLEN,OINTERP,OGLOBLON,PILO1,PILO2,POLO,&
                                KO,KINLO,KP,PLOP  )  
!   ###########################################################################
!
!!****  *HORIBL_SURF_COEF* - 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_surf_par,  ONLY : xundef
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE modi_abor1_sfx
!
IMPLICIT NONE
!
!*      0.1. Declaration of arguments
!         
INTEGER,                   INTENT(IN)  :: KOLEN   ! size of output array
LOGICAL, DIMENSION(:), INTENT(IN)  :: OINTERP ! .true. where physical value is needed
!
LOGICAL, INTENT(IN) :: OGLOBLON  ! True if the map is circular
REAL, INTENT(IN) :: PILO1     ! Longitude of the first data point
REAL, INTENT(IN) :: PILO2     ! Longitude of the last data point
REAL, DIMENSION(:), INTENT(IN) :: POLO
INTEGER, DIMENSION(:,:), INTENT(IN) :: KO
!
INTEGER, DIMENSION(:), INTENT(IN) :: KINLO     ! Extended KINLO
!
! Posiiton in the array of the twelwe surrounding points
INTEGER, DIMENSION(:,:), INTENT(OUT) :: KP  
  ! Latitudes and longitudes of the surrounding points
REAL, DIMENSION(:,:), INTENT(OUT) :: PLOP  
!
!*      0.2. Declaration of local variables
!            
 ! Variables used to perform the interpolation
REAL                               :: ZIDLO    ! Delta longitude
INTEGER, DIMENSION(:), ALLOCATABLE :: IOFS     ! Offset of each parallel in the array
!
INTEGER :: IINLA
INTEGER :: JL, IS1, IS2, JI
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!------------------------------------------------------------------------------
!
!*     1. DETERMINATION  of the latitude of the poles (depending of the latitude
!         -------------                                 of the first data point)
!
IF (lhook) CALL dr_hook('HORIBL_SURF_COEF',0,zhook_handle)
!
!------------------------------------------------------------------------------
!
!*     3.   LOOP OVER ALL THE POINTS OF THE OUTPUT GRID
!           -------------------------------------------
!
iinla = SIZE(kinlo)
ALLOCATE (iofs(iinla))
iofs(1) = 1
IF (ogloblon) iofs(1)=iofs(1)+2
DO jl = 2,iinla
  iofs(jl) = iofs(jl-1) + kinlo(jl-1)
  IF (ogloblon) iofs(jl) = iofs(jl) + 4
END DO
!
DO jl = 1, kolen
  !
  IF (.NOT. ointerp(jl)) cycle
  !
  ! 3.1.2. northern
  zidlo = (pilo2 - pilo1) / (kinlo(ko(jl,1)))
  kp(jl,1) = int((polo(jl) - pilo1) / zidlo)
  kp(jl,2) = kp(jl,1) + 1
  IF (.NOT.ogloblon) kp(jl,2) = min(kinlo(ko(jl,1))-1, kp(jl,2))
  plop(jl,1) = pilo1 + kp(jl,1) * zidlo
  plop(jl,2) = plop(jl,1) + zidlo
!
      ! 3.1.3. north
  zidlo = (pilo2 - pilo1) / (kinlo(ko(jl,2) ))
  kp(jl,4)   = int((polo(jl) - pilo1) / zidlo)
  kp(jl,3) = kp(jl,4) - 1
  kp(jl,5) = kp(jl,4) + 1
  kp(jl,6) = kp(jl,4) + 2
  IF (.NOT.ogloblon) THEN
    kp(jl,3) = max(0,kp(jl,3))
    kp(jl,5:6) = min(kinlo(ko(jl,2))-1,kp(jl,5:6))
  ENDIF
  plop(jl,4) = pilo1 + kp(jl,4) * zidlo
  plop(jl,3) = plop(jl,4) - zidlo
  plop(jl,5) = plop(jl,4) + zidlo
  plop(jl,6) = plop(jl,4) + 2*zidlo
!
      ! 3.1.4. south
  zidlo = (pilo2 - pilo1) / (kinlo(ko(jl,3) ))
  kp(jl,8) = int((polo(jl) - pilo1) / zidlo)
  kp(jl,7)  = kp(jl,8) - 1
  kp(jl,9)  = kp(jl,8) + 1
  kp(jl,10) = kp(jl,8) + 2
  IF (.NOT.ogloblon) THEN
    kp(jl,7)  = max(0,kp(jl,7))
    kp(jl,9:10) = min(kinlo(ko(jl,3))-1,kp(jl,9:10))
  ENDIF
  plop(jl,8)  = pilo1 + kp(jl,8) * zidlo
  plop(jl,7)  = plop(jl,8) - zidlo
  plop(jl,9)  = plop(jl,8) + zidlo
  plop(jl,10) = plop(jl,8) + 2*zidlo
!
      ! 3.1.5. southern
  zidlo = (pilo2 - pilo1) / (kinlo(ko(jl,4)))
  kp(jl,11) = int((polo(jl) - pilo1) / zidlo) 
  kp(jl,12) = kp(jl,11) + 1 
  IF (.NOT.ogloblon) kp(jl,12) = min(kinlo(ko(jl,4))-1,kp(jl,12))
  plop(jl,11) = pilo1 + kp(jl,11) * zidlo
  plop(jl,12) = plop(jl,11) + zidlo
!
      ! 3.1.6. check position of points
  IF (ogloblon) THEN 
    is1 = -2
    is2 = 1
  ELSE
    is1 = 0
    is2 = -1
  ENDIF
  !
  IF ((kp(jl,1) <is1) .OR. (kp(jl,2) >kinlo(ko(jl,1))+is2) .OR. &
      (kp(jl,3) <is1) .OR. (kp(jl,6) >kinlo(ko(jl,2))+is2) .OR. &
      (kp(jl,7) <is1) .OR. (kp(jl,10)>kinlo(ko(jl,3))+is2) .OR. &
      (kp(jl,11)<is1) .OR. (kp(jl,12)>kinlo(ko(jl,4))+is2)) THEN  
    CALL abor1_sfx('HORIBLE_SURF: INPUT DOMAIN SMALLER THAN OUTPUT ONE - LONGITUDE GLOBAL')
  END IF
!
      ! 3.1.7. add parallel offset
  kp(jl,1:2)  = kp(jl,1:2) + iofs(ko(jl,1))
  kp(jl,3:6)  = kp(jl,3:6) + iofs(ko(jl,2))
  kp(jl,7:10) = kp(jl,7:10) + iofs(ko(jl,3))
  kp(jl,11:12)= kp(jl,11:12) + iofs(ko(jl,4))
!
END DO
!
!
IF (lhook) CALL dr_hook('HORIBL_SURF_COEF',1,zhook_handle)
!
!------------------------------------------------------------------------------
!
!
END SUBROUTINE horibl_surf_coef