ecearth3/sources/nemo-3.6/NEMO/OPA_SRC/OBS/obs_grd_bruteforce.h90

   SUBROUTINE obs_grd_bruteforce( kpi, kpj, kpiglo, kpjglo,       &
      &                            kldi, klei, kldj, klej,         &
      &                            kmyproc, ktotproc,              &
      &                            pglam, pgphi, pmask,            &
      &                            kobs, plam, pphi, kobsi, kobsj, &
      &                                   kproc)
      !!----------------------------------------------------------------------
      !!                ***  ROUTINE obs_grd_bruteforce ***
      !!
      !! ** Purpose : Search gridpoints to find the grid box containing
      !!              the observations
      !!
      !! ** Method  : Call to linquad
      !!
      !! ** Action  : Return kproc holding the observation and kiobsi,kobsj
      !!              valid on kproc=kmyproc processor only.
      !!   
      !! History :
      !!        !  2001-11  (N. Daget, A. Weaver)
      !!        !  2006-03  (A. Weaver) NEMOVAR migration.
      !!        !  2006-05  (K. Mogensen) Moved to to separate routine.
      !!        !  2007-10  (A. Vidard) Bug fix in wrap around checks; cleanup
      !!----------------------------------------------------------------------

      !! * Arguments
      INTEGER, INTENT(IN) :: kpi                ! Number of local longitudes
      INTEGER, INTENT(IN) :: kpj                ! Number of local latitudes
      INTEGER, INTENT(IN) :: kpiglo             ! Number of global longitudes
      INTEGER, INTENT(IN) :: kpjglo             ! Number of global latitudes
      INTEGER, INTENT(IN) :: kldi               ! Start of inner domain in i
      INTEGER, INTENT(IN) :: klei               ! End of inner domain in i
      INTEGER, INTENT(IN) :: kldj               ! Start of inner domain in j
      INTEGER, INTENT(IN) :: klej               ! End of inner domain in j
      INTEGER, INTENT(IN) :: kmyproc            ! Processor number for MPP
      INTEGER, INTENT(IN) :: ktotproc           ! Total number of processors
      REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) :: &
         & pglam,   &               ! Grid point longitude
         & pgphi,   &               ! Grid point latitude
         & pmask                    ! Grid point mask
      INTEGER,INTENT(IN) :: kobs                ! Size of the observation arrays
      REAL(KIND=wp), DIMENSION(kobs), INTENT(IN) :: &
         & plam, &                  ! Longitude of obsrvations 
         & pphi                     ! Latitude of observations
      INTEGER, DIMENSION(kobs), INTENT(OUT) :: &
         & kobsi, &                 ! I-index of observations 
         & kobsj, &                 ! J-index of observations 
         & kproc                    ! Processor number of observations
  
      !! * Local declarations
      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
         & zplam, zpphi
      REAL(wp) :: zlammax
      REAL(wp) :: zlam
      INTEGER :: ji
      INTEGER :: jj
      INTEGER :: jk
      INTEGER :: jo
      INTEGER :: jlon
      INTEGER :: jlat
      INTEGER :: joffset
      INTEGER :: jostride
      REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: &
         & zlamg, &
         & zphig, &
         & zmskg, &
         & zphitmax,&
         & zphitmin,&
         & zlamtmax,&
         & zlamtmin
      LOGICAL, DIMENSION(:,:), ALLOCATABLE :: &
         & llinvalidcell
      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
         & zlamtm,  &
         & zphitm

      !-----------------------------------------------------------------------
      ! Define grid setup for grid search
      !-----------------------------------------------------------------------
      IF (ln_grid_global) THEN
         jlon     = kpiglo
         jlat     = kpjglo
         joffset  = kmyproc
         jostride = ktotproc
      ELSE
         jlon     = kpi
         jlat     = kpj
         joffset  = 0
         jostride = 1
      ENDIF
      !-----------------------------------------------------------------------
      ! Set up data for grid search
      !-----------------------------------------------------------------------
      ALLOCATE( &
         & zlamg(jlon,jlat),             &
         & zphig(jlon,jlat),             &
         & zmskg(jlon,jlat),             &
         & zphitmax(jlon-1,jlat-1),      &
         & zphitmin(jlon-1,jlat-1),      &
         & zlamtmax(jlon-1,jlat-1),      &
         & zlamtmin(jlon-1,jlat-1),      &
         & llinvalidcell(jlon-1,jlat-1), &
         & zlamtm(4,jlon-1,jlat-1),      &
         & zphitm(4,jlon-1,jlat-1)       &
         & )
      !-----------------------------------------------------------------------
      ! Copy data to local arrays
      !-----------------------------------------------------------------------
      IF (ln_grid_global) THEN
         zlamg(:,:) = -1.e+10
         zphig(:,:) = -1.e+10
         zmskg(:,:) = -1.e+10
         DO jj = kldj, klej
            DO ji = kldi, klei
               zlamg(mig(ji),mjg(jj)) = pglam(ji,jj)
               zphig(mig(ji),mjg(jj)) = pgphi(ji,jj)
               zmskg(mig(ji),mjg(jj)) = pmask(ji,jj)
            END DO
         END DO
         CALL mpp_global_max( zlamg )
         CALL mpp_global_max( zphig )
         CALL mpp_global_max( zmskg )
      ELSE
         DO jj = 1, jlat
            DO ji = 1, jlon
               zlamg(ji,jj) = pglam(ji,jj)
               zphig(ji,jj) = pgphi(ji,jj)
               zmskg(ji,jj) = pmask(ji,jj)
            END DO
         END DO
      ENDIF
      !-----------------------------------------------------------------------
      ! Copy longitudes and latitudes
      !-----------------------------------------------------------------------
      ALLOCATE( &
         & zplam(kobs), &
         & zpphi(kobs)  &
         & )
      DO jo = 1, kobs
         zplam(jo) = plam(jo)
         zpphi(jo) = pphi(jo)
      END DO
      !-----------------------------------------------------------------------
      ! Set default values for output
      !-----------------------------------------------------------------------
      kproc(:) = -1
      kobsi(:) = -1
      kobsj(:) = -1
      !-----------------------------------------------------------------------
      ! Copy grid positions to temporary arrays and renormalize to 0 to 360.
      !-----------------------------------------------------------------------
      DO jj = 1, jlat-1
         DO ji = 1, jlon-1
            zlamtm(1,ji,jj) = zlamg(ji  ,jj  )
            zphitm(1,ji,jj) = zphig(ji  ,jj  )
            zlamtm(2,ji,jj) = zlamg(ji+1,jj  )
            zphitm(2,ji,jj) = zphig(ji+1,jj  )
            zlamtm(3,ji,jj) = zlamg(ji+1,jj+1)
            zphitm(3,ji,jj) = zphig(ji+1,jj+1)
            zlamtm(4,ji,jj) = zlamg(ji  ,jj+1)
            zphitm(4,ji,jj) = zphig(ji  ,jj+1)
         END DO
      END DO
      WHERE ( zlamtm(:,:,:) < 0.0_wp )
         zlamtm(:,:,:) = zlamtm(:,:,:) + 360.0_wp
      END WHERE
      WHERE ( zlamtm(:,:,:) > 360.0_wp )
         zlamtm(:,:,:) = zlamtm(:,:,:) - 360.0_wp
      END WHERE
      !-----------------------------------------------------------------------
      ! Handle case of the wraparound; beware, not working with orca180
      !-----------------------------------------------------------------------
      DO jj = 1, jlat-1
         DO ji = 1, jlon-1
            zlammax = MAXVAL( zlamtm(:,ji,jj) )
            WHERE (zlammax - zlamtm(:, ji, jj) > 180 ) &
               & zlamtm(:,ji,jj) = zlamtm(:,ji,jj) + 360._wp
            zphitmax(ji,jj) = MAXVAL(zphitm(:,ji,jj))
            zphitmin(ji,jj) = MINVAL(zphitm(:,ji,jj))
            zlamtmax(ji,jj) = MAXVAL(zlamtm(:,ji,jj))
            zlamtmin(ji,jj) = MINVAL(zlamtm(:,ji,jj))
         END DO
      END DO
      !-----------------------------------------------------------------------
      ! Search for boxes with only land points mark them invalid
      !-----------------------------------------------------------------------
      llinvalidcell(:,:) = .FALSE.
      DO jj = 1, jlat-1
         DO ji = 1, jlon-1
            llinvalidcell(ji,jj) =               &
               & zmskg(ji  ,jj  ) == 0.0_wp .AND. &
               & zmskg(ji+1,jj  ) == 0.0_wp .AND. &
               & zmskg(ji+1,jj+1) == 0.0_wp .AND. &
               & zmskg(ji  ,jj+1) == 0.0_wp
         END DO
      END DO

      !------------------------------------------------------------------------
      ! Master loop for grid search
      !------------------------------------------------------------------------

      DO jo = 1+joffset, kobs, jostride

         !---------------------------------------------------------------------
         ! Ensure that all observation longtiudes are between 0 and 360
         !---------------------------------------------------------------------

         IF ( zplam(jo) <   0.0_wp ) zplam(jo) = zplam(jo) + 360.0_wp
         IF ( zplam(jo) > 360.0_wp ) zplam(jo) = zplam(jo) - 360.0_wp

         !---------------------------------------------------------------------
         ! Find observations which are on within 1e-6 of a grid point
         !---------------------------------------------------------------------

         gridloop: DO jj = 1, jlat-1
            DO ji = 1, jlon-1
               IF ( ABS( zphig(ji,jj) - zpphi(jo) ) < 1e-6 )  THEN
                  zlam = zlamg(ji,jj)
                  IF ( zlam <   0.0_wp ) zlam = zlam + 360.0_wp
                  IF ( zlam > 360.0_wp ) zlam = zlam - 360.0_wp
                  IF ( ABS( zlam - zplam(jo) ) < 1e-6 ) THEN
                     IF ( llinvalidcell(ji,jj) ) THEN
                        kproc(jo) = kmyproc + 1000000
                        kobsi(jo) = ji + 1
                        kobsj(jo) = jj + 1
                        CYCLE
                     ELSE
                        kproc(jo) = kmyproc
                        kobsi(jo) = ji + 1
                        kobsj(jo) = jj + 1
                        EXIT gridloop
                     ENDIF
                  ENDIF
               ENDIF
            END DO
         END DO gridloop
         
         !---------------------------------------------------------------------
         ! Ensure that all observation longtiudes are between -180 and 180
         !---------------------------------------------------------------------

         IF ( zplam(jo) > 180 ) zplam(jo) = zplam(jo) - 360.0_wp

         !---------------------------------------------------------------------
         ! Do coordinate search using brute force.
         ! - For land points kproc is set to number of the processor + 1000000
         !   and we continue the search.
         ! - For ocean points kproc is set to the number of the processor 
         !   and we stop the search.
         !---------------------------------------------------------------------

         IF ( kproc(jo) == -1 ) THEN

            ! Normal case
            gridpoints : DO jj = 1, jlat-1
               DO ji = 1, jlon-1
                  
                  IF ( ( zplam(jo) > zlamtmax(ji,jj) ) .OR. &
                     & ( zplam(jo) < zlamtmin(ji,jj) ) ) CYCLE
                  
                  IF ( ABS( zpphi(jo) ) < 85 ) THEN
                     IF ( ( zpphi(jo) > zphitmax(ji,jj) ) .OR. &
                        & ( zpphi(jo) < zphitmin(ji,jj) ) ) CYCLE
                  ENDIF
                  
                  IF ( linquad( zplam(jo), zpphi(jo), &
                     &          zlamtm(:,ji,jj), zphitm(:,ji,jj) ) ) THEN
                     IF ( llinvalidcell(ji,jj) ) THEN
                        kproc(jo) = kmyproc + 1000000
                        kobsi(jo) = ji + 1
                        kobsj(jo) = jj + 1
                        CYCLE
                     ELSE
                        kproc(jo) = kmyproc
                        kobsi(jo) = ji + 1
                        kobsj(jo) = jj + 1
                        EXIT gridpoints
                     ENDIF
                  ENDIF
                  
               END DO
            END DO gridpoints

         ENDIF

         ! In case of failure retry for obs. longtiude + 360.
         IF ( kproc(jo) == -1 ) THEN
            gridpoints_greenwich : DO jj = 1, jlat-1
               DO ji = 1, jlon-1

                  IF ( ( zplam(jo)+360.0_wp > zlamtmax(ji,jj) ) .OR. &
                     & ( zplam(jo)+360.0_wp < zlamtmin(ji,jj) ) ) CYCLE

                  IF ( ABS( zpphi(jo) ) < 85 ) THEN
                     IF ( ( zpphi(jo) > zphitmax(ji,jj) ) .OR. &
                        & ( zpphi(jo) < zphitmin(ji,jj) ) ) CYCLE
                  ENDIF

                  IF ( linquad( zplam(jo)+360.0_wp, zpphi(jo), &
                     &          zlamtm(:,ji,jj), zphitm(:,ji,jj) ) ) THEN
                     IF ( llinvalidcell(ji,jj) ) THEN
                        kproc(jo) = kmyproc + 1000000
                        kobsi(jo) = ji + 1
                        kobsj(jo) = jj + 1
                        CYCLE
                     ELSE
                        kproc(jo) = kmyproc
                        kobsi(jo) = ji + 1
                        kobsj(jo) = jj + 1
                        EXIT gridpoints_greenwich
                     ENDIF
                  ENDIF

               END DO
            END DO gridpoints_greenwich

         ENDIF
      END DO

      !----------------------------------------------------------------------
      ! Synchronize kproc on all processors
      !----------------------------------------------------------------------
      IF ( ln_grid_global ) THEN
         CALL obs_mpp_max_integer( kproc, kobs )
         CALL obs_mpp_max_integer( kobsi, kobs )
         CALL obs_mpp_max_integer( kobsj, kobs )
      ELSE
         CALL obs_mpp_find_obs_proc( kproc, kobsi, kobsj, kobs )
      ENDIF

      WHERE( kproc(:) >= 1000000 )
         kproc(:) = kproc(:) - 1000000
      END WHERE

      DEALLOCATE( &
         & zlamg,         &
         & zphig,         &
         & zmskg,         &
         & zphitmax,      &
         & zphitmin,      &
         & zlamtmax,      &
         & zlamtmin,      &
         & llinvalidcell, &
         & zlamtm,        &
         & zphitm,        &
         & zplam,         &
         & zpphi          &
         & )

   END SUBROUTINE obs_grd_bruteforce