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

Context Navigation

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/OBS/obs_grid.F90 @ 4428

Last change on this file since 4428 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.

Property svn:keywords set to Id

File size: 47.3 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: