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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

obs_prep.F90 in branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

Context Navigation

source: branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 9186

Last change on this file since 9186 was 9186, checked in by dford, 6 years ago
Initial implementation of 3D biogeochemistry observation operator.
File size: 61.6 KB

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

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

Download in other formats: