Context Navigation

obs_prep.F90 @ 2287

Last change on this file since 2287 was 2287, checked in by smasson, 14 years ago
update licence of all NEMO files...
Property svn:keywords set to `Id`
File size: 80.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_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
30	IMPLICIT NONE
31
32	!! * Routine accessibility
33	PRIVATE
34
35	PUBLIC &
36	& obs_pre_pro, & ! First level check and screening of profiles
37	& obs_pre_sla, & ! First level check and screening of SLA data
38	& obs_pre_sst, & ! First level check and screening of SLA data
39	& obs_pre_seaice, & ! First level check and screening of sea ice data
40	& obs_pre_vel, & ! First level check and screening of velocity profiles
41	& calc_month_len ! Calculate the number of days in the months of a year
42
43	!!----------------------------------------------------------------------
44	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
45	!! $Id$
46	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
47	!!----------------------------------------------------------------------
48
49	CONTAINS
50
51	SUBROUTINE obs_pre_pro( profdata, prodatqc, ld_t3d, ld_s3d, ld_nea, &
52	& kdailyavtypes )
53	!!----------------------------------------------------------------------
54	!! * ROUTINE obs_pre_pro *
55	!!
56	!! ** Purpose : First level check and screening of T and S profiles
57	!!
58	!! ** Method : First level check and screening of T and S profiles
59	!!
60	!! ** Action :
61	!!
62	!! References :
63	!!
64	!! History :
65	!! ! 2007-01 (K. Mogensen) Merge of obs_pre_t3d and obs_pre_s3d
66	!! ! 2007-03 (K. Mogensen) General handling of profiles
67	!! ! 2007-06 (K. Mogensen et al) Reject obs. near land.
68	!!----------------------------------------------------------------------
69	!! * Modules used
70	USE domstp ! Domain: set the time-step
71	USE par_oce ! Ocean parameters
72	USE dom_oce, ONLY : & ! Geographical information
73	& glamt, &
74	& gphit, &
75	& gdept_0, &
76	& tmask, &
77	& nproc
78	!! * Arguments
79	TYPE(obs_prof), INTENT(INOUT) :: profdata ! Full set of profile data
80	TYPE(obs_prof), INTENT(INOUT) :: prodatqc ! Subset of profile data not failing screening
81	LOGICAL, INTENT(IN) :: ld_t3d ! Switch for temperature
82	LOGICAL, INTENT(IN) :: ld_s3d ! Switch for salinity
83	LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land
84	INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
85	& kdailyavtypes! Types for daily averages
86	!! * Local declarations
87	INTEGER :: iyea0 ! Initial date
88	INTEGER :: imon0 ! - (year, month, day, hour, minute)
89	INTEGER :: iday0
90	INTEGER :: ihou0
91	INTEGER :: imin0
92	INTEGER :: icycle ! Current assimilation cycle
93	! Counters for observations that
94	INTEGER :: iotdobs ! - outside time domain
95	INTEGER :: iosdtobs ! - outside space domain (temperature)
96	INTEGER :: iosdsobs ! - outside space domain (salinity)
97	INTEGER :: ilantobs ! - within a model land cell (temperature)
98	INTEGER :: ilansobs ! - within a model land cell (salinity)
99	INTEGER :: inlatobs ! - close to land (temperature)
100	INTEGER :: inlasobs ! - close to land (salinity)
101	INTEGER :: igrdobs ! - fail the grid search
102	! Global counters for observations that
103	INTEGER :: iotdobsmpp ! - outside time domain
104	INTEGER :: iosdtobsmpp ! - outside space domain (temperature)
105	INTEGER :: iosdsobsmpp ! - outside space domain (salinity)
106	INTEGER :: ilantobsmpp ! - within a model land cell (temperature)
107	INTEGER :: ilansobsmpp ! - within a model land cell (salinity)
108	INTEGER :: inlatobsmpp ! - close to land (temperature)
109	INTEGER :: inlasobsmpp ! - close to land (salinity)
110	INTEGER :: igrdobsmpp ! - fail the grid search
111	TYPE(obs_prof_valid) :: llvalid ! Profile selection
112	TYPE(obs_prof_valid), DIMENSION(profdata%nvar) :: &
113	& llvvalid ! T,S selection
114	INTEGER :: jvar ! Variable loop variable
115	INTEGER :: jobs ! Obs. loop variable
116	INTEGER :: jstp ! Time loop variable
117	INTEGER :: inrc ! Time index variable
118
119	IF(lwp) WRITE(numout,*)'obs_pre_pro : Preparing the profile observations...'
120
121	! Initial date initialization (year, month, day, hour, minute)
122	iyea0 = ndate0 / 10000
123	imon0 = ( ndate0 - iyea0 * 10000 ) / 100
124	iday0 = ndate0 - iyea0 * 10000 - imon0 * 100
125	ihou0 = 0
126	imin0 = 0
127
128	icycle = no ! Assimilation cycle
129
130	! Diagnotics counters for various failures.
131
132	iotdobs = 0
133	igrdobs = 0
134	iosdtobs = 0
135	iosdsobs = 0
136	ilantobs = 0
137	ilansobs = 0
138	inlatobs = 0
139	inlasobs = 0
140
141	! -----------------------------------------------------------------------
142	! Find time coordinate for profiles
143	! -----------------------------------------------------------------------
144
145	IF ( PRESENT(kdailyavtypes) ) THEN
146	CALL obs_coo_tim_prof( icycle, &
147	& iyea0, imon0, iday0, ihou0, imin0, &
148	& profdata%nprof, profdata%nyea, profdata%nmon, &
149	& profdata%nday, profdata%nhou, profdata%nmin, &
150	& profdata%ntyp, profdata%nqc, profdata%mstp, &
151	& iotdobs, kdailyavtypes = kdailyavtypes )
152	ELSE
153	CALL obs_coo_tim_prof( icycle, &
154	& iyea0, imon0, iday0, ihou0, imin0, &
155	& profdata%nprof, profdata%nyea, profdata%nmon, &
156	& profdata%nday, profdata%nhou, profdata%nmin, &
157	& profdata%ntyp, profdata%nqc, profdata%mstp, &
158	& iotdobs )
159	ENDIF
160	CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
161
162	! -----------------------------------------------------------------------
163	! Check for profiles failing the grid search
164	! -----------------------------------------------------------------------
165
166	CALL obs_coo_grd( profdata%nprof, profdata%mi, profdata%mj, &
167	& profdata%nqc, igrdobs )
168
169	CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
170
171	! -----------------------------------------------------------------------
172	! Reject all observations for profiles with nqc > 10
173	! -----------------------------------------------------------------------
174
175	CALL obs_pro_rej( profdata )
176
177	! -----------------------------------------------------------------------
178	! Check for land points. This includes points below the model
179	! bathymetry so this is done for every point in the profile
180	! -----------------------------------------------------------------------
181
182	! Temperature
183
184	CALL obs_coo_spc_3d( profdata%nprof, profdata%nvprot(1), &
185	& profdata%npvsta(:,1), profdata%npvend(:,1), &
186	& jpi, jpj, &
187	& jpk, &
188	& profdata%mi, profdata%mj, &
189	& profdata%var(1)%mvk, &
190	& profdata%rlam, profdata%rphi, &
191	& profdata%var(1)%vdep, &
192	& glamt, gphit, &
193	& gdept_0, tmask, &
194	& profdata%nqc, profdata%var(1)%nvqc, &
195	& iosdtobs, ilantobs, &
196	& inlatobs, ld_nea )
197
198	CALL obs_mpp_sum_integer( iosdtobs, iosdtobsmpp )
199	CALL obs_mpp_sum_integer( ilantobs, ilantobsmpp )
200	CALL obs_mpp_sum_integer( inlatobs, inlatobsmpp )
201
202	! Salinity
203
204	CALL obs_coo_spc_3d( profdata%nprof, profdata%nvprot(2), &
205	& profdata%npvsta(:,2), profdata%npvend(:,2), &
206	& jpi, jpj, &
207	& jpk, &
208	& profdata%mi, profdata%mj, &
209	& profdata%var(2)%mvk, &
210	& profdata%rlam, profdata%rphi, &
211	& profdata%var(2)%vdep, &
212	& glamt, gphit, &
213	& gdept_0, tmask, &
214	& profdata%nqc, profdata%var(2)%nvqc, &
215	& iosdsobs, ilansobs, &
216	& inlasobs, ld_nea )
217
218	CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
219	CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
220	CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
221
222	! -----------------------------------------------------------------------
223	! Copy useful data from the profdata data structure to
224	! the prodatqc data structure
225	! -----------------------------------------------------------------------
226
227	! Allocate the selection arrays
228
229	ALLOCATE( llvalid%luse(profdata%nprof) )
230	DO jvar = 1,profdata%nvar
231	ALLOCATE( llvvalid(jvar)%luse(profdata%nvprot(jvar)) )
232	END DO
233
234	! We want all data which has qc flags <= 10
235
236	llvalid%luse(:) = ( profdata%nqc(:) <= 10 )
237	DO jvar = 1,profdata%nvar
238	llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= 10 )
239	END DO
240
241	! The actual copying
242
243	CALL obs_prof_compress( profdata, prodatqc, .TRUE., numout, &
244	& lvalid=llvalid, lvvalid=llvvalid )
245
246	! Dellocate the selection arrays
247	DEALLOCATE( llvalid%luse )
248	DO jvar = 1,profdata%nvar
249	DEALLOCATE( llvvalid(jvar)%luse )
250	END DO
251
252	! -----------------------------------------------------------------------
253	! Print information about what observations are left after qc
254	! -----------------------------------------------------------------------
255
256	! Update the total observation counter array
257
258	IF(lwp) THEN
259	WRITE(numout,*)
260	WRITE(numout,*) 'obs_pre_pro :'
261	WRITE(numout,*) '~~~~~~~~~~~'
262	WRITE(numout,*)
263	WRITE(numout,*) ' Profiles outside time domain = ', &
264	& iotdobsmpp
265	WRITE(numout,*) ' Remaining profiles that failed grid search = ', &
266	& igrdobsmpp
267	WRITE(numout,*) ' Remaining T data outside space domain = ', &
268	& iosdtobsmpp
269	WRITE(numout,*) ' Remaining T data at land points = ', &
270	& ilantobsmpp
271	IF (ld_nea) THEN
272	WRITE(numout,*) ' Remaining T data near land points (removed) = ',&
273	& inlatobsmpp
274	ELSE
275	WRITE(numout,*) ' Remaining T data near land points (kept) = ',&
276	& inlatobsmpp
277	ENDIF
278	WRITE(numout,*) ' T data accepted = ', &
279	& prodatqc%nvprotmpp(1)
280	WRITE(numout,*) ' Remaining S data outside space domain = ', &
281	& iosdsobsmpp
282	WRITE(numout,*) ' Remaining S data at land points = ', &
283	& ilansobsmpp
284	IF (ld_nea) THEN
285	WRITE(numout,*) ' Remaining S data near land points (removed) = ',&
286	& inlasobsmpp
287	ELSE
288	WRITE(numout,*) ' Remaining S data near land points (kept) = ',&
289	& inlasobsmpp
290	ENDIF
291	WRITE(numout,*) ' S data accepted = ', &
292	& prodatqc%nvprotmpp(2)
293
294	WRITE(numout,*)
295	WRITE(numout,*) ' Number of observations per time step :'
296	WRITE(numout,*)
297	WRITE(numout,997)
298	WRITE(numout,998)
299	ENDIF
300
301	DO jobs = 1, prodatqc%nprof
302	inrc = prodatqc%mstp(jobs) + 2 - nit000
303	prodatqc%npstp(inrc) = prodatqc%npstp(inrc) + 1
304	DO jvar = 1, prodatqc%nvar
305	IF ( prodatqc%npvend(jobs,jvar) > 0 ) THEN
306	prodatqc%nvstp(inrc,jvar) = prodatqc%nvstp(inrc,jvar) + &
307	& ( prodatqc%npvend(jobs,jvar) - &
308	& prodatqc%npvsta(jobs,jvar) + 1 )
309	ENDIF
310	END DO
311	END DO
312
313
314	CALL obs_mpp_sum_integers( prodatqc%npstp, prodatqc%npstpmpp, &
315	& nitend - nit000 + 2 )
316	DO jvar = 1, prodatqc%nvar
317	CALL obs_mpp_sum_integers( prodatqc%nvstp(:,jvar), &
318	& prodatqc%nvstpmpp(:,jvar), &
319	& nitend - nit000 + 2 )
320	END DO
321
322	IF ( lwp ) THEN
323	DO jstp = nit000 - 1, nitend
324	inrc = jstp - nit000 + 2
325	WRITE(numout,999) jstp, prodatqc%npstpmpp(inrc), &
326	& prodatqc%nvstpmpp(inrc,1), &
327	& prodatqc%nvstpmpp(inrc,2)
328	END DO
329	ENDIF
330
331	997 FORMAT(10X,'Time step',5X,'Profiles',5X,'Temperature',5X,'Salinity')
332	998 FORMAT(10X,'---------',5X,'--------',5X,'-----------',5X,'--------')
333	999 FORMAT(10X,I9,5X,I8,5X,I11,5X,I8)
334
335	END SUBROUTINE obs_pre_pro
336
337	SUBROUTINE obs_pre_sla( sladata, sladatqc, ld_sla, ld_nea )
338	!!----------------------------------------------------------------------
339	!! * ROUTINE obs_pre_sla *
340	!!
341	!! ** Purpose : First level check and screening of SLA observations
342	!!
343	!! ** Method : First level check and screening of SLA observations
344	!!
345	!! ** Action :
346	!!
347	!! References :
348	!!
349	!! History :
350	!! ! 2007-03 (A. Weaver, K. Mogensen) Original
351	!! ! 2007-06 (K. Mogensen et al) Reject obs. near land.
352	!!----------------------------------------------------------------------
353	!! * Modules used
354	USE domstp ! Domain: set the time-step
355	USE par_oce ! Ocean parameters
356	USE dom_oce, ONLY : & ! Geographical information
357	& glamt, &
358	& gphit, &
359	& tmask, &
360	& nproc
361	!! * Arguments
362	TYPE(obs_surf), INTENT(INOUT) :: sladata ! Full set of SLA data
363	TYPE(obs_surf), INTENT(INOUT) :: sladatqc ! Subset of SLA data not failing screening
364	LOGICAL, INTENT(IN) :: ld_sla ! Switch for SLA data
365	LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land
366	!! * Local declarations
367	INTEGER :: iyea0 ! Initial date
368	INTEGER :: imon0 ! - (year, month, day, hour, minute)
369	INTEGER :: iday0
370	INTEGER :: ihou0
371	INTEGER :: imin0
372	INTEGER :: icycle ! Current assimilation cycle
373	! Counters for observations that
374	INTEGER :: iotdobs ! - outside time domain
375	INTEGER :: iosdsobs ! - outside space domain
376	INTEGER :: ilansobs ! - within a model land cell
377	INTEGER :: inlasobs ! - close to land
378	INTEGER :: igrdobs ! - fail the grid search
379	! Global counters for observations that
380	INTEGER :: iotdobsmpp ! - outside time domain
381	INTEGER :: iosdsobsmpp ! - outside space domain
382	INTEGER :: ilansobsmpp ! - within a model land cell
383	INTEGER :: inlasobsmpp ! - close to land
384	INTEGER :: igrdobsmpp ! - fail the grid search
385	LOGICAL, DIMENSION(:), ALLOCATABLE :: &
386	& llvalid ! SLA data selection
387	INTEGER :: jobs ! Obs. loop variable
388	INTEGER :: jstp ! Time loop variable
389	INTEGER :: inrc ! Time index variable
390
391	IF(lwp) WRITE(numout,*)'obs_pre_sla : Preparing the SLA observations...'
392
393	! Initial date initialization (year, month, day, hour, minute)
394	iyea0 = ndate0 / 10000
395	imon0 = ( ndate0 - iyea0 * 10000 ) / 100
396	iday0 = ndate0 - iyea0 * 10000 - imon0 * 100
397	ihou0 = 0
398	imin0 = 0
399
400	icycle = no ! Assimilation cycle
401
402	! Diagnotics counters for various failures.
403
404	iotdobs = 0
405	igrdobs = 0
406	iosdsobs = 0
407	ilansobs = 0
408	inlasobs = 0
409
410	! -----------------------------------------------------------------------
411	! Find time coordinate for SLA data
412	! -----------------------------------------------------------------------
413
414	CALL obs_coo_tim( icycle, &
415	& iyea0, imon0, iday0, ihou0, imin0, &
416	& sladata%nsurf, sladata%nyea, sladata%nmon, &
417	& sladata%nday, sladata%nhou, sladata%nmin, &
418	& sladata%nqc, sladata%mstp, iotdobs )
419
420	CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
421
422	! -----------------------------------------------------------------------
423	! Check for SLA data failing the grid search
424	! -----------------------------------------------------------------------
425
426	CALL obs_coo_grd( sladata%nsurf, sladata%mi, sladata%mj, &
427	& sladata%nqc, igrdobs )
428
429	CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
430
431	! -----------------------------------------------------------------------
432	! Check for land points.
433	! -----------------------------------------------------------------------
434
435	CALL obs_coo_spc_2d( sladata%nsurf, &
436	& jpi, jpj, &
437	& sladata%mi, sladata%mj, &
438	& sladata%rlam, sladata%rphi, &
439	& glamt, gphit, &
440	& tmask(:,:,1), sladata%nqc, &
441	& iosdsobs, ilansobs, &
442	& inlasobs, ld_nea )
443
444	CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
445	CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
446	CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
447
448	! -----------------------------------------------------------------------
449	! Copy useful data from the sladata data structure to
450	! the sladatqc data structure
451	! -----------------------------------------------------------------------
452
453	! Allocate the selection arrays
454
455	ALLOCATE( llvalid(sladata%nsurf) )
456
457	! We want all data which has qc flags <= 10
458
459	llvalid(:) = ( sladata%nqc(:) <= 10 )
460
461	! The actual copying
462
463	CALL obs_surf_compress( sladata, sladatqc, .TRUE., numout, &
464	& lvalid=llvalid )
465
466	! Dellocate the selection arrays
467	DEALLOCATE( llvalid )
468
469	! -----------------------------------------------------------------------
470	! Print information about what observations are left after qc
471	! -----------------------------------------------------------------------
472
473	! Update the total observation counter array
474
475	IF(lwp) THEN
476	WRITE(numout,*)
477	WRITE(numout,*) 'obs_pre_sla :'
478	WRITE(numout,*) '~~~~~~~~~~~'
479	WRITE(numout,*)
480	WRITE(numout,*) ' SLA data outside time domain = ', &
481	& iotdobsmpp
482	WRITE(numout,*) ' Remaining SLA data that failed grid search = ', &
483	& igrdobsmpp
484	WRITE(numout,*) ' Remaining SLA data outside space domain = ', &
485	& iosdsobsmpp
486	WRITE(numout,*) ' Remaining SLA data at land points = ', &
487	& ilansobsmpp
488	IF (ld_nea) THEN
489	WRITE(numout,*) ' Remaining SLA data near land points (removed) = ', &
490	& inlasobsmpp
491	ELSE
492	WRITE(numout,*) ' Remaining SLA data near land points (kept) = ', &
493	& inlasobsmpp
494	ENDIF
495	WRITE(numout,*) ' SLA data accepted = ', &
496	& sladatqc%nsurfmpp
497
498	WRITE(numout,*)
499	WRITE(numout,*) ' Number of observations per time step :'
500	WRITE(numout,*)
501	WRITE(numout,1997)
502	WRITE(numout,1998)
503	ENDIF
504
505	DO jobs = 1, sladatqc%nsurf
506	inrc = sladatqc%mstp(jobs) + 2 - nit000
507	sladatqc%nsstp(inrc) = sladatqc%nsstp(inrc) + 1
508	END DO
509
510	CALL obs_mpp_sum_integers( sladatqc%nsstp, sladatqc%nsstpmpp, &
511	& nitend - nit000 + 2 )
512
513	IF ( lwp ) THEN
514	DO jstp = nit000 - 1, nitend
515	inrc = jstp - nit000 + 2
516	WRITE(numout,1999) jstp, sladatqc%nsstpmpp(inrc)
517	END DO
518	ENDIF
519
520	1997 FORMAT(10X,'Time step',5X,'Sea level anomaly')
521	1998 FORMAT(10X,'---------',5X,'-----------------')
522	1999 FORMAT(10X,I9,5X,I17)
523
524	END SUBROUTINE obs_pre_sla
525
526	SUBROUTINE obs_pre_sst( sstdata, sstdatqc, ld_sst, ld_nea )
527	!!----------------------------------------------------------------------
528	!! * ROUTINE obs_pre_sst *
529	!!
530	!! ** Purpose : First level check and screening of SST observations
531	!!
532	!! ** Method : First level check and screening of SST observations
533	!!
534	!! ** Action :
535	!!
536	!! References :
537	!!
538	!! History :
539	!! ! 2007-03 (S. Ricci) SST data preparation
540	!!----------------------------------------------------------------------
541	!! * Modules used
542	USE domstp ! Domain: set the time-step
543	USE par_oce ! Ocean parameters
544	USE dom_oce, ONLY : & ! Geographical information
545	& glamt, &
546	& gphit, &
547	& tmask, &
548	& nproc
549	!! * Arguments
550	TYPE(obs_surf), INTENT(INOUT) :: sstdata ! Full set of SST data
551	TYPE(obs_surf), INTENT(INOUT) :: sstdatqc ! Subset of SST data not failing screening
552	LOGICAL :: ld_sst ! Switch for SST data
553	LOGICAL :: ld_nea ! Switch for rejecting observation near land
554	!! * Local declarations
555	INTEGER :: iyea0 ! Initial date
556	INTEGER :: imon0 ! - (year, month, day, hour, minute)
557	INTEGER :: iday0
558	INTEGER :: ihou0
559	INTEGER :: imin0
560	INTEGER :: icycle ! Current assimilation cycle
561	! Counters for observations that
562	INTEGER :: iotdobs ! - outside time domain
563	INTEGER :: iosdsobs ! - outside space domain
564	INTEGER :: ilansobs ! - within a model land cell
565	INTEGER :: inlasobs ! - close to land
566	INTEGER :: igrdobs ! - fail the grid search
567	! Global counters for observations that
568	INTEGER :: iotdobsmpp ! - outside time domain
569	INTEGER :: iosdsobsmpp ! - outside space domain
570	INTEGER :: ilansobsmpp ! - within a model land cell
571	INTEGER :: inlasobsmpp ! - close to land
572	INTEGER :: igrdobsmpp ! - fail the grid search
573	LOGICAL, DIMENSION(:), ALLOCATABLE :: &
574	& llvalid ! SST data selection
575	INTEGER :: jobs ! Obs. loop variable
576	INTEGER :: jstp ! Time loop variable
577	INTEGER :: inrc ! Time index variable
578
579	IF(lwp) WRITE(numout,*)'obs_pre_sst : Preparing the SST observations...'
580
581	! Initial date initialization (year, month, day, hour, minute)
582	iyea0 = ndate0 / 10000
583	imon0 = ( ndate0 - iyea0 * 10000 ) / 100
584	iday0 = ndate0 - iyea0 * 10000 - imon0 * 100
585	ihou0 = 0
586	imin0 = 0
587
588	icycle = no ! Assimilation cycle
589
590	! Diagnotics counters for various failures.
591
592	iotdobs = 0
593	igrdobs = 0
594	iosdsobs = 0
595	ilansobs = 0
596	inlasobs = 0
597
598	! -----------------------------------------------------------------------
599	! Find time coordinate for SST data
600	! -----------------------------------------------------------------------
601
602	CALL obs_coo_tim( icycle, &
603	& iyea0, imon0, iday0, ihou0, imin0, &
604	& sstdata%nsurf, sstdata%nyea, sstdata%nmon, &
605	& sstdata%nday, sstdata%nhou, sstdata%nmin, &
606	& sstdata%nqc, sstdata%mstp, iotdobs )
607	CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
608	! -----------------------------------------------------------------------
609	! Check for SST data failing the grid search
610	! -----------------------------------------------------------------------
611
612	CALL obs_coo_grd( sstdata%nsurf, sstdata%mi, sstdata%mj, &
613	& sstdata%nqc, igrdobs )
614	CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
615
616	! -----------------------------------------------------------------------
617	! Check for land points.
618	! -----------------------------------------------------------------------
619
620	CALL obs_coo_spc_2d( sstdata%nsurf, &
621	& jpi, jpj, &
622	& sstdata%mi, sstdata%mj, &
623	& sstdata%rlam, sstdata%rphi, &
624	& glamt, gphit, &
625	& tmask(:,:,1), sstdata%nqc, &
626	& iosdsobs, ilansobs, &
627	& inlasobs, ld_nea )
628
629	CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
630	CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
631	CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
632
633	! -----------------------------------------------------------------------
634	! Copy useful data from the sstdata data structure to
635	! the sstdatqc data structure
636	! -----------------------------------------------------------------------
637
638	! Allocate the selection arrays
639
640	ALLOCATE( llvalid(sstdata%nsurf) )
641
642	! We want all data which has qc flags <= 0
643
644	llvalid(:) = ( sstdata%nqc(:) <= 10 )
645
646	! The actual copying
647
648	CALL obs_surf_compress( sstdata, sstdatqc, .TRUE., numout, &
649	& lvalid=llvalid )
650
651	! Dellocate the selection arrays
652	DEALLOCATE( llvalid )
653
654	! -----------------------------------------------------------------------
655	! Print information about what observations are left after qc
656	! -----------------------------------------------------------------------
657
658	! Update the total observation counter array
659
660	IF(lwp) THEN
661	WRITE(numout,*)
662	WRITE(numout,*) 'obs_pre_sst :'
663	WRITE(numout,*) '~~~~~~~~~~~'
664	WRITE(numout,*)
665	WRITE(numout,*) ' SST data outside time domain = ', &
666	& iotdobsmpp
667	WRITE(numout,*) ' Remaining SST data that failed grid search = ', &
668	& igrdobsmpp
669	WRITE(numout,*) ' Remaining SST data outside space domain = ', &
670	& iosdsobsmpp
671	WRITE(numout,*) ' Remaining SST data at land points = ', &
672	& ilansobsmpp
673	IF (ld_nea) THEN
674	WRITE(numout,*) ' Remaining SST data near land points (removed) = ', &
675	& inlasobsmpp
676	ELSE
677	WRITE(numout,*) ' Remaining SST data near land points (kept) = ', &
678	& inlasobsmpp
679	ENDIF
680	WRITE(numout,*) ' SST data accepted = ', &
681	& sstdatqc%nsurfmpp
682
683	WRITE(numout,*)
684	WRITE(numout,*) ' Number of observations per time step :'
685	WRITE(numout,*)
686	WRITE(numout,1997)
687	WRITE(numout,1998)
688	ENDIF
689
690	DO jobs = 1, sstdatqc%nsurf
691	inrc = sstdatqc%mstp(jobs) + 2 - nit000
692	sstdatqc%nsstp(inrc) = sstdatqc%nsstp(inrc) + 1
693	END DO
694
695	CALL obs_mpp_sum_integers( sstdatqc%nsstp, sstdatqc%nsstpmpp, &
696	& nitend - nit000 + 2 )
697
698	IF ( lwp ) THEN
699	DO jstp = nit000 - 1, nitend
700	inrc = jstp - nit000 + 2
701	WRITE(numout,1999) jstp, sstdatqc%nsstpmpp(inrc)
702	END DO
703	ENDIF
704
705	1997 FORMAT(10X,'Time step',5X,'Sea surface temperature')
706	1998 FORMAT(10X,'---------',5X,'-----------------')
707	1999 FORMAT(10X,I9,5X,I17)
708
709	END SUBROUTINE obs_pre_sst
710
711	SUBROUTINE obs_pre_seaice( seaicedata, seaicedatqc, ld_seaice, ld_nea )
712	!!----------------------------------------------------------------------
713	!! * ROUTINE obs_pre_seaice *
714	!!
715	!! ** Purpose : First level check and screening of Sea Ice observations
716	!!
717	!! ** Method : First level check and screening of Sea Ice observations
718	!!
719	!! ** Action :
720	!!
721	!! References :
722	!!
723	!! History :
724	!! ! 2007-11 (D. Lea) based on obs_pre_sst
725	!!----------------------------------------------------------------------
726	!! * Modules used
727	USE domstp ! Domain: set the time-step
728	USE par_oce ! Ocean parameters
729	USE dom_oce, ONLY : & ! Geographical information
730	& glamt, &
731	& gphit, &
732	& tmask, &
733	& nproc
734	!! * Arguments
735	TYPE(obs_surf), INTENT(INOUT) :: seaicedata ! Full set of Sea Ice data
736	TYPE(obs_surf), INTENT(INOUT) :: seaicedatqc ! Subset of sea ice data not failing screening
737	LOGICAL :: ld_seaice ! Switch for sea ice data
738	LOGICAL :: ld_nea ! Switch for rejecting observation near land
739	!! * Local declarations
740	INTEGER :: iyea0 ! Initial date
741	INTEGER :: imon0 ! - (year, month, day, hour, minute)
742	INTEGER :: iday0
743	INTEGER :: ihou0
744	INTEGER :: imin0
745	INTEGER :: icycle ! Current assimilation cycle
746	! Counters for observations that
747	INTEGER :: iotdobs ! - outside time domain
748	INTEGER :: iosdsobs ! - outside space domain
749	INTEGER :: ilansobs ! - within a model land cell
750	INTEGER :: inlasobs ! - close to land
751	INTEGER :: igrdobs ! - fail the grid search
752	! Global counters for observations that
753	INTEGER :: iotdobsmpp ! - outside time domain
754	INTEGER :: iosdsobsmpp ! - outside space domain
755	INTEGER :: ilansobsmpp ! - within a model land cell
756	INTEGER :: inlasobsmpp ! - close to land
757	INTEGER :: igrdobsmpp ! - fail the grid search
758	LOGICAL, DIMENSION(:), ALLOCATABLE :: &
759	& llvalid ! data selection
760	INTEGER :: jobs ! Obs. loop variable
761	INTEGER :: jstp ! Time loop variable
762	INTEGER :: inrc ! Time index variable
763
764	IF (lwp) WRITE(numout,*)'obs_pre_seaice : Preparing the sea ice observations...'
765
766	! Initial date initialization (year, month, day, hour, minute)
767	iyea0 = ndate0 / 10000
768	imon0 = ( ndate0 - iyea0 * 10000 ) / 100
769	iday0 = ndate0 - iyea0 * 10000 - imon0 * 100
770	ihou0 = 0
771	imin0 = 0
772
773	icycle = no ! Assimilation cycle
774
775	! Diagnotics counters for various failures.
776
777	iotdobs = 0
778	igrdobs = 0
779	iosdsobs = 0
780	ilansobs = 0
781	inlasobs = 0
782
783	! -----------------------------------------------------------------------
784	! Find time coordinate for sea ice data
785	! -----------------------------------------------------------------------
786
787	CALL obs_coo_tim( icycle, &
788	& iyea0, imon0, iday0, ihou0, imin0, &
789	& seaicedata%nsurf, seaicedata%nyea, seaicedata%nmon, &
790	& seaicedata%nday, seaicedata%nhou, seaicedata%nmin, &
791	& seaicedata%nqc, seaicedata%mstp, iotdobs )
792	CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
793	! -----------------------------------------------------------------------
794	! Check for sea ice data failing the grid search
795	! -----------------------------------------------------------------------
796
797	CALL obs_coo_grd( seaicedata%nsurf, seaicedata%mi, seaicedata%mj, &
798	& seaicedata%nqc, igrdobs )
799	CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
800
801	! -----------------------------------------------------------------------
802	! Check for land points.
803	! -----------------------------------------------------------------------
804
805	CALL obs_coo_spc_2d( seaicedata%nsurf, &
806	& jpi, jpj, &
807	& seaicedata%mi, seaicedata%mj, &
808	& seaicedata%rlam, seaicedata%rphi, &
809	& glamt, gphit, &
810	& tmask(:,:,1), seaicedata%nqc, &
811	& iosdsobs, ilansobs, &
812	& inlasobs, ld_nea )
813
814	CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
815	CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
816	CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
817
818	! -----------------------------------------------------------------------
819	! Copy useful data from the seaicedata data structure to
820	! the seaicedatqc data structure
821	! -----------------------------------------------------------------------
822
823	! Allocate the selection arrays
824
825	ALLOCATE( llvalid(seaicedata%nsurf) )
826
827	! We want all data which has qc flags <= 0
828
829	llvalid(:) = ( seaicedata%nqc(:) <= 10 )
830
831	! The actual copying
832
833	CALL obs_surf_compress( seaicedata, seaicedatqc, .TRUE., numout, &
834	& lvalid=llvalid )
835
836	! Dellocate the selection arrays
837	DEALLOCATE( llvalid )
838
839	! -----------------------------------------------------------------------
840	! Print information about what observations are left after qc
841	! -----------------------------------------------------------------------
842
843	! Update the total observation counter array
844
845	IF(lwp) THEN
846	WRITE(numout,*)
847	WRITE(numout,*) 'obs_pre_seaice :'
848	WRITE(numout,*) '~~~~~~~~~~~'
849	WRITE(numout,*)
850	WRITE(numout,*) ' Sea ice data outside time domain = ', &
851	& iotdobsmpp
852	WRITE(numout,*) ' Remaining sea ice data that failed grid search = ', &
853	& igrdobsmpp
854	WRITE(numout,*) ' Remaining sea ice data outside space domain = ', &
855	& iosdsobsmpp
856	WRITE(numout,*) ' Remaining sea ice data at land points = ', &
857	& ilansobsmpp
858	IF (ld_nea) THEN
859	WRITE(numout,*) ' Remaining sea ice data near land points (removed) = ', &
860	& inlasobsmpp
861	ELSE
862	WRITE(numout,*) ' Remaining sea ice data near land points (kept) = ', &
863	& inlasobsmpp
864	ENDIF
865	WRITE(numout,*) ' Sea ice data accepted = ', &
866	& seaicedatqc%nsurfmpp
867
868	WRITE(numout,*)
869	WRITE(numout,*) ' Number of observations per time step :'
870	WRITE(numout,*)
871	WRITE(numout,1997)
872	WRITE(numout,1998)
873	ENDIF
874
875	DO jobs = 1, seaicedatqc%nsurf
876	inrc = seaicedatqc%mstp(jobs) + 2 - nit000
877	seaicedatqc%nsstp(inrc) = seaicedatqc%nsstp(inrc) + 1
878	END DO
879
880	CALL obs_mpp_sum_integers( seaicedatqc%nsstp, seaicedatqc%nsstpmpp, &
881	& nitend - nit000 + 2 )
882
883	IF ( lwp ) THEN
884	DO jstp = nit000 - 1, nitend
885	inrc = jstp - nit000 + 2
886	WRITE(numout,1999) jstp, seaicedatqc%nsstpmpp(inrc)
887	END DO
888	ENDIF
889
890	1997 FORMAT(10X,'Time step',5X,'Sea ice data ')
891	1998 FORMAT(10X,'---------',5X,'-----------------')
892	1999 FORMAT(10X,I9,5X,I17)
893
894	END SUBROUTINE obs_pre_seaice
895
896	SUBROUTINE obs_pre_vel( profdata, prodatqc, ld_vel3d, ld_nea, ld_dailyav )
897	!!----------------------------------------------------------------------
898	!! * ROUTINE obs_pre_taovel *
899	!!
900	!! ** Purpose : First level check and screening of U and V profiles
901	!!
902	!! ** Method : First level check and screening of U and V profiles
903	!!
904	!! History :
905	!! ! 2007-06 (K. Mogensen) original : T and S profile data
906	!! ! 2008-09 (M. Valdivieso) : TAO velocity data
907	!! ! 2009-01 (K. Mogensen) : New feedback strictuer
908	!!
909	!!----------------------------------------------------------------------
910	!! * Modules used
911	USE domstp ! Domain: set the time-step
912	USE par_oce ! Ocean parameters
913	USE dom_oce, ONLY : & ! Geographical information
914	& glamt, glamu, glamv, &
915	& gphit, gphiu, gphiv, &
916	& gdept_0, &
917	& tmask, umask, vmask, &
918	& nproc
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_0, 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_0, 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_0
1710
1711	!! * Arguments
1712	INTEGER, INTENT(IN) :: kprofno ! Number of profiles
1713	INTEGER, INTENT(IN) :: kobsno ! Total number of observations
1714	INTEGER, INTENT(IN) :: kpi ! Number of grid points in (i,j,k)
1715	INTEGER, INTENT(IN) :: kpj
1716	INTEGER, INTENT(IN) :: kpk
1717	INTEGER, DIMENSION(kprofno), INTENT(IN) :: &
1718	& kpstart, & ! Start of individual profiles
1719	& kpend ! End of individual profiles
1720	INTEGER, DIMENSION(kprofno), INTENT(IN) :: &
1721	& kobsi, & ! Observation (i,j) coordinates
1722	& kobsj
1723	INTEGER, DIMENSION(kobsno), INTENT(IN) :: &
1724	& kobsk ! Observation k coordinate
1725	REAL(KIND=wp), DIMENSION(kprofno), INTENT(IN) :: &
1726	& pobslam, & ! Observation (lon,lat) coordinates
1727	& pobsphi
1728	REAL(KIND=wp), DIMENSION(kobsno), INTENT(INOUT) :: &
1729	& pobsdep ! Observation depths
1730	REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) :: &
1731	& plam, pphi ! Model (lon,lat) coordinates
1732	REAL(KIND=wp), DIMENSION(kpk), INTENT(IN) :: &
1733	& pdep ! Model depth coordinates
1734	REAL(KIND=wp), DIMENSION(kpi,kpj,kpk), INTENT(IN) :: &
1735	& pmask ! Land mask array
1736	INTEGER, DIMENSION(kprofno), INTENT(INOUT) :: &
1737	& kpobsqc ! Profile quality control
1738	INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
1739	& kobsqc ! Observation quality control
1740	INTEGER, INTENT(INOUT) :: kosdobs ! Observations outside space domain
1741	INTEGER, INTENT(INOUT) :: klanobs ! Observations within a model land cell
1742	INTEGER, INTENT(INOUT) :: knlaobs ! Observations near land
1743	LOGICAL, INTENT(IN) :: ld_nea ! Flag observations near land
1744	!! * Local declarations
1745	REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: &
1746	& zgmsk ! Grid mask
1747	REAL(KIND=wp), DIMENSION(2,2,kprofno) :: &
1748	& zglam, & ! Model longitude at grid points
1749	& zgphi ! Model latitude at grid points
1750	INTEGER, DIMENSION(2,2,kprofno) :: &
1751	& igrdi, & ! Grid i,j
1752	& igrdj
1753	LOGICAL :: lgridobs ! Is observation on a model grid point.
1754	INTEGER :: iig, ijg ! i,j of observation on model grid point.
1755	INTEGER :: jobs, jobsp, jk, ji, jj
1756
1757	! Get grid point indices
1758
1759	DO jobs = 1, kprofno
1760
1761	! For invalid points use 2,2
1762
1763	IF ( kpobsqc(jobs) >= 10 ) THEN
1764
1765	igrdi(1,1,jobs) = 1
1766	igrdj(1,1,jobs) = 1
1767	igrdi(1,2,jobs) = 1
1768	igrdj(1,2,jobs) = 2
1769	igrdi(2,1,jobs) = 2
1770	igrdj(2,1,jobs) = 1
1771	igrdi(2,2,jobs) = 2
1772	igrdj(2,2,jobs) = 2
1773
1774	ELSE
1775
1776	igrdi(1,1,jobs) = kobsi(jobs)-1
1777	igrdj(1,1,jobs) = kobsj(jobs)-1
1778	igrdi(1,2,jobs) = kobsi(jobs)-1
1779	igrdj(1,2,jobs) = kobsj(jobs)
1780	igrdi(2,1,jobs) = kobsi(jobs)
1781	igrdj(2,1,jobs) = kobsj(jobs)-1
1782	igrdi(2,2,jobs) = kobsi(jobs)
1783	igrdj(2,2,jobs) = kobsj(jobs)
1784
1785	ENDIF
1786
1787	END DO
1788
1789	CALL obs_int_comm_3d( 2, 2, kprofno, kpk, igrdi, igrdj, pmask, zgmsk )
1790	CALL obs_int_comm_2d( 2, 2, kprofno, igrdi, igrdj, plam, zglam )
1791	CALL obs_int_comm_2d( 2, 2, kprofno, igrdi, igrdj, pphi, zgphi )
1792
1793	DO jobs = 1, kprofno
1794
1795	! Skip bad profiles
1796	IF ( kpobsqc(jobs) >= 10 ) CYCLE
1797
1798	! Check if this observation is on a grid point
1799
1800	lgridobs = .FALSE.
1801	iig = -1
1802	ijg = -1
1803	DO jj = 1, 2
1804	DO ji = 1, 2
1805	IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) &
1806	& .AND. &
1807	& ( ABS( zglam(ji,jj,jobs) - pobslam(jobs) ) < 1.0e-6_wp ) &
1808	& ) THEN
1809	lgridobs = .TRUE.
1810	iig = ji
1811	ijg = jj
1812	ENDIF
1813	END DO
1814	END DO
1815
1816	! Reject observations
1817
1818	DO jobsp = kpstart(jobs), kpend(jobs)
1819
1820	! Flag if the observation falls outside the model spatial domain
1821	IF ( ( pobslam(jobs) < -180. ) &
1822	& .OR. ( pobslam(jobs) > 180. ) &
1823	& .OR. ( pobsphi(jobs) < -90. ) &
1824	& .OR. ( pobsphi(jobs) > 90. ) &
1825	& .OR. ( pobsdep(jobsp) < 0.0 ) &
1826	& .OR. ( pobsdep(jobsp) > gdepw_0(kpk) ) ) THEN
1827	kobsqc(jobsp) = kobsqc(jobsp) + 11
1828	kosdobs = kosdobs + 1
1829	CYCLE
1830	ENDIF
1831
1832	! Flag if the observation falls with a model land cell
1833	IF ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) &
1834	& == 0.0_wp ) THEN
1835	kobsqc(jobsp) = kobsqc(jobsp) + 12
1836	klanobs = klanobs + 1
1837	CYCLE
1838	ENDIF
1839
1840	! For observations on the grid reject them if their are at
1841	! a masked point
1842
1843	IF (lgridobs) THEN
1844	IF (zgmsk(iig,ijg,kobsk(jobsp)-1,jobs) == 0.0_wp ) THEN
1845	kobsqc(jobsp) = kobsqc(jobsp) + 12
1846	klanobs = klanobs + 1
1847	CYCLE
1848	ENDIF
1849	ENDIF
1850
1851	! Flag if the observation falls is close to land
1852	IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == &
1853	& 0.0_wp) THEN
1854	IF (ld_nea) kobsqc(jobsp) = kobsqc(jobsp) + 14
1855	knlaobs = knlaobs + 1
1856	ENDIF
1857
1858	! Set observation depth equal to that of the first model depth
1859	IF ( pobsdep(jobsp) <= pdep(1) ) THEN
1860	pobsdep(jobsp) = pdep(1)
1861	ENDIF
1862
1863	END DO
1864	END DO
1865
1866	END SUBROUTINE obs_coo_spc_3d
1867
1868	SUBROUTINE obs_pro_rej( profdata )
1869	!!----------------------------------------------------------------------
1870	!! * ROUTINE obs_pro_rej *
1871	!!
1872	!! ** Purpose : Reject all data within a rejected profile
1873	!!
1874	!! ** Method :
1875	!!
1876	!! ** Action :
1877	!!
1878	!! References :
1879	!!
1880	!! History :
1881	!! ! 2007-10 (K. Mogensen) Original code
1882	!!----------------------------------------------------------------------
1883	!! * Modules used
1884	!! * Arguments
1885	TYPE(obs_prof), INTENT(INOUT) :: profdata ! Profile data
1886	!! * Local declarations
1887	INTEGER :: jprof
1888	INTEGER :: jvar
1889	INTEGER :: jobs
1890
1891	! Loop over profiles
1892
1893	DO jprof = 1, profdata%nprof
1894
1895	IF ( profdata%nqc(jprof) > 10 ) THEN
1896
1897	DO jvar = 1, profdata%nvar
1898
1899	DO jobs = profdata%npvsta(jprof,jvar), &
1900	& profdata%npvend(jprof,jvar)
1901
1902	profdata%var(jvar)%nvqc(jobs) = &
1903	& profdata%var(jvar)%nvqc(jobs) + 26
1904
1905	END DO
1906
1907	END DO
1908
1909	ENDIF
1910
1911	END DO
1912
1913	END SUBROUTINE obs_pro_rej
1914
1915	SUBROUTINE obs_uv_rej( profdata, knumu, knumv )
1916	!!----------------------------------------------------------------------
1917	!! * ROUTINE obs_uv_rej *
1918	!!
1919	!! ** Purpose : Reject u if v is rejected and vice versa
1920	!!
1921	!! ** Method :
1922	!!
1923	!! ** Action :
1924	!!
1925	!! References :
1926	!!
1927	!! History :
1928	!! ! 2009-2 (K. Mogensen) Original code
1929	!!----------------------------------------------------------------------
1930	!! * Modules used
1931	!! * Arguments
1932	TYPE(obs_prof), INTENT(INOUT) :: profdata ! Profile data
1933	INTEGER, INTENT(INOUT) :: knumu ! Number of u rejected
1934	INTEGER, INTENT(INOUT) :: knumv ! Number of v rejected
1935	!! * Local declarations
1936	INTEGER :: jprof
1937	INTEGER :: jvar
1938	INTEGER :: jobs
1939
1940	! Loop over profiles
1941
1942	DO jprof = 1, profdata%nprof
1943
1944	IF ( ( profdata%npvsta(jprof,1) /= profdata%npvsta(jprof,2) ) .OR. &
1945	& ( profdata%npvend(jprof,1) /= profdata%npvend(jprof,2) ) ) THEN
1946
1947	CALL ctl_stop('U,V profiles inconsistent in obs_uv_rej')
1948	RETURN
1949
1950	ENDIF
1951
1952	DO jobs = profdata%npvsta(jprof,1), profdata%npvend(jprof,1)
1953
1954	IF ( ( profdata%var(1)%nvqc(jobs) > 10 ) .AND. &
1955	& ( profdata%var(2)%nvqc(jobs) <= 10) ) THEN
1956	profdata%var(2)%nvqc(jobs) = profdata%var(2)%nvqc(jobs) + 42
1957	knumv = knumv + 1
1958	ENDIF
1959	IF ( ( profdata%var(2)%nvqc(jobs) > 10 ) .AND. &
1960	& ( profdata%var(1)%nvqc(jobs) <= 10) ) THEN
1961	profdata%var(1)%nvqc(jobs) = profdata%var(1)%nvqc(jobs) + 42
1962	knumu = knumu + 1
1963	ENDIF
1964
1965	END DO
1966
1967	END DO
1968
1969	END SUBROUTINE obs_uv_rej
1970
1971	END MODULE obs_prep

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

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

Context Navigation

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 2287

Download in other formats: