get_mesh_index_gauss.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_gauss(KNBLINES,KGRID_PAR,KSSO,PGRID_PAR,PLAT,PLON,&
                                      kindex,kissox,kissoy,pvalue,pnodata)
!     ###############################################################
!
!!**** *GET_MESH_INDEX_GAUSS* get the grid mesh where point (lat,lon) is located
!!
!!    PURPOSE
!!    -------
!!
!!    AUTHOR
!!    ------
!!
!!    V. Masson         Meteo-France
!!
!!    MODIFICATION
!!    ------------
!!
!!    Original    12/09/95
!!
!----------------------------------------------------------------------------
!
!*    0.     DECLARATION
!            -----------
!
USE eggangles,ONLY : lola, angle_domain
!
USE modd_get_mesh_index_gauss, ONLY : xxcen, xycen, xyinf, xysup, xxinf, xxsup,       &
                                      nnlati, nnlopa, xlapo, xlopo, xcodil, xsints,   &
                                      xlon, xlat, xcost, xsintc, xcosn, xsinn, xlonp, &
                                      xlatp, xcosp, xsinp, xpi, x1, x2, xdr,          &
                                      nfracdx, nfracgx, nfracdy, nfacty, xxdif, xydif,&
                                      lrotstretch
!
USE mode_gridtype_gauss
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*    0.1    Declaration of arguments
!            ------------------------
!
INTEGER,                       INTENT(IN)   :: knblines
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  (degrees)
REAL,    DIMENSION(:),         INTENT(IN)   :: plon      ! longitude of the point (degrees)
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
!            ------------------------------------
!
REAL, DIMENSION(SIZE(PLAT))       :: zx, zy       ! pseudo longitude of input point
REAL, DIMENSION(SIZE(PLAT))       :: zvalue
REAL :: zpc2
REAL :: znodata
!
INTEGER, DIMENSION(SIZE(PLAT))    :: icj
INTEGER :: ilgrid, isize  ! number of grid points
INTEGER :: ifactx, isizex, isizey
INTEGER :: ji, jj, jl       ! loop counter in x
INTEGER :: jgrid, igrid0    ! loop counter on grid  points
!
REAL(KIND=JPRB) :: zhook_handle
!----------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('GET_MESH_INDEX_GAUSS_1',0,zhook_handle)
!
IF (present(pvalue) .AND. present(pnodata)) THEN
  zvalue(:) = pvalue(:)
  znodata = pnodata
ELSE
  zvalue(:) = 1
  znodata = 0
ENDIF
!
kindex(:,:) = 0
kissox(:,:) = 0
kissoy(:,:) = 0
!
IF (.NOT. ALLOCATED(nnlopa)) THEN
  !
  !*    1.     Gets parameters of the projection
  !            ---------------------------------
  !
  xpi = 4.*atan(1.)
  xdr = xpi / 180.
  !
  CALL get_gridtype_gauss(pgrid_par,nnlati,kl=ilgrid)
  !
  ALLOCATE(nnlopa(0:nnlati))
  ALLOCATE(xxcen(ilgrid))
  ALLOCATE(xycen(ilgrid))
  CALL get_gridtype_gauss(pgrid_par,nnlati,xlapo,xlopo,xcodil,nnlopa(1:nnlati), &
                          ilgrid,plat_xy=xycen,plon_xy=xxcen             )
  nnlopa(0) = 0
  !
  zpc2  = xcodil*xcodil
  x1 = 1.0-zpc2
  x2 = 1.0+zpc2
  !
  xlonp = angle_domain(xlopo,dom='0+',unit='D') * xdr
  xlatp = xlapo * xdr
  !
  xcosp = cos(xlatp)
  xsinp = sin(xlatp)
  !
  !*    2.     Limits of grid meshes in x and y
  !            --------------------------------
  !
  ALLOCATE(xxinf(ilgrid))
  ALLOCATE(xyinf(ilgrid))
  ALLOCATE(xxsup(ilgrid))
  ALLOCATE(xysup(ilgrid))
  ALLOCATE(xxdif(ilgrid))
  ALLOCATE(xydif(ilgrid))  
  !
  CALL gauss_grid_limits(nnlati,nnlopa(1:nnlati),xxinf,xxsup,xyinf,xysup)
  DO jj=1,ilgrid
    xxdif(jj) = 1. / (xxsup(jj) - xxinf(jj))
    xydif(jj) = 1. / (xysup(jj) - xyinf(jj))
  ENDDO
  !
  ifactx = floor(sqrt(float(nnlati))) + 1
  isizex = floor(float(nnlati) / ifactx)
  !
  ALLOCATE(nfracdx(0:ifactx))
  ALLOCATE(nfracgx(0:ifactx))
  !
  nfracdx(0) = 0
  nfracgx(0) = 1
  nfracdx(1) = 1
  nfracgx(1) = 1
  nfracdx(ifactx) = nnlati
  nfracgx(ifactx) = sum(nnlopa(:))
  DO jj=2,ifactx-1
    nfracdx(jj) = 1 + (jj-1) * isizex
    nfracgx(jj) = nfracgx(jj-1) + sum(nnlopa(nfracdx(jj-1):(jj-1)*isizex))
  ENDDO
  !
  !
  ALLOCATE(nfacty(nnlati))
  nfacty(:) = floor(sqrt(float(nnlopa(1:nnlati))))+1
  !
  ALLOCATE(nfracdy(nnlati,0:maxval(nfacty)))
  !
  DO jj=1,nnlati
    isizey = floor(float(nnlopa(jj)) / nfacty(jj))
    nfracdy(jj,0) = 0
    nfracdy(jj,1) = 1
    nfracdy(jj,nfacty(jj)) = nnlopa(jj)
    DO ji=2,nfacty(jj)-1
      nfracdy(jj,ji) = 1 + (ji-1) * isizey
    ENDDO
  ENDDO
  !
  !*    3.     Find if rotated pole and/or stretching to improve CPU time
  !            ----------------------------------------------------------
  !
  lrotstretch = .true.
  IF (xcodil==1.0.AND.xlapo==90.0.AND.xlopo==0.0) lrotstretch = .false.
  !
ENDIF
!
isize = SIZE(plat)/knblines
!
IF (ALLOCATED(xlon)) THEN
  IF ( isize/=SIZE(xlon) .OR. knblines/=SIZE(xlat) ) THEN
    DEALLOCATE(xlon)
    DEALLOCATE(xlat)
    DEALLOCATE(xcost)
    DEALLOCATE(xsintc)
    DEALLOCATE(xsints)
    DEALLOCATE(xcosn)
    DEALLOCATE(xsinn)
  ENDIF
ENDIF
!
IF (.NOT.ALLOCATED(xlon)) THEN
  !
  ALLOCATE(xlon(isize))
  ALLOCATE(xlat(knblines))
  !
  ALLOCATE(xcost(SIZE(xlat)))
  ALLOCATE(xsintc(SIZE(xlat)))
  ALLOCATE(xsints(SIZE(xlat)))
  ALLOCATE(xcosn(SIZE(xlon)))
  ALLOCATE(xsinn(SIZE(xlon)))
  !
  xlon(:) = angle_domain(plon(1:isize),dom='0+',unit='D') * xdr
  xcosn(:) = cos(xlon(:)-xlonp)
  xsinn(:) = sin(xlon(:)-xlonp)
  !
ENDIF
!
IF (lhook) CALL dr_hook('GET_MESH_INDEX_GAUSS_1',1,zhook_handle)
IF (lhook) CALL dr_hook('GET_MESH_INDEX_GAUSS_2',0,zhook_handle)

!
!*    3.     Projection of input points into pseudo coordinates
!            --------------------------------------------------
!
DO jj=1,knblines
  xlat(jj) = plat(jj*isize) * xdr
  xsintc(jj) = sin(xlat(jj)) * xcosp
  xsints(jj) = sin(xlat(jj)) * xsinp
  xcost(jj) = cos(xlat(jj))
ENDDO
!
IF (lrotstretch) THEN
  CALL xy_gauss(xcodil,isize,znodata,zvalue,zy,zx)
ELSE
  zx(:) = plon(:)
  zy(:) = plat(:) 
ENDIF
!
IF (lhook) CALL dr_hook('GET_MESH_INDEX_GAUSS_2',1,zhook_handle)
IF (lhook) CALL dr_hook('GET_MESH_INDEX_GAUSS_3',0,zhook_handle)
!
!-------------------------------------------------------------------------------
!
!*    5.     Localisation of the data points on (x,y) grid
!            ---------------------------------------------
!
icj(:) = 0
!
ifactx = SIZE(nfracdx)
!$OMP PARALLEL DO PRIVATE(JJ,JI,JGRID)
DO jl=1,SIZE(plat)
  !
 IF (zvalue(jl)==znodata) cycle
  !
  fracx: &  
  DO jj=1,ifactx
    !
    IF (zy(jl)>=xyinf(nfracgx(jj))) THEN
      !
      jgrid = nfracgx(jj-1)  
      !
      DO ji=nfracdx(jj-1)+1,nfracdx(jj)-1
        !
        jgrid = jgrid + nnlopa(ji-1)
        !
        IF (zy(jl)>=xyinf(jgrid)) THEN
          !
          icj(jl) = ji
          EXIT fracx
          !
        ENDIF
        !
      ENDDO
      !
      icj(jl) = nfracdx(jj)
      !
      exit fracx
      !
    ENDIF
    !
  ENDDO fracx
  !
ENDDO
!$OMP END PARALLEL DO
!
IF (lhook) CALL dr_hook('GET_MESH_INDEX_GAUSS_3',1,zhook_handle)
IF (lhook) CALL dr_hook('GET_MESH_INDEX_GAUSS_4',0,zhook_handle)
!
!$OMP PARALLEL DO PRIVATE(IGRID0,JGRID,JI,JJ)
DO jl=1,SIZE(plat)
  !
  IF (zvalue(jl)==znodata) cycle
  !
  igrid0 = 0
  IF (icj(jl)/=1) igrid0 = sum(nnlopa(1:icj(jl)-1))
  !
  !* loop on grid points: longitude
  fracy: &
  DO jj = 1,nfacty(icj(jl))
    !
    IF (zx(jl)<xxsup(igrid0+nfracdy(icj(jl),jj))) THEN
      !
      jgrid = igrid0 + nfracdy(icj(jl),jj-1)
      !
      DO ji=nfracdy(icj(jl),jj-1)+1,nfracdy(icj(jl),jj)
        !
        jgrid = jgrid + 1
        IF (zx(jl)<=xxcen(jgrid)-180. .AND. zx(jl)<xxsup(jgrid)-360.) zx(jl) = zx(jl) + 360.
    !* imput point is in this grid mesh
        IF (zx(jl)>=xxinf(jgrid) .AND. zx(jl)<xxsup(jgrid)) THEN
          !
          kindex(1,jl) = jgrid
          !
          EXIT fracy
          !
        ENDIF
        !
      ENDDO 
      !
    END IF
    !
  END DO fracy
  !
END DO
!$OMP END PARALLEL DO 
!
IF (lhook) CALL dr_hook('GET_MESH_INDEX_GAUSS_4',1,zhook_handle)
IF (lhook) CALL dr_hook('GET_MESH_INDEX_GAUSS_5',0,zhook_handle)
!
!*    6.     Localisation of the data points on in the subgrid of this mesh
!            --------------------------------------------------------------
!
IF (ksso/=0) THEN
!$OMP PARALLEL DO 
  DO jl=1,SIZE(plat)
    IF (kindex(1,jl)/=0) THEN
      kissox(1,jl) = 1 + int( float(ksso) * (zx(jl)-xxinf(kindex(1,jl)))/(xxsup(kindex(1,jl))-xxinf(kindex(1,jl))) )
      kissoy(1,jl) = 1 + int( float(ksso) * (zy(jl)-xyinf(kindex(1,jl)))/(xysup(kindex(1,jl))-xyinf(kindex(1,jl))) ) 
    ENDIF 
  ENDDO
!$OMP END PARALLEL DO
ENDIF
!
IF (lhook) CALL dr_hook('GET_MESH_INDEX_GAUSS_5',1,zhook_handle)
!-------------------------------------------------------------------------------
!
END SUBROUTINE get_mesh_index_gauss