bilin_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 bilin_coef (KLUOUT,PX1,PY1,PX2,PY2,PCX,PCY,KCI,KCJ)
!     #########################################################################
!
!!****  *BILIN_COEFEAR * - subroutine to interpolate surface FIELD
!!
!!    PURPOSE
!!    -------
!!
!!
!!**  METHOD
!!    ------
!!
!!     Interpolation is bilinear, and uses 9 grid points, located in the
!!   center of model 1 grid mesh, and at the boundaries of this grid
!!   mesh (2 X limits, 2 Y limits and 4 corners).
!!     This implies that the grid mesh values located around the model 1
!!   grid mesh are not used directly. The values at the boundaries of the
!!   grid mesh are defined by the average between the middle point
!!   (this grid mesh value), and the one in the considered direction.
!!   So the eight grid meshes around the considered grid mesh are used
!!   equally.
!!     This is important to note that these average values are erased
!!   and replaced by zero if they are at the limit of any grid
!!   mesh with the zero value. This allows to insure zero value in model 2
!!   grid meshes where there was not the considered class in corresponding
!!   model 1 grid mesh, and to insure continuity of the FIELD type
!!   at such boundaries.
!!
!!
!!    The arrays and array index are defined on the following (model1) grid:
!!
!!
!!        XFIELD                    XFIELD                    XFIELD
!!          *                         *                         *
!!       i-1,j+1                    i,j+1                    i+1,j+1
!!
!!
!!
!!                   ZFIELD_XY     ZFIELD_Y    ZFIELD_XY
!!                       *            *            *
!!                     i,j+1        i,j+1       i+1,j+1
!!
!!
!!
!!        XFIELD      ZFIELD_X      XFIELD      ZFIELD_X      XFIELD
!!          *            *            *            *            *
!!        i-1,j         i,j          i,j         i+1,j        i+1,j
!!
!!
!!
!!                    ZFIELD_XY     ZFIELD_Y    ZFIELD_XY
!!                       *            *            *
!!                      i,j          i,j         i+1,j
!!
!!
!!
!!        XFIELD                    XFIELD                    XFIELD
!!          *                         *                         *
!!       i-1,j-1                    i,j-1                    i+1,j-1
!!
!!
!!
!!
!!
!!    AUTHOR
!!    ------
!!
!!       V. Masson     * METEO-FRANCE *
!!
!!    MODIFICATIONS
!!    -------------
!!
!!      Original     01/2004
! TD&DD: added OpenMP directives

!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE modd_surfex_mpi, ONLY : nrank
USE modd_surf_par,   ONLY : xundef
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*       0.1   Declarations of dummy arguments :
!
INTEGER,              INTENT(IN)  :: KLUOUT  ! output listing logical unit
REAL, DIMENSION(:),   INTENT(IN)  :: PX1     ! X coordinate of the regular input grid
REAL, DIMENSION(:),   INTENT(IN)  :: PY1     ! Y coordinate of the regular input grid
REAL, DIMENSION(:),   INTENT(IN)  :: PX2     ! X coordinate of all points of output grid
REAL, DIMENSION(:),   INTENT(IN)  :: PY2     ! Y coordinate of all points of output grid
REAL, DIMENSION(:,:), INTENT(OUT) :: PCX, PCY
INTEGER, DIMENSION(:), INTENT(OUT):: KCI, KCJ

!
!*       0.2    Declarations of local variables for print on FM file
!
! 
REAL, DIMENSION (SIZE(PX1)+1)        :: ZX       ! X coordinate of left   limit of input meshes
REAL, DIMENSION (SIZE(PY1)+1)        :: ZY       ! Y coordinate of bottom limit of input meshes
!
!
INTEGER                              :: IIU       ! model 1 X size
INTEGER                              :: IJU       ! model 1 Y size
!
INTEGER                              :: IDMIN, IDMAX
INTEGER                              :: JI        ! grid 1 x index
INTEGER                              :: JJ        ! grid 1 y index
INTEGER                              :: IIN       ! loop counter on all input points
!
INTEGER                              :: JL        ! grid 2 index
!
REAL                                 :: ZEPS=1.e-3
REAL(KIND=JPRB) :: ZHOOK_HANDLE, ZHOOK_HANDLE_OMP
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('BILIN_COEF_1',0,zhook_handle)

iiu=SIZE(px1)
iju=SIZE(py1)
!
!*       4.     Coordinates of points between input grid points
!               -----------------------------------------------
!
IF (iiu>1) THEN
  zx(1)     = 1.5*px1(1)  -0.5*px1(2)
  zx(iiu+1) = 1.5*px1(iiu)-0.5*px1(iiu-1)
  DO jj = 2,iiu
    zx(jj) = 0.5*(px1(jj-1)+px1(jj))
  ENDDO  
ELSE
  zx(1)     = px1(1) - 1.e6 ! uniform field in X direction if only 1 point is
  zx(2)     = px1(1) + 1.e6 ! available. Arbitrary mesh length of 2000km assumed
END IF
!
IF (iju>1) THEN
  zy(1)     = 1.5*py1(1)  -0.5*py1(2)
  zy(iju+1) = 1.5*py1(iju)-0.5*py1(iju-1)
  DO jj = 2,iju
    zy(jj) = 0.5*(py1(jj-1)+py1(jj))
  ENDDO
ELSE
  zy(1)     = py1(1) - 1.e6 ! uniform field in Y direction if only 1 point is
  zy(2)     = py1(1) + 1.e6 ! available. Arbitrary mesh length of 2000km assumed
END IF
!
!-------------------------------------------------------------------------------
!
!*       5.     Interpolation
!               -------------
!
!* loop on all output grid points
!
kci(:) = 1
kcj(:) = 1
!
IF (lhook) CALL dr_hook('BILIN_COEF_1',1,zhook_handle)
!
!$OMP PARALLEL PRIVATE(ZHOOK_HANDLE_OMP)
IF (lhook) CALL dr_hook('BILIN_COEF_2',0,zhook_handle_omp)
!$OMP DO PRIVATE(JL,JJ)
DO jl=1,SIZE(px2)
  DO jj=SIZE(zx),1,-1
    IF (zx(jj)<=px2(jl)) THEN
      kci(jl) = jj
      EXIT
    ENDIF
  ENDDO
  DO jj=SIZE(zy),1,-1
    IF (zy(jj)<=py2(jl)) THEN
      kcj(jl) = jj
      EXIT
    ENDIF
  ENDDO
ENDDO
!$OMP END DO
IF (lhook) CALL dr_hook('BILIN_COEF_2',1,zhook_handle_omp)
!$OMP END PARALLEL
!
IF (lhook) CALL dr_hook('BILIN_COEF_3',0,zhook_handle)

DO jl=1,SIZE(px2)
!
!* interpolation weights in X direction
!
  ji=kci(jl)
  ji=max(min(ji,iiu),0)
  IF (px1(ji)<=px2(jl)) THEN
    pcx(jl,2) = (px2(jl)-zx(ji+1))/(px1(ji)-zx(ji+1))
    pcx(jl,2) = max(min(pcx(jl,2),1.),0.)
    pcx(jl,3) = 1. - pcx(jl,2)
    pcx(jl,1) = 0.
  ELSE
    pcx(jl,2) = (px2(jl)-zx(ji))/(px1(ji)-zx(ji))
    pcx(jl,2) = max(min(pcx(jl,2),1.),0.)
    pcx(jl,1) = 1. - pcx(jl,2)
    pcx(jl,3) = 0.
  END IF
!
!  interpolation weights in Y direction
!
  jj=kcj(jl)
  jj=max(min(jj,iju),0)
  IF (py1(jj)<=py2(jl)) THEN
    pcy(jl,2) = (py2(jl)-zy(jj+1))/(py1(jj)-zy(jj+1))
    pcy(jl,2) = max(min(pcy(jl,2),1.),0.)
    pcy(jl,3) = 1. - pcy(jl,2)
    pcy(jl,1) = 0.
  ELSE
    pcy(jl,2) = (py2(jl)-zy(jj))/(py1(jj)-zy(jj))
    pcy(jl,2) = max(min(pcy(jl,2),1.),0.)
    pcy(jl,1) = 1. - pcy(jl,2)
    pcy(jl,3) = 0.
  END IF
!
END DO
!
IF (lhook) CALL dr_hook('BILIN_COEF_3',1,zhook_handle)
!-------------------------------------------------------------------------------
!
END SUBROUTINE bilin_coef