horibl_surf_value.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_value(KILEN,KOLEN,PAROUT,OINTERP,PARIN,KLSMIN,&
                                 POLO,POLA,PLA,PLOP,&
                                 KMASKIN,KLSMOUT )  
!   ###########################################################################
!
!!****  *HORIBL_SURF_VALUE* - 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 PARIN 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)  :: KILEN
INTEGER,                   INTENT(IN)  :: KOLEN   ! size of output array
REAL,    DIMENSION(:), INTENT(OUT) :: PAROUT  ! output array
LOGICAL, DIMENSION(:), INTENT(IN)  :: OINTERP ! .true. where physical value is needed
REAL, DIMENSION(:,:), INTENT(IN) :: PARIN
INTEGER, DIMENSION(:,:), INTENT(IN) :: KLSMIN
REAL, DIMENSION(:), INTENT(IN) :: POLO, POLA
REAL, DIMENSION(:,:), INTENT(IN) :: PLA
REAL, DIMENSION(:,:), INTENT(IN) :: PLOP
INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL  :: KMASKIN  ! input land/sea mask
INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL  :: KLSMOUT ! output land/sea mask
!
!*      0.2. Declaration of local variables
!            
! Weights of the latitudes and of the points
REAL, DIMENSION(4) :: ZWV
REAL, DIMENSION(12) :: ZWP  
  ! Land/sea mask coefficient for each point : 0 -> point not taken in account,
  !                                            1 -> point taken in account
REAL, DIMENSION(12) :: ZLSMP
REAL, DIMENSION(4) :: ZLSMV
REAL :: ZLSMTOT  
 ! Variables implied in the extension procedure
LOGICAL :: LDLSM     ! Specify if land/sea mask is present or not
 ! Loop counters
INTEGER :: JL, JI    ! Dummy counter
!
!------------------------------------------------------------------------------
REAL, DIMENSION(3) :: ZP
REAL :: ZMAX, ZT     ! Max of 12 surrounding values
REAL :: ZMIN      ! Min of 12 surrounding values
INTEGER :: JL2    ! Dummy counter
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_VALUE_1',0,zhook_handle)
!
ldlsm = .false.
IF (PRESENT(kmaskin) .AND. PRESENT(klsmout)) ldlsm = .true.
!
!------------------------------------------------------------------------------
!
!*     3.   LOOP OVER ALL THE POINTS OF THE OUTPUT GRID
!           -------------------------------------------
!
parout(:) = xundef
!
IF (lhook) CALL dr_hook('HORIBL_SURF_VALUE_1',1,zhook_handle)
!$OMP PARALLEL PRIVATE(ZHOOK_HANDLE_OMP)
IF (lhook) CALL dr_hook('HORIBL_SURF_VALUE_2',0,zhook_handle_omp)
!$OMP DO SCHEDULE(DYNAMIC,1) PRIVATE(JL,ZLSMP,ZLSMV,ZLSMTOT,ZP,ZT,&
!$OMP   ZWV,ZWP,JI,ZMIN,ZMAX,JL2)
DO jl = 1, kolen
!
  IF (.NOT. ointerp(jl)) cycle
!
!*  3.2 Land / Sea mask
!
  zlsmp(:) = 1.
  zlsmv(:) = 1.
  !
  IF (ldlsm) THEN
    !
    DO ji = 1,12
      IF (klsmin(jl,ji).NE.klsmout(jl)) zlsmp(ji) = 0.
    ENDDO
    !
    zlsmv(1) = zlsmp(1) + zlsmp(2)
    zlsmv(2) = zlsmp(3) + zlsmp(4) + zlsmp(5) + zlsmp(6)
    zlsmv(3) = zlsmp(7) + zlsmp(8) + zlsmp(9) + zlsmp(10)
    zlsmv(4) = zlsmp(11) + zlsmp(12)
    !
    zlsmv(:) = min(zlsmv(:),1.)
    zlsmtot = min(sum(zlsmv),1.)
    !
    IF (zlsmv(1) < 1.e-3) zlsmp(1:2) = 1.
    IF (zlsmv(2) < 1.e-3) zlsmp(3:6) = 1.
    IF (zlsmv(3) < 1.e-3) zlsmp(7:10) = 1.
    IF (zlsmv(4) < 1.e-3) zlsmp(11:12) = 1.
    IF (zlsmtot < 1.e-3) zlsmv(:) = 1.
    !
  ENDIF
!
!*  3.3 Weight of points
!
  zp(1) = pla(jl,1) - pla(jl,2)
  zp(2) = pla(jl,1) - pla(jl,3)
  zp(3) = pla(jl,2) - pla(jl,3)
  !
  zt = pola(jl) - pla(jl,1)
  zwv(1) = zlsmv(1) * (1.+zlsmv(2)*zt/zp(1)) * (1.+zlsmv(3)*zt/zp(2)) * (1.+zlsmv(4)*zt/(pla(jl,1)-pla(jl,4)))  
  zt = pola(jl) - pla(jl,2)
  zwv(2) = zlsmv(2) * (1.-zlsmv(1)*zt/zp(1)) * (1.+zlsmv(3)*zt/zp(3)) * (1.+zlsmv(4)*zt/(pla(jl,2)-pla(jl,4))) 
  zt = pola(jl) - pla(jl,3)
  zwv(3) = zlsmv(3) * (1.-zlsmv(1)*zt/zp(2)) * (1.-zlsmv(2)*zt/zp(3)) * (1.+zlsmv(4)*zt/(pla(jl,3)-pla(jl,4))) 
  zwv(4) = 1. - zwv(1) - zwv(2) - zwv(3)
!
  ! 3.3.1 northern
  zwp(1)  = zlsmp(1) * (1.+zlsmp(2) *(polo(jl) -plop(jl,1))/(plop(jl,1) -plop(jl,2)))
  zwp(2)  = 1. - zwp(1)
  ! 3.3.4. southern
  zwp(11) = zlsmp(11)* (1.+zlsmp(12)*(polo(jl)-plop(jl,11))/(plop(jl,11)-plop(jl,12)))
  zwp(12) = 1. - zwp(11)

  ! 3.3.2. north
  zp(1) = plop(jl,3) - plop(jl,4)
  zp(2) = plop(jl,3) - plop(jl,5)
  zp(3) = plop(jl,4) - plop(jl,5) 
  !
  zt = polo(jl) - plop(jl,3)
  zwp(3) = zlsmp(3) * (1.+zlsmp(4)*zt/zp(1)) * (1.+zlsmp(5)*zt/zp(2)) * (1.+zlsmp(6)*zt/(plop(jl,3)-plop(jl,6)))  
  zt = polo(jl) - plop(jl,4)
  zwp(4) = zlsmp(4) * (1.-zlsmp(3)*zt/zp(1)) * (1.+zlsmp(5)*zt/zp(3)) * (1.+zlsmp(6)*zt/(plop(jl,4)-plop(jl,6)))
  zt = polo(jl) - plop(jl,5)
  zwp(5) = zlsmp(5) * (1.-zlsmp(3)*zt/zp(2)) * (1.-zlsmp(4)*zt/zp(3)) * (1.+zlsmp(6)*zt/(plop(jl,5)-plop(jl,6)))  
  zwp(6) = 1. - zwp(3) - zwp(4) - zwp(5)
!
  ! 3.3.3. south
  zp(1) = plop(jl,7) - plop(jl,8)
  zp(2) = plop(jl,7) - plop(jl,9)
  zp(3) = plop(jl,8) - plop(jl,9)
  !
  zt = polo(jl) - plop(jl,7)
  zwp(7)  = zlsmp(7) * (1.+zlsmp(8)*zt/zp(1)) * (1.+zlsmp(9)*zt/zp(2)) * (1.+zlsmp(10)*zt/(plop(jl,7)-plop(jl,10)))  
  zt = polo(jl) - plop(jl,8)
  zwp(8)  = zlsmp(8) * (1.-zlsmp(7)*zt/zp(1)) * (1.+zlsmp(9)*zt/zp(3)) * (1.+zlsmp(10)*zt/(plop(jl,8)-plop(jl,10)))  
  zt = polo(jl) - plop(jl,9)
  zwp(9)  = zlsmp(9) * (1.-zlsmp(7)*zt/zp(2)) * (1.-zlsmp(8)*zt/zp(3)) * (1.+zlsmp(10)*zt/(plop(jl,9)-plop(jl,10)))  
  zwp(10) = 1. - zwp(7) - zwp(8) - zwp(9)
!
! In order to exclude undef values from computation of PAROUT, 
! weights w2, w6, w10, w12 and wss which can be numerically very low
! because they are residual, are set to 0
!
  IF (abs(zwp(2)) < 1.e-10) zwp(2) =0.
  IF (abs(zwp(6)) < 1.e-10) zwp(6) =0.
  IF (abs(zwp(10))< 1.e-10) zwp(10)=0.
  IF (abs(zwp(12))< 1.e-10) zwp(12)=0.
  IF (abs(zwv(4)) < 1.e-10) zwv(4) =0.
!
      ! 3.3.5. longitude weight x latitude weight
  zwp(1:2) = zwp(1:2) * zwv(1)
  zwp(3:6) = zwp(3:6) * zwv(2)
  zwp(7:10) = zwp(7:10) * zwv(3)
  zwp(11:12) = zwp(11:12) * zwv(4)
  WHERE (zwp(:)<1.e-4) zwp(:) = 0.
  zwp(:) = zwp(:) / sum(zwp)
!
  parout(jl) = 0.
  DO ji = 1,12
    parout(jl) = parout(jl) + zwp(ji) * parin(jl,ji)
  ENDDO
!
! For surface fields, the interpoalted value is bounded 
! by the min max values of the initial field
  
  IF (PRESENT(kmaskin)) THEN

    zmin=xundef
    zmax=xundef

    DO jl2=1,12
      IF (parin(jl,jl2)==xundef) cycle

      IF ((zmax==xundef)) THEN
        zmax=parin(jl,jl2)
        zmin=parin(jl,jl2)
      ELSE
        zmax=max(zmax,parin(jl,jl2))
        zmin=min(zmin,parin(jl,jl2))
      ENDIF

    END DO
    
    parout(jl) = max(min(parout(jl),zmax),zmin)

  ENDIF

END DO
!$OMP END DO
IF (lhook) CALL dr_hook('HORIBL_SURF_VALUE_2',1,zhook_handle_omp)
!$OMP END PARALLEL
!
IF (lhook) CALL dr_hook('HORIBL_SURF_VALUE_3',0,zhook_handle)
WHERE(abs(parout-xundef)<1.e-6) parout=xundef
IF (lhook) CALL dr_hook('HORIBL_SURF_VALUE_3',1,zhook_handle)
!
!------------------------------------------------------------------------------
!
!
END SUBROUTINE horibl_surf_value