Context Navigation

obs_prep.F90 @ 7932

Last change on this file since 7932 was 7932, checked in by jwhile, 8 years ago
Bug fixes
File size: 60.8 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_var1, ld_var2, &
258	& kpi, kpj, kpk, &
259	& zmask1, pglam1, pgphi1, zmask2, pglam2, pgphi2, &
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, INTENT(IN) :: ld_var1 ! Observed variables switches
287	LOGICAL, INTENT(IN) :: ld_var2
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) :: &
294	& zmask1, &
295	& zmask2
296	REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
297	& pglam1, &
298	& pglam2, &
299	& pgphi1, &
300	& pgphi2
301	INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff ! cut off for QC value
302
303	!! * Local declarations
304	INTEGER :: iqc_cutoff = 255 ! cut off for QC value
305	INTEGER :: iyea0 ! Initial date
306	INTEGER :: imon0 ! - (year, month, day, hour, minute)
307	INTEGER :: iday0
308	INTEGER :: ihou0
309	INTEGER :: imin0
310	INTEGER :: icycle ! Current assimilation cycle
311	! Counters for observations that are
312	INTEGER :: iotdobs ! - outside time domain
313	INTEGER :: iosdv1obs ! - outside space domain (variable 1)
314	INTEGER :: iosdv2obs ! - outside space domain (variable 2)
315	INTEGER :: ilanv1obs ! - within a model land cell (variable 1)
316	INTEGER :: ilanv2obs ! - within a model land cell (variable 2)
317	INTEGER :: inlav1obs ! - close to land (variable 1)
318	INTEGER :: inlav2obs ! - close to land (variable 2)
319	INTEGER :: ibdyv1obs ! - boundary (variable 1)
320	INTEGER :: ibdyv2obs ! - boundary (variable 2)
321	INTEGER :: igrdobs ! - fail the grid search
322	INTEGER :: iuvchku ! - reject u if v rejected and vice versa
323	INTEGER :: iuvchkv !
324	! Global counters for observations that are
325	INTEGER :: iotdobsmpp ! - outside time domain
326	INTEGER :: iosdv1obsmpp ! - outside space domain (variable 1)
327	INTEGER :: iosdv2obsmpp ! - outside space domain (variable 2)
328	INTEGER :: ilanv1obsmpp ! - within a model land cell (variable 1)
329	INTEGER :: ilanv2obsmpp ! - within a model land cell (variable 2)
330	INTEGER :: inlav1obsmpp ! - close to land (variable 1)
331	INTEGER :: inlav2obsmpp ! - close to land (variable 2)
332	INTEGER :: ibdyv1obsmpp ! - boundary (variable 1)
333	INTEGER :: ibdyv2obsmpp ! - boundary (variable 2)
334	INTEGER :: igrdobsmpp ! - fail the grid search
335	INTEGER :: iuvchkumpp ! - reject var1 if var2 rejected and vice versa
336	INTEGER :: iuvchkvmpp !
337	TYPE(obs_prof_valid) :: llvalid ! Profile selection
338	TYPE(obs_prof_valid), DIMENSION(profdata%nvar) :: &
339	& llvvalid ! var1,var2 selection
340	INTEGER :: jvar ! Variable loop variable
341	INTEGER :: jobs ! Obs. loop variable
342	INTEGER :: jstp ! Time loop variable
343	INTEGER :: inrc ! Time index variable
344
345	IF(lwp) WRITE(numout,*)'obs_pre_prof: Preparing the profile data...'
346	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
347
348	! Initial date initialization (year, month, day, hour, minute)
349	iyea0 = ndate0 / 10000
350	imon0 = ( ndate0 - iyea0 * 10000 ) / 100
351	iday0 = ndate0 - iyea0 * 10000 - imon0 * 100
352	ihou0 = 0
353	imin0 = 0
354
355	icycle = no ! Assimilation cycle
356
357	! Diagnotics counters for various failures.
358
359	iotdobs = 0
360	igrdobs = 0
361	iosdv1obs = 0
362	iosdv2obs = 0
363	ilanv1obs = 0
364	ilanv2obs = 0
365	inlav1obs = 0
366	inlav2obs = 0
367	ibdyv1obs = 0
368	ibdyv2obs = 0
369	iuvchku = 0
370	iuvchkv = 0
371
372
373	! Set QC cutoff to optional value if provided
374	IF ( PRESENT(kqc_cutoff) ) iqc_cutoff=kqc_cutoff
375
376	! -----------------------------------------------------------------------
377	! Find time coordinate for profiles
378	! -----------------------------------------------------------------------
379
380	IF ( PRESENT(kdailyavtypes) ) THEN
381	CALL obs_coo_tim_prof( icycle, &
382	& iyea0, imon0, iday0, ihou0, imin0, &
383	& profdata%nprof, profdata%nyea, profdata%nmon, &
384	& profdata%nday, profdata%nhou, profdata%nmin, &
385	& profdata%ntyp, profdata%nqc, profdata%mstp, &
386	& iotdobs, kdailyavtypes = kdailyavtypes, &
387	& kqc_cutoff = iqc_cutoff )
388	ELSE
389	CALL obs_coo_tim_prof( icycle, &
390	& iyea0, imon0, iday0, ihou0, imin0, &
391	& profdata%nprof, profdata%nyea, profdata%nmon, &
392	& profdata%nday, profdata%nhou, profdata%nmin, &
393	& profdata%ntyp, profdata%nqc, profdata%mstp, &
394	& iotdobs, kqc_cutoff = iqc_cutoff )
395	ENDIF
396
397	CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
398
399	! -----------------------------------------------------------------------
400	! Check for profiles failing the grid search
401	! -----------------------------------------------------------------------
402
403	CALL obs_coo_grd( profdata%nprof, profdata%mi(:,1), profdata%mj(:,1), &
404	& profdata%nqc, igrdobs )
405	CALL obs_coo_grd( profdata%nprof, profdata%mi(:,2), profdata%mj(:,2), &
406	& profdata%nqc, igrdobs )
407
408	CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
409
410	! -----------------------------------------------------------------------
411	! Reject all observations for profiles with nqc > iqc_cutoff
412	! -----------------------------------------------------------------------
413
414	CALL obs_pro_rej( profdata, kqc_cutoff = iqc_cutoff )
415
416	! -----------------------------------------------------------------------
417	! Check for land points. This includes points below the model
418	! bathymetry so this is done for every point in the profile
419	! -----------------------------------------------------------------------
420
421	! Variable 1
422	CALL obs_coo_spc_3d( profdata%nprof, profdata%nvprot(1), &
423	& profdata%npvsta(:,1), profdata%npvend(:,1), &
424	& jpi, jpj, &
425	& jpk, &
426	& profdata%mi, profdata%mj, &
427	& profdata%var(1)%mvk, &
428	& profdata%rlam, profdata%rphi, &
429	& profdata%var(1)%vdep, &
430	& pglam1, pgphi1, &
431	& gdept_1d, zmask1, &
432	& profdata%nqc, profdata%var(1)%nvqc, &
433	& iosdv1obs, ilanv1obs, &
434	& inlav1obs, ld_nea, &
435	& ibdyv1obs, ld_bound_reject, &
436	& iqc_cutoff )
437
438	CALL obs_mpp_sum_integer( iosdv1obs, iosdv1obsmpp )
439	CALL obs_mpp_sum_integer( ilanv1obs, ilanv1obsmpp )
440	CALL obs_mpp_sum_integer( inlav1obs, inlav1obsmpp )
441	CALL obs_mpp_sum_integer( ibdyv1obs, ibdyv1obsmpp )
442
443	! Variable 2
444	CALL obs_coo_spc_3d( profdata%nprof, profdata%nvprot(2), &
445	& profdata%npvsta(:,2), profdata%npvend(:,2), &
446	& jpi, jpj, &
447	& jpk, &
448	& profdata%mi, profdata%mj, &
449	& profdata%var(2)%mvk, &
450	& profdata%rlam, profdata%rphi, &
451	& profdata%var(2)%vdep, &
452	& pglam2, pgphi2, &
453	& gdept_1d, zmask2, &
454	& profdata%nqc, profdata%var(2)%nvqc, &
455	& iosdv2obs, ilanv2obs, &
456	& inlav2obs, ld_nea, &
457	& ibdyv2obs, ld_bound_reject, &
458	& iqc_cutoff )
459
460	CALL obs_mpp_sum_integer( iosdv2obs, iosdv2obsmpp )
461	CALL obs_mpp_sum_integer( ilanv2obs, ilanv2obsmpp )
462	CALL obs_mpp_sum_integer( inlav2obs, inlav2obsmpp )
463	CALL obs_mpp_sum_integer( ibdyv2obs, ibdyv2obsmpp )
464
465	! -----------------------------------------------------------------------
466	! Reject u if v is rejected and vice versa
467	! -----------------------------------------------------------------------
468
469	IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN
470	CALL obs_uv_rej( profdata, iuvchku, iuvchkv, iqc_cutoff )
471	CALL obs_mpp_sum_integer( iuvchku, iuvchkumpp )
472	CALL obs_mpp_sum_integer( iuvchkv, iuvchkvmpp )
473	ENDIF
474
475	! -----------------------------------------------------------------------
476	! Copy useful data from the profdata data structure to
477	! the prodatqc data structure
478	! -----------------------------------------------------------------------
479
480	! Allocate the selection arrays
481
482	ALLOCATE( llvalid%luse(profdata%nprof) )
483	DO jvar = 1,profdata%nvar
484	ALLOCATE( llvvalid(jvar)%luse(profdata%nvprot(jvar)) )
485	END DO
486
487	! We want all data which has qc flags <= iqc_cutoff
488
489	llvalid%luse(:) = ( profdata%nqc(:) <= iqc_cutoff )
490	DO jvar = 1,profdata%nvar
491	llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= iqc_cutoff )
492	END DO
493
494	! The actual copying
495
496	CALL obs_prof_compress( profdata, prodatqc, .TRUE., numout, &
497	& lvalid=llvalid, lvvalid=llvvalid )
498
499	! Dellocate the selection arrays
500	DEALLOCATE( llvalid%luse )
501	DO jvar = 1,profdata%nvar
502	DEALLOCATE( llvvalid(jvar)%luse )
503	END DO
504
505	! -----------------------------------------------------------------------
506	! Print information about what observations are left after qc
507	! -----------------------------------------------------------------------
508
509	! Update the total observation counter array
510
511	IF(lwp) THEN
512
513	WRITE(numout,*)
514	WRITE(numout,*) ' Profiles outside time domain = ', &
515	& iotdobsmpp
516	WRITE(numout,*) ' Remaining profiles that failed grid search = ', &
517	& igrdobsmpp
518	WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data outside space domain = ', &
519	& iosdv1obsmpp
520	WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data at land points = ', &
521	& ilanv1obsmpp
522	IF (ld_nea) THEN
523	WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data near land points (removed) = ',&
524	& inlav1obsmpp
525	ELSE
526	WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data near land points (kept) = ',&
527	& inlav1obsmpp
528	ENDIF
529	IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN
530	WRITE(numout,*) ' U observation rejected since V rejected = ', &
531	& iuvchku
532	ENDIF
533	WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data near open boundary (removed) = ',&
534	& ibdyv1obsmpp
535	WRITE(numout,*) ' '//prodatqc%cvars(1)//' data accepted = ', &
536	& prodatqc%nvprotmpp(1)
537	WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data outside space domain = ', &
538	& iosdv2obsmpp
539	WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data at land points = ', &
540	& ilanv2obsmpp
541	IF (ld_nea) THEN
542	WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near land points (removed) = ',&
543	& inlav2obsmpp
544	ELSE
545	WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near land points (kept) = ',&
546	& inlav2obsmpp
547	ENDIF
548	IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN
549	WRITE(numout,*) ' V observation rejected since U rejected = ', &
550	& iuvchkv
551	ENDIF
552	WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near open boundary (removed) = ',&
553	& ibdyv2obsmpp
554	WRITE(numout,*) ' '//prodatqc%cvars(2)//' data accepted = ', &
555	& prodatqc%nvprotmpp(2)
556
557	WRITE(numout,*)
558	WRITE(numout,*) ' Number of observations per time step :'
559	WRITE(numout,*)
560	WRITE(numout,'(10X,A,5X,A,5X,A,A)')'Time step','Profiles', &
561	& ' '//prodatqc%cvars(1)//' ', &
562	& ' '//prodatqc%cvars(2)//' '
563	WRITE(numout,998)
564	ENDIF
565
566	DO jobs = 1, prodatqc%nprof
567	inrc = prodatqc%mstp(jobs) + 2 - nit000
568	prodatqc%npstp(inrc) = prodatqc%npstp(inrc) + 1
569	DO jvar = 1, prodatqc%nvar
570	IF ( prodatqc%npvend(jobs,jvar) > 0 ) THEN
571	prodatqc%nvstp(inrc,jvar) = prodatqc%nvstp(inrc,jvar) + &
572	& ( prodatqc%npvend(jobs,jvar) - &
573	& prodatqc%npvsta(jobs,jvar) + 1 )
574	ENDIF
575	END DO
576	END DO
577
578
579	CALL obs_mpp_sum_integers( prodatqc%npstp, prodatqc%npstpmpp, &
580	& nitend - nit000 + 2 )
581	DO jvar = 1, prodatqc%nvar
582	CALL obs_mpp_sum_integers( prodatqc%nvstp(:,jvar), &
583	& prodatqc%nvstpmpp(:,jvar), &
584	& nitend - nit000 + 2 )
585	END DO
586
587	IF ( lwp ) THEN
588	DO jstp = nit000 - 1, nitend
589	inrc = jstp - nit000 + 2
590	WRITE(numout,999) jstp, prodatqc%npstpmpp(inrc), &
591	& prodatqc%nvstpmpp(inrc,1), &
592	& prodatqc%nvstpmpp(inrc,2)
593	END DO
594	ENDIF
595
596	998 FORMAT(10X,'---------',5X,'--------',5X,'-----------',5X,'----------------')
597	999 FORMAT(10X,I9,5X,I8,5X,I11,5X,I8)
598
599	END SUBROUTINE obs_pre_prof
600
601	SUBROUTINE obs_coo_tim( kcycle, &
602	& kyea0, kmon0, kday0, khou0, kmin0, &
603	& kobsno, &
604	& kobsyea, kobsmon, kobsday, kobshou, kobsmin, &
605	& kobsqc, kobsstp, kotdobs )
606	!!----------------------------------------------------------------------
607	!! * ROUTINE obs_coo_tim *
608	!!
609	!! ** Purpose : Compute the number of time steps to the observation time.
610	!!
611	!! ** Method : For time coordinates ( yea_obs, mon_obs, day_obs,
612	!! hou_obs, min_obs ), this routine locates the time step
613	!! that is closest to this time.
614	!!
615	!! ** Action :
616	!!
617	!! References :
618	!!
619	!! History :
620	!! ! 1997-07 (A. Weaver) Original
621	!! ! 2006-08 (A. Weaver) NEMOVAR migration
622	!! ! 2006-10 (A. Weaver) Cleanup
623	!! ! 2007-01 (K. Mogensen) Rewritten with loop
624	!! ! 2010-05 (D. Lea) Fix in leap year calculation for NEMO vn3.2
625	!!----------------------------------------------------------------------
626	!! * Modules used
627	USE dom_oce, ONLY : & ! Geographical information
628	& rdt
629	USE phycst, ONLY : & ! Physical constants
630	& rday, &
631	& rmmss, &
632	& rhhmm
633	!! * Arguments
634	INTEGER, INTENT(IN) :: kcycle ! Current cycle
635	INTEGER, INTENT(IN) :: kyea0 ! Initial date coordinates
636	INTEGER, INTENT(IN) :: kmon0
637	INTEGER, INTENT(IN) :: kday0
638	INTEGER, INTENT(IN) :: khou0
639	INTEGER, INTENT(IN) :: kmin0
640	INTEGER, INTENT(IN) :: kobsno ! Number of observations
641	INTEGER, INTENT(INOUT) :: kotdobs ! Number of observations failing time check
642	INTEGER, DIMENSION(kobsno), INTENT(IN ) :: &
643	& kobsyea, & ! Observation time coordinates
644	& kobsmon, &
645	& kobsday, &
646	& kobshou, &
647	& kobsmin
648	INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
649	& kobsqc ! Quality control flag
650	INTEGER, DIMENSION(kobsno), INTENT(OUT) :: &
651	& kobsstp ! Number of time steps up to the
652	! observation time
653
654	!! * Local declarations
655	INTEGER :: jyea
656	INTEGER :: jmon
657	INTEGER :: jday
658	INTEGER :: jobs
659	INTEGER :: iyeastr
660	INTEGER :: iyeaend
661	INTEGER :: imonstr
662	INTEGER :: imonend
663	INTEGER :: idaystr
664	INTEGER :: idayend
665	INTEGER :: iskip
666	INTEGER :: idaystp
667	REAL(KIND=wp) :: zminstp
668	REAL(KIND=wp) :: zhoustp
669	REAL(KIND=wp) :: zobsstp
670	INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year
671
672	!-----------------------------------------------------------------------
673	! Initialization
674	!-----------------------------------------------------------------------
675
676	! Intialize the number of time steps per day
677	idaystp = NINT( rday / rdt )
678
679	!---------------------------------------------------------------------
680	! Locate the model time coordinates for interpolation
681	!---------------------------------------------------------------------
682
683	DO jobs = 1, kobsno
684
685	! Initialize the time step counter
686	kobsstp(jobs) = nit000 - 1
687
688	! Flag if observation date is less than the initial date
689
690	IF ( ( kobsyea(jobs) < kyea0 ) &
691	& .OR. ( ( kobsyea(jobs) == kyea0 ) &
692	& .AND. ( kobsmon(jobs) < kmon0 ) ) &
693	& .OR. ( ( kobsyea(jobs) == kyea0 ) &
694	& .AND. ( kobsmon(jobs) == kmon0 ) &
695	& .AND. ( kobsday(jobs) < kday0 ) ) &
696	& .OR. ( ( kobsyea(jobs) == kyea0 ) &
697	& .AND. ( kobsmon(jobs) == kmon0 ) &
698	& .AND. ( kobsday(jobs) == kday0 ) &
699	& .AND. ( kobshou(jobs) < khou0 ) ) &
700	& .OR. ( ( kobsyea(jobs) == kyea0 ) &
701	& .AND. ( kobsmon(jobs) == kmon0 ) &
702	& .AND. ( kobsday(jobs) == kday0 ) &
703	& .AND. ( kobshou(jobs) == khou0 ) &
704	& .AND. ( kobsmin(jobs) <= kmin0 ) ) ) THEN
705	kobsstp(jobs) = -1
706	kobsqc(jobs) = kobsqc(jobs) + 11
707	kotdobs = kotdobs + 1
708	CYCLE
709	ENDIF
710
711	! Compute the number of time steps to the observation day
712	iyeastr = kyea0
713	iyeaend = kobsyea(jobs)
714
715	!---------------------------------------------------------------------
716	! Year loop
717	!---------------------------------------------------------------------
718	DO jyea = iyeastr, iyeaend
719
720	CALL calc_month_len( jyea, imonth_len )
721
722	imonstr = 1
723	IF ( jyea == kyea0 ) imonstr = kmon0
724	imonend = 12
725	IF ( jyea == kobsyea(jobs) ) imonend = kobsmon(jobs)
726
727	! Month loop
728	DO jmon = imonstr, imonend
729
730	idaystr = 1
731	IF ( ( jmon == kmon0 ) &
732	& .AND. ( jyea == kyea0 ) ) idaystr = kday0
733	idayend = imonth_len(jmon)
734	IF ( ( jmon == kobsmon(jobs) ) &
735	& .AND. ( jyea == kobsyea(jobs) ) ) idayend = kobsday(jobs) - 1
736
737	! Day loop
738	DO jday = idaystr, idayend
739	kobsstp(jobs) = kobsstp(jobs) + idaystp
740	END DO
741
742	END DO
743
744	END DO
745
746	! Add in the number of time steps to the observation minute
747	zminstp = rmmss / rdt
748	zhoustp = rhhmm * zminstp
749
750	zobsstp = REAL( kobsmin(jobs) - kmin0, KIND=wp ) * zminstp &
751	& + REAL( kobshou(jobs) - khou0, KIND=wp ) * zhoustp
752	kobsstp(jobs) = kobsstp(jobs) + NINT( zobsstp )
753
754	! Flag if observation step outside the time window
755	IF ( ( kobsstp(jobs) < ( nit000 - 1 ) ) &
756	& .OR.( kobsstp(jobs) > nitend ) ) THEN
757	kobsqc(jobs) = kobsqc(jobs) + 12
758	kotdobs = kotdobs + 1
759	CYCLE
760	ENDIF
761
762	END DO
763
764	END SUBROUTINE obs_coo_tim
765
766	SUBROUTINE calc_month_len( iyear, imonth_len )
767	!!----------------------------------------------------------------------
768	!! * ROUTINE calc_month_len *
769	!!
770	!! ** Purpose : Compute the number of days in a months given a year.
771	!!
772	!! ** Method :
773	!!
774	!! ** Action :
775	!!
776	!! History :
777	!! ! 10-05 (D. Lea) New routine based on day_init
778	!!----------------------------------------------------------------------
779
780	INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year
781	INTEGER :: iyear !: year
782
783	! length of the month of the current year (from nleapy, read in namelist)
784	IF ( nleapy < 2 ) THEN
785	imonth_len(:) = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)
786	IF ( nleapy == 1 ) THEN ! we are using calendar with leap years
787	IF ( MOD(iyear, 4) == 0 .AND. ( MOD(iyear, 400) == 0 .OR. MOD(iyear, 100) /= 0 ) ) THEN
788	imonth_len(2) = 29
789	ENDIF
790	ENDIF
791	ELSE
792	imonth_len(:) = nleapy ! all months with nleapy days per year
793	ENDIF
794
795	END SUBROUTINE
796
797	SUBROUTINE obs_coo_tim_prof( kcycle, &
798	& kyea0, kmon0, kday0, khou0, kmin0, &
799	& kobsno, &
800	& kobsyea, kobsmon, kobsday, kobshou, kobsmin, &
801	& ktyp, kobsqc, kobsstp, kotdobs, kdailyavtypes, &
802	& kqc_cutoff )
803	!!----------------------------------------------------------------------
804	!! * ROUTINE obs_coo_tim *
805	!!
806	!! ** Purpose : Compute the number of time steps to the observation time.
807	!!
808	!! ** Method : For time coordinates ( yea_obs, mon_obs, day_obs,
809	!! hou_obs, min_obs ), this routine locates the time step
810	!! that is closest to this time.
811	!!
812	!! ** Action :
813	!!
814	!! References :
815	!!
816	!! History :
817	!! ! 1997-07 (A. Weaver) Original
818	!! ! 2006-08 (A. Weaver) NEMOVAR migration
819	!! ! 2006-10 (A. Weaver) Cleanup
820	!! ! 2007-01 (K. Mogensen) Rewritten with loop
821	!!----------------------------------------------------------------------
822	!! * Modules used
823	!! * Arguments
824	INTEGER, INTENT(IN) :: kcycle ! Current cycle
825	INTEGER, INTENT(IN) :: kyea0 ! Initial date coordinates
826	INTEGER, INTENT(IN) :: kmon0
827	INTEGER, INTENT(IN) :: kday0
828	INTEGER, INTENT(IN) :: khou0
829	INTEGER, INTENT(IN) :: kmin0
830	INTEGER, INTENT(IN) :: kobsno ! Number of observations
831	INTEGER, INTENT(INOUT) :: kotdobs ! Number of observations failing time check
832	INTEGER, DIMENSION(kobsno), INTENT(IN ) :: &
833	& kobsyea, & ! Observation time coordinates
834	& kobsmon, &
835	& kobsday, &
836	& kobshou, &
837	& kobsmin, &
838	& ktyp ! Observation type.
839	INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
840	& kobsqc ! Quality control flag
841	INTEGER, DIMENSION(kobsno), INTENT(OUT) :: &
842	& kobsstp ! Number of time steps up to the
843	! observation time
844	INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
845	& kdailyavtypes ! Types for daily averages
846	INTEGER, OPTIONAL, INTENT(IN) :: kqc_cutoff ! QC cutoff value
847
848	!! * Local declarations
849	INTEGER :: jobs
850	INTEGER :: iqc_cutoff=255
851
852	!-----------------------------------------------------------------------
853	! Call standard obs_coo_tim
854	!-----------------------------------------------------------------------
855
856	CALL obs_coo_tim( kcycle, &
857	& kyea0, kmon0, kday0, khou0, kmin0, &
858	& kobsno, &
859	& kobsyea, kobsmon, kobsday, kobshou, kobsmin, &
860	& kobsqc, kobsstp, kotdobs )
861
862	!------------------------------------------------------------------------
863	! Always reject daily averaged data (e.g. MRB data (820)) at initial time
864	!------------------------------------------------------------------------
865
866	IF ( PRESENT(kdailyavtypes) ) THEN
867	DO jobs = 1, kobsno
868
869	IF ( kobsqc(jobs) <= iqc_cutoff ) THEN
870
871	IF ( ( kobsstp(jobs) == (nit000 - 1) ).AND.&
872	& ( ANY (kdailyavtypes(:) == ktyp(jobs)) ) ) THEN
873	kobsqc(jobs) = IBSET(kobsqc(jobs),13)
874	kotdobs = kotdobs + 1
875	CYCLE
876	ENDIF
877
878	ENDIF
879	END DO
880	ENDIF
881
882
883	END SUBROUTINE obs_coo_tim_prof
884
885	SUBROUTINE obs_coo_grd( kobsno, kobsi, kobsj, kobsqc, kgrdobs )
886	!!----------------------------------------------------------------------
887	!! * ROUTINE obs_coo_grd *
888	!!
889	!! ** Purpose : Verify that the grid search has not failed
890	!!
891	!! ** Method : The previously computed i,j indeces are checked
892	!!
893	!! ** Action :
894	!!
895	!! References :
896	!!
897	!! History :
898	!! ! 2007-01 (K. Mogensen) Original
899	!!----------------------------------------------------------------------
900	!! * Arguments
901	INTEGER, INTENT(IN) :: kobsno ! Number of observations
902	INTEGER, DIMENSION(kobsno), INTENT(IN ) :: &
903	& kobsi, & ! i,j indeces previously computed
904	& kobsj
905	INTEGER, INTENT(INOUT) :: kgrdobs ! Number of observations failing the check
906	INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
907	& kobsqc ! Quality control flag
908
909	!! * Local declarations
910	INTEGER :: jobs ! Loop variable
911
912	! Flag if the grid search failed
913
914	DO jobs = 1, kobsno
915	IF ( ( kobsi(jobs) <= 0 ) .AND. ( kobsj(jobs) <= 0 ) ) THEN
916	kobsqc(jobs) = IBSET(kobsqc(jobs),12)
917	kgrdobs = kgrdobs + 1
918	ENDIF
919	END DO
920
921	END SUBROUTINE obs_coo_grd
922
923	SUBROUTINE obs_coo_spc_2d( kobsno, kpi, kpj, &
924	& kobsi, kobsj, pobslam, pobsphi, &
925	& plam, pphi, pmask, &
926	& kobsqc, kosdobs, klanobs, &
927	& knlaobs,ld_nea, &
928	& kbdyobs,ld_bound_reject, &
929	& kqc_cutoff )
930	!!----------------------------------------------------------------------
931	!! * ROUTINE obs_coo_spc_2d *
932	!!
933	!! ** Purpose : Check for points outside the domain and land points
934	!!
935	!! ** Method : Remove the observations that are outside the model space
936	!! and time domain or located within model land cells.
937	!!
938	!! ** Action :
939	!!
940	!! History :
941	!! ! 2007-03 (A. Weaver, K. Mogensen) Original
942	!! ! 2007-06 (K. Mogensen et al) Reject obs. near land.
943	!!----------------------------------------------------------------------
944	!! * Modules used
945
946	!! * Arguments
947	INTEGER, INTENT(IN) :: kobsno ! Total number of observations
948	INTEGER, INTENT(IN) :: kpi ! Number of grid points in (i,j)
949	INTEGER, INTENT(IN) :: kpj
950	INTEGER, DIMENSION(kobsno), INTENT(IN) :: &
951	& kobsi, & ! Observation (i,j) coordinates
952	& kobsj
953	REAL(KIND=wp), DIMENSION(kobsno), INTENT(IN) :: &
954	& pobslam, & ! Observation (lon,lat) coordinates
955	& pobsphi
956	REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) :: &
957	& plam, pphi ! Model (lon,lat) coordinates
958	REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) :: &
959	& pmask ! Land mask array
960	INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
961	& kobsqc ! Observation quality control
962	INTEGER, INTENT(INOUT) :: kosdobs ! Observations outside space domain
963	INTEGER, INTENT(INOUT) :: klanobs ! Observations within a model land cell
964	INTEGER, INTENT(INOUT) :: knlaobs ! Observations near land
965	INTEGER, INTENT(INOUT) :: kbdyobs ! Observations near boundary
966	LOGICAL, INTENT(IN) :: ld_nea ! Flag observations near land
967	LOGICAL, INTENT(IN) :: ld_bound_reject ! Flag observations near open boundary
968	INTEGER, INTENT(IN) :: kqc_cutoff ! Cutoff QC value
969
970	!! * Local declarations
971	REAL(KIND=wp), DIMENSION(2,2,kobsno) :: &
972	& zgmsk ! Grid mask
973	#if defined key_bdy
974	REAL(KIND=wp), DIMENSION(2,2,kobsno) :: &
975	& zbmsk ! Boundary mask
976	REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask
977	#endif
978	REAL(KIND=wp), DIMENSION(2,2,kobsno) :: &
979	& zglam, & ! Model longitude at grid points
980	& zgphi ! Model latitude at grid points
981	INTEGER, DIMENSION(2,2,kobsno) :: &
982	& igrdi, & ! Grid i,j
983	& igrdj
984	LOGICAL :: lgridobs ! Is observation on a model grid point.
985	INTEGER :: iig, ijg ! i,j of observation on model grid point.
986	INTEGER :: jobs, ji, jj
987
988	! Get grid point indices
989
990	DO jobs = 1, kobsno
991
992	! For invalid points use 2,2
993
994	IF ( kobsqc(jobs) >= kqc_cutoff ) THEN
995
996	igrdi(1,1,jobs) = 1
997	igrdj(1,1,jobs) = 1
998	igrdi(1,2,jobs) = 1
999	igrdj(1,2,jobs) = 2
1000	igrdi(2,1,jobs) = 2
1001	igrdj(2,1,jobs) = 1
1002	igrdi(2,2,jobs) = 2
1003	igrdj(2,2,jobs) = 2
1004
1005	ELSE
1006
1007	igrdi(1,1,jobs) = kobsi(jobs)-1
1008	igrdj(1,1,jobs) = kobsj(jobs)-1
1009	igrdi(1,2,jobs) = kobsi(jobs)-1
1010	igrdj(1,2,jobs) = kobsj(jobs)
1011	igrdi(2,1,jobs) = kobsi(jobs)
1012	igrdj(2,1,jobs) = kobsj(jobs)-1
1013	igrdi(2,2,jobs) = kobsi(jobs)
1014	igrdj(2,2,jobs) = kobsj(jobs)
1015
1016	ENDIF
1017
1018	END DO
1019
1020	#if defined key_bdy
1021	! Create a mask grid points in boundary rim
1022	IF (ld_bound_reject) THEN
1023	zbdymask(:,:) = 1.0_wp
1024	DO ji = 1, nb_bdy
1025	DO jj = 1, idx_bdy(ji)%nblen(1)
1026	zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp
1027	ENDDO
1028	ENDDO
1029
1030	CALL obs_int_comm_2d( 2, 2, kobsno, igrdi, igrdj, zbdymask, zbmsk )
1031	ENDIF
1032	#endif
1033
1034	CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, pmask, zgmsk )
1035	CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, plam, zglam )
1036	CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, pphi, zgphi )
1037
1038	DO jobs = 1, kobsno
1039
1040	! Skip bad observations
1041	IF ( kobsqc(jobs) >= kqc_cutoff ) CYCLE
1042
1043	! Flag if the observation falls outside the model spatial domain
1044	IF ( ( pobslam(jobs) < -180. ) &
1045	& .OR. ( pobslam(jobs) > 180. ) &
1046	& .OR. ( pobsphi(jobs) < -90. ) &
1047	& .OR. ( pobsphi(jobs) > 90. ) ) THEN
1048	kobsqc(jobs) = IBSET(kobsqc(jobs),11)
1049	kosdobs = kosdobs + 1
1050	CYCLE
1051	ENDIF
1052
1053	! Flag if the observation falls with a model land cell
1054	IF ( SUM( zgmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN
1055	kobsqc(jobs) = IBSET(kobsqc(jobs),10)
1056	klanobs = klanobs + 1
1057	CYCLE
1058	ENDIF
1059
1060	! Check if this observation is on a grid point
1061
1062	lgridobs = .FALSE.
1063	iig = -1
1064	ijg = -1
1065	DO jj = 1, 2
1066	DO ji = 1, 2
1067	IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) &
1068	& .AND. &
1069	& ( ABS( zglam(ji,jj,jobs) - pobslam(jobs) ) < 1.0e-6_wp ) &
1070	& ) THEN
1071	lgridobs = .TRUE.
1072	iig = ji
1073	ijg = jj
1074	ENDIF
1075	END DO
1076	END DO
1077
1078	IF (lgridobs) THEN
1079	IF (zgmsk(iig,ijg,jobs) == 0.0_wp ) THEN
1080	kobsqc(jobs) = IBSET(kobsqc(jobs),10)
1081	klanobs = klanobs + 1
1082	CYCLE
1083	ENDIF
1084	ENDIF
1085
1086
1087	! Flag if the observation falls is close to land
1088	IF ( MINVAL( zgmsk(1:2,1:2,jobs) ) == 0.0_wp) THEN
1089	knlaobs = knlaobs + 1
1090	IF (ld_nea) THEN
1091	kobsqc(jobs) = IBSET(kobsqc(jobs),9)
1092	CYCLE
1093	ENDIF
1094	ENDIF
1095
1096	#if defined key_bdy
1097	! Flag if the observation falls close to the boundary rim
1098	IF (ld_bound_reject) THEN
1099	IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN
1100	kobsqc(jobs) = kobsqc(jobs) + 15
1101	kbdyobs = kbdyobs + 1
1102	CYCLE
1103	ENDIF
1104	! for observations on the grid...
1105	IF (lgridobs) THEN
1106	IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN
1107	kobsqc(jobs) = kobsqc(jobs) + 15
1108	kbdyobs = kbdyobs + 1
1109	CYCLE
1110	ENDIF
1111	ENDIF
1112	ENDIF
1113	#endif
1114
1115	END DO
1116
1117	END SUBROUTINE obs_coo_spc_2d
1118
1119	SUBROUTINE obs_coo_spc_3d( kprofno, kobsno, kpstart, kpend, &
1120	& kpi, kpj, kpk, &
1121	& kobsi, kobsj, kobsk, &
1122	& pobslam, pobsphi, pobsdep, &
1123	& plam, pphi, pdep, pmask, &
1124	& kpobsqc, kobsqc, kosdobs, &
1125	& klanobs, knlaobs, ld_nea, &
1126	& kbdyobs, ld_bound_reject, &
1127	& kqc_cutoff )
1128	!!----------------------------------------------------------------------
1129	!! * ROUTINE obs_coo_spc_3d *
1130	!!
1131	!! ** Purpose : Check for points outside the domain and land points
1132	!! Reset depth of observation above highest model level
1133	!! to the value of highest model level
1134	!!
1135	!! ** Method : Remove the observations that are outside the model space
1136	!! and time domain or located within model land cells.
1137	!!
1138	!! NB. T and S profile observations lying between the ocean
1139	!! surface and the depth of the first model T point are
1140	!! assigned a depth equal to that of the first model T pt.
1141	!!
1142	!! ** Action :
1143	!!
1144	!! History :
1145	!! ! 2007-01 (K. Mogensen) Rewrite of parts of obs_scr
1146	!! ! 2007-06 (K. Mogensen et al) Reject obs. near land.
1147	!!----------------------------------------------------------------------
1148	!! * Modules used
1149	USE dom_oce, ONLY : & ! Geographical information
1150	& gdepw_1d, &
1151	& gdepw_0, &
1152	#if defined key_vvl
1153	& gdepw_n, &
1154	& gdept_n, &
1155	#endif
1156	& ln_zco, &
1157	& ln_zps, &
1158	& lk_vvl
1159
1160	!! * Arguments
1161	INTEGER, INTENT(IN) :: kprofno ! Number of profiles
1162	INTEGER, INTENT(IN) :: kobsno ! Total number of observations
1163	INTEGER, INTENT(IN) :: kpi ! Number of grid points in (i,j,k)
1164	INTEGER, INTENT(IN) :: kpj
1165	INTEGER, INTENT(IN) :: kpk
1166	INTEGER, DIMENSION(kprofno), INTENT(IN) :: &
1167	& kpstart, & ! Start of individual profiles
1168	& kpend ! End of individual profiles
1169	INTEGER, DIMENSION(kprofno), INTENT(IN) :: &
1170	& kobsi, & ! Observation (i,j) coordinates
1171	& kobsj
1172	INTEGER, DIMENSION(kobsno), INTENT(IN) :: &
1173	& kobsk ! Observation k coordinate
1174	REAL(KIND=wp), DIMENSION(kprofno), INTENT(IN) :: &
1175	& pobslam, & ! Observation (lon,lat) coordinates
1176	& pobsphi
1177	REAL(KIND=wp), DIMENSION(kobsno), INTENT(INOUT) :: &
1178	& pobsdep ! Observation depths
1179	REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) :: &
1180	& plam, pphi ! Model (lon,lat) coordinates
1181	REAL(KIND=wp), DIMENSION(kpk), INTENT(IN) :: &
1182	& pdep ! Model depth coordinates
1183	REAL(KIND=wp), DIMENSION(kpi,kpj,kpk), INTENT(IN) :: &
1184	& pmask ! Land mask array
1185	INTEGER, DIMENSION(kprofno), INTENT(INOUT) :: &
1186	& kpobsqc ! Profile quality control
1187	INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
1188	& kobsqc ! Observation quality control
1189	INTEGER, INTENT(INOUT) :: kosdobs ! Observations outside space domain
1190	INTEGER, INTENT(INOUT) :: klanobs ! Observations within a model land cell
1191	INTEGER, INTENT(INOUT) :: knlaobs ! Observations near land
1192	INTEGER, INTENT(INOUT) :: kbdyobs ! Observations near boundary
1193	LOGICAL, INTENT(IN) :: ld_nea ! Flag observations near land
1194	LOGICAL, INTENT(IN) :: ld_bound_reject ! Flag observations near open boundary
1195	INTEGER, INTENT(IN) :: kqc_cutoff ! Cutoff QC value
1196
1197	!! * Local declarations
1198	REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: &
1199	& zgmsk ! Grid mask
1200	#if defined key_bdy
1201	REAL(KIND=wp), DIMENSION(2,2,kprofno) :: &
1202	& zbmsk ! Boundary mask
1203	REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask
1204	#endif
1205	REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: &
1206	& zgdepw
1207	REAL(KIND=wp), DIMENSION(2,2,kprofno) :: &
1208	& zglam, & ! Model longitude at grid points
1209	& zgphi ! Model latitude at grid points
1210	INTEGER, DIMENSION(2,2,kprofno) :: &
1211	& igrdi, & ! Grid i,j
1212	& igrdj
1213	LOGICAL :: lgridobs ! Is observation on a model grid point.
1214	LOGICAL :: ll_next_to_land ! Is a profile next to land
1215	INTEGER :: iig, ijg ! i,j of observation on model grid point.
1216	INTEGER :: jobs, jobsp, jk, ji, jj
1217
1218	! Get grid point indices
1219
1220	DO jobs = 1, kprofno
1221
1222	! For invalid points use 2,2
1223
1224	IF ( kpobsqc(jobs) >= kqc_cutoff ) THEN
1225
1226	igrdi(1,1,jobs) = 1
1227	igrdj(1,1,jobs) = 1
1228	igrdi(1,2,jobs) = 1
1229	igrdj(1,2,jobs) = 2
1230	igrdi(2,1,jobs) = 2
1231	igrdj(2,1,jobs) = 1
1232	igrdi(2,2,jobs) = 2
1233	igrdj(2,2,jobs) = 2
1234
1235	ELSE
1236
1237	igrdi(1,1,jobs) = kobsi(jobs)-1
1238	igrdj(1,1,jobs) = kobsj(jobs)-1
1239	igrdi(1,2,jobs) = kobsi(jobs)-1
1240	igrdj(1,2,jobs) = kobsj(jobs)
1241	igrdi(2,1,jobs) = kobsi(jobs)
1242	igrdj(2,1,jobs) = kobsj(jobs)-1
1243	igrdi(2,2,jobs) = kobsi(jobs)
1244	igrdj(2,2,jobs) = kobsj(jobs)
1245
1246	ENDIF
1247
1248	END DO
1249
1250	#if defined key_bdy
1251	! Create a mask grid points in boundary rim
1252	IF (ld_bound_reject) THEN
1253	zbdymask(:,:) = 1.0_wp
1254	DO ji = 1, nb_bdy
1255	DO jj = 1, idx_bdy(ji)%nblen(1)
1256	zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp
1257	ENDDO
1258	ENDDO
1259	ENDIF
1260
1261	CALL obs_int_comm_2d( 2, 2, kprofno, igrdi, igrdj, zbdymask, zbmsk )
1262	#endif
1263
1264	CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, pmask, zgmsk )
1265	CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, plam, zglam )
1266	CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, pphi, zgphi )
1267	! Need to know the bathy depth for each observation for sco
1268	CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, fsdepw(:,:,:), &
1269	& zgdepw )
1270
1271	DO jobs = 1, kprofno
1272
1273	! Skip bad profiles
1274	IF ( kpobsqc(jobs) >= kqc_cutoff ) CYCLE
1275
1276	! Check if this observation is on a grid point
1277
1278	lgridobs = .FALSE.
1279	iig = -1
1280	ijg = -1
1281	DO jj = 1, 2
1282	DO ji = 1, 2
1283	IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) &
1284	& .AND. &
1285	& ( ABS( zglam(ji,jj,jobs) - pobslam(jobs) ) < 1.0e-6_wp ) &
1286	& ) THEN
1287	lgridobs = .TRUE.
1288	iig = ji
1289	ijg = jj
1290	ENDIF
1291	END DO
1292	END DO
1293
1294	! Check if next to land
1295	IF ( ANY( zgmsk(1:2,1:2,1,jobs) == 0.0_wp ) ) THEN
1296	ll_next_to_land=.TRUE.
1297	ELSE
1298	ll_next_to_land=.FALSE.
1299	ENDIF
1300
1301	! Reject observations
1302
1303	DO jobsp = kpstart(jobs), kpend(jobs)
1304
1305	! Flag if the observation falls outside the model spatial domain
1306	IF ( ( pobslam(jobs) < -180. ) &
1307	& .OR. ( pobslam(jobs) > 180. ) &
1308	& .OR. ( pobsphi(jobs) < -90. ) &
1309	& .OR. ( pobsphi(jobs) > 90. ) &
1310	& .OR. ( pobsdep(jobsp) < 0.0 ) &
1311	& .OR. ( pobsdep(jobsp) > gdepw_1d(kpk)) ) THEN
1312	kobsqc(jobsp) = IBSET(kobsqc(jobsp),11)
1313	kosdobs = kosdobs + 1
1314	CYCLE
1315	ENDIF
1316
1317	! To check if an observations falls within land there are two cases:
1318	! 1: z-coordibnates, where the check uses the mask
1319	! 2: terrain following (eg s-coordinates),
1320	! where we use the depth of the bottom cell to mask observations
1321
1322	IF( (.NOT. lk_vvl) .AND. ( ln_zps .OR. ln_zco ) ) THEN !(CASE 1)
1323
1324	! Flag if the observation falls with a model land cell
1325	IF ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) &
1326	& == 0.0_wp ) THEN
1327	kobsqc(jobsp) = IBSET(kobsqc(jobsp),10)
1328	klanobs = klanobs + 1
1329	CYCLE
1330	ENDIF
1331
1332	! Flag if the observation is close to land
1333	IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == &
1334	& 0.0_wp) THEN
1335	knlaobs = knlaobs + 1
1336	IF (ld_nea) THEN
1337	kobsqc(jobsp) = IBSET(kobsqc(jobsp),10)
1338	ENDIF
1339	ENDIF
1340
1341	ELSE ! Case 2
1342	! Flag if the observation is deeper than the bathymetry
1343	! Or if it is within the mask
1344	IF ( ANY( zgdepw(1:2,1:2,kpk,jobs) < pobsdep(jobsp) ) &
1345	& .OR. &
1346	& ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) &
1347	& == 0.0_wp) ) THEN
1348	kobsqc(jobsp) = IBSET(kobsqc(jobsp),10)
1349	klanobs = klanobs + 1
1350	CYCLE
1351	ENDIF
1352
1353	! Flag if the observation is close to land
1354	IF ( ll_next_to_land ) THEN
1355	knlaobs = knlaobs + 1
1356	IF (ld_nea) THEN
1357	kobsqc(jobsp) = IBSET(kobsqc(jobsp),10)
1358	ENDIF
1359	ENDIF
1360
1361	ENDIF
1362
1363	! For observations on the grid reject them if their are at
1364	! a masked point
1365
1366	IF (lgridobs) THEN
1367	IF (zgmsk(iig,ijg,kobsk(jobsp)-1,jobs) == 0.0_wp ) THEN
1368	kobsqc(jobsp) = kobsqc(jobsp) + 12
1369	klanobs = klanobs + 1
1370	CYCLE
1371	ENDIF
1372	ENDIF
1373
1374	! Set observation depth equal to that of the first model depth
1375	IF ( pobsdep(jobsp) <= pdep(1) ) THEN
1376	pobsdep(jobsp) = pdep(1)
1377	ENDIF
1378
1379	#if defined key_bdy
1380	! Flag if the observation falls close to the boundary rim
1381	IF (ld_bound_reject) THEN
1382	IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN
1383	kobsqc(jobsp) = kobsqc(jobsp) + 15
1384	kbdyobs = kbdyobs + 1
1385	CYCLE
1386	ENDIF
1387	! for observations on the grid...
1388	IF (lgridobs) THEN
1389	IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN
1390	kobsqc(jobsp) = kobsqc(jobsp) + 15
1391	kbdyobs = kbdyobs + 1
1392	CYCLE
1393	ENDIF
1394	ENDIF
1395	ENDIF
1396	#endif
1397
1398	END DO
1399	END DO
1400
1401	END SUBROUTINE obs_coo_spc_3d
1402
1403	SUBROUTINE obs_pro_rej( profdata, kqc_cutoff )
1404	!!----------------------------------------------------------------------
1405	!! * ROUTINE obs_pro_rej *
1406	!!
1407	!! ** Purpose : Reject all data within a rejected profile
1408	!!
1409	!! ** Method :
1410	!!
1411	!! ** Action :
1412	!!
1413	!! References :
1414	!!
1415	!! History :
1416	!! ! 2007-10 (K. Mogensen) Original code
1417	!!----------------------------------------------------------------------
1418	!! * Modules used
1419	!! * Arguments
1420	TYPE(obs_prof), INTENT(INOUT) :: profdata ! Profile data
1421	INTEGER, INTENT(IN) :: kqc_cutoff ! QC cutoff value
1422
1423	!! * Local declarations
1424	INTEGER :: jprof
1425	INTEGER :: jvar
1426	INTEGER :: jobs
1427
1428	! Loop over profiles
1429
1430	DO jprof = 1, profdata%nprof
1431
1432	IF ( profdata%nqc(jprof) > kqc_cutoff ) THEN
1433
1434	DO jvar = 1, profdata%nvar
1435
1436	DO jobs = profdata%npvsta(jprof,jvar), &
1437	& profdata%npvend(jprof,jvar)
1438
1439	profdata%var(jvar)%nvqc(jobs) = &
1440	& IBSET(profdata%var(jvar)%nvqc(jobs),14)
1441
1442	END DO
1443
1444	END DO
1445
1446	ENDIF
1447
1448	END DO
1449
1450	END SUBROUTINE obs_pro_rej
1451
1452	SUBROUTINE obs_uv_rej( profdata, knumu, knumv, kqc_cutoff )
1453	!!----------------------------------------------------------------------
1454	!! * ROUTINE obs_uv_rej *
1455	!!
1456	!! ** Purpose : Reject u if v is rejected and vice versa
1457	!!
1458	!! ** Method :
1459	!!
1460	!! ** Action :
1461	!!
1462	!! References :
1463	!!
1464	!! History :
1465	!! ! 2009-2 (K. Mogensen) Original code
1466	!!----------------------------------------------------------------------
1467	!! * Modules used
1468	!! * Arguments
1469	TYPE(obs_prof), INTENT(INOUT) :: profdata ! Profile data
1470	INTEGER, INTENT(INOUT) :: knumu ! Number of u rejected
1471	INTEGER, INTENT(INOUT) :: knumv ! Number of v rejected
1472	INTEGER, INTENT(IN) :: kqc_cutoff ! QC cutoff value
1473
1474	!! * Local declarations
1475	INTEGER :: jprof
1476	INTEGER :: jvar
1477	INTEGER :: jobs
1478
1479	! Loop over profiles
1480
1481	DO jprof = 1, profdata%nprof
1482
1483	IF ( ( profdata%npvsta(jprof,1) /= profdata%npvsta(jprof,2) ) .OR. &
1484	& ( profdata%npvend(jprof,1) /= profdata%npvend(jprof,2) ) ) THEN
1485
1486	CALL ctl_stop('U,V profiles inconsistent in obs_uv_rej')
1487	RETURN
1488
1489	ENDIF
1490
1491	DO jobs = profdata%npvsta(jprof,1), profdata%npvend(jprof,1)
1492
1493	IF ( ( profdata%var(1)%nvqc(jobs) > kqc_cutoff ) .AND. &
1494	& ( profdata%var(2)%nvqc(jobs) <= kqc_cutoff) ) THEN
1495	profdata%var(2)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15)
1496	knumv = knumv + 1
1497	ENDIF
1498	IF ( ( profdata%var(2)%nvqc(jobs) > kqc_cutoff ) .AND. &
1499	& ( profdata%var(1)%nvqc(jobs) <= kqc_cutoff) ) THEN
1500	profdata%var(1)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15)
1501	knumu = knumu + 1
1502	ENDIF
1503
1504	END DO
1505
1506	END DO
1507
1508	END SUBROUTINE obs_uv_rej
1509
1510	END MODULE obs_prep

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

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

Context Navigation

source: branches/UKMO/obs_oper_do_not_assim_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 7932

Download in other formats: