Forked from
NEMO Workspace / Nemo
968 commits behind, 63 commits ahead of the upstream repository.
-
Sebastien Masson authorede0e6ed97
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
obs_grd_bruteforce.h90 13.90 KiB
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,nn_hls),mjg(jj,nn_hls)) = pglam(ji,jj)
zphig(mig(ji,nn_hls),mjg(jj,nn_hls)) = pgphi(ji,jj)
zmskg(mig(ji,nn_hls),mjg(jj,nn_hls)) = 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, 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