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(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_surfex_omp, ONLY : nblocktot
USE modd_surfex_mpi, ONLY : nrank
USE modd_get_mesh_index_lonlat_reg, ONLY : xlonlim, xlatlim, nlat, nlon, xlon0,&
                                           nfracdlat, nfracdlon
USE mode_gridtype_lonlat_reg
!
USE modd_point_overlay, ONLY : novmx
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*    0.1    Declaration of arguments
!            ------------------------
!
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
!
INTEGER :: IFACTLON, ISIZELON, IFACTLAT, ISIZELAT, ISIZE_OMP
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE, ZHOOK_HANDLE_OMP
!
!----------------------------------------------------------------------------
!
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)
  !
  ifactlon = floor(sqrt(float(nlon+1))) + 1
  isizelon = floor(float(nlon+1) / ifactlon)
  ALLOCATE(nfracdlon(ifactlon+1))
  DO jj=1,ifactlon
    nfracdlon(jj) = 1 + (jj-1) * isizelon
  ENDDO
  nfracdlon(ifactlon+1) = nlon+1
  !
  ifactlat = floor(sqrt(float(nlat+1))) + 1
  isizelat = floor(float(nlat+1) / ifactlat)
  ALLOCATE(nfracdlat(ifactlat+1))
  DO jj=1,ifactlat
    nfracdlat(jj) = 1 + (jj-1) * isizelat
  ENDDO
  nfracdlat(ifactlat+1) = nlat+1
  !
END IF
IF (lhook) CALL dr_hook('GET_MESH_INDEX_LONLAT_REG_1',1,zhook_handle)
!----------------------------------------------------------------------------
!
!*    3.     Reshifts the longitudes with respect to projection reference point
!            ------------------------------------------------------------------
!
!
!$OMP PARALLEL PRIVATE(ZHOOK_HANDLE_OMP)
IF (lhook) CALL dr_hook('GET_MESH_INDEX_LONLAT_REG_2',0,zhook_handle_omp)
!$OMP DO PRIVATE(JL)
DO jl = 1,SIZE(plon)
  IF (plon(jl)>=xlon0+360.) THEN
    zlon(jl) = plon(jl) - 360.
  ELSEIF (plon(jl)<=xlon0-360.) THEN
    zlon(jl) = plon(jl) + 360.
  ELSE
    zlon(jl) = plon(jl)
  ENDIF
ENDDO
!$OMP END DO 
IF (lhook) CALL dr_hook('GET_MESH_INDEX_LONLAT_REG_2',1,zhook_handle_omp)
!$OMP END PARALLEL
!
!ZLON(JL) = PLON(JL)+NINT((XLON0-PLON(JL))/360.)*360.
!
!----------------------------------------------------------------------------
!
!*    4.     Localisation of the data points on (x,y) grid
!            ---------------------------------------------
!
IF (lhook) CALL dr_hook('GET_MESH_INDEX_LONLAT_REG_2b',0,zhook_handle)
!
!IF (SIZE(PLAT)/=NLON*NLAT) THEN
!  KINDEX = 0
!  KISSOX = 0
!  KISSOY = 0
!END IF
!
ifactlat = SIZE(nfracdlat)-1
ifactlon = SIZE(nfracdlon)-1
!
IF (lhook) CALL dr_hook('GET_MESH_INDEX_LONLAT_REG_2b',1,zhook_handle)
!
isize_omp = max(1,SIZE(plat)/nblocktot)
!
!$OMP PARALLEL PRIVATE(ZHOOK_HANDLE_OMP)
IF (lhook) CALL dr_hook('GET_MESH_INDEX_LONLAT_REG_3',0,zhook_handle_omp)
!$OMP DO SCHEDULE(STATIC,ISIZE_OMP) PRIVATE(JL,JJ,JI)
DO jl=1,SIZE(plat)
  !
  ici(jl) = 0
  icj(jl) = 0
  !
  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) ) cycle
  !
  fraclat: &
  DO jj = ifactlat,1,-1
    IF (plat(jl)>xlatlim(nfracdlat(jj))) THEN
      DO ji = nfracdlat(jj+1),nfracdlat(jj)+1,-1
        IF (plat(jl)>=xlatlim(ji)) THEN
          icj(jl) = ji
          EXIT fraclat
        ENDIF
      ENDDO
      icj(jl) = nfracdlat(jj)
      EXIT fraclat
    ENDIF
  ENDDO fraclat
  !
  fraclon: &
  DO jj = ifactlon,1,-1
    IF (zlon(jl)>xlonlim(nfracdlon(jj))) THEN
      DO ji = nfracdlon(jj+1),nfracdlon(jj)+1,-1
        IF (zlon(jl)>=xlonlim(ji)) THEN
          ici(jl) = ji
          EXIT fraclon
        ENDIF
      ENDDO
      ici(jl) = nfracdlon(jj)
      EXIT fraclon
    ENDIF
  ENDDO fraclon
  !  
ENDDO
!$OMP END DO
IF (lhook) CALL dr_hook('GET_MESH_INDEX_LONLAT_REG_3',1,zhook_handle_omp)
!$OMP END PARALLEL
!
!$OMP PARALLEL PRIVATE(ZHOOK_HANDLE_OMP)
IF (lhook) CALL dr_hook('GET_MESH_INDEX_LONLAT_REG_4',0,zhook_handle_omp)
!$OMP DO PRIVATE(JL,JI,JJ)
DO jl=1,SIZE(plat)
  !
  IF (novmx>1) kindex(2:novmx,jl) = 0
  !
  IF (zvalue(jl)==znodata) THEN
    !
    kindex(1,jl) = 0
    ! 
  ELSEIF ( ici(jl)==0 .OR. icj(jl)==0 ) THEN    
    !
    kindex(1,jl) = 0
    !
    IF (ksso/=0) THEN
      kissox(1,jl) = 0
      kissoy(1,jl) = 0
    ENDIF
    !
  ELSE

    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

  ENDIF

END DO
!$OMP END DO
IF (lhook) CALL dr_hook('GET_MESH_INDEX_LONLAT_REG_4',1,zhook_handle_omp)
!$OMP END PARALLEL
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE get_mesh_index_lonlat_reg