get_mesh_index_lonlat_reg.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_lonlat_reg(KGRID_PAR,KSSO,PGRID_PAR,PLAT,PLON,&
                        kindex,kissox,kissoy,pvalue,pnodata)
!     ###############################################################
!
!!**** *GET_MESH_INDEX_LONLAT_REG* get the grid mesh where point (lat,lon) is located
!!
!!    PURPOSE
!!    -------
!!
!!    AUTHOR
!!    ------
!!
!!    V. Masson         Meteo-France
!!
!!    MODIFICATION
!!    ------------
!!
!!    Original    12/09/95
!!
!----------------------------------------------------------------------------
!
!*    0.     DECLARATION
!            -----------
!
USE modd_get_mesh_index_lonlat_reg, ONLY : xlonlim, xlatlim, nlat, nlon, xlon0
USE mode_gridtype_lonlat_reg
!
!
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                           :: ji       ! loop counter in x
INTEGER                           :: jj       ! loop counter in y
INTEGER                           :: jl       ! loop counter on input points
!
REAL    :: zlonmin ! minimum longitude (degrees)
REAL    :: zlonmax ! maximum longitude (degrees)
REAL    :: zlatmin ! minimum latitude  (degrees)
REAL    :: zlatmax ! maximum latitude  (degrees)
REAL    :: zdlon   ! longitude grid size
REAL    :: zdlat   ! latitude  grid size
!
REAL :: znodata
!
REAL, DIMENSION(SIZE(PLAT))       :: zvalue
!
REAL, DIMENSION(SIZE(PLON)) :: zlon
!
INTEGER, DIMENSION(SIZE(PLAT))    :: ici, icj
!
REAL(KIND=JPRB) :: zhook_handle
!
!----------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('GET_MESH_INDEX_LONLAT_REG_1',0,zhook_handle)
!
IF (present(pvalue) .AND. present(pnodata)) THEN
  zvalue(:) = pvalue(:)
  znodata = pnodata
ELSE
  zvalue(:) = 1
  znodata = 0
ENDIF
!
IF (.NOT. ALLOCATED(xlatlim)) THEN
!
!*    1.     Uncode parameters of the grid
!            -----------------------------
!
  CALL get_gridtype_lonlat_reg(pgrid_par,zlonmin,zlonmax, &
                                 zlatmin,zlatmax,nlon,nlat  )  
!
!----------------------------------------------------------------------------
!
!*    2.     Limits of grid meshes
!            ---------------------
!
  zdlon = (zlonmax-zlonmin) / float(nlon)
  zdlat = (zlatmax-zlatmin) / float(nlat)
!
  ALLOCATE(xlonlim(nlon+1))
  DO ji=1,nlon+1
    xlonlim(ji) = zlonmin + float(ji-1)*zdlon
  END DO

  ALLOCATE(xlatlim(nlat+1))
  DO ji=1,nlat+1
    xlatlim(ji) = zlatmin + float(ji-1)*zdlat
  END DO
!
  xlon0 = 0.5*(zlonmin+zlonmax)
!
END IF
IF (lhook) CALL dr_hook('GET_MESH_INDEX_LONLAT_REG_1',1,zhook_handle)
IF (lhook) CALL dr_hook('GET_MESH_INDEX_LONLAT_REG_2',0,zhook_handle)
!----------------------------------------------------------------------------
!
!*    3.     Reshifts the longitudes with respect to projection reference point
!            ------------------------------------------------------------------
!
!
zlon(:) = plon(:)+nint((xlon0-plon(:))/360.)*360.
!
!----------------------------------------------------------------------------
!
!*    4.     Localisation of the data points on (x,y) grid
!            ---------------------------------------------
!
IF (SIZE(plat)/=nlon*nlat) THEN
  kindex = 0
  kissox = 0
  kissoy = 0
END IF
!
IF (lhook) CALL dr_hook('GET_MESH_INDEX_LONLAT_REG_2',1,zhook_handle)
IF (lhook) CALL dr_hook('GET_MESH_INDEX_LONLAT_REG_3',0,zhook_handle)
!
ici(:) = 0
icj(:) = 0
!$OMP PARALLEL DO PRIVATE(JL,JJ)
DO jl=1,SIZE(plat)
  !
  IF (zvalue(jl)==znodata) cycle
  ! 
  DO jj=SIZE(xlonlim),1,-1
    IF (xlonlim(jj)<=zlon(jl)) THEN
      ici(jl) = jj
      EXIT
    ENDIF
  ENDDO
  DO jj=SIZE(xlatlim),1,-1
    IF (xlatlim(jj)<=plat(jl)) THEN
      icj(jl) = jj
      EXIT
    ENDIF
  ENDDO
  !
ENDDO
!$OMP END PARALLEL DO
!
IF (lhook) CALL dr_hook('GET_MESH_INDEX_LONLAT_REG_3',1,zhook_handle)
IF (lhook) CALL dr_hook('GET_MESH_INDEX_LONLAT_REG_4',0,zhook_handle)
!
kindex(:,:) = 0

DO jl=1,SIZE(plat)
  !
  IF (zvalue(jl)==znodata) cycle
  ! 
  IF (     zlon(jl)<xlonlim(1) .OR. zlon(jl)>=xlonlim(nlon+1) &
        .OR. plat(jl)<xlatlim(1) .OR. plat(jl)>=xlatlim(nlat+1) ) THEN

    IF (ksso/=0) THEN
      kissox(1,jl) = 0
      kissoy(1,jl) = 0
    END IF
    cycle

  END IF

  ji = ici(jl)
  jj = icj(jl)
  kindex(1,jl) = (jj-1) * nlon + ji
!
!
!*    6.     Localisation of the data points on in the subgrid of this mesh
!            --------------------------------------------------------------
!
  IF (ksso/=0) THEN
    kissox(1,jl) = 1 + int( float(ksso) * (zlon(jl)-xlonlim(ji))/(xlonlim(ji+1)-xlonlim(ji)) )
    kissoy(1,jl) = 1 + int( float(ksso) * (plat(jl)-xlatlim(jj))/(xlatlim(jj+1)-xlatlim(jj)) )
  END IF
END DO
IF (lhook) CALL dr_hook('GET_MESH_INDEX_LONLAT_REG_4',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE get_mesh_index_lonlat_reg