source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/OBS/obs_grd_bruteforce.h90 @ 2358

Last change on this file since 2358 was 2358, checked in by rblod, 10 years ago

Changes to be able to compile v3_3_beta with key_agrif,see ticket #753 ; just compilation fixes, I was to scared to try to run AGRIF

File size: 13.9 KB
Line 
1   SUBROUTINE obs_grd_bruteforce( kpi, kpj, kpiglo, kpjglo,       &
2      &                            kldi, klei, kldj, klej,         &
3      &                            kmyproc, ktotproc,              &
4      &                            pglam, pgphi, pmask,            &
5      &                            kobs, plam, pphi, kobsi, kobsj, &
6      &                                   kproc)
7      !!----------------------------------------------------------------------
8      !!                ***  ROUTINE obs_grd_bruteforce ***
9      !!
10      !! ** Purpose : Search gridpoints to find the grid box containing
11      !!              the observations
12      !!
13      !! ** Method  : Call to linquad
14      !!
15      !! ** Action  : Return kproc holding the observation and kiobsi,kobsj
16      !!              valid on kproc=kmyproc processor only.
17      !!   
18      !! History :
19      !!        !  2001-11  (N. Daget, A. Weaver)
20      !!        !  2006-03  (A. Weaver) NEMOVAR migration.
21      !!        !  2006-05  (K. Mogensen) Moved to to separate routine.
22      !!        !  2007-10  (A. Vidard) Bug fix in wrap around checks; cleanup
23      !!----------------------------------------------------------------------
24
25      !! * Arguments
26      INTEGER, INTENT(IN) :: kpi                ! Number of local longitudes
27      INTEGER, INTENT(IN) :: kpj                ! Number of local latitudes
28      INTEGER, INTENT(IN) :: kpiglo             ! Number of global longitudes
29      INTEGER, INTENT(IN) :: kpjglo             ! Number of global latitudes
30      INTEGER, INTENT(IN) :: kldi               ! Start of inner domain in i
31      INTEGER, INTENT(IN) :: klei               ! End of inner domain in i
32      INTEGER, INTENT(IN) :: kldj               ! Start of inner domain in j
33      INTEGER, INTENT(IN) :: klej               ! End of inner domain in j
34      INTEGER, INTENT(IN) :: kmyproc            ! Processor number for MPP
35      INTEGER, INTENT(IN) :: ktotproc           ! Total number of processors
36      REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) :: &
37         & pglam,   &               ! Grid point longitude
38         & pgphi,   &               ! Grid point latitude
39         & pmask                    ! Grid point mask
40      INTEGER,INTENT(IN) :: kobs                ! Size of the observation arrays
41      REAL(KIND=wp), DIMENSION(kobs), INTENT(IN) :: &
42         & plam, &                  ! Longitude of obsrvations
43         & pphi                     ! Latitude of observations
44      INTEGER, DIMENSION(kobs), INTENT(OUT) :: &
45         & kobsi, &                 ! I-index of observations
46         & kobsj, &                 ! J-index of observations
47         & kproc                    ! Processor number of observations
48 
49      !! * Local declarations
50      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
51         & zplam, zpphi
52      REAL(wp) :: zlammax
53      REAL(wp) :: zlam
54      INTEGER :: ji
55      INTEGER :: jj
56      INTEGER :: jk
57      INTEGER :: jo
58      INTEGER :: jlon
59      INTEGER :: jlat
60      INTEGER :: joffset
61      INTEGER :: jostride
62      REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: &
63         & zlamg, &
64         & zphig, &
65         & zmskg, &
66         & zphitmax,&
67         & zphitmin,&
68         & zlamtmax,&
69         & zlamtmin
70      LOGICAL, DIMENSION(:,:), ALLOCATABLE :: &
71         & llinvalidcell
72      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
73         & zlamtm,  &
74         & zphitm
75
76      !-----------------------------------------------------------------------
77      ! Define grid setup for grid search
78      !-----------------------------------------------------------------------
79      IF (ln_grid_global) THEN
80         jlon     = kpiglo
81         jlat     = kpjglo
82         joffset  = kmyproc
83         jostride = ktotproc
84      ELSE
85         jlon     = kpi
86         jlat     = kpj
87         joffset  = 0
88         jostride = 1
89      ENDIF
90      !-----------------------------------------------------------------------
91      ! Set up data for grid search
92      !-----------------------------------------------------------------------
93      ALLOCATE( &
94         & zlamg(jlon,jlat),             &
95         & zphig(jlon,jlat),             &
96         & zmskg(jlon,jlat),             &
97         & zphitmax(jlon-1,jlat-1),      &
98         & zphitmin(jlon-1,jlat-1),      &
99         & zlamtmax(jlon-1,jlat-1),      &
100         & zlamtmin(jlon-1,jlat-1),      &
101         & llinvalidcell(jlon-1,jlat-1), &
102         & zlamtm(4,jlon-1,jlat-1),      &
103         & zphitm(4,jlon-1,jlat-1)       &
104         & )
105      !-----------------------------------------------------------------------
106      ! Copy data to local arrays
107      !-----------------------------------------------------------------------
108      IF (ln_grid_global) THEN
109         zlamg(:,:) = -1.e+10
110         zphig(:,:) = -1.e+10
111         zmskg(:,:) = -1.e+10
112         DO jj = kldj, klej
113            DO ji = kldi, klei
114               zlamg(mig(ji),mjg(jj)) = pglam(ji,jj)
115               zphig(mig(ji),mjg(jj)) = pgphi(ji,jj)
116               zmskg(mig(ji),mjg(jj)) = pmask(ji,jj)
117            END DO
118         END DO
119         CALL mpp_global_max( zlamg )
120         CALL mpp_global_max( zphig )
121         CALL mpp_global_max( zmskg )
122      ELSE
123         DO jj = 1, jlat
124            DO ji = 1, jlon
125               zlamg(ji,jj) = pglam(ji,jj)
126               zphig(ji,jj) = pgphi(ji,jj)
127               zmskg(ji,jj) = pmask(ji,jj)
128            END DO
129         END DO
130      ENDIF
131      !-----------------------------------------------------------------------
132      ! Copy longitudes and latitudes
133      !-----------------------------------------------------------------------
134      ALLOCATE( &
135         & zplam(kobs), &
136         & zpphi(kobs)  &
137         & )
138      DO jo = 1, kobs
139         zplam(jo) = plam(jo)
140         zpphi(jo) = pphi(jo)
141      END DO
142      !-----------------------------------------------------------------------
143      ! Set default values for output
144      !-----------------------------------------------------------------------
145      kproc(:) = -1
146      kobsi(:) = -1
147      kobsj(:) = -1
148      !-----------------------------------------------------------------------
149      ! Copy grid positions to temporary arrays and renormalize to 0 to 360.
150      !-----------------------------------------------------------------------
151      DO jj = 1, jlat-1
152         DO ji = 1, jlon-1
153            zlamtm(1,ji,jj) = zlamg(ji  ,jj  )
154            zphitm(1,ji,jj) = zphig(ji  ,jj  )
155            zlamtm(2,ji,jj) = zlamg(ji+1,jj  )
156            zphitm(2,ji,jj) = zphig(ji+1,jj  )
157            zlamtm(3,ji,jj) = zlamg(ji+1,jj+1)
158            zphitm(3,ji,jj) = zphig(ji+1,jj+1)
159            zlamtm(4,ji,jj) = zlamg(ji  ,jj+1)
160            zphitm(4,ji,jj) = zphig(ji  ,jj+1)
161         END DO
162      END DO
163      WHERE ( zlamtm(:,:,:) < 0.0_wp )
164         zlamtm(:,:,:) = zlamtm(:,:,:) + 360.0_wp
165      END WHERE
166      WHERE ( zlamtm(:,:,:) > 360.0_wp )
167         zlamtm(:,:,:) = zlamtm(:,:,:) - 360.0_wp
168      END WHERE
169      !-----------------------------------------------------------------------
170      ! Handle case of the wraparound; beware, not working with orca180
171      !-----------------------------------------------------------------------
172      DO jj = 1, jlat-1
173         DO ji = 1, jlon-1
174            zlammax = MAXVAL( zlamtm(:,ji,jj) )
175            WHERE (zlammax - zlamtm(:, ji, jj) > 180 ) &
176               & zlamtm(:,ji,jj) = zlamtm(:,ji,jj) + 360._wp
177            zphitmax(ji,jj) = MAXVAL(zphitm(:,ji,jj))
178            zphitmin(ji,jj) = MINVAL(zphitm(:,ji,jj))
179            zlamtmax(ji,jj) = MAXVAL(zlamtm(:,ji,jj))
180            zlamtmin(ji,jj) = MINVAL(zlamtm(:,ji,jj))
181         END DO
182      END DO
183      !-----------------------------------------------------------------------
184      ! Search for boxes with only land points mark them invalid
185      !-----------------------------------------------------------------------
186      llinvalidcell(:,:) = .FALSE.
187      DO jj = 1, jlat-1
188         DO ji = 1, jlon-1
189            llinvalidcell(ji,jj) =               &
190               & zmskg(ji  ,jj  ) == 0.0_wp .AND. &
191               & zmskg(ji+1,jj  ) == 0.0_wp .AND. &
192               & zmskg(ji+1,jj+1) == 0.0_wp .AND. &
193               & zmskg(ji  ,jj+1) == 0.0_wp
194         END DO
195      END DO
196
197      !------------------------------------------------------------------------
198      ! Master loop for grid search
199      !------------------------------------------------------------------------
200
201      DO jo = 1+joffset, kobs, jostride
202
203         !---------------------------------------------------------------------
204         ! Ensure that all observation longtiudes are between 0 and 360
205         !---------------------------------------------------------------------
206
207         IF ( zplam(jo) <   0.0_wp ) zplam(jo) = zplam(jo) + 360.0_wp
208         IF ( zplam(jo) > 360.0_wp ) zplam(jo) = zplam(jo) - 360.0_wp
209
210         !---------------------------------------------------------------------
211         ! Find observations which are on within 1e-6 of a grid point
212         !---------------------------------------------------------------------
213
214         gridloop: DO jj = 1, jlat-1
215            DO ji = 1, jlon-1
216               IF ( ABS( zphig(ji,jj) - zpphi(jo) ) < 1e-6 )  THEN
217                  zlam = zlamg(ji,jj)
218                  IF ( zlam <   0.0_wp ) zlam = zlam + 360.0_wp
219                  IF ( zlam > 360.0_wp ) zlam = zlam - 360.0_wp
220                  IF ( ABS( zlam - zplam(jo) ) < 1e-6 ) THEN
221                     IF ( llinvalidcell(ji,jj) ) THEN
222                        kproc(jo) = kmyproc + 1000000
223                        kobsi(jo) = ji + 1
224                        kobsj(jo) = jj + 1
225                        CYCLE
226                     ELSE
227                        kproc(jo) = kmyproc
228                        kobsi(jo) = ji + 1
229                        kobsj(jo) = jj + 1
230                        EXIT gridloop
231                     ENDIF
232                  ENDIF
233               ENDIF
234            END DO
235         END DO gridloop
236         
237         !---------------------------------------------------------------------
238         ! Ensure that all observation longtiudes are between -180 and 180
239         !---------------------------------------------------------------------
240
241         IF ( zplam(jo) > 180 ) zplam(jo) = zplam(jo) - 360.0_wp
242
243         !---------------------------------------------------------------------
244         ! Do coordinate search using brute force.
245         ! - For land points kproc is set to number of the processor + 1000000
246         !   and we continue the search.
247         ! - For ocean points kproc is set to the number of the processor
248         !   and we stop the search.
249         !---------------------------------------------------------------------
250
251         IF ( kproc(jo) == -1 ) THEN
252
253            ! Normal case
254            gridpoints : DO jj = 1, jlat-1
255               DO ji = 1, jlon-1
256                 
257                  IF ( ( zplam(jo) > zlamtmax(ji,jj) ) .OR. &
258                     & ( zplam(jo) < zlamtmin(ji,jj) ) ) CYCLE
259                 
260                  IF ( ABS( zpphi(jo) ) < 85 ) THEN
261                     IF ( ( zpphi(jo) > zphitmax(ji,jj) ) .OR. &
262                        & ( zpphi(jo) < zphitmin(ji,jj) ) ) CYCLE
263                  ENDIF
264                 
265                  IF ( linquad( zplam(jo), zpphi(jo), &
266                     &          zlamtm(:,ji,jj), zphitm(:,ji,jj) ) ) THEN
267                     IF ( llinvalidcell(ji,jj) ) THEN
268                        kproc(jo) = kmyproc + 1000000
269                        kobsi(jo) = ji + 1
270                        kobsj(jo) = jj + 1
271                        CYCLE
272                     ELSE
273                        kproc(jo) = kmyproc
274                        kobsi(jo) = ji + 1
275                        kobsj(jo) = jj + 1
276                        EXIT gridpoints
277                     ENDIF
278                  ENDIF
279                 
280               END DO
281            END DO gridpoints
282
283         ENDIF
284
285         ! In case of failure retry for obs. longtiude + 360.
286         IF ( kproc(jo) == -1 ) THEN
287            gridpoints_greenwich : DO jj = 1, jlat-1
288               DO ji = 1, jlon-1
289
290                  IF ( ( zplam(jo)+360.0_wp > zlamtmax(ji,jj) ) .OR. &
291                     & ( zplam(jo)+360.0_wp < zlamtmin(ji,jj) ) ) CYCLE
292
293                  IF ( ABS( zpphi(jo) ) < 85 ) THEN
294                     IF ( ( zpphi(jo) > zphitmax(ji,jj) ) .OR. &
295                        & ( zpphi(jo) < zphitmin(ji,jj) ) ) CYCLE
296                  ENDIF
297
298                  IF ( linquad( zplam(jo)+360.0_wp, zpphi(jo), &
299                     &          zlamtm(:,ji,jj), zphitm(:,ji,jj) ) ) THEN
300                     IF ( llinvalidcell(ji,jj) ) THEN
301                        kproc(jo) = kmyproc + 1000000
302                        kobsi(jo) = ji + 1
303                        kobsj(jo) = jj + 1
304                        CYCLE
305                     ELSE
306                        kproc(jo) = kmyproc
307                        kobsi(jo) = ji + 1
308                        kobsj(jo) = jj + 1
309                        EXIT gridpoints_greenwich
310                     ENDIF
311                  ENDIF
312
313               END DO
314            END DO gridpoints_greenwich
315
316         ENDIF
317      END DO
318
319      !----------------------------------------------------------------------
320      ! Synchronize kproc on all processors
321      !----------------------------------------------------------------------
322      IF ( ln_grid_global ) THEN
323         CALL obs_mpp_max_integer( kproc, kobs )
324         CALL obs_mpp_max_integer( kobsi, kobs )
325         CALL obs_mpp_max_integer( kobsj, kobs )
326      ELSE
327         CALL obs_mpp_find_obs_proc( kproc, kobsi, kobsj, kobs )
328      ENDIF
329
330      WHERE( kproc(:) >= 1000000 )
331         kproc(:) = kproc(:) - 1000000
332      END WHERE
333
334      DEALLOCATE( &
335         & zlamg,         &
336         & zphig,         &
337         & zmskg,         &
338         & zphitmax,      &
339         & zphitmin,      &
340         & zlamtmax,      &
341         & zlamtmin,      &
342         & llinvalidcell, &
343         & zlamtm,        &
344         & zphitm,        &
345         & zplam,         &
346         & zpphi          &
347         & )
348
349   END SUBROUTINE obs_grd_bruteforce
Note: See TracBrowser for help on using the repository browser.