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

Context Navigation

source: branches/UKMO/dev_r5518_obs_oper_update_sit/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 11932

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

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

Download in other formats: