Context Navigation

obs_prep.F90 @ 15540

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

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

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

Context Navigation

source: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/OBS/obs_prep.F90 @ 15540

Download in other formats: