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 @ 2287

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

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