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

Context Navigation

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

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

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

Download in other formats: