MODULE domngb
   !!======================================================================
   !!                       ***  MODULE  domngb  ***
   !! Grid search:  find the closest grid point from a given on/lat position
   !!======================================================================
   !! History : 3.2  !  2009-11  (S. Masson)  Original code
   !!----------------------------------------------------------------------

   !!----------------------------------------------------------------------
   !!   dom_ngb       : find the closest grid point from a given lon/lat position
   !!----------------------------------------------------------------------
   USE dom_oce        ! ocean space and time domain
   USE in_out_manager ! I/O manager
   USE lib_mpp        ! for mppsum
   USE wrk_nemo       ! Memory allocation
   USE timing         ! Timing

   IMPLICIT NONE
   PRIVATE

   PUBLIC   dom_ngb   ! routine called in iom.F90 module

   !!----------------------------------------------------------------------
   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
   !! $Id$ 
   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------
CONTAINS

   SUBROUTINE dom_ngb( plon, plat, kii, kjj, cdgrid, kkk )
      !!----------------------------------------------------------------------
      !!                    ***  ROUTINE dom_ngb  ***
      !!
      !! ** Purpose :   find the closest grid point from a given lon/lat position
      !!
      !! ** Method  :   look for minimum distance in cylindrical projection 
      !!                -> not good if located at too high latitude...
      !!----------------------------------------------------------------------
      REAL(wp)        , INTENT(in   ) ::   plon, plat   ! longitude,latitude of the point
      INTEGER         , INTENT(  out) ::   kii, kjj     ! i-,j-index of the closes grid point
      INTEGER         , INTENT(in   ), OPTIONAL :: kkk  ! k-index of the mask level used
      CHARACTER(len=1), INTENT(in   ) ::   cdgrid       ! grid name 'T', 'U', 'V', 'W'
      !
      INTEGER :: ik, jj, ji         ! working level
      INTEGER , DIMENSION(2) ::   iloc
      REAL(wp)               ::   zlon, zmini
      REAL(wp), POINTER, DIMENSION(:,:) ::  zglam, zgphi, zmask, zdist
      !!--------------------------------------------------------------------
      !
      IF( nn_timing == 1 )  CALL timing_start('dom_ngb')
      !
      CALL wrk_alloc( jpi,jpj,   zglam, zgphi, zmask, zdist )
      !
!$OMP PARALLEL DO schedule(static) private(jj, ji)
      DO jj = 1, jpj
         DO ji = 1, jpi
            zmask(ji,jj) = 0._wp
         END DO
      END DO
      ik = 1
      IF ( PRESENT(kkk) ) ik=kkk
      SELECT CASE( cdgrid )
      CASE( 'U' )
!$OMP PARALLEL DO schedule(static) private(jj, ji)
      DO jj = 1, jpj
         DO ji = 1, jpi
            zglam(ji,jj) = glamu(ji,jj) ; zgphi(ji,jj) = gphiu(ji,jj)
         END DO
      END DO
      zmask(nldi:nlei,nldj:nlej) = umask(nldi:nlei,nldj:nlej,ik)
      CASE( 'V' )  
!$OMP PARALLEL DO schedule(static) private(jj, ji)
      DO jj = 1, jpj
         DO ji = 1, jpi
            zglam(ji,jj) = glamv(ji,jj) ; zgphi(ji,jj) = gphiv(ji,jj)
         END DO
      END DO
      zmask(nldi:nlei,nldj:nlej) = vmask(nldi:nlei,nldj:nlej,ik)
      CASE( 'F' )  
!$OMP PARALLEL DO schedule(static) private(jj, ji)
      DO jj = 1, jpj
         DO ji = 1, jpi
            zglam(ji,jj) = glamf(ji,jj) ; zgphi(ji,jj) = gphif(ji,jj)
         END DO
      END DO
      zmask(nldi:nlei,nldj:nlej) = fmask(nldi:nlei,nldj:nlej,ik)
      CASE DEFAULT 
!$OMP PARALLEL DO schedule(static) private(jj, ji)
      DO jj = 1, jpj
         DO ji = 1, jpi
            zglam(ji,jj) = glamt(ji,jj) ; zgphi(ji,jj) = gphit(ji,jj)
         END DO
      END DO
      zmask(nldi:nlei,nldj:nlej) = tmask(nldi:nlei,nldj:nlej,ik)
      END SELECT

      IF (jphgr_msh /= 2 .AND. jphgr_msh /= 3) THEN
         zlon       = MOD( plon       + 720., 360. )                                     ! plon between    0 and 360
!$OMP PARALLEL DO schedule(static) private(jj, ji, zglam, zlon)
         DO jj = 1, jpj
            DO ji = 1, jpi
               zglam(ji,jj) = MOD( zglam(ji,jj) + 720., 360. )                                     ! glam between    0 and 360
               IF( zlon > 270. )   zlon = zlon - 360.                                          ! zlon between  -90 and 270
               IF( zlon <  90. ) THEN
                 IF( zglam(ji,jj) > 180. ) zglam(ji,jj) = zglam(ji,jj) - 360.   ! glam between -180 and 180
               END IF
               zglam(ji,jj) = zglam(ji,jj) - zlon
            END DO
         END DO
      ELSE
!$OMP PARALLEL DO schedule(static) private(jj, ji)
         DO jj = 1, jpj
            DO ji = 1, jpi
               zglam(ji,jj) = zglam(ji,jj) - plon
            END DO
         END DO
      END IF
!$OMP PARALLEL DO schedule(static) private(jj, ji)
      DO jj = 1, jpj
         DO ji = 1, jpi
            zgphi(ji,jj) = zgphi(ji,jj) - plat
            zdist(ji,jj) = zglam(ji,jj) * zglam(ji,jj) + zgphi(ji,jj) * zgphi(ji,jj)
         END DO
      END DO
      IF( lk_mpp ) THEN  
         CALL mpp_minloc( zdist(:,:), zmask, zmini, kii, kjj)
      ELSE
         iloc(:) = MINLOC( zdist(:,:), mask = zmask(:,:) == 1.e0 )
         kii = iloc(1) + nimpp - 1
         kjj = iloc(2) + njmpp - 1
      ENDIF
      !
      CALL wrk_dealloc( jpi,jpj,   zglam, zgphi, zmask, zdist )
      !
      IF( nn_timing == 1 )  CALL timing_stop('dom_ngb')
      !
   END SUBROUTINE dom_ngb

   !!======================================================================
END MODULE domngb