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 trunk/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/OBS/obs_grid.F90 @ 2715

Last change on this file since 2715 was 2715, checked in by rblod, 13 years ago

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