Context Navigation

obs_prep.F90 @ 7713

Last change on this file since 7713 was 7713, checked in by dford, 7 years ago
Add observation operator code for surface log10(chlorophyll), SPM, pCO2 and fCO2, for use with FABM-ERSEM, HadOCC and MEDUSA. See internal Met Office NEMO ticket 660.
File size: 123.5 KB

Line
1	MODULE obs_prep
2	!!=====================================================================
3	!! * MODULE obs_prep *
4	!! Observation diagnostics: Prepare observation arrays: screening,
5	!! sorting, coordinate search
6	!!=====================================================================
7
8	!!---------------------------------------------------------------------
9	!! obs_pre_pro : First level check and screening of T/S profiles
10	!! obs_pre_sla : First level check and screening of SLA observations
11	!! obs_pre_sst : First level check and screening of SLA observations
12	!! obs_pre_seaice : First level check and screening of sea ice observations
13	!! obs_pre_vel : First level check and screening of velocity obs.
14	!! obs_pre_logchl : First level check and screening of logchl obs.
15	!! obs_pre_spm : First level check and screening of spm obs.
16	!! obs_pre_fco2 : First level check and screening of fco2 obs.
17	!! obs_pre_pco2 : First level check and screening of pco2 obs.
18	!! obs_scr : Basic screening of the observations
19	!! obs_coo_tim : Compute number of time steps to the observation time
20	!! obs_sor : Sort the observation arrays
21	!!---------------------------------------------------------------------
22	!! * Modules used
23	USE par_kind, ONLY : & ! Precision variables
24	& wp
25	USE in_out_manager ! I/O manager
26	USE obs_profiles_def ! Definitions for storage arrays for profiles
27	USE obs_surf_def ! Definitions for storage arrays for surface data
28	USE obs_mpp, ONLY : & ! MPP support routines for observation diagnostics
29	& obs_mpp_sum_integer, &
30	& obs_mpp_sum_integers
31	USE obs_inter_sup ! Interpolation support
32	USE obs_oper ! Observation operators
33	#if defined key_bdy
34	USE bdy_oce, ONLY : & ! Boundary information
35	idx_bdy, nb_bdy
36	#endif
37	USE lib_mpp, ONLY : &
38	& ctl_warn, ctl_stop
39
40	IMPLICIT NONE
41
42	!! * Routine accessibility
43	PRIVATE
44
45	PUBLIC &
46	& obs_pre_pro, & ! First level check and screening of profiles
47	& obs_pre_sla, & ! First level check and screening of SLA data
48	& obs_pre_sst, & ! First level check and screening of SLA data
49	& obs_pre_seaice, & ! First level check and screening of sea ice data
50	& obs_pre_vel, & ! First level check and screening of velocity profiles
51	& obs_pre_logchl, & ! First level check and screening of logchl data
52	& obs_pre_spm, & ! First level check and screening of spm data
53	& obs_pre_fco2, & ! First level check and screening of fco2 data
54	& obs_pre_pco2, & ! First level check and screening of pco2 data
55	& calc_month_len ! Calculate the number of days in the months of a year
56
57	LOGICAL, PUBLIC :: ln_bound_reject !: Remove obs near open boundaries
58
59	!!----------------------------------------------------------------------
60	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
61	!! $Id$
62	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
63	!!----------------------------------------------------------------------
64
65	!! * Substitutions
66	# include "domzgr_substitute.h90"
67
68	CONTAINS
69
70	SUBROUTINE obs_pre_pro( profdata, prodatqc, ld_t3d, ld_s3d, ld_nea, &
71	& kdailyavtypes )
72	!!----------------------------------------------------------------------
73	!! * ROUTINE obs_pre_pro *
74	!!
75	!! ** Purpose : First level check and screening of T and S profiles
76	!!
77	!! ** Method : First level check and screening of T and S profiles
78	!!
79	!! ** Action :
80	!!
81	!! References :
82	!!
83	!! History :
84	!! ! 2007-01 (K. Mogensen) Merge of obs_pre_t3d and obs_pre_s3d
85	!! ! 2007-03 (K. Mogensen) General handling of profiles
86	!! ! 2007-06 (K. Mogensen et al) Reject obs. near land.
87	!!----------------------------------------------------------------------
88	!! * Modules used
89	USE domstp ! Domain: set the time-step
90	USE par_oce ! Ocean parameters
91	USE dom_oce, ONLY : & ! Geographical information
92	& glamt, &
93	& gphit, &
94	& gdept_1d,&
95	#if defined key_vvl
96	& gdepw_n, &
97	& gdept_n, &
98	#else
99	& gdepw_1d, &
100	& gdept_1d, &
101	#endif
102	& tmask, &
103	& ln_zco, &
104	& ln_zps, &
105	& nproc
106	!! * Arguments
107	TYPE(obs_prof), INTENT(INOUT) :: profdata ! Full set of profile data
108	TYPE(obs_prof), INTENT(INOUT) :: prodatqc ! Subset of profile data not failing screening
109	LOGICAL, INTENT(IN) :: ld_t3d ! Switch for temperature
110	LOGICAL, INTENT(IN) :: ld_s3d ! Switch for salinity
111	LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land
112	INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
113	& kdailyavtypes! Types for daily averages
114	!! * Local declarations
115	INTEGER :: iyea0 ! Initial date
116	INTEGER :: imon0 ! - (year, month, day, hour, minute)
117	INTEGER :: iday0
118	INTEGER :: ihou0
119	INTEGER :: imin0
120	INTEGER :: icycle ! Current assimilation cycle
121	! Counters for observations that
122	INTEGER :: iotdobs ! - outside time domain
123	INTEGER :: iosdtobs ! - outside space domain (temperature)
124	INTEGER :: iosdsobs ! - outside space domain (salinity)
125	INTEGER :: ilantobs ! - within a model land cell (temperature)
126	INTEGER :: ilansobs ! - within a model land cell (salinity)
127	INTEGER :: ibdytobs ! - boundary (temperature)
128	INTEGER :: ibdysobs ! - boundary (salinity)
129	INTEGER :: inlatobs ! - close to land (temperature)
130	INTEGER :: inlasobs ! - close to land (salinity)
131	INTEGER :: igrdobs ! - fail the grid search
132	! Global counters for observations that
133	INTEGER :: iotdobsmpp ! - outside time domain
134	INTEGER :: iosdtobsmpp ! - outside space domain (temperature)
135	INTEGER :: iosdsobsmpp ! - outside space domain (salinity)
136	INTEGER :: ilantobsmpp ! - within a model land cell (temperature)
137	INTEGER :: ilansobsmpp ! - within a model land cell (salinity)
138	INTEGER :: ibdytobsmpp ! - boundary (temperature)
139	INTEGER :: ibdysobsmpp ! - boundary (salinity)
140	INTEGER :: inlatobsmpp ! - close to land (temperature)
141	INTEGER :: inlasobsmpp ! - close to land (salinity)
142	INTEGER :: igrdobsmpp ! - fail the grid search
143	TYPE(obs_prof_valid) :: llvalid ! Profile selection
144	TYPE(obs_prof_valid), DIMENSION(profdata%nvar) :: &
145	& llvvalid ! T,S selection
146	INTEGER :: jvar ! Variable loop variable
147	INTEGER :: jobs ! Obs. loop variable
148	INTEGER :: jstp ! Time loop variable
149	INTEGER :: inrc ! Time index variable
150
151	IF(lwp) WRITE(numout,*)'obs_pre_pro : Preparing the profile observations...'
152
153	! Initial date initialization (year, month, day, hour, minute)
154	iyea0 = ndate0 / 10000
155	imon0 = ( ndate0 - iyea0 * 10000 ) / 100
156	iday0 = ndate0 - iyea0 * 10000 - imon0 * 100
157	ihou0 = 0
158	imin0 = 0
159
160	icycle = no ! Assimilation cycle
161
162	! Diagnotics counters for various failures.
163
164	iotdobs = 0
165	igrdobs = 0
166	iosdtobs = 0
167	iosdsobs = 0
168	ilantobs = 0
169	ilansobs = 0
170	inlatobs = 0
171	inlasobs = 0
172	ibdytobs = 0
173	ibdysobs = 0
174
175	! -----------------------------------------------------------------------
176	! Find time coordinate for profiles
177	! -----------------------------------------------------------------------
178
179	IF ( PRESENT(kdailyavtypes) ) THEN
180	CALL obs_coo_tim_prof( icycle, &
181	& iyea0, imon0, iday0, ihou0, imin0, &
182	& profdata%nprof, profdata%nyea, profdata%nmon, &
183	& profdata%nday, profdata%nhou, profdata%nmin, &
184	& profdata%ntyp, profdata%nqc, profdata%mstp, &
185	& iotdobs, kdailyavtypes = kdailyavtypes )
186	ELSE
187	CALL obs_coo_tim_prof( icycle, &
188	& iyea0, imon0, iday0, ihou0, imin0, &
189	& profdata%nprof, profdata%nyea, profdata%nmon, &
190	& profdata%nday, profdata%nhou, profdata%nmin, &
191	& profdata%ntyp, profdata%nqc, profdata%mstp, &
192	& iotdobs )
193	ENDIF
194	CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
195
196	! -----------------------------------------------------------------------
197	! Check for profiles failing the grid search
198	! -----------------------------------------------------------------------
199
200	CALL obs_coo_grd( profdata%nprof, profdata%mi, profdata%mj, &
201	& profdata%nqc, igrdobs )
202
203	CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
204
205	! -----------------------------------------------------------------------
206	! Reject all observations for profiles with nqc > 10
207	! -----------------------------------------------------------------------
208
209	CALL obs_pro_rej( profdata )
210
211	! -----------------------------------------------------------------------
212	! Check for land points. This includes points below the model
213	! bathymetry so this is done for every point in the profile
214	! -----------------------------------------------------------------------
215
216	! Temperature
217
218	CALL obs_coo_spc_3d( profdata%nprof, profdata%nvprot(1), &
219	& profdata%npvsta(:,1), profdata%npvend(:,1), &
220	& jpi, jpj, &
221	& jpk, &
222	& profdata%mi, profdata%mj, &
223	& profdata%var(1)%mvk, &
224	& profdata%rlam, profdata%rphi, &
225	& profdata%var(1)%vdep, &
226	& glamt, gphit, &
227	& gdept_1d, tmask, &
228	& profdata%nqc, profdata%var(1)%nvqc, &
229	& iosdtobs, ilantobs, &
230	& inlatobs, ld_nea, &
231	& ibdytobs, ln_bound_reject )
232
233
234	CALL obs_mpp_sum_integer( iosdtobs, iosdtobsmpp )
235	CALL obs_mpp_sum_integer( ilantobs, ilantobsmpp )
236	CALL obs_mpp_sum_integer( inlatobs, inlatobsmpp )
237	CALL obs_mpp_sum_integer( ibdytobs, ibdytobsmpp )
238
239	! Salinity
240
241	CALL obs_coo_spc_3d( profdata%nprof, profdata%nvprot(2), &
242	& profdata%npvsta(:,2), profdata%npvend(:,2), &
243	& jpi, jpj, &
244	& jpk, &
245	& profdata%mi, profdata%mj, &
246	& profdata%var(2)%mvk, &
247	& profdata%rlam, profdata%rphi, &
248	& profdata%var(2)%vdep, &
249	& glamt, gphit, &
250	& gdept_1d, tmask, &
251	& profdata%nqc, profdata%var(2)%nvqc, &
252	& iosdsobs, ilansobs, &
253	& inlasobs, ld_nea, &
254	& ibdysobs, ln_bound_reject )
255
256	CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
257	CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
258	CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
259	CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp )
260
261	! -----------------------------------------------------------------------
262	! Copy useful data from the profdata data structure to
263	! the prodatqc data structure
264	! -----------------------------------------------------------------------
265
266	! Allocate the selection arrays
267
268	ALLOCATE( llvalid%luse(profdata%nprof) )
269	DO jvar = 1,profdata%nvar
270	ALLOCATE( llvvalid(jvar)%luse(profdata%nvprot(jvar)) )
271	END DO
272
273	! We want all data which has qc flags <= 10
274
275	llvalid%luse(:) = ( profdata%nqc(:) <= 10 )
276	DO jvar = 1,profdata%nvar
277	llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= 10 )
278	END DO
279
280	! The actual copying
281
282	CALL obs_prof_compress( profdata, prodatqc, .TRUE., numout, &
283	& lvalid=llvalid, lvvalid=llvvalid )
284
285	! Dellocate the selection arrays
286	DEALLOCATE( llvalid%luse )
287	DO jvar = 1,profdata%nvar
288	DEALLOCATE( llvvalid(jvar)%luse )
289	END DO
290
291	! -----------------------------------------------------------------------
292	! Print information about what observations are left after qc
293	! -----------------------------------------------------------------------
294
295	! Update the total observation counter array
296
297	IF(lwp) THEN
298	WRITE(numout,*)
299	WRITE(numout,*) 'obs_pre_pro :'
300	WRITE(numout,*) '~~~~~~~~~~~'
301	WRITE(numout,*)
302	WRITE(numout,*) ' Profiles outside time domain = ', &
303	& iotdobsmpp
304	WRITE(numout,*) ' Remaining profiles that failed grid search = ', &
305	& igrdobsmpp
306	WRITE(numout,*) ' Remaining T data outside space domain = ', &
307	& iosdtobsmpp
308	WRITE(numout,*) ' Remaining T data at land points = ', &
309	& ilantobsmpp
310	IF (ld_nea) THEN
311	WRITE(numout,*) ' Remaining T data near land points (removed) = ',&
312	& inlatobsmpp
313	ELSE
314	WRITE(numout,*) ' Remaining T data near land points (kept) = ',&
315	& inlatobsmpp
316	ENDIF
317	WRITE(numout,*) ' Remaining T data near open boundary (removed) = ',&
318	& ibdytobsmpp
319	WRITE(numout,*) ' T data accepted = ', &
320	& prodatqc%nvprotmpp(1)
321	WRITE(numout,*) ' Remaining S data outside space domain = ', &
322	& iosdsobsmpp
323	WRITE(numout,*) ' Remaining S data at land points = ', &
324	& ilansobsmpp
325	IF (ld_nea) THEN
326	WRITE(numout,*) ' Remaining S data near land points (removed) = ',&
327	& inlasobsmpp
328	ELSE
329	WRITE(numout,*) ' Remaining S data near land points (kept) = ',&
330	& inlasobsmpp
331	ENDIF
332	WRITE(numout,*) ' Remaining S data near open boundary (removed) = ',&
333	& ibdysobsmpp
334	WRITE(numout,*) ' S data accepted = ', &
335	& prodatqc%nvprotmpp(2)
336
337	WRITE(numout,*)
338	WRITE(numout,*) ' Number of observations per time step :'
339	WRITE(numout,*)
340	WRITE(numout,997)
341	WRITE(numout,998)
342	ENDIF
343
344	DO jobs = 1, prodatqc%nprof
345	inrc = prodatqc%mstp(jobs) + 2 - nit000
346	prodatqc%npstp(inrc) = prodatqc%npstp(inrc) + 1
347	DO jvar = 1, prodatqc%nvar
348	IF ( prodatqc%npvend(jobs,jvar) > 0 ) THEN
349	prodatqc%nvstp(inrc,jvar) = prodatqc%nvstp(inrc,jvar) + &
350	& ( prodatqc%npvend(jobs,jvar) - &
351	& prodatqc%npvsta(jobs,jvar) + 1 )
352	ENDIF
353	END DO
354	END DO
355
356
357	CALL obs_mpp_sum_integers( prodatqc%npstp, prodatqc%npstpmpp, &
358	& nitend - nit000 + 2 )
359	DO jvar = 1, prodatqc%nvar
360	CALL obs_mpp_sum_integers( prodatqc%nvstp(:,jvar), &
361	& prodatqc%nvstpmpp(:,jvar), &
362	& nitend - nit000 + 2 )
363	END DO
364
365	IF ( lwp ) THEN
366	DO jstp = nit000 - 1, nitend
367	inrc = jstp - nit000 + 2
368	WRITE(numout,999) jstp, prodatqc%npstpmpp(inrc), &
369	& prodatqc%nvstpmpp(inrc,1), &
370	& prodatqc%nvstpmpp(inrc,2)
371	END DO
372	ENDIF
373
374	997 FORMAT(10X,'Time step',5X,'Profiles',5X,'Temperature',5X,'Salinity')
375	998 FORMAT(10X,'---------',5X,'--------',5X,'-----------',5X,'--------')
376	999 FORMAT(10X,I9,5X,I8,5X,I11,5X,I8)
377
378	END SUBROUTINE obs_pre_pro
379
380	SUBROUTINE obs_pre_sla( sladata, sladatqc, ld_sla, ld_nea )
381	!!----------------------------------------------------------------------
382	!! * ROUTINE obs_pre_sla *
383	!!
384	!! ** Purpose : First level check and screening of SLA observations
385	!!
386	!! ** Method : First level check and screening of SLA observations
387	!!
388	!! ** Action :
389	!!
390	!! References :
391	!!
392	!! History :
393	!! ! 2007-03 (A. Weaver, K. Mogensen) Original
394	!! ! 2007-06 (K. Mogensen et al) Reject obs. near land.
395	!!----------------------------------------------------------------------
396	!! * Modules used
397	USE domstp ! Domain: set the time-step
398	USE par_oce ! Ocean parameters
399	USE dom_oce, ONLY : & ! Geographical information
400	& glamt, &
401	& gphit, &
402	& tmask
403	!! * Arguments
404	TYPE(obs_surf), INTENT(INOUT) :: sladata ! Full set of SLA data
405	TYPE(obs_surf), INTENT(INOUT) :: sladatqc ! Subset of SLA data not failing screening
406	LOGICAL, INTENT(IN) :: ld_sla ! Switch for SLA data
407	LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land
408	!! * Local declarations
409	INTEGER :: iyea0 ! Initial date
410	INTEGER :: imon0 ! - (year, month, day, hour, minute)
411	INTEGER :: iday0
412	INTEGER :: ihou0
413	INTEGER :: imin0
414	INTEGER :: icycle ! Current assimilation cycle
415	! Counters for observations that
416	INTEGER :: iotdobs ! - outside time domain
417	INTEGER :: iosdsobs ! - outside space domain
418	INTEGER :: ilansobs ! - within a model land cell
419	INTEGER :: inlasobs ! - close to land
420	INTEGER :: igrdobs ! - fail the grid search
421	INTEGER :: ibdysobs ! - close to open boundary
422	! Global counters for observations that
423	INTEGER :: iotdobsmpp ! - outside time domain
424	INTEGER :: iosdsobsmpp ! - outside space domain
425	INTEGER :: ilansobsmpp ! - within a model land cell
426	INTEGER :: inlasobsmpp ! - close to land
427	INTEGER :: igrdobsmpp ! - fail the grid search
428	INTEGER :: ibdysobsmpp ! - close to open boundary
429	LOGICAL, DIMENSION(:), ALLOCATABLE :: &
430	& llvalid ! SLA data selection
431	INTEGER :: jobs ! Obs. loop variable
432	INTEGER :: jstp ! Time loop variable
433	INTEGER :: inrc ! Time index variable
434	INTEGER :: irec ! Record index
435
436	IF(lwp) WRITE(numout,*)'obs_pre_sla : Preparing the SLA observations...'
437
438	! Initial date initialization (year, month, day, hour, minute)
439	iyea0 = ndate0 / 10000
440	imon0 = ( ndate0 - iyea0 * 10000 ) / 100
441	iday0 = ndate0 - iyea0 * 10000 - imon0 * 100
442	ihou0 = 0
443	imin0 = 0
444
445	icycle = no ! Assimilation cycle
446
447	! Diagnotics counters for various failures.
448
449	iotdobs = 0
450	igrdobs = 0
451	iosdsobs = 0
452	ilansobs = 0
453	inlasobs = 0
454	ibdysobs = 0
455
456	! -----------------------------------------------------------------------
457	! Find time coordinate for SLA data
458	! -----------------------------------------------------------------------
459
460	CALL obs_coo_tim( icycle, &
461	& iyea0, imon0, iday0, ihou0, imin0, &
462	& sladata%nsurf, sladata%nyea, sladata%nmon, &
463	& sladata%nday, sladata%nhou, sladata%nmin, &
464	& sladata%nqc, sladata%mstp, iotdobs )
465
466	CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
467
468	! -----------------------------------------------------------------------
469	! Check for SLA data failing the grid search
470	! -----------------------------------------------------------------------
471
472	CALL obs_coo_grd( sladata%nsurf, sladata%mi, sladata%mj, &
473	& sladata%nqc, igrdobs )
474
475	CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
476
477	! -----------------------------------------------------------------------
478	! Check for land points.
479	! -----------------------------------------------------------------------
480
481	CALL obs_coo_spc_2d( sladata%nsurf, &
482	& jpi, jpj, &
483	& sladata%mi, sladata%mj, &
484	& sladata%rlam, sladata%rphi, &
485	& glamt, gphit, &
486	& tmask(:,:,1), sladata%nqc, &
487	& iosdsobs, ilansobs, &
488	& inlasobs, ld_nea, &
489	& ibdysobs, ln_bound_reject )
490
491	CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
492	CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
493	CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
494	CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp )
495
496	! -----------------------------------------------------------------------
497	! Copy useful data from the sladata data structure to
498	! the sladatqc data structure
499	! -----------------------------------------------------------------------
500
501	! Allocate the selection arrays
502
503	ALLOCATE( llvalid(sladata%nsurf) )
504
505	! We want all data which has qc flags <= 10
506
507	llvalid(:) = ( sladata%nqc(:) <= 10 )
508
509	! The actual copying
510
511	CALL obs_surf_compress( sladata, sladatqc, .TRUE., numout, &
512	& lvalid=llvalid )
513
514	! Dellocate the selection arrays
515	DEALLOCATE( llvalid )
516
517	! -----------------------------------------------------------------------
518	! Print information about what observations are left after qc
519	! -----------------------------------------------------------------------
520
521	! Update the total observation counter array
522
523	IF(lwp) THEN
524	WRITE(numout,*)
525	WRITE(numout,*) 'obs_pre_sla :'
526	WRITE(numout,*) '~~~~~~~~~~~'
527	WRITE(numout,*)
528	WRITE(numout,*) ' SLA data outside time domain = ', &
529	& iotdobsmpp
530	WRITE(numout,*) ' Remaining SLA data that failed grid search = ', &
531	& igrdobsmpp
532	WRITE(numout,*) ' Remaining SLA data outside space domain = ', &
533	& iosdsobsmpp
534	WRITE(numout,*) ' Remaining SLA data at land points = ', &
535	& ilansobsmpp
536	IF (ld_nea) THEN
537	WRITE(numout,*) ' Remaining SLA data near land points (removed) = ', &
538	& inlasobsmpp
539	ELSE
540	WRITE(numout,*) ' Remaining SLA data near land points (kept) = ', &
541	& inlasobsmpp
542	ENDIF
543	WRITE(numout,*) ' Remaining SLA data near open boundary (removed) = ', &
544	& ibdysobsmpp
545	WRITE(numout,*) ' SLA data accepted = ', &
546	& sladatqc%nsurfmpp
547
548	WRITE(numout,*)
549	WRITE(numout,*) ' Number of observations per time step :'
550	WRITE(numout,*)
551	WRITE(numout,1997)
552	WRITE(numout,1998)
553	ENDIF
554
555	DO jobs = 1, sladatqc%nsurf
556	inrc = sladatqc%mstp(jobs) + 2 - nit000
557	sladatqc%nsstp(inrc) = sladatqc%nsstp(inrc) + 1
558	END DO
559
560	CALL obs_mpp_sum_integers( sladatqc%nsstp, sladatqc%nsstpmpp, &
561	& nitend - nit000 + 2 )
562
563	IF ( lwp ) THEN
564	DO jstp = nit000 - 1, nitend
565	inrc = jstp - nit000 + 2
566	WRITE(numout,1999) jstp, sladatqc%nsstpmpp(inrc)
567	END DO
568	ENDIF
569
570	!---------------------------------------------------------
571	! Record handling
572	!---------------------------------------------------------
573	! First count the number of records
574	sladatqc%nrec = 0
575	DO jstp = nit000 - 1, nitend
576	inrc = jstp - nit000 + 2
577	IF ( sladatqc%nsstpmpp(inrc) > 0 ) THEN
578	sladatqc%nrec = sladatqc%nrec + 1
579	ENDIF
580	END DO
581	! Allocate record data
582	ALLOCATE( &
583	& sladatqc%mrecstp(sladatqc%nrec) &
584	& )
585	! Finally save the time step corresponding to record rank
586	irec = 0
587	DO jstp = nit000 - 1, nitend
588	inrc = jstp - nit000 + 2
589	IF ( sladatqc%nsstpmpp(inrc) > 0 ) THEN
590	irec = irec + 1
591	sladatqc%mrecstp(irec) = inrc
592	ENDIF
593	IF ( lwp ) THEN
594	WRITE(numout,1999) inrc, sladatqc%nsstpmpp(inrc)
595	ENDIF
596	END DO
597
598	! Print record information
599	IF( lwp ) THEN
600	WRITE(numout,*)
601	WRITE(numout,2000)
602	WRITE(numout,2001)
603	DO irec = 1, sladatqc%nrec
604	WRITE(numout,1999) irec, sladatqc%mrecstp(irec)
605	END DO
606	ENDIF
607
608
609	1997 FORMAT(10X,'Time step',5X,'Sea level anomaly')
610	1998 FORMAT(10X,'---------',5X,'-----------------')
611	1999 FORMAT(10X,I9,5X,I17)
612	2000 FORMAT(15X,'Record',10X,'Time step')
613	2001 FORMAT(15X,'------',10X,'---------')
614
615	END SUBROUTINE obs_pre_sla
616
617	SUBROUTINE obs_pre_sst( sstdata, sstdatqc, ld_sst, ld_nea )
618	!!----------------------------------------------------------------------
619	!! * ROUTINE obs_pre_sst *
620	!!
621	!! ** Purpose : First level check and screening of SST observations
622	!!
623	!! ** Method : First level check and screening of SST observations
624	!!
625	!! ** Action :
626	!!
627	!! References :
628	!!
629	!! History :
630	!! ! 2007-03 (S. Ricci) SST data preparation
631	!!----------------------------------------------------------------------
632	!! * Modules used
633	USE domstp ! Domain: set the time-step
634	USE par_oce ! Ocean parameters
635	USE dom_oce, ONLY : & ! Geographical information
636	& glamt, &
637	& gphit, &
638	& tmask
639	!! * Arguments
640	TYPE(obs_surf), INTENT(INOUT) :: sstdata ! Full set of SST data
641	TYPE(obs_surf), INTENT(INOUT) :: sstdatqc ! Subset of SST data not failing screening
642	LOGICAL :: ld_sst ! Switch for SST data
643	LOGICAL :: ld_nea ! Switch for rejecting observation near land
644	!! * Local declarations
645	INTEGER :: iyea0 ! Initial date
646	INTEGER :: imon0 ! - (year, month, day, hour, minute)
647	INTEGER :: iday0
648	INTEGER :: ihou0
649	INTEGER :: imin0
650	INTEGER :: icycle ! Current assimilation cycle
651	! Counters for observations that
652	INTEGER :: iotdobs ! - outside time domain
653	INTEGER :: iosdsobs ! - outside space domain
654	INTEGER :: ilansobs ! - within a model land cell
655	INTEGER :: inlasobs ! - close to land
656	INTEGER :: igrdobs ! - fail the grid search
657	INTEGER :: ibdysobs ! - close to open boundary
658	! Global counters for observations that
659	INTEGER :: iotdobsmpp ! - outside time domain
660	INTEGER :: iosdsobsmpp ! - outside space domain
661	INTEGER :: ilansobsmpp ! - within a model land cell
662	INTEGER :: inlasobsmpp ! - close to land
663	INTEGER :: igrdobsmpp ! - fail the grid search
664	INTEGER :: ibdysobsmpp ! - close to open boundary
665	LOGICAL, DIMENSION(:), ALLOCATABLE :: &
666	& llvalid ! SST data selection
667	INTEGER :: jobs ! Obs. loop variable
668	INTEGER :: jstp ! Time loop variable
669	INTEGER :: inrc ! Time index variable
670	INTEGER :: irec ! Record index
671
672	IF(lwp) WRITE(numout,*)'obs_pre_sst : Preparing the SST observations...'
673
674	! Initial date initialization (year, month, day, hour, minute)
675	iyea0 = ndate0 / 10000
676	imon0 = ( ndate0 - iyea0 * 10000 ) / 100
677	iday0 = ndate0 - iyea0 * 10000 - imon0 * 100
678	ihou0 = 0
679	imin0 = 0
680
681	icycle = no ! Assimilation cycle
682
683	! Diagnotics counters for various failures.
684
685	iotdobs = 0
686	igrdobs = 0
687	iosdsobs = 0
688	ilansobs = 0
689	inlasobs = 0
690	ibdysobs = 0
691
692	! -----------------------------------------------------------------------
693	! Find time coordinate for SST data
694	! -----------------------------------------------------------------------
695
696	CALL obs_coo_tim( icycle, &
697	& iyea0, imon0, iday0, ihou0, imin0, &
698	& sstdata%nsurf, sstdata%nyea, sstdata%nmon, &
699	& sstdata%nday, sstdata%nhou, sstdata%nmin, &
700	& sstdata%nqc, sstdata%mstp, iotdobs )
701	CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
702	! -----------------------------------------------------------------------
703	! Check for SST data failing the grid search
704	! -----------------------------------------------------------------------
705
706	CALL obs_coo_grd( sstdata%nsurf, sstdata%mi, sstdata%mj, &
707	& sstdata%nqc, igrdobs )
708	CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
709
710	! -----------------------------------------------------------------------
711	! Check for land points.
712	! -----------------------------------------------------------------------
713
714	CALL obs_coo_spc_2d( sstdata%nsurf, &
715	& jpi, jpj, &
716	& sstdata%mi, sstdata%mj, &
717	& sstdata%rlam, sstdata%rphi, &
718	& glamt, gphit, &
719	& tmask(:,:,1), sstdata%nqc, &
720	& iosdsobs, ilansobs, &
721	& inlasobs, ld_nea, &
722	& ibdysobs, ln_bound_reject )
723
724	CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
725	CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
726	CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
727	CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp )
728
729	! -----------------------------------------------------------------------
730	! Copy useful data from the sstdata data structure to
731	! the sstdatqc data structure
732	! -----------------------------------------------------------------------
733
734	! Allocate the selection arrays
735
736	ALLOCATE( llvalid(sstdata%nsurf) )
737
738	! We want all data which has qc flags <= 0
739
740	llvalid(:) = ( sstdata%nqc(:) <= 10 )
741
742	! The actual copying
743
744	CALL obs_surf_compress( sstdata, sstdatqc, .TRUE., numout, &
745	& lvalid=llvalid )
746
747	! Dellocate the selection arrays
748	DEALLOCATE( llvalid )
749
750	! -----------------------------------------------------------------------
751	! Print information about what observations are left after qc
752	! -----------------------------------------------------------------------
753
754	! Update the total observation counter array
755
756	IF(lwp) THEN
757	WRITE(numout,*)
758	WRITE(numout,*) 'obs_pre_sst :'
759	WRITE(numout,*) '~~~~~~~~~~~'
760	WRITE(numout,*)
761	WRITE(numout,*) ' SST data outside time domain = ', &
762	& iotdobsmpp
763	WRITE(numout,*) ' Remaining SST data that failed grid search = ', &
764	& igrdobsmpp
765	WRITE(numout,*) ' Remaining SST data outside space domain = ', &
766	& iosdsobsmpp
767	WRITE(numout,*) ' Remaining SST data at land points = ', &
768	& ilansobsmpp
769	IF (ld_nea) THEN
770	WRITE(numout,*) ' Remaining SST data near land points (removed) = ', &
771	& inlasobsmpp
772	ELSE
773	WRITE(numout,*) ' Remaining SST data near land points (kept) = ', &
774	& inlasobsmpp
775	ENDIF
776	WRITE(numout,*) ' Remaining SST data near open boundary (removed) = ', &
777	& ibdysobsmpp
778	WRITE(numout,*) ' SST data accepted = ', &
779	& sstdatqc%nsurfmpp
780
781	WRITE(numout,*)
782	WRITE(numout,*) ' Number of observations per time step :'
783	WRITE(numout,*)
784	WRITE(numout,1997)
785	WRITE(numout,1998)
786	ENDIF
787
788	DO jobs = 1, sstdatqc%nsurf
789	inrc = sstdatqc%mstp(jobs) + 2 - nit000
790	sstdatqc%nsstp(inrc) = sstdatqc%nsstp(inrc) + 1
791	END DO
792
793	CALL obs_mpp_sum_integers( sstdatqc%nsstp, sstdatqc%nsstpmpp, &
794	& nitend - nit000 + 2 )
795
796	IF ( lwp ) THEN
797	DO jstp = nit000 - 1, nitend
798	inrc = jstp - nit000 + 2
799	WRITE(numout,1999) jstp, sstdatqc%nsstpmpp(inrc)
800	END DO
801	ENDIF
802
803	!---------------------------------------------------------
804	! Record handling
805	!---------------------------------------------------------
806	! First count the number of records
807	sstdatqc%nrec = 0
808	DO jstp = nit000 - 1, nitend
809	inrc = jstp - nit000 + 2
810	IF ( sstdatqc%nsstpmpp(inrc) > 0 ) THEN
811	sstdatqc%nrec = sstdatqc%nrec + 1
812	ENDIF
813	END DO
814	! Allocate record data
815	ALLOCATE( &
816	& sstdatqc%mrecstp(sstdatqc%nrec) &
817	& )
818	! Finally save the time step corresponding to record rank
819	irec = 0
820	DO jstp = nit000 - 1, nitend
821	inrc = jstp - nit000 + 2
822	IF ( sstdatqc%nsstpmpp(inrc) > 0 ) THEN
823	irec = irec + 1
824	sstdatqc%mrecstp(irec) = inrc
825	ENDIF
826	IF ( lwp ) THEN
827	WRITE(numout,1999) jstp, sstdatqc%nsstpmpp(inrc)
828	ENDIF
829	END DO
830
831	! Print record information
832	IF( lwp ) THEN
833	WRITE(numout,*)
834	WRITE(numout,2000)
835	WRITE(numout,2001)
836	DO irec = 1, sstdatqc%nrec
837	WRITE(numout,1999) irec, sstdatqc%mrecstp(irec) - 1
838	END DO
839	ENDIF
840
841	1997 FORMAT(10X,'Time step',5X,'Sea surface temperature')
842	1998 FORMAT(10X,'---------',5X,'-----------------')
843	1999 FORMAT(10X,I9,5X,I17)
844	2000 FORMAT(15X,'Record',10X,'Time step')
845	2001 FORMAT(15X,'------',10X,'---------')
846
847	END SUBROUTINE obs_pre_sst
848
849	SUBROUTINE obs_pre_seaice( seaicedata, seaicedatqc, ld_seaice, ld_nea )
850	!!----------------------------------------------------------------------
851	!! * ROUTINE obs_pre_seaice *
852	!!
853	!! ** Purpose : First level check and screening of Sea Ice observations
854	!!
855	!! ** Method : First level check and screening of Sea Ice observations
856	!!
857	!! ** Action :
858	!!
859	!! References :
860	!!
861	!! History :
862	!! ! 2007-11 (D. Lea) based on obs_pre_sst
863	!!----------------------------------------------------------------------
864	!! * Modules used
865	USE domstp ! Domain: set the time-step
866	USE par_oce ! Ocean parameters
867	USE dom_oce, ONLY : & ! Geographical information
868	& glamt, &
869	& gphit, &
870	& tmask
871	!! * Arguments
872	TYPE(obs_surf), INTENT(INOUT) :: seaicedata ! Full set of Sea Ice data
873	TYPE(obs_surf), INTENT(INOUT) :: seaicedatqc ! Subset of sea ice data not failing screening
874	LOGICAL :: ld_seaice ! Switch for sea ice data
875	LOGICAL :: ld_nea ! Switch for rejecting observation near land
876	!! * Local declarations
877	INTEGER :: iyea0 ! Initial date
878	INTEGER :: imon0 ! - (year, month, day, hour, minute)
879	INTEGER :: iday0
880	INTEGER :: ihou0
881	INTEGER :: imin0
882	INTEGER :: icycle ! Current assimilation cycle
883	! Counters for observations that
884	INTEGER :: iotdobs ! - outside time domain
885	INTEGER :: iosdsobs ! - outside space domain
886	INTEGER :: ilansobs ! - within a model land cell
887	INTEGER :: inlasobs ! - close to land
888	INTEGER :: igrdobs ! - fail the grid search
889	INTEGER :: ibdysobs ! - close to open boundary
890	! Global counters for observations that
891	INTEGER :: iotdobsmpp ! - outside time domain
892	INTEGER :: iosdsobsmpp ! - outside space domain
893	INTEGER :: ilansobsmpp ! - within a model land cell
894	INTEGER :: inlasobsmpp ! - close to land
895	INTEGER :: igrdobsmpp ! - fail the grid search
896	INTEGER :: ibdysobsmpp ! - close to open boundary
897	LOGICAL, DIMENSION(:), ALLOCATABLE :: &
898	& llvalid ! data selection
899	INTEGER :: jobs ! Obs. loop variable
900	INTEGER :: jstp ! Time loop variable
901	INTEGER :: inrc ! Time index variable
902	INTEGER :: irec ! Record index
903
904	IF (lwp) WRITE(numout,*)'obs_pre_seaice : Preparing the sea ice observations...'
905
906	! Initial date initialization (year, month, day, hour, minute)
907	iyea0 = ndate0 / 10000
908	imon0 = ( ndate0 - iyea0 * 10000 ) / 100
909	iday0 = ndate0 - iyea0 * 10000 - imon0 * 100
910	ihou0 = 0
911	imin0 = 0
912
913	icycle = no ! Assimilation cycle
914
915	! Diagnotics counters for various failures.
916
917	iotdobs = 0
918	igrdobs = 0
919	iosdsobs = 0
920	ilansobs = 0
921	inlasobs = 0
922	ibdysobs = 0
923
924	! -----------------------------------------------------------------------
925	! Find time coordinate for sea ice data
926	! -----------------------------------------------------------------------
927
928	CALL obs_coo_tim( icycle, &
929	& iyea0, imon0, iday0, ihou0, imin0, &
930	& seaicedata%nsurf, seaicedata%nyea, seaicedata%nmon, &
931	& seaicedata%nday, seaicedata%nhou, seaicedata%nmin, &
932	& seaicedata%nqc, seaicedata%mstp, iotdobs )
933	CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
934	! -----------------------------------------------------------------------
935	! Check for sea ice data failing the grid search
936	! -----------------------------------------------------------------------
937
938	CALL obs_coo_grd( seaicedata%nsurf, seaicedata%mi, seaicedata%mj, &
939	& seaicedata%nqc, igrdobs )
940	CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
941
942	! -----------------------------------------------------------------------
943	! Check for land points.
944	! -----------------------------------------------------------------------
945
946	CALL obs_coo_spc_2d( seaicedata%nsurf, &
947	& jpi, jpj, &
948	& seaicedata%mi, seaicedata%mj, &
949	& seaicedata%rlam, seaicedata%rphi, &
950	& glamt, gphit, &
951	& tmask(:,:,1), seaicedata%nqc, &
952	& iosdsobs, ilansobs, &
953	& inlasobs, ld_nea, &
954	& ibdysobs, ln_bound_reject )
955
956	CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
957	CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
958	CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
959	CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp )
960
961	! -----------------------------------------------------------------------
962	! Copy useful data from the seaicedata data structure to
963	! the seaicedatqc data structure
964	! -----------------------------------------------------------------------
965
966	! Allocate the selection arrays
967
968	ALLOCATE( llvalid(seaicedata%nsurf) )
969
970	! We want all data which has qc flags <= 0
971
972	llvalid(:) = ( seaicedata%nqc(:) <= 10 )
973
974	! The actual copying
975
976	CALL obs_surf_compress( seaicedata, seaicedatqc, .TRUE., numout, &
977	& lvalid=llvalid )
978
979	! Dellocate the selection arrays
980	DEALLOCATE( llvalid )
981
982	! -----------------------------------------------------------------------
983	! Print information about what observations are left after qc
984	! -----------------------------------------------------------------------
985
986	! Update the total observation counter array
987
988	IF(lwp) THEN
989	WRITE(numout,*)
990	WRITE(numout,*) 'obs_pre_seaice :'
991	WRITE(numout,*) '~~~~~~~~~~~'
992	WRITE(numout,*)
993	WRITE(numout,*) ' Sea ice data outside time domain = ', &
994	& iotdobsmpp
995	WRITE(numout,*) ' Remaining sea ice data that failed grid search = ', &
996	& igrdobsmpp
997	WRITE(numout,*) ' Remaining sea ice data outside space domain = ', &
998	& iosdsobsmpp
999	WRITE(numout,*) ' Remaining sea ice data at land points = ', &
1000	& ilansobsmpp
1001	IF (ld_nea) THEN
1002	WRITE(numout,*) ' Remaining sea ice data near land points (removed) = ', &
1003	& inlasobsmpp
1004	ELSE
1005	WRITE(numout,*) ' Remaining sea ice data near land points (kept) = ', &
1006	& inlasobsmpp
1007	ENDIF
1008	WRITE(numout,*) ' Remaining sea ice data near open boundary (removed) = ', &
1009	& ibdysobsmpp
1010	WRITE(numout,*) ' Sea ice data accepted = ', &
1011	& seaicedatqc%nsurfmpp
1012
1013	WRITE(numout,*)
1014	WRITE(numout,*) ' Number of observations per time step :'
1015	WRITE(numout,*)
1016	WRITE(numout,1997)
1017	WRITE(numout,1998)
1018	ENDIF
1019
1020	DO jobs = 1, seaicedatqc%nsurf
1021	inrc = seaicedatqc%mstp(jobs) + 2 - nit000
1022	seaicedatqc%nsstp(inrc) = seaicedatqc%nsstp(inrc) + 1
1023	END DO
1024
1025	CALL obs_mpp_sum_integers( seaicedatqc%nsstp, seaicedatqc%nsstpmpp, &
1026	& nitend - nit000 + 2 )
1027
1028	IF ( lwp ) THEN
1029	DO jstp = nit000 - 1, nitend
1030	inrc = jstp - nit000 + 2
1031	WRITE(numout,1999) jstp, seaicedatqc%nsstpmpp(inrc)
1032	END DO
1033	ENDIF
1034
1035	!---------------------------------------------------------
1036	! Record handling
1037	!---------------------------------------------------------
1038	! First count the number of records
1039	seaicedatqc%nrec = 0
1040	DO jstp = nit000 - 1, nitend
1041	inrc = jstp - nit000 + 2
1042	IF ( seaicedatqc%nsstpmpp(inrc) > 0 ) THEN
1043	seaicedatqc%nrec = seaicedatqc%nrec + 1
1044	ENDIF
1045	END DO
1046	! Allocate record data
1047	ALLOCATE( &
1048	& seaicedatqc%mrecstp(seaicedatqc%nrec) &
1049	& )
1050	! Finally save the time step corresponding to record rank
1051	irec = 0
1052	DO jstp = nit000 - 1, nitend
1053	inrc = jstp - nit000 + 2
1054	IF ( seaicedatqc%nsstpmpp(inrc) > 0 ) THEN
1055	irec = irec + 1
1056	seaicedatqc%mrecstp(irec) = inrc
1057	ENDIF
1058	IF ( lwp ) THEN
1059	WRITE(numout,1999) inrc, seaicedatqc%nsstpmpp(inrc)
1060	ENDIF
1061	END DO
1062
1063	! Print record information
1064	IF( lwp ) THEN
1065	WRITE(numout,*)
1066	WRITE(numout,2000)
1067	WRITE(numout,2001)
1068	DO irec = 1, seaicedatqc%nrec
1069	WRITE(numout,1999) irec, seaicedatqc%mrecstp(irec)
1070	END DO
1071	ENDIF
1072
1073	1997 FORMAT(10X,'Time step',5X,'Sea ice data ')
1074	1998 FORMAT(10X,'---------',5X,'-----------------')
1075	1999 FORMAT(10X,I9,5X,I17)
1076	2000 FORMAT(15X,'Record',10X,'Time step')
1077	2001 FORMAT(15X,'------',10X,'---------')
1078
1079	END SUBROUTINE obs_pre_seaice
1080
1081	SUBROUTINE obs_pre_vel( profdata, prodatqc, ld_vel3d, ld_nea, ld_dailyav )
1082	!!----------------------------------------------------------------------
1083	!! * ROUTINE obs_pre_taovel *
1084	!!
1085	!! ** Purpose : First level check and screening of U and V profiles
1086	!!
1087	!! ** Method : First level check and screening of U and V profiles
1088	!!
1089	!! History :
1090	!! ! 2007-06 (K. Mogensen) original : T and S profile data
1091	!! ! 2008-09 (M. Valdivieso) : TAO velocity data
1092	!! ! 2009-01 (K. Mogensen) : New feedback strictuer
1093	!!
1094	!!----------------------------------------------------------------------
1095	!! * Modules used
1096	USE domstp ! Domain: set the time-step
1097	USE par_oce ! Ocean parameters
1098	USE dom_oce, ONLY : & ! Geographical information
1099	& glamt, glamu, glamv, &
1100	& gphit, gphiu, gphiv, &
1101	& gdept_1d, &
1102	& tmask, umask, vmask
1103	!! * Arguments
1104	TYPE(obs_prof), INTENT(INOUT) :: profdata ! Full set of profile data
1105	TYPE(obs_prof), INTENT(INOUT) :: prodatqc ! Subset of profile data not failing screening
1106	LOGICAL, INTENT(IN) :: ld_vel3d ! Switch for zonal and meridional velocity components
1107	LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land
1108	LOGICAL, INTENT(IN) :: ld_dailyav ! Switch for daily average data
1109	!! * Local declarations
1110	INTEGER :: iyea0 ! Initial date
1111	INTEGER :: imon0 ! - (year, month, day, hour, minute)
1112	INTEGER :: iday0
1113	INTEGER :: ihou0
1114	INTEGER :: imin0
1115	INTEGER :: icycle ! Current assimilation cycle
1116	! Counters for observations that
1117	INTEGER :: iotdobs ! - outside time domain
1118	INTEGER :: iosduobs ! - outside space domain (zonal velocity component)
1119	INTEGER :: iosdvobs ! - outside space domain (meridional velocity component)
1120	INTEGER :: ilanuobs ! - within a model land cell (zonal velocity component)
1121	INTEGER :: ilanvobs ! - within a model land cell (meridional velocity component)
1122	INTEGER :: inlauobs ! - close to land (zonal velocity component)
1123	INTEGER :: inlavobs ! - close to land (meridional velocity component)
1124	INTEGER :: igrdobs ! - fail the grid search
1125	INTEGER :: ibdyuobs ! - close to open boundary
1126	INTEGER :: ibdyvobs ! - close to open boundary
1127	INTEGER :: iuvchku ! - reject u if v rejected and vice versa
1128	INTEGER :: iuvchkv !
1129	! Global counters for observations that
1130	INTEGER :: iotdobsmpp ! - outside time domain
1131	INTEGER :: iosduobsmpp ! - outside space domain (zonal velocity component)
1132	INTEGER :: iosdvobsmpp ! - outside space domain (meridional velocity component)
1133	INTEGER :: ilanuobsmpp ! - within a model land cell (zonal velocity component)
1134	INTEGER :: ilanvobsmpp ! - within a model land cell (meridional velocity component)
1135	INTEGER :: inlauobsmpp ! - close to land (zonal velocity component)
1136	INTEGER :: inlavobsmpp ! - close to land (meridional velocity component)
1137	INTEGER :: igrdobsmpp ! - fail the grid search
1138	INTEGER :: ibdyuobsmpp ! - close to open boundary
1139	INTEGER :: ibdyvobsmpp ! - close to open boundary
1140	INTEGER :: iuvchkumpp ! - reject u if v rejected and vice versa
1141	INTEGER :: iuvchkvmpp !
1142	TYPE(obs_prof_valid) :: llvalid ! Profile selection
1143	TYPE(obs_prof_valid), DIMENSION(profdata%nvar) :: &
1144	& llvvalid ! U,V selection
1145	INTEGER :: jvar ! Variable loop variable
1146	INTEGER :: jobs ! Obs. loop variable
1147	INTEGER :: jstp ! Time loop variable
1148	INTEGER :: inrc ! Time index variable
1149
1150	IF(lwp) WRITE(numout,*)'obs_pre_vel: Preparing the velocity profile data'
1151
1152	! Initial date initialization (year, month, day, hour, minute)
1153	iyea0 = ndate0 / 10000
1154	imon0 = ( ndate0 - iyea0 * 10000 ) / 100
1155	iday0 = ndate0 - iyea0 * 10000 - imon0 * 100
1156	ihou0 = 0
1157	imin0 = 0
1158
1159	icycle = no ! Assimilation cycle
1160
1161	! Diagnotics counters for various failures.
1162
1163	iotdobs = 0
1164	igrdobs = 0
1165	iosduobs = 0
1166	iosdvobs = 0
1167	ilanuobs = 0
1168	ilanvobs = 0
1169	inlauobs = 0
1170	inlavobs = 0
1171	ibdyuobs = 0
1172	ibdyvobs = 0
1173	iuvchku = 0
1174	iuvchkv = 0
1175
1176	! -----------------------------------------------------------------------
1177	! Find time coordinate for profiles
1178	! -----------------------------------------------------------------------
1179
1180	CALL obs_coo_tim_prof( icycle, &
1181	& iyea0, imon0, iday0, ihou0, imin0, &
1182	& profdata%nprof, profdata%nyea, profdata%nmon, &
1183	& profdata%nday, profdata%nhou, profdata%nmin, &
1184	& profdata%ntyp, profdata%nqc, profdata%mstp, &
1185	& iotdobs, ld_dailyav = ld_dailyav )
1186
1187	CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
1188
1189	! -----------------------------------------------------------------------
1190	! Check for profiles failing the grid search
1191	! -----------------------------------------------------------------------
1192
1193	CALL obs_coo_grd( profdata%nprof, profdata%mi(:,1), profdata%mj(:,1), &
1194	& profdata%nqc, igrdobs )
1195	CALL obs_coo_grd( profdata%nprof, profdata%mi(:,2), profdata%mj(:,2), &
1196	& profdata%nqc, igrdobs )
1197
1198	CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
1199
1200	! -----------------------------------------------------------------------
1201	! Reject all observations for profiles with nqc > 10
1202	! -----------------------------------------------------------------------
1203
1204	CALL obs_pro_rej( profdata )
1205
1206	! -----------------------------------------------------------------------
1207	! Check for land points. This includes points below the model
1208	! bathymetry so this is done for every point in the profile
1209	! -----------------------------------------------------------------------
1210
1211	! Zonal Velocity Component
1212
1213	CALL obs_coo_spc_3d( profdata%nprof, profdata%nvprot(1), &
1214	& profdata%npvsta(:,1), profdata%npvend(:,1), &
1215	& jpi, jpj, &
1216	& jpk, &
1217	& profdata%mi, profdata%mj, &
1218	& profdata%var(1)%mvk, &
1219	& profdata%rlam, profdata%rphi, &
1220	& profdata%var(1)%vdep, &
1221	& glamu, gphiu, &
1222	& gdept_1d, umask, &
1223	& profdata%nqc, profdata%var(1)%nvqc, &
1224	& iosduobs, ilanuobs, &
1225	& inlauobs, ld_nea, &
1226	& ibdyuobs, ln_bound_reject )
1227
1228	CALL obs_mpp_sum_integer( iosduobs, iosduobsmpp )
1229	CALL obs_mpp_sum_integer( ilanuobs, ilanuobsmpp )
1230	CALL obs_mpp_sum_integer( inlauobs, inlauobsmpp )
1231	CALL obs_mpp_sum_integer( ibdyuobs, ibdyuobsmpp )
1232
1233	! Meridional Velocity Component
1234
1235	CALL obs_coo_spc_3d( profdata%nprof, profdata%nvprot(2), &
1236	& profdata%npvsta(:,2), profdata%npvend(:,2), &
1237	& jpi, jpj, &
1238	& jpk, &
1239	& profdata%mi, profdata%mj, &
1240	& profdata%var(2)%mvk, &
1241	& profdata%rlam, profdata%rphi, &
1242	& profdata%var(2)%vdep, &
1243	& glamv, gphiv, &
1244	& gdept_1d, vmask, &
1245	& profdata%nqc, profdata%var(2)%nvqc, &
1246	& iosdvobs, ilanvobs, &
1247	& inlavobs, ld_nea, &
1248	& ibdyvobs, ln_bound_reject )
1249
1250	CALL obs_mpp_sum_integer( iosdvobs, iosdvobsmpp )
1251	CALL obs_mpp_sum_integer( ilanvobs, ilanvobsmpp )
1252	CALL obs_mpp_sum_integer( inlavobs, inlavobsmpp )
1253	CALL obs_mpp_sum_integer( ibdyvobs, ibdyvobsmpp )
1254
1255	! -----------------------------------------------------------------------
1256	! Reject u if v is rejected and vice versa
1257	! -----------------------------------------------------------------------
1258
1259	CALL obs_uv_rej( profdata, iuvchku, iuvchkv )
1260	CALL obs_mpp_sum_integer( iuvchku, iuvchkumpp )
1261	CALL obs_mpp_sum_integer( iuvchkv, iuvchkvmpp )
1262
1263	! -----------------------------------------------------------------------
1264	! Copy useful data from the profdata data structure to
1265	! the prodatqc data structure
1266	! -----------------------------------------------------------------------
1267
1268	! Allocate the selection arrays
1269
1270	ALLOCATE( llvalid%luse(profdata%nprof) )
1271	DO jvar = 1,profdata%nvar
1272	ALLOCATE( llvvalid(jvar)%luse(profdata%nvprot(jvar)) )
1273	END DO
1274
1275	! We want all data which has qc flags = 0
1276
1277	llvalid%luse(:) = ( profdata%nqc(:) <= 10 )
1278	DO jvar = 1,profdata%nvar
1279	llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= 10 )
1280	END DO
1281
1282	! The actual copying
1283
1284	CALL obs_prof_compress( profdata, prodatqc, .TRUE., numout, &
1285	& lvalid=llvalid, lvvalid=llvvalid )
1286
1287	! Dellocate the selection arrays
1288	DEALLOCATE( llvalid%luse )
1289	DO jvar = 1,profdata%nvar
1290	DEALLOCATE( llvvalid(jvar)%luse )
1291	END DO
1292
1293	! -----------------------------------------------------------------------
1294	! Print information about what observations are left after qc
1295	! -----------------------------------------------------------------------
1296
1297	! Update the total observation counter array
1298
1299	IF(lwp) THEN
1300	WRITE(numout,*)
1301	WRITE(numout,*) 'obs_pre_vel :'
1302	WRITE(numout,*) '~~~~~~~~~~~'
1303	WRITE(numout,*)
1304	WRITE(numout,*) ' Profiles outside time domain = ', &
1305	& iotdobsmpp
1306	WRITE(numout,*) ' Remaining profiles that failed grid search = ', &
1307	& igrdobsmpp
1308	WRITE(numout,*) ' Remaining U data outside space domain = ', &
1309	& iosduobsmpp
1310	WRITE(numout,*) ' Remaining U data at land points = ', &
1311	& ilanuobsmpp
1312	IF (ld_nea) THEN
1313	WRITE(numout,*) ' Remaining U data near land points (removed) = ',&
1314	& inlauobsmpp
1315	ELSE
1316	WRITE(numout,*) ' Remaining U data near land points (kept) = ',&
1317	& inlauobsmpp
1318	ENDIF
1319	WRITE(numout,*) ' Remaining U data near open boundary (removed) = ', &
1320	& ibdyuobsmpp
1321	WRITE(numout,*) ' U observation rejected since V rejected = ', &
1322	& iuvchku
1323	WRITE(numout,*) ' U data accepted = ', &
1324	& prodatqc%nvprotmpp(1)
1325	WRITE(numout,*) ' Remaining V data outside space domain = ', &
1326	& iosdvobsmpp
1327	WRITE(numout,*) ' Remaining V data at land points = ', &
1328	& ilanvobsmpp
1329	IF (ld_nea) THEN
1330	WRITE(numout,*) ' Remaining V data near land points (removed) = ',&
1331	& inlavobsmpp
1332	ELSE
1333	WRITE(numout,*) ' Remaining V data near land points (kept) = ',&
1334	& inlavobsmpp
1335	ENDIF
1336	WRITE(numout,*) ' Remaining V data near open boundary (removed) = ', &
1337	& ibdyvobsmpp
1338	WRITE(numout,*) ' V observation rejected since U rejected = ', &
1339	& iuvchkv
1340	WRITE(numout,*) ' V data accepted = ', &
1341	& prodatqc%nvprotmpp(2)
1342
1343	WRITE(numout,*)
1344	WRITE(numout,*) ' Number of observations per time step :'
1345	WRITE(numout,*)
1346	WRITE(numout,997)
1347	WRITE(numout,998)
1348	ENDIF
1349
1350	DO jobs = 1, prodatqc%nprof
1351	inrc = prodatqc%mstp(jobs) + 2 - nit000
1352	prodatqc%npstp(inrc) = prodatqc%npstp(inrc) + 1
1353	DO jvar = 1, prodatqc%nvar
1354	IF ( prodatqc%npvend(jobs,jvar) > 0 ) THEN
1355	prodatqc%nvstp(inrc,jvar) = prodatqc%nvstp(inrc,jvar) + &
1356	& ( prodatqc%npvend(jobs,jvar) - &
1357	& prodatqc%npvsta(jobs,jvar) + 1 )
1358	ENDIF
1359	END DO
1360	END DO
1361
1362
1363	CALL obs_mpp_sum_integers( prodatqc%npstp, prodatqc%npstpmpp, &
1364	& nitend - nit000 + 2 )
1365	DO jvar = 1, prodatqc%nvar
1366	CALL obs_mpp_sum_integers( prodatqc%nvstp(:,jvar), &
1367	& prodatqc%nvstpmpp(:,jvar), &
1368	& nitend - nit000 + 2 )
1369	END DO
1370
1371	IF ( lwp ) THEN
1372	DO jstp = nit000 - 1, nitend
1373	inrc = jstp - nit000 + 2
1374	WRITE(numout,999) jstp, prodatqc%npstpmpp(inrc), &
1375	& prodatqc%nvstpmpp(inrc,1), &
1376	& prodatqc%nvstpmpp(inrc,2)
1377	END DO
1378	ENDIF
1379
1380	997 FORMAT(10X,'Time step',5X,'Profiles',5X,'Zonal Comp.',5X,'Meridional Comp.')
1381	998 FORMAT(10X,'---------',5X,'--------',5X,'-----------',5X,'----------------')
1382	999 FORMAT(10X,I9,5X,I8,5X,I11,5X,I8)
1383
1384	END SUBROUTINE obs_pre_vel
1385
1386	SUBROUTINE obs_pre_logchl( logchldata, logchldatqc, ld_logchl, ld_nea )
1387	!!----------------------------------------------------------------------
1388	!! * ROUTINE obs_pre_logchl *
1389	!!
1390	!! ** Purpose : First level check and screening of logchl observations
1391	!!
1392	!! ** Method : First level check and screening of logchl observations
1393	!!
1394	!! ** Action :
1395	!!
1396	!! References :
1397	!!
1398	!! History :
1399	!!----------------------------------------------------------------------
1400	!! * Modules used
1401	USE domstp ! Domain: set the time-step
1402	USE par_oce ! Ocean parameters
1403	USE dom_oce, ONLY : & ! Geographical information
1404	& glamt, &
1405	& gphit, &
1406	& tmask
1407	!! * Arguments
1408	TYPE(obs_surf), INTENT(INOUT) :: logchldata ! Full set of logchl data
1409	TYPE(obs_surf), INTENT(INOUT) :: logchldatqc ! Subset of logchl data not failing screening
1410	LOGICAL :: ld_logchl ! Switch for logchl data
1411	LOGICAL :: ld_nea ! Switch for rejecting observation near land
1412	!! * Local declarations
1413	INTEGER :: iyea0 ! Initial date
1414	INTEGER :: imon0 ! - (year, month, day, hour, minute)
1415	INTEGER :: iday0
1416	INTEGER :: ihou0
1417	INTEGER :: imin0
1418	INTEGER :: icycle ! Current assimilation cycle
1419	! Counters for observations that
1420	INTEGER :: iotdobs ! - outside time domain
1421	INTEGER :: iosdsobs ! - outside space domain
1422	INTEGER :: ilansobs ! - within a model land cell
1423	INTEGER :: inlasobs ! - close to land
1424	INTEGER :: igrdobs ! - fail the grid search
1425	INTEGER :: ibdysobs ! - close to open boundary
1426	! Global counters for observations that
1427	INTEGER :: iotdobsmpp ! - outside time domain
1428	INTEGER :: iosdsobsmpp ! - outside space domain
1429	INTEGER :: ilansobsmpp ! - within a model land cell
1430	INTEGER :: inlasobsmpp ! - close to land
1431	INTEGER :: igrdobsmpp ! - fail the grid search
1432	INTEGER :: ibdysobsmpp ! - close to open boundary
1433	LOGICAL, DIMENSION(:), ALLOCATABLE :: &
1434	& llvalid ! data selection
1435	INTEGER :: jobs ! Obs. loop variable
1436	INTEGER :: jstp ! Time loop variable
1437	INTEGER :: inrc ! Time index variable
1438	INTEGER :: irec ! Record index
1439
1440	IF (lwp) WRITE(numout,*)'obs_pre_logchl : Preparing the logchl observations...'
1441
1442	! Initial date initialization (year, month, day, hour, minute)
1443	iyea0 = ndate0 / 10000
1444	imon0 = ( ndate0 - iyea0 * 10000 ) / 100
1445	iday0 = ndate0 - iyea0 * 10000 - imon0 * 100
1446	ihou0 = 0
1447	imin0 = 0
1448
1449	icycle = no ! Assimilation cycle
1450
1451	! Diagnostics counters for various failures.
1452
1453	iotdobs = 0
1454	igrdobs = 0
1455	iosdsobs = 0
1456	ilansobs = 0
1457	inlasobs = 0
1458	ibdysobs = 0
1459
1460	! -----------------------------------------------------------------------
1461	! Find time coordinate for logchl data
1462	! -----------------------------------------------------------------------
1463
1464	CALL obs_coo_tim( icycle, &
1465	& iyea0, imon0, iday0, ihou0, imin0, &
1466	& logchldata%nsurf, logchldata%nyea, logchldata%nmon, &
1467	& logchldata%nday, logchldata%nhou, logchldata%nmin, &
1468	& logchldata%nqc, logchldata%mstp, iotdobs )
1469	CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
1470	! -----------------------------------------------------------------------
1471	! Check for logchl data failing the grid search
1472	! -----------------------------------------------------------------------
1473
1474	CALL obs_coo_grd( logchldata%nsurf, logchldata%mi, logchldata%mj, &
1475	& logchldata%nqc, igrdobs )
1476	CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
1477
1478	! -----------------------------------------------------------------------
1479	! Check for land points.
1480	! -----------------------------------------------------------------------
1481
1482	CALL obs_coo_spc_2d( logchldata%nsurf, &
1483	& jpi, jpj, &
1484	& logchldata%mi, logchldata%mj, &
1485	& logchldata%rlam, logchldata%rphi, &
1486	& glamt, gphit, &
1487	& tmask(:,:,1), logchldata%nqc, &
1488	& iosdsobs, ilansobs, &
1489	& inlasobs, ld_nea, &
1490	& ibdysobs, ln_bound_reject )
1491
1492	CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
1493	CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
1494	CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
1495	CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp )
1496
1497	! -----------------------------------------------------------------------
1498	! Copy useful data from the logchldata data structure to
1499	! the logchldatqc data structure
1500	! -----------------------------------------------------------------------
1501
1502	! Allocate the selection arrays
1503
1504	ALLOCATE( llvalid(logchldata%nsurf) )
1505
1506	! We want all data which has qc flags <= 0
1507
1508	llvalid(:) = ( logchldata%nqc(:) <= 10 )
1509
1510	! The actual copying
1511
1512	CALL obs_surf_compress( logchldata, logchldatqc, .TRUE., numout, &
1513	& lvalid=llvalid )
1514
1515	! Dellocate the selection arrays
1516	DEALLOCATE( llvalid )
1517
1518	! -----------------------------------------------------------------------
1519	! Print information about what observations are left after qc
1520	! -----------------------------------------------------------------------
1521
1522	! Update the total observation counter array
1523
1524	IF(lwp) THEN
1525	WRITE(numout,*)
1526	WRITE(numout,*) 'obs_pre_logchl :'
1527	WRITE(numout,*) '~~~~~~~~~~~'
1528	WRITE(numout,*)
1529	WRITE(numout,*) ' logchl data outside time domain = ', &
1530	& iotdobsmpp
1531	WRITE(numout,*) ' Remaining logchl data that failed grid search = ', &
1532	& igrdobsmpp
1533	WRITE(numout,*) ' Remaining logchl data outside space domain = ', &
1534	& iosdsobsmpp
1535	WRITE(numout,*) ' Remaining logchl data at land points = ', &
1536	& ilansobsmpp
1537	IF (ld_nea) THEN
1538	WRITE(numout,*) ' Remaining logchl data near land points (removed) = ', &
1539	& inlasobsmpp
1540	ELSE
1541	WRITE(numout,*) ' Remaining logchl data near land points (kept) = ', &
1542	& inlasobsmpp
1543	ENDIF
1544	WRITE(numout,*) ' Remaining logchl data near open boundary (removed) = ', &
1545	& ibdysobsmpp
1546	WRITE(numout,*) ' logchl data accepted = ', &
1547	& logchldatqc%nsurfmpp
1548
1549	WRITE(numout,*)
1550	WRITE(numout,*) ' Number of observations per time step :'
1551	WRITE(numout,*)
1552	WRITE(numout,1997)
1553	WRITE(numout,1998)
1554	ENDIF
1555
1556	DO jobs = 1, logchldatqc%nsurf
1557	inrc = logchldatqc%mstp(jobs) + 2 - nit000
1558	logchldatqc%nsstp(inrc) = logchldatqc%nsstp(inrc) + 1
1559	END DO
1560
1561	CALL obs_mpp_sum_integers( logchldatqc%nsstp, logchldatqc%nsstpmpp, &
1562	& nitend - nit000 + 2 )
1563
1564	IF ( lwp ) THEN
1565	DO jstp = nit000 - 1, nitend
1566	inrc = jstp - nit000 + 2
1567	WRITE(numout,1999) jstp, logchldatqc%nsstpmpp(inrc)
1568	END DO
1569	ENDIF
1570
1571	1997 FORMAT(10X,'Time step',5X,'logchl data')
1572	1998 FORMAT(10X,'---------',5X,'------------')
1573	1999 FORMAT(10X,I9,5X,I17)
1574
1575	END SUBROUTINE obs_pre_logchl
1576
1577	SUBROUTINE obs_pre_spm( spmdata, spmdatqc, ld_spm, ld_nea )
1578	!!----------------------------------------------------------------------
1579	!! * ROUTINE obs_pre_spm *
1580	!!
1581	!! ** Purpose : First level check and screening of spm observations
1582	!!
1583	!! ** Method : First level check and screening of spm observations
1584	!!
1585	!! ** Action :
1586	!!
1587	!! References :
1588	!!
1589	!! History :
1590	!!----------------------------------------------------------------------
1591	!! * Modules used
1592	USE domstp ! Domain: set the time-step
1593	USE par_oce ! Ocean parameters
1594	USE dom_oce, ONLY : & ! Geographical information
1595	& glamt, &
1596	& gphit, &
1597	& tmask
1598	!! * Arguments
1599	TYPE(obs_surf), INTENT(INOUT) :: spmdata ! Full set of spm data
1600	TYPE(obs_surf), INTENT(INOUT) :: spmdatqc ! Subset of spm data not failing screening
1601	LOGICAL :: ld_spm ! Switch for spm data
1602	LOGICAL :: ld_nea ! Switch for rejecting observation near land
1603	!! * Local declarations
1604	INTEGER :: iyea0 ! Initial date
1605	INTEGER :: imon0 ! - (year, month, day, hour, minute)
1606	INTEGER :: iday0
1607	INTEGER :: ihou0
1608	INTEGER :: imin0
1609	INTEGER :: icycle ! Current assimilation cycle
1610	! Counters for observations that
1611	INTEGER :: iotdobs ! - outside time domain
1612	INTEGER :: iosdsobs ! - outside space domain
1613	INTEGER :: ilansobs ! - within a model land cell
1614	INTEGER :: inlasobs ! - close to land
1615	INTEGER :: igrdobs ! - fail the grid search
1616	INTEGER :: ibdysobs ! - close to open boundary
1617	! Global counters for observations that
1618	INTEGER :: iotdobsmpp ! - outside time domain
1619	INTEGER :: iosdsobsmpp ! - outside space domain
1620	INTEGER :: ilansobsmpp ! - within a model land cell
1621	INTEGER :: inlasobsmpp ! - close to land
1622	INTEGER :: igrdobsmpp ! - fail the grid search
1623	INTEGER :: ibdysobsmpp ! - close to open boundary
1624	LOGICAL, DIMENSION(:), ALLOCATABLE :: &
1625	& llvalid ! data selection
1626	INTEGER :: jobs ! Obs. loop variable
1627	INTEGER :: jstp ! Time loop variable
1628	INTEGER :: inrc ! Time index variable
1629	INTEGER :: irec ! Record index
1630
1631	IF (lwp) WRITE(numout,*)'obs_pre_spm : Preparing the spm observations...'
1632
1633	! Initial date initialization (year, month, day, hour, minute)
1634	iyea0 = ndate0 / 10000
1635	imon0 = ( ndate0 - iyea0 * 10000 ) / 100
1636	iday0 = ndate0 - iyea0 * 10000 - imon0 * 100
1637	ihou0 = 0
1638	imin0 = 0
1639
1640	icycle = no ! Assimilation cycle
1641
1642	! Diagnostics counters for various failures.
1643
1644	iotdobs = 0
1645	igrdobs = 0
1646	iosdsobs = 0
1647	ilansobs = 0
1648	inlasobs = 0
1649	ibdysobs = 0
1650
1651	! -----------------------------------------------------------------------
1652	! Find time coordinate for spm data
1653	! -----------------------------------------------------------------------
1654
1655	CALL obs_coo_tim( icycle, &
1656	& iyea0, imon0, iday0, ihou0, imin0, &
1657	& spmdata%nsurf, spmdata%nyea, spmdata%nmon, &
1658	& spmdata%nday, spmdata%nhou, spmdata%nmin, &
1659	& spmdata%nqc, spmdata%mstp, iotdobs )
1660	CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
1661	! -----------------------------------------------------------------------
1662	! Check for spm data failing the grid search
1663	! -----------------------------------------------------------------------
1664
1665	CALL obs_coo_grd( spmdata%nsurf, spmdata%mi, spmdata%mj, &
1666	& spmdata%nqc, igrdobs )
1667	CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
1668
1669	! -----------------------------------------------------------------------
1670	! Check for land points.
1671	! -----------------------------------------------------------------------
1672
1673	CALL obs_coo_spc_2d( spmdata%nsurf, &
1674	& jpi, jpj, &
1675	& spmdata%mi, spmdata%mj, &
1676	& spmdata%rlam, spmdata%rphi, &
1677	& glamt, gphit, &
1678	& tmask(:,:,1), spmdata%nqc, &
1679	& iosdsobs, ilansobs, &
1680	& inlasobs, ld_nea, &
1681	& ibdysobs, ln_bound_reject )
1682
1683	CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
1684	CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
1685	CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
1686	CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp )
1687
1688	! -----------------------------------------------------------------------
1689	! Copy useful data from the spmdata data structure to
1690	! the spmdatqc data structure
1691	! -----------------------------------------------------------------------
1692
1693	! Allocate the selection arrays
1694
1695	ALLOCATE( llvalid(spmdata%nsurf) )
1696
1697	! We want all data which has qc flags <= 0
1698
1699	llvalid(:) = ( spmdata%nqc(:) <= 10 )
1700
1701	! The actual copying
1702
1703	CALL obs_surf_compress( spmdata, spmdatqc, .TRUE., numout, &
1704	& lvalid=llvalid )
1705
1706	! Dellocate the selection arrays
1707	DEALLOCATE( llvalid )
1708
1709	! -----------------------------------------------------------------------
1710	! Print information about what observations are left after qc
1711	! -----------------------------------------------------------------------
1712
1713	! Update the total observation counter array
1714
1715	IF(lwp) THEN
1716	WRITE(numout,*)
1717	WRITE(numout,*) 'obs_pre_spm :'
1718	WRITE(numout,*) '~~~~~~~~~~~'
1719	WRITE(numout,*)
1720	WRITE(numout,*) ' spm data outside time domain = ', &
1721	& iotdobsmpp
1722	WRITE(numout,*) ' Remaining spm data that failed grid search = ', &
1723	& igrdobsmpp
1724	WRITE(numout,*) ' Remaining spm data outside space domain = ', &
1725	& iosdsobsmpp
1726	WRITE(numout,*) ' Remaining spm data at land points = ', &
1727	& ilansobsmpp
1728	IF (ld_nea) THEN
1729	WRITE(numout,*) ' Remaining spm data near land points (removed) = ', &
1730	& inlasobsmpp
1731	ELSE
1732	WRITE(numout,*) ' Remaining spm data near land points (kept) = ', &
1733	& inlasobsmpp
1734	ENDIF
1735	WRITE(numout,*) ' Remaining spm data near open boundary (removed) = ', &
1736	& ibdysobsmpp
1737	WRITE(numout,*) ' spm data accepted = ', &
1738	& spmdatqc%nsurfmpp
1739
1740	WRITE(numout,*)
1741	WRITE(numout,*) ' Number of observations per time step :'
1742	WRITE(numout,*)
1743	WRITE(numout,1997)
1744	WRITE(numout,1998)
1745	ENDIF
1746
1747	DO jobs = 1, spmdatqc%nsurf
1748	inrc = spmdatqc%mstp(jobs) + 2 - nit000
1749	spmdatqc%nsstp(inrc) = spmdatqc%nsstp(inrc) + 1
1750	END DO
1751
1752	CALL obs_mpp_sum_integers( spmdatqc%nsstp, spmdatqc%nsstpmpp, &
1753	& nitend - nit000 + 2 )
1754
1755	IF ( lwp ) THEN
1756	DO jstp = nit000 - 1, nitend
1757	inrc = jstp - nit000 + 2
1758	WRITE(numout,1999) jstp, spmdatqc%nsstpmpp(inrc)
1759	END DO
1760	ENDIF
1761
1762	1997 FORMAT(10X,'Time step',5X,'spm data')
1763	1998 FORMAT(10X,'---------',5X,'------------')
1764	1999 FORMAT(10X,I9,5X,I17)
1765
1766	END SUBROUTINE obs_pre_spm
1767
1768	SUBROUTINE obs_pre_fco2( fco2data, fco2datqc, ld_fco2, ld_nea )
1769	!!----------------------------------------------------------------------
1770	!! * ROUTINE obs_pre_fco2 *
1771	!!
1772	!! ** Purpose : First level check and screening of fco2 observations
1773	!!
1774	!! ** Method : First level check and screening of fco2 observations
1775	!!
1776	!! ** Action :
1777	!!
1778	!! References :
1779	!!
1780	!! History :
1781	!!----------------------------------------------------------------------
1782	!! * Modules used
1783	USE domstp ! Domain: set the time-step
1784	USE par_oce ! Ocean parameters
1785	USE dom_oce, ONLY : & ! Geographical information
1786	& glamt, &
1787	& gphit, &
1788	& tmask
1789	!! * Arguments
1790	TYPE(obs_surf), INTENT(INOUT) :: fco2data ! Full set of fco2 data
1791	TYPE(obs_surf), INTENT(INOUT) :: fco2datqc ! Subset of fco2 data not failing screening
1792	LOGICAL :: ld_fco2 ! Switch for fco2 data
1793	LOGICAL :: ld_nea ! Switch for rejecting observation near land
1794	!! * Local declarations
1795	INTEGER :: iyea0 ! Initial date
1796	INTEGER :: imon0 ! - (year, month, day, hour, minute)
1797	INTEGER :: iday0
1798	INTEGER :: ihou0
1799	INTEGER :: imin0
1800	INTEGER :: icycle ! Current assimilation cycle
1801	! Counters for observations that
1802	INTEGER :: iotdobs ! - outside time domain
1803	INTEGER :: iosdsobs ! - outside space domain
1804	INTEGER :: ilansobs ! - within a model land cell
1805	INTEGER :: inlasobs ! - close to land
1806	INTEGER :: igrdobs ! - fail the grid search
1807	INTEGER :: ibdysobs ! - close to open boundary
1808	! Global counters for observations that
1809	INTEGER :: iotdobsmpp ! - outside time domain
1810	INTEGER :: iosdsobsmpp ! - outside space domain
1811	INTEGER :: ilansobsmpp ! - within a model land cell
1812	INTEGER :: inlasobsmpp ! - close to land
1813	INTEGER :: igrdobsmpp ! - fail the grid search
1814	INTEGER :: ibdysobsmpp ! - close to open boundary
1815	LOGICAL, DIMENSION(:), ALLOCATABLE :: &
1816	& llvalid ! data selection
1817	INTEGER :: jobs ! Obs. loop variable
1818	INTEGER :: jstp ! Time loop variable
1819	INTEGER :: inrc ! Time index variable
1820	INTEGER :: irec ! Record index
1821
1822	IF (lwp) WRITE(numout,*)'obs_pre_fco2 : Preparing the fco2 observations...'
1823
1824	! Initial date initialization (year, month, day, hour, minute)
1825	iyea0 = ndate0 / 10000
1826	imon0 = ( ndate0 - iyea0 * 10000 ) / 100
1827	iday0 = ndate0 - iyea0 * 10000 - imon0 * 100
1828	ihou0 = 0
1829	imin0 = 0
1830
1831	icycle = no ! Assimilation cycle
1832
1833	! Diagnostics counters for various failures.
1834
1835	iotdobs = 0
1836	igrdobs = 0
1837	iosdsobs = 0
1838	ilansobs = 0
1839	inlasobs = 0
1840	ibdysobs = 0
1841
1842	! -----------------------------------------------------------------------
1843	! Find time coordinate for fco2 data
1844	! -----------------------------------------------------------------------
1845
1846	CALL obs_coo_tim( icycle, &
1847	& iyea0, imon0, iday0, ihou0, imin0, &
1848	& fco2data%nsurf, fco2data%nyea, fco2data%nmon, &
1849	& fco2data%nday, fco2data%nhou, fco2data%nmin, &
1850	& fco2data%nqc, fco2data%mstp, iotdobs )
1851	CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
1852	! -----------------------------------------------------------------------
1853	! Check for fco2 data failing the grid search
1854	! -----------------------------------------------------------------------
1855
1856	CALL obs_coo_grd( fco2data%nsurf, fco2data%mi, fco2data%mj, &
1857	& fco2data%nqc, igrdobs )
1858	CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
1859
1860	! -----------------------------------------------------------------------
1861	! Check for land points.
1862	! -----------------------------------------------------------------------
1863
1864	CALL obs_coo_spc_2d( fco2data%nsurf, &
1865	& jpi, jpj, &
1866	& fco2data%mi, fco2data%mj, &
1867	& fco2data%rlam, fco2data%rphi, &
1868	& glamt, gphit, &
1869	& tmask(:,:,1), fco2data%nqc, &
1870	& iosdsobs, ilansobs, &
1871	& inlasobs, ld_nea, &
1872	& ibdysobs, ln_bound_reject )
1873
1874	CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
1875	CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
1876	CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
1877	CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp )
1878
1879	! -----------------------------------------------------------------------
1880	! Copy useful data from the fco2data data structure to
1881	! the fco2datqc data structure
1882	! -----------------------------------------------------------------------
1883
1884	! Allocate the selection arrays
1885
1886	ALLOCATE( llvalid(fco2data%nsurf) )
1887
1888	! We want all data which has qc flags <= 0
1889
1890	llvalid(:) = ( fco2data%nqc(:) <= 10 )
1891
1892	! The actual copying
1893
1894	CALL obs_surf_compress( fco2data, fco2datqc, .TRUE., numout, &
1895	& lvalid=llvalid )
1896
1897	! Dellocate the selection arrays
1898	DEALLOCATE( llvalid )
1899
1900	! -----------------------------------------------------------------------
1901	! Print information about what observations are left after qc
1902	! -----------------------------------------------------------------------
1903
1904	! Update the total observation counter array
1905
1906	IF(lwp) THEN
1907	WRITE(numout,*)
1908	WRITE(numout,*) 'obs_pre_fco2 :'
1909	WRITE(numout,*) '~~~~~~~~~~~'
1910	WRITE(numout,*)
1911	WRITE(numout,*) ' fco2 data outside time domain = ', &
1912	& iotdobsmpp
1913	WRITE(numout,*) ' Remaining fco2 data that failed grid search = ', &
1914	& igrdobsmpp
1915	WRITE(numout,*) ' Remaining fco2 data outside space domain = ', &
1916	& iosdsobsmpp
1917	WRITE(numout,*) ' Remaining fco2 data at land points = ', &
1918	& ilansobsmpp
1919	IF (ld_nea) THEN
1920	WRITE(numout,*) ' Remaining fco2 data near land points (removed) = ', &
1921	& inlasobsmpp
1922	ELSE
1923	WRITE(numout,*) ' Remaining fco2 data near land points (kept) = ', &
1924	& inlasobsmpp
1925	ENDIF
1926	WRITE(numout,*) ' Remaining fco2 data near open boundary (removed) = ', &
1927	& ibdysobsmpp
1928	WRITE(numout,*) ' fco2 data accepted = ', &
1929	& fco2datqc%nsurfmpp
1930
1931	WRITE(numout,*)
1932	WRITE(numout,*) ' Number of observations per time step :'
1933	WRITE(numout,*)
1934	WRITE(numout,1997)
1935	WRITE(numout,1998)
1936	ENDIF
1937
1938	DO jobs = 1, fco2datqc%nsurf
1939	inrc = fco2datqc%mstp(jobs) + 2 - nit000
1940	fco2datqc%nsstp(inrc) = fco2datqc%nsstp(inrc) + 1
1941	END DO
1942
1943	CALL obs_mpp_sum_integers( fco2datqc%nsstp, fco2datqc%nsstpmpp, &
1944	& nitend - nit000 + 2 )
1945
1946	IF ( lwp ) THEN
1947	DO jstp = nit000 - 1, nitend
1948	inrc = jstp - nit000 + 2
1949	WRITE(numout,1999) jstp, fco2datqc%nsstpmpp(inrc)
1950	END DO
1951	ENDIF
1952
1953	1997 FORMAT(10X,'Time step',5X,'fco2 data')
1954	1998 FORMAT(10X,'---------',5X,'------------')
1955	1999 FORMAT(10X,I9,5X,I17)
1956
1957	END SUBROUTINE obs_pre_fco2
1958
1959	SUBROUTINE obs_pre_pco2( pco2data, pco2datqc, ld_pco2, ld_nea )
1960	!!----------------------------------------------------------------------
1961	!! * ROUTINE obs_pre_pco2 *
1962	!!
1963	!! ** Purpose : First level check and screening of pco2 observations
1964	!!
1965	!! ** Method : First level check and screening of pco2 observations
1966	!!
1967	!! ** Action :
1968	!!
1969	!! References :
1970	!!
1971	!! History :
1972	!!----------------------------------------------------------------------
1973	!! * Modules used
1974	USE domstp ! Domain: set the time-step
1975	USE par_oce ! Ocean parameters
1976	USE dom_oce, ONLY : & ! Geographical information
1977	& glamt, &
1978	& gphit, &
1979	& tmask
1980	!! * Arguments
1981	TYPE(obs_surf), INTENT(INOUT) :: pco2data ! Full set of pco2 data
1982	TYPE(obs_surf), INTENT(INOUT) :: pco2datqc ! Subset of pco2 data not failing screening
1983	LOGICAL :: ld_pco2 ! Switch for pco2 data
1984	LOGICAL :: ld_nea ! Switch for rejecting observation near land
1985	!! * Local declarations
1986	INTEGER :: iyea0 ! Initial date
1987	INTEGER :: imon0 ! - (year, month, day, hour, minute)
1988	INTEGER :: iday0
1989	INTEGER :: ihou0
1990	INTEGER :: imin0
1991	INTEGER :: icycle ! Current assimilation cycle
1992	! Counters for observations that
1993	INTEGER :: iotdobs ! - outside time domain
1994	INTEGER :: iosdsobs ! - outside space domain
1995	INTEGER :: ilansobs ! - within a model land cell
1996	INTEGER :: inlasobs ! - close to land
1997	INTEGER :: igrdobs ! - fail the grid search
1998	INTEGER :: ibdysobs ! - close to open boundary
1999	! Global counters for observations that
2000	INTEGER :: iotdobsmpp ! - outside time domain
2001	INTEGER :: iosdsobsmpp ! - outside space domain
2002	INTEGER :: ilansobsmpp ! - within a model land cell
2003	INTEGER :: inlasobsmpp ! - close to land
2004	INTEGER :: igrdobsmpp ! - fail the grid search
2005	INTEGER :: ibdysobsmpp ! - close to open boundary
2006	LOGICAL, DIMENSION(:), ALLOCATABLE :: &
2007	& llvalid ! data selection
2008	INTEGER :: jobs ! Obs. loop variable
2009	INTEGER :: jstp ! Time loop variable
2010	INTEGER :: inrc ! Time index variable
2011	INTEGER :: irec ! Record index
2012
2013	IF (lwp) WRITE(numout,*)'obs_pre_pco2 : Preparing the pco2 observations...'
2014
2015	! Initial date initialization (year, month, day, hour, minute)
2016	iyea0 = ndate0 / 10000
2017	imon0 = ( ndate0 - iyea0 * 10000 ) / 100
2018	iday0 = ndate0 - iyea0 * 10000 - imon0 * 100
2019	ihou0 = 0
2020	imin0 = 0
2021
2022	icycle = no ! Assimilation cycle
2023
2024	! Diagnostics counters for various failures.
2025
2026	iotdobs = 0
2027	igrdobs = 0
2028	iosdsobs = 0
2029	ilansobs = 0
2030	inlasobs = 0
2031	ibdysobs = 0
2032
2033	! -----------------------------------------------------------------------
2034	! Find time coordinate for pco2 data
2035	! -----------------------------------------------------------------------
2036
2037	CALL obs_coo_tim( icycle, &
2038	& iyea0, imon0, iday0, ihou0, imin0, &
2039	& pco2data%nsurf, pco2data%nyea, pco2data%nmon, &
2040	& pco2data%nday, pco2data%nhou, pco2data%nmin, &
2041	& pco2data%nqc, pco2data%mstp, iotdobs )
2042	CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
2043	! -----------------------------------------------------------------------
2044	! Check for pco2 data failing the grid search
2045	! -----------------------------------------------------------------------
2046
2047	CALL obs_coo_grd( pco2data%nsurf, pco2data%mi, pco2data%mj, &
2048	& pco2data%nqc, igrdobs )
2049	CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
2050
2051	! -----------------------------------------------------------------------
2052	! Check for land points.
2053	! -----------------------------------------------------------------------
2054
2055	CALL obs_coo_spc_2d( pco2data%nsurf, &
2056	& jpi, jpj, &
2057	& pco2data%mi, pco2data%mj, &
2058	& pco2data%rlam, pco2data%rphi, &
2059	& glamt, gphit, &
2060	& tmask(:,:,1), pco2data%nqc, &
2061	& iosdsobs, ilansobs, &
2062	& inlasobs, ld_nea, &
2063	& ibdysobs, ln_bound_reject )
2064
2065	CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
2066	CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
2067	CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
2068	CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp )
2069
2070	! -----------------------------------------------------------------------
2071	! Copy useful data from the pco2data data structure to
2072	! the pco2datqc data structure
2073	! -----------------------------------------------------------------------
2074
2075	! Allocate the selection arrays
2076
2077	ALLOCATE( llvalid(pco2data%nsurf) )
2078
2079	! We want all data which has qc flags <= 0
2080
2081	llvalid(:) = ( pco2data%nqc(:) <= 10 )
2082
2083	! The actual copying
2084
2085	CALL obs_surf_compress( pco2data, pco2datqc, .TRUE., numout, &
2086	& lvalid=llvalid )
2087
2088	! Dellocate the selection arrays
2089	DEALLOCATE( llvalid )
2090
2091	! -----------------------------------------------------------------------
2092	! Print information about what observations are left after qc
2093	! -----------------------------------------------------------------------
2094
2095	! Update the total observation counter array
2096
2097	IF(lwp) THEN
2098	WRITE(numout,*)
2099	WRITE(numout,*) 'obs_pre_pco2 :'
2100	WRITE(numout,*) '~~~~~~~~~~~'
2101	WRITE(numout,*)
2102	WRITE(numout,*) ' pco2 data outside time domain = ', &
2103	& iotdobsmpp
2104	WRITE(numout,*) ' Remaining pco2 data that failed grid search = ', &
2105	& igrdobsmpp
2106	WRITE(numout,*) ' Remaining pco2 data outside space domain = ', &
2107	& iosdsobsmpp
2108	WRITE(numout,*) ' Remaining pco2 data at land points = ', &
2109	& ilansobsmpp
2110	IF (ld_nea) THEN
2111	WRITE(numout,*) ' Remaining pco2 data near land points (removed) = ', &
2112	& inlasobsmpp
2113	ELSE
2114	WRITE(numout,*) ' Remaining pco2 data near land points (kept) = ', &
2115	& inlasobsmpp
2116	ENDIF
2117	WRITE(numout,*) ' Remaining pco2 data near open boundary (removed) = ', &
2118	& ibdysobsmpp
2119	WRITE(numout,*) ' pco2 data accepted = ', &
2120	& pco2datqc%nsurfmpp
2121
2122	WRITE(numout,*)
2123	WRITE(numout,*) ' Number of observations per time step :'
2124	WRITE(numout,*)
2125	WRITE(numout,1997)
2126	WRITE(numout,1998)
2127	ENDIF
2128
2129	DO jobs = 1, pco2datqc%nsurf
2130	inrc = pco2datqc%mstp(jobs) + 2 - nit000
2131	pco2datqc%nsstp(inrc) = pco2datqc%nsstp(inrc) + 1
2132	END DO
2133
2134	CALL obs_mpp_sum_integers( pco2datqc%nsstp, pco2datqc%nsstpmpp, &
2135	& nitend - nit000 + 2 )
2136
2137	IF ( lwp ) THEN
2138	DO jstp = nit000 - 1, nitend
2139	inrc = jstp - nit000 + 2
2140	WRITE(numout,1999) jstp, pco2datqc%nsstpmpp(inrc)
2141	END DO
2142	ENDIF
2143
2144	1997 FORMAT(10X,'Time step',5X,'pco2 data')
2145	1998 FORMAT(10X,'---------',5X,'------------')
2146	1999 FORMAT(10X,I9,5X,I17)
2147
2148	END SUBROUTINE obs_pre_pco2
2149
2150	SUBROUTINE obs_coo_tim( kcycle, &
2151	& kyea0, kmon0, kday0, khou0, kmin0, &
2152	& kobsno, &
2153	& kobsyea, kobsmon, kobsday, kobshou, kobsmin, &
2154	& kobsqc, kobsstp, kotdobs )
2155	!!----------------------------------------------------------------------
2156	!! * ROUTINE obs_coo_tim *
2157	!!
2158	!! ** Purpose : Compute the number of time steps to the observation time.
2159	!!
2160	!! ** Method : For time coordinates ( yea_obs, mon_obs, day_obs,
2161	!! hou_obs, min_obs ), this routine locates the time step
2162	!! that is closest to this time.
2163	!!
2164	!! ** Action :
2165	!!
2166	!! References :
2167	!!
2168	!! History :
2169	!! ! 1997-07 (A. Weaver) Original
2170	!! ! 2006-08 (A. Weaver) NEMOVAR migration
2171	!! ! 2006-10 (A. Weaver) Cleanup
2172	!! ! 2007-01 (K. Mogensen) Rewritten with loop
2173	!! ! 2010-05 (D. Lea) Fix in leap year calculation for NEMO vn3.2
2174	!!----------------------------------------------------------------------
2175	!! * Modules used
2176	USE dom_oce, ONLY : & ! Geographical information
2177	& rdt
2178	USE phycst, ONLY : & ! Physical constants
2179	& rday, &
2180	& rmmss, &
2181	& rhhmm
2182	!! * Arguments
2183	INTEGER, INTENT(IN) :: kcycle ! Current cycle
2184	INTEGER, INTENT(IN) :: kyea0 ! Initial date coordinates
2185	INTEGER, INTENT(IN) :: kmon0
2186	INTEGER, INTENT(IN) :: kday0
2187	INTEGER, INTENT(IN) :: khou0
2188	INTEGER, INTENT(IN) :: kmin0
2189	INTEGER, INTENT(IN) :: kobsno ! Number of observations
2190	INTEGER, INTENT(INOUT) :: kotdobs ! Number of observations failing time check
2191	INTEGER, DIMENSION(kobsno), INTENT(IN ) :: &
2192	& kobsyea, & ! Observation time coordinates
2193	& kobsmon, &
2194	& kobsday, &
2195	& kobshou, &
2196	& kobsmin
2197	INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
2198	& kobsqc ! Quality control flag
2199	INTEGER, DIMENSION(kobsno), INTENT(OUT) :: &
2200	& kobsstp ! Number of time steps up to the
2201	! observation time
2202
2203	!! * Local declarations
2204	INTEGER :: jyea
2205	INTEGER :: jmon
2206	INTEGER :: jday
2207	INTEGER :: jobs
2208	INTEGER :: iyeastr
2209	INTEGER :: iyeaend
2210	INTEGER :: imonstr
2211	INTEGER :: imonend
2212	INTEGER :: idaystr
2213	INTEGER :: idayend
2214	INTEGER :: iskip
2215	INTEGER :: idaystp
2216	REAL(KIND=wp) :: zminstp
2217	REAL(KIND=wp) :: zhoustp
2218	REAL(KIND=wp) :: zobsstp
2219	INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year
2220
2221	!-----------------------------------------------------------------------
2222	! Initialization
2223	!-----------------------------------------------------------------------
2224
2225	! Intialize the number of time steps per day
2226	idaystp = NINT( rday / rdt )
2227
2228	!---------------------------------------------------------------------
2229	! Locate the model time coordinates for interpolation
2230	!---------------------------------------------------------------------
2231
2232	DO jobs = 1, kobsno
2233
2234	! Initialize the time step counter
2235	kobsstp(jobs) = nit000 - 1
2236
2237	! Flag if observation date is less than the initial date
2238
2239	IF ( ( kobsyea(jobs) < kyea0 ) &
2240	& .OR. ( ( kobsyea(jobs) == kyea0 ) &
2241	& .AND. ( kobsmon(jobs) < kmon0 ) ) &
2242	& .OR. ( ( kobsyea(jobs) == kyea0 ) &
2243	& .AND. ( kobsmon(jobs) == kmon0 ) &
2244	& .AND. ( kobsday(jobs) < kday0 ) ) &
2245	& .OR. ( ( kobsyea(jobs) == kyea0 ) &
2246	& .AND. ( kobsmon(jobs) == kmon0 ) &
2247	& .AND. ( kobsday(jobs) == kday0 ) &
2248	& .AND. ( kobshou(jobs) < khou0 ) ) &
2249	& .OR. ( ( kobsyea(jobs) == kyea0 ) &
2250	& .AND. ( kobsmon(jobs) == kmon0 ) &
2251	& .AND. ( kobsday(jobs) == kday0 ) &
2252	& .AND. ( kobshou(jobs) == khou0 ) &
2253	& .AND. ( kobsmin(jobs) <= kmin0 ) ) ) THEN
2254	kobsstp(jobs) = -1
2255	kobsqc(jobs) = kobsqc(jobs) + 11
2256	kotdobs = kotdobs + 1
2257	CYCLE
2258	ENDIF
2259
2260	! Compute the number of time steps to the observation day
2261	iyeastr = kyea0
2262	iyeaend = kobsyea(jobs)
2263
2264	!---------------------------------------------------------------------
2265	! Year loop
2266	!---------------------------------------------------------------------
2267	DO jyea = iyeastr, iyeaend
2268
2269	CALL calc_month_len( jyea, imonth_len )
2270
2271	imonstr = 1
2272	IF ( jyea == kyea0 ) imonstr = kmon0
2273	imonend = 12
2274	IF ( jyea == kobsyea(jobs) ) imonend = kobsmon(jobs)
2275
2276	! Month loop
2277	DO jmon = imonstr, imonend
2278
2279	idaystr = 1
2280	IF ( ( jmon == kmon0 ) &
2281	& .AND. ( jyea == kyea0 ) ) idaystr = kday0
2282	idayend = imonth_len(jmon)
2283	IF ( ( jmon == kobsmon(jobs) ) &
2284	& .AND. ( jyea == kobsyea(jobs) ) ) idayend = kobsday(jobs) - 1
2285
2286	! Day loop
2287	DO jday = idaystr, idayend
2288	kobsstp(jobs) = kobsstp(jobs) + idaystp
2289	END DO
2290
2291	END DO
2292
2293	END DO
2294
2295	! Add in the number of time steps to the observation minute
2296	zminstp = rmmss / rdt
2297	zhoustp = rhhmm * zminstp
2298
2299	zobsstp = REAL( kobsmin(jobs) - kmin0, KIND=wp ) * zminstp &
2300	& + REAL( kobshou(jobs) - khou0, KIND=wp ) * zhoustp
2301	kobsstp(jobs) = kobsstp(jobs) + NINT( zobsstp )
2302
2303	! Flag if observation step outside the time window
2304	IF ( ( kobsstp(jobs) < ( nit000 - 1 ) ) &
2305	& .OR.( kobsstp(jobs) > nitend ) ) THEN
2306	kobsqc(jobs) = kobsqc(jobs) + 12
2307	kotdobs = kotdobs + 1
2308	CYCLE
2309	ENDIF
2310
2311	END DO
2312
2313	END SUBROUTINE obs_coo_tim
2314
2315	SUBROUTINE calc_month_len( iyear, imonth_len )
2316	!!----------------------------------------------------------------------
2317	!! * ROUTINE calc_month_len *
2318	!!
2319	!! ** Purpose : Compute the number of days in a months given a year.
2320	!!
2321	!! ** Method :
2322	!!
2323	!! ** Action :
2324	!!
2325	!! History :
2326	!! ! 10-05 (D. Lea) New routine based on day_init
2327	!!----------------------------------------------------------------------
2328
2329	INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year
2330	INTEGER :: iyear !: year
2331
2332	! length of the month of the current year (from nleapy, read in namelist)
2333	IF ( nleapy < 2 ) THEN
2334	imonth_len(:) = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)
2335	IF ( nleapy == 1 ) THEN ! we are using calendar with leap years
2336	IF ( MOD(iyear, 4) == 0 .AND. ( MOD(iyear, 400) == 0 .OR. MOD(iyear, 100) /= 0 ) ) THEN
2337	imonth_len(2) = 29
2338	ENDIF
2339	ENDIF
2340	ELSE
2341	imonth_len(:) = nleapy ! all months with nleapy days per year
2342	ENDIF
2343
2344	END SUBROUTINE
2345
2346	SUBROUTINE obs_coo_tim_prof( kcycle, &
2347	& kyea0, kmon0, kday0, khou0, kmin0, &
2348	& kobsno, &
2349	& kobsyea, kobsmon, kobsday, kobshou, kobsmin, &
2350	& ktyp, kobsqc, kobsstp, kotdobs, kdailyavtypes, &
2351	& ld_dailyav )
2352	!!----------------------------------------------------------------------
2353	!! * ROUTINE obs_coo_tim *
2354	!!
2355	!! ** Purpose : Compute the number of time steps to the observation time.
2356	!!
2357	!! ** Method : For time coordinates ( yea_obs, mon_obs, day_obs,
2358	!! hou_obs, min_obs ), this routine locates the time step
2359	!! that is closest to this time.
2360	!!
2361	!! ** Action :
2362	!!
2363	!! References :
2364	!!
2365	!! History :
2366	!! ! 1997-07 (A. Weaver) Original
2367	!! ! 2006-08 (A. Weaver) NEMOVAR migration
2368	!! ! 2006-10 (A. Weaver) Cleanup
2369	!! ! 2007-01 (K. Mogensen) Rewritten with loop
2370	!!----------------------------------------------------------------------
2371	!! * Modules used
2372	!! * Arguments
2373	INTEGER, INTENT(IN) :: kcycle ! Current cycle
2374	INTEGER, INTENT(IN) :: kyea0 ! Initial date coordinates
2375	INTEGER, INTENT(IN) :: kmon0
2376	INTEGER, INTENT(IN) :: kday0
2377	INTEGER, INTENT(IN) :: khou0
2378	INTEGER, INTENT(IN) :: kmin0
2379	INTEGER, INTENT(IN) :: kobsno ! Number of observations
2380	INTEGER, INTENT(INOUT) :: kotdobs ! Number of observations failing time check
2381	INTEGER, DIMENSION(kobsno), INTENT(IN ) :: &
2382	& kobsyea, & ! Observation time coordinates
2383	& kobsmon, &
2384	& kobsday, &
2385	& kobshou, &
2386	& kobsmin, &
2387	& ktyp ! Observation type.
2388	INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
2389	& kobsqc ! Quality control flag
2390	INTEGER, DIMENSION(kobsno), INTENT(OUT) :: &
2391	& kobsstp ! Number of time steps up to the
2392	! observation time
2393	INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
2394	& kdailyavtypes ! Types for daily averages
2395	LOGICAL, OPTIONAL :: ld_dailyav ! All types are daily averages
2396	!! * Local declarations
2397	INTEGER :: jobs
2398
2399	!-----------------------------------------------------------------------
2400	! Call standard obs_coo_tim
2401	!-----------------------------------------------------------------------
2402
2403	CALL obs_coo_tim( kcycle, &
2404	& kyea0, kmon0, kday0, khou0, kmin0, &
2405	& kobsno, &
2406	& kobsyea, kobsmon, kobsday, kobshou, kobsmin, &
2407	& kobsqc, kobsstp, kotdobs )
2408
2409	!------------------------------------------------------------------------
2410	! Always reject daily averaged data (e.g. MRB data (820)) at initial time
2411	!------------------------------------------------------------------------
2412
2413	IF ( PRESENT(kdailyavtypes) ) THEN
2414	DO jobs = 1, kobsno
2415
2416	IF ( kobsqc(jobs) <= 10 ) THEN
2417
2418	IF ( ( kobsstp(jobs) == (nit000 - 1) ).AND.&
2419	& ( ANY (kdailyavtypes(:) == ktyp(jobs)) ) ) THEN
2420	kobsqc(jobs) = kobsqc(jobs) + 14
2421	kotdobs = kotdobs + 1
2422	CYCLE
2423	ENDIF
2424
2425	ENDIF
2426	END DO
2427	ENDIF
2428
2429	!------------------------------------------------------------------------
2430	! If ld_dailyav is set then all data assumed to be daily averaged
2431	!------------------------------------------------------------------------
2432
2433	IF ( PRESENT( ld_dailyav) ) THEN
2434	IF (ld_dailyav) THEN
2435	DO jobs = 1, kobsno
2436
2437	IF ( kobsqc(jobs) <= 10 ) THEN
2438
2439	IF ( kobsstp(jobs) == (nit000 - 1) ) THEN
2440	kobsqc(jobs) = kobsqc(jobs) + 14
2441	kotdobs = kotdobs + 1
2442	CYCLE
2443	ENDIF
2444
2445	ENDIF
2446	END DO
2447	ENDIF
2448	ENDIF
2449
2450	END SUBROUTINE obs_coo_tim_prof
2451
2452	SUBROUTINE obs_coo_grd( kobsno, kobsi, kobsj, kobsqc, kgrdobs )
2453	!!----------------------------------------------------------------------
2454	!! * ROUTINE obs_coo_grd *
2455	!!
2456	!! ** Purpose : Verify that the grid search has not failed
2457	!!
2458	!! ** Method : The previously computed i,j indeces are checked
2459	!!
2460	!! ** Action :
2461	!!
2462	!! References :
2463	!!
2464	!! History :
2465	!! ! 2007-01 (K. Mogensen) Original
2466	!!----------------------------------------------------------------------
2467	!! * Arguments
2468	INTEGER, INTENT(IN) :: kobsno ! Number of observations
2469	INTEGER, DIMENSION(kobsno), INTENT(IN ) :: &
2470	& kobsi, & ! i,j indeces previously computed
2471	& kobsj
2472	INTEGER, INTENT(INOUT) :: kgrdobs ! Number of observations failing the check
2473	INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
2474	& kobsqc ! Quality control flag
2475
2476	!! * Local declarations
2477	INTEGER :: jobs ! Loop variable
2478
2479	! Flag if the grid search failed
2480
2481	DO jobs = 1, kobsno
2482	IF ( ( kobsi(jobs) <= 0 ) .AND. ( kobsj(jobs) <= 0 ) ) THEN
2483	kobsqc(jobs) = kobsqc(jobs) + 18
2484	kgrdobs = kgrdobs + 1
2485	ENDIF
2486	END DO
2487
2488	END SUBROUTINE obs_coo_grd
2489
2490	SUBROUTINE obs_coo_spc_2d( kobsno, kpi, kpj, &
2491	& kobsi, kobsj, pobslam, pobsphi, &
2492	& plam, pphi, pmask, &
2493	& kobsqc, kosdobs, klanobs, &
2494	& knlaobs,ld_nea, &
2495	& kbdyobs,ld_bound_reject )
2496	!!----------------------------------------------------------------------
2497	!! * ROUTINE obs_coo_spc_2d *
2498	!!
2499	!! ** Purpose : Check for points outside the domain and land points
2500	!!
2501	!! ** Method : Remove the observations that are outside the model space
2502	!! and time domain or located within model land cells.
2503	!!
2504	!! ** Action :
2505	!!
2506	!! History :
2507	!! ! 2007-03 (A. Weaver, K. Mogensen) Original
2508	!! ! 2007-06 (K. Mogensen et al) Reject obs. near land.
2509	!!----------------------------------------------------------------------
2510	!! * Modules used
2511
2512	!! * Arguments
2513	INTEGER, INTENT(IN) :: kobsno ! Total number of observations
2514	INTEGER, INTENT(IN) :: kpi ! Number of grid points in (i,j)
2515	INTEGER, INTENT(IN) :: kpj
2516	INTEGER, DIMENSION(kobsno), INTENT(IN) :: &
2517	& kobsi, & ! Observation (i,j) coordinates
2518	& kobsj
2519	REAL(KIND=wp), DIMENSION(kobsno), INTENT(IN) :: &
2520	& pobslam, & ! Observation (lon,lat) coordinates
2521	& pobsphi
2522	REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) :: &
2523	& plam, pphi ! Model (lon,lat) coordinates
2524	REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) :: &
2525	& pmask ! Land mask array
2526	INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
2527	& kobsqc ! Observation quality control
2528	INTEGER, INTENT(INOUT) :: kosdobs ! Observations outside space domain
2529	INTEGER, INTENT(INOUT) :: klanobs ! Observations within a model land cell
2530	INTEGER, INTENT(INOUT) :: knlaobs ! Observations near land
2531	INTEGER, INTENT(INOUT) :: kbdyobs ! Observations near boundary
2532	LOGICAL, INTENT(IN) :: ld_nea ! Flag observations near land
2533	LOGICAL, INTENT(IN) :: ld_bound_reject ! Flag observations near open boundary
2534	!! * Local declarations
2535	REAL(KIND=wp), DIMENSION(2,2,kobsno) :: &
2536	& zgmsk ! Grid mask
2537	#if defined key_bdy
2538	REAL(KIND=wp), DIMENSION(2,2,kobsno) :: &
2539	& zbmsk ! Boundary mask
2540	REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask
2541	#endif
2542	REAL(KIND=wp), DIMENSION(2,2,kobsno) :: &
2543	& zglam, & ! Model longitude at grid points
2544	& zgphi ! Model latitude at grid points
2545	INTEGER, DIMENSION(2,2,kobsno) :: &
2546	& igrdi, & ! Grid i,j
2547	& igrdj
2548	LOGICAL :: lgridobs ! Is observation on a model grid point.
2549	INTEGER :: iig, ijg ! i,j of observation on model grid point.
2550	INTEGER :: jobs, ji, jj
2551
2552	! Get grid point indices
2553
2554	DO jobs = 1, kobsno
2555
2556	! For invalid points use 2,2
2557
2558	IF ( kobsqc(jobs) >= 10 ) THEN
2559
2560	igrdi(1,1,jobs) = 1
2561	igrdj(1,1,jobs) = 1
2562	igrdi(1,2,jobs) = 1
2563	igrdj(1,2,jobs) = 2
2564	igrdi(2,1,jobs) = 2
2565	igrdj(2,1,jobs) = 1
2566	igrdi(2,2,jobs) = 2
2567	igrdj(2,2,jobs) = 2
2568
2569	ELSE
2570
2571	igrdi(1,1,jobs) = kobsi(jobs)-1
2572	igrdj(1,1,jobs) = kobsj(jobs)-1
2573	igrdi(1,2,jobs) = kobsi(jobs)-1
2574	igrdj(1,2,jobs) = kobsj(jobs)
2575	igrdi(2,1,jobs) = kobsi(jobs)
2576	igrdj(2,1,jobs) = kobsj(jobs)-1
2577	igrdi(2,2,jobs) = kobsi(jobs)
2578	igrdj(2,2,jobs) = kobsj(jobs)
2579
2580	ENDIF
2581
2582	END DO
2583
2584	#if defined key_bdy
2585	! Create a mask grid points in boundary rim
2586	IF (ld_bound_reject) THEN
2587	zbdymask(:,:) = 1.0_wp
2588	DO ji = 1, nb_bdy
2589	DO jj = 1, idx_bdy(ji)%nblen(1)
2590	zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp
2591	ENDDO
2592	ENDDO
2593
2594	CALL obs_int_comm_2d( 2, 2, kobsno, igrdi, igrdj, zbdymask, zbmsk )
2595	ENDIF
2596	#endif
2597
2598	CALL obs_int_comm_2d( 2, 2, kobsno, igrdi, igrdj, pmask, zgmsk )
2599	CALL obs_int_comm_2d( 2, 2, kobsno, igrdi, igrdj, plam, zglam )
2600	CALL obs_int_comm_2d( 2, 2, kobsno, igrdi, igrdj, pphi, zgphi )
2601
2602	DO jobs = 1, kobsno
2603
2604	! Skip bad observations
2605	IF ( kobsqc(jobs) >= 10 ) CYCLE
2606
2607	! Flag if the observation falls outside the model spatial domain
2608	IF ( ( pobslam(jobs) < -180. ) &
2609	& .OR. ( pobslam(jobs) > 180. ) &
2610	& .OR. ( pobsphi(jobs) < -90. ) &
2611	& .OR. ( pobsphi(jobs) > 90. ) ) THEN
2612	kobsqc(jobs) = kobsqc(jobs) + 11
2613	kosdobs = kosdobs + 1
2614	CYCLE
2615	ENDIF
2616
2617	! Flag if the observation falls with a model land cell
2618	IF ( SUM( zgmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN
2619	kobsqc(jobs) = kobsqc(jobs) + 12
2620	klanobs = klanobs + 1
2621	CYCLE
2622	ENDIF
2623
2624	! Check if this observation is on a grid point
2625
2626	lgridobs = .FALSE.
2627	iig = -1
2628	ijg = -1
2629	DO jj = 1, 2
2630	DO ji = 1, 2
2631	IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) &
2632	& .AND. &
2633	& ( ABS( zglam(ji,jj,jobs) - pobslam(jobs) ) < 1.0e-6_wp ) &
2634	& ) THEN
2635	lgridobs = .TRUE.
2636	iig = ji
2637	ijg = jj
2638	ENDIF
2639	END DO
2640	END DO
2641
2642	! Flag if the observation falls is close to land
2643	IF ( MINVAL( zgmsk(1:2,1:2,jobs) ) == 0.0_wp) THEN
2644	knlaobs = knlaobs + 1
2645	IF (ld_nea) THEN
2646	kobsqc(jobs) = kobsqc(jobs) + 14
2647	CYCLE
2648	ENDIF
2649	ENDIF
2650
2651	#if defined key_bdy
2652	! Flag if the observation falls close to the boundary rim
2653	IF (ld_bound_reject) THEN
2654	IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN
2655	kobsqc(jobs) = kobsqc(jobs) + 15
2656	kbdyobs = kbdyobs + 1
2657	CYCLE
2658	ENDIF
2659	! for observations on the grid...
2660	IF (lgridobs) THEN
2661	IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN
2662	kobsqc(jobs) = kobsqc(jobs) + 15
2663	kbdyobs = kbdyobs + 1
2664	CYCLE
2665	ENDIF
2666	ENDIF
2667	ENDIF
2668	#endif
2669
2670	END DO
2671
2672	END SUBROUTINE obs_coo_spc_2d
2673
2674	SUBROUTINE obs_coo_spc_3d( kprofno, kobsno, kpstart, kpend, &
2675	& kpi, kpj, kpk, &
2676	& kobsi, kobsj, kobsk, &
2677	& pobslam, pobsphi, pobsdep, &
2678	& plam, pphi, pdep, pmask, &
2679	& kpobsqc, kobsqc, kosdobs, &
2680	& klanobs, knlaobs, ld_nea, &
2681	& kbdyobs, ld_bound_reject )
2682	!!----------------------------------------------------------------------
2683	!! * ROUTINE obs_coo_spc_3d *
2684	!!
2685	!! ** Purpose : Check for points outside the domain and land points
2686	!! Reset depth of observation above highest model level
2687	!! to the value of highest model level
2688	!!
2689	!! ** Method : Remove the observations that are outside the model space
2690	!! and time domain or located within model land cells.
2691	!!
2692	!! NB. T and S profile observations lying between the ocean
2693	!! surface and the depth of the first model T point are
2694	!! assigned a depth equal to that of the first model T pt.
2695	!!
2696	!! ** Action :
2697	!!
2698	!! History :
2699	!! ! 2007-01 (K. Mogensen) Rewrite of parts of obs_scr
2700	!! ! 2007-06 (K. Mogensen et al) Reject obs. near land.
2701	!!----------------------------------------------------------------------
2702	!! * Modules used
2703	USE dom_oce, ONLY : & ! Geographical information
2704	& gdepw_1d, &
2705	& gdepw_0, &
2706	#if defined key_vvl
2707	& gdepw_n, &
2708	& gdept_n, &
2709	#endif
2710	& ln_zco, &
2711	& ln_zps, &
2712	& lk_vvl
2713
2714	!! * Arguments
2715	INTEGER, INTENT(IN) :: kprofno ! Number of profiles
2716	INTEGER, INTENT(IN) :: kobsno ! Total number of observations
2717	INTEGER, INTENT(IN) :: kpi ! Number of grid points in (i,j,k)
2718	INTEGER, INTENT(IN) :: kpj
2719	INTEGER, INTENT(IN) :: kpk
2720	INTEGER, DIMENSION(kprofno), INTENT(IN) :: &
2721	& kpstart, & ! Start of individual profiles
2722	& kpend ! End of individual profiles
2723	INTEGER, DIMENSION(kprofno), INTENT(IN) :: &
2724	& kobsi, & ! Observation (i,j) coordinates
2725	& kobsj
2726	INTEGER, DIMENSION(kobsno), INTENT(IN) :: &
2727	& kobsk ! Observation k coordinate
2728	REAL(KIND=wp), DIMENSION(kprofno), INTENT(IN) :: &
2729	& pobslam, & ! Observation (lon,lat) coordinates
2730	& pobsphi
2731	REAL(KIND=wp), DIMENSION(kobsno), INTENT(INOUT) :: &
2732	& pobsdep ! Observation depths
2733	REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) :: &
2734	& plam, pphi ! Model (lon,lat) coordinates
2735	REAL(KIND=wp), DIMENSION(kpk), INTENT(IN) :: &
2736	& pdep ! Model depth coordinates
2737	REAL(KIND=wp), DIMENSION(kpi,kpj,kpk), INTENT(IN) :: &
2738	& pmask ! Land mask array
2739	INTEGER, DIMENSION(kprofno), INTENT(INOUT) :: &
2740	& kpobsqc ! Profile quality control
2741	INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
2742	& kobsqc ! Observation quality control
2743	INTEGER, INTENT(INOUT) :: kosdobs ! Observations outside space domain
2744	INTEGER, INTENT(INOUT) :: klanobs ! Observations within a model land cell
2745	INTEGER, INTENT(INOUT) :: knlaobs ! Observations near land
2746	INTEGER, INTENT(INOUT) :: kbdyobs ! Observations near boundary
2747	LOGICAL, INTENT(IN) :: ld_nea ! Flag observations near land
2748	LOGICAL, INTENT(IN) :: ld_bound_reject ! Flag observations near open boundary
2749	!! * Local declarations
2750	REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: &
2751	& zgmsk ! Grid mask
2752	#if defined key_bdy
2753	REAL(KIND=wp), DIMENSION(2,2,kprofno) :: &
2754	& zbmsk ! Boundary mask
2755	REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask
2756	#endif
2757	REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: &
2758	& zgdepw
2759	REAL(KIND=wp), DIMENSION(2,2,kprofno) :: &
2760	& zglam, & ! Model longitude at grid points
2761	& zgphi ! Model latitude at grid points
2762	INTEGER, DIMENSION(2,2,kprofno) :: &
2763	& igrdi, & ! Grid i,j
2764	& igrdj
2765	LOGICAL :: lgridobs ! Is observation on a model grid point.
2766	LOGICAL :: ll_next_to_land ! Is a profile next to land
2767	INTEGER :: iig, ijg ! i,j of observation on model grid point.
2768	INTEGER :: jobs, jobsp, jk, ji, jj
2769
2770	! Get grid point indices
2771
2772	DO jobs = 1, kprofno
2773
2774	! For invalid points use 2,2
2775
2776	IF ( kpobsqc(jobs) >= 10 ) THEN
2777
2778	igrdi(1,1,jobs) = 1
2779	igrdj(1,1,jobs) = 1
2780	igrdi(1,2,jobs) = 1
2781	igrdj(1,2,jobs) = 2
2782	igrdi(2,1,jobs) = 2
2783	igrdj(2,1,jobs) = 1
2784	igrdi(2,2,jobs) = 2
2785	igrdj(2,2,jobs) = 2
2786
2787	ELSE
2788
2789	igrdi(1,1,jobs) = kobsi(jobs)-1
2790	igrdj(1,1,jobs) = kobsj(jobs)-1
2791	igrdi(1,2,jobs) = kobsi(jobs)-1
2792	igrdj(1,2,jobs) = kobsj(jobs)
2793	igrdi(2,1,jobs) = kobsi(jobs)
2794	igrdj(2,1,jobs) = kobsj(jobs)-1
2795	igrdi(2,2,jobs) = kobsi(jobs)
2796	igrdj(2,2,jobs) = kobsj(jobs)
2797
2798	ENDIF
2799
2800	END DO
2801
2802	#if defined key_bdy
2803	! Create a mask grid points in boundary rim
2804	IF (ld_bound_reject) THEN
2805	zbdymask(:,:) = 1.0_wp
2806	DO ji = 1, nb_bdy
2807	DO jj = 1, idx_bdy(ji)%nblen(1)
2808	zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp
2809	ENDDO
2810	ENDDO
2811	ENDIF
2812
2813	CALL obs_int_comm_2d( 2, 2, kprofno, igrdi, igrdj, zbdymask, zbmsk )
2814	#endif
2815
2816	CALL obs_int_comm_3d( 2, 2, kprofno, kpk, igrdi, igrdj, pmask, zgmsk )
2817	CALL obs_int_comm_2d( 2, 2, kprofno, igrdi, igrdj, plam, zglam )
2818	CALL obs_int_comm_2d( 2, 2, kprofno, igrdi, igrdj, pphi, zgphi )
2819	! Need to know the bathy depth for each observation for sco
2820	CALL obs_int_comm_3d( 2, 2, kprofno, kpk, igrdi, igrdj, fsdepw(:,:,:), &
2821	& zgdepw )
2822
2823	DO jobs = 1, kprofno
2824
2825	! Skip bad profiles
2826	IF ( kpobsqc(jobs) >= 10 ) CYCLE
2827
2828	! Check if this observation is on a grid point
2829
2830	lgridobs = .FALSE.
2831	iig = -1
2832	ijg = -1
2833	DO jj = 1, 2
2834	DO ji = 1, 2
2835	IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) &
2836	& .AND. &
2837	& ( ABS( zglam(ji,jj,jobs) - pobslam(jobs) ) < 1.0e-6_wp ) &
2838	& ) THEN
2839	lgridobs = .TRUE.
2840	iig = ji
2841	ijg = jj
2842	ENDIF
2843	END DO
2844	END DO
2845
2846	! Check if next to land
2847	IF ( ANY( zgmsk(1:2,1:2,1,jobs) == 0.0_wp ) ) THEN
2848	ll_next_to_land=.TRUE.
2849	ELSE
2850	ll_next_to_land=.FALSE.
2851	ENDIF
2852
2853	! Reject observations
2854
2855	DO jobsp = kpstart(jobs), kpend(jobs)
2856
2857	! Flag if the observation falls outside the model spatial domain
2858	IF ( ( pobslam(jobs) < -180. ) &
2859	& .OR. ( pobslam(jobs) > 180. ) &
2860	& .OR. ( pobsphi(jobs) < -90. ) &
2861	& .OR. ( pobsphi(jobs) > 90. ) &
2862	& .OR. ( pobsdep(jobsp) < 0.0 ) &
2863	& .OR. ( pobsdep(jobsp) > gdepw_1d(kpk)) ) THEN
2864	kobsqc(jobsp) = kobsqc(jobsp) + 11
2865	kosdobs = kosdobs + 1
2866	CYCLE
2867	ENDIF
2868
2869	! To check if an observations falls within land there are two cases:
2870	! 1: z-coordibnates, where the check uses the mask
2871	! 2: terrain following (eg s-coordinates),
2872	! where we use the depth of the bottom cell to mask observations
2873
2874	IF( (.NOT. lk_vvl) .AND. ( ln_zps .OR. ln_zco ) ) THEN !(CASE 1)
2875
2876	! Flag if the observation falls with a model land cell
2877	IF ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) &
2878	& == 0.0_wp ) THEN
2879	kobsqc(jobsp) = kobsqc(jobsp) + 12
2880	klanobs = klanobs + 1
2881	CYCLE
2882	ENDIF
2883
2884	! Flag if the observation is close to land
2885	IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == &
2886	& 0.0_wp) THEN
2887	knlaobs = knlaobs + 1
2888	IF (ld_nea) THEN
2889	kobsqc(jobsp) = kobsqc(jobsp) + 14
2890	ENDIF
2891	ENDIF
2892
2893	ELSE ! Case 2
2894	! Flag if the observation is deeper than the bathymetry
2895	! Or if it is within the mask
2896	IF ( ANY( zgdepw(1:2,1:2,kpk,jobs) < pobsdep(jobsp) ) &
2897	& .OR. &
2898	& ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) &
2899	& == 0.0_wp) ) THEN
2900	kobsqc(jobsp) = kobsqc(jobsp) + 12
2901	klanobs = klanobs + 1
2902	CYCLE
2903	ENDIF
2904
2905	! Flag if the observation is close to land
2906	IF ( ll_next_to_land ) THEN
2907	knlaobs = knlaobs + 1
2908	IF (ld_nea) THEN
2909	kobsqc(jobsp) = kobsqc(jobsp) + 14
2910	ENDIF
2911	ENDIF
2912
2913	ENDIF
2914
2915	! For observations on the grid reject them if their are at
2916	! a masked point
2917
2918	IF (lgridobs) THEN
2919	IF (zgmsk(iig,ijg,kobsk(jobsp)-1,jobs) == 0.0_wp ) THEN
2920	kobsqc(jobsp) = kobsqc(jobsp) + 12
2921	klanobs = klanobs + 1
2922	CYCLE
2923	ENDIF
2924	ENDIF
2925
2926	! Set observation depth equal to that of the first model depth
2927	IF ( pobsdep(jobsp) <= pdep(1) ) THEN
2928	pobsdep(jobsp) = pdep(1)
2929	ENDIF
2930
2931	#if defined key_bdy
2932	! Flag if the observation falls close to the boundary rim
2933	IF (ld_bound_reject) THEN
2934	IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN
2935	kobsqc(jobsp) = kobsqc(jobsp) + 15
2936	kbdyobs = kbdyobs + 1
2937	CYCLE
2938	ENDIF
2939	! for observations on the grid...
2940	IF (lgridobs) THEN
2941	IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN
2942	kobsqc(jobsp) = kobsqc(jobsp) + 15
2943	kbdyobs = kbdyobs + 1
2944	CYCLE
2945	ENDIF
2946	ENDIF
2947	ENDIF
2948	#endif
2949
2950	END DO
2951	END DO
2952
2953	END SUBROUTINE obs_coo_spc_3d
2954
2955	SUBROUTINE obs_pro_rej( profdata )
2956	!!----------------------------------------------------------------------
2957	!! * ROUTINE obs_pro_rej *
2958	!!
2959	!! ** Purpose : Reject all data within a rejected profile
2960	!!
2961	!! ** Method :
2962	!!
2963	!! ** Action :
2964	!!
2965	!! References :
2966	!!
2967	!! History :
2968	!! ! 2007-10 (K. Mogensen) Original code
2969	!!----------------------------------------------------------------------
2970	!! * Modules used
2971	!! * Arguments
2972	TYPE(obs_prof), INTENT(INOUT) :: profdata ! Profile data
2973	!! * Local declarations
2974	INTEGER :: jprof
2975	INTEGER :: jvar
2976	INTEGER :: jobs
2977
2978	! Loop over profiles
2979
2980	DO jprof = 1, profdata%nprof
2981
2982	IF ( profdata%nqc(jprof) > 10 ) THEN
2983
2984	DO jvar = 1, profdata%nvar
2985
2986	DO jobs = profdata%npvsta(jprof,jvar), &
2987	& profdata%npvend(jprof,jvar)
2988
2989	profdata%var(jvar)%nvqc(jobs) = &
2990	& profdata%var(jvar)%nvqc(jobs) + 26
2991
2992	END DO
2993
2994	END DO
2995
2996	ENDIF
2997
2998	END DO
2999
3000	END SUBROUTINE obs_pro_rej
3001
3002	SUBROUTINE obs_uv_rej( profdata, knumu, knumv )
3003	!!----------------------------------------------------------------------
3004	!! * ROUTINE obs_uv_rej *
3005	!!
3006	!! ** Purpose : Reject u if v is rejected and vice versa
3007	!!
3008	!! ** Method :
3009	!!
3010	!! ** Action :
3011	!!
3012	!! References :
3013	!!
3014	!! History :
3015	!! ! 2009-2 (K. Mogensen) Original code
3016	!!----------------------------------------------------------------------
3017	!! * Modules used
3018	!! * Arguments
3019	TYPE(obs_prof), INTENT(INOUT) :: profdata ! Profile data
3020	INTEGER, INTENT(INOUT) :: knumu ! Number of u rejected
3021	INTEGER, INTENT(INOUT) :: knumv ! Number of v rejected
3022	!! * Local declarations
3023	INTEGER :: jprof
3024	INTEGER :: jvar
3025	INTEGER :: jobs
3026
3027	! Loop over profiles
3028
3029	DO jprof = 1, profdata%nprof
3030
3031	IF ( ( profdata%npvsta(jprof,1) /= profdata%npvsta(jprof,2) ) .OR. &
3032	& ( profdata%npvend(jprof,1) /= profdata%npvend(jprof,2) ) ) THEN
3033
3034	CALL ctl_stop('U,V profiles inconsistent in obs_uv_rej')
3035	RETURN
3036
3037	ENDIF
3038
3039	DO jobs = profdata%npvsta(jprof,1), profdata%npvend(jprof,1)
3040
3041	IF ( ( profdata%var(1)%nvqc(jobs) > 10 ) .AND. &
3042	& ( profdata%var(2)%nvqc(jobs) <= 10) ) THEN
3043	profdata%var(2)%nvqc(jobs) = profdata%var(2)%nvqc(jobs) + 42
3044	knumv = knumv + 1
3045	ENDIF
3046	IF ( ( profdata%var(2)%nvqc(jobs) > 10 ) .AND. &
3047	& ( profdata%var(1)%nvqc(jobs) <= 10) ) THEN
3048	profdata%var(1)%nvqc(jobs) = profdata%var(1)%nvqc(jobs) + 42
3049	knumu = knumu + 1
3050	ENDIF
3051
3052	END DO
3053
3054	END DO
3055
3056	END SUBROUTINE obs_uv_rej
3057
3058	END MODULE obs_prep

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

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

Context Navigation

source: branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 7713

Download in other formats: