get_mesh_index_lonlatval.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_lonlatval(KGRID_PAR,KSSO,PGRID_PAR,PLAT,PLON,KINDEX,KISSOX,KISSOY)
!     ###############################################################
!
!!**** *GET_MESH_INDEX_LONLATVAL* 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_lonlatval, ONLY : xxlim, xylim, xx_min, xx_max, xy_min, &
                                      xy_max, xdx, xdy, xxlims, nxids, xdx_max,&
                                      nfracd  
USE mode_gridtype_lonlatval
!
USE modd_point_overlay
!
!
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
!
!*    0.2    Declaration of other local variables
!            ------------------------------------
!
REAL  :: xlon0
REAL  :: zvalx
!
REAL, DIMENSION(SIZE(PLON)) :: zlon
!
REAL, DIMENSION(:), ALLOCATABLE   :: zx       ! X Lambert   coordinate
REAL, DIMENSION(:), ALLOCATABLE   :: zy       ! Y Lambert   coordinate
REAL, DIMENSION(:), ALLOCATABLE   :: zxlim       ! X Lambert   coordinate
!
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
!
LOGICAL, DIMENSION(SIZE(PLAT)) :: gmask
!
REAL(KIND=JPRB) :: zhook_handle
!----------------------------------------------------------------------------
!
!
IF (lhook) CALL dr_hook('GET_MESH_INDEX_LONLATVAL',0,zhook_handle)
IF (.NOT. ALLOCATED(xxlim)) THEN
!
!*    1.     Uncode parameters of the grid
!            -----------------------------
!
  CALL get_gridtype_lonlatval(pgrid_par,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_lonlatval(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 )
  
END IF
!
xlon0 = 0.5*(xx_min+xx_max)
!
!*    3.     Projection
!            ----------
!
  CALL get_gridtype_lonlatval(pgrid_par)
!
  zlon(:) = plon(:)+nint((xlon0-plon(:))/360.)*360.
!
gmask(:) = .false.
DO jl=1,SIZE(plat)
  IF (     zlon(jl)<xx_min .OR. zlon(jl)>xx_max  &
      .OR. plat(jl)<xy_min .OR. plat(jl)>xy_max ) gmask(jl) = .true.
ENDDO

!*    5.     Localisation of the data points on (x,y) grid
!            ---------------------------------------------
!
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 (gmask(jl)) cycle
  !
  frac: &
  DO jj=ifact,1,-1
    !
    IF (zlon(jl)>xxlims(nfracd(jj))) THEN
      !
      DO ji = nfracd(jj+1),nfracd(jj),-1
        IF (zlon(jl)>xxlims(ji)) THEN
          ici(jl,2) = ji
          EXIT
        ENDIF
      ENDDO
      !
      DO ji = ici(jl,2),0,-1
        IF (zlon(jl)>=xxlims(ji)+xdx_max) THEN
          ici(jl,1) = ji+1
          EXIT
        ENDIF
      ENDDO
      !
      EXIT frac
      !
    ENDIF 
    !
  ENDDO frac
  !
ENDDO
!$OMP END PARALLEL DO
!
DO jl=1,SIZE(plat)
  !
  IF (gmask(jl)) cycle
  !
  icpt = 0
  DO ji=ici(jl,1),ici(jl,2)
    !
    IF (plat(jl)>xylim(nxids(ji)) .AND. plat(jl)<xylim(nxids(ji))+xdy(nxids(ji)) &
    .AND. zlon(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) * (zlon(jl)-xxlim(nxids(ji)))/xdx(nxids(ji)) )   
        kissoy(icpt,jl) = 1 + int( float(ksso) * (plat(jl)-xylim(nxids(ji)))/xdy(nxids(ji)) ) 
      ENDIF     
      !
      IF (icpt==novmx) EXIT
      !
    ENDIF 
    !
  ENDDO
  !
ENDDO
!
IF (lhook) CALL dr_hook('GET_MESH_INDEX_LONLATVAL',1,zhook_handle)
!-------------------------------------------------------------------------------
!
END SUBROUTINE get_mesh_index_lonlatval