get_mesh_index_ign.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 get_mesh_index_ign(KGRID_PAR,KSSO,PGRID_PAR,PLAT,PLON,&
                                kindex,kissox,kissoy,pvalue,pnodata)
!     ###############################################################
!
!!**** *GET_MESH_INDEX_IGN* get the grid mesh where point (lat,lon) is located
!!
!!    PURPOSE
!!    -------
!!
!!    AUTHOR
!!    ------
!!
!!    E. Martin         Meteo-France
!!
!!    MODIFICATION
!!    ------------
!!
!!    Original    10/2007  
!!
!----------------------------------------------------------------------------
!
!*    0.     DECLARATION
!            -----------
!
USE modd_get_mesh_index_ign, ONLY : xxlim, xylim, xx_min, xx_max, xy_min, &
                                    xy_max, xdx, xdy, xxlims, nxids, xdx_max,&
                                    nfracd
USE modd_point_overlay
!
USE mode_gridtype_ign
!
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*    0.1    Declaration of arguments
!            ------------------------
!
INTEGER,                       INTENT(IN)   :: kgrid_par ! size of PGRID_PAR
INTEGER,                       INTENT(IN)   :: ksso      ! number of subgrid mesh in each direction
REAL,    DIMENSION(:),        INTENT(IN)    :: pgrid_par ! grid parameters
REAL,    DIMENSION(:),        INTENT(IN)    :: plat      ! latitude of the point
REAL,    DIMENSION(:),        INTENT(IN)    :: plon      ! longitude of the point
INTEGER, DIMENSION(:,:),      INTENT(OUT)   :: kindex    ! index of the grid mesh where the point is
INTEGER, DIMENSION(:,:),      INTENT(OUT)   :: kissox    ! X index of the subgrid mesh
INTEGER, DIMENSION(:,:),      INTENT(OUT)   :: kissoy    ! Y index of the subgrid mesh
!
REAL, DIMENSION(:), OPTIONAL, INTENT(IN)    :: pvalue  ! value of the point to add
REAL, OPTIONAL, INTENT(IN) :: pnodata
!
!*    0.2    Declaration of other local variables
!            ------------------------------------
!
INTEGER                           :: ilambert ! Lambert type
!
REAL, DIMENSION(:), ALLOCATABLE   :: zx       ! X Lambert   coordinate
REAL, DIMENSION(:), ALLOCATABLE   :: zy       ! Y Lambert   coordinate
!
REAL, DIMENSION(:), ALLOCATABLE   :: zxlim       ! X Lambert   coordinate
!
REAL, DIMENSION(SIZE(PLAT))       :: zvalue
!
LOGICAL, DIMENSION(SIZE(PLAT)) :: gmask
!
REAL  :: zvalx
!
REAL :: znodata
!
INTEGER :: isize, ifact
INTEGER                           :: il, icpt       ! Grid dimension
INTEGER                           :: jl       ! loop counter in lambert grid
INTEGER                           :: ji, jj       ! loop counter on input points
INTEGER, DIMENSION(SIZE(PLAT),2) :: ici
INTEGER, DIMENSION(1)             :: idx0
REAL(KIND=JPRB) :: zhook_handle
!----------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('GET_MESH_INDEX_IGN_1',0,zhook_handle)
!
IF (present(pvalue) .AND. present(pnodata)) THEN
  zvalue(:) = pvalue(:)
  znodata = pnodata
ELSE
  zvalue(:) = 1
  znodata = 0
ENDIF
!
IF (.NOT. ALLOCATED(xxlim)) THEN
!
!*    1.     Uncode parameters of the grid
!            -----------------------------
!
  CALL get_gridtype_ign(pgrid_par,klambert=ilambert,kl=il)
!
  ifact = floor(sqrt(float(il)))
  isize = floor(float(il) / ifact)
  ALLOCATE(nfracd(ifact+1))
  nfracd(1) = 1
  nfracd(ifact+1) = il
  DO jj=2,ifact
    nfracd(jj) = 1 + (jj-1)*isize
  ENDDO
!
  ALLOCATE(zx(il))
  ALLOCATE(zy(il))
  ALLOCATE(xdx(il))
  ALLOCATE(xdy(il))
  ALLOCATE(xxlim(il))
  ALLOCATE(xylim(il))

  ALLOCATE(xxlims(0:il))
  ALLOCATE(nxids(il))

  ALLOCATE(zxlim(il))  
!
  CALL get_gridtype_ign(pgrid_par,px=zx,py=zy,pdx=xdx,pdy=xdy)
!
!*    2.     Limits of grid meshes in x and y
!            --------------------------------
!
  xxlim(:)=zx(:)-xdx(:)/2.
  xylim(:)=zy(:)-xdy(:)/2.
!
  xx_min = minval(xxlim)
  xx_max = maxval(xxlim+xdx)
  xy_min = minval(xylim)
  xy_max = maxval(xylim+xdy)

  xdx_max = minval(xdx)

  zxlim(:) = xxlim(:)

  zvalx = maxval(zxlim) + 1.
  DO ji=1,il
    idx0 = minloc(zxlim) 
    xxlims(ji) = zxlim(idx0(1))
    nxids(ji) = idx0(1)
    zxlim(idx0(1)) = zvalx
  ENDDO
  xxlims(0) = xxlims(1) - xdx_max -1.

  DEALLOCATE(zxlim)
  DEALLOCATE(zx )
  DEALLOCATE(zy )
  !
ENDIF
IF (lhook) CALL dr_hook('GET_MESH_INDEX_IGN_1',1,zhook_handle)
IF (lhook) CALL dr_hook('GET_MESH_INDEX_IGN_2',0,zhook_handle)
!
!*    3.     Projection
!            ----------
!
ALLOCATE(zx(SIZE(plat)))
ALLOCATE(zy(SIZE(plat)))
!
  CALL get_gridtype_ign(pgrid_par,klambert=ilambert)
!
!       write (*,*),'APPEL  XY_IGN'
  CALL xy_ign(ilambert,zx,zy,plat,plon)
!   write(*,*) '    X Y SORTIE : ', ZX(1),ZY(1)
!
gmask(:) = .false.
DO jl=1,SIZE(plat)
  IF (     zx(jl)<xx_min .OR. zx(jl)>xx_max  &
      .OR. zy(jl)<xy_min .OR. zy(jl)>xy_max ) gmask(jl) = .true.
ENDDO
!
!*    5.     Localisation of the data points on (x,y) grid
!            ---------------------------------------------
!
IF (lhook) CALL dr_hook('GET_MESH_INDEX_IGN_2',1,zhook_handle)
IF (lhook) CALL dr_hook('GET_MESH_INDEX_IGN_3',0,zhook_handle)
!
ifact = SIZE(nfracd) - 1
!
kindex(:,:)=0
!
kissox(:,:) = 0
kissoy(:,:) = 0
!
ici(:,:) = 0
!$OMP PARALLEL DO PRIVATE(JL,JI,JJ)
DO jl=1,SIZE(plat)
  !
  IF (zvalue(jl)==znodata) cycle  
  IF (gmask(jl)) cycle
  !
  frac: &
  DO jj=ifact,1,-1
    !
    IF (zx(jl)>xxlims(nfracd(jj))) THEN
      !
      DO ji = nfracd(jj+1),nfracd(jj),-1
        IF (zx(jl)>xxlims(ji)) THEN
          ici(jl,2) = ji
          EXIT
        ENDIF
      ENDDO
      !
      DO ji = ici(jl,2),0,-1
        IF (zx(jl)>=xxlims(ji)+xdx_max) THEN
          ici(jl,1) = ji+1
          EXIT
        ENDIF
      ENDDO
      !
      EXIT frac
      !
    ENDIF 
    !
  ENDDO frac
  !
ENDDO
!$OMP END PARALLEL DO
!
IF (lhook) CALL dr_hook('GET_MESH_INDEX_IGN_3',1,zhook_handle)
IF (lhook) CALL dr_hook('GET_MESH_INDEX_IGN_4',0,zhook_handle)
!
DO jl=1,SIZE(plat)
  !
  IF (zvalue(jl)==znodata) cycle  
  IF (gmask(jl)) cycle
  !
  icpt = 0
  DO ji=ici(jl,1),ici(jl,2)
    !
    IF (zy(jl)>xylim(nxids(ji)) .AND. zy(jl)<xylim(nxids(ji))+xdy(nxids(ji)) &
    .AND. zx(jl)<xxlims(ji)+xdx(nxids(ji))) THEN
      !
      icpt = icpt + 1
      !
      kindex(icpt,jl) = nxids(ji)
      !
      IF (ksso/=0) THEN
        kissox(icpt,jl) = 1 + int( float(ksso) * (zx(jl)-xxlim(nxids(ji)))/xdx(nxids(ji)) )   
        kissoy(icpt,jl) = 1 + int( float(ksso) * (zy(jl)-xylim(nxids(ji)))/xdy(nxids(ji)) ) 
      ENDIF     
      !
      IF (icpt==novmx) EXIT
      !
    ENDIF 
    !
  ENDDO
  !
ENDDO
!
IF (lhook) CALL dr_hook('GET_MESH_INDEX_IGN_4',1,zhook_handle)
!-------------------------------------------------------------------------------
DEALLOCATE(zx )
DEALLOCATE(zy )
!-------------------------------------------------------------------------------
!
END SUBROUTINE get_mesh_index_ign