horibl_surf_init.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_init(PILA1,PILO1,PILA2,PILO2,KINLA,KINLO,KOLEN,&
                           PXOUT,PYOUT,OINTERP,OGLOBLON,OGLOBN,OGLOBS,&
                           KO,KINLO_OUT,POLA,POLO,PILO1_OUT,&
                           PILO2_OUT,PLA,PILATARRAY )  
!   ###########################################################################
!
!!****  *HORIBL_SURF_INIT* - 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
!!
!!   MODIFICATKONS
!!   -------------
!!
!!   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. DECLARATKONS
!       ---------------
!
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
!         
REAL,                      INTENT(IN)  :: PILA1   ! Lat. (y) of first input point KDGSA
REAL,                      INTENT(IN)  :: PILO1   ! Lon. (x) of first input point 
REAL,                      INTENT(IN)  :: PILA2   ! Lat. (y) of last input point KDGEN
REAL,                      INTENT(IN)  :: PILO2   ! Lon. (x) of last input point
INTEGER,                   INTENT(IN)  :: KINLA   ! Number of parallels
INTEGER, DIMENSION(:), INTENT(IN)  :: KINLO   ! Number of point along a parallel
INTEGER,                   INTENT(IN)  :: KOLEN   ! size of output array
REAL,    DIMENSION(:), INTENT(IN)  :: PXOUT   ! X (lon.) of output points
REAL,    DIMENSION(:), INTENT(IN)  :: PYOUT   ! Y (lat.) of output points
LOGICAL, DIMENSION(:), INTENT(IN)  :: OINTERP ! .true. where physical value is needed
! 
LOGICAL, INTENT(OUT)  :: OGLOBLON  ! True if the map is circular
LOGICAL, INTENT(OUT)  :: OGLOBN    ! True if the map has the north pole
LOGICAL, INTENT(OUT)  :: OGLOBS    ! True if the map has the south pole
  ! Number of the surrounding latitudes
INTEGER, DIMENSION(:,:), INTENT(OUT) :: KO
INTEGER, DIMENSION(:), INTENT(OUT) :: KINLO_OUT     ! Extended KINLO
REAL, INTENT(OUT) :: PILO1_OUT     ! Longitude of the first data point
REAL, INTENT(OUT) :: PILO2_OUT     ! Longitude of the last data point
! Variables used to perform the interpolation
REAL, DIMENSION(:), INTENT(OUT) :: POLA     ! Latitude of the output point
REAL, DIMENSION(:), INTENT(OUT) :: POLO     ! Longitude of the output point
  ! Latitudes and longitudes of the surrounding points
REAL, DIMENSION(:,:), INTENT(OUT) :: PLA
REAL,    DIMENSION(:), INTENT(IN), OPTIONAL  :: PILATARRAY! latitudes array
!
!*      0.2. Declaration of local variables
!
REAL, DIMENSION(:), ALLOCATABLE    :: ZIDLAT   ! Deltai latitude
REAL                               :: ZIDLA    ! Delta latitude
REAL                               :: ZSOUTHPOLE! south pole latitude (-90 or  90)
REAL                               :: ZNORTHPOLE! north pole latitude ( 90 or -90)
!
! Variables implied in the extension procedure
!
INTEGER                            :: IOFFSET   ! Offset in map
INTEGER                            :: IINLA     ! Number of parallel
 ! Loop counters
INTEGER                            :: JOPOS     ! Output position
INTEGER                            :: JL, JL2   ! Dummy counter
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!------------------------------------------------------------------------------
!
!*     1. DETERMINATKON  of the latitude of the poles (depending of the latitude
!         -------------                                 of the first data point)
!
IF (lhook) CALL dr_hook('HORIBL_SURF_INIT',0,zhook_handle)
!
IF (pila1>0.) THEN
  zsouthpole= 90.
  znorthpole=-90.
ELSE
  zsouthpole=-90.
  znorthpole= 90.
