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

Context Navigation

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

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

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

Download in other formats: