New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
obs_grid.F90 in branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/OBS/obs_grid.F90 @ 2281

Last change on this file since 2281 was 2281, checked in by smasson, 14 years ago

set proper svn properties to all files...
  • Property svn:keywords set to Id
File size: 46.7 KB
Line 
1MODULE obs_grid
2   !!======================================================================
3   !!                       ***  MODULE  obs_grid  ***
4   !! Observation diagnostics: Various tools for grid searching etc.
5   !!======================================================================
6   !!----------------------------------------------------------------------
7   !!   obs_grid_search   : Find i,j on the ORCA grid from lat,lon
8   !!   obs_level_search  : Find level from depth
9   !!   obs_zlevel_search : Find depth level from observed depth
10   !!   obs_tlevel_search : Find temperature level from observed temp
11   !!   obs_rlevel_search : Find density level from observed density
12   !!----------------------------------------------------------------------
13   !! * Modules used   
14   USE par_kind, ONLY : &          ! Precision variables
15      & wp 
16   USE par_oce, ONLY :  &          ! Ocean parameters
17      & jpk,     &
18      & jpni,    &
19      & jpnj,    &
20      & jpnij
21   USE dom_oce                     ! Ocean space and time domain variables
22   USE obs_mpp, ONLY : &           ! MPP support routines for observation diagnostics
23      & obs_mpp_find_obs_proc, &
24      & mpp_global_max,        &
25      & obs_mpp_max_integer
26   USE phycst, ONLY : &            ! Physical constants
27      & rad
28   USE obs_utils, ONLY : &         ! Observation operator utility functions
29      & grt_cir_dis, &
30      & chkerr
31   USE in_out_manager              ! Printing support
32   USE netcdf
33   USE obs_const, ONLY :      &
34      & obfillflt                  ! Fillvalue   
35   
36   IMPLICIT NONE
37
38   !! * Routine accessibility
39   PUBLIC  obs_grid_setup,      & ! Setup grid searching
40      &    obs_grid_search,     & ! Find i, j on the ORCA grid from lat, lon
41      &    obs_grid_deallocate, & ! Deallocate the look up table
42      &    obs_level_search       ! Find level from depth
43
44   PRIVATE linquad,                    & ! Determine whether a point lies within a cell
45      &    maxdist,                    & ! Find the maximum distance between 2 pts in a cell
46      &    obs_grid_search_bruteforce, & ! Find i, j on the ORCA grid from lat, lon
47      &    obs_grid_search_lookup        ! Find i, j on the ORCA grid from lat, lon quicker
48
49   !!* Module variables
50
51   !! Default values
52   REAL, PUBLIC :: grid_search_res = 0.5    ! Resolution of grid
53   INTEGER, PRIVATE :: gsearch_nlons_def    ! Num of longitudes
54   INTEGER, PRIVATE :: gsearch_nlats_def    ! Num of latitudes
55   REAL(wp), PRIVATE :: gsearch_lonmin_def  ! Min longitude
56   REAL(wp), PRIVATE :: gsearch_latmin_def  ! Min latitude
57   REAL(wp), PRIVATE :: gsearch_dlon_def    ! Lon spacing
58   REAL(wp), PRIVATE :: gsearch_dlat_def    ! Lat spacing
59   !! Variable versions
60   INTEGER, PRIVATE :: nlons     ! Num of longitudes
61   INTEGER, PRIVATE :: nlats     ! Num of latitudes
62   REAL(wp), PRIVATE :: lonmin ! Min longitude
63   REAL(wp), PRIVATE :: latmin ! Min latitude
64   REAL(wp), PRIVATE :: dlon     ! Lon spacing
65   REAL(wp), PRIVATE :: dlat     ! Lat spacing
66   
67   INTEGER, PRIVATE :: maxxdiff, maxydiff ! Max diffs between model points
68   INTEGER, PRIVATE :: limxdiff, limydiff
69   
70   ! Data storage
71   REAL(wp), PRIVATE, DIMENSION(:,:), ALLOCATABLE :: &
72      & lons,  &
73      & lats
74   INTEGER, PRIVATE, DIMENSION(:,:), ALLOCATABLE :: &
75      & ixpos, &
76      & iypos, &
77      & iproc   
78
79   ! Switches
80   LOGICAL, PUBLIC :: ln_grid_search_lookup  ! Use lookup table to speed up grid search
81   LOGICAL, PUBLIC :: ln_grid_global         ! Use global distribution of observations
82   CHARACTER(LEN=44), PUBLIC :: &
83      & grid_search_file    ! file name head for grid search lookup
84
85CONTAINS
86
87   SUBROUTINE obs_grid_search( kobsin, plam, pphi, kobsi, kobsj, kproc, &
88      &                        cdgrid )
89      !!----------------------------------------------------------------------
90      !!                ***  ROUTINE obs_grid_search ***
91      !!
92      !! ** Purpose : Search local gridpoints to find the grid box containing
93      !!              the observations calls either
94      !!              obs_grid_search_bruteforce - the original brute force search
95      !!                     or
96      !!              obs_grid_search_lookup - uses a lookup table to do a fast
97      !!search
98      !!History :
99      !!        !  2007-12  (D. Lea)
100      !!------------------------------------------------------------------------
101
102      !! * Arguments
103      INTEGER :: &
104         & kobsin                     ! Size of the observation arrays
105      REAL(KIND=wp), DIMENSION(kobsin), INTENT(IN) :: &
106         & plam, &                  ! Longitude of obsrvations
107         & pphi                     ! Latitude of observations
108      INTEGER, DIMENSION(kobsin), INTENT(OUT) :: &
109         & kobsi, &                 ! I-index of observations
110         & kobsj, &                 ! J-index of observations
111         & kproc                    ! Processor number of observations
112      CHARACTER(LEN=1) :: &
113         & cdgrid                   ! Grid to search
114
115      IF(kobsin > 0) THEN
116
117         IF ( ln_grid_search_lookup .AND. ( cdgrid == 'T' ) ) THEN
118            CALL obs_grid_search_lookup( kobsin, plam, pphi, &
119               &                         kobsi, kobsj, kproc )
120         ELSE
121            IF ( cdgrid == 'T' ) THEN
122               CALL obs_grid_search_bruteforce( jpi, jpj, jpiglo, jpjglo, &
123                  &                             nldi, nlei,nldj,  nlej,   &
124                  &                             nproc, jpnij,             &
125                  &                             glamt, gphit, tmask,      &
126                  &                             kobsin, plam, pphi,       &
127                  &                             kobsi, kobsj, kproc )
128            ELSEIF ( cdgrid == 'U' ) THEN
129               CALL obs_grid_search_bruteforce( jpi, jpj, jpiglo, jpjglo, &
130                  &                             nldi, nlei,nldj,  nlej,   &
131                  &                             nproc, jpnij,             &
132                  &                             glamu, gphiu, umask,      &
133                  &                             kobsin, plam, pphi,       &
134                  &                             kobsi, kobsj, kproc )
135            ELSEIF ( cdgrid == 'V' ) THEN
136               CALL obs_grid_search_bruteforce( jpi, jpj, jpiglo, jpjglo, &
137                  &                             nldi, nlei,nldj,  nlej,   &
138                  &                             nproc, jpnij,             &
139                  &                             glamv, gphiv, vmask,      &
140                  &                             kobsin, plam, pphi,       &
141                  &                             kobsi, kobsj, kproc )
142            ELSEIF ( cdgrid == 'F' ) THEN
143               CALL obs_grid_search_bruteforce( jpi, jpj, jpiglo, jpjglo, &
144                  &                             nldi, nlei,nldj,  nlej,   &
145                  &                             nproc, jpnij,             &
146                  &                             glamf, gphif, fmask,      &
147                  &                             kobsin, plam, pphi,       &
148                  &                             kobsi, kobsj, kproc )
149            ELSE
150               CALL ctl_stop( 'Grid not supported' )
151            ENDIF
152         ENDIF
153         
154      ENDIF
155
156   END SUBROUTINE obs_grid_search
157
158#include "obs_grid_search_bruteforce.h90"
159   
160   SUBROUTINE obs_grid_search_lookup( kobs, plam, pphi, kobsi, kobsj, kproc )
161      !!----------------------------------------------------------------------
162      !!                ***  ROUTINE obs_grid_search_lookup ***
163      !!
164      !! ** Purpose : Search local gridpoints to find the grid box containing
165      !!              the observations (much faster then obs_grid_search_bruteforce)
166      !!
167      !! ** Method  : Call to linquad
168      !!
169      !! ** Action  : Return kproc holding the observation and kiobsi,kobsj
170      !!              valid on kproc=nproc processor only.
171      !!   
172      !! History :
173      !!        !  2007-12 (D. Lea) new routine based on obs_grid_search
174      !!! updated with fixes from new version of obs_grid_search_bruteforce
175      !!! speeded up where points are not near a "difficult" region like an edge
176      !!----------------------------------------------------------------------
177
178      !! * Arguments
179      INTEGER :: kobs                     ! Size of the observation arrays
180      REAL(KIND=wp), DIMENSION(kobs), INTENT(IN) :: &
181         & plam, &                  ! Longitude of obsrvations
182         & pphi                     ! Latitude of observations
183      INTEGER, DIMENSION(kobs), INTENT(OUT) :: &
184         & kobsi, &                 ! I-index of observations
185         & kobsj, &                 ! J-index of observations
186         & kproc                    ! Processor number of observations
187 
188      !! * Local declarations
189      REAL(KIND=wp), DIMENSION(:), ALLOCATABLE :: &
190         & zplam
191      REAL(wp) :: zlammax
192      REAL(wp) :: zlam
193      INTEGER :: ji
194      INTEGER :: jj
195      INTEGER :: jk
196      INTEGER :: jo
197      INTEGER :: isx
198      INTEGER :: isy
199      INTEGER :: jimin
200      INTEGER :: jimax
201      INTEGER :: jjmin
202      INTEGER :: jjmax
203      INTEGER :: jojimin
204      INTEGER :: jojimax
205      INTEGER :: jojjmin
206      INTEGER :: jojjmax
207      INTEGER :: ipx1
208      INTEGER :: ipy1
209      INTEGER :: ip
210      INTEGER :: jp
211      INTEGER :: ipx
212      INTEGER :: ipy
213      INTEGER :: ipmx
214      INTEGER :: jlon
215      INTEGER :: jlat
216      INTEGER :: joffset
217      INTEGER :: jostride
218      REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: &
219         & zlamg, &
220         & zphig, &
221         & zmskg, &
222         & zphitmax,&
223         & zphitmin,&
224         & zlamtmax,&
225         & zlamtmin
226      LOGICAL, DIMENSION(:,:), ALLOCATABLE :: &
227         & llinvalidcell
228      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
229         & zlamtm,  &
230         & zphitm
231      LOGICAL :: llfourflag
232      INTEGER :: ifourflagcountt
233      INTEGER :: ifourflagcountf
234      INTEGER, DIMENSION(5) :: ifourflagcountr
235
236      !-----------------------------------------------------------------------
237      ! Define grid for grid search
238      !-----------------------------------------------------------------------
239      IF (ln_grid_global) THEN
240         jlon     = jpiglo
241         jlat     = jpjglo
242         joffset  = nproc
243         jostride = jpnij
244      ELSE
245         jlon     = jpi
246         jlat     = jpj
247         joffset  = 0
248         jostride = 1
249      ENDIF
250      !-----------------------------------------------------------------------
251      ! Set up data for grid search
252      !-----------------------------------------------------------------------
253      ALLOCATE( &
254         & zlamg(jlon,jlat),             &
255         & zphig(jlon,jlat),             &
256         & zmskg(jlon,jlat),             &
257         & zphitmax(jlon-1,jlat-1),      &
258         & zphitmin(jlon-1,jlat-1),      &
259         & zlamtmax(jlon-1,jlat-1),      &
260         & zlamtmin(jlon-1,jlat-1),      &
261         & llinvalidcell(jlon-1,jlat-1), &
262         & zlamtm(4,jlon-1,jlat-1),      &
263         & zphitm(4,jlon-1,jlat-1)       &
264         & )
265      !-----------------------------------------------------------------------
266      ! Copy data to local arrays
267      !-----------------------------------------------------------------------
268      IF (ln_grid_global) THEN
269         zlamg(:,:) = -1.e+10
270         zphig(:,:) = -1.e+10
271         zmskg(:,:) = -1.e+10
272         ! Add various grids here.
273         DO jj = nldj, nlej
274            DO ji = nldi, nlei
275               zlamg(mig(ji),mjg(jj)) = glamt(ji,jj)
276               zphig(mig(ji),mjg(jj)) = gphit(ji,jj)
277               zmskg(mig(ji),mjg(jj)) = tmask(ji,jj,1)
278            END DO
279         END DO
280         CALL mpp_global_max( zlamg )
281         CALL mpp_global_max( zphig )
282         CALL mpp_global_max( zmskg )
283      ELSE
284         ! Add various grids here.
285         DO jj = 1, jlat
286            DO ji = 1, jlon
287               zlamg(ji,jj) = glamt(ji,jj)
288               zphig(ji,jj) = gphit(ji,jj)
289               zmskg(ji,jj) = tmask(ji,jj,1)
290            END DO
291         END DO
292      ENDIF
293      !-----------------------------------------------------------------------
294      ! Copy longitudes
295      !-----------------------------------------------------------------------
296      ALLOCATE( &
297         & zplam(kobs) &
298         & )
299      DO jo = 1, kobs
300         zplam(jo) = plam(jo)
301      END DO
302      !-----------------------------------------------------------------------
303      ! Set default values for output
304      !-----------------------------------------------------------------------
305      kproc(:) = -1
306      kobsi(:) = -1
307      kobsj(:) = -1
308      !-----------------------------------------------------------------------
309      ! Copy grid positions to temporary arrays and renormalize to 0 to 360.
310      !-----------------------------------------------------------------------
311      DO jj = 1, jlat-1
312         DO ji = 1, jlon-1
313            zlamtm(1,ji,jj) = zlamg(ji  ,jj  )
314            zphitm(1,ji,jj) = zphig(ji  ,jj  )
315            zlamtm(2,ji,jj) = zlamg(ji+1,jj  )
316            zphitm(2,ji,jj) = zphig(ji+1,jj  )
317            zlamtm(3,ji,jj) = zlamg(ji+1,jj+1)
318            zphitm(3,ji,jj) = zphig(ji+1,jj+1)
319            zlamtm(4,ji,jj) = zlamg(ji  ,jj+1)
320            zphitm(4,ji,jj) = zphig(ji  ,jj+1)
321         END DO
322      END DO
323      WHERE ( zlamtm(:,:,:) < 0.0_wp )
324         zlamtm(:,:,:) = zlamtm(:,:,:) + 360.0_wp
325      END WHERE
326      WHERE ( zlamtm(:,:,:) > 360.0_wp )
327         zlamtm(:,:,:) = zlamtm(:,:,:) - 360.0_wp
328      END WHERE
329      !-----------------------------------------------------------------------
330      ! Handle case of the wraparound; beware, not working with orca180
331      !-----------------------------------------------------------------------
332      DO jj = 1, jlat-1
333         DO ji = 1, jlon-1
334            zlammax = MAXVAL( zlamtm(:,ji,jj) )
335            WHERE (zlammax - zlamtm(:, ji, jj) > 180 ) & 
336               & zlamtm(:,ji,jj) = zlamtm(:,ji,jj) + 360._wp
337            zphitmax(ji,jj) = MAXVAL(zphitm(:,ji,jj))
338            zphitmin(ji,jj) = MINVAL(zphitm(:,ji,jj))
339            zlamtmax(ji,jj) = MAXVAL(zlamtm(:,ji,jj))
340            zlamtmin(ji,jj) = MINVAL(zlamtm(:,ji,jj))
341         END DO
342      END DO
343      !-----------------------------------------------------------------------
344      ! Search for boxes with only land points mark them invalid
345      !-----------------------------------------------------------------------
346      llinvalidcell(:,:) = .FALSE.
347      DO jj = 1, jlat-1
348         DO ji = 1, jlon-1
349            llinvalidcell(ji,jj) =               &
350               & zmskg(ji  ,jj  ) == 0.0_wp .AND. &
351               & zmskg(ji+1,jj  ) == 0.0_wp .AND. &
352               & zmskg(ji+1,jj+1) == 0.0_wp .AND. &
353               & zmskg(ji  ,jj+1) == 0.0_wp
354         END DO
355      END DO
356
357      if(lwp) WRITE(numout,*) 'obs_grid_search_lookup do coordinate search using lookup table'
358
359      !-----------------------------------------------------------------------
360      ! Do coordinate search using lookup table with local searches.
361      ! - For land points kproc is set to number of the processor + 1000000
362      !   and we continue the search.
363      ! - For ocean points kproc is set to the number of the processor
364      !   and we stop the search.
365      !-----------------------------------------------------------------------
366      ifourflagcountt = 0
367      ifourflagcountf = 0
368      ifourflagcountr(:) = 0
369
370      !------------------------------------------------------------------------
371      ! Master loop for grid search
372      !------------------------------------------------------------------------
373         
374      gpkobs: DO jo = 1+joffset, kobs, jostride
375         ! Normal case
376         !        specify 4 points which surround the lat lon of interest
377         !                          x i,j+1  x i+1, j+1
378         !
379         !
380         !                             * lon,lat
381         !                          x i,j    x i+1,j
382
383         ! bottom corner point
384         ipx1 = INT( ( zplam(jo)  - lonmin ) / dlon + 1.0 ) 
385         ipy1 = INT( ( pphi (jo)  - latmin ) / dlat + 1.0 ) 
386         
387         ipx = ipx1 + 1
388         ipy = ipy1 + 1
389
390         ! flag for searching around four points separately
391         ! default to false
392         llfourflag = .FALSE.
393         
394         ! check for point fully outside of region
395         IF ( (ipx1 > nlons) .OR. (ipy1 > nlats) .OR. &
396            & (ipx < 1) .OR. (ipy < 1) ) THEN
397            CYCLE
398         ENDIF
399         ! check wrap around
400         IF ( (ipx > nlons) .OR. (ipy > nlats) .OR. &
401            & (ipx1 < 1) .OR. (ipy1 < 1) ) THEN
402            llfourflag=.TRUE.
403            ifourflagcountr(1) = ifourflagcountr(1) + 1
404         ENDIF
405
406         IF (MAXVAL(ixpos(ipx1:ipx,ipy1:ipy)) == -1) CYCLE! cycle if no lookup points found
407         
408         jimin = 0
409         jimax = 0
410         jjmin = 0
411         jjmax = 0
412         
413         IF (.NOT. llfourflag) THEN 
414
415            ! calculate points range
416            ! define a square region encompassing the four corner points
417            ! do I need the -1 points?
418
419            jojimin = MINVAL(ixpos(ipx1:ipx,ipy1:ipy)) - 1
420            jojimax = MAXVAL(ixpos(ipx1:ipx,ipy1:ipy)) + 1
421            jojjmin = MINVAL(iypos(ipx1:ipx,ipy1:ipy)) - 1
422            jojjmax = MAXVAL(iypos(ipx1:ipx,ipy1:ipy)) + 1
423
424            jimin = jojimin - 1
425            jimax = jojimax + 1
426            jjmin = jojjmin - 1
427            jjmax = jojjmax + 1
428           
429            IF ( jojimin < 0 .OR. jojjmin < 0) THEN
430               llfourflag = .TRUE.
431               ifourflagcountr(2) = ifourflagcountr(2) + 1
432            ENDIF
433            IF ( jojimax - jojimin > maxxdiff) THEN
434               llfourflag = .TRUE.
435               ifourflagcountr(3) = ifourflagcountr(3) + 1
436            ENDIF
437            IF ( jojjmax - jojjmin > maxydiff) THEN
438               llfourflag = .TRUE.
439               ifourflagcountr(4) = ifourflagcountr(4) + 1
440            ENDIF
441           
442         ENDIF
443
444         ipmx = 0
445         IF (llfourflag) ipmx = 1
446
447         IF (llfourflag) THEN
448            ifourflagcountt = ifourflagcountt + 1
449         ELSE
450            ifourflagcountf = ifourflagcountf + 1
451         ENDIF
452
453         gridpointsn : DO ip = 0, ipmx
454            DO jp = 0, ipmx
455               
456               IF ( kproc(jo) /= -1 ) EXIT gridpointsn
457       
458               ipx = ipx1 + ip
459               ipy = ipy1 + jp
460               
461               IF (llfourflag) THEN
462
463                  ! deal with wrap around
464                  IF ( ipx > nlons ) ipx = 1
465                  IF ( ipy > nlats ) ipy = 1
466                  IF ( ipx < 1     ) ipx = nlons
467                  IF ( ipy < 1     ) ipy = nlats
468
469                  ! get i,j
470                  isx = ixpos(ipx,ipy)
471                  isy = iypos(ipx,ipy)
472                 
473                  ! estimate appropriate search region (use max/min values)
474                  jimin = isx - maxxdiff - 1
475                  jimax = isx + maxxdiff + 1
476                  jjmin = isy - maxydiff - 1
477                  jjmax = isy + maxydiff + 1
478
479               ENDIF
480
481               IF ( jimin < 1      ) jimin = 1
482               IF ( jimax > jlon-1 ) jimax = jlon-1
483               IF ( jjmin < 1      ) jjmin = 1
484               IF ( jjmax > jlat-1 ) jjmax = jlat-1
485               
486               !---------------------------------------------------------------
487               ! Ensure that all observation longtiudes are between 0 and 360
488               !---------------------------------------------------------------
489
490               IF ( zplam(jo) <   0.0_wp ) zplam(jo) = zplam(jo) + 360.0_wp
491               IF ( zplam(jo) > 360.0_wp ) zplam(jo) = zplam(jo) - 360.0_wp
492     
493               !---------------------------------------------------------------
494               ! Find observations which are on within 1e-6 of a grid point
495               !---------------------------------------------------------------
496
497               gridloop: DO jj = jjmin, jjmax
498                  DO ji = jimin, jimax
499                     IF ( ABS( zphig(ji,jj) - pphi(jo) ) < 1e-6 )  THEN
500                        zlam = zlamg(ji,jj)
501                        IF ( zlam <   0.0_wp ) zlam = zlam + 360.0_wp
502                        IF ( zlam > 360.0_wp ) zlam = zlam - 360.0_wp
503                        IF ( ABS( zlam - zplam(jo) ) < 1e-6 ) THEN
504                           IF ( llinvalidcell(ji,jj) ) THEN
505                              kproc(jo) = nproc + 1000000
506                              kobsi(jo) = ji + 1
507                              kobsj(jo) = jj + 1
508                              CYCLE
509                           ELSE
510                              kproc(jo) = nproc
511                              kobsi(jo) = ji + 1
512                              kobsj(jo) = jj + 1
513                              EXIT gridloop
514                           ENDIF
515                        ENDIF
516                     ENDIF
517                  END DO
518               END DO gridloop
519
520               !---------------------------------------------------------------
521               ! Ensure that all observation longtiudes are between -180/180
522               !---------------------------------------------------------------
523
524               IF ( zplam(jo) > 180 ) zplam(jo) = zplam(jo) - 360.0_wp
525
526               IF ( kproc(jo) == -1 ) THEN
527                 
528                  ! Normal case
529                  gridpoints : DO jj = jjmin, jjmax
530                     DO ji = jimin, jimax
531
532
533                        IF ( ( zplam(jo) > zlamtmax(ji,jj) ) .OR. &
534                           & ( zplam(jo) < zlamtmin(ji,jj) ) ) CYCLE
535                       
536                        IF ( ABS( pphi(jo) ) < 85 ) THEN
537                           IF ( ( pphi(jo) > zphitmax(ji,jj) ) .OR. &
538                              & ( pphi(jo) < zphitmin(ji,jj) ) ) CYCLE
539                        ENDIF
540                       
541                        IF ( linquad( zplam(jo), pphi(jo), &
542                           &          zlamtm(:,ji,jj), zphitm(:,ji,jj) ) ) THEN
543                           IF ( llinvalidcell(ji,jj) ) THEN
544                              kproc(jo) = nproc + 1000000
545                              kobsi(jo) = ji + 1
546                              kobsj(jo) = jj + 1
547                              CYCLE
548                           ELSE
549                              kproc(jo) = nproc
550                              kobsi(jo) = ji + 1
551                              kobsj(jo) = jj + 1
552                              EXIT gridpoints
553                           ENDIF
554                        ENDIF
555                       
556                     END DO
557                  END DO gridpoints
558               ENDIF
559
560               ! In case of failure retry for obs. longtiude + 360.
561               IF ( kproc(jo) == -1 ) THEN
562                  gridpoints_greenwich : DO jj = jjmin, jjmax
563                     DO ji = jimin, jimax
564                       
565                        IF ( ( zplam(jo)+360.0_wp > zlamtmax(ji,jj) ) .OR. &
566                           & ( zplam(jo)+360.0_wp < zlamtmin(ji,jj) ) ) CYCLE
567
568                        IF ( ABS( pphi(jo) ) < 85 ) THEN
569                           IF ( ( pphi(jo) > zphitmax(ji,jj) ) .OR. &
570                              & ( pphi(jo) < zphitmin(ji,jj) ) ) CYCLE
571                        ENDIF
572
573                        IF ( linquad( zplam(jo)+360.0_wp, pphi(jo), &
574                           &          zlamtm(:,ji,jj), zphitm(:,ji,jj) ) ) THEN
575                           IF ( llinvalidcell(ji,jj) ) THEN
576                              kproc(jo) = nproc + 1000000
577                              kobsi(jo) = ji + 1
578                              kobsj(jo) = jj + 1
579                              CYCLE
580                           ELSE
581                              kproc(jo) = nproc
582                              kobsi(jo) = ji + 1
583                              kobsj(jo) = jj + 1
584                              EXIT gridpoints_greenwich
585                           ENDIF
586                        ENDIF
587                       
588                     END DO
589                  END DO gridpoints_greenwich
590                 
591               ENDIF   ! kproc
592               
593            END DO
594         END DO gridpointsn
595      END DO gpkobs  ! kobs
596
597      !----------------------------------------------------------------------
598      ! Synchronize kproc on all processors
599      !----------------------------------------------------------------------
600      IF ( ln_grid_global ) THEN
601         CALL obs_mpp_max_integer( kproc, kobs )
602         CALL obs_mpp_max_integer( kobsi, kobs )
603         CALL obs_mpp_max_integer( kobsj, kobs )
604      ELSE
605         CALL obs_mpp_find_obs_proc( kproc, kobsi, kobsj, kobs )
606      ENDIF
607
608      WHERE( kproc(:) >= 1000000 )
609         kproc(:) = kproc(:) - 1000000
610      END WHERE
611
612      DEALLOCATE( &
613         & zlamg,         &
614         & zphig,         &
615         & zmskg,         &
616         & zphitmax,      &
617         & zphitmin,      &
618         & zlamtmax,      &
619         & zlamtmin,      &
620         & llinvalidcell, &
621         & zlamtm,        &
622         & zphitm,        &
623         & zplam          &
624         & )
625     
626   END SUBROUTINE obs_grid_search_lookup
627
628
629   SUBROUTINE obs_grid_setup
630      !!----------------------------------------------------------------------
631      !!                ***  ROUTINE obs_grid_setup ***
632      !!
633      !! ** Purpose : Setup a lookup table to reduce the searching required
634      !!              for converting lat lons to grid point location
635      !!              produces or reads in a preexisting file for use in
636      !!              obs_grid_search_lookup_local
637      !!
638      !! ** Method : calls obs_grid_search_bruteforce_local with a array
639      !!             of lats and lons
640      !!
641      !! History :
642      !!        !  2007-12 (D. Lea) new routine
643      !!----------------------------------------------------------------------
644     
645      !! * Local declarations
646      CHARACTER(LEN=15), PARAMETER :: &
647         & cpname = 'obs_grid_setup'
648      CHARACTER(LEN=40) :: cfname     
649      INTEGER :: ji
650      INTEGER :: jj
651      INTEGER :: jk
652      INTEGER :: jo
653      INTEGER :: idfile, idny, idnx, idxpos, idypos
654      INTEGER :: idlat, idlon, fileexist
655      INTEGER, DIMENSION(2) :: incdim
656      CHARACTER(LEN=20) :: datestr=" ",timestr=" "
657      REAL(wp) :: tmpx1, tmpx2, tmpy1, tmpy2
658      REAL(wp) :: meanxdiff, meanydiff
659      REAL(wp) :: meanxdiff1, meanydiff1
660      REAL(wp) :: meanxdiff2, meanydiff2
661      INTEGER :: numx1, numx2, numy1, numy2, df
662      INTEGER :: jimin, jimax, jjmin, jjmax
663      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
664         & lonsi,     &
665         & latsi
666      INTEGER, DIMENSION(:,:), ALLOCATABLE :: & 
667         & ixposi,    &
668         & iyposi,    & 
669         & iproci   
670      INTEGER, PARAMETER :: histsize=90
671      INTEGER, DIMENSION(histsize) :: &
672         & histx1, histx2, histy1, histy2
673      REAL, DIMENSION(histsize) :: &
674         & fhistx1, fhistx2, fhisty1, fhisty2
675      REAL(wp) :: histtol
676     
677      IF (ln_grid_search_lookup) THEN
678         
679         WRITE(numout,*) 'Calling obs_grid_setup'
680         
681#ifdef key_datetime_out
682         IF (lwp) THEN
683            CALL DATE_AND_TIME(datestr,timestr)
684            WRITE(numout,*) 'obs_grid_setup begin date_and_time ',datestr,' ',timestr
685         ENDIF
686#endif
687         
688         IF(lwp) WRITE(numout,*)
689         IF(lwp) WRITE(numout,*)'Grid search resolution : ', grid_search_res
690         
691         gsearch_nlons_def  = NINT( 360.0_wp / grid_search_res ) 
692         gsearch_nlats_def  = NINT( 180.0_wp / grid_search_res )
693         gsearch_lonmin_def = -180.0_wp + 0.5_wp * grid_search_res
694         gsearch_latmin_def =  -90.0_wp + 0.5_wp * grid_search_res
695         gsearch_dlon_def   = grid_search_res
696         gsearch_dlat_def   = grid_search_res
697         
698         IF (lwp) THEN
699            WRITE(numout,*)'Grid search gsearch_nlons_def  = ',gsearch_nlons_def
700            WRITE(numout,*)'Grid search gsearch_nlats_def  = ',gsearch_nlats_def
701            WRITE(numout,*)'Grid search gsearch_lonmin_def = ',gsearch_lonmin_def
702            WRITE(numout,*)'Grid search gsearch_latmin_def = ',gsearch_latmin_def
703            WRITE(numout,*)'Grid search gsearch_dlon_def   = ',gsearch_dlon_def
704            WRITE(numout,*)'Grid search gsearch_dlat_def   = ',gsearch_dlat_def
705         ENDIF
706
707         IF ( ln_grid_global ) THEN
708            WRITE(cfname, FMT="(A,'_',A)") &
709               &          TRIM(grid_search_file), 'global.nc'
710         ELSE
711            WRITE(cfname, FMT="(A,'_',I4.4,'of',I4.4,'by',I4.4,'.nc')") &
712               &          TRIM(grid_search_file), nproc, jpni, jpnj
713         ENDIF
714
715         fileexist=nf90_open( TRIM( cfname ), nf90_nowrite, &
716            &                  idfile )
717         
718         IF ( fileexist == nf90_noerr ) THEN
719           
720            ! read data
721            ! initially assume size is as defined (to be fixed)
722           
723            WRITE(numout,*) 'Reading: ',cfname
724           
725            CALL chkerr( nf90_open( TRIM( cfname ), nf90_nowrite, idfile ), &
726               &         cpname, __LINE__ )
727            CALL chkerr( nf90_get_att( idfile, nf90_global, 'maxxdiff', maxxdiff ), &
728               &         cpname, __LINE__ )       
729            CALL chkerr( nf90_get_att( idfile, nf90_global, 'maxydiff', maxydiff ), &
730               &         cpname, __LINE__ )       
731            CALL chkerr( nf90_get_att( idfile, nf90_global, 'dlon', dlon ), &
732            &         cpname, __LINE__ )       
733            CALL chkerr( nf90_get_att( idfile, nf90_global, 'dlat', dlat ), &
734               &         cpname, __LINE__ )       
735            CALL chkerr( nf90_get_att( idfile, nf90_global, 'lonmin', lonmin ), &
736               &         cpname, __LINE__ )       
737            CALL chkerr( nf90_get_att( idfile, nf90_global, 'latmin', latmin ), &
738               &         cpname, __LINE__ )       
739           
740            CALL chkerr( nf90_inq_dimid(idfile, 'nx'  , idnx), &
741               &         cpname, __LINE__ )
742            CALL chkerr( nf90_inquire_dimension( idfile, idnx, len = nlons ),     &
743               &         cpname, __LINE__ ) 
744            CALL chkerr( nf90_inq_dimid(idfile, 'ny'  , idny), &
745               &         cpname, __LINE__ )
746            CALL chkerr( nf90_inquire_dimension( idfile, idny, len = nlats ),     &
747               &         cpname, __LINE__ ) 
748           
749            ALLOCATE( &
750               & lons(nlons,nlats),  &
751               & lats(nlons,nlats),  &
752               & ixpos(nlons,nlats), &
753               & iypos(nlons,nlats), &
754               & iproc(nlons,nlats)  &
755               & )
756           
757            CALL chkerr( nf90_inq_varid( idfile, 'XPOS', idxpos ), & 
758               &         cpname, __LINE__ )
759            CALL chkerr( nf90_get_var  ( idfile, idxpos, ixpos),   &
760               &         cpname, __LINE__ )
761            CALL chkerr( nf90_inq_varid( idfile, 'YPOS', idypos ), & 
762               &         cpname, __LINE__ )
763            CALL chkerr( nf90_get_var  ( idfile, idypos, iypos),   &
764               &         cpname, __LINE__ )
765           
766            CALL chkerr( nf90_close( idfile ), cpname, __LINE__ )
767           
768            ! setup arrays
769           
770            DO ji = 1, nlons
771               DO jj = 1, nlats
772                  lons(ji,jj) = lonmin + (ji-1) * dlon
773                  lats(ji,jj) = latmin + (jj-1) * dlat
774               END DO
775            END DO
776           
777            ! if we are not reading the file we need to create it
778            ! create new obs grid search lookup file
779           
780         ELSE 
781           
782            ! call obs_grid_search
783           
784            IF (lwp) THEN
785               WRITE(numout,*) 'creating: ',cfname
786               WRITE(numout,*) 'calling obs_grid_search: ',nlons*nlats
787            ENDIF
788
789            ! set parameters from default values
790            nlons  = gsearch_nlons_def
791            nlats  = gsearch_nlats_def
792            lonmin = gsearch_lonmin_def
793            latmin = gsearch_latmin_def
794            dlon   = gsearch_dlon_def
795            dlat   = gsearch_dlat_def
796           
797            ! setup arrays
798           
799            ALLOCATE( &
800               & lonsi(nlons,nlats),   &
801               & latsi(nlons,nlats),   &
802               & ixposi(nlons,nlats),  &
803               & iyposi(nlons,nlats),  &
804               & iproci(nlons,nlats)   &
805               & )
806         
807            DO ji = 1, nlons
808               DO jj = 1, nlats
809                  lonsi(ji,jj) = lonmin + (ji-1) * dlon
810                  latsi(ji,jj) = latmin + (jj-1) * dlat
811               END DO
812            END DO
813           
814            CALL obs_grid_search_bruteforce( jpi, jpj, jpiglo, jpjglo,  &
815               &                             nldi, nlei,nldj,  nlej,    &
816               &                             nproc, jpnij,              &
817               &                             glamt, gphit, tmask,       &
818               &                             nlons*nlats, lonsi, latsi, &
819               &                             ixposi, iyposi, iproci )
820           
821            ! minimise file size by removing regions with no data from xypos file
822            ! should be able to just use xpos (ypos will have the same areas of missing data)
823         
824            jimin=1
825            jimax=nlons
826            jjmin=1
827            jjmax=nlats
828
829            minlon_xpos: DO ji= 1, nlons
830               IF (COUNT(ixposi(ji,:) >= 0) > 0) THEN
831                  jimin=ji
832                  EXIT minlon_xpos
833               ENDIF
834            END DO minlon_xpos
835
836            maxlon_xpos: DO ji= nlons, 1, -1
837               IF (COUNT(ixposi(ji,:) >= 0) > 0) THEN
838                  jimax=ji
839                  EXIT maxlon_xpos
840               ENDIF
841            END DO maxlon_xpos
842
843            minlat_xpos: DO jj= 1, nlats
844               IF (COUNT(ixposi(:,jj) >= 0) > 0) THEN
845                  jjmin=jj
846                  EXIT minlat_xpos
847               ENDIF
848            END DO minlat_xpos
849
850            maxlat_xpos: DO jj= nlats, 1, -1
851               IF (COUNT(ixposi(:,jj) >= 0) > 0) THEN
852                  jjmax=jj
853                  EXIT maxlat_xpos
854               ENDIF
855            END DO maxlat_xpos
856
857            lonmin = lonsi(jimin,jjmin)
858            latmin = latsi(jimin,jjmin)
859            nlons  = jimax-jimin+1
860            nlats  = jjmax-jjmin+1
861
862            ! construct new arrays
863
864            ALLOCATE( &
865               & lons(nlons,nlats),    &
866               & lats(nlons,nlats),    &
867               & ixpos(nlons,nlats),   &
868               & iypos(nlons,nlats),   &
869               & iproc(nlons,nlats)    &
870               & )
871
872            lons(:,:) = lonsi(jimin:jimax,jjmin:jjmax)
873            lats(:,:) = latsi(jimin:jimax,jjmin:jjmax)
874            ixpos(:,:) = ixposi(jimin:jimax,jjmin:jjmax)
875            iypos(:,:) = iyposi(jimin:jimax,jjmin:jjmax)
876            iproc(:,:) = iproci(jimin:jimax,jjmin:jjmax)
877
878            DEALLOCATE(lonsi,latsi,ixposi,iyposi,iproci)
879
880            ! calculate (estimate) maxxdiff, maxydiff
881            ! this is to help define the search area for obs_grid_search_lookup
882
883            maxxdiff = 1
884            maxydiff = 1
885
886            tmpx1 = 0
887            tmpx2 = 0
888            tmpy1 = 0
889            tmpy2 = 0 
890
891            numx1 = 0
892            numx2 = 0
893            numy1 = 0
894            numy2 = 0
895
896            ! calculate the mean absolute xdiff and ydiff
897            ! also calculate a histogram
898            ! note the reason why looking for xdiff and ydiff in both directions
899            ! is to allow for rotated grids
900
901            DO ji = 1, nlons-1
902               DO jj = 1, nlats-1
903                  IF ( ixpos(ji,jj) > 0 .AND. iypos(ji,jj) > 0 ) THEN
904                     IF ( ixpos(ji+1,jj) > 0 ) THEN
905                        df = ABS( ixpos(ji+1,jj) - ixpos(ji,jj) )
906                        tmpx1 = tmpx1+df
907                        numx1 = numx1+1
908                        IF ( df < histsize ) histx1(df+1) = histx1(df+1) + 1
909                     ENDIF
910                     IF ( ixpos(ji,jj+1) > 0 ) THEN
911                        df = ABS( ixpos(ji,jj+1) - ixpos(ji,jj) )
912                        tmpx2 = tmpx2 + df
913                        numx2 = numx2 + 1
914                        IF ( df < histsize ) histx2(df+1) = histx2(df+1) + 1
915                     ENDIF
916                     IF (iypos(ji+1,jj) > 0) THEN
917                        df = ABS( iypos(ji+1,jj) - iypos(ji,jj) )
918                        tmpy1 = tmpy1 + df
919                        numy1 = numy1 + 1
920                        IF ( df < histsize ) histy1(df+1) = histy1(df+1) + 1
921                     ENDIF
922                     IF ( iypos(ji,jj+1) > 0 ) THEN
923                        df = ABS( iypos(ji,jj+1) - iypos(ji,jj) )
924                        tmpy2 = tmpy2 + df
925                        numy2 = numy2 + 1
926                        IF ( df < histsize ) histy2(df+1) = histy2(df+1) + 1
927                     ENDIF
928                  ENDIF
929               END DO
930            END DO
931
932            IF (lwp) THEN
933               WRITE(numout,*) 'histograms'
934               WRITE(numout,*) '0   1   2   3   4   5   6   7   8   9   10 ...'
935               WRITE(numout,*) 'histx1'
936               WRITE(numout,*) histx1
937               WRITE(numout,*) 'histx2'
938               WRITE(numout,*) histx2
939               WRITE(numout,*) 'histy1'
940               WRITE(numout,*) histy1
941               WRITE(numout,*) 'histy2'
942               WRITE(numout,*) histy2
943            ENDIF
944
945            meanxdiff1 = tmpx1 / numx1
946            meanydiff1 = tmpy1 / numy1
947            meanxdiff2 = tmpx2 / numx2
948            meanydiff2 = tmpy2 / numy2
949
950            meanxdiff = MAXVAL((/ meanxdiff1, meanxdiff2 /))
951            meanydiff = MAXVAL((/ meanydiff1, meanydiff2 /))
952
953            IF (lwp) THEN
954               WRITE(numout,*) tmpx1, tmpx2, tmpy1, tmpy2
955               WRITE(numout,*) numx1, numx2, numy1, numy2
956               WRITE(numout,*) 'meanxdiff: ',meanxdiff, meanxdiff1, meanxdiff2
957               WRITE(numout,*) 'meanydiff: ',meanydiff, meanydiff1, meanydiff2
958            ENDIF
959
960            tmpx1 = 0
961            tmpx2 = 0
962            tmpy1 = 0
963            tmpy2 = 0
964
965            numx1 = 0
966            numx2 = 0
967            numy1 = 0
968            numy2 = 0
969
970            histx1(:) = 0
971            histx2(:) = 0
972            histy1(:) = 0
973            histy2(:) = 0
974
975            limxdiff = meanxdiff * 4! limit the difference to avoid picking up wraparound
976            limydiff = meanydiff * 4
977
978            DO ji = 1, nlons-1
979               DO jj = 1, nlats-1
980                  IF ( ixpos(ji,jj) > 0 .AND. iypos(ji,jj) > 0 ) THEN
981
982                     IF ( ixpos(ji+1,jj) > 0 ) THEN
983                        df = ABS( ixpos(ji+1,jj)-ixpos(ji,jj) )
984                        tmpx1 = df
985                        IF ( df < limxdiff ) numx1 = numx1+1
986                        IF ( df < histsize ) histx1(df+1) = histx1(df+1) + 1
987                     ENDIF
988                     IF ( ixpos(ji,jj+1) > 0 ) THEN
989                        df = ABS( ixpos(ji,jj+1) - ixpos(ji,jj) )
990                        tmpx2 = df
991                        IF ( df < limxdiff ) numx2 = numx2 + 1
992                        IF ( df < histsize ) histx2(df+1) = histx2(df+1) + 1
993                     ENDIF
994                     IF (iypos(ji+1,jj) > 0) THEN
995                        df = ABS( iypos(ji+1,jj) - iypos(ji,jj) )
996                        tmpy1 = df
997                        IF ( df < limxdiff ) numy1 = numy1 + 1
998                        IF ( df < histsize ) histy1(df+1) = histy1(df+1) + 1
999                     ENDIF
1000                     IF (iypos(ji,jj+1) > 0) THEN
1001                        df = ABS( iypos(ji,jj+1) - iypos(ji,jj) )
1002                        tmpy2 = df
1003                        IF ( df < limxdiff ) numy2 = numy2+1
1004                        IF ( df < histsize ) histy2(df+1) = histy2(df+1)+1
1005                     ENDIF
1006
1007                     IF ( maxxdiff < tmpx1 .AND. tmpx1 < limxdiff ) &
1008                        & maxxdiff = tmpx1
1009                     IF ( maxxdiff < tmpx2 .AND. tmpx2 < limxdiff ) &
1010                        & maxxdiff = tmpx2
1011                     IF ( maxydiff < tmpy1 .AND. tmpy1 < limydiff ) &
1012                        & maxydiff = tmpy1
1013                     IF ( maxydiff < tmpy2 .AND. tmpy2 < limydiff ) &
1014                        & maxydiff = tmpy2
1015
1016                  ENDIF
1017               END DO
1018            END DO
1019
1020            ! cumulative histograms
1021
1022            DO ji = 1, histsize - 1
1023               histx1(ji+1) = histx1(ji+1) + histx1(ji)
1024               histx2(ji+1) = histx2(ji+1) + histx2(ji)
1025               histy1(ji+1) = histy1(ji+1) + histy1(ji)
1026               histy2(ji+1) = histy2(ji+1) + histy2(ji)
1027            END DO
1028
1029            fhistx1(:) = histx1(:) * 1.0 / numx1
1030            fhistx2(:) = histx2(:) * 1.0 / numx2
1031            fhisty1(:) = histy1(:) * 1.0 / numy1
1032            fhisty2(:) = histy2(:) * 1.0 / numy2
1033
1034            ! output new histograms
1035
1036            IF (lwp) THEN
1037               WRITE(numout,*) 'cumulative histograms'
1038               WRITE(numout,*) '0   1   2   3   4   5   6   7   8   9   10 ...'
1039               WRITE(numout,*) 'fhistx1'
1040               WRITE(numout,*) fhistx1
1041               WRITE(numout,*) 'fhistx2'
1042               WRITE(numout,*) fhistx2
1043               WRITE(numout,*) 'fhisty1'
1044               WRITE(numout,*) fhisty1
1045               WRITE(numout,*) 'fhisty2'
1046               WRITE(numout,*) fhisty2
1047            ENDIF
1048
1049            ! calculate maxxdiff and maxydiff based on cumulative histograms
1050            ! where > 0.999 of points are
1051
1052            ! maxval just converts 1x1 vector return from maxloc to a scalar
1053
1054            histtol = 0.999
1055            tmpx1 = MAXVAL( MAXLOC( fhistx1(:), mask = ( fhistx1(:) <= histtol ) ) )
1056            tmpx2 = MAXVAL( MAXLOC( fhistx2(:), mask = ( fhistx2(:) <= histtol ) ) )
1057            tmpy1 = MAXVAL( MAXLOC( fhisty1(:), mask = ( fhisty1(:) <= histtol ) ) )
1058            tmpy2 = MAXVAL( MAXLOC( fhisty2(:), mask = ( fhisty2(:) <= histtol ) ) )
1059
1060            maxxdiff = MAXVAL( (/ tmpx1, tmpx2 /) ) + 1
1061            maxydiff = MAXVAL( (/ tmpy1, tmpy2 /) ) + 1
1062
1063            ! Write out data
1064
1065            IF ( ( .NOT. ln_grid_global ) .OR. &
1066               & ( ( ln_grid_global ) .AND. ( nproc==0 ) ) ) THEN
1067
1068               CALL chkerr( nf90_create (TRIM(cfname), nf90_clobber, idfile), &
1069                  &         cpname, __LINE__ )
1070               CALL chkerr( nf90_put_att( idfile, nf90_global, 'title',       &
1071                  &         'Mapping file from lon/lat to model grid point' ),&
1072                  &         cpname,__LINE__ ) 
1073               CALL chkerr( nf90_put_att( idfile, nf90_global, 'maxxdiff',    &
1074                  &                       maxxdiff ),                         &
1075                  &         cpname,__LINE__ ) 
1076               CALL chkerr( nf90_put_att( idfile, nf90_global, 'maxydiff',    &
1077                  &                       maxydiff ),                         &
1078                  &         cpname,__LINE__ ) 
1079               CALL chkerr( nf90_put_att( idfile, nf90_global, 'dlon', dlon ),&
1080                  &         cpname,__LINE__ ) 
1081               CALL chkerr( nf90_put_att( idfile, nf90_global, 'dlat', dlat ),&
1082                  &         cpname,__LINE__ ) 
1083               CALL chkerr( nf90_put_att( idfile, nf90_global, 'lonmin',      &
1084                  &                       lonmin ),                           &
1085                  &         cpname,__LINE__ ) 
1086               CALL chkerr( nf90_put_att( idfile, nf90_global, 'latmin',      &
1087                  &                       latmin ),                           &
1088                  &         cpname,__LINE__ ) 
1089
1090               CALL chkerr( nf90_def_dim(idfile, 'nx'  , nlons, idnx),        &
1091                  &         cpname,__LINE__ )
1092               CALL chkerr( nf90_def_dim(idfile, 'ny'  , nlats, idny),        &
1093                  &         cpname,__LINE__ )
1094
1095               incdim(1) = idnx
1096               incdim(2) = idny
1097               
1098               CALL chkerr( nf90_def_var( idfile, 'LON', nf90_float, incdim,  &
1099                  &                       idlon ),                            &
1100                  &         cpname, __LINE__ )
1101               CALL chkerr( nf90_put_att( idfile, idlon, 'long_name',         &
1102                  &                       'longitude' ),                      &
1103                  &         cpname, __LINE__ )
1104               
1105               CALL chkerr( nf90_def_var( idfile, 'LAT', nf90_float, incdim,  &
1106                  &                       idlat ),                            &
1107                  &         cpname, __LINE__ )
1108               CALL chkerr( nf90_put_att( idfile, idlat, 'long_name',         &
1109                  &                       'latitude' ),                       &
1110                  &         cpname, __LINE__ )
1111
1112               CALL chkerr( nf90_def_var( idfile, 'XPOS', nf90_int, incdim,   &
1113                  &                       idxpos ),                           &
1114                  &         cpname, __LINE__ )
1115               CALL chkerr( nf90_put_att( idfile, idxpos, 'long_name',        &
1116                  &                       'x position' ),                     &
1117                  &         cpname, __LINE__ )
1118               CALL chkerr( nf90_put_att( idfile, idxpos, '_FillValue', -1 ), &
1119                  &         cpname, __LINE__ )
1120
1121               CALL chkerr( nf90_def_var( idfile, 'YPOS', nf90_int, incdim,   &
1122                  &                       idypos ),                           &
1123                  &         cpname, __LINE__ )
1124               CALL chkerr( nf90_put_att( idfile, idypos, 'long_name',        &
1125                  &                       'y position' ),                     &
1126                  &         cpname, __LINE__ )
1127               CALL chkerr( nf90_put_att( idfile, idypos, '_FillValue', -1 ), &
1128                  &         cpname, __LINE__ )
1129
1130               CALL chkerr( nf90_enddef( idfile ), cpname, __LINE__ )
1131               
1132               CALL chkerr( nf90_put_var( idfile, idlon, lons),               &
1133                  &         cpname, __LINE__ )
1134               CALL chkerr( nf90_put_var( idfile, idlat, lats),               &
1135                  &         cpname, __LINE__ )
1136               CALL chkerr( nf90_put_var( idfile, idxpos, ixpos),             &
1137                  &         cpname, __LINE__ )
1138               CALL chkerr( nf90_put_var( idfile, idypos, iypos),             &
1139                  &         cpname, __LINE__ )
1140               
1141               CALL chkerr( nf90_close( idfile ), cpname, __LINE__ )
1142               
1143               ! should also output max i, max j spacing for use in
1144               ! obs_grid_search_lookup
1145               
1146            ENDIF
1147
1148         ENDIF
1149
1150      ENDIF
1151
1152   END SUBROUTINE obs_grid_setup
1153   
1154   SUBROUTINE obs_grid_deallocate( )
1155      !!----------------------------------------------------------------------
1156      !!                ***  ROUTINE obs_grid_setup ***
1157      !!
1158      !! ** Purpose : Deallocate arrays setup by obs_grid_setup
1159      !!
1160      !! History :
1161      !!        !  2007-12 (D. Lea) new routine
1162      !!-----------------------------------------------------------------------
1163
1164      IF (ln_grid_search_lookup) THEN
1165         DEALLOCATE( lons, lats, ixpos, iypos, iproc )
1166      ENDIF
1167     
1168   END SUBROUTINE obs_grid_deallocate
1169
1170#include "obs_level_search.h90"
1171
1172#include "linquad.h90"
1173
1174#include "maxdist.h90"
1175
1176#include "find_obs_proc.h90"
1177
1178END MODULE obs_grid
1179
Note: See TracBrowser for help on using the repository browser.