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

Context Navigation

source: branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 5704

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

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

Download in other formats: