bilin.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 (KLUOUT,PX1,PY1,PFIELD1,PX2,PY2,PFIELD2,OINTERP)
!     #########################################################################
!
!!****  *BILINEAR * - 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 modi_hor_extrapol_surf_cheap
!
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)  :: pfield1 ! FIELD on 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) :: pfield2 ! FIELD on model 2
LOGICAL, DIMENSION(:),INTENT(IN)  :: ointerp ! .true. where physical value is needed
!
!
!*       0.2    Declarations of local variables for print on FM file
!
! 
REAL, DIMENSION (:,:), ALLOCATABLE   :: zw       ! weight. 1 if defined, 
                                                 !         0 if FIELD=XUNDEF
REAL, DIMENSION (:,:), ALLOCATABLE   :: zw_field ! weight. * FIELD
REAL, DIMENSION (:,:), ALLOCATABLE   :: zfield_x ! FIELD at mesh interface
REAL, DIMENSION (:,:), ALLOCATABLE   :: zfield_y ! FIELD at mesh interface
REAL, DIMENSION (:,:), ALLOCATABLE   :: zfield_xy! FIELD at mesh corner
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, DIMENSION(SIZE(PFIELD2))    :: ici, icj
REAL                                 :: zc1_x    ! coefficient for left   points
REAL                                 :: zc2_x    ! coefficient for middle points
REAL                                 :: zc3_x    ! coefficient for right  points
REAL                                 :: zc1_y    ! coefficient for bottom points
REAL                                 :: zc2_y    ! coefficient for middle points
REAL                                 :: zc3_y    ! coefficient for top    points
!
INTEGER                              :: iiu       ! model 1 X size
INTEGER                              :: iju       ! model 1 Y size
!
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
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('BILIN_1',0,zhook_handle)

iiu=SIZE(pfield1,1)
iju=SIZE(pfield1,2)
!
!* weighting factor
!
ALLOCATE(zw(iiu,iju))
WHERE (pfield1/=xundef)
  zw=1.
ELSEWHERE
  zw=0.
END WHERE
!
!* weighted FIELD
!
ALLOCATE(zw_field(iiu,iju))
zw_field=zw*pfield1
!
!-------------------------------------------------------------------------------
!
!*       1.     FIELD type at grid mesh interfaces (in X directions)
!               ----------------------------------
!
ALLOCATE(zfield_x(iiu+1,iju))
!
!*       1.1    Standard case
!               -------------
!
IF (iiu>1)                                                      &
zfield_x(2:iiu  ,:) =  (zw_field(1:iiu-1,:)+zw_field(2:iiu,:))  &
                    / max(zw(1:iiu-1,:)+zw(2:iiu,:),zeps)  
zfield_x(1      ,:) =         zw_field(1  ,:)
zfield_x(  iiu+1,:) =         zw_field(iiu,:)
!
!
!*       1.2    FIELD type value is 0 in the grid mesh
!               --------------------------------------
!
WHERE ( pfield1(1:iiu,:)==0.)
  zfield_x(1:iiu,:)  = 0.
END WHERE
!
WHERE ( pfield1(1:iiu,:)==0.)
  zfield_x(2:iiu+1,:)  = 0.
END WHERE
!
!-------------------------------------------------------------------------------
!
!*       2.     FIELD type at grid mesh interfaces (in X directions)
!               ----------------------------------
!
ALLOCATE(zfield_y(iiu,iju+1))
!
!*       2.1    Standard case
!               -------------
!
IF (iju>1)                                                     &
zfield_y(:,2:iju  ) =  (zw_field(:,1:iju-1)+zw_field(:,2:iju)) &
                    / max(zw(:,1:iju-1)+zw(:,2:iju),zeps)  
zfield_y(:,1      ) =       zw_field(:,1      )
zfield_y(:,  iju+1) =       zw_field(:,  iju  )
!
!
!*       2.3    FIELD type value is 0 in the grid mesh
!               --------------------------------------
!
WHERE ( pfield1(:,1:iju)==0.)
  zfield_y(:,1:iju)  = 0.
END WHERE
!
WHERE ( pfield1(:,1:iju)==0.)
  zfield_y(:,2:iju+1)  = 0.
END WHERE
!
!-------------------------------------------------------------------------------
!
!*       3.     FIELD type at grid mesh corners
!               -------------------------------
!
ALLOCATE(zfield_xy(iiu+1,iju+1))
!
!*       3.1    Standard case
!               -------------
!

