get_near_meshes_ign.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_near_meshes_ign(KGRID_PAR,KL,PGRID_PAR,KNEAR_NBR,KNEAR)
!     ##############################################################
!
!!**** *GET_NEAR_MESHES_IGN* get the near grid mesh indices
!!
!!    PURPOSE
!!    -------
!!
!!    METHOD
!!    ------
!!   
!!    REFERENCE
!!    ---------
!!
!!    AUTHOR
!!    ------
!!
!!    V. Masson         Meteo-France
!!
!!    MODIFICATION
!!    ------------
!!
!!    Original    03/2004
!!
!
!     Modifié par Renaud Lestrigant (02/2016) : changement complet de l'algo
!     de recherche des plus proches voisins.
!     Récupération d'un code sur Internet et adaptation locale.
!     (http://jblevins.org/mirror/amiller/qsort.f90)
!
!----------------------------------------------------------------------------
!
!*    0.     DECLARATION
!            -----------
!
USE mode_gridtype_ign
!
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)    :: kl        ! number of points
INTEGER,                         INTENT(IN)    :: knear_nbr ! number of nearest points wanted
REAL,    DIMENSION(KGRID_PAR),   INTENT(IN)    :: pgrid_par ! grid parameters
INTEGER, DIMENSION(:,:),POINTER :: knear    ! near mesh indices
!
!*    0.2    Declaration of other local variables
!            ------------------------------------
!
REAL,DIMENSION(KL)  :: zx
REAL,DIMENSION(KL)  :: zy
REAL,DIMENSION(KL)  :: zdx
REAL,DIMENSION(KL)  :: zdy
REAL,DIMENSION(KL)  :: zdis
INTEGER,DIMENSION(KL) :: indzdis

REAL :: zmaxvaldis

INTEGER :: isize, jp
REAL(KIND=JPRB) :: zhook_handle
!
!----------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('GET_NEAR_MESHES_IGN_1',0,zhook_handle)
!
 CALL get_gridtype_ign(pgrid_par,px=zx,py=zy,pdx=zdx,pdy=zdy)
!
knear(:,:) = 0
!
isize = min(knear_nbr,kl)
!
! calcul de la distance de tous les points 2 à 2
!
!
IF (lhook) CALL dr_hook('GET_NEAR_MESHES_IGN_1',1,zhook_handle)
IF (lhook) CALL dr_hook('GET_NEAR_MESHES_IGN_2',0,zhook_handle)
!
!$OMP PARALLEL DO PRIVATE(JP,ZDIS,ZMAXVALDIS,INDZDIS)
!
DO jp=1,kl
  !
    ! distance du point JP à tous les autres points
    zdis(:) = sqrt((zx(:)-zx(jp))**2 + (zy(:)-zy(jp))**2)
    ! distance maximale entre JP et les autres moints
    zmaxvaldis = 2. * maxval(zdis)
    zdis(jp) = zmaxvaldis

    CALL quick_sort(zdis, indzdis)
    knear(jp,1:isize) = indzdis(1:isize)
  !
ENDDO
!$OMP END PARALLEL DO
!
IF (lhook) CALL dr_hook('GET_NEAR_MESHES_IGN_2',1,zhook_handle)
!
 CONTAINS
!
RECURSIVE SUBROUTINE quick_sort(PLIST, KORDER)

! Quick sort routine from:
! Brainerd, W.S., Goldberg, C.H. & Adams, J.C. (1990) "Programmer's Guide to
! Fortran 90", McGraw-Hill  ISBN 0-07-000248-7, pages 149-150.
! Modified by Alan Miller to include an associated integer array which gives
! the positions of the elements in the original order.
!
IMPLICIT NONE
!
REAL, DIMENSION (:), INTENT(INOUT)  :: plist
INTEGER, DIMENSION (:), INTENT(OUT)  :: korder
!
! Local variable
INTEGER :: ji

DO ji = 1, SIZE(plist)
  korder(ji) = ji
END DO

 CALL quick_sort_1(1, SIZE(plist), plist, korder)

END SUBROUTINE quick_sort


RECURSIVE SUBROUTINE quick_sort_1(KLEFT_END, KRIGHT_END, PLIST1, KORDER1)

INTEGER, INTENT(IN) :: kleft_end, kright_end
REAL, DIMENSION (:), INTENT(INOUT)  :: plist1
INTEGER, DIMENSION (:), INTENT(INOUT)  :: korder1
!     Local variables
INTEGER             :: ji, jj, itemp
REAL                :: zref, ztemp
INTEGER, PARAMETER  :: imax_simple_sort_size = 6

IF (kright_end < kleft_end + imax_simple_sort_size) THEN
  ! Use interchange sort for small PLISTs
  CALL interchange_sort(kleft_end, kright_end, plist1, korder1)
  !
ELSE
  !
  ! Use partition ("quick") sort
  ! valeur au centre du tableau
  zref = plist1((kleft_end + kright_end)/2)
  ji = kleft_end - 1
  jj = kright_end + 1

  DO
    ! Scan PLIST from left end until element >= ZREF is found
    DO
      ji = ji + 1
      IF (plist1(ji) >= zref) EXIT
    END DO
    ! Scan PLIST from right end until element <= ZREF is found
    DO
      jj = jj - 1
      IF (plist1(jj) <= zref) EXIT
    END DO


    IF (ji < jj) THEN
      ! Swap two out-of-order elements
      ztemp = plist1(ji)
      plist1(ji) = plist1(jj)
      plist1(jj) = ztemp
      itemp = korder1(ji)
      korder1(ji) = korder1(jj)
      korder1(jj) = itemp
    ELSE IF (ji == jj) THEN
      ji = ji + 1
      EXIT
    ELSE
      EXIT
    END IF
  END DO

  IF (kleft_end < jj) CALL quick_sort_1(kleft_end, jj, plist1, korder1)
  IF (ji < kright_end) CALL quick_sort_1(ji, kright_end,plist1,korder1)
END IF

END SUBROUTINE quick_sort_1


SUBROUTINE interchange_sort(KLEFT_END, KRIGHT_END, PLIST2, KORDER2)

INTEGER, INTENT(IN) :: kleft_end, kright_end
REAL, DIMENSION (:), INTENT(INOUT)  :: plist2
INTEGER, DIMENSION (:), INTENT(INOUT)  :: korder2
!     Local variables
INTEGER             :: ji, jj, itemp
REAL                :: ztemp

! boucle sur tous les points
DO ji = kleft_end, kright_end - 1
  !
  ! boucle sur les points suivants le point JI
  DO jj = ji+1, kright_end
    !
    ! si la distance de JI au point est plus grande que celle de JJ
    IF (plist2(ji) > plist2(jj)) THEN
      ! distance de JI au point (la plus grande)
      ztemp = plist2(ji)
      ! le point JJ est déplacé à l'indice JI dans le tableau 
      plist2(ji) = plist2(jj)
      ! le point JI est déplacé à l'indice JJ dans le tableau
      plist2(jj) = ztemp
      ! indice du point JI dans le tableau
      itemp = korder2(ji)
      ! l'indice du point JJ est mis à la place JI
      korder2(ji) = korder2(jj)
      ! l'indice du point JI est mis à la place JJ
      korder2(jj) = itemp
    END IF
    !
  END DO
  !
END DO

END SUBROUTINE interchange_sort
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE get_near_meshes_ign