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

Context Navigation

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 4460

Last change on this file since 4460 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.

Property svn:keywords set to Id

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

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

Download in other formats: