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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

obs_prep.F90 in branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

Context Navigation

source: branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 9205

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

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

Download in other formats: