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 NEMO/trunk/src/OCE/OBS – NEMO

Context Navigation

source: NEMO/trunk/src/OCE/OBS/obs_prep.F90 @ 14587

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

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

Download in other formats: