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 NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/OBS – NEMO

source: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/OBS/obs_grid.F90 @ 14894

Last change on this file since 14894 was 14894, checked in by sparonuz, 3 years ago

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