Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

obs_prep.F90 @ 15799

Last change on this file since 15799 was 15799, checked in by dford, 2 years ago
More generic interface and structure for OBS code. See Met Office utils tickets 471 and 530.
File size: 62.1 KB

Line
1	MODULE obs_prep
2	!!=====================================================================
3	!! * MODULE obs_prep *
4	!! Observation diagnostics: Prepare observation arrays: screening,
5	!! sorting, coordinate search
6	!!=====================================================================
7
8	!!---------------------------------------------------------------------
9	!! obs_pre_prof : First level check and screening of profile observations
10	!! obs_pre_surf : First level check and screening of surface observations
11	!! obs_scr : Basic screening of the observations
12	!! obs_coo_tim : Compute number of time steps to the observation time
13	!! obs_sor : Sort the observation arrays
14	!!---------------------------------------------------------------------
15	USE par_kind, ONLY : wp ! Precision variables
16	USE in_out_manager ! I/O manager
17	USE obs_profiles_def ! Definitions for storage arrays for profiles
18	USE obs_surf_def ! Definitions for storage arrays for surface data
19	USE obs_mpp, ONLY : & ! MPP support routines for observation diagnostics
20	& obs_mpp_sum_integer, &
21	& obs_mpp_sum_integers
22	USE obs_inter_sup ! Interpolation support
23	USE obs_oper ! Observation operators
24	USE obs_group_def, ONLY : & ! Observation variable information
25	& cobsname_uvel, &
26	& cobsname_vvel, &
27	& imaxavtypes
28	USE lib_mpp, ONLY : ctl_warn, ctl_stop
29	USE bdy_oce, ONLY : & ! Boundary information
30	idx_bdy, nb_bdy, ln_bdy
31
32	IMPLICIT NONE
33	PRIVATE
34
35	PUBLIC obs_pre_prof ! First level check and screening of profile obs
36	PUBLIC obs_pre_surf ! First level check and screening of surface obs
37	PUBLIC calc_month_len ! Calculate the number of days in the months of a year
38
39	!!----------------------------------------------------------------------
40	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
41	!! $Id$
42	!! Software governed by the CeCILL license (see ./LICENSE)
43	!!----------------------------------------------------------------------
44
45
46	CONTAINS
47
48	SUBROUTINE obs_pre_surf( surfdata, surfdataqc, &
49	& kpi, kpj, &
50	& zmask, pglam, pgphi, &
51	& ld_nea, ld_bound_reject, &
52	& kqc_cutoff )
53	!!----------------------------------------------------------------------
54	!! * ROUTINE obs_pre_surf *
55	!!
56	!! ** Purpose : First level check and screening of surface observations
57	!!
58	!! ** Method : First level check and screening of surface observations
59	!!
60	!! ** Action :
61	!!
62	!! References :
63	!!
64	!! History :
65	!! ! 2007-03 (A. Weaver, K. Mogensen) Original
66	!! ! 2007-06 (K. Mogensen et al) Reject obs. near land.
67	!! ! 2015-02 (M. Martin) Combined routine for surface types.
68	!!----------------------------------------------------------------------
69	!! * Modules used
70	USE par_oce ! Ocean parameters
71	USE dom_oce, ONLY : glamt, gphit, tmask, nproc ! Geographical information
72	!! * Arguments
73	TYPE(obs_surf), INTENT(INOUT) :: surfdata ! Full set of surface data
74	TYPE(obs_surf), INTENT(INOUT) :: surfdataqc ! Subset of surface data not failing screening
75	INTEGER, INTENT(IN) :: kpi, kpj ! Local domain sizes
76	REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,surfdata%nvar) :: &
77	& zmask
78	REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,surfdata%nvar) :: &
79	& pglam, &
80	& pgphi
81	LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land
82	LOGICAL, INTENT(IN) :: ld_bound_reject ! Switch for rejecting obs near the boundary
83	INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff ! cut off for QC value
84	!! * Local declarations
85	INTEGER :: iqc_cutoff = 255 ! cut off for QC value
86	INTEGER :: iyea0 ! Initial date
87	INTEGER :: imon0 ! - (year, month, day, hour, minute)
88	INTEGER :: iday0
89	INTEGER :: ihou0
90	INTEGER :: imin0
91	INTEGER :: icycle ! Current assimilation cycle
92	! Counters for observations that are
93	INTEGER :: iotdobs ! - outside time domain
94	INTEGER, DIMENSION(surfdata%nvar) :: iosdsobs ! - outside space domain
95	INTEGER, DIMENSION(surfdata%nvar) :: ilansobs ! - within a model land cell
96	INTEGER, DIMENSION(surfdata%nvar) :: inlasobs ! - close to land
97	INTEGER, DIMENSION(surfdata%nvar) :: ibdysobs ! - close to open boundary
98	INTEGER :: igrdobs ! - fail the grid search
99	! Global counters for observations that
100	INTEGER :: iotdobsmpp ! - outside time domain
101	INTEGER, DIMENSION(surfdata%nvar) :: iosdsobsmpp ! - outside space domain
102	INTEGER, DIMENSION(surfdata%nvar) :: ilansobsmpp ! - within a model land cell
103	INTEGER, DIMENSION(surfdata%nvar) :: inlasobsmpp ! - close to land
104	INTEGER, DIMENSION(surfdata%nvar) :: ibdysobsmpp ! - close to open boundary
105	INTEGER :: igrdobsmpp ! - fail the grid search
106
107	LOGICAL, DIMENSION(:), ALLOCATABLE :: &
108	& llvalid ! SLA data selection
109	INTEGER :: jobs ! Obs. loop counter
110	INTEGER :: jvar ! Variable loop counter
111	INTEGER :: jstp ! Time loop variable
112	INTEGER :: inrc ! Time index variable
113	CHARACTER(LEN=256) :: cout1 ! Diagnostic output line
114	!!----------------------------------------------------------------------
115
116	IF(lwp) WRITE(numout,*) 'obs_pre_surf : Preparing the surface observations...'
117	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
118
119	! Initial date initialization (year, month, day, hour, minute)
120	iyea0 = ndate0 / 10000
121	imon0 = ( ndate0 - iyea0 * 10000 ) / 100
122	iday0 = ndate0 - iyea0 * 10000 - imon0 * 100
123	ihou0 = nn_time0 / 100
124	imin0 = ( nn_time0 - ihou0 * 100 )
125
126	icycle = nn_no ! Assimilation cycle
127
128	! Diagnostic counters for various failures.
129
130	iotdobs = 0
131	igrdobs = 0
132	iosdsobs(:) = 0
133	ilansobs(:) = 0
134	inlasobs(:) = 0
135	ibdysobs(:) = 0
136
137	! Set QC cutoff to optional value if provided
138	IF ( PRESENT(kqc_cutoff) ) iqc_cutoff = kqc_cutoff
139
140	! -----------------------------------------------------------------------
141	! Find time coordinate for surface data
142	! -----------------------------------------------------------------------
143
144	CALL obs_coo_tim( icycle, &
145	& iyea0, imon0, iday0, ihou0, imin0, &
146	& surfdata%nsurf, surfdata%nyea, surfdata%nmon, &
147	& surfdata%nday, surfdata%nhou, surfdata%nmin, &
148	& surfdata%nqc, surfdata%mstp, iotdobs )
149
150	CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
151
152	! -----------------------------------------------------------------------
153	! Check for surface data failing the grid search
154	! -----------------------------------------------------------------------
155
156	DO jvar = 1, surfdata%nvar
157	CALL obs_coo_grd( surfdata%nsurf, surfdata%mi(:,jvar), &
158	& surfdata%mj(:,jvar), surfdata%nqc, igrdobs )
159	END DO
160
161	CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
162
163	! -----------------------------------------------------------------------
164	! Check for land points.
165	! -----------------------------------------------------------------------
166
167	DO jvar = 1, surfdata%nvar
168	CALL obs_coo_spc_2d( surfdata%nsurf, &
169	& jpi, jpj, &
170	& surfdata%mi(:,jvar), surfdata%mj(:,jvar), &
171	& surfdata%rlam, surfdata%rphi, &
172	& pglam(:,:,jvar), pgphi(:,:,jvar), &
173	& zmask(:,:,jvar), surfdata%nqc, &
174	& iosdsobs(jvar), ilansobs(jvar), &
175	& inlasobs(jvar), ld_nea, &
176	& ibdysobs(jvar), ld_bound_reject, &
177	& iqc_cutoff )
178	CALL obs_mpp_sum_integer( iosdsobs(jvar), iosdsobsmpp(jvar) )
179	CALL obs_mpp_sum_integer( ilansobs(jvar), ilansobsmpp(jvar) )
180	CALL obs_mpp_sum_integer( inlasobs(jvar), inlasobsmpp(jvar) )
181	CALL obs_mpp_sum_integer( ibdysobs(jvar), ibdysobsmpp(jvar) )
182	END DO
183
184	! -----------------------------------------------------------------------
185	! Copy useful data from the surfdata data structure to
186	! the surfdataqc data structure
187	! -----------------------------------------------------------------------
188
189	! Allocate the selection arrays
190
191	ALLOCATE( llvalid(surfdata%nsurf) )
192
193	! We want all data which has qc flags <= iqc_cutoff
194
195	llvalid(:) = ( surfdata%nqc(:) <= iqc_cutoff )
196
197	! The actual copying
198
199	CALL obs_surf_compress( surfdata, surfdataqc, .TRUE., numout, &
200	& lvalid=llvalid )
201
202	! Dellocate the selection arrays
203	DEALLOCATE( llvalid )
204
205	! -----------------------------------------------------------------------
206	! Print information about what observations are left after qc
207	! -----------------------------------------------------------------------
208
209	! Update the total observation counter array
210
211	IF(lwp) THEN
212	DO jvar = 1, surfdataqc%nvar
213	IF ( jvar == 1 ) THEN
214	cout1 = TRIM(surfdataqc%cvars(1))
215	ELSE
216	WRITE(cout1,'(A,A1,A)') TRIM(cout1), '/', TRIM(surfdataqc%cvars(jvar))
217	ENDIF
218	END DO
219
220	WRITE(numout,*)
221	WRITE(numout,*) ' '//TRIM(cout1)//' data outside time domain = ', &
222	& iotdobsmpp
223	WRITE(numout,*) ' Remaining '//TRIM(cout1)//' data that failed grid search = ', &
224	& igrdobsmpp
225
226	DO jvar = 1, surfdataqc%nvar
227	WRITE(numout,*) ' Remaining '//surfdataqc%cvars(jvar)//' data outside space domain = ', &
228	& iosdsobsmpp(jvar)
229	WRITE(numout,*) ' Remaining '//surfdataqc%cvars(jvar)//' data at land points = ', &
230	& ilansobsmpp(jvar)
231	IF (ld_nea) THEN
232	WRITE(numout,*) ' Remaining '//surfdataqc%cvars(jvar)//' data near land points (removed) = ', &
233	& inlasobsmpp(jvar)
234	ELSE
235	WRITE(numout,*) ' Remaining '//surfdataqc%cvars(jvar)//' data near land points (kept) = ', &
236	& inlasobsmpp(jvar)
237	ENDIF
238	WRITE(numout,*) ' Remaining '//surfdataqc%cvars(jvar)//' data near open boundary (removed) = ', &
239	& ibdysobsmpp(jvar)
240	END DO
241	WRITE(numout,*) ' '//TRIM(cout1)//' data accepted = ', &
242	& surfdataqc%nsurfmpp
243
244	WRITE(numout,*)
245	WRITE(numout,*) ' Number of observations per time step :'
246	WRITE(numout,*)
247	WRITE(numout,'(10X,A,10X,A)') 'Time step', TRIM(cout1)
248	WRITE(numout,'(10X,A,5X,A)') '---------', '-----------------'
249	CALL FLUSH(numout)
250	ENDIF
251
252	DO jobs = 1, surfdataqc%nsurf
253	inrc = surfdataqc%mstp(jobs) + 2 - nit000
254	surfdataqc%nsstp(inrc) = surfdataqc%nsstp(inrc) + 1
255	END DO
256
257	CALL obs_mpp_sum_integers( surfdataqc%nsstp, surfdataqc%nsstpmpp, &
258	& nitend - nit000 + 2 )
259
260	IF ( lwp ) THEN
261	DO jstp = nit000 - 1, nitend
262	inrc = jstp - nit000 + 2
263	WRITE(numout,1999) jstp, surfdataqc%nsstpmpp(inrc)
264	CALL FLUSH(numout)
265	END DO
266	ENDIF
267
268	1999 FORMAT(10X,I9,5X,I17)
269
270	END SUBROUTINE obs_pre_surf
271
272
273	SUBROUTINE obs_pre_prof( profdata, prodatqc, ld_var, &
274	& kpi, kpj, kpk, &
275	& zmask, pglam, pgphi, &
276	& ld_nea, ld_bound_reject, kdailyavtypes, kqc_cutoff )
277
278	!!----------------------------------------------------------------------
279	!! * ROUTINE obs_pre_prof *
280	!!
281	!! ** Purpose : First level check and screening of profiles
282	!!
283	!! ** Method : First level check and screening of profiles
284	!!
285	!! History :
286	!! ! 2007-06 (K. Mogensen) original : T and S profile data
287	!! ! 2008-09 (M. Valdivieso) : TAO velocity data
288	!! ! 2009-01 (K. Mogensen) : New feedback stricture
289	!! ! 2015-02 (M. Martin) : Combined profile routine.
290	!!
291	!!----------------------------------------------------------------------
292	!! * Modules used
293	USE par_oce ! Ocean parameters
294	USE dom_oce, ONLY : & ! Geographical information
295	& gdept_1d, &
296	& nproc
297
298	!! * Arguments
299	TYPE(obs_prof), INTENT(INOUT) :: profdata ! Full set of profile data
300	TYPE(obs_prof), INTENT(INOUT) :: prodatqc ! Subset of profile data not failing screening
301	LOGICAL, DIMENSION(profdata%nvar), INTENT(IN) :: &
302	& ld_var ! Observed variables switches
303	LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land
304	LOGICAL, INTENT(IN) :: ld_bound_reject ! Switch for rejecting observations near the boundary
305	INTEGER, INTENT(IN) :: kpi, kpj, kpk ! Local domain sizes
306	INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
307	& kdailyavtypes ! Types for daily averages
308	REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,kpk,profdata%nvar) :: &
309	& zmask
310	REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,profdata%nvar) :: &
311	& pglam, &
312	& pgphi
313	INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff ! cut off for QC value
314
315	!! * Local declarations
316	INTEGER :: iqc_cutoff = 255 ! cut off for QC value
317	INTEGER :: iyea0 ! Initial date
318	INTEGER :: imon0 ! - (year, month, day, hour, minute)
319	INTEGER :: iday0
320	INTEGER :: ihou0
321	INTEGER :: imin0
322	INTEGER :: icycle ! Current assimilation cycle
323	! Counters for observations that are
324	INTEGER :: iotdobs ! - outside time domain
325	INTEGER, DIMENSION(profdata%nvar) :: iosdvobs ! - outside space domain
326	INTEGER, DIMENSION(profdata%nvar) :: ilanvobs ! - within a model land cell
327	INTEGER, DIMENSION(profdata%nvar) :: inlavobs ! - close to land
328	INTEGER, DIMENSION(profdata%nvar) :: ibdyvobs ! - boundary
329	INTEGER :: igrdobs ! - fail the grid search
330	INTEGER :: iuvchku ! - reject UVEL if VVEL rejected
331	INTEGER :: iuvchkv ! - reject VVEL if UVEL rejected
332	! Global counters for observations that are
333	INTEGER :: iotdobsmpp ! - outside time domain
334	INTEGER, DIMENSION(profdata%nvar) :: iosdvobsmpp ! - outside space domain
335	INTEGER, DIMENSION(profdata%nvar) :: ilanvobsmpp ! - within a model land cell
336	INTEGER, DIMENSION(profdata%nvar) :: inlavobsmpp ! - close to land
337	INTEGER, DIMENSION(profdata%nvar) :: ibdyvobsmpp ! - boundary
338	INTEGER :: igrdobsmpp ! - fail the grid search
339	INTEGER :: iuvchkumpp ! - reject UVEL if VVEL rejected
340	INTEGER :: iuvchkvmpp ! - reject VVEL if UVEL rejected
341	INTEGER :: iuvar ! - UVEL index
342	INTEGER :: ivvar ! - VVEL index
343	TYPE(obs_prof_valid) :: llvalid ! Profile selection
344	TYPE(obs_prof_valid), DIMENSION(profdata%nvar) :: &
345	& llvvalid ! var selection
346	INTEGER :: jvar ! Variable loop variable
347	INTEGER :: jobs ! Obs. loop variable
348	INTEGER :: jstp ! Time loop variable
349	INTEGER :: inrc ! Time index variable
350	CHARACTER(LEN=256) :: cout1 ! Diagnostic output line
351	CHARACTER(LEN=256) :: cout2 ! Diagnostic output line
352	!!----------------------------------------------------------------------
353
354	IF(lwp) WRITE(numout,*)'obs_pre_prof: Preparing the profile data...'
355	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
356
357	! Initial date initialization (year, month, day, hour, minute)
358	iyea0 = ndate0 / 10000
359	imon0 = ( ndate0 - iyea0 * 10000 ) / 100
360	iday0 = ndate0 - iyea0 * 10000 - imon0 * 100
361	ihou0 = nn_time0 / 100
362	imin0 = ( nn_time0 - ihou0 * 100 )
363
364	icycle = nn_no ! Assimilation cycle
365
366	! Diagnostic counters for various failures.
367
368	iotdobs = 0
369	igrdobs = 0
370	iosdvobs(:) = 0
371	ilanvobs(:) = 0
372	inlavobs(:) = 0
373	ibdyvobs(:) = 0
374	iuvchku = 0
375	iuvchkv = 0
376
377
378	! Set QC cutoff to optional value if provided
379	IF ( PRESENT(kqc_cutoff) ) iqc_cutoff = kqc_cutoff
380
381	! -----------------------------------------------------------------------
382	! Find time coordinate for profiles
383	! -----------------------------------------------------------------------
384
385	IF ( PRESENT(kdailyavtypes) ) THEN
386	CALL obs_coo_tim_prof( icycle, &
387	& iyea0, imon0, iday0, ihou0, imin0, &
388	& profdata%nprof, profdata%nyea, profdata%nmon, &
389	& profdata%nday, profdata%nhou, profdata%nmin, &
390	& profdata%ntyp, profdata%nqc, profdata%mstp, &
391	& iotdobs, kdailyavtypes = kdailyavtypes, &
392	& kqc_cutoff = iqc_cutoff )
393	ELSE
394	CALL obs_coo_tim_prof( icycle, &
395	& iyea0, imon0, iday0, ihou0, imin0, &
396	& profdata%nprof, profdata%nyea, profdata%nmon, &
397	& profdata%nday, profdata%nhou, profdata%nmin, &
398	& profdata%ntyp, profdata%nqc, profdata%mstp, &
399	& iotdobs, kqc_cutoff = iqc_cutoff )
400	ENDIF
401
402	CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
403
404	! -----------------------------------------------------------------------
405	! Check for profiles failing the grid search
406	! -----------------------------------------------------------------------
407
408	DO jvar = 1, profdata%nvar
409	CALL obs_coo_grd( profdata%nprof, profdata%mi(:,jvar), profdata%mj(:,jvar), &
410	& profdata%nqc, igrdobs )
411	END DO
412
413	CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
414
415	! -----------------------------------------------------------------------
416	! Reject all observations for profiles with nqc > iqc_cutoff
417	! -----------------------------------------------------------------------
418
419	CALL obs_pro_rej( profdata, kqc_cutoff = iqc_cutoff )
420
421	! -----------------------------------------------------------------------
422	! Check for land points. This includes points below the model
423	! bathymetry so this is done for every point in the profile
424	! -----------------------------------------------------------------------
425
426	DO jvar = 1, profdata%nvar
427	CALL obs_coo_spc_3d( profdata%nprof, profdata%nvprot(jvar), &
428	& profdata%npvsta(:,jvar), profdata%npvend(:,jvar), &
429	& jpi, jpj, &
430	& jpk, &
431	& profdata%mi, profdata%mj, &
432	& profdata%var(jvar)%mvk, &
433	& profdata%rlam, profdata%rphi, &
434	& profdata%var(jvar)%vdep, &
435	& pglam(:,:,jvar), pgphi(:,:,jvar), &
436	& gdept_1d, zmask(:,:,:,jvar), &
437	& profdata%nqc, profdata%var(jvar)%nvqc, &
438	& iosdvobs(jvar), ilanvobs(jvar), &
439	& inlavobs(jvar), ld_nea, &
440	& ibdyvobs(jvar), ld_bound_reject, &
441	& iqc_cutoff )
442
443	CALL obs_mpp_sum_integer( iosdvobs(jvar), iosdvobsmpp(jvar) )
444	CALL obs_mpp_sum_integer( ilanvobs(jvar), ilanvobsmpp(jvar) )
445	CALL obs_mpp_sum_integer( inlavobs(jvar), inlavobsmpp(jvar) )
446	CALL obs_mpp_sum_integer( ibdyvobs(jvar), ibdyvobsmpp(jvar) )
447	END DO
448
449	! -----------------------------------------------------------------------
450	! Reject u if v is rejected and vice versa
451	! -----------------------------------------------------------------------
452
453	iuvar = 0
454	ivvar = 0
455	DO jvar = 1,profdata%nvar
456	IF ( TRIM(profdata%cvars(jvar)) == cobsname_uvel ) THEN
457	iuvar = jvar
458	ELSEIF ( TRIM(profdata%cvars(jvar)) == cobsname_vvel ) THEN
459	ivvar = jvar
460	ENDIF
461	END DO
462	IF ( (iuvar > 0) .AND. (ivvar > 0) ) THEN
463	CALL obs_uv_rej( profdata, iuvchku, iuvchkv, iqc_cutoff, iuvar, ivvar )
464	CALL obs_mpp_sum_integer( iuvchku, iuvchkumpp )
465	CALL obs_mpp_sum_integer( iuvchkv, iuvchkvmpp )
466	ENDIF
467
468	! -----------------------------------------------------------------------
469	! Copy useful data from the profdata data structure to
470	! the prodatqc data structure
471	! -----------------------------------------------------------------------
472
473	! Allocate the selection arrays
474
475	ALLOCATE( llvalid%luse(profdata%nprof) )
476	DO jvar = 1,profdata%nvar
477	ALLOCATE( llvvalid(jvar)%luse(profdata%nvprot(jvar)) )
478	END DO
479
480	! We want all data which has qc flags <= iqc_cutoff
481
482	llvalid%luse(:) = ( profdata%nqc(:) <= iqc_cutoff )
483	DO jvar = 1,profdata%nvar
484	llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= iqc_cutoff )
485	END DO
486
487	! The actual copying
488
489	CALL obs_prof_compress( profdata, prodatqc, .TRUE., numout, &
490	& lvalid=llvalid, lvvalid=llvvalid )
491
492	! Dellocate the selection arrays
493	DEALLOCATE( llvalid%luse )
494	DO jvar = 1,profdata%nvar
495	DEALLOCATE( llvvalid(jvar)%luse )
496	END DO
497
498	! -----------------------------------------------------------------------
499	! Print information about what observations are left after qc
500	! -----------------------------------------------------------------------
501
502	! Update the total observation counter array
503
504	IF(lwp) THEN
505
506	WRITE(numout,*)
507	WRITE(numout,*) ' Profiles outside time domain = ', &
508	& iotdobsmpp
509	WRITE(numout,*) ' Remaining profiles that failed grid search = ', &
510	& igrdobsmpp
511	DO jvar = 1, profdata%nvar
512	WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data outside space domain = ', &
513	& iosdvobsmpp(jvar)
514	WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data at land points = ', &
515	& ilanvobsmpp(jvar)
516	IF (ld_nea) THEN
517	WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data near land points (removed) = ', &
518	& inlavobsmpp(jvar)
519	ELSE
520	WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data near land points (kept) = ', &
521	& inlavobsmpp(jvar)
522	ENDIF
523	IF ( (iuvar > 0) .AND. (ivvar > 0) ) THEN
524	IF ( TRIM(profdata%cvars(jvar)) == cobsname_uvel ) THEN
525	WRITE(numout,*) ' U observation rejected since V rejected = ', &
526	& iuvchku
527	ELSE IF ( TRIM(profdata%cvars(jvar)) == cobsname_vvel ) THEN
528	WRITE(numout,*) ' V observation rejected since U rejected = ', &
529	& iuvchkv
530	ENDIF
531	ENDIF
532	WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data near open boundary (removed) = ', &
533	& ibdyvobsmpp(jvar)
534	WRITE(numout,*) ' '//prodatqc%cvars(jvar)//' data accepted = ', &
535	& prodatqc%nvprotmpp(jvar)
536	END DO
537
538	WRITE(numout,*)
539	WRITE(numout,*) ' Number of observations per time step :'
540	WRITE(numout,*)
541	WRITE(cout1,'(10X,A9,5X,A8)') 'Time step', 'Profiles'
542	WRITE(cout2,'(10X,A9,5X,A8)') '---------', '--------'
543	DO jvar = 1, prodatqc%nvar
544	WRITE(cout1,'(A,5X,A11)') TRIM(cout1), TRIM(prodatqc%cvars(jvar))
545	WRITE(cout2,'(A,5X,A11)') TRIM(cout2), '-----------'
546	END DO
547	WRITE(numout,*) cout1
548	WRITE(numout,*) cout2
549	ENDIF
550
551	DO jobs = 1, prodatqc%nprof
552	inrc = prodatqc%mstp(jobs) + 2 - nit000
553	prodatqc%npstp(inrc) = prodatqc%npstp(inrc) + 1
554	DO jvar = 1, prodatqc%nvar
555	IF ( prodatqc%npvend(jobs,jvar) > 0 ) THEN
556	prodatqc%nvstp(inrc,jvar) = prodatqc%nvstp(inrc,jvar) + &
557	& ( prodatqc%npvend(jobs,jvar) - &
558	& prodatqc%npvsta(jobs,jvar) + 1 )
559	ENDIF
560	END DO
561	END DO
562
563
564	CALL obs_mpp_sum_integers( prodatqc%npstp, prodatqc%npstpmpp, &
565	& nitend - nit000 + 2 )
566	DO jvar = 1, prodatqc%nvar
567	CALL obs_mpp_sum_integers( prodatqc%nvstp(:,jvar), &
568	& prodatqc%nvstpmpp(:,jvar), &
569	& nitend - nit000 + 2 )
570	END DO
571
572	IF ( lwp ) THEN
573	DO jstp = nit000 - 1, nitend
574	inrc = jstp - nit000 + 2
575	WRITE(cout1,'(10X,I9,5X,I8)') jstp, prodatqc%npstpmpp(inrc)
576	DO jvar = 1, prodatqc%nvar
577	WRITE(cout1,'(A,5X,I11)') TRIM(cout1), prodatqc%nvstpmpp(inrc,jvar)
578	END DO
579	WRITE(numout,*) cout1
580	END DO
581	ENDIF
582
583	END SUBROUTINE obs_pre_prof
584
585	SUBROUTINE obs_coo_tim( kcycle, &
586	& kyea0, kmon0, kday0, khou0, kmin0, &
587	& kobsno, &
588	& kobsyea, kobsmon, kobsday, kobshou, kobsmin, &
589	& kobsqc, kobsstp, kotdobs )
590	!!----------------------------------------------------------------------
591	!! * ROUTINE obs_coo_tim *
592	!!
593	!! ** Purpose : Compute the number of time steps to the observation time.
594	!!
595	!! ** Method : For time coordinates ( yea_obs, mon_obs, day_obs,
596	!! hou_obs, min_obs ), this routine locates the time step
597	!! that is closest to this time.
598	!!
599	!! ** Action :
600	!!
601	!! References :
602	!!
603	!! History :
604	!! ! 1997-07 (A. Weaver) Original
605	!! ! 2006-08 (A. Weaver) NEMOVAR migration
606	!! ! 2006-10 (A. Weaver) Cleanup
607	!! ! 2007-01 (K. Mogensen) Rewritten with loop
608	!! ! 2010-05 (D. Lea) Fix in leap year calculation for NEMO vn3.2
609	!!----------------------------------------------------------------------
610	!! * Modules used
611	USE dom_oce, ONLY : & ! Geographical information
612	& rdt
613	USE phycst, ONLY : & ! Physical constants
614	& rday, &
615	& rmmss, &
616	& rhhmm
617	!! * Arguments
618	INTEGER, INTENT(IN) :: kcycle ! Current cycle
619	INTEGER, INTENT(IN) :: kyea0 ! Initial date coordinates
620	INTEGER, INTENT(IN) :: kmon0
621	INTEGER, INTENT(IN) :: kday0
622	INTEGER, INTENT(IN) :: khou0
623	INTEGER, INTENT(IN) :: kmin0
624	INTEGER, INTENT(IN) :: kobsno ! Number of observations
625	INTEGER, INTENT(INOUT) :: kotdobs ! Number of observations failing time check
626	INTEGER, DIMENSION(kobsno), INTENT(IN ) :: &
627	& kobsyea, & ! Observation time coordinates
628	& kobsmon, &
629	& kobsday, &
630	& kobshou, &
631	& kobsmin
632	INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
633	& kobsqc ! Quality control flag
634	INTEGER, DIMENSION(kobsno), INTENT(OUT) :: &
635	& kobsstp ! Number of time steps up to the
636	! observation time
637
638	!! * Local declarations
639	INTEGER :: jyea
640	INTEGER :: jmon
641	INTEGER :: jday
642	INTEGER :: jobs
643	INTEGER :: iyeastr
644	INTEGER :: iyeaend
645	INTEGER :: imonstr
646	INTEGER :: imonend
647	INTEGER :: idaystr
648	INTEGER :: idayend
649	INTEGER :: iskip
650	INTEGER :: idaystp
651	REAL(KIND=wp) :: zminstp
652	REAL(KIND=wp) :: zhoustp
653	REAL(KIND=wp) :: zobsstp
654	INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year
655
656	!-----------------------------------------------------------------------
657	! Initialization
658	!-----------------------------------------------------------------------
659
660	! Intialize the number of time steps per day
661	idaystp = NINT( rday / rdt )
662
663	!---------------------------------------------------------------------
664	! Locate the model time coordinates for interpolation
665	!---------------------------------------------------------------------
666
667	DO jobs = 1, kobsno
668
669	! Initialize the time step counter
670	kobsstp(jobs) = nit000 - 1
671
672	! Flag if observation date is less than the initial date
673
674	IF ( ( kobsyea(jobs) < kyea0 ) &
675	& .OR. ( ( kobsyea(jobs) == kyea0 ) &
676	& .AND. ( kobsmon(jobs) < kmon0 ) ) &
677	& .OR. ( ( kobsyea(jobs) == kyea0 ) &
678	& .AND. ( kobsmon(jobs) == kmon0 ) &
679	& .AND. ( kobsday(jobs) < kday0 ) ) &
680	& .OR. ( ( kobsyea(jobs) == kyea0 ) &
681	& .AND. ( kobsmon(jobs) == kmon0 ) &
682	& .AND. ( kobsday(jobs) == kday0 ) &
683	& .AND. ( kobshou(jobs) < khou0 ) ) &
684	& .OR. ( ( kobsyea(jobs) == kyea0 ) &
685	& .AND. ( kobsmon(jobs) == kmon0 ) &
686	& .AND. ( kobsday(jobs) == kday0 ) &
687	& .AND. ( kobshou(jobs) == khou0 ) &
688	& .AND. ( kobsmin(jobs) <= kmin0 ) ) ) THEN
689	kobsstp(jobs) = -1
690	kobsqc(jobs) = IBSET(kobsqc(jobs),13)
691	kotdobs = kotdobs + 1
692	CYCLE
693	ENDIF
694
695	! Compute the number of time steps to the observation day
696	iyeastr = kyea0
697	iyeaend = kobsyea(jobs)
698
699	!---------------------------------------------------------------------
700	! Year loop
701	!---------------------------------------------------------------------
702	DO jyea = iyeastr, iyeaend
703
704	CALL calc_month_len( jyea, imonth_len )
705
706	imonstr = 1
707	IF ( jyea == kyea0 ) imonstr = kmon0
708	imonend = 12
709	IF ( jyea == kobsyea(jobs) ) imonend = kobsmon(jobs)
710
711	! Month loop
712	DO jmon = imonstr, imonend
713
714	idaystr = 1
715	IF ( ( jmon == kmon0 ) &
716	& .AND. ( jyea == kyea0 ) ) idaystr = kday0
717	idayend = imonth_len(jmon)
718	IF ( ( jmon == kobsmon(jobs) ) &
719	& .AND. ( jyea == kobsyea(jobs) ) ) idayend = kobsday(jobs) - 1
720
721	! Day loop
722	DO jday = idaystr, idayend
723	kobsstp(jobs) = kobsstp(jobs) + idaystp
724	END DO
725
726	END DO
727
728	END DO
729
730	! Add in the number of time steps to the observation minute
731	zminstp = rmmss / rdt
732	zhoustp = rhhmm * zminstp
733
734	zobsstp = REAL( kobsmin(jobs) - kmin0, KIND=wp ) * zminstp &
735	& + REAL( kobshou(jobs) - khou0, KIND=wp ) * zhoustp
736	kobsstp(jobs) = kobsstp(jobs) + NINT( zobsstp )
737
738	! Flag if observation step outside the time window
739	IF ( ( kobsstp(jobs) < ( nit000 - 1 ) ) &
740	& .OR.( kobsstp(jobs) > nitend ) ) THEN
741	kobsqc(jobs) = IBSET(kobsqc(jobs),13)
742	kotdobs = kotdobs + 1
743	CYCLE
744	ENDIF
745
746	END DO
747
748	END SUBROUTINE obs_coo_tim
749
750	SUBROUTINE calc_month_len( iyear, imonth_len )
751	!!----------------------------------------------------------------------
752	!! * ROUTINE calc_month_len *
753	!!
754	!! ** Purpose : Compute the number of days in a months given a year.
755	!!
756	!! ** Method :
757	!!
758	!! ** Action :
759	!!
760	!! History :
761	!! ! 10-05 (D. Lea) New routine based on day_init
762	!!----------------------------------------------------------------------
763
764	INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year
765	INTEGER :: iyear !: year
766
767	! length of the month of the current year (from nleapy, read in namelist)
768	IF ( nleapy < 2 ) THEN
769	imonth_len(:) = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)
770	IF ( nleapy == 1 ) THEN ! we are using calendar with leap years
771	IF ( MOD(iyear, 4) == 0 .AND. ( MOD(iyear, 400) == 0 .OR. MOD(iyear, 100) /= 0 ) ) THEN
772	imonth_len(2) = 29
773	ENDIF
774	ENDIF
775	ELSE
776	imonth_len(:) = nleapy ! all months with nleapy days per year
777	ENDIF
778
779	END SUBROUTINE
780
781	SUBROUTINE obs_coo_tim_prof( kcycle, &
782	& kyea0, kmon0, kday0, khou0, kmin0, &
783	& kobsno, &
784	& kobsyea, kobsmon, kobsday, kobshou, kobsmin, &
785	& ktyp, kobsqc, kobsstp, kotdobs, kdailyavtypes, &
786	& kqc_cutoff )
787	!!----------------------------------------------------------------------
788	!! * ROUTINE obs_coo_tim *
789	!!
790	!! ** Purpose : Compute the number of time steps to the observation time.
791	!!
792	!! ** Method : For time coordinates ( yea_obs, mon_obs, day_obs,
793	!! hou_obs, min_obs ), this routine locates the time step
794	!! that is closest to this time.
795	!!
796	!! ** Action :
797	!!
798	!! References :
799	!!
800	!! History :
801	!! ! 1997-07 (A. Weaver) Original
802	!! ! 2006-08 (A. Weaver) NEMOVAR migration
803	!! ! 2006-10 (A. Weaver) Cleanup
804	!! ! 2007-01 (K. Mogensen) Rewritten with loop
805	!!----------------------------------------------------------------------
806	!! * Modules used
807	!! * Arguments
808	INTEGER, INTENT(IN) :: kcycle ! Current cycle
809	INTEGER, INTENT(IN) :: kyea0 ! Initial date coordinates
810	INTEGER, INTENT(IN) :: kmon0
811	INTEGER, INTENT(IN) :: kday0
812	INTEGER, INTENT(IN) :: khou0
813	INTEGER, INTENT(IN) :: kmin0
814	INTEGER, INTENT(IN) :: kobsno ! Number of observations
815	INTEGER, INTENT(INOUT) :: kotdobs ! Number of observations failing time check
816	INTEGER, DIMENSION(kobsno), INTENT(IN ) :: &
817	& kobsyea, & ! Observation time coordinates
818	& kobsmon, &
819	& kobsday, &
820	& kobshou, &
821	& kobsmin, &
822	& ktyp ! Observation type.
823	INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
824	& kobsqc ! Quality control flag
825	INTEGER, DIMENSION(kobsno), INTENT(OUT) :: &
826	& kobsstp ! Number of time steps up to the
827	! observation time
828	INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
829	& kdailyavtypes ! Types for daily averages
830	INTEGER, OPTIONAL, INTENT(IN) :: kqc_cutoff ! QC cutoff value
831
832	!! * Local declarations
833	INTEGER :: jobs
834	INTEGER :: iqc_cutoff = 255
835
836	!-----------------------------------------------------------------------
837	! Call standard obs_coo_tim
838	!-----------------------------------------------------------------------
839
840	CALL obs_coo_tim( kcycle, &
841	& kyea0, kmon0, kday0, khou0, kmin0, &
842	& kobsno, &
843	& kobsyea, kobsmon, kobsday, kobshou, kobsmin, &
844	& kobsqc, kobsstp, kotdobs )
845
846	!------------------------------------------------------------------------
847	! Always reject daily averaged data (e.g. MRB data (820)) at initial time
848	!------------------------------------------------------------------------
849
850	IF ( PRESENT(kdailyavtypes) ) THEN
851	DO jobs = 1, kobsno
852
853	IF ( kobsqc(jobs) <= iqc_cutoff ) THEN
854
855	IF ( ( kobsstp(jobs) == (nit000 - 1) ).AND.&
856	& ( ANY (kdailyavtypes(:) == ktyp(jobs)) ) ) THEN
857	kobsqc(jobs) = IBSET(kobsqc(jobs),13)
858	kotdobs = kotdobs + 1
859	CYCLE
860	ENDIF
861
862	ENDIF
863	END DO
864	ENDIF
865
866
867	END SUBROUTINE obs_coo_tim_prof
868
869	SUBROUTINE obs_coo_grd( kobsno, kobsi, kobsj, kobsqc, kgrdobs )
870	!!----------------------------------------------------------------------
871	!! * ROUTINE obs_coo_grd *
872	!!
873	!! ** Purpose : Verify that the grid search has not failed
874	!!
875	!! ** Method : The previously computed i,j indeces are checked
876	!!
877	!! ** Action :
878	!!
879	!! References :
880	!!
881	!! History :
882	!! ! 2007-01 (K. Mogensen) Original
883	!!----------------------------------------------------------------------
884	INTEGER, INTENT(IN) :: kobsno ! Number of observations
885	INTEGER, DIMENSION(kobsno), INTENT(IN ) :: &
886	& kobsi, & ! i,j indeces previously computed
887	& kobsj
888	INTEGER, INTENT(INOUT) :: kgrdobs ! Number of observations failing the check
889	INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
890	& kobsqc ! Quality control flag
891
892	!! * Local declarations
893	INTEGER :: jobs ! Loop variable
894
895	! Flag if the grid search failed
896
897	DO jobs = 1, kobsno
898	IF ( ( kobsi(jobs) <= 0 ) .AND. ( kobsj(jobs) <= 0 ) ) THEN
899	kobsqc(jobs) = IBSET(kobsqc(jobs),12)
900	kgrdobs = kgrdobs + 1
901	ENDIF
902	END DO
903
904	END SUBROUTINE obs_coo_grd
905
906	SUBROUTINE obs_coo_spc_2d( kobsno, kpi, kpj, &
907	& kobsi, kobsj, pobslam, pobsphi, &
908	& plam, pphi, pmask, &
909	& kobsqc, kosdobs, klanobs, &
910	& knlaobs,ld_nea, &
911	& kbdyobs,ld_bound_reject, &
912	& kqc_cutoff )
913	!!----------------------------------------------------------------------
914	!! * ROUTINE obs_coo_spc_2d *
915	!!
916	!! ** Purpose : Check for points outside the domain and land points
917	!!
918	!! ** Method : Remove the observations that are outside the model space
919	!! and time domain or located within model land cells.
920	!!
921	!! ** Action :
922	!!
923	!! History : 2007-03 (A. Weaver, K. Mogensen) Original
924	!! ! 2007-06 (K. Mogensen et al) Reject obs. near land.
925	!!----------------------------------------------------------------------
926	INTEGER , INTENT(in ) :: kobsno ! Total number of observations
927	INTEGER , INTENT(in ) :: kpi , kpj ! Number of grid points in (i,j)
928	INTEGER , INTENT(in ), DIMENSION(kobsno) :: kobsi , kobsj ! Observation (i,j) coordinates
929	REAL(wp), INTENT(in ), DIMENSION(kobsno) :: pobslam, pobsphi ! Observation (lon,lat) coordinates
930	REAL(wp), INTENT(in ), DIMENSION(kpi,kpj) :: plam , pphi ! Model (lon,lat) coordinates
931	REAL(wp), INTENT(in ), DIMENSION(kpi,kpj) :: pmask ! Land mask array
932	INTEGER , INTENT(inout), DIMENSION(kobsno) :: kobsqc ! Observation quality control
933	INTEGER , INTENT(inout) :: kosdobs ! Observations outside space domain
934	INTEGER , INTENT(inout) :: klanobs ! Observations within a model land cell
935	INTEGER , INTENT(inout) :: knlaobs ! Observations near land
936	INTEGER , INTENT(inout) :: kbdyobs ! Observations near boundary
937	LOGICAL , INTENT(in ) :: ld_nea ! Flag observations near land
938	LOGICAL , INTENT(in ) :: ld_bound_reject ! Flag observations near open boundary
939	INTEGER , INTENT(in ) :: kqc_cutoff ! Cutoff QC value
940	!
941	REAL(KIND=wp), DIMENSION(2,2,kobsno) :: zgmsk ! Grid mask
942	REAL(KIND=wp), DIMENSION(2,2,kobsno) :: zbmsk ! Boundary mask
943	REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask
944	REAL(KIND=wp), DIMENSION(2,2,kobsno) :: zglam, zgphi ! Model Lon/lat at grid points
945	INTEGER , DIMENSION(2,2,kobsno) :: igrdi, igrdj ! Grid i,j
946	LOGICAL :: lgridobs ! Is observation on a model grid point.
947	INTEGER :: iig, ijg ! i,j of observation on model grid point.
948	INTEGER :: jobs, ji, jj
949	!!----------------------------------------------------------------------
950
951	! Get grid point indices
952
953	DO jobs = 1, kobsno
954
955	! For invalid points use 2,2
956
957	IF ( kobsqc(jobs) >= kqc_cutoff ) THEN
958
959	igrdi(1,1,jobs) = 1
960	igrdj(1,1,jobs) = 1
961	igrdi(1,2,jobs) = 1
962	igrdj(1,2,jobs) = 2
963	igrdi(2,1,jobs) = 2
964	igrdj(2,1,jobs) = 1
965	igrdi(2,2,jobs) = 2
966	igrdj(2,2,jobs) = 2
967
968	ELSE
969
970	igrdi(1,1,jobs) = kobsi(jobs)-1
971	igrdj(1,1,jobs) = kobsj(jobs)-1
972	igrdi(1,2,jobs) = kobsi(jobs)-1
973	igrdj(1,2,jobs) = kobsj(jobs)
974	igrdi(2,1,jobs) = kobsi(jobs)
975	igrdj(2,1,jobs) = kobsj(jobs)-1
976	igrdi(2,2,jobs) = kobsi(jobs)
977	igrdj(2,2,jobs) = kobsj(jobs)
978
979	ENDIF
980
981	END DO
982
983	IF (ln_bdy) THEN
984	! Create a mask grid points in boundary rim
985	IF (ld_bound_reject) THEN
986	zbdymask(:,:) = 1.0_wp
987	DO ji = 1, nb_bdy
988	DO jj = 1, idx_bdy(ji)%nblen(1)
989	zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp
990	ENDDO
991	ENDDO
992
993	CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, zbdymask, zbmsk )
994	ENDIF
995	ENDIF
996
997
998	CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, pmask, zgmsk )
999	CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, plam, zglam )
1000	CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, pphi, zgphi )
1001
1002	DO jobs = 1, kobsno
1003
1004	! Skip bad observations
1005	IF ( kobsqc(jobs) >= kqc_cutoff ) CYCLE
1006
1007	! Flag if the observation falls outside the model spatial domain
1008	IF ( ( pobslam(jobs) < -180. ) &
1009	& .OR. ( pobslam(jobs) > 180. ) &
1010	& .OR. ( pobsphi(jobs) < -90. ) &
1011	& .OR. ( pobsphi(jobs) > 90. ) ) THEN
1012	kobsqc(jobs) = IBSET(kobsqc(jobs),11)
1013	kosdobs = kosdobs + 1
1014	CYCLE
1015	ENDIF
1016
1017	! Flag if the observation falls with a model land cell
1018	IF ( SUM( zgmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN
1019	kobsqc(jobs) = IBSET(kobsqc(jobs),10)
1020	klanobs = klanobs + 1
1021	CYCLE
1022	ENDIF
1023
1024	! Check if this observation is on a grid point
1025
1026	lgridobs = .FALSE.
1027	iig = -1
1028	ijg = -1
1029	DO jj = 1, 2
1030	DO ji = 1, 2
1031	IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) &
1032	& .AND. &
1033	& ( ABS( MOD( zglam(ji,jj,jobs) - pobslam(jobs),360.0) ) &
1034	& < 1.0e-6_wp ) ) THEN
1035	lgridobs = .TRUE.
1036	iig = ji
1037	ijg = jj
1038	ENDIF
1039	END DO
1040	END DO
1041
1042	! For observations on the grid reject them if their are at
1043	! a masked point
1044
1045	IF (lgridobs) THEN
1046	IF (zgmsk(iig,ijg,jobs) == 0.0_wp ) THEN
1047	kobsqc(jobs) = IBSET(kobsqc(jobs),10)
1048	klanobs = klanobs + 1
1049	CYCLE
1050	ENDIF
1051	ENDIF
1052
1053	! Flag if the observation falls is close to land
1054	IF ( MINVAL( zgmsk(1:2,1:2,jobs) ) == 0.0_wp) THEN
1055	knlaobs = knlaobs + 1
1056	IF (ld_nea) THEN
1057	kobsqc(jobs) = IBSET(kobsqc(jobs),9)
1058	CYCLE
1059	ENDIF
1060	ENDIF
1061
1062	IF (ln_bdy) THEN
1063	! Flag if the observation falls close to the boundary rim
1064	IF (ld_bound_reject) THEN
1065	IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN
1066	kobsqc(jobs) = IBSET(kobsqc(jobs),8)
1067	kbdyobs = kbdyobs + 1
1068	CYCLE
1069	ENDIF
1070	! for observations on the grid...
1071	IF (lgridobs) THEN
1072	IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN
1073	kobsqc(jobs) = IBSET(kobsqc(jobs),8)
1074	kbdyobs = kbdyobs + 1
1075	CYCLE
1076	ENDIF
1077	ENDIF
1078	ENDIF
1079	ENDIF
1080	!
1081	END DO
1082	!
1083	END SUBROUTINE obs_coo_spc_2d
1084
1085
1086	SUBROUTINE obs_coo_spc_3d( kprofno, kobsno, kpstart, kpend, &
1087	& kpi, kpj, kpk, &
1088	& kobsi, kobsj, kobsk, &
1089	& pobslam, pobsphi, pobsdep, &
1090	& plam, pphi, pdep, pmask, &
1091	& kpobsqc, kobsqc, kosdobs, &
1092	& klanobs, knlaobs, ld_nea, &
1093	& kbdyobs, ld_bound_reject, &
1094	& kqc_cutoff )
1095	!!----------------------------------------------------------------------
1096	!! * ROUTINE obs_coo_spc_3d *
1097	!!
1098	!! ** Purpose : Check for points outside the domain and land points
1099	!! Reset depth of observation above highest model level
1100	!! to the value of highest model level
1101	!!
1102	!! ** Method : Remove the observations that are outside the model space
1103	!! and time domain or located within model land cells.
1104	!!
1105	!! NB. T and S profile observations lying between the ocean
1106	!! surface and the depth of the first model T point are
1107	!! assigned a depth equal to that of the first model T pt.
1108	!!
1109	!! ** Action :
1110	!!
1111	!! History :
1112	!! ! 2007-01 (K. Mogensen) Rewrite of parts of obs_scr
1113	!! ! 2007-06 (K. Mogensen et al) Reject obs. near land.
1114	!!----------------------------------------------------------------------
1115	!! * Modules used
1116	USE dom_oce, ONLY : & ! Geographical information
1117	& gdepw_1d, &
1118	& gdepw_0, &
1119	& gdepw_n, &
1120	& gdept_n, &
1121	& ln_zco, &
1122	& ln_zps, &
1123	& mbkt
1124
1125	!! * Arguments
1126	INTEGER, INTENT(IN) :: kprofno ! Number of profiles
1127	INTEGER, INTENT(IN) :: kobsno ! Total number of observations
1128	INTEGER, INTENT(IN) :: kpi ! Number of grid points in (i,j,k)
1129	INTEGER, INTENT(IN) :: kpj
1130	INTEGER, INTENT(IN) :: kpk
1131	INTEGER, DIMENSION(kprofno), INTENT(IN) :: &
1132	& kpstart, & ! Start of individual profiles
1133	& kpend ! End of individual profiles
1134	INTEGER, DIMENSION(kprofno), INTENT(IN) :: &
1135	& kobsi, & ! Observation (i,j) coordinates
1136	& kobsj
1137	INTEGER, DIMENSION(kobsno), INTENT(IN) :: &
1138	& kobsk ! Observation k coordinate
1139	REAL(KIND=wp), DIMENSION(kprofno), INTENT(IN) :: &
1140	& pobslam, & ! Observation (lon,lat) coordinates
1141	& pobsphi
1142	REAL(KIND=wp), DIMENSION(kobsno), INTENT(INOUT) :: &
1143	& pobsdep ! Observation depths
1144	REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) :: &
1145	& plam, pphi ! Model (lon,lat) coordinates
1146	REAL(KIND=wp), DIMENSION(kpk), INTENT(IN) :: &
1147	& pdep ! Model depth coordinates
1148	REAL(KIND=wp), DIMENSION(kpi,kpj,kpk), INTENT(IN) :: &
1149	& pmask ! Land mask array
1150	INTEGER, DIMENSION(kprofno), INTENT(INOUT) :: &
1151	& kpobsqc ! Profile quality control
1152	INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
1153	& kobsqc ! Observation quality control
1154	INTEGER, INTENT(INOUT) :: kosdobs ! Observations outside space domain
1155	INTEGER, INTENT(INOUT) :: klanobs ! Observations within a model land cell
1156	INTEGER, INTENT(INOUT) :: knlaobs ! Observations near land
1157	INTEGER, INTENT(INOUT) :: kbdyobs ! Observations near boundary
1158	LOGICAL, INTENT(IN) :: ld_nea ! Flag observations near land
1159	LOGICAL, INTENT(IN) :: ld_bound_reject ! Flag observations near open boundary
1160	INTEGER, INTENT(IN) :: kqc_cutoff ! Cutoff QC value
1161
1162	!! * Local declarations
1163	REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: &
1164	& zgmsk ! Grid mask
1165	REAL(KIND=wp), DIMENSION(2,2,kprofno) :: &
1166	& zbmsk ! Boundary mask
1167	REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask
1168	REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: &
1169	& zgdept, &
1170	& zgdepw
1171	REAL(KIND=wp), DIMENSION(2,2,kprofno) :: &
1172	& zglam, & ! Model longitude at grid points
1173	& zgphi, & ! Model latitude at grid points
1174	& zbathy ! Index of deepest wet level at grid points
1175	INTEGER, DIMENSION(2,2,kprofno) :: &
1176	& igrdi, & ! Grid i,j
1177	& igrdj
1178	LOGICAL :: lgridobs ! Is observation on a model grid point.
1179	LOGICAL :: ll_next_to_land ! Is a profile next to land
1180	INTEGER :: iig, ijg ! i,j of observation on model grid point.
1181	INTEGER :: jobs, jobsp, jk, ji, jj
1182	REAL(KIND=wp) :: maxdept, maxdepw
1183	!!----------------------------------------------------------------------
1184
1185	! Get grid point indices
1186
1187	DO jobs = 1, kprofno
1188
1189	! For invalid points use 2,2
1190
1191	IF ( kpobsqc(jobs) >= kqc_cutoff ) THEN
1192
1193	igrdi(1,1,jobs) = 1
1194	igrdj(1,1,jobs) = 1
1195	igrdi(1,2,jobs) = 1
1196	igrdj(1,2,jobs) = 2
1197	igrdi(2,1,jobs) = 2
1198	igrdj(2,1,jobs) = 1
1199	igrdi(2,2,jobs) = 2
1200	igrdj(2,2,jobs) = 2
1201
1202	ELSE
1203
1204	igrdi(1,1,jobs) = kobsi(jobs)-1
1205	igrdj(1,1,jobs) = kobsj(jobs)-1
1206	igrdi(1,2,jobs) = kobsi(jobs)-1
1207	igrdj(1,2,jobs) = kobsj(jobs)
1208	igrdi(2,1,jobs) = kobsi(jobs)
1209	igrdj(2,1,jobs) = kobsj(jobs)-1
1210	igrdi(2,2,jobs) = kobsi(jobs)
1211	igrdj(2,2,jobs) = kobsj(jobs)
1212
1213	ENDIF
1214
1215	END DO
1216
1217	IF (ln_bdy) THEN
1218	! Create a mask grid points in boundary rim
1219	IF (ld_bound_reject) THEN
1220	zbdymask(:,:) = 1.0_wp
1221	DO ji = 1, nb_bdy
1222	DO jj = 1, idx_bdy(ji)%nblen(1)
1223	zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp
1224	ENDDO
1225	ENDDO
1226	ENDIF
1227
1228	CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, zbdymask, zbmsk )
1229	ENDIF
1230
1231	CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, pmask, zgmsk )
1232	CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, plam, zglam )
1233	CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, pphi, zgphi )
1234	CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, REAL(mbkt), zbathy )
1235	! Need to know the bathy depth for each observation for sco
1236	CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, gdepw_n(:,:,:), &
1237	& zgdepw )
1238	CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, gdept_n(:,:,:), &
1239	& zgdept )
1240
1241	DO jobs = 1, kprofno
1242
1243	! Skip bad profiles
1244	IF ( kpobsqc(jobs) >= kqc_cutoff ) CYCLE
1245
1246	! Check if this observation is on a grid point
1247
1248	lgridobs = .FALSE.
1249	iig = -1
1250	ijg = -1
1251	DO jj = 1, 2
1252	DO ji = 1, 2
1253	IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) &
1254	& .AND. &
1255	& ( ABS( MOD( zglam(ji,jj,jobs) - pobslam(jobs),360.0) ) < 1.0e-6_wp ) &
1256	& ) THEN
1257	lgridobs = .TRUE.
1258	iig = ji
1259	ijg = jj
1260	ENDIF
1261	END DO
1262	END DO
1263
1264	! Check if next to land
1265	IF ( ANY( zgmsk(1:2,1:2,1,jobs) == 0.0_wp ) ) THEN
1266	ll_next_to_land=.TRUE.
1267	ELSE
1268	ll_next_to_land=.FALSE.
1269	ENDIF
1270
1271	! Reject observations
1272
1273	DO jobsp = kpstart(jobs), kpend(jobs)
1274
1275	! Calculate max T and W depths of 2x2 grid
1276	maxdept = zgdept(1,1,NINT(zbathy(1,1,jobs)),jobs)
1277	maxdepw = zgdepw(1,1,NINT(zbathy(1,1,jobs))+1,jobs)
1278	DO jj = 1, 2
1279	DO ji = 1, 2
1280	IF ( zgdept(ji,jj,NINT(zbathy(ji,jj,jobs)),jobs) > maxdept ) THEN
1281	maxdept = zgdept(ji,jj,NINT(zbathy(ji,jj,jobs)),jobs)
1282	END IF
1283	IF ( zgdepw(ji,jj,NINT(zbathy(ji,jj,jobs))+1,jobs) > maxdepw ) THEN
1284	maxdepw = zgdepw(ji,jj,NINT(zbathy(ji,jj,jobs))+1,jobs)
1285	END IF
1286	END DO
1287	END DO
1288
1289	! Flag if the observation falls outside the model spatial domain
1290	IF ( ( pobslam(jobs) < -180. ) &
1291	& .OR. ( pobslam(jobs) > 180. ) &
1292	& .OR. ( pobsphi(jobs) < -90. ) &
1293	& .OR. ( pobsphi(jobs) > 90. ) &
1294	& .OR. ( pobsdep(jobsp) < 0.0 ) &
1295	& .OR. ( pobsdep(jobsp) > gdepw_1d(kpk) ) ) THEN
1296	kobsqc(jobsp) = IBSET(kobsqc(jobsp),11)
1297	kosdobs = kosdobs + 1
1298	CYCLE
1299	ENDIF
1300
1301	! To check if an observations falls within land:
1302
1303	! Flag if the observation is deeper than the bathymetry
1304	! Or if it is within the mask
1305	IF ( ALL( zgdepw(1:2,1:2,kpk,jobs) < pobsdep(jobsp) ) &
1306	& .OR. &
1307	& ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) &
1308	& == 0.0_wp) ) THEN
1309	kobsqc(jobsp) = IBSET(kobsqc(jobsp),10)
1310	klanobs = klanobs + 1
1311	CYCLE
1312	ENDIF
1313
1314	! Flag if the observation is close to land
1315	IF ( ll_next_to_land ) THEN
1316	knlaobs = knlaobs + 1
1317	IF (ld_nea) THEN
1318	kobsqc(jobsp) = IBSET(kobsqc(jobsp),10)
1319	ENDIF
1320	ENDIF
1321
1322	! For observations on the grid reject them if their are at
1323	! a masked point
1324
1325	IF (lgridobs) THEN
1326	IF (zgmsk(iig,ijg,kobsk(jobsp)-1,jobs) == 0.0_wp ) THEN
1327	kobsqc(jobsp) = IBSET(kobsqc(jobs),10)
1328	klanobs = klanobs + 1
1329	CYCLE
1330	ENDIF
1331	ENDIF
1332
1333	! Flag if the observation falls is close to land
1334	IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == &
1335	& 0.0_wp) THEN
1336	IF (ld_nea) kobsqc(jobsp) = kobsqc(jobsp) + 14
1337	knlaobs = knlaobs + 1
1338	ENDIF
1339
1340	! Set observation depth equal to that of the first model depth
1341	IF ( pobsdep(jobsp) <= pdep(1) ) THEN
1342	pobsdep(jobsp) = pdep(1)
1343	ENDIF
1344	!IF ( pobsdep(jobsp) < MINVAL(zgdept(1:2,1:2,1,jobs) ) ) THEN
1345	! pobsdep(jobsp) = MINVAL(zgdept(1:2,1:2,1,jobs))
1346	!ENDIF
1347
1348	! Set observation depth equal to that of the last wet T-point
1349	!IF ( ( pobsdep(jobsp) > maxdept ) .AND. &
1350	! & ( pobsdep(jobsp) < maxdepw ) ) THEN
1351	! pobsdep(jobsp) = maxdept
1352	!END IF
1353
1354	IF (ln_bdy) THEN
1355	! Flag if the observation falls close to the boundary rim
1356	IF (ld_bound_reject) THEN
1357	IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN
1358	kobsqc(jobsp) = IBSET(kobsqc(jobs),8)
1359	kbdyobs = kbdyobs + 1
1360	CYCLE
1361	ENDIF
1362	! for observations on the grid...
1363	IF (lgridobs) THEN
1364	IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN
1365	kobsqc(jobsp) = IBSET(kobsqc(jobs),8)
1366	kbdyobs = kbdyobs + 1
1367	CYCLE
1368	ENDIF
1369	ENDIF
1370	ENDIF
1371	ENDIF
1372	!
1373	END DO
1374	END DO
1375	!
1376	END SUBROUTINE obs_coo_spc_3d
1377
1378
1379	SUBROUTINE obs_pro_rej( profdata, kqc_cutoff )
1380	!!----------------------------------------------------------------------
1381	!! * ROUTINE obs_pro_rej *
1382	!!
1383	!! ** Purpose : Reject all data within a rejected profile
1384	!!
1385	!! ** Method :
1386	!!
1387	!! ** Action :
1388	!!
1389	!! References :
1390	!!
1391	!! History : 2007-10 (K. Mogensen) Original code
1392	!!----------------------------------------------------------------------
1393	TYPE(obs_prof), INTENT(inout) :: profdata ! Profile data
1394	INTEGER , INTENT(in ) :: kqc_cutoff ! QC cutoff value
1395	!
1396	INTEGER :: jprof
1397	INTEGER :: jvar
1398	INTEGER :: jobs
1399	!!----------------------------------------------------------------------
1400
1401	! Loop over profiles
1402
1403	DO jprof = 1, profdata%nprof
1404
1405	IF ( profdata%nqc(jprof) > kqc_cutoff ) THEN
1406
1407	DO jvar = 1, profdata%nvar
1408
1409	DO jobs = profdata%npvsta(jprof,jvar), &
1410	& profdata%npvend(jprof,jvar)
1411
1412	profdata%var(jvar)%nvqc(jobs) = &
1413	& IBSET(profdata%var(jvar)%nvqc(jobs),14)
1414
1415	END DO
1416
1417	END DO
1418
1419	ENDIF
1420
1421	END DO
1422	!
1423	END SUBROUTINE obs_pro_rej
1424
1425
1426	SUBROUTINE obs_uv_rej( profdata, knumu, knumv, kqc_cutoff, kuvar, kvvar )
1427	!!----------------------------------------------------------------------
1428	!! * ROUTINE obs_uv_rej *
1429	!!
1430	!! ** Purpose : Reject u if v is rejected and vice versa
1431	!!
1432	!! ** Method :
1433	!!
1434	!! ** Action :
1435	!!
1436	!! References :
1437	!!
1438	!! History : 2009-2 (K. Mogensen) Original code
1439	!!----------------------------------------------------------------------
1440	TYPE(obs_prof), INTENT(INOUT) :: profdata ! Profile data
1441	INTEGER, INTENT(INOUT) :: knumu ! Number of u rejected
1442	INTEGER, INTENT(INOUT) :: knumv ! Number of v rejected
1443	INTEGER, INTENT(IN) :: kqc_cutoff ! QC cutoff value
1444	INTEGER, INTENT(IN) :: kuvar ! Index of u
1445	INTEGER, INTENT(IN) :: kvvar ! Index of v
1446	!
1447	INTEGER :: jprof
1448	INTEGER :: jvar
1449	INTEGER :: jobs
1450	!!----------------------------------------------------------------------
1451
1452	DO jprof = 1, profdata%nprof !== Loop over profiles ==!
1453	!
1454	IF ( ( profdata%npvsta(jprof,kuvar) /= profdata%npvsta(jprof,kvvar) ) .OR. &
1455	& ( profdata%npvend(jprof,kuvar) /= profdata%npvend(jprof,kvvar) ) ) THEN
1456	!
1457	CALL ctl_stop('U,V profiles inconsistent in obs_uv_rej')
1458	RETURN
1459	!
1460	ENDIF
1461	!
1462	DO jobs = profdata%npvsta(jprof,kuvar), profdata%npvend(jprof,kuvar)
1463	!
1464	IF ( ( profdata%var(kuvar)%nvqc(jobs) > kqc_cutoff ) .AND. &
1465	& ( profdata%var(kvvar)%nvqc(jobs) <= kqc_cutoff) ) THEN
1466	profdata%var(kvvar)%nvqc(jobs) = IBSET(profdata%var(kuvar)%nvqc(jobs),15)
1467	knumv = knumv + 1
1468	ENDIF
1469	IF ( ( profdata%var(kvvar)%nvqc(jobs) > kqc_cutoff ) .AND. &
1470	& ( profdata%var(kuvar)%nvqc(jobs) <= kqc_cutoff) ) THEN
1471	profdata%var(kuvar)%nvqc(jobs) = IBSET(profdata%var(kuvar)%nvqc(jobs),15)
1472	knumu = knumu + 1
1473	ENDIF
1474	!
1475	END DO
1476	!
1477	END DO
1478	!
1479	END SUBROUTINE obs_uv_rej
1480
1481	!!=====================================================================
1482	END MODULE obs_prep

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

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

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.4_FOAM_package/src/OCE/OBS/obs_prep.F90 @ 15799

Download in other formats: