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

Context Navigation

source: branches/UKMO/dev_r5518_obs_oper_update_utils335/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 12498

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

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

Download in other formats: