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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

obs_prep.F90 in branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

Context Navigation

source: branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 6854

Last change on this file since 6854 was 6854, checked in by dford, 8 years ago
Initial implementation of observation operator for LogChl?.
File size: 90.1 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_scr : Basic screening of the observations
16	!! obs_coo_tim : Compute number of time steps to the observation time
17	!! obs_sor : Sort the observation arrays
18	!!---------------------------------------------------------------------
19	!! * Modules used
20	USE par_kind, ONLY : & ! Precision variables
21	& wp
22	USE in_out_manager ! I/O manager
23	USE obs_profiles_def ! Definitions for storage arrays for profiles
24	USE obs_surf_def ! Definitions for storage arrays for surface data
25	USE obs_mpp, ONLY : & ! MPP support routines for observation diagnostics
26	& obs_mpp_sum_integer, &
27	& obs_mpp_sum_integers
28	USE obs_inter_sup ! Interpolation support
29	USE obs_oper ! Observation operators
30	USE lib_mpp, ONLY : &
31	& ctl_warn, ctl_stop
32
33	IMPLICIT NONE
34
35	!! * Routine accessibility
36	PRIVATE
37
38	PUBLIC &
39	& obs_pre_pro, & ! First level check and screening of profiles
40	& obs_pre_sla, & ! First level check and screening of SLA data
41	& obs_pre_sst, & ! First level check and screening of SLA data
42	& obs_pre_seaice, & ! First level check and screening of sea ice data
43	& obs_pre_vel, & ! First level check and screening of velocity profiles
44	& obs_pre_logchl, & ! First level check and screening of logchl data
45	& calc_month_len ! Calculate the number of days in the months of a year
46
47	!!----------------------------------------------------------------------
48	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
49	!! $Id$
50	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
51	!!----------------------------------------------------------------------
52
53	!! * Substitutions
54	# include "domzgr_substitute.h90"
55
56	CONTAINS
57
58	SUBROUTINE obs_pre_pro( profdata, prodatqc, ld_t3d, ld_s3d, ld_nea, &
59	& kdailyavtypes )
60	!!----------------------------------------------------------------------
61	!! * ROUTINE obs_pre_pro *
62	!!
63	!! ** Purpose : First level check and screening of T and S profiles
64	!!
65	!! ** Method : First level check and screening of T and S profiles
66	!!
67	!! ** Action :
68	!!
69	!! References :
70	!!
71	!! History :
72	!! ! 2007-01 (K. Mogensen) Merge of obs_pre_t3d and obs_pre_s3d
73	!! ! 2007-03 (K. Mogensen) General handling of profiles
74	!! ! 2007-06 (K. Mogensen et al) Reject obs. near land.
75	!!----------------------------------------------------------------------
76	!! * Modules used
77	USE domstp ! Domain: set the time-step
78	USE par_oce ! Ocean parameters
79	USE dom_oce, ONLY : & ! Geographical information
80	& glamt, &
81	& gphit, &
82	& gdept_1d,&
83	& tmask
84	!! * Arguments
85	TYPE(obs_prof), INTENT(INOUT) :: profdata ! Full set of profile data
86	TYPE(obs_prof), INTENT(INOUT) :: prodatqc ! Subset of profile data not failing screening
87	LOGICAL, INTENT(IN) :: ld_t3d ! Switch for temperature
88	LOGICAL, INTENT(IN) :: ld_s3d ! Switch for salinity
89	LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land
90	INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
91	& kdailyavtypes! Types for daily averages
92	!! * Local declarations
93	INTEGER :: iyea0 ! Initial date
94	INTEGER :: imon0 ! - (year, month, day, hour, minute)
95	INTEGER :: iday0
96	INTEGER :: ihou0
97	INTEGER :: imin0
98	INTEGER :: icycle ! Current assimilation cycle
99	! Counters for observations that
100	INTEGER :: iotdobs ! - outside time domain
101	INTEGER :: iosdtobs ! - outside space domain (temperature)
102	INTEGER :: iosdsobs ! - outside space domain (salinity)
103	INTEGER :: ilantobs ! - within a model land cell (temperature)
104	INTEGER :: ilansobs ! - within a model land cell (salinity)
105	INTEGER :: inlatobs ! - close to land (temperature)
106	INTEGER :: inlasobs ! - close to land (salinity)
107	INTEGER :: igrdobs ! - fail the grid search
108	! Global counters for observations that
109	INTEGER :: iotdobsmpp ! - outside time domain
110	INTEGER :: iosdtobsmpp ! - outside space domain (temperature)
111	INTEGER :: iosdsobsmpp ! - outside space domain (salinity)
112	INTEGER :: ilantobsmpp ! - within a model land cell (temperature)
113	INTEGER :: ilansobsmpp ! - within a model land cell (salinity)
114	INTEGER :: inlatobsmpp ! - close to land (temperature)
115	INTEGER :: inlasobsmpp ! - close to land (salinity)
116	INTEGER :: igrdobsmpp ! - fail the grid search
117	TYPE(obs_prof_valid) :: llvalid ! Profile selection
118	TYPE(obs_prof_valid), DIMENSION(profdata%nvar) :: &
119	& llvvalid ! T,S selection
120	INTEGER :: jvar ! Variable loop variable
121	INTEGER :: jobs ! Obs. loop variable
122	INTEGER :: jstp ! Time loop variable
123	INTEGER :: inrc ! Time index variable
124
125	IF(lwp) WRITE(numout,*)'obs_pre_pro : Preparing the profile observations...'
126
127	! Initial date initialization (year, month, day, hour, minute)
128	iyea0 = ndate0 / 10000
129	imon0 = ( ndate0 - iyea0 * 10000 ) / 100
130	iday0 = ndate0 - iyea0 * 10000 - imon0 * 100
131	ihou0 = 0
132	imin0 = 0
133
134	icycle = no ! Assimilation cycle
135
136	! Diagnotics counters for various failures.
137
138	iotdobs = 0
139	igrdobs = 0
140	iosdtobs = 0
141	iosdsobs = 0
142	ilantobs = 0
143	ilansobs = 0
144	inlatobs = 0
145	inlasobs = 0
146
147	! -----------------------------------------------------------------------
148	! Find time coordinate for profiles
149	! -----------------------------------------------------------------------
150
151	IF ( PRESENT(kdailyavtypes) ) THEN
152	CALL obs_coo_tim_prof( icycle, &
153	& iyea0, imon0, iday0, ihou0, imin0, &
154	& profdata%nprof, profdata%nyea, profdata%nmon, &
155	& profdata%nday, profdata%nhou, profdata%nmin, &
156	& profdata%ntyp, profdata%nqc, profdata%mstp, &
157	& iotdobs, kdailyavtypes = kdailyavtypes )
158	ELSE
159	CALL obs_coo_tim_prof( icycle, &
160	& iyea0, imon0, iday0, ihou0, imin0, &
161	& profdata%nprof, profdata%nyea, profdata%nmon, &
162	& profdata%nday, profdata%nhou, profdata%nmin, &
163	& profdata%ntyp, profdata%nqc, profdata%mstp, &
164	& iotdobs )
165	ENDIF
166	CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
167
168	! -----------------------------------------------------------------------
169	! Check for profiles failing the grid search
170	! -----------------------------------------------------------------------
171
172	CALL obs_coo_grd( profdata%nprof, profdata%mi, profdata%mj, &
173	& profdata%nqc, igrdobs )
174
175	CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
176
177	! -----------------------------------------------------------------------
178	! Reject all observations for profiles with nqc > 10
179	! -----------------------------------------------------------------------
180
181	CALL obs_pro_rej( profdata )
182
183	! -----------------------------------------------------------------------
184	! Check for land points. This includes points below the model
185	! bathymetry so this is done for every point in the profile
186	! -----------------------------------------------------------------------
187
188	! Temperature
189
190	CALL obs_coo_spc_3d( profdata%nprof, profdata%nvprot(1), &
191	& profdata%npvsta(:,1), profdata%npvend(:,1), &
192	& jpi, jpj, &
193	& jpk, &
194	& profdata%mi, profdata%mj, &
195	& profdata%var(1)%mvk, &
196	& profdata%rlam, profdata%rphi, &
197	& profdata%var(1)%vdep, &
198	& glamt, gphit, &
199	& gdept_1d, tmask, &
200	& profdata%nqc, profdata%var(1)%nvqc, &
201	& iosdtobs, ilantobs, &
202	& inlatobs, ld_nea )
203
204	CALL obs_mpp_sum_integer( iosdtobs, iosdtobsmpp )
205	CALL obs_mpp_sum_integer( ilantobs, ilantobsmpp )
206	CALL obs_mpp_sum_integer( inlatobs, inlatobsmpp )
207
208	! Salinity
209
210	CALL obs_coo_spc_3d( profdata%nprof, profdata%nvprot(2), &
211	& profdata%npvsta(:,2), profdata%npvend(:,2), &
212	& jpi, jpj, &
213	& jpk, &
214	& profdata%mi, profdata%mj, &
215	& profdata%var(2)%mvk, &
216	& profdata%rlam, profdata%rphi, &
217	& profdata%var(2)%vdep, &
218	& glamt, gphit, &
219	& gdept_1d, tmask, &
220	& profdata%nqc, profdata%var(2)%nvqc, &
221	& iosdsobs, ilansobs, &
222	& inlasobs, ld_nea )
223
224	CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
225	CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
226	CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
227
228	! -----------------------------------------------------------------------
229	! Copy useful data from the profdata data structure to
230	! the prodatqc data structure
231	! -----------------------------------------------------------------------
232
233	! Allocate the selection arrays
234
235	ALLOCATE( llvalid%luse(profdata%nprof) )
236	DO jvar = 1,profdata%nvar
237	ALLOCATE( llvvalid(jvar)%luse(profdata%nvprot(jvar)) )
238	END DO
239
240	! We want all data which has qc flags <= 10
241
242	llvalid%luse(:) = ( profdata%nqc(:) <= 10 )
243	DO jvar = 1,profdata%nvar
244	llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= 10 )
245	END DO
246
247	! The actual copying
248
249	CALL obs_prof_compress( profdata, prodatqc, .TRUE., numout, &
250	& lvalid=llvalid, lvvalid=llvvalid )
251
252	! Dellocate the selection arrays
253	DEALLOCATE( llvalid%luse )
254	DO jvar = 1,profdata%nvar
255	DEALLOCATE( llvvalid(jvar)%luse )
256	END DO
257
258	! -----------------------------------------------------------------------
259	! Print information about what observations are left after qc
260	! -----------------------------------------------------------------------
261
262	! Update the total observation counter array
263
264	IF(lwp) THEN
265	WRITE(numout,*)
266	WRITE(numout,*) 'obs_pre_pro :'
267	WRITE(numout,*) '~~~~~~~~~~~'
268	WRITE(numout,*)
269	WRITE(numout,*) ' Profiles outside time domain = ', &
270	& iotdobsmpp
271	WRITE(numout,*) ' Remaining profiles that failed grid search = ', &
272	& igrdobsmpp
273	WRITE(numout,*) ' Remaining T data outside space domain = ', &
274	& iosdtobsmpp
275	WRITE(numout,*) ' Remaining T data at land points = ', &
276	& ilantobsmpp
277	IF (ld_nea) THEN
278	WRITE(numout,*) ' Remaining T data near land points (removed) = ',&
279	& inlatobsmpp
280	ELSE
281	WRITE(numout,*) ' Remaining T data near land points (kept) = ',&
282	& inlatobsmpp
283	ENDIF
284	WRITE(numout,*) ' T data accepted = ', &
285	& prodatqc%nvprotmpp(1)
286	WRITE(numout,*) ' Remaining S data outside space domain = ', &
287	& iosdsobsmpp
288	WRITE(numout,*) ' Remaining S data at land points = ', &
289	& ilansobsmpp
290	IF (ld_nea) THEN
291	WRITE(numout,*) ' Remaining S data near land points (removed) = ',&
292	& inlasobsmpp
293	ELSE
294	WRITE(numout,*) ' Remaining S data near land points (kept) = ',&
295	& inlasobsmpp
296	ENDIF
297	WRITE(numout,*) ' S data accepted = ', &
298	& prodatqc%nvprotmpp(2)
299
300	WRITE(numout,*)
301	WRITE(numout,*) ' Number of observations per time step :'
302	WRITE(numout,*)
303	WRITE(numout,997)
304	WRITE(numout,998)
305	ENDIF
306
307	DO jobs = 1, prodatqc%nprof
308	inrc = prodatqc%mstp(jobs) + 2 - nit000
309	prodatqc%npstp(inrc) = prodatqc%npstp(inrc) + 1
310	DO jvar = 1, prodatqc%nvar
311	IF ( prodatqc%npvend(jobs,jvar) > 0 ) THEN
312	prodatqc%nvstp(inrc,jvar) = prodatqc%nvstp(inrc,jvar) + &
313	& ( prodatqc%npvend(jobs,jvar) - &
314	& prodatqc%npvsta(jobs,jvar) + 1 )
315	ENDIF
316	END DO
317	END DO
318
319
320	CALL obs_mpp_sum_integers( prodatqc%npstp, prodatqc%npstpmpp, &
321	& nitend - nit000 + 2 )
322	DO jvar = 1, prodatqc%nvar
323	CALL obs_mpp_sum_integers( prodatqc%nvstp(:,jvar), &
324	& prodatqc%nvstpmpp(:,jvar), &
325	& nitend - nit000 + 2 )
326	END DO
327
328	IF ( lwp ) THEN
329	DO jstp = nit000 - 1, nitend
330	inrc = jstp - nit000 + 2
331	WRITE(numout,999) jstp, prodatqc%npstpmpp(inrc), &
332	& prodatqc%nvstpmpp(inrc,1), &
333	& prodatqc%nvstpmpp(inrc,2)
334	END DO
335	ENDIF
336
337	997 FORMAT(10X,'Time step',5X,'Profiles',5X,'Temperature',5X,'Salinity')
338	998 FORMAT(10X,'---------',5X,'--------',5X,'-----------',5X,'--------')
339	999 FORMAT(10X,I9,5X,I8,5X,I11,5X,I8)
340
341	END SUBROUTINE obs_pre_pro
342
343	SUBROUTINE obs_pre_sla( sladata, sladatqc, ld_sla, ld_nea )
344	!!----------------------------------------------------------------------
345	!! * ROUTINE obs_pre_sla *
346	!!
347	!! ** Purpose : First level check and screening of SLA observations
348	!!
349	!! ** Method : First level check and screening of SLA observations
350	!!
351	!! ** Action :
352	!!
353	!! References :
354	!!
355	!! History :
356	!! ! 2007-03 (A. Weaver, K. Mogensen) Original
357	!! ! 2007-06 (K. Mogensen et al) Reject obs. near land.
358	!!----------------------------------------------------------------------
359	!! * Modules used
360	USE domstp ! Domain: set the time-step
361	USE par_oce ! Ocean parameters
362	USE dom_oce, ONLY : & ! Geographical information
363	& glamt, &
364	& gphit, &
365	& tmask
366	!! * Arguments
367	TYPE(obs_surf), INTENT(INOUT) :: sladata ! Full set of SLA data
368	TYPE(obs_surf), INTENT(INOUT) :: sladatqc ! Subset of SLA data not failing screening
369	LOGICAL, INTENT(IN) :: ld_sla ! Switch for SLA data
370	LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land
371	!! * Local declarations
372	INTEGER :: iyea0 ! Initial date
373	INTEGER :: imon0 ! - (year, month, day, hour, minute)
374	INTEGER :: iday0
375	INTEGER :: ihou0
376	INTEGER :: imin0
377	INTEGER :: icycle ! Current assimilation cycle
378	! Counters for observations that
379	INTEGER :: iotdobs ! - outside time domain
380	INTEGER :: iosdsobs ! - outside space domain
381	INTEGER :: ilansobs ! - within a model land cell
382	INTEGER :: inlasobs ! - close to land
383	INTEGER :: igrdobs ! - fail the grid search
384	! Global counters for observations that
385	INTEGER :: iotdobsmpp ! - outside time domain
386	INTEGER :: iosdsobsmpp ! - outside space domain
387	INTEGER :: ilansobsmpp ! - within a model land cell
388	INTEGER :: inlasobsmpp ! - close to land
389	INTEGER :: igrdobsmpp ! - fail the grid search
390	LOGICAL, DIMENSION(:), ALLOCATABLE :: &
391	& llvalid ! SLA data selection
392	INTEGER :: jobs ! Obs. loop variable
393	INTEGER :: jstp ! Time loop variable
394	INTEGER :: inrc ! Time index variable
395
396	IF(lwp) WRITE(numout,*)'obs_pre_sla : Preparing the SLA observations...'
397
398	! Initial date initialization (year, month, day, hour, minute)
399	iyea0 = ndate0 / 10000
400	imon0 = ( ndate0 - iyea0 * 10000 ) / 100
401	iday0 = ndate0 - iyea0 * 10000 - imon0 * 100
402	ihou0 = 0
403	imin0 = 0
404
405	icycle = no ! Assimilation cycle
406
407	! Diagnotics counters for various failures.
408
409	iotdobs = 0
410	igrdobs = 0
411	iosdsobs = 0
412	ilansobs = 0
413	inlasobs = 0
414
415	! -----------------------------------------------------------------------
416	! Find time coordinate for SLA data
417	! -----------------------------------------------------------------------
418
419	CALL obs_coo_tim( icycle, &
420	& iyea0, imon0, iday0, ihou0, imin0, &
421	& sladata%nsurf, sladata%nyea, sladata%nmon, &
422	& sladata%nday, sladata%nhou, sladata%nmin, &
423	& sladata%nqc, sladata%mstp, iotdobs )
424
425	CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
426
427	! -----------------------------------------------------------------------
428	! Check for SLA data failing the grid search
429	! -----------------------------------------------------------------------
430
431	CALL obs_coo_grd( sladata%nsurf, sladata%mi, sladata%mj, &
432	& sladata%nqc, igrdobs )
433
434	CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
435
436	! -----------------------------------------------------------------------
437	! Check for land points.
438	! -----------------------------------------------------------------------
439
440	CALL obs_coo_spc_2d( sladata%nsurf, &
441	& jpi, jpj, &
442	& sladata%mi, sladata%mj, &
443	& sladata%rlam, sladata%rphi, &
444	& glamt, gphit, &
445	& tmask(:,:,1), sladata%nqc, &
446	& iosdsobs, ilansobs, &
447	& inlasobs, ld_nea )
448
449	CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
450	CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
451	CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
452
453	! -----------------------------------------------------------------------
454	! Copy useful data from the sladata data structure to
455	! the sladatqc data structure
456	! -----------------------------------------------------------------------
457
458	! Allocate the selection arrays
459
460	ALLOCATE( llvalid(sladata%nsurf) )
461
462	! We want all data which has qc flags <= 10
463
464	llvalid(:) = ( sladata%nqc(:) <= 10 )
465
466	! The actual copying
467
468	CALL obs_surf_compress( sladata, sladatqc, .TRUE., numout, &
469	& lvalid=llvalid )
470
471	! Dellocate the selection arrays
472	DEALLOCATE( llvalid )
473
474	! -----------------------------------------------------------------------
475	! Print information about what observations are left after qc
476	! -----------------------------------------------------------------------
477
478	! Update the total observation counter array
479
480	IF(lwp) THEN
481	WRITE(numout,*)
482	WRITE(numout,*) 'obs_pre_sla :'
483	WRITE(numout,*) '~~~~~~~~~~~'
484	WRITE(numout,*)
485	WRITE(numout,*) ' SLA data outside time domain = ', &
486	& iotdobsmpp
487	WRITE(numout,*) ' Remaining SLA data that failed grid search = ', &
488	& igrdobsmpp
489	WRITE(numout,*) ' Remaining SLA data outside space domain = ', &
490	& iosdsobsmpp
491	WRITE(numout,*) ' Remaining SLA data at land points = ', &
492	& ilansobsmpp
493	IF (ld_nea) THEN
494	WRITE(numout,*) ' Remaining SLA data near land points (removed) = ', &
495	& inlasobsmpp
496	ELSE
497	WRITE(numout,*) ' Remaining SLA data near land points (kept) = ', &
498	& inlasobsmpp
499	ENDIF
500	WRITE(numout,*) ' SLA data accepted = ', &
501	& sladatqc%nsurfmpp
502
503	WRITE(numout,*)
504	WRITE(numout,*) ' Number of observations per time step :'
505	WRITE(numout,*)
506	WRITE(numout,1997)
507	WRITE(numout,1998)
508	ENDIF
509
510	DO jobs = 1, sladatqc%nsurf
511	inrc = sladatqc%mstp(jobs) + 2 - nit000
512	sladatqc%nsstp(inrc) = sladatqc%nsstp(inrc) + 1
513	END DO
514
515	CALL obs_mpp_sum_integers( sladatqc%nsstp, sladatqc%nsstpmpp, &
516	& nitend - nit000 + 2 )
517
518	IF ( lwp ) THEN
519	DO jstp = nit000 - 1, nitend
520	inrc = jstp - nit000 + 2
521	WRITE(numout,1999) jstp, sladatqc%nsstpmpp(inrc)
522	END DO
523	ENDIF
524
525	1997 FORMAT(10X,'Time step',5X,'Sea level anomaly')
526	1998 FORMAT(10X,'---------',5X,'-----------------')
527	1999 FORMAT(10X,I9,5X,I17)
528
529	END SUBROUTINE obs_pre_sla
530
531	SUBROUTINE obs_pre_sst( sstdata, sstdatqc, ld_sst, ld_nea )
532	!!----------------------------------------------------------------------
533	!! * ROUTINE obs_pre_sst *
534	!!
535	!! ** Purpose : First level check and screening of SST observations
536	!!
537	!! ** Method : First level check and screening of SST observations
538	!!
539	!! ** Action :
540	!!
541	!! References :
542	!!
543	!! History :
544	!! ! 2007-03 (S. Ricci) SST data preparation
545	!!----------------------------------------------------------------------
546	!! * Modules used
547	USE domstp ! Domain: set the time-step
548	USE par_oce ! Ocean parameters
549	USE dom_oce, ONLY : & ! Geographical information
550	& glamt, &
551	& gphit, &
552	& tmask
553	!! * Arguments
554	TYPE(obs_surf), INTENT(INOUT) :: sstdata ! Full set of SST data
555	TYPE(obs_surf), INTENT(INOUT) :: sstdatqc ! Subset of SST data not failing screening
556	LOGICAL :: ld_sst ! Switch for SST data
557	LOGICAL :: ld_nea ! Switch for rejecting observation near land
558	!! * Local declarations
559	INTEGER :: iyea0 ! Initial date
560	INTEGER :: imon0 ! - (year, month, day, hour, minute)
561	INTEGER :: iday0
562	INTEGER :: ihou0
563	INTEGER :: imin0
564	INTEGER :: icycle ! Current assimilation cycle
565	! Counters for observations that
566	INTEGER :: iotdobs ! - outside time domain
567	INTEGER :: iosdsobs ! - outside space domain
568	INTEGER :: ilansobs ! - within a model land cell
569	INTEGER :: inlasobs ! - close to land
570	INTEGER :: igrdobs ! - fail the grid search
571	! Global counters for observations that
572	INTEGER :: iotdobsmpp ! - outside time domain
573	INTEGER :: iosdsobsmpp ! - outside space domain
574	INTEGER :: ilansobsmpp ! - within a model land cell
575	INTEGER :: inlasobsmpp ! - close to land
576	INTEGER :: igrdobsmpp ! - fail the grid search
577	LOGICAL, DIMENSION(:), ALLOCATABLE :: &
578	& llvalid ! SST data selection
579	INTEGER :: jobs ! Obs. loop variable
580	INTEGER :: jstp ! Time loop variable
581	INTEGER :: inrc ! Time index variable
582
583	IF(lwp) WRITE(numout,*)'obs_pre_sst : Preparing the SST observations...'
584
585	! Initial date initialization (year, month, day, hour, minute)
586	iyea0 = ndate0 / 10000
587	imon0 = ( ndate0 - iyea0 * 10000 ) / 100
588	iday0 = ndate0 - iyea0 * 10000 - imon0 * 100
589	ihou0 = 0
590	imin0 = 0
591
592	icycle = no ! Assimilation cycle
593
594	! Diagnotics counters for various failures.
595
596	iotdobs = 0
597	igrdobs = 0
598	iosdsobs = 0
599	ilansobs = 0
600	inlasobs = 0
601
602	! -----------------------------------------------------------------------
603	! Find time coordinate for SST data
604	! -----------------------------------------------------------------------
605
606	CALL obs_coo_tim( icycle, &
607	& iyea0, imon0, iday0, ihou0, imin0, &
608	& sstdata%nsurf, sstdata%nyea, sstdata%nmon, &
609	& sstdata%nday, sstdata%nhou, sstdata%nmin, &
610	& sstdata%nqc, sstdata%mstp, iotdobs )
611	CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
612	! -----------------------------------------------------------------------
613	! Check for SST data failing the grid search
614	! -----------------------------------------------------------------------
615
616	CALL obs_coo_grd( sstdata%nsurf, sstdata%mi, sstdata%mj, &
617	& sstdata%nqc, igrdobs )
618	CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
619
620	! -----------------------------------------------------------------------
621	! Check for land points.
622	! -----------------------------------------------------------------------
623
624	CALL obs_coo_spc_2d( sstdata%nsurf, &
625	& jpi, jpj, &
626	& sstdata%mi, sstdata%mj, &
627	& sstdata%rlam, sstdata%rphi, &
628	& glamt, gphit, &
629	& tmask(:,:,1), sstdata%nqc, &
630	& iosdsobs, ilansobs, &
631	& inlasobs, ld_nea )
632
633	CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
634	CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
635	CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
636
637	! -----------------------------------------------------------------------
638	! Copy useful data from the sstdata data structure to
639	! the sstdatqc data structure
640	! -----------------------------------------------------------------------
641
642	! Allocate the selection arrays
643
644	ALLOCATE( llvalid(sstdata%nsurf) )
645
646	! We want all data which has qc flags <= 0
647
648	llvalid(:) = ( sstdata%nqc(:) <= 10 )
649
650	! The actual copying
651
652	CALL obs_surf_compress( sstdata, sstdatqc, .TRUE., numout, &
653	& lvalid=llvalid )
654
655	! Dellocate the selection arrays
656	DEALLOCATE( llvalid )
657
658	! -----------------------------------------------------------------------
659	! Print information about what observations are left after qc
660	! -----------------------------------------------------------------------
661
662	! Update the total observation counter array
663
664	IF(lwp) THEN
665	WRITE(numout,*)
666	WRITE(numout,*) 'obs_pre_sst :'
667	WRITE(numout,*) '~~~~~~~~~~~'
668	WRITE(numout,*)
669	WRITE(numout,*) ' SST data outside time domain = ', &
670	& iotdobsmpp
671	WRITE(numout,*) ' Remaining SST data that failed grid search = ', &
672	& igrdobsmpp
673	WRITE(numout,*) ' Remaining SST data outside space domain = ', &
674	& iosdsobsmpp
675	WRITE(numout,*) ' Remaining SST data at land points = ', &
676	& ilansobsmpp
677	IF (ld_nea) THEN
678	WRITE(numout,*) ' Remaining SST data near land points (removed) = ', &
679	& inlasobsmpp
680	ELSE
681	WRITE(numout,*) ' Remaining SST data near land points (kept) = ', &
682	& inlasobsmpp
683	ENDIF
684	WRITE(numout,*) ' SST data accepted = ', &
685	& sstdatqc%nsurfmpp
686
687	WRITE(numout,*)
688	WRITE(numout,*) ' Number of observations per time step :'
689	WRITE(numout,*)
690	WRITE(numout,1997)
691	WRITE(numout,1998)
692	ENDIF
693
694	DO jobs = 1, sstdatqc%nsurf
695	inrc = sstdatqc%mstp(jobs) + 2 - nit000
696	sstdatqc%nsstp(inrc) = sstdatqc%nsstp(inrc) + 1
697	END DO
698
699	CALL obs_mpp_sum_integers( sstdatqc%nsstp, sstdatqc%nsstpmpp, &
700	& nitend - nit000 + 2 )
701
702	IF ( lwp ) THEN
703	DO jstp = nit000 - 1, nitend
704	inrc = jstp - nit000 + 2
705	WRITE(numout,1999) jstp, sstdatqc%nsstpmpp(inrc)
706	END DO
707	ENDIF
708
709	1997 FORMAT(10X,'Time step',5X,'Sea surface temperature')
710	1998 FORMAT(10X,'---------',5X,'-----------------')
711	1999 FORMAT(10X,I9,5X,I17)
712
713	END SUBROUTINE obs_pre_sst
714
715	SUBROUTINE obs_pre_seaice( seaicedata, seaicedatqc, ld_seaice, ld_nea )
716	!!----------------------------------------------------------------------
717	!! * ROUTINE obs_pre_seaice *
718	!!
719	!! ** Purpose : First level check and screening of Sea Ice observations
720	!!
721	!! ** Method : First level check and screening of Sea Ice observations
722	!!
723	!! ** Action :
724	!!
725	!! References :
726	!!
727	!! History :
728	!! ! 2007-11 (D. Lea) based on obs_pre_sst
729	!!----------------------------------------------------------------------
730	!! * Modules used
731	USE domstp ! Domain: set the time-step
732	USE par_oce ! Ocean parameters
733	USE dom_oce, ONLY : & ! Geographical information
734	& glamt, &
735	& gphit, &
736	& tmask
737	!! * Arguments
738	TYPE(obs_surf), INTENT(INOUT) :: seaicedata ! Full set of Sea Ice data
739	TYPE(obs_surf), INTENT(INOUT) :: seaicedatqc ! Subset of sea ice data not failing screening
740	LOGICAL :: ld_seaice ! Switch for sea ice data
741	LOGICAL :: ld_nea ! Switch for rejecting observation near land
742	!! * Local declarations
743	INTEGER :: iyea0 ! Initial date
744	INTEGER :: imon0 ! - (year, month, day, hour, minute)
745	INTEGER :: iday0
746	INTEGER :: ihou0
747	INTEGER :: imin0
748	INTEGER :: icycle ! Current assimilation cycle
749	! Counters for observations that
750	INTEGER :: iotdobs ! - outside time domain
751	INTEGER :: iosdsobs ! - outside space domain
752	INTEGER :: ilansobs ! - within a model land cell
753	INTEGER :: inlasobs ! - close to land
754	INTEGER :: igrdobs ! - fail the grid search
755	! Global counters for observations that
756	INTEGER :: iotdobsmpp ! - outside time domain
757	INTEGER :: iosdsobsmpp ! - outside space domain
758	INTEGER :: ilansobsmpp ! - within a model land cell
759	INTEGER :: inlasobsmpp ! - close to land
760	INTEGER :: igrdobsmpp ! - fail the grid search
761	LOGICAL, DIMENSION(:), ALLOCATABLE :: &
762	& llvalid ! data selection
763	INTEGER :: jobs ! Obs. loop variable
764	INTEGER :: jstp ! Time loop variable
765	INTEGER :: inrc ! Time index variable
766
767	IF (lwp) WRITE(numout,*)'obs_pre_seaice : Preparing the sea ice observations...'
768
769	! Initial date initialization (year, month, day, hour, minute)
770	iyea0 = ndate0 / 10000
771	imon0 = ( ndate0 - iyea0 * 10000 ) / 100
772	iday0 = ndate0 - iyea0 * 10000 - imon0 * 100
773	ihou0 = 0
774	imin0 = 0
775
776	icycle = no ! Assimilation cycle
777
778	! Diagnotics counters for various failures.
779
780	iotdobs = 0
781	igrdobs = 0
782	iosdsobs = 0
783	ilansobs = 0
784	inlasobs = 0
785
786	! -----------------------------------------------------------------------
787	! Find time coordinate for sea ice data
788	! -----------------------------------------------------------------------
789
790	CALL obs_coo_tim( icycle, &
791	& iyea0, imon0, iday0, ihou0, imin0, &
792	& seaicedata%nsurf, seaicedata%nyea, seaicedata%nmon, &
793	& seaicedata%nday, seaicedata%nhou, seaicedata%nmin, &
794	& seaicedata%nqc, seaicedata%mstp, iotdobs )
795	CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
796	! -----------------------------------------------------------------------
797	! Check for sea ice data failing the grid search
798	! -----------------------------------------------------------------------
799
800	CALL obs_coo_grd( seaicedata%nsurf, seaicedata%mi, seaicedata%mj, &
801	& seaicedata%nqc, igrdobs )
802	CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
803
804	! -----------------------------------------------------------------------
805	! Check for land points.
806	! -----------------------------------------------------------------------
807
808	CALL obs_coo_spc_2d( seaicedata%nsurf, &
809	& jpi, jpj, &
810	& seaicedata%mi, seaicedata%mj, &
811	& seaicedata%rlam, seaicedata%rphi, &
812	& glamt, gphit, &
813	& tmask(:,:,1), seaicedata%nqc, &
814	& iosdsobs, ilansobs, &
815	& inlasobs, ld_nea )
816
817	CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
818	CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
819	CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
820
821	! -----------------------------------------------------------------------
822	! Copy useful data from the seaicedata data structure to
823	! the seaicedatqc data structure
824	! -----------------------------------------------------------------------
825
826	! Allocate the selection arrays
827
828	ALLOCATE( llvalid(seaicedata%nsurf) )
829
830	! We want all data which has qc flags <= 0
831
832	llvalid(:) = ( seaicedata%nqc(:) <= 10 )
833
834	! The actual copying
835
836	CALL obs_surf_compress( seaicedata, seaicedatqc, .TRUE., numout, &
837	& lvalid=llvalid )
838
839	! Dellocate the selection arrays
840	DEALLOCATE( llvalid )
841
842	! -----------------------------------------------------------------------
843	! Print information about what observations are left after qc
844	! -----------------------------------------------------------------------
845
846	! Update the total observation counter array
847
848	IF(lwp) THEN
849	WRITE(numout,*)
850	WRITE(numout,*) 'obs_pre_seaice :'
851	WRITE(numout,*) '~~~~~~~~~~~'
852	WRITE(numout,*)
853	WRITE(numout,*) ' Sea ice data outside time domain = ', &
854	& iotdobsmpp
855	WRITE(numout,*) ' Remaining sea ice data that failed grid search = ', &
856	& igrdobsmpp
857	WRITE(numout,*) ' Remaining sea ice data outside space domain = ', &
858	& iosdsobsmpp
859	WRITE(numout,*) ' Remaining sea ice data at land points = ', &
860	& ilansobsmpp
861	IF (ld_nea) THEN
862	WRITE(numout,*) ' Remaining sea ice data near land points (removed) = ', &
863	& inlasobsmpp
864	ELSE
865	WRITE(numout,*) ' Remaining sea ice data near land points (kept) = ', &
866	& inlasobsmpp
867	ENDIF
868	WRITE(numout,*) ' Sea ice data accepted = ', &
869	& seaicedatqc%nsurfmpp
870
871	WRITE(numout,*)
872	WRITE(numout,*) ' Number of observations per time step :'
873	WRITE(numout,*)
874	WRITE(numout,1997)
875	WRITE(numout,1998)
876	ENDIF
877
878	DO jobs = 1, seaicedatqc%nsurf
879	inrc = seaicedatqc%mstp(jobs) + 2 - nit000
880	seaicedatqc%nsstp(inrc) = seaicedatqc%nsstp(inrc) + 1
881	END DO
882
883	CALL obs_mpp_sum_integers( seaicedatqc%nsstp, seaicedatqc%nsstpmpp, &
884	& nitend - nit000 + 2 )
885
886	IF ( lwp ) THEN
887	DO jstp = nit000 - 1, nitend
888	inrc = jstp - nit000 + 2
889	WRITE(numout,1999) jstp, seaicedatqc%nsstpmpp(inrc)
890	END DO
891	ENDIF
892
893	1997 FORMAT(10X,'Time step',5X,'Sea ice data ')
894	1998 FORMAT(10X,'---------',5X,'-----------------')
895	1999 FORMAT(10X,I9,5X,I17)
896
897	END SUBROUTINE obs_pre_seaice
898
899	SUBROUTINE obs_pre_vel( profdata, prodatqc, ld_vel3d, ld_nea, ld_dailyav )
900	!!----------------------------------------------------------------------
901	!! * ROUTINE obs_pre_taovel *
902	!!
903	!! ** Purpose : First level check and screening of U and V profiles
904	!!
905	!! ** Method : First level check and screening of U and V profiles
906	!!
907	!! History :
908	!! ! 2007-06 (K. Mogensen) original : T and S profile data
909	!! ! 2008-09 (M. Valdivieso) : TAO velocity data
910	!! ! 2009-01 (K. Mogensen) : New feedback strictuer
911	!!
912	!!----------------------------------------------------------------------
913	!! * Modules used
914	USE domstp ! Domain: set the time-step
915	USE par_oce ! Ocean parameters
916	USE dom_oce, ONLY : & ! Geographical information
917	& glamt, glamu, glamv, &
918	& gphit, gphiu, gphiv, &
919	& gdept_1d, &
920	& tmask, umask, vmask
921	!! * Arguments
922	TYPE(obs_prof), INTENT(INOUT) :: profdata ! Full set of profile data
923	TYPE(obs_prof), INTENT(INOUT) :: prodatqc ! Subset of profile data not failing screening
924	LOGICAL, INTENT(IN) :: ld_vel3d ! Switch for zonal and meridional velocity components
925	LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land
926	LOGICAL, INTENT(IN) :: ld_dailyav ! Switch for daily average data
927	!! * Local declarations
928	INTEGER :: iyea0 ! Initial date
929	INTEGER :: imon0 ! - (year, month, day, hour, minute)
930	INTEGER :: iday0
931	INTEGER :: ihou0
932	INTEGER :: imin0
933	INTEGER :: icycle ! Current assimilation cycle
934	! Counters for observations that
935	INTEGER :: iotdobs ! - outside time domain
936	INTEGER :: iosduobs ! - outside space domain (zonal velocity component)
937	INTEGER :: iosdvobs ! - outside space domain (meridional velocity component)
938	INTEGER :: ilanuobs ! - within a model land cell (zonal velocity component)
939	INTEGER :: ilanvobs ! - within a model land cell (meridional velocity component)
940	INTEGER :: inlauobs ! - close to land (zonal velocity component)
941	INTEGER :: inlavobs ! - close to land (meridional velocity component)
942	INTEGER :: igrdobs ! - fail the grid search
943	INTEGER :: iuvchku ! - reject u if v rejected and vice versa
944	INTEGER :: iuvchkv !
945	! Global counters for observations that
946	INTEGER :: iotdobsmpp ! - outside time domain
947	INTEGER :: iosduobsmpp ! - outside space domain (zonal velocity component)
948	INTEGER :: iosdvobsmpp ! - outside space domain (meridional velocity component)
949	INTEGER :: ilanuobsmpp ! - within a model land cell (zonal velocity component)
950	INTEGER :: ilanvobsmpp ! - within a model land cell (meridional velocity component)
951	INTEGER :: inlauobsmpp ! - close to land (zonal velocity component)
952	INTEGER :: inlavobsmpp ! - close to land (meridional velocity component)
953	INTEGER :: igrdobsmpp ! - fail the grid search
954	INTEGER :: iuvchkumpp ! - reject u if v rejected and vice versa
955	INTEGER :: iuvchkvmpp !
956	TYPE(obs_prof_valid) :: llvalid ! Profile selection
957	TYPE(obs_prof_valid), DIMENSION(profdata%nvar) :: &
958	& llvvalid ! U,V selection
959	INTEGER :: jvar ! Variable loop variable
960	INTEGER :: jobs ! Obs. loop variable
961	INTEGER :: jstp ! Time loop variable
962	INTEGER :: inrc ! Time index variable
963
964	IF(lwp) WRITE(numout,*)'obs_pre_vel: Preparing the velocity profile data'
965
966	! Initial date initialization (year, month, day, hour, minute)
967	iyea0 = ndate0 / 10000
968	imon0 = ( ndate0 - iyea0 * 10000 ) / 100
969	iday0 = ndate0 - iyea0 * 10000 - imon0 * 100
970	ihou0 = 0
971	imin0 = 0
972
973	icycle = no ! Assimilation cycle
974
975	! Diagnotics counters for various failures.
976
977	iotdobs = 0
978	igrdobs = 0
979	iosduobs = 0
980	iosdvobs = 0
981	ilanuobs = 0
982	ilanvobs = 0
983	inlauobs = 0
984	inlavobs = 0
985	iuvchku = 0
986	iuvchkv = 0
987
988	! -----------------------------------------------------------------------
989	! Find time coordinate for profiles
990	! -----------------------------------------------------------------------
991
992	CALL obs_coo_tim_prof( icycle, &
993	& iyea0, imon0, iday0, ihou0, imin0, &
994	& profdata%nprof, profdata%nyea, profdata%nmon, &
995	& profdata%nday, profdata%nhou, profdata%nmin, &
996	& profdata%ntyp, profdata%nqc, profdata%mstp, &
997	& iotdobs, ld_dailyav = ld_dailyav )
998
999	CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
1000
1001	! -----------------------------------------------------------------------
1002	! Check for profiles failing the grid search
1003	! -----------------------------------------------------------------------
1004
1005	CALL obs_coo_grd( profdata%nprof, profdata%mi(:,1), profdata%mj(:,1), &
1006	& profdata%nqc, igrdobs )
1007	CALL obs_coo_grd( profdata%nprof, profdata%mi(:,2), profdata%mj(:,2), &
1008	& profdata%nqc, igrdobs )
1009
1010	CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
1011
1012	! -----------------------------------------------------------------------
1013	! Reject all observations for profiles with nqc > 10
1014	! -----------------------------------------------------------------------
1015
1016	CALL obs_pro_rej( profdata )
1017
1018	! -----------------------------------------------------------------------
1019	! Check for land points. This includes points below the model
1020	! bathymetry so this is done for every point in the profile
1021	! -----------------------------------------------------------------------
1022
1023	! Zonal Velocity Component
1024
1025	CALL obs_coo_spc_3d( profdata%nprof, profdata%nvprot(1), &
1026	& profdata%npvsta(:,1), profdata%npvend(:,1), &
1027	& jpi, jpj, &
1028	& jpk, &
1029	& profdata%mi, profdata%mj, &
1030	& profdata%var(1)%mvk, &
1031	& profdata%rlam, profdata%rphi, &
1032	& profdata%var(1)%vdep, &
1033	& glamu, gphiu, &
1034	& gdept_1d, umask, &
1035	& profdata%nqc, profdata%var(1)%nvqc, &
1036	& iosduobs, ilanuobs, &
1037	& inlauobs, ld_nea )
1038
1039	CALL obs_mpp_sum_integer( iosduobs, iosduobsmpp )
1040	CALL obs_mpp_sum_integer( ilanuobs, ilanuobsmpp )
1041	CALL obs_mpp_sum_integer( inlauobs, inlauobsmpp )
1042
1043	! Meridional Velocity Component
1044
1045	CALL obs_coo_spc_3d( profdata%nprof, profdata%nvprot(2), &
1046	& profdata%npvsta(:,2), profdata%npvend(:,2), &
1047	& jpi, jpj, &
1048	& jpk, &
1049	& profdata%mi, profdata%mj, &
1050	& profdata%var(2)%mvk, &
1051	& profdata%rlam, profdata%rphi, &
1052	& profdata%var(2)%vdep, &
1053	& glamv, gphiv, &
1054	& gdept_1d, vmask, &
1055	& profdata%nqc, profdata%var(2)%nvqc, &
1056	& iosdvobs, ilanvobs, &
1057	& inlavobs, ld_nea )
1058
1059	CALL obs_mpp_sum_integer( iosdvobs, iosdvobsmpp )
1060	CALL obs_mpp_sum_integer( ilanvobs, ilanvobsmpp )
1061	CALL obs_mpp_sum_integer( inlavobs, inlavobsmpp )
1062
1063	! -----------------------------------------------------------------------
1064	! Reject u if v is rejected and vice versa
1065	! -----------------------------------------------------------------------
1066
1067	CALL obs_uv_rej( profdata, iuvchku, iuvchkv )
1068	CALL obs_mpp_sum_integer( iuvchku, iuvchkumpp )
1069	CALL obs_mpp_sum_integer( iuvchkv, iuvchkvmpp )
1070
1071	! -----------------------------------------------------------------------
1072	! Copy useful data from the profdata data structure to
1073	! the prodatqc data structure
1074	! -----------------------------------------------------------------------
1075
1076	! Allocate the selection arrays
1077
1078	ALLOCATE( llvalid%luse(profdata%nprof) )
1079	DO jvar = 1,profdata%nvar
1080	ALLOCATE( llvvalid(jvar)%luse(profdata%nvprot(jvar)) )
1081	END DO
1082
1083	! We want all data which has qc flags = 0
1084
1085	llvalid%luse(:) = ( profdata%nqc(:) <= 10 )
1086	DO jvar = 1,profdata%nvar
1087	llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= 10 )
1088	END DO
1089
1090	! The actual copying
1091
1092	CALL obs_prof_compress( profdata, prodatqc, .TRUE., numout, &
1093	& lvalid=llvalid, lvvalid=llvvalid )
1094
1095	! Dellocate the selection arrays
1096	DEALLOCATE( llvalid%luse )
1097	DO jvar = 1,profdata%nvar
1098	DEALLOCATE( llvvalid(jvar)%luse )
1099	END DO
1100
1101	! -----------------------------------------------------------------------
1102	! Print information about what observations are left after qc
1103	! -----------------------------------------------------------------------
1104
1105	! Update the total observation counter array
1106
1107	IF(lwp) THEN
1108	WRITE(numout,*)
1109	WRITE(numout,*) 'obs_pre_vel :'
1110	WRITE(numout,*) '~~~~~~~~~~~'
1111	WRITE(numout,*)
1112	WRITE(numout,*) ' Profiles outside time domain = ', &
1113	& iotdobsmpp
1114	WRITE(numout,*) ' Remaining profiles that failed grid search = ', &
1115	& igrdobsmpp
1116	WRITE(numout,*) ' Remaining U data outside space domain = ', &
1117	& iosduobsmpp
1118	WRITE(numout,*) ' Remaining U data at land points = ', &
1119	& ilanuobsmpp
1120	IF (ld_nea) THEN
1121	WRITE(numout,*) ' Remaining U data near land points (removed) = ',&
1122	& inlauobsmpp
1123	ELSE
1124	WRITE(numout,*) ' Remaining U data near land points (kept) = ',&
1125	& inlauobsmpp
1126	ENDIF
1127	WRITE(numout,*) ' U observation rejected since V rejected = ', &
1128	& iuvchku
1129	WRITE(numout,*) ' U data accepted = ', &
1130	& prodatqc%nvprotmpp(1)
1131	WRITE(numout,*) ' Remaining V data outside space domain = ', &
1132	& iosdvobsmpp
1133	WRITE(numout,*) ' Remaining V data at land points = ', &
1134	& ilanvobsmpp
1135	IF (ld_nea) THEN
1136	WRITE(numout,*) ' Remaining V data near land points (removed) = ',&
1137	& inlavobsmpp
1138	ELSE
1139	WRITE(numout,*) ' Remaining V data near land points (kept) = ',&
1140	& inlavobsmpp
1141	ENDIF
1142	WRITE(numout,*) ' V observation rejected since U rejected = ', &
1143	& iuvchkv
1144	WRITE(numout,*) ' V data accepted = ', &
1145	& prodatqc%nvprotmpp(2)
1146
1147	WRITE(numout,*)
1148	WRITE(numout,*) ' Number of observations per time step :'
1149	WRITE(numout,*)
1150	WRITE(numout,997)
1151	WRITE(numout,998)
1152	ENDIF
1153
1154	DO jobs = 1, prodatqc%nprof
1155	inrc = prodatqc%mstp(jobs) + 2 - nit000
1156	prodatqc%npstp(inrc) = prodatqc%npstp(inrc) + 1
1157	DO jvar = 1, prodatqc%nvar
1158	IF ( prodatqc%npvend(jobs,jvar) > 0 ) THEN
1159	prodatqc%nvstp(inrc,jvar) = prodatqc%nvstp(inrc,jvar) + &
1160	& ( prodatqc%npvend(jobs,jvar) - &
1161	& prodatqc%npvsta(jobs,jvar) + 1 )
1162	ENDIF
1163	END DO
1164	END DO
1165
1166
1167	CALL obs_mpp_sum_integers( prodatqc%npstp, prodatqc%npstpmpp, &
1168	& nitend - nit000 + 2 )
1169	DO jvar = 1, prodatqc%nvar
1170	CALL obs_mpp_sum_integers( prodatqc%nvstp(:,jvar), &
1171	& prodatqc%nvstpmpp(:,jvar), &
1172	& nitend - nit000 + 2 )
1173	END DO
1174
1175	IF ( lwp ) THEN
1176	DO jstp = nit000 - 1, nitend
1177	inrc = jstp - nit000 + 2
1178	WRITE(numout,999) jstp, prodatqc%npstpmpp(inrc), &
1179	& prodatqc%nvstpmpp(inrc,1), &
1180	& prodatqc%nvstpmpp(inrc,2)
1181	END DO
1182	ENDIF
1183
1184	997 FORMAT(10X,'Time step',5X,'Profiles',5X,'Zonal Comp.',5X,'Meridional Comp.')
1185	998 FORMAT(10X,'---------',5X,'--------',5X,'-----------',5X,'----------------')
1186	999 FORMAT(10X,I9,5X,I8,5X,I11,5X,I8)
1187
1188	END SUBROUTINE obs_pre_vel
1189
1190	SUBROUTINE obs_pre_logchl( logchldata, logchldatqc, ld_logchl, ld_nea )
1191	!!----------------------------------------------------------------------
1192	!! * ROUTINE obs_pre_logchl *
1193	!!
1194	!! ** Purpose : First level check and screening of logchl observations
1195	!!
1196	!! ** Method : First level check and screening of logchl observations
1197	!!
1198	!! ** Action :
1199	!!
1200	!! References :
1201	!!
1202	!! History :
1203	!!----------------------------------------------------------------------
1204	!! * Modules used
1205	USE domstp ! Domain: set the time-step
1206	USE par_oce ! Ocean parameters
1207	USE dom_oce, ONLY : & ! Geographical information
1208	& glamt, &
1209	& gphit, &
1210	& tmask
1211	!! * Arguments
1212	TYPE(obs_surf), INTENT(INOUT) :: logchldata ! Full set of logchl data
1213	TYPE(obs_surf), INTENT(INOUT) :: logchldatqc ! Subset of logchl data not failing screening
1214	LOGICAL :: ld_logchl ! Switch for logchl data
1215	LOGICAL :: ld_nea ! Switch for rejecting observation near land
1216	!! * Local declarations
1217	INTEGER :: iyea0 ! Initial date
1218	INTEGER :: imon0 ! - (year, month, day, hour, minute)
1219	INTEGER :: iday0
1220	INTEGER :: ihou0
1221	INTEGER :: imin0
1222	INTEGER :: icycle ! Current assimilation cycle
1223	! Counters for observations that
1224	INTEGER :: iotdobs ! - outside time domain
1225	INTEGER :: iosdsobs ! - outside space domain
1226	INTEGER :: ilansobs ! - within a model land cell
1227	INTEGER :: inlasobs ! - close to land
1228	INTEGER :: igrdobs ! - fail the grid search
1229	! Global counters for observations that
1230	INTEGER :: iotdobsmpp ! - outside time domain
1231	INTEGER :: iosdsobsmpp ! - outside space domain
1232	INTEGER :: ilansobsmpp ! - within a model land cell
1233	INTEGER :: inlasobsmpp ! - close to land
1234	INTEGER :: igrdobsmpp ! - fail the grid search
1235	LOGICAL, DIMENSION(:), ALLOCATABLE :: &
1236	& llvalid ! data selection
1237	INTEGER :: jobs ! Obs. loop variable
1238	INTEGER :: jstp ! Time loop variable
1239	INTEGER :: inrc ! Time index variable
1240
1241	IF (lwp) WRITE(numout,*)'obs_pre_logchl : Preparing the logchl observations...'
1242
1243	! Initial date initialization (year, month, day, hour, minute)
1244	iyea0 = ndate0 / 10000
1245	imon0 = ( ndate0 - iyea0 * 10000 ) / 100
1246	iday0 = ndate0 - iyea0 * 10000 - imon0 * 100
1247	ihou0 = 0
1248	imin0 = 0
1249
1250	icycle = no ! Assimilation cycle
1251
1252	! Diagnostics counters for various failures.
1253
1254	iotdobs = 0
1255	igrdobs = 0
1256	iosdsobs = 0
1257	ilansobs = 0
1258	inlasobs = 0
1259
1260	! -----------------------------------------------------------------------
1261	! Find time coordinate for logchl data
1262	! -----------------------------------------------------------------------
1263
1264	CALL obs_coo_tim( icycle, &
1265	& iyea0, imon0, iday0, ihou0, imin0, &
1266	& logchldata%nsurf, logchldata%nyea, logchldata%nmon, &
1267	& logchldata%nday, logchldata%nhou, logchldata%nmin, &
1268	& logchldata%nqc, logchldata%mstp, iotdobs )
1269	CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
1270	! -----------------------------------------------------------------------
1271	! Check for logchl data failing the grid search
1272	! -----------------------------------------------------------------------
1273
1274	CALL obs_coo_grd( logchldata%nsurf, logchldata%mi, logchldata%mj, &
1275	& logchldata%nqc, igrdobs )
1276	CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
1277
1278	! -----------------------------------------------------------------------
1279	! Check for land points.
1280	! -----------------------------------------------------------------------
1281
1282	CALL obs_coo_spc_2d( logchldata%nsurf, &
1283	& jpi, jpj, &
1284	& logchldata%mi, logchldata%mj, &
1285	& logchldata%rlam, logchldata%rphi, &
1286	& glamt, gphit, &
1287	& tmask(:,:,1), logchldata%nqc, &
1288	& iosdsobs, ilansobs, &
1289	& inlasobs, ld_nea )
1290
1291	CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
1292	CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
1293	CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
1294
1295	! -----------------------------------------------------------------------
1296	! Copy useful data from the logchldata data structure to
1297	! the logchldatqc data structure
1298	! -----------------------------------------------------------------------
1299
1300	! Allocate the selection arrays
1301
1302	ALLOCATE( llvalid(logchldata%nsurf) )
1303
1304	! We want all data which has qc flags <= 0
1305
1306	llvalid(:) = ( logchldata%nqc(:) <= 10 )
1307
1308	! The actual copying
1309
1310	CALL obs_surf_compress( logchldata, logchldatqc, .TRUE., numout, &
1311	& lvalid=llvalid )
1312
1313	! Dellocate the selection arrays
1314	DEALLOCATE( llvalid )
1315
1316	! -----------------------------------------------------------------------
1317	! Print information about what observations are left after qc
1318	! -----------------------------------------------------------------------
1319
1320	! Update the total observation counter array
1321
1322	IF(lwp) THEN
1323	WRITE(numout,*)
1324	WRITE(numout,*) 'obs_pre_logchl :'
1325	WRITE(numout,*) '~~~~~~~~~~~'
1326	WRITE(numout,*)
1327	WRITE(numout,*) ' logchl data outside time domain = ', &
1328	& iotdobsmpp
1329	WRITE(numout,*) ' Remaining logchl data that failed grid search = ', &
1330	& igrdobsmpp
1331	WRITE(numout,*) ' Remaining logchl data outside space domain = ', &
1332	& iosdsobsmpp
1333	WRITE(numout,*) ' Remaining logchl data at land points = ', &
1334	& ilansobsmpp
1335	IF (ld_nea) THEN
1336	WRITE(numout,*) ' Remaining logchl data near land points (removed) = ', &
1337	& inlasobsmpp
1338	ELSE
1339	WRITE(numout,*) ' Remaining logchl data near land points (kept) = ', &
1340	& inlasobsmpp
1341	ENDIF
1342	WRITE(numout,*) ' logchl data accepted = ', &
1343	& logchldatqc%nsurfmpp
1344
1345	WRITE(numout,*)
1346	WRITE(numout,*) ' Number of observations per time step :'
1347	WRITE(numout,*)
1348	WRITE(numout,1997)
1349	WRITE(numout,1998)
1350	ENDIF
1351
1352	DO jobs = 1, logchldatqc%nsurf
1353	inrc = logchldatqc%mstp(jobs) + 2 - nit000
1354	logchldatqc%nsstp(inrc) = logchldatqc%nsstp(inrc) + 1
1355	END DO
1356
1357	CALL obs_mpp_sum_integers( logchldatqc%nsstp, logchldatqc%nsstpmpp, &
1358	& nitend - nit000 + 2 )
1359
1360	IF ( lwp ) THEN
1361	DO jstp = nit000 - 1, nitend
1362	inrc = jstp - nit000 + 2
1363	WRITE(numout,1999) jstp, logchldatqc%nsstpmpp(inrc)
1364	END DO
1365	ENDIF
1366
1367	1997 FORMAT(10X,'Time step',5X,'logchl data')
1368	1998 FORMAT(10X,'---------',5X,'------------')
1369	1999 FORMAT(10X,I9,5X,I17)
1370
1371	END SUBROUTINE obs_pre_logchl
1372
1373	SUBROUTINE obs_coo_tim( kcycle, &
1374	& kyea0, kmon0, kday0, khou0, kmin0, &
1375	& kobsno, &
1376	& kobsyea, kobsmon, kobsday, kobshou, kobsmin, &
1377	& kobsqc, kobsstp, kotdobs )
1378	!!----------------------------------------------------------------------
1379	!! * ROUTINE obs_coo_tim *
1380	!!
1381	!! ** Purpose : Compute the number of time steps to the observation time.
1382	!!
1383	!! ** Method : For time coordinates ( yea_obs, mon_obs, day_obs,
1384	!! hou_obs, min_obs ), this routine locates the time step
1385	!! that is closest to this time.
1386	!!
1387	!! ** Action :
1388	!!
1389	!! References :
1390	!!
1391	!! History :
1392	!! ! 1997-07 (A. Weaver) Original
1393	!! ! 2006-08 (A. Weaver) NEMOVAR migration
1394	!! ! 2006-10 (A. Weaver) Cleanup
1395	!! ! 2007-01 (K. Mogensen) Rewritten with loop
1396	!! ! 2010-05 (D. Lea) Fix in leap year calculation for NEMO vn3.2
1397	!!----------------------------------------------------------------------
1398	!! * Modules used
1399	USE dom_oce, ONLY : & ! Geographical information
1400	& rdt
1401	USE phycst, ONLY : & ! Physical constants
1402	& rday, &
1403	& rmmss, &
1404	& rhhmm
1405	!! * Arguments
1406	INTEGER, INTENT(IN) :: kcycle ! Current cycle
1407	INTEGER, INTENT(IN) :: kyea0 ! Initial date coordinates
1408	INTEGER, INTENT(IN) :: kmon0
1409	INTEGER, INTENT(IN) :: kday0
1410	INTEGER, INTENT(IN) :: khou0
1411	INTEGER, INTENT(IN) :: kmin0
1412	INTEGER, INTENT(IN) :: kobsno ! Number of observations
1413	INTEGER, INTENT(INOUT) :: kotdobs ! Number of observations failing time check
1414	INTEGER, DIMENSION(kobsno), INTENT(IN ) :: &
1415	& kobsyea, & ! Observation time coordinates
1416	& kobsmon, &
1417	& kobsday, &
1418	& kobshou, &
1419	& kobsmin
1420	INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
1421	& kobsqc ! Quality control flag
1422	INTEGER, DIMENSION(kobsno), INTENT(OUT) :: &
1423	& kobsstp ! Number of time steps up to the
1424	! observation time
1425
1426	!! * Local declarations
1427	INTEGER :: jyea
1428	INTEGER :: jmon
1429	INTEGER :: jday
1430	INTEGER :: jobs
1431	INTEGER :: iyeastr
1432	INTEGER :: iyeaend
1433	INTEGER :: imonstr
1434	INTEGER :: imonend
1435	INTEGER :: idaystr
1436	INTEGER :: idayend
1437	INTEGER :: iskip
1438	INTEGER :: idaystp
1439	REAL(KIND=wp) :: zminstp
1440	REAL(KIND=wp) :: zhoustp
1441	REAL(KIND=wp) :: zobsstp
1442	INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year
1443
1444	!-----------------------------------------------------------------------
1445	! Initialization
1446	!-----------------------------------------------------------------------
1447
1448	! Intialize the number of time steps per day
1449	idaystp = NINT( rday / rdt )
1450
1451	!---------------------------------------------------------------------
1452	! Locate the model time coordinates for interpolation
1453	!---------------------------------------------------------------------
1454
1455	DO jobs = 1, kobsno
1456
1457	! Initialize the time step counter
1458	kobsstp(jobs) = nit000 - 1
1459
1460	! Flag if observation date is less than the initial date
1461
1462	IF ( ( kobsyea(jobs) < kyea0 ) &
1463	& .OR. ( ( kobsyea(jobs) == kyea0 ) &
1464	& .AND. ( kobsmon(jobs) < kmon0 ) ) &
1465	& .OR. ( ( kobsyea(jobs) == kyea0 ) &
1466	& .AND. ( kobsmon(jobs) == kmon0 ) &
1467	& .AND. ( kobsday(jobs) < kday0 ) ) &
1468	& .OR. ( ( kobsyea(jobs) == kyea0 ) &
1469	& .AND. ( kobsmon(jobs) == kmon0 ) &
1470	& .AND. ( kobsday(jobs) == kday0 ) &
1471	& .AND. ( kobshou(jobs) < khou0 ) ) &
1472	& .OR. ( ( kobsyea(jobs) == kyea0 ) &
1473	& .AND. ( kobsmon(jobs) == kmon0 ) &
1474	& .AND. ( kobsday(jobs) == kday0 ) &
1475	& .AND. ( kobshou(jobs) == khou0 ) &
1476	& .AND. ( kobsmin(jobs) <= kmin0 ) ) ) THEN
1477	kobsstp(jobs) = -1
1478	kobsqc(jobs) = kobsqc(jobs) + 11
1479	kotdobs = kotdobs + 1
1480	CYCLE
1481	ENDIF
1482
1483	! Compute the number of time steps to the observation day
1484	iyeastr = kyea0
1485	iyeaend = kobsyea(jobs)
1486
1487	!---------------------------------------------------------------------
1488	! Year loop
1489	!---------------------------------------------------------------------
1490	DO jyea = iyeastr, iyeaend
1491
1492	CALL calc_month_len( jyea, imonth_len )
1493
1494	imonstr = 1
1495	IF ( jyea == kyea0 ) imonstr = kmon0
1496	imonend = 12
1497	IF ( jyea == kobsyea(jobs) ) imonend = kobsmon(jobs)
1498
1499	! Month loop
1500	DO jmon = imonstr, imonend
1501
1502	idaystr = 1
1503	IF ( ( jmon == kmon0 ) &
1504	& .AND. ( jyea == kyea0 ) ) idaystr = kday0
1505	idayend = imonth_len(jmon)
1506	IF ( ( jmon == kobsmon(jobs) ) &
1507	& .AND. ( jyea == kobsyea(jobs) ) ) idayend = kobsday(jobs) - 1
1508
1509	! Day loop
1510	DO jday = idaystr, idayend
1511	kobsstp(jobs) = kobsstp(jobs) + idaystp
1512	END DO
1513
1514	END DO
1515
1516	END DO
1517
1518	! Add in the number of time steps to the observation minute
1519	zminstp = rmmss / rdt
1520	zhoustp = rhhmm * zminstp
1521
1522	zobsstp = REAL( kobsmin(jobs) - kmin0, KIND=wp ) * zminstp &
1523	& + REAL( kobshou(jobs) - khou0, KIND=wp ) * zhoustp
1524	kobsstp(jobs) = kobsstp(jobs) + NINT( zobsstp )
1525
1526	! Flag if observation step outside the time window
1527	IF ( ( kobsstp(jobs) < ( nit000 - 1 ) ) &
1528	& .OR.( kobsstp(jobs) > nitend ) ) THEN
1529	kobsqc(jobs) = kobsqc(jobs) + 12
1530	kotdobs = kotdobs + 1
1531	CYCLE
1532	ENDIF
1533
1534	END DO
1535
1536	END SUBROUTINE obs_coo_tim
1537
1538	SUBROUTINE calc_month_len( iyear, imonth_len )
1539	!!----------------------------------------------------------------------
1540	!! * ROUTINE calc_month_len *
1541	!!
1542	!! ** Purpose : Compute the number of days in a months given a year.
1543	!!
1544	!! ** Method :
1545	!!
1546	!! ** Action :
1547	!!
1548	!! History :
1549	!! ! 10-05 (D. Lea) New routine based on day_init
1550	!!----------------------------------------------------------------------
1551
1552	INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year
1553	INTEGER :: iyear !: year
1554
1555	! length of the month of the current year (from nleapy, read in namelist)
1556	IF ( nleapy < 2 ) THEN
1557	imonth_len(:) = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)
1558	IF ( nleapy == 1 ) THEN ! we are using calendar with leap years
1559	IF ( MOD(iyear, 4) == 0 .AND. ( MOD(iyear, 400) == 0 .OR. MOD(iyear, 100) /= 0 ) ) THEN
1560	imonth_len(2) = 29
1561	ENDIF
1562	ENDIF
1563	ELSE
1564	imonth_len(:) = nleapy ! all months with nleapy days per year
1565	ENDIF
1566
1567	END SUBROUTINE
1568
1569	SUBROUTINE obs_coo_tim_prof( kcycle, &
1570	& kyea0, kmon0, kday0, khou0, kmin0, &
1571	& kobsno, &
1572	& kobsyea, kobsmon, kobsday, kobshou, kobsmin, &
1573	& ktyp, kobsqc, kobsstp, kotdobs, kdailyavtypes, &
1574	& ld_dailyav )
1575	!!----------------------------------------------------------------------
1576	!! * ROUTINE obs_coo_tim *
1577	!!
1578	!! ** Purpose : Compute the number of time steps to the observation time.
1579	!!
1580	!! ** Method : For time coordinates ( yea_obs, mon_obs, day_obs,
1581	!! hou_obs, min_obs ), this routine locates the time step
1582	!! that is closest to this time.
1583	!!
1584	!! ** Action :
1585	!!
1586	!! References :
1587	!!
1588	!! History :
1589	!! ! 1997-07 (A. Weaver) Original
1590	!! ! 2006-08 (A. Weaver) NEMOVAR migration
1591	!! ! 2006-10 (A. Weaver) Cleanup
1592	!! ! 2007-01 (K. Mogensen) Rewritten with loop
1593	!!----------------------------------------------------------------------
1594	!! * Modules used
1595	!! * Arguments
1596	INTEGER, INTENT(IN) :: kcycle ! Current cycle
1597	INTEGER, INTENT(IN) :: kyea0 ! Initial date coordinates
1598	INTEGER, INTENT(IN) :: kmon0
1599	INTEGER, INTENT(IN) :: kday0
1600	INTEGER, INTENT(IN) :: khou0
1601	INTEGER, INTENT(IN) :: kmin0
1602	INTEGER, INTENT(IN) :: kobsno ! Number of observations
1603	INTEGER, INTENT(INOUT) :: kotdobs ! Number of observations failing time check
1604	INTEGER, DIMENSION(kobsno), INTENT(IN ) :: &
1605	& kobsyea, & ! Observation time coordinates
1606	& kobsmon, &
1607	& kobsday, &
1608	& kobshou, &
1609	& kobsmin, &
1610	& ktyp ! Observation type.
1611	INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
1612	& kobsqc ! Quality control flag
1613	INTEGER, DIMENSION(kobsno), INTENT(OUT) :: &
1614	& kobsstp ! Number of time steps up to the
1615	! observation time
1616	INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
1617	& kdailyavtypes ! Types for daily averages
1618	LOGICAL, OPTIONAL :: ld_dailyav ! All types are daily averages
1619	!! * Local declarations
1620	INTEGER :: jobs
1621
1622	!-----------------------------------------------------------------------
1623	! Call standard obs_coo_tim
1624	!-----------------------------------------------------------------------
1625
1626	CALL obs_coo_tim( kcycle, &
1627	& kyea0, kmon0, kday0, khou0, kmin0, &
1628	& kobsno, &
1629	& kobsyea, kobsmon, kobsday, kobshou, kobsmin, &
1630	& kobsqc, kobsstp, kotdobs )
1631
1632	!------------------------------------------------------------------------
1633	! Always reject daily averaged data (e.g. MRB data (820)) at initial time
1634	!------------------------------------------------------------------------
1635
1636	IF ( PRESENT(kdailyavtypes) ) THEN
1637	DO jobs = 1, kobsno
1638
1639	IF ( kobsqc(jobs) <= 10 ) THEN
1640
1641	IF ( ( kobsstp(jobs) == (nit000 - 1) ).AND.&
1642	& ( ANY (kdailyavtypes(:) == ktyp(jobs)) ) ) THEN
1643	kobsqc(jobs) = kobsqc(jobs) + 14
1644	kotdobs = kotdobs + 1
1645	CYCLE
1646	ENDIF
1647
1648	ENDIF
1649	END DO
1650	ENDIF
1651
1652	!------------------------------------------------------------------------
1653	! If ld_dailyav is set then all data assumed to be daily averaged
1654	!------------------------------------------------------------------------
1655
1656	IF ( PRESENT( ld_dailyav) ) THEN
1657	IF (ld_dailyav) THEN
1658	DO jobs = 1, kobsno
1659
1660	IF ( kobsqc(jobs) <= 10 ) THEN
1661
1662	IF ( kobsstp(jobs) == (nit000 - 1) ) THEN
1663	kobsqc(jobs) = kobsqc(jobs) + 14
1664	kotdobs = kotdobs + 1
1665	CYCLE
1666	ENDIF
1667
1668	ENDIF
1669	END DO
1670	ENDIF
1671	ENDIF
1672
1673	END SUBROUTINE obs_coo_tim_prof
1674
1675	SUBROUTINE obs_coo_grd( kobsno, kobsi, kobsj, kobsqc, kgrdobs )
1676	!!----------------------------------------------------------------------
1677	!! * ROUTINE obs_coo_grd *
1678	!!
1679	!! ** Purpose : Verify that the grid search has not failed
1680	!!
1681	!! ** Method : The previously computed i,j indeces are checked
1682	!!
1683	!! ** Action :
1684	!!
1685	!! References :
1686	!!
1687	!! History :
1688	!! ! 2007-01 (K. Mogensen) Original
1689	!!----------------------------------------------------------------------
1690	!! * Arguments
1691	INTEGER, INTENT(IN) :: kobsno ! Number of observations
1692	INTEGER, DIMENSION(kobsno), INTENT(IN ) :: &
1693	& kobsi, & ! i,j indeces previously computed
1694	& kobsj
1695	INTEGER, INTENT(INOUT) :: kgrdobs ! Number of observations failing the check
1696	INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
1697	& kobsqc ! Quality control flag
1698
1699	!! * Local declarations
1700	INTEGER :: jobs ! Loop variable
1701
1702	! Flag if the grid search failed
1703
1704	DO jobs = 1, kobsno
1705	IF ( ( kobsi(jobs) <= 0 ) .AND. ( kobsj(jobs) <= 0 ) ) THEN
1706	kobsqc(jobs) = kobsqc(jobs) + 18
1707	kgrdobs = kgrdobs + 1
1708	ENDIF
1709	END DO
1710
1711	END SUBROUTINE obs_coo_grd
1712
1713	SUBROUTINE obs_coo_spc_2d( kobsno, kpi, kpj, &
1714	& kobsi, kobsj, pobslam, pobsphi, &
1715	& plam, pphi, pmask, &
1716	& kobsqc, kosdobs, klanobs, &
1717	& knlaobs,ld_nea )
1718	!!----------------------------------------------------------------------
1719	!! * ROUTINE obs_coo_spc_2d *
1720	!!
1721	!! ** Purpose : Check for points outside the domain and land points
1722	!!
1723	!! ** Method : Remove the observations that are outside the model space
1724	!! and time domain or located within model land cells.
1725	!!
1726	!! ** Action :
1727	!!
1728	!! History :
1729	!! ! 2007-03 (A. Weaver, K. Mogensen) Original
1730	!! ! 2007-06 (K. Mogensen et al) Reject obs. near land.
1731	!!----------------------------------------------------------------------
1732	!! * Modules used
1733
1734	!! * Arguments
1735	INTEGER, INTENT(IN) :: kobsno ! Total number of observations
1736	INTEGER, INTENT(IN) :: kpi ! Number of grid points in (i,j)
1737	INTEGER, INTENT(IN) :: kpj
1738	INTEGER, DIMENSION(kobsno), INTENT(IN) :: &
1739	& kobsi, & ! Observation (i,j) coordinates
1740	& kobsj
1741	REAL(KIND=wp), DIMENSION(kobsno), INTENT(IN) :: &
1742	& pobslam, & ! Observation (lon,lat) coordinates
1743	& pobsphi
1744	REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) :: &
1745	& plam, pphi ! Model (lon,lat) coordinates
1746	REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) :: &
1747	& pmask ! Land mask array
1748	INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
1749	& kobsqc ! Observation quality control
1750	INTEGER, INTENT(INOUT) :: kosdobs ! Observations outside space domain
1751	INTEGER, INTENT(INOUT) :: klanobs ! Observations within a model land cell
1752	INTEGER, INTENT(INOUT) :: knlaobs ! Observations near land
1753	LOGICAL, INTENT(IN) :: ld_nea ! Flag observations near land
1754	!! * Local declarations
1755	REAL(KIND=wp), DIMENSION(2,2,kobsno) :: &
1756	& zgmsk ! Grid mask
1757	REAL(KIND=wp), DIMENSION(2,2,kobsno) :: &
1758	& zglam, & ! Model longitude at grid points
1759	& zgphi ! Model latitude at grid points
1760	INTEGER, DIMENSION(2,2,kobsno) :: &
1761	& igrdi, & ! Grid i,j
1762	& igrdj
1763	LOGICAL :: lgridobs ! Is observation on a model grid point.
1764	INTEGER :: iig, ijg ! i,j of observation on model grid point.
1765	INTEGER :: jobs, ji, jj
1766
1767	! Get grid point indices
1768
1769	DO jobs = 1, kobsno
1770
1771	! For invalid points use 2,2
1772
1773	IF ( kobsqc(jobs) >= 10 ) THEN
1774
1775	igrdi(1,1,jobs) = 1
1776	igrdj(1,1,jobs) = 1
1777	igrdi(1,2,jobs) = 1
1778	igrdj(1,2,jobs) = 2
1779	igrdi(2,1,jobs) = 2
1780	igrdj(2,1,jobs) = 1
1781	igrdi(2,2,jobs) = 2
1782	igrdj(2,2,jobs) = 2
1783
1784	ELSE
1785
1786	igrdi(1,1,jobs) = kobsi(jobs)-1
1787	igrdj(1,1,jobs) = kobsj(jobs)-1
1788	igrdi(1,2,jobs) = kobsi(jobs)-1
1789	igrdj(1,2,jobs) = kobsj(jobs)
1790	igrdi(2,1,jobs) = kobsi(jobs)
1791	igrdj(2,1,jobs) = kobsj(jobs)-1
1792	igrdi(2,2,jobs) = kobsi(jobs)
1793	igrdj(2,2,jobs) = kobsj(jobs)
1794
1795	ENDIF
1796
1797	END DO
1798
1799	CALL obs_int_comm_2d( 2, 2, kobsno, igrdi, igrdj, pmask, zgmsk )
1800	CALL obs_int_comm_2d( 2, 2, kobsno, igrdi, igrdj, plam, zglam )
1801	CALL obs_int_comm_2d( 2, 2, kobsno, igrdi, igrdj, pphi, zgphi )
1802
1803	DO jobs = 1, kobsno
1804
1805	! Skip bad observations
1806	IF ( kobsqc(jobs) >= 10 ) CYCLE
1807
1808	! Flag if the observation falls outside the model spatial domain
1809	IF ( ( pobslam(jobs) < -180. ) &
1810	& .OR. ( pobslam(jobs) > 180. ) &
1811	& .OR. ( pobsphi(jobs) < -90. ) &
1812	& .OR. ( pobsphi(jobs) > 90. ) ) THEN
1813	kobsqc(jobs) = kobsqc(jobs) + 11
1814	kosdobs = kosdobs + 1
1815	CYCLE
1816	ENDIF
1817
1818	! Flag if the observation falls with a model land cell
1819	IF ( SUM( zgmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN
1820	kobsqc(jobs) = kobsqc(jobs) + 12
1821	klanobs = klanobs + 1
1822	CYCLE
1823	ENDIF
1824
1825	! Check if this observation is on a grid point
1826
1827	lgridobs = .FALSE.
1828	iig = -1
1829	ijg = -1
1830	DO jj = 1, 2
1831	DO ji = 1, 2
1832	IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) &
1833	& .AND. &
1834	& ( ABS( zglam(ji,jj,jobs) - pobslam(jobs) ) < 1.0e-6_wp ) &
1835	& ) THEN
1836	lgridobs = .TRUE.
1837	iig = ji
1838	ijg = jj
1839	ENDIF
1840	END DO
1841	END DO
1842
1843	! For observations on the grid reject them if their are at
1844	! a masked point
1845
1846	IF (lgridobs) THEN
1847	IF (zgmsk(iig,ijg,jobs) == 0.0_wp ) THEN
1848	kobsqc(jobs) = kobsqc(jobs) + 12
1849	klanobs = klanobs + 1
1850	CYCLE
1851	ENDIF
1852	ENDIF
1853
1854	! Flag if the observation falls is close to land
1855	IF ( MINVAL( zgmsk(1:2,1:2,jobs) ) == 0.0_wp) THEN
1856	IF (ld_nea) kobsqc(jobs) = kobsqc(jobs) + 14
1857	knlaobs = knlaobs + 1
1858	CYCLE
1859	ENDIF
1860
1861	END DO
1862
1863	END SUBROUTINE obs_coo_spc_2d
1864
1865	SUBROUTINE obs_coo_spc_3d( kprofno, kobsno, kpstart, kpend, &
1866	& kpi, kpj, kpk, &
1867	& kobsi, kobsj, kobsk, &
1868	& pobslam, pobsphi, pobsdep, &
1869	& plam, pphi, pdep, pmask, &
1870	& kpobsqc, kobsqc, kosdobs, &
1871	& klanobs, knlaobs, ld_nea )
1872	!!----------------------------------------------------------------------
1873	!! * ROUTINE obs_coo_spc_3d *
1874	!!
1875	!! ** Purpose : Check for points outside the domain and land points
1876	!! Reset depth of observation above highest model level
1877	!! to the value of highest model level
1878	!!
1879	!! ** Method : Remove the observations that are outside the model space
1880	!! and time domain or located within model land cells.
1881	!!
1882	!! NB. T and S profile observations lying between the ocean
1883	!! surface and the depth of the first model T point are
1884	!! assigned a depth equal to that of the first model T pt.
1885	!!
1886	!! ** Action :
1887	!!
1888	!! History :
1889	!! ! 2007-01 (K. Mogensen) Rewrite of parts of obs_scr
1890	!! ! 2007-06 (K. Mogensen et al) Reject obs. near land.
1891	!!----------------------------------------------------------------------
1892	!! * Modules used
1893	USE dom_oce, ONLY : & ! Geographical information
1894	& gdepw_1d, &
1895	& gdepw_0, &
1896	#if defined key_vvl
1897	& gdepw_n, &
1898	& gdept_n, &
1899	#endif
1900	& ln_zco, &
1901	& ln_zps, &
1902	& lk_vvl
1903
1904	!! * Arguments
1905	INTEGER, INTENT(IN) :: kprofno ! Number of profiles
1906	INTEGER, INTENT(IN) :: kobsno ! Total number of observations
1907	INTEGER, INTENT(IN) :: kpi ! Number of grid points in (i,j,k)
1908	INTEGER, INTENT(IN) :: kpj
1909	INTEGER, INTENT(IN) :: kpk
1910	INTEGER, DIMENSION(kprofno), INTENT(IN) :: &
1911	& kpstart, & ! Start of individual profiles
1912	& kpend ! End of individual profiles
1913	INTEGER, DIMENSION(kprofno), INTENT(IN) :: &
1914	& kobsi, & ! Observation (i,j) coordinates
1915	& kobsj
1916	INTEGER, DIMENSION(kobsno), INTENT(IN) :: &
1917	& kobsk ! Observation k coordinate
1918	REAL(KIND=wp), DIMENSION(kprofno), INTENT(IN) :: &
1919	& pobslam, & ! Observation (lon,lat) coordinates
1920	& pobsphi
1921	REAL(KIND=wp), DIMENSION(kobsno), INTENT(INOUT) :: &
1922	& pobsdep ! Observation depths
1923	REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) :: &
1924	& plam, pphi ! Model (lon,lat) coordinates
1925	REAL(KIND=wp), DIMENSION(kpk), INTENT(IN) :: &
1926	& pdep ! Model depth coordinates
1927	REAL(KIND=wp), DIMENSION(kpi,kpj,kpk), INTENT(IN) :: &
1928	& pmask ! Land mask array
1929	INTEGER, DIMENSION(kprofno), INTENT(INOUT) :: &
1930	& kpobsqc ! Profile quality control
1931	INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
1932	& kobsqc ! Observation quality control
1933	INTEGER, INTENT(INOUT) :: kosdobs ! Observations outside space domain
1934	INTEGER, INTENT(INOUT) :: klanobs ! Observations within a model land cell
1935	INTEGER, INTENT(INOUT) :: knlaobs ! Observations near land
1936	LOGICAL, INTENT(IN) :: ld_nea ! Flag observations near land
1937	!! * Local declarations
1938	REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: &
1939	& zgmsk ! Grid mask
1940	REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: &
1941	& zgdepw
1942	REAL(KIND=wp), DIMENSION(2,2,kprofno) :: &
1943	& zglam, & ! Model longitude at grid points
1944	& zgphi ! Model latitude at grid points
1945	INTEGER, DIMENSION(2,2,kprofno) :: &
1946	& igrdi, & ! Grid i,j
1947	& igrdj
1948	LOGICAL :: lgridobs ! Is observation on a model grid point.
1949	LOGICAL :: ll_next_to_land ! Is a profile next to land
1950	INTEGER :: iig, ijg ! i,j of observation on model grid point.
1951	INTEGER :: jobs, jobsp, jk, ji, jj
1952
1953	! Get grid point indices
1954
1955	DO jobs = 1, kprofno
1956
1957	! For invalid points use 2,2
1958
1959	IF ( kpobsqc(jobs) >= 10 ) THEN
1960
1961	igrdi(1,1,jobs) = 1
1962	igrdj(1,1,jobs) = 1
1963	igrdi(1,2,jobs) = 1
1964	igrdj(1,2,jobs) = 2
1965	igrdi(2,1,jobs) = 2
1966	igrdj(2,1,jobs) = 1
1967	igrdi(2,2,jobs) = 2
1968	igrdj(2,2,jobs) = 2
1969
1970	ELSE
1971
1972	igrdi(1,1,jobs) = kobsi(jobs)-1
1973	igrdj(1,1,jobs) = kobsj(jobs)-1
1974	igrdi(1,2,jobs) = kobsi(jobs)-1
1975	igrdj(1,2,jobs) = kobsj(jobs)
1976	igrdi(2,1,jobs) = kobsi(jobs)
1977	igrdj(2,1,jobs) = kobsj(jobs)-1
1978	igrdi(2,2,jobs) = kobsi(jobs)
1979	igrdj(2,2,jobs) = kobsj(jobs)
1980
1981	ENDIF
1982
1983	END DO
1984
1985	CALL obs_int_comm_3d( 2, 2, kprofno, kpk, igrdi, igrdj, pmask, zgmsk )
1986	CALL obs_int_comm_2d( 2, 2, kprofno, igrdi, igrdj, plam, zglam )
1987	CALL obs_int_comm_2d( 2, 2, kprofno, igrdi, igrdj, pphi, zgphi )
1988	! Need to know the bathy depth for each observation for sco
1989	CALL obs_int_comm_3d( 2, 2, kprofno, kpk, igrdi, igrdj, fsdepw(:,:,:), &
1990	& zgdepw )
1991
1992	DO jobs = 1, kprofno
1993
1994	! Skip bad profiles
1995	IF ( kpobsqc(jobs) >= 10 ) CYCLE
1996
1997	! Check if this observation is on a grid point
1998
1999	lgridobs = .FALSE.
2000	iig = -1
2001	ijg = -1
2002	DO jj = 1, 2
2003	DO ji = 1, 2
2004	IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) &
2005	& .AND. &
2006	& ( ABS( zglam(ji,jj,jobs) - pobslam(jobs) ) < 1.0e-6_wp ) &
2007	& ) THEN
2008	lgridobs = .TRUE.
2009	iig = ji
2010	ijg = jj
2011	ENDIF
2012	END DO
2013	END DO
2014
2015	! Check if next to land
2016	IF ( ANY( zgmsk(1:2,1:2,1,jobs) == 0.0_wp ) ) THEN
2017	ll_next_to_land=.TRUE.
2018	ELSE
2019	ll_next_to_land=.FALSE.
2020	ENDIF
2021
2022	! Reject observations
2023
2024	DO jobsp = kpstart(jobs), kpend(jobs)
2025
2026	! Flag if the observation falls outside the model spatial domain
2027	IF ( ( pobslam(jobs) < -180. ) &
2028	& .OR. ( pobslam(jobs) > 180. ) &
2029	& .OR. ( pobsphi(jobs) < -90. ) &
2030	& .OR. ( pobsphi(jobs) > 90. ) &
2031	& .OR. ( pobsdep(jobsp) < 0.0 ) &
2032	& .OR. ( pobsdep(jobsp) > gdepw_1d(kpk)) ) THEN
2033	kobsqc(jobsp) = kobsqc(jobsp) + 11
2034	kosdobs = kosdobs + 1
2035	CYCLE
2036	ENDIF
2037
2038	! To check if an observations falls within land there are two cases:
2039	! 1: z-coordibnates, where the check uses the mask
2040	! 2: terrain following (eg s-coordinates),
2041	! where we use the depth of the bottom cell to mask observations
2042
2043	IF( (.NOT. lk_vvl) .AND. ( ln_zps .OR. ln_zco ) ) THEN !(CASE 1)
2044
2045	! Flag if the observation falls with a model land cell
2046	IF ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) &
2047	& == 0.0_wp ) THEN
2048	kobsqc(jobsp) = kobsqc(jobsp) + 12
2049	klanobs = klanobs + 1
2050	CYCLE
2051	ENDIF
2052
2053	! Flag if the observation is close to land
2054	IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == &
2055	& 0.0_wp) THEN
2056	knlaobs = knlaobs + 1
2057	IF (ld_nea) THEN
2058	kobsqc(jobsp) = kobsqc(jobsp) + 14
2059	ENDIF
2060	ENDIF
2061
2062	ELSE ! Case 2
2063
2064	! Flag if the observation is deeper than the bathymetry
2065	! Or if it is within the mask
2066	IF ( ANY( zgdepw(1:2,1:2,kpk,jobs) < pobsdep(jobsp) ) &
2067	& .OR. &
2068	& ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) &
2069	& == 0.0_wp) ) THEN
2070	kobsqc(jobsp) = kobsqc(jobsp) + 12
2071	klanobs = klanobs + 1
2072	CYCLE
2073	ENDIF
2074
2075	! Flag if the observation is close to land
2076	IF ( ll_next_to_land ) THEN
2077	knlaobs = knlaobs + 1
2078	IF (ld_nea) THEN
2079	kobsqc(jobsp) = kobsqc(jobsp) + 14
2080	ENDIF
2081	ENDIF
2082
2083	ENDIF
2084
2085	! For observations on the grid reject them if their are at
2086	! a masked point
2087
2088	IF (lgridobs) THEN
2089	IF (zgmsk(iig,ijg,kobsk(jobsp)-1,jobs) == 0.0_wp ) THEN
2090	kobsqc(jobsp) = kobsqc(jobsp) + 12
2091	klanobs = klanobs + 1
2092	CYCLE
2093	ENDIF
2094	ENDIF
2095
2096	! Flag if the observation falls is close to land
2097	IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == &
2098	& 0.0_wp) THEN
2099	IF (ld_nea) kobsqc(jobsp) = kobsqc(jobsp) + 14
2100	knlaobs = knlaobs + 1
2101	ENDIF
2102
2103	! Set observation depth equal to that of the first model depth
2104	IF ( pobsdep(jobsp) <= pdep(1) ) THEN
2105	pobsdep(jobsp) = pdep(1)
2106	ENDIF
2107
2108	END DO
2109	END DO
2110
2111	END SUBROUTINE obs_coo_spc_3d
2112
2113	SUBROUTINE obs_pro_rej( profdata )
2114	!!----------------------------------------------------------------------
2115	!! * ROUTINE obs_pro_rej *
2116	!!
2117	!! ** Purpose : Reject all data within a rejected profile
2118	!!
2119	!! ** Method :
2120	!!
2121	!! ** Action :
2122	!!
2123	!! References :
2124	!!
2125	!! History :
2126	!! ! 2007-10 (K. Mogensen) Original code
2127	!!----------------------------------------------------------------------
2128	!! * Modules used
2129	!! * Arguments
2130	TYPE(obs_prof), INTENT(INOUT) :: profdata ! Profile data
2131	!! * Local declarations
2132	INTEGER :: jprof
2133	INTEGER :: jvar
2134	INTEGER :: jobs
2135
2136	! Loop over profiles
2137
2138	DO jprof = 1, profdata%nprof
2139
2140	IF ( profdata%nqc(jprof) > 10 ) THEN
2141
2142	DO jvar = 1, profdata%nvar
2143
2144	DO jobs = profdata%npvsta(jprof,jvar), &
2145	& profdata%npvend(jprof,jvar)
2146
2147	profdata%var(jvar)%nvqc(jobs) = &
2148	& profdata%var(jvar)%nvqc(jobs) + 26
2149
2150	END DO
2151
2152	END DO
2153
2154	ENDIF
2155
2156	END DO
2157
2158	END SUBROUTINE obs_pro_rej
2159
2160	SUBROUTINE obs_uv_rej( profdata, knumu, knumv )
2161	!!----------------------------------------------------------------------
2162	!! * ROUTINE obs_uv_rej *
2163	!!
2164	!! ** Purpose : Reject u if v is rejected and vice versa
2165	!!
2166	!! ** Method :
2167	!!
2168	!! ** Action :
2169	!!
2170	!! References :
2171	!!
2172	!! History :
2173	!! ! 2009-2 (K. Mogensen) Original code
2174	!!----------------------------------------------------------------------
2175	!! * Modules used
2176	!! * Arguments
2177	TYPE(obs_prof), INTENT(INOUT) :: profdata ! Profile data
2178	INTEGER, INTENT(INOUT) :: knumu ! Number of u rejected
2179	INTEGER, INTENT(INOUT) :: knumv ! Number of v rejected
2180	!! * Local declarations
2181	INTEGER :: jprof
2182	INTEGER :: jvar
2183	INTEGER :: jobs
2184
2185	! Loop over profiles
2186
2187	DO jprof = 1, profdata%nprof
2188
2189	IF ( ( profdata%npvsta(jprof,1) /= profdata%npvsta(jprof,2) ) .OR. &
2190	& ( profdata%npvend(jprof,1) /= profdata%npvend(jprof,2) ) ) THEN
2191
2192	CALL ctl_stop('U,V profiles inconsistent in obs_uv_rej')
2193	RETURN
2194
2195	ENDIF
2196
2197	DO jobs = profdata%npvsta(jprof,1), profdata%npvend(jprof,1)
2198
2199	IF ( ( profdata%var(1)%nvqc(jobs) > 10 ) .AND. &
2200	& ( profdata%var(2)%nvqc(jobs) <= 10) ) THEN
2201	profdata%var(2)%nvqc(jobs) = profdata%var(2)%nvqc(jobs) + 42
2202	knumv = knumv + 1
2203	ENDIF
2204	IF ( ( profdata%var(2)%nvqc(jobs) > 10 ) .AND. &
2205	& ( profdata%var(1)%nvqc(jobs) <= 10) ) THEN
2206	profdata%var(1)%nvqc(jobs) = profdata%var(1)%nvqc(jobs) + 42
2207	knumu = knumu + 1
2208	ENDIF
2209
2210	END DO
2211
2212	END DO
2213
2214	END SUBROUTINE obs_uv_rej
2215
2216	END MODULE obs_prep

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

Download in other formats: