get_mesh_index_conf_proj.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_conf_proj(KGRID_PAR,KSSO,PGRID_PAR,PLAT,PLON,KINDEX,KISSOX,KISSOY)
!     ###############################################################
!
!!**** *GET_MESH_INDEX_CONF_PROJ* get the grid mesh where point (lat,lon) is located
!!
!!    PURPOSE
!!    -------
!!
!!    AUTHOR
!!    ------
!!
!!    V. Masson         Meteo-France
!!
!!    MODIFICATION
!!    ------------
!!
!!    Original    12/09/95
!!    J.Escobar  22/10/2011 : reintroduce optimisation for JI/JJ number of lines computation
!!     J. Escobar  06/2013  : modif for REAL*4
!!
!----------------------------------------------------------------------------
!
!*    0.     DECLARATION
!            -----------
!
USE modd_get_mesh_index_conf_proj, ONLY : xlat0, xlon0, xrpk, xbeta,    &
                                            xlator, xlonor, nimax, njmax, &
                                            xxlim, xylim  
USE mode_gridtype_conf_proj
!
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
USE modd_csts ,ONLY : xsurf_epsilon
!
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, DIMENSION(SIZE(PLON))       :: zlon     ! longitude
!
REAL, DIMENSION(:), ALLOCATABLE   :: zx       ! X conformal coordinate
REAL, DIMENSION(:), ALLOCATABLE   :: zy       ! Y conformal coordinate
REAL, DIMENSION(:), ALLOCATABLE   :: zdx      ! X grid mesh size
REAL, DIMENSION(:), ALLOCATABLE   :: zdy      ! Y grid mesh size
REAL                              :: zdxlim   ! X grid mesh size
REAL                              :: zdylim   ! Y grid mesh size
!
INTEGER                           :: ji       ! loop counter in x
INTEGER                           :: jj       ! loop counter in y
INTEGER                           :: jl       ! loop counter on input points
REAL(KIND=JPRB) :: zhook_handle
!----------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('GET_MESH_INDEX_CONF_PROJ',0,zhook_handle)
IF (.NOT. ALLOCATED(xxlim)) THEN
!
!*    1.     Uncode parameters of the grid
!            -----------------------------
!
  CALL get_gridtype_conf_proj(pgrid_par,xlat0,xlon0,xrpk,xbeta,&
                                xlator,xlonor,nimax,njmax        )  
!
  ALLOCATE(zx(nimax*njmax))
  ALLOCATE(zy(nimax*njmax))
  ALLOCATE(zdx(nimax*njmax))
  ALLOCATE(zdy(nimax*njmax))
!
  CALL get_gridtype_conf_proj(pgrid_par,                       &
                                px=zx,py=zy,pdx=zdx,pdy=zdy      )  
!
!*    2.     Limits of grid meshes in x and y
!            --------------------------------
!
  ALLOCATE(xxlim(nimax+1))
  xxlim(1) = zx(1) - zdx(1)/2.
  DO ji=2,nimax
    jl = ji
    xxlim(ji) = zx(jl) - (zdx(jl-1)+zdx(jl))/4.
  END DO
  xxlim(nimax+1) = zx(nimax) + zdx(nimax)/2.

  ALLOCATE(xylim(njmax+1))
  xylim(1) = zy(1) - zdy(1)/2.
  DO jj=2,njmax
    jl = 1 + (jj-1) * nimax
    xylim(jj) = zy(jl) - (zdy(jl-nimax)+zdy(jl))/4.
  END DO
  xylim(njmax+1) = zy(1+(njmax-1)*nimax) + zdy(1+(njmax-1)*nimax)/2.
!
!
  DEALLOCATE(zx )
  DEALLOCATE(zy )
  DEALLOCATE(zdx)
  DEALLOCATE(zdy)
END IF
!
zdxlim = xxlim(2) - xxlim(1)
zdylim = xylim(2) - xylim(1)
!
!*    2.     Reshifts the longitudes with respect to projection reference point
!            ------------------------------------------------------------------
!
zlon(:) = plon(:)+nint((xlon0-plon(:)+180.0*xsurf_epsilon)/360.)*360.
!
!*    3.     Projection
!            ----------
!
ALLOCATE(zx(SIZE(plat)))
ALLOCATE(zy(SIZE(plat)))
!
 CALL xy_conf_proj(xlat0,xlon0,xrpk,xbeta,xlator,xlonor, &
                    zx,zy,plat,zlon                       )  
!
!
!
!*    5.     Localisation of the data points on (x,y) grid
!            ---------------------------------------------
!
kindex(:,:) = 0
!
DO jl=1,SIZE(plon)
  IF (     zx(jl)<xxlim(1) .OR. zx(jl)>=xxlim(nimax+1) &
        .OR. zy(jl)<xylim(1) .OR. zy(jl)>=xylim(njmax+1) ) THEN  
    kindex(1,jl) = 0
    IF (ksso/=0) THEN
      kissox(1,jl) = 0
      kissoy(1,jl) = 0
    END IF
    cycle
  END IF
  ji = min(int( (zx(jl) - xxlim(1))/zdxlim+1),nimax)
  jj = min(int( (zy(jl) - xylim(1))/zdylim+1),njmax)

  kindex(1,jl) = (jj-1) * nimax + 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) * (zx(jl)-xxlim(ji))/(xxlim(ji+1)-xxlim(ji)) )
    kissoy(1,jl) = 1 + int( float(ksso) * (zy(jl)-xylim(jj))/(xylim(jj+1)-xylim(jj)) )
  END IF
END DO
!
!-------------------------------------------------------------------------------
DEALLOCATE(zx )
DEALLOCATE(zy )
IF (lhook) CALL dr_hook('GET_MESH_INDEX_CONF_PROJ',1,zhook_handle)
!-------------------------------------------------------------------------------
!
END SUBROUTINE get_mesh_index_conf_proj