END IF
!
!------------------------------------------------------------------------------
!
!*     2. EXTEND DATA GRID
!         ----------------
!
!*    2.1 Alias input data
!
pilo1_out = pilo1
pilo2_out = pilo2
!
!*   2.2 Center input domain in order to have Lo1 < Lo 2
!
IF (pilo2_out < 0.)    pilo2_out = pilo2_out + 360.
IF (pilo1_out < 0.)    pilo1_out = pilo1_out + 360.
IF (pilo2_out < pilo1_out) pilo1_out = pilo1_out - 360.
!
!*   2.3 Extend one point (needed for reduced grids)
!
! Longitude coordinate of points are found by :
!                      i
!  Lon(i) = Lon1 + ------------- . (Lon2 - Lon1)
!                   Npts(Lat)-1
! Where i goes from 0 to Npts(Lat)-1. The result of this is that the last point of 
! each parallel is located at Lon2. This is not the case for reduced grid where the 
! position of the last point depends upon the number of points of the parallel. For
! reduced grid, the right formula to use is the following :
!                       i
!  Lon(i) = Lon1 + ----------- . (Lon2' - Lon1)
!                   Npts(Lat)
! Where Lon2' = Lon1 + 2.PI.
!
!                                              Lon2 - Lon1
! This can be generalized with Lon2' = Lon2 + -------------
!                                              Nptsmax - 1
!
jopos = maxval(kinlo(1:kinla))
pilo2_out = pilo1_out + (pilo2_out - pilo1_out) * jopos / (jopos - 1.)
!
!* 2.4 Test if the input is global or partially global
!
! Note that we must have a global map to make extension around the poles
oglobn   = .false.
oglobs   = .false.
ogloblon = .false.
IF (pilo2_out-360.>pilo1_out-1.e-3) ogloblon = .true.
zidla = (pila2 - pila1) / (kinla - 1)
IF (PRESENT(pilatarray)) THEN
  ALLOCATE(zidlat(kinla+1))
  DO jl=2,kinla
    zidlat(jl) = pilatarray(jl)-pilatarray(jl-1)
  END DO
  zidlat(1)=zidlat(2)
  zidlat(kinla+1) = zidlat(kinla)
ENDIF
IF ((pila1-zidla>= 90.) .OR. (pila1-zidla<=-90.)) oglobs=ogloblon
IF ((pila2+zidla>= 90.) .OR. (pila2+zidla<=-90.)) oglobn=ogloblon
! Aladin case (input PILA2, PILO2 are in meters) no extension
IF ( pila2 > 100. ) THEN
  oglobn   = .false.
  oglobs   = .false.
  ogloblon = .false.
END IF
!
!
!*  2.7  Compute the resulting parameters of the map
!
iinla = kinla
IF (oglobs) iinla = iinla + 2
IF (oglobn) iinla = iinla + 2
!
ioffset = 0
IF (oglobs) THEN
  kinlo_out(ioffset+1) = kinlo(2)
  kinlo_out(ioffset+2) = kinlo(1)
  ioffset = ioffset + 2
END IF
kinlo_out(ioffset+1:ioffset+kinla) = kinlo(1:kinla)
ioffset = ioffset + kinla
IF (oglobn) THEN
  kinlo_out(ioffset+1) = kinlo(kinla)
  kinlo_out(ioffset+2) = kinlo(kinla-1)
  ioffset = ioffset + 2
END IF
!
!------------------------------------------------------------------------------
!
!*     3.   LOOP OVER ALL THE POINTS OF THE OUTPUT GRID
!           -------------------------------------------
!
pola(:) = 0.
polo(:) = 0.
!
DO jl = 1, kolen
  !
  IF (.NOT. ointerp(jl)) cycle
  !
  pola(jl)  = pyout(jl)
  !
  polo(jl)  = pxout(jl)
  IF (polo(jl) < pilo1_out) polo(jl) = polo(jl) + 360.
  IF (polo(jl) > pilo2_out) polo(jl) = polo(jl) - 360.
  !
  ! 3.1.1. find positions of latitudes
  IF (PRESENT(pilatarray)) THEN
    !
    DO jl2 = 1,kinla
      IF((pola(jl)>=pilatarray(jl2)-zidlat(jl2)/2..AND.pola(jl)<pilatarray(jl2)+zidlat(jl2+1)/2.).OR.&
         (pola(jl)<=pilatarray(jl2)-zidlat(jl2)/2..AND.pola(jl)>pilatarray(jl2)+zidlat(jl2+1)/2.)) THEN
        ko(jl,3) = jl2
        EXIT
      ELSEIF (pola(jl)>maxval(pilatarray(:))) THEN
        ko(jl,3) = maxloc(pilatarray,1)
      ELSEIF (pola(jl)<minval(pilatarray(:))) THEN
        ko(jl,3) = minloc(pilatarray,1)
      ENDIF
    ENDDO
    pla(jl,3) = pilatarray(ko(jl,3))
    !
  ELSE
    !
    ko(jl,3)  = nint( (pola(jl)-pila1)/zidla - 0.5 ) ! because of the zero
    IF ( ko(jl,3)<-1) CALL abor1_sfx('HORIBLE_SURF_INIT: INPUT DOMAIN SMALLER THAN OUTPUT ONE - LATITUDE')
    !
    pla(jl,3) = pila1 + ko(jl,3) * zidla
    !
    ko(jl,3) = ko(jl,3) + 1
    !
  ENDIF
  !
  IF (PRESENT(pilatarray)) THEN
    !
    IF (ko(jl,3)==1) THEN
      pla(jl,4) = pilatarray(1) - zidlat(1)
    ELSE
      pla(jl,4) = pilatarray(ko(jl,3)-1)
    ENDIF
    IF (ko(jl,3)>=kinla) THEN
      pla(jl,2) = pilatarray(kinla) + zidlat(kinla)
    ELSE
      pla(jl,2) = pilatarray(ko(jl,3)+1)
    ENDIF
    IF (ko(jl,3)>=kinla) THEN
      pla(jl,1) = pilatarray(kinla) + 2.*zidlat(kinla)
    ELSE
      pla(jl,1) = pilatarray(ko(jl,3)+2)
    ENDIF
    !
  ELSE
    !
    pla(jl,1) = pla(jl,3) + 2*zidla
    pla(jl,2) = pla(jl,3) + zidla
    pla(jl,4) = pla(jl,3) - zidla
    !
  ENDIF
  !
  IF (oglobs) ko(jl,3) = ko(jl,3) + 2
  !
  ko(jl,1) = ko(jl,3) + 2
  ko(jl,2) = ko(jl,3) + 1
  ko(jl,4) = ko(jl,3) - 1
  !
  IF (.NOT.oglobs) THEN
    ko(jl,1:2) = min(ko(jl,1:2),iinla)
    ko(jl,3:4) = max(ko(jl,3:4),1)    
  ENDIF  
  !
  ! extra latitudes are computed symetrically compared to the poles
  !
  IF (oglobs) THEN
    IF (ko(jl,3)==2) THEN
      pla(jl,4) = 2. * zsouthpole - pla(jl,1)
      pla(jl,3) = 2. * zsouthpole - pla(jl,2)
    ELSEIF (ko(jl,3)==3) THEN
      pla(jl,4) = 2. * zsouthpole - pla(jl,3)
    END IF
  ENDIF
  IF (oglobn) THEN
    IF (ko(jl,3)==iinla-2) THEN
      pla(jl,1) = 2. * znorthpole - pla(jl,4)
      pla(jl,2) = 2. * znorthpole - pla(jl,3)
    ELSEIF (ko(jl,3)==iinla-3) THEN
     pla(jl,1) = 2. * znorthpole - pla(jl,2)
    END IF
  ENDIF
  !
  IF ((ko(jl,4)<1).OR.any((ko(jl,:)>iinla))) THEN  
    CALL abor1_sfx('HORIBLE_SURF_INIT: INPUT DOMAIN SMALLER THAN OUTPUT ONE - LATITUDE')
  END IF
  !
END DO
!
IF (lhook) CALL dr_hook('HORIBL_SURF_INIT',1,zhook_handle)
!
!------------------------------------------------------------------------------
!
!
END SUBROUTINE horibl_surf_init