IF (iiu>1 .AND. iju>1)                                       &
zfield_xy(2:iiu  ,2:iju  ) =  (  zw_field(1:iiu-1,1:iju-1)   &
                                 + zw_field(1:iiu-1,2:iju  )   &
                                 + zw_field(2:iiu  ,1:iju-1)   &
                                 + zw_field(2:iiu  ,2:iju  ) ) &
                           / max(  zw(1:iiu-1,1:iju-1)   &
                                 + zw(1:iiu-1,2:iju  )   &
                                 + zw(2:iiu  ,1:iju-1)   &
                                 + zw(2:iiu  ,2:iju  ) , zeps)  
!
IF (iju>1)                                                   &
zfield_xy(1      ,2:iju  ) = (  zw_field(1       ,1:iju-1)   &
                                + zw_field(1       ,2:iju  ) ) &
                          / max(  zw(1       ,1:iju-1)   &
                                + zw(1       ,2:iju  ) , zeps)  
!
IF (iju>1)                                                   &
zfield_xy(  iiu+1,2:iju  ) = (  zw_field(iiu     ,1:iju-1)   &
                                + zw_field(iiu     ,2:iju  ) ) &
                          / max(  zw(iiu     ,1:iju-1)   &
                                + zw(iiu     ,2:iju  ) , zeps)  
!
IF (iiu>1)                                                   &
zfield_xy(2:iiu  ,1      ) = (   zw_field(1:iiu-1,1      )   &
                                 + zw_field(2:iiu  ,1      ) ) &
                          / max(   zw(1:iiu-1,1      )   &
                                 + zw(2:iiu  ,1      ) , zeps)  
!
IF (iiu>1)                                                   &
zfield_xy(2:iiu  ,iju+1  ) = (   zw_field(1:iiu-1,iju    )   &
                                 + zw_field(2:iiu  ,iju    ) ) &
                          / max(   zw(1:iiu-1,iju    )   &
                                 + zw(2:iiu  ,iju    ) , zeps)  
!
zfield_xy(1      ,1      ) =       (   pfield1(1      ,1      ) )
zfield_xy(iiu+1  ,1      ) =       (   pfield1(iiu    ,1      ) )
zfield_xy(1      ,iju+1  ) =       (   pfield1(1      ,iju    ) )
zfield_xy(iiu+1  ,iju+1  ) =       (   pfield1(iiu    ,iju    ) )
!
!*       3.2    FIELD type value is 0 in one grid mesh
!               --------------------------------------
!
WHERE ( pfield1(1:iiu,1:iju)==0.)
  zfield_xy(1:iiu,1:iju)  = 0.
END WHERE
!
WHERE ( pfield1(1:iiu,1:iju)==0.)
  zfield_xy(1:iiu,2:iju+1)  = 0.
END WHERE
!
WHERE ( pfield1(1:iiu,1:iju)==0.)
  zfield_xy(2:iiu+1,1:iju)  = 0.
END WHERE
!
WHERE ( pfield1(1:iiu,1:iju)==0.)
  zfield_xy(2:iiu+1,2:iju+1)  = 0.
END WHERE
!
!-------------------------------------------------------------------------------
!
!*       4.     Coordinates of points between input grid points
!               -----------------------------------------------
!
IF (iiu>1) THEN
  zx(iiu+1) = 1.5*px1(iiu)    -0.5*px1(iiu-1)
  zx(2:iiu) = 0.5*px1(1:iiu-1)+0.5*px1(2:iiu)
  zx(1)     = 1.5*px1(1)      -0.5*px1(2)
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(iju+1) = 1.5*py1(iju)    -0.5*py1(iju-1)
  zy(2:iju) = 0.5*py1(1:iju-1)+0.5*py1(2:iju)
  zy(1)     = 1.5*py1(1)      -0.5*py1(2)
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
!               -------------
!
pfield2(:) = xundef
!
!* loop on all output grid points
!
IF (lhook) CALL dr_hook('BILIN_1',1,zhook_handle)
IF (lhook) CALL dr_hook('BILIN_2',0,zhook_handle)
!
ici(:) = 1
icj(:) = 1
!$OMP PARALLEL DO PRIVATE(JL,JJ)
DO jl=1,SIZE(pfield2)
  DO jj=SIZE(zx),1,-1
    IF (zx(jj)<=px2(jl)) THEN
      ici(jl) = jj
      EXIT
    ENDIF
  ENDDO
  DO jj=SIZE(zy),1,-1
    IF (zy(jj)<=py2(jl)) THEN
      icj(jl) = jj
      EXIT
    ENDIF
  ENDDO
ENDDO
!$OMP END PARALLEL DO

IF (lhook) CALL dr_hook('BILIN_2',1,zhook_handle)
IF (lhook) CALL dr_hook('BILIN_3',0,zhook_handle)

DO jl=1,SIZE(pfield2,1)
!
!* interpolation weights in X direction
!
  ji=ici(jl)
  ji=max(min(ji,iiu),0)
  IF (px1(ji)<=px2(jl)) THEN
    zc2_x = (px2(jl)-zx(ji+1))/(px1(ji)-zx(ji+1))
    zc2_x = max(min(zc2_x,1.),0.)
    zc3_x = 1. - zc2_x
    zc1_x = 0.
  ELSE
    zc2_x = (px2(jl)-zx(ji))/(px1(ji)-zx(ji))
    zc2_x = max(min(zc2_x,1.),0.)
    zc1_x = 1. - zc2_x
    zc3_x = 0.
  END IF
!
!  interpolation weights in Y direction
!
  jj=icj(jl)
  jj=max(min(jj,iju),0)
  IF (py1(jj)<=py2(jl)) THEN
    zc2_y = (py2(jl)-zy(jj+1))/(py1(jj)-zy(jj+1))
    zc2_y = max(min(zc2_y,1.),0.)
    zc3_y = 1. - zc2_y
    zc1_y = 0.
  ELSE
    zc2_y = (py2(jl)-zy(jj))/(py1(jj)-zy(jj))
    zc2_y = max(min(zc2_y,1.),0.)
    zc1_y = 1. - zc2_y
    zc3_y = 0.
  END IF
!
!* interpolation
!
  IF(pfield1(ji,jj) /= xundef)                       &
       pfield2(jl) =                                   &
            zc1_y * (   zc1_x * zfield_xy(ji  ,jj  )   &
                      + zc2_x * zfield_y(ji  ,jj  )   &
                      + zc3_x * zfield_xy(ji+1,jj  ) ) &
          + zc2_y * (   zc1_x * zfield_x(ji  ,jj  )   &
                      + zc2_x * pfield1(ji  ,jj  )   &
                      + zc3_x * zfield_x(ji+1,jj  ) ) &
          + zc3_y * (   zc1_x * zfield_xy(ji  ,jj+1)   &
                      + zc2_x * zfield_y(ji  ,jj+1)   &
                      + zc3_x * zfield_xy(ji+1,jj+1) )  

END DO
!
IF (lhook) CALL dr_hook('BILIN_3',1,zhook_handle)
IF (lhook) CALL dr_hook('BILIN_4',0,zhook_handle)
!
!-------------------------------------------------------------------------------
!
DEALLOCATE(zw)
DEALLOCATE(zw_field)
DEALLOCATE(zfield_x )
DEALLOCATE(zfield_y )
DEALLOCATE(zfield_xy)
!
!-------------------------------------------------------------------------------
!
WHERE(abs(pfield2-xundef)<1.e-6) pfield2=xundef
!
!------------------------------------------------------------------------------
!
!*     6.   EXTRAPOLATIONS IF SOME POINTS WERE NOT INTERPOLATED
!           ---------------------------------------------------
!
!* no missing point
IF (count(pfield2(:)==xundef .AND. ointerp(:))==0 .AND. lhook) CALL dr_hook('BILIN',1,zhook_handle)
IF (count(pfield2(:)==xundef .AND. ointerp(:))==0) RETURN
!
!* no data point
IF (count(pfield1(:,:)/=xundef)==0 .AND. lhook) CALL dr_hook('BILIN',1,zhook_handle)
IF (count(pfield1(:,:)/=xundef)==0) RETURN
!
WRITE(kluout,*) ' Remaining horizontal extrapolations'
WRITE(kluout,*) ' Total number of input data     : ',count(pfield1(:,:)/=xundef)
WRITE(kluout,*) ' Number of points to interpolate: ',count(pfield2(:)==xundef .AND. ointerp(:))
!
!* input grid coordinates
!
 CALL hor_extrapol_surf_cheap(kluout,'XY  ',py1,px1,pfield1,py2,px2,pfield2,ointerp)
!
IF (lhook) CALL dr_hook('BILIN_4',1,zhook_handle)
!-------------------------------------------------------------------------------
!
END SUBROUTINE bilin