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

Context Navigation

source: branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/OBS/obs_grid.F90 @ 10246

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

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

Download in other formats: