Context Navigation

obs_grid.F90 @ 15540

Last change on this file since 15540 was 15540, checked in by sparonuz, 2 years ago
Mixed precision version, tested up to 30 years on ORCA2.
Property svn:keywords set to `Id`
File size: 47.0 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(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(dp), PRIVATE :: gsearch_lonmin_def ! Min longitude
58	REAL(dp), PRIVATE :: gsearch_latmin_def ! Min latitude
59	REAL(dp), PRIVATE :: gsearch_dlon_def ! Lon spacing
60	REAL(dp), 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(dp), PRIVATE :: lonmin ! Min longitude
65	REAL(dp), PRIVATE :: latmin ! Min latitude
66	REAL(dp), PRIVATE :: dlon ! Lon spacing
67	REAL(dp), 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(dp), PRIVATE, DIMENSION(:,:), ALLOCATABLE :: lons, lats
74
75
76
77	INTEGER, PRIVATE, DIMENSION(:,:), ALLOCATABLE :: &
78	& ixpos, &
79	& iypos, &
80	& iprocn
81
82	! Switches
83	LOGICAL, PUBLIC :: ln_grid_search_lookup ! Use lookup table to speed up grid search
84	LOGICAL, PUBLIC :: ln_grid_global ! Use global distribution of observations
85	CHARACTER(LEN=44), PUBLIC :: &
86	& cn_gridsearchfile ! file name head for grid search lookup
87
88
89	# include "single_precision_substitute.h90"
90
91
92	!!----------------------------------------------------------------------
93	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
94	!! $Id$
95	!! Software governed by the CeCILL license (see ./LICENSE)
96	!!----------------------------------------------------------------------
97
98	CONTAINS
99
100	SUBROUTINE obs_grid_search( kobsin, plam, pphi, kobsi, kobsj, kproc, &
101	& cdgrid )
102	!!----------------------------------------------------------------------
103	!! * ROUTINE obs_grid_search *
104	!!
105	!! ** Purpose : Search local gridpoints to find the grid box containing
106	!! the observations calls either
107	!! obs_grd_bruteforce - the original brute force search
108	!! or
109	!! obs_grd_lookup - uses a lookup table to do a fast
110	!!search
111	!!History :
112	!! ! 2007-12 (D. Lea)
113	!!------------------------------------------------------------------------
114
115	!! * Arguments
116	INTEGER :: &
117	& kobsin ! Size of the observation arrays
118	REAL(KIND=dp), DIMENSION(kobsin), INTENT(IN) :: &
119	& plam, & ! Longitude of obsrvations
120	& pphi ! Latitude of observations
121	INTEGER, DIMENSION(kobsin), INTENT(OUT) :: &
122	& kobsi, & ! I-index of observations
123	& kobsj, & ! J-index of observations
124	& kproc ! Processor number of observations
125	CHARACTER(LEN=1) , INTENT(IN) :: &
126	& cdgrid ! Grid to search
127
128	IF(kobsin > 0) THEN
129
130	IF ( ln_grid_search_lookup .AND. ( cdgrid == 'T' ) ) THEN
131	CALL obs_grd_lookup( kobsin, plam, pphi, &
132	& kobsi, kobsj, kproc )
133	ELSE
134	IF ( cdgrid == 'T' ) THEN
135	CALL obs_grd_bruteforce( jpi, jpj, jpiglo, jpjglo, &
136	& 1, jpi, 1, jpj, &
137	& narea-1, jpnij, &
138	& glamt, gphit, tmask, &
139	& kobsin, plam, pphi, &
140	& kobsi, kobsj, kproc )
141	ELSEIF ( cdgrid == 'U' ) THEN
142	CALL obs_grd_bruteforce( jpi, jpj, jpiglo, jpjglo, &
143	& 1, jpi, 1, jpj, &
144	& narea-1, jpnij, &
145	& glamu, gphiu, umask, &
146	& kobsin, plam, pphi, &
147	& kobsi, kobsj, kproc )
148	ELSEIF ( cdgrid == 'V' ) THEN
149	CALL obs_grd_bruteforce( jpi, jpj, jpiglo, jpjglo, &
150	& 1, jpi, 1, jpj, &
151	& narea-1, jpnij, &
152	& glamv, gphiv, vmask, &
153	& kobsin, plam, pphi, &
154	& kobsi, kobsj, kproc )
155	ELSEIF ( cdgrid == 'F' ) THEN
156	CALL obs_grd_bruteforce( jpi, jpj, jpiglo, jpjglo, &
157	& 1, jpi, 1, jpj, &
158	& narea-1, jpnij, &
159	& glamf, gphif, fmask, &
160	& kobsin, plam, pphi, &
161	& kobsi, kobsj, kproc )
162	ELSE
163	CALL ctl_stop( 'Grid not supported' )
164	ENDIF
165	ENDIF
166
167	ENDIF
168
169	END SUBROUTINE obs_grid_search
170
171	#include "obs_grd_bruteforce.h90"
172
173	SUBROUTINE obs_grd_lookup( kobs, plam, pphi, kobsi, kobsj, kproc )
174	!!----------------------------------------------------------------------
175	!! * ROUTINE obs_grid_lookup *
176	!!
177	!! ** Purpose : Search local gridpoints to find the grid box containing
178	!! the observations (much faster then obs_grd_bruteforce)
179	!!
180	!! ** Method : Call to linquad
181	!!
182	!! ** Action : Return kproc holding the observation and kiobsi,kobsj
183	!! valid on kproc=narea-1 processor only.
184	!!
185	!! History :
186	!! ! 2007-12 (D. Lea) new routine based on obs_grid_search
187	!!! updated with fixes from new version of obs_grid_search_bruteforce
188	!!! speeded up where points are not near a "difficult" region like an edge
189	!!----------------------------------------------------------------------
190
191	!! * Arguments
192	INTEGER :: kobs ! Size of the observation arrays
193	REAL(KIND=dp), DIMENSION(kobs), INTENT(IN) :: &
194	& plam, & ! Longitude of obsrvations
195	& pphi ! Latitude of observations
196	INTEGER, DIMENSION(kobs), INTENT(OUT) :: &
197	& kobsi, & ! I-index of observations
198	& kobsj, & ! J-index of observations
199	& kproc ! Processor number of observations
200
201	!! * Local declarations
202	REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: zplam
203	REAL(dp) :: zlammax
204	REAL(dp) :: 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=dp), DIMENSION(:,:), ALLOCATABLE :: &
231	& zlamg, &
232	& zphig, &
233	& zmskg, &
234	& zphitmax,&
235	& zphitmin,&
236	& zlamtmax,&
237	& zlamtmin
238	LOGICAL, DIMENSION(:,:), ALLOCATABLE :: &
239	& llinvalidcell
240	REAL(dp), 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(dp) :: tmpx1, tmpx2, tmpy1, tmpy2
672	REAL(dp) :: meanxdiff, meanydiff
673	REAL(dp) :: meanxdiff1, meanydiff1
674	REAL(dp) :: meanxdiff2, meanydiff2
675	INTEGER :: numx1, numx2, numy1, numy2, df
676	INTEGER :: jimin, jimax, jjmin, jjmax
677	REAL(dp), DIMENSION(:,:), ALLOCATABLE :: lonsi, latsi
678	INTEGER, DIMENSION(:,:), ALLOCATABLE :: &
679	& ixposi, &
680	& iyposi, &
681	& iproci
682	INTEGER, PARAMETER :: histsize=90
683	INTEGER, DIMENSION(histsize) :: &
684	& histx1, histx2, histy1, histy2
685	REAL(dp), DIMENSION(histsize) :: &
686	& fhistx1, fhistx2, fhisty1, fhisty2
687	REAL(dp) :: histtol
688	CHARACTER(LEN=26) :: clfmt ! writing format
689	INTEGER :: idg ! number of digits
690
691	IF (ln_grid_search_lookup) THEN
692
693	WRITE(numout,*) 'Calling obs_grid_setup'
694
695	IF(lwp) WRITE(numout,*)
696	IF(lwp) WRITE(numout,*)'Grid search resolution : ', rn_gridsearchres
697
698	gsearch_nlons_def = NINT( 360.0_wp / rn_gridsearchres )
699	gsearch_nlats_def = NINT( 180.0_wp / rn_gridsearchres )
700	gsearch_lonmin_def = -180.0_wp + 0.5_wp * rn_gridsearchres
701	gsearch_latmin_def = -90.0_wp + 0.5_wp * rn_gridsearchres
702	gsearch_dlon_def = rn_gridsearchres
703	gsearch_dlat_def = rn_gridsearchres
704
705	IF (lwp) THEN
706	WRITE(numout,*)'Grid search gsearch_nlons_def = ',gsearch_nlons_def
707	WRITE(numout,*)'Grid search gsearch_nlats_def = ',gsearch_nlats_def
708	WRITE(numout,*)'Grid search gsearch_lonmin_def = ',gsearch_lonmin_def
709	WRITE(numout,*)'Grid search gsearch_latmin_def = ',gsearch_latmin_def
710	WRITE(numout,*)'Grid search gsearch_dlon_def = ',gsearch_dlon_def
711	WRITE(numout,*)'Grid search gsearch_dlat_def = ',gsearch_dlat_def
712	ENDIF
713
714	IF ( ln_grid_global ) THEN
715	WRITE(cfname, FMT="(A,'_',A)") TRIM(cn_gridsearchfile), 'global.nc'
716	ELSE
717	idg = MAX( INT(LOG10(REAL(jpnij,wp))) + 1, 4 ) ! how many digits to we need to write? min=4, max=9
718	! define the following format: "(a,a,ix.x,a,ix.x,a,ix.x,a)"
719	WRITE(clfmt, "('(a,a,i', i1, '.', i1',a,i', i1, '.', i1',a,i', i1, '.', i1',a)')") idg, idg, idg, idg, idg, idg
720	WRITE(cfname, clfmt ) TRIM(cn_gridsearchfile),'_', narea-1,'of', jpni,'by', jpnj,'.nc'
721	ENDIF
722
723	fileexist=nf90_open( TRIM( cfname ), nf90_nowrite, &
724	& idfile )
725
726	IF ( fileexist == nf90_noerr ) THEN
727
728	! read data
729	! initially assume size is as defined (to be fixed)
730
731	WRITE(numout,*) 'Reading: ',cfname
732
733	CALL chkerr( nf90_open( TRIM( cfname ), nf90_nowrite, idfile ), &
734	& cpname, __LINE__ )
735	CALL chkerr( nf90_get_att( idfile, nf90_global, 'maxxdiff', maxxdiff ), &
736	& cpname, __LINE__ )
737	CALL chkerr( nf90_get_att( idfile, nf90_global, 'maxydiff', maxydiff ), &
738	& cpname, __LINE__ )
739	CALL chkerr( nf90_get_att( idfile, nf90_global, 'dlon', dlon ), &
740	& cpname, __LINE__ )
741	CALL chkerr( nf90_get_att( idfile, nf90_global, 'dlat', dlat ), &
742	& cpname, __LINE__ )
743	CALL chkerr( nf90_get_att( idfile, nf90_global, 'lonmin', lonmin ), &
744	& cpname, __LINE__ )
745	CALL chkerr( nf90_get_att( idfile, nf90_global, 'latmin', latmin ), &
746	& cpname, __LINE__ )
747
748	CALL chkerr( nf90_inq_dimid(idfile, 'nx' , idnx), &
749	& cpname, __LINE__ )
750	CALL chkerr( nf90_inquire_dimension( idfile, idnx, len = nlons ), &
751	& cpname, __LINE__ )
752	CALL chkerr( nf90_inq_dimid(idfile, 'ny' , idny), &
753	& cpname, __LINE__ )
754	CALL chkerr( nf90_inquire_dimension( idfile, idny, len = nlats ), &
755	& cpname, __LINE__ )
756
757	ALLOCATE( &
758	& lons(nlons,nlats), &
759	& lats(nlons,nlats), &
760	& ixpos(nlons,nlats), &
761	& iypos(nlons,nlats), &
762	& iprocn(nlons,nlats) &
763	& )
764
765	CALL chkerr( nf90_inq_varid( idfile, 'XPOS', idxpos ), &
766	& cpname, __LINE__ )
767	CALL chkerr( nf90_get_var ( idfile, idxpos, ixpos), &
768	& cpname, __LINE__ )
769	CALL chkerr( nf90_inq_varid( idfile, 'YPOS', idypos ), &
770	& cpname, __LINE__ )
771	CALL chkerr( nf90_get_var ( idfile, idypos, iypos), &
772	& cpname, __LINE__ )
773
774	CALL chkerr( nf90_close( idfile ), cpname, __LINE__ )
775
776	! setup arrays
777
778	DO ji = 1, nlons
779	DO jj = 1, nlats
780	lons(ji,jj) = lonmin + (ji-1) * dlon
781	lats(ji,jj) = latmin + (jj-1) * dlat
782	END DO
783	END DO
784
785	! if we are not reading the file we need to create it
786	! create new obs grid search lookup file
787
788	ELSE
789
790	! call obs_grid_search
791
792	IF (lwp) THEN
793	WRITE(numout,*) 'creating: ',cfname
794	WRITE(numout,) 'calling obs_grid_search: ',nlonsnlats
795	ENDIF
796
797	! set parameters from default values
798	nlons = gsearch_nlons_def
799	nlats = gsearch_nlats_def
800	lonmin = gsearch_lonmin_def
801	latmin = gsearch_latmin_def
802	dlon = gsearch_dlon_def
803	dlat = gsearch_dlat_def
804
805	! setup arrays
806
807	ALLOCATE( &
808	& lonsi(nlons,nlats), &
809	& latsi(nlons,nlats), &
810	& ixposi(nlons,nlats), &
811	& iyposi(nlons,nlats), &
812	& iproci(nlons,nlats) &
813	& )
814
815	DO ji = 1, nlons
816	DO jj = 1, nlats
817	lonsi(ji,jj) = lonmin + (ji-1) * dlon
818	latsi(ji,jj) = latmin + (jj-1) * dlat
819	END DO
820	END DO
821
822	CALL obs_grd_bruteforce( jpi, jpj, jpiglo, jpjglo, &
823	& 1, jpi, 1, jpj, &
824	& narea-1, jpnij, &
825	& glamt, gphit, tmask, &
826	& nlons*nlats, lonsi, latsi, &
827	& ixposi, iyposi, iproci )
828
829	! minimise file size by removing regions with no data from xypos file
830	! should be able to just use xpos (ypos will have the same areas of missing data)
831
832	jimin=1
833	jimax=nlons
834	jjmin=1
835	jjmax=nlats
836
837	minlon_xpos: DO ji= 1, nlons
838	IF (COUNT(ixposi(ji,:) >= 0) > 0) THEN
839	jimin=ji
840	EXIT minlon_xpos
841	ENDIF
842	END DO minlon_xpos
843
844	maxlon_xpos: DO ji= nlons, 1, -1
845	IF (COUNT(ixposi(ji,:) >= 0) > 0) THEN
846	jimax=ji
847	EXIT maxlon_xpos
848	ENDIF
849	END DO maxlon_xpos
850
851	minlat_xpos: DO jj= 1, nlats
852	IF (COUNT(ixposi(:,jj) >= 0) > 0) THEN
853	jjmin=jj
854	EXIT minlat_xpos
855	ENDIF
856	END DO minlat_xpos
857
858	maxlat_xpos: DO jj= nlats, 1, -1
859	IF (COUNT(ixposi(:,jj) >= 0) > 0) THEN
860	jjmax=jj
861	EXIT maxlat_xpos
862	ENDIF
863	END DO maxlat_xpos
864
865	lonmin = lonsi(jimin,jjmin)
866	latmin = latsi(jimin,jjmin)
867	nlons = jimax-jimin+1
868	nlats = jjmax-jjmin+1
869
870	! construct new arrays
871
872	ALLOCATE( &
873	& lons(nlons,nlats), &
874	& lats(nlons,nlats), &
875	& ixpos(nlons,nlats), &
876	& iypos(nlons,nlats), &
877	& iprocn(nlons,nlats) &
878	& )
879
880	lons(:,:) = lonsi(jimin:jimax,jjmin:jjmax)
881	lats(:,:) = latsi(jimin:jimax,jjmin:jjmax)
882	ixpos(:,:) = ixposi(jimin:jimax,jjmin:jjmax)
883	iypos(:,:) = iyposi(jimin:jimax,jjmin:jjmax)
884	iprocn(:,:) = iproci(jimin:jimax,jjmin:jjmax)
885
886	DEALLOCATE(lonsi,latsi,ixposi,iyposi,iproci)
887
888	! calculate (estimate) maxxdiff, maxydiff
889	! this is to help define the search area for obs_grid_search_lookup
890
891	maxxdiff = 1
892	maxydiff = 1
893
894	tmpx1 = 0
895	tmpx2 = 0
896	tmpy1 = 0
897	tmpy2 = 0
898
899	numx1 = 0
900	numx2 = 0
901	numy1 = 0
902	numy2 = 0
903
904	! calculate the mean absolute xdiff and ydiff
905	! also calculate a histogram
906	! note the reason why looking for xdiff and ydiff in both directions
907	! is to allow for rotated grids
908
909	DO ji = 1, nlons-1
910	DO jj = 1, nlats-1
911	IF ( ixpos(ji,jj) > 0 .AND. iypos(ji,jj) > 0 ) THEN
912	IF ( ixpos(ji+1,jj) > 0 ) THEN
913	df = ABS( ixpos(ji+1,jj) - ixpos(ji,jj) )
914	tmpx1 = tmpx1+df
915	numx1 = numx1+1
916	IF ( df < histsize ) histx1(df+1) = histx1(df+1) + 1
917	ENDIF
918	IF ( ixpos(ji,jj+1) > 0 ) THEN
919	df = ABS( ixpos(ji,jj+1) - ixpos(ji,jj) )
920	tmpx2 = tmpx2 + df
921	numx2 = numx2 + 1
922	IF ( df < histsize ) histx2(df+1) = histx2(df+1) + 1
923	ENDIF
924	IF (iypos(ji+1,jj) > 0) THEN
925	df = ABS( iypos(ji+1,jj) - iypos(ji,jj) )
926	tmpy1 = tmpy1 + df
927	numy1 = numy1 + 1
928	IF ( df < histsize ) histy1(df+1) = histy1(df+1) + 1
929	ENDIF
930	IF ( iypos(ji,jj+1) > 0 ) THEN
931	df = ABS( iypos(ji,jj+1) - iypos(ji,jj) )
932	tmpy2 = tmpy2 + df
933	numy2 = numy2 + 1
934	IF ( df < histsize ) histy2(df+1) = histy2(df+1) + 1
935	ENDIF
936	ENDIF
937	END DO
938	END DO
939
940	IF (lwp) THEN
941	WRITE(numout,*) 'histograms'
942	WRITE(numout,*) '0 1 2 3 4 5 6 7 8 9 10 ...'
943	WRITE(numout,*) 'histx1'
944	WRITE(numout,*) histx1
945	WRITE(numout,*) 'histx2'
946	WRITE(numout,*) histx2
947	WRITE(numout,*) 'histy1'
948	WRITE(numout,*) histy1
949	WRITE(numout,*) 'histy2'
950	WRITE(numout,*) histy2
951	ENDIF
952
953	meanxdiff1 = tmpx1 / numx1
954	meanydiff1 = tmpy1 / numy1
955	meanxdiff2 = tmpx2 / numx2
956	meanydiff2 = tmpy2 / numy2
957
958	meanxdiff = MAXVAL((/ meanxdiff1, meanxdiff2 /))
959	meanydiff = MAXVAL((/ meanydiff1, meanydiff2 /))
960
961	IF (lwp) THEN
962	WRITE(numout,*) tmpx1, tmpx2, tmpy1, tmpy2
963	WRITE(numout,*) numx1, numx2, numy1, numy2
964	WRITE(numout,*) 'meanxdiff: ',meanxdiff, meanxdiff1, meanxdiff2
965	WRITE(numout,*) 'meanydiff: ',meanydiff, meanydiff1, meanydiff2
966	ENDIF
967
968	tmpx1 = 0
969	tmpx2 = 0
970	tmpy1 = 0
971	tmpy2 = 0
972
973	numx1 = 0
974	numx2 = 0
975	numy1 = 0
976	numy2 = 0
977
978	histx1(:) = 0
979	histx2(:) = 0
980	histy1(:) = 0
981	histy2(:) = 0
982
983	limxdiff = meanxdiff * 4! limit the difference to avoid picking up wraparound
984	limydiff = meanydiff * 4
985
986	DO ji = 1, nlons-1
987	DO jj = 1, nlats-1
988	IF ( ixpos(ji,jj) > 0 .AND. iypos(ji,jj) > 0 ) THEN
989
990	IF ( ixpos(ji+1,jj) > 0 ) THEN
991	df = ABS( ixpos(ji+1,jj)-ixpos(ji,jj) )
992	tmpx1 = df
993	IF ( df < limxdiff ) numx1 = numx1+1
994	IF ( df < histsize ) histx1(df+1) = histx1(df+1) + 1
995	ENDIF
996	IF ( ixpos(ji,jj+1) > 0 ) THEN
997	df = ABS( ixpos(ji,jj+1) - ixpos(ji,jj) )
998	tmpx2 = df
999	IF ( df < limxdiff ) numx2 = numx2 + 1
1000	IF ( df < histsize ) histx2(df+1) = histx2(df+1) + 1
1001	ENDIF
1002	IF (iypos(ji+1,jj) > 0) THEN
1003	df = ABS( iypos(ji+1,jj) - iypos(ji,jj) )
1004	tmpy1 = df
1005	IF ( df < limydiff ) numy1 = numy1 + 1
1006	IF ( df < histsize ) histy1(df+1) = histy1(df+1) + 1
1007	ENDIF
1008	IF (iypos(ji,jj+1) > 0) THEN
1009	df = ABS( iypos(ji,jj+1) - iypos(ji,jj) )
1010	tmpy2 = df
1011	IF ( df < limydiff ) numy2 = numy2+1
1012	IF ( df < histsize ) histy2(df+1) = histy2(df+1)+1
1013	ENDIF
1014
1015	IF ( maxxdiff < tmpx1 .AND. tmpx1 < limxdiff ) &
1016	& maxxdiff = tmpx1
1017	IF ( maxxdiff < tmpx2 .AND. tmpx2 < limxdiff ) &
1018	& maxxdiff = tmpx2
1019	IF ( maxydiff < tmpy1 .AND. tmpy1 < limydiff ) &
1020	& maxydiff = tmpy1
1021	IF ( maxydiff < tmpy2 .AND. tmpy2 < limydiff ) &
1022	& maxydiff = tmpy2
1023
1024	ENDIF
1025	END DO
1026	END DO
1027
1028	! cumulative histograms
1029
1030	DO ji = 1, histsize - 1
1031	histx1(ji+1) = histx1(ji+1) + histx1(ji)
1032	histx2(ji+1) = histx2(ji+1) + histx2(ji)
1033	histy1(ji+1) = histy1(ji+1) + histy1(ji)
1034	histy2(ji+1) = histy2(ji+1) + histy2(ji)
1035	END DO
1036
1037	fhistx1(:) = histx1(:) * 1.0 / numx1
1038	fhistx2(:) = histx2(:) * 1.0 / numx2
1039	fhisty1(:) = histy1(:) * 1.0 / numy1
1040	fhisty2(:) = histy2(:) * 1.0 / numy2
1041
1042	! output new histograms
1043
1044	IF (lwp) THEN
1045	WRITE(numout,*) 'cumulative histograms'
1046	WRITE(numout,*) '0 1 2 3 4 5 6 7 8 9 10 ...'
1047	WRITE(numout,*) 'fhistx1'
1048	WRITE(numout,*) fhistx1
1049	WRITE(numout,*) 'fhistx2'
1050	WRITE(numout,*) fhistx2
1051	WRITE(numout,*) 'fhisty1'
1052	WRITE(numout,*) fhisty1
1053	WRITE(numout,*) 'fhisty2'
1054	WRITE(numout,*) fhisty2
1055	ENDIF
1056
1057	! calculate maxxdiff and maxydiff based on cumulative histograms
1058	! where > 0.999 of points are
1059
1060	! maxval just converts 1x1 vector return from maxloc to a scalar
1061
1062	histtol = 0.999
1063	tmpx1 = MAXVAL( MAXLOC( fhistx1(:), mask = ( fhistx1(:) <= histtol ) ) )
1064	tmpx2 = MAXVAL( MAXLOC( fhistx2(:), mask = ( fhistx2(:) <= histtol ) ) )
1065	tmpy1 = MAXVAL( MAXLOC( fhisty1(:), mask = ( fhisty1(:) <= histtol ) ) )
1066	tmpy2 = MAXVAL( MAXLOC( fhisty2(:), mask = ( fhisty2(:) <= histtol ) ) )
1067
1068	maxxdiff = MAXVAL( (/ tmpx1, tmpx2 /) ) + 1
1069	maxydiff = MAXVAL( (/ tmpy1, tmpy2 /) ) + 1
1070
1071	! Write out data
1072
1073	IF ( ( .NOT. ln_grid_global ) .OR. &
1074	& ( ( ln_grid_global ) .AND. ( narea-1==0 ) ) ) THEN
1075
1076	CALL chkerr( nf90_create (TRIM(cfname), nf90_clobber, idfile), &
1077	& cpname, __LINE__ )
1078	CALL chkerr( nf90_put_att( idfile, nf90_global, 'title', &
1079	& 'Mapping file from lon/lat to model grid point' ),&
1080	& cpname,__LINE__ )
1081	CALL chkerr( nf90_put_att( idfile, nf90_global, 'maxxdiff', &
1082	& maxxdiff ), &
1083	& cpname,__LINE__ )
1084	CALL chkerr( nf90_put_att( idfile, nf90_global, 'maxydiff', &
1085	& maxydiff ), &
1086	& cpname,__LINE__ )
1087	CALL chkerr( nf90_put_att( idfile, nf90_global, 'dlon', dlon ),&
1088	& cpname,__LINE__ )
1089	CALL chkerr( nf90_put_att( idfile, nf90_global, 'dlat', dlat ),&
1090	& cpname,__LINE__ )
1091	CALL chkerr( nf90_put_att( idfile, nf90_global, 'lonmin', &
1092	& lonmin ), &
1093	& cpname,__LINE__ )
1094	CALL chkerr( nf90_put_att( idfile, nf90_global, 'latmin', &
1095	& latmin ), &
1096	& cpname,__LINE__ )
1097
1098	CALL chkerr( nf90_def_dim(idfile, 'nx' , nlons, idnx), &
1099	& cpname,__LINE__ )
1100	CALL chkerr( nf90_def_dim(idfile, 'ny' , nlats, idny), &
1101	& cpname,__LINE__ )
1102
1103	incdim(1) = idnx
1104	incdim(2) = idny
1105
1106	CALL chkerr( nf90_def_var( idfile, 'LON', nf90_float, incdim, &
1107	& idlon ), &
1108	& cpname, __LINE__ )
1109	CALL chkerr( nf90_put_att( idfile, idlon, 'long_name', &
1110	& 'longitude' ), &
1111	& cpname, __LINE__ )
1112
1113	CALL chkerr( nf90_def_var( idfile, 'LAT', nf90_float, incdim, &
1114	& idlat ), &
1115	& cpname, __LINE__ )
1116	CALL chkerr( nf90_put_att( idfile, idlat, 'long_name', &
1117	& 'latitude' ), &
1118	& cpname, __LINE__ )
1119
1120	CALL chkerr( nf90_def_var( idfile, 'XPOS', nf90_int, incdim, &
1121	& idxpos ), &
1122	& cpname, __LINE__ )
1123	CALL chkerr( nf90_put_att( idfile, idxpos, 'long_name', &
1124	& 'x position' ), &
1125	& cpname, __LINE__ )
1126	CALL chkerr( nf90_put_att( idfile, idxpos, '_FillValue', -1 ), &
1127	& cpname, __LINE__ )
1128
1129	CALL chkerr( nf90_def_var( idfile, 'YPOS', nf90_int, incdim, &
1130	& idypos ), &
1131	& cpname, __LINE__ )
1132	CALL chkerr( nf90_put_att( idfile, idypos, 'long_name', &
1133	& 'y position' ), &
1134	& cpname, __LINE__ )
1135	CALL chkerr( nf90_put_att( idfile, idypos, '_FillValue', -1 ), &
1136	& cpname, __LINE__ )
1137
1138	CALL chkerr( nf90_enddef( idfile ), cpname, __LINE__ )
1139
1140	CALL chkerr( nf90_put_var( idfile, idlon, lons), &
1141	& cpname, __LINE__ )
1142	CALL chkerr( nf90_put_var( idfile, idlat, lats), &
1143	& cpname, __LINE__ )
1144	CALL chkerr( nf90_put_var( idfile, idxpos, ixpos), &
1145	& cpname, __LINE__ )
1146	CALL chkerr( nf90_put_var( idfile, idypos, iypos), &
1147	& cpname, __LINE__ )
1148
1149	CALL chkerr( nf90_close( idfile ), cpname, __LINE__ )
1150
1151	! should also output max i, max j spacing for use in
1152	! obs_grid_search_lookup
1153
1154	ENDIF
1155
1156	ENDIF
1157
1158	ENDIF
1159
1160	END SUBROUTINE obs_grid_setup
1161
1162	SUBROUTINE obs_grid_deallocate( )
1163	!!----------------------------------------------------------------------
1164	!! * ROUTINE obs_grid_setup *
1165	!!
1166	!! ** Purpose : Deallocate arrays setup by obs_grid_setup
1167	!!
1168	!! History :
1169	!! ! 2007-12 (D. Lea) new routine
1170	!!-----------------------------------------------------------------------
1171
1172	IF (ln_grid_search_lookup) THEN
1173	DEALLOCATE( lons, lats, ixpos, iypos, iprocn )
1174	ENDIF
1175
1176	END SUBROUTINE obs_grid_deallocate
1177
1178	#include "obs_level_search.h90"
1179
1180	#include "linquad.h90"
1181
1182	#include "maxdist.h90"
1183
1184	#include "find_obs_proc.h90"
1185
1186	END MODULE obs_grid
1187

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

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

Download in other formats: