Context Navigation

obs_prep.F90 @ 7773

Last change on this file since 7773 was 7773, checked in by mattmartin, 7 years ago

Committing updates after doing the following:

merging the branch dev_r4650_general_vert_coord_obsoper@7763 into this branch
updating it so that the following OBS changes were implemented correctly on top of the simplification changes:
- generalised vertical coordinate for profile obs. This was done so that is now the default option.
- sst bias correction implemented with the new simplified obs code.
- included the biogeochemical obs types int he new simplified obs code.
- included the changes to exclude obs in the boundary for limited area models
- included other changes for the efficiency of the obs operator to remove global arrays.

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

Download in other formats: