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

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

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

Context Navigation

source: branches/UKMO/dev_rev5518_OBS_DoNotAssim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 7841

Last change on this file since 7841 was 7841, checked in by jwhile, 7 years ago
Added "Do not Assimlate" funtionality to OBS code
Property svn:keywords set to `Id`
File size: 82.8 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_pro : First level check and screening of T/S profiles
10	!! obs_pre_sla : First level check and screening of SLA observations
11	!! obs_pre_sst : First level check and screening of SLA observations
12	!! obs_pre_seaice : First level check and screening of sea ice observations
13	!! obs_pre_vel : First level check and screening of velocity obs.
14	!! obs_scr : Basic screening of the observations
15	!! obs_coo_tim : Compute number of time steps to the observation time
16	!! obs_sor : Sort the observation arrays
17	!!---------------------------------------------------------------------
18	!! * Modules used
19	USE par_kind, ONLY : & ! Precision variables
20	& wp
21	USE in_out_manager ! I/O manager
22	USE obs_profiles_def ! Definitions for storage arrays for profiles
23	USE obs_surf_def ! Definitions for storage arrays for surface data
24	USE obs_mpp, ONLY : & ! MPP support routines for observation diagnostics
25	& obs_mpp_sum_integer, &
26	& obs_mpp_sum_integers
27	USE obs_inter_sup ! Interpolation support
28	USE obs_oper ! Observation operators
29	USE lib_mpp, ONLY : &
30	& ctl_warn, ctl_stop
31
32	IMPLICIT NONE
33
34	!! * Routine accessibility
35	PRIVATE
36
37	PUBLIC &
38	& obs_pre_pro, & ! First level check and screening of profiles
39	& obs_pre_sla, & ! First level check and screening of SLA data
40	& obs_pre_sst, & ! First level check and screening of SLA data
41	& obs_pre_seaice, & ! First level check and screening of sea ice data
42	& obs_pre_vel, & ! First level check and screening of velocity profiles
43	& calc_month_len ! Calculate the number of days in the months of a year
44
45	!!----------------------------------------------------------------------
46	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
47	!! $Id$
48	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
49	!!----------------------------------------------------------------------
50
51	CONTAINS
52
53	SUBROUTINE obs_pre_pro( profdata, prodatqc, ld_t3d, ld_s3d, ld_nea, &
54	& kdailyavtypes, kqc_cutoff )
55	!!----------------------------------------------------------------------
56	!! * ROUTINE obs_pre_pro *
57	!!
58	!! ** Purpose : First level check and screening of T and S profiles
59	!!
60	!! ** Method : First level check and screening of T and S profiles
61	!!
62	!! ** Action :
63	!!
64	!! References :
65	!!
66	!! History :
67	!! ! 2007-01 (K. Mogensen) Merge of obs_pre_t3d and obs_pre_s3d
68	!! ! 2007-03 (K. Mogensen) General handling of profiles
69	!! ! 2007-06 (K. Mogensen et al) Reject obs. near land.
70	!!----------------------------------------------------------------------
71	!! * Modules used
72	USE domstp ! Domain: set the time-step
73	USE par_oce ! Ocean parameters
74	USE dom_oce, ONLY : & ! Geographical information
75	& glamt, &
76	& gphit, &
77	& gdept_1d,&
78	& tmask, &
79	& nproc
80	!! * Arguments
81	TYPE(obs_prof), INTENT(INOUT) :: profdata ! Full set of profile data
82	TYPE(obs_prof), INTENT(INOUT) :: prodatqc ! Subset of profile data not failing screening
83	LOGICAL, INTENT(IN) :: ld_t3d ! Switch for temperature
84	LOGICAL, INTENT(IN) :: ld_s3d ! Switch for salinity
85	LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land
86	INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
87	& kdailyavtypes! Types for daily averages
88	INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff ! cut off for QC value
89	!! * Local declarations
90	INTEGER :: iqc_cutoff = 255 ! cut off for QC value
91	INTEGER :: iyea0 ! Initial date
92	INTEGER :: imon0 ! - (year, month, day, hour, minute)
93	INTEGER :: iday0
94	INTEGER :: ihou0
95	INTEGER :: imin0
96	INTEGER :: icycle ! Current assimilation cycle
97	! Counters for observations that
98	INTEGER :: iotdobs ! - outside time domain
99	INTEGER :: iosdtobs ! - outside space domain (temperature)
100	INTEGER :: iosdsobs ! - outside space domain (salinity)
101	INTEGER :: ilantobs ! - within a model land cell (temperature)
102	INTEGER :: ilansobs ! - within a model land cell (salinity)
103	INTEGER :: inlatobs ! - close to land (temperature)
104	INTEGER :: inlasobs ! - close to land (salinity)
105	INTEGER :: igrdobs ! - fail the grid search
106	! Global counters for observations that
107	INTEGER :: iotdobsmpp ! - outside time domain
108	INTEGER :: iosdtobsmpp ! - outside space domain (temperature)
109	INTEGER :: iosdsobsmpp ! - outside space domain (salinity)
110	INTEGER :: ilantobsmpp ! - within a model land cell (temperature)
111	INTEGER :: ilansobsmpp ! - within a model land cell (salinity)
112	INTEGER :: inlatobsmpp ! - close to land (temperature)
113	INTEGER :: inlasobsmpp ! - close to land (salinity)
114	INTEGER :: igrdobsmpp ! - fail the grid search
115	TYPE(obs_prof_valid) :: llvalid ! Profile selection
116	TYPE(obs_prof_valid), DIMENSION(profdata%nvar) :: &
117	& llvvalid ! T,S selection
118	INTEGER :: jvar ! Variable loop variable
119	INTEGER :: jobs ! Obs. loop variable
120	INTEGER :: jstp ! Time loop variable
121	INTEGER :: inrc ! Time index variable
122
123	IF(lwp) WRITE(numout,*)'obs_pre_pro : Preparing the profile observations...'
124
125	! Initial date initialization (year, month, day, hour, minute)
126	iyea0 = ndate0 / 10000
127	imon0 = ( ndate0 - iyea0 * 10000 ) / 100
128	iday0 = ndate0 - iyea0 * 10000 - imon0 * 100
129	ihou0 = 0
130	imin0 = 0
131
132	icycle = no ! Assimilation cycle
133
134	! Diagnotics counters for various failures.
135
136	iotdobs = 0
137	igrdobs = 0
138	iosdtobs = 0
139	iosdsobs = 0
140	ilantobs = 0
141	ilansobs = 0
142	inlatobs = 0
143	inlasobs = 0
144
145	! Set QC cutoff to optional value if provided
146	IF ( PRESENT(kqc_cutoff) ) iqc_cutoff=kqc_cutoff
147
148	! -----------------------------------------------------------------------
149	! Find time coordinate for profiles
150	! -----------------------------------------------------------------------
151
152	IF ( PRESENT(kdailyavtypes) ) THEN
153	CALL obs_coo_tim_prof( icycle, &
154	& iyea0, imon0, iday0, ihou0, imin0, &
155	& profdata%nprof, profdata%nyea, profdata%nmon, &
156	& profdata%nday, profdata%nhou, profdata%nmin, &
157	& profdata%ntyp, profdata%nqc, profdata%mstp, &
158	& iotdobs, kdailyavtypes = kdailyavtypes, &
159	& kqc_cutoff = iqc_cutoff )
160	ELSE
161	CALL obs_coo_tim_prof( icycle, &
162	& iyea0, imon0, iday0, ihou0, imin0, &
163	& profdata%nprof, profdata%nyea, profdata%nmon, &
164	& profdata%nday, profdata%nhou, profdata%nmin, &
165	& profdata%ntyp, profdata%nqc, profdata%mstp, &
166	& iotdobs, kqc_cutoff = iqc_cutoff )
167	ENDIF
168	CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
169
170	! -----------------------------------------------------------------------
171	! Check for profiles failing the grid search
172	! -----------------------------------------------------------------------
173
174	CALL obs_coo_grd( profdata%nprof, profdata%mi, profdata%mj, &
175	& profdata%nqc, igrdobs )
176
177	CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
178
179	! -----------------------------------------------------------------------
180	! Reject all observations for profiles with nqc > iqc_cutoff
181	! -----------------------------------------------------------------------
182
183	CALL obs_pro_rej( profdata, iqc_cutoff )
184
185	! -----------------------------------------------------------------------
186	! Check for land points. This includes points below the model
187	! bathymetry so this is done for every point in the profile
188	! -----------------------------------------------------------------------
189
190	! Temperature
191
192	CALL obs_coo_spc_3d( profdata%nprof, profdata%nvprot(1), &
193	& profdata%npvsta(:,1), profdata%npvend(:,1), &
194	& jpi, jpj, &
195	& jpk, &
196	& profdata%mi, profdata%mj, &
197	& profdata%var(1)%mvk, &
198	& profdata%rlam, profdata%rphi, &
199	& profdata%var(1)%vdep, &
200	& glamt, gphit, &
201	& gdept_1d, tmask, &
202	& profdata%nqc, profdata%var(1)%nvqc, &
203	& iosdtobs, ilantobs, &
204	& inlatobs, ld_nea, &
205	& iqc_cutoff )
206
207	CALL obs_mpp_sum_integer( iosdtobs, iosdtobsmpp )
208	CALL obs_mpp_sum_integer( ilantobs, ilantobsmpp )
209	CALL obs_mpp_sum_integer( inlatobs, inlatobsmpp )
210
211	! Salinity
212
213	CALL obs_coo_spc_3d( profdata%nprof, profdata%nvprot(2), &
214	& profdata%npvsta(:,2), profdata%npvend(:,2), &
215	& jpi, jpj, &
216	& jpk, &
217	& profdata%mi, profdata%mj, &
218	& profdata%var(2)%mvk, &
219	& profdata%rlam, profdata%rphi, &
220	& profdata%var(2)%vdep, &
221	& glamt, gphit, &
222	& gdept_1d, tmask, &
223	& profdata%nqc, profdata%var(2)%nvqc, &
224	& iosdsobs, ilansobs, &
225	& inlasobs, ld_nea, &
226	& iqc_cutoff )
227
228	CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
229	CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
230	CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
231
232	! -----------------------------------------------------------------------
233	! Copy useful data from the profdata data structure to
234	! the prodatqc data structure
235	! -----------------------------------------------------------------------
236
237	! Allocate the selection arrays
238
239	ALLOCATE( llvalid%luse(profdata%nprof) )
240	DO jvar = 1,profdata%nvar
241	ALLOCATE( llvvalid(jvar)%luse(profdata%nvprot(jvar)) )
242	END DO
243
244	! We want all data which has qc flags <= iqc_cutoff
245
246	llvalid%luse(:) = ( profdata%nqc(:) <= iqc_cutoff )
247	DO jvar = 1,profdata%nvar
248	llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= iqc_cutoff )
249	END DO
250
251	! The actual copying
252
253	CALL obs_prof_compress( profdata, prodatqc, .TRUE., numout, &
254	& lvalid=llvalid, lvvalid=llvvalid )
255
256	! Dellocate the selection arrays
257	DEALLOCATE( llvalid%luse )
258	DO jvar = 1,profdata%nvar
259	DEALLOCATE( llvvalid(jvar)%luse )
260	END DO
261
262	! -----------------------------------------------------------------------
263	! Print information about what observations are left after qc
264	! -----------------------------------------------------------------------
265
266	! Update the total observation counter array
267
268	IF(lwp) THEN
269	WRITE(numout,*)
270	WRITE(numout,*) 'obs_pre_pro :'
271	WRITE(numout,*) '~~~~~~~~~~~'
272	WRITE(numout,*)
273	WRITE(numout,*) ' Profiles outside time domain = ', &
274	& iotdobsmpp
275	WRITE(numout,*) ' Remaining profiles that failed grid search = ', &
276	& igrdobsmpp
277	WRITE(numout,*) ' Remaining T data outside space domain = ', &
278	& iosdtobsmpp
279	WRITE(numout,*) ' Remaining T data at land points = ', &
280	& ilantobsmpp
281	IF (ld_nea) THEN
282	WRITE(numout,*) ' Remaining T data near land points (removed) = ',&
283	& inlatobsmpp
284	ELSE
285	WRITE(numout,*) ' Remaining T data near land points (kept) = ',&
286	& inlatobsmpp
287	ENDIF
288	WRITE(numout,*) ' T data accepted = ', &
289	& prodatqc%nvprotmpp(1)
290	WRITE(numout,*) ' Remaining S data outside space domain = ', &
291	& iosdsobsmpp
292	WRITE(numout,*) ' Remaining S data at land points = ', &
293	& ilansobsmpp
294	IF (ld_nea) THEN
295	WRITE(numout,*) ' Remaining S data near land points (removed) = ',&
296	& inlasobsmpp
297	ELSE
298	WRITE(numout,*) ' Remaining S data near land points (kept) = ',&
299	& inlasobsmpp
300	ENDIF
301	WRITE(numout,*) ' S data accepted = ', &
302	& prodatqc%nvprotmpp(2)
303
304	WRITE(numout,*)
305	WRITE(numout,*) ' Number of observations per time step :'
306	WRITE(numout,*)
307	WRITE(numout,997)
308	WRITE(numout,998)
309	ENDIF
310
311	DO jobs = 1, prodatqc%nprof
312	inrc = prodatqc%mstp(jobs) + 2 - nit000
313	prodatqc%npstp(inrc) = prodatqc%npstp(inrc) + 1
314	DO jvar = 1, prodatqc%nvar
315	IF ( prodatqc%npvend(jobs,jvar) > 0 ) THEN
316	prodatqc%nvstp(inrc,jvar) = prodatqc%nvstp(inrc,jvar) + &
317	& ( prodatqc%npvend(jobs,jvar) - &
318	& prodatqc%npvsta(jobs,jvar) + 1 )
319	ENDIF
320	END DO
321	END DO
322
323
324	CALL obs_mpp_sum_integers( prodatqc%npstp, prodatqc%npstpmpp, &
325	& nitend - nit000 + 2 )
326	DO jvar = 1, prodatqc%nvar
327	CALL obs_mpp_sum_integers( prodatqc%nvstp(:,jvar), &
328	& prodatqc%nvstpmpp(:,jvar), &
329	& nitend - nit000 + 2 )
330	END DO
331
332	IF ( lwp ) THEN
333	DO jstp = nit000 - 1, nitend
334	inrc = jstp - nit000 + 2
335	WRITE(numout,999) jstp, prodatqc%npstpmpp(inrc), &
336	& prodatqc%nvstpmpp(inrc,1), &
337	& prodatqc%nvstpmpp(inrc,2)
338	END DO
339	ENDIF
340
341	997 FORMAT(10X,'Time step',5X,'Profiles',5X,'Temperature',5X,'Salinity')
342	998 FORMAT(10X,'---------',5X,'--------',5X,'-----------',5X,'--------')
343	999 FORMAT(10X,I9,5X,I8,5X,I11,5X,I8)
344
345	END SUBROUTINE obs_pre_pro
346
347	SUBROUTINE obs_pre_sla( sladata, sladatqc, ld_sla, ld_nea, kqc_cutoff )
348	!!----------------------------------------------------------------------
349	!! * ROUTINE obs_pre_sla *
350	!!
351	!! ** Purpose : First level check and screening of SLA observations
352	!!
353	!! ** Method : First level check and screening of SLA observations
354	!!
355	!! ** Action :
356	!!
357	!! References :
358	!!
359	!! History :
360	!! ! 2007-03 (A. Weaver, K. Mogensen) Original
361	!! ! 2007-06 (K. Mogensen et al) Reject obs. near land.
362	!!----------------------------------------------------------------------
363	!! * Modules used
364	USE domstp ! Domain: set the time-step
365	USE par_oce ! Ocean parameters
366	USE dom_oce, ONLY : & ! Geographical information
367	& glamt, &
368	& gphit, &
369	& tmask, &
370	& nproc
371	!! * Arguments
372	TYPE(obs_surf), INTENT(INOUT) :: sladata ! Full set of SLA data
373	TYPE(obs_surf), INTENT(INOUT) :: sladatqc ! Subset of SLA data not failing screening
374	LOGICAL, INTENT(IN) :: ld_sla ! Switch for SLA data
375	LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land
376	INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff ! cut off for QC value
377	!! * Local declarations
378	INTEGER :: iqc_cutoff = 255 ! cut off for QC value
379	INTEGER :: iyea0 ! Initial date
380	INTEGER :: imon0 ! - (year, month, day, hour, minute)
381	INTEGER :: iday0
382	INTEGER :: ihou0
383	INTEGER :: imin0
384	INTEGER :: icycle ! Current assimilation cycle
385	! Counters for observations that
386	INTEGER :: iotdobs ! - outside time domain
387	INTEGER :: iosdsobs ! - outside space domain
388	INTEGER :: ilansobs ! - within a model land cell
389	INTEGER :: inlasobs ! - close to land
390	INTEGER :: igrdobs ! - fail the grid search
391	! Global counters for observations that
392	INTEGER :: iotdobsmpp ! - outside time domain
393	INTEGER :: iosdsobsmpp ! - outside space domain
394	INTEGER :: ilansobsmpp ! - within a model land cell
395	INTEGER :: inlasobsmpp ! - close to land
396	INTEGER :: igrdobsmpp ! - fail the grid search
397	LOGICAL, DIMENSION(:), ALLOCATABLE :: &
398	& llvalid ! SLA data selection
399	INTEGER :: jobs ! Obs. loop variable
400	INTEGER :: jstp ! Time loop variable
401	INTEGER :: inrc ! Time index variable
402
403	IF(lwp) WRITE(numout,*)'obs_pre_sla : Preparing the SLA observations...'
404
405	! Initial date initialization (year, month, day, hour, minute)
406	iyea0 = ndate0 / 10000
407	imon0 = ( ndate0 - iyea0 * 10000 ) / 100
408	iday0 = ndate0 - iyea0 * 10000 - imon0 * 100
409	ihou0 = 0
410	imin0 = 0
411
412	icycle = no ! Assimilation cycle
413
414	! Diagnotics counters for various failures.
415
416	iotdobs = 0
417	igrdobs = 0
418	iosdsobs = 0
419	ilansobs = 0
420	inlasobs = 0
421
422	! Set QC cutoff to optional value if provided
423
424	IF ( PRESENT(kqc_cutoff) ) iqc_cutoff=kqc_cutoff
425
426	! -----------------------------------------------------------------------
427	! Find time coordinate for SLA data
428	! -----------------------------------------------------------------------
429
430	CALL obs_coo_tim( icycle, &
431	& iyea0, imon0, iday0, ihou0, imin0, &
432	& sladata%nsurf, sladata%nyea, sladata%nmon, &
433	& sladata%nday, sladata%nhou, sladata%nmin, &
434	& sladata%nqc, sladata%mstp, iotdobs )
435
436	CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
437
438	! -----------------------------------------------------------------------
439	! Check for SLA data failing the grid search
440	! -----------------------------------------------------------------------
441
442	CALL obs_coo_grd( sladata%nsurf, sladata%mi, sladata%mj, &
443	& sladata%nqc, igrdobs )
444
445	CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
446
447	! -----------------------------------------------------------------------
448	! Check for land points.
449	! -----------------------------------------------------------------------
450
451	CALL obs_coo_spc_2d( sladata%nsurf, &
452	& jpi, jpj, &
453	& sladata%mi, sladata%mj, &
454	& sladata%rlam, sladata%rphi, &
455	& glamt, gphit, &
456	& tmask(:,:,1), sladata%nqc, &
457	& iosdsobs, ilansobs, &
458	& inlasobs, ld_nea, &
459	& iqc_cutoff )
460
461	CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
462	CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
463	CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
464
465	! -----------------------------------------------------------------------
466	! Copy useful data from the sladata data structure to
467	! the sladatqc data structure
468	! -----------------------------------------------------------------------
469
470	! Allocate the selection arrays
471
472	ALLOCATE( llvalid(sladata%nsurf) )
473
474	! We want all data which has qc flags <= iqc_cutoff
475
476	llvalid(:) = ( sladata%nqc(:) <= iqc_cutoff )
477
478	! The actual copying
479
480	CALL obs_surf_compress( sladata, sladatqc, .TRUE., numout, &
481	& lvalid=llvalid )
482
483	! Dellocate the selection arrays
484	DEALLOCATE( llvalid )
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	WRITE(numout,*)
494	WRITE(numout,*) 'obs_pre_sla :'
495	WRITE(numout,*) '~~~~~~~~~~~'
496	WRITE(numout,*)
497	WRITE(numout,*) ' SLA data outside time domain = ', &
498	& iotdobsmpp
499	WRITE(numout,*) ' Remaining SLA data that failed grid search = ', &
500	& igrdobsmpp
501	WRITE(numout,*) ' Remaining SLA data outside space domain = ', &
502	& iosdsobsmpp
503	WRITE(numout,*) ' Remaining SLA data at land points = ', &
504	& ilansobsmpp
505	IF (ld_nea) THEN
506	WRITE(numout,*) ' Remaining SLA data near land points (removed) = ', &
507	& inlasobsmpp
508	ELSE
509	WRITE(numout,*) ' Remaining SLA data near land points (kept) = ', &
510	& inlasobsmpp
511	ENDIF
512	WRITE(numout,*) ' SLA data accepted = ', &
513	& sladatqc%nsurfmpp
514
515	WRITE(numout,*)
516	WRITE(numout,*) ' Number of observations per time step :'
517	WRITE(numout,*)
518	WRITE(numout,1997)
519	WRITE(numout,1998)
520	ENDIF
521
522	DO jobs = 1, sladatqc%nsurf
523	inrc = sladatqc%mstp(jobs) + 2 - nit000
524	sladatqc%nsstp(inrc) = sladatqc%nsstp(inrc) + 1
525	END DO
526
527	CALL obs_mpp_sum_integers( sladatqc%nsstp, sladatqc%nsstpmpp, &
528	& nitend - nit000 + 2 )
529
530	IF ( lwp ) THEN
531	DO jstp = nit000 - 1, nitend
532	inrc = jstp - nit000 + 2
533	WRITE(numout,1999) jstp, sladatqc%nsstpmpp(inrc)
534	END DO
535	ENDIF
536
537	1997 FORMAT(10X,'Time step',5X,'Sea level anomaly')
538	1998 FORMAT(10X,'---------',5X,'-----------------')
539	1999 FORMAT(10X,I9,5X,I17)
540
541	END SUBROUTINE obs_pre_sla
542
543	SUBROUTINE obs_pre_sst( sstdata, sstdatqc, ld_sst, ld_nea, kqc_cutoff )
544	!!----------------------------------------------------------------------
545	!! * ROUTINE obs_pre_sst *
546	!!
547	!! ** Purpose : First level check and screening of SST observations
548	!!
549	!! ** Method : First level check and screening of SST observations
550	!!
551	!! ** Action :
552	!!
553	!! References :
554	!!
555	!! History :
556	!! ! 2007-03 (S. Ricci) SST data preparation
557	!!----------------------------------------------------------------------
558	!! * Modules used
559	USE domstp ! Domain: set the time-step
560	USE par_oce ! Ocean parameters
561	USE dom_oce, ONLY : & ! Geographical information
562	& glamt, &
563	& gphit, &
564	& tmask, &
565	& nproc
566	!! * Arguments
567	TYPE(obs_surf), INTENT(INOUT) :: sstdata ! Full set of SST data
568	TYPE(obs_surf), INTENT(INOUT) :: sstdatqc ! Subset of SST data not failing screening
569	LOGICAL :: ld_sst ! Switch for SST data
570	LOGICAL :: ld_nea ! Switch for rejecting observation near land
571	INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff ! cut off for QC value
572	!! * Local declarations
573	INTEGER :: iqc_cutoff = 255 ! cut off for QC value
574	INTEGER :: iyea0 ! Initial date
575	INTEGER :: imon0 ! - (year, month, day, hour, minute)
576	INTEGER :: iday0
577	INTEGER :: ihou0
578	INTEGER :: imin0
579	INTEGER :: icycle ! Current assimilation cycle
580	! Counters for observations that
581	INTEGER :: iotdobs ! - outside time domain
582	INTEGER :: iosdsobs ! - outside space domain
583	INTEGER :: ilansobs ! - within a model land cell
584	INTEGER :: inlasobs ! - close to land
585	INTEGER :: igrdobs ! - fail the grid search
586	! Global counters for observations that
587	INTEGER :: iotdobsmpp ! - outside time domain
588	INTEGER :: iosdsobsmpp ! - outside space domain
589	INTEGER :: ilansobsmpp ! - within a model land cell
590	INTEGER :: inlasobsmpp ! - close to land
591	INTEGER :: igrdobsmpp ! - fail the grid search
592	LOGICAL, DIMENSION(:), ALLOCATABLE :: &
593	& llvalid ! SST data selection
594	INTEGER :: jobs ! Obs. loop variable
595	INTEGER :: jstp ! Time loop variable
596	INTEGER :: inrc ! Time index variable
597
598	IF(lwp) WRITE(numout,*)'obs_pre_sst : Preparing the SST observations...'
599
600	! Initial date initialization (year, month, day, hour, minute)
601	iyea0 = ndate0 / 10000
602	imon0 = ( ndate0 - iyea0 * 10000 ) / 100
603	iday0 = ndate0 - iyea0 * 10000 - imon0 * 100
604	ihou0 = 0
605	imin0 = 0
606
607	icycle = no ! Assimilation cycle
608
609	! Diagnotics counters for various failures.
610
611	iotdobs = 0
612	igrdobs = 0
613	iosdsobs = 0
614	ilansobs = 0
615	inlasobs = 0
616
617	! Set QC cutoff to optional value if provided
618
619	IF ( PRESENT(kqc_cutoff) ) iqc_cutoff=kqc_cutoff
620
621	! -----------------------------------------------------------------------
622	! Find time coordinate for SST data
623	! -----------------------------------------------------------------------
624
625	CALL obs_coo_tim( icycle, &
626	& iyea0, imon0, iday0, ihou0, imin0, &
627	& sstdata%nsurf, sstdata%nyea, sstdata%nmon, &
628	& sstdata%nday, sstdata%nhou, sstdata%nmin, &
629	& sstdata%nqc, sstdata%mstp, iotdobs )
630	CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
631	! -----------------------------------------------------------------------
632	! Check for SST data failing the grid search
633	! -----------------------------------------------------------------------
634
635	CALL obs_coo_grd( sstdata%nsurf, sstdata%mi, sstdata%mj, &
636	& sstdata%nqc, igrdobs )
637	CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
638
639	! -----------------------------------------------------------------------
640	! Check for land points.
641	! -----------------------------------------------------------------------
642
643	CALL obs_coo_spc_2d( sstdata%nsurf, &
644	& jpi, jpj, &
645	& sstdata%mi, sstdata%mj, &
646	& sstdata%rlam, sstdata%rphi, &
647	& glamt, gphit, &
648	& tmask(:,:,1), sstdata%nqc, &
649	& iosdsobs, ilansobs, &
650	& inlasobs, ld_nea, &
651	& iqc_cutoff )
652
653	CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
654	CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
655	CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
656
657	! -----------------------------------------------------------------------
658	! Copy useful data from the sstdata data structure to
659	! the sstdatqc data structure
660	! -----------------------------------------------------------------------
661
662	! Allocate the selection arrays
663
664	ALLOCATE( llvalid(sstdata%nsurf) )
665
666	! We want all data which has qc flags <= iqc_cutoff
667
668	llvalid(:) = ( sstdata%nqc(:) <= iqc_cutoff )
669
670	! The actual copying
671
672	CALL obs_surf_compress( sstdata, sstdatqc, .TRUE., numout, &
673	& lvalid=llvalid )
674
675	! Dellocate the selection arrays
676	DEALLOCATE( llvalid )
677
678	! -----------------------------------------------------------------------
679	! Print information about what observations are left after qc
680	! -----------------------------------------------------------------------
681
682	! Update the total observation counter array
683
684	IF(lwp) THEN
685	WRITE(numout,*)
686	WRITE(numout,*) 'obs_pre_sst :'
687	WRITE(numout,*) '~~~~~~~~~~~'
688	WRITE(numout,*)
689	WRITE(numout,*) ' SST data outside time domain = ', &
690	& iotdobsmpp
691	WRITE(numout,*) ' Remaining SST data that failed grid search = ', &
692	& igrdobsmpp
693	WRITE(numout,*) ' Remaining SST data outside space domain = ', &
694	& iosdsobsmpp
695	WRITE(numout,*) ' Remaining SST data at land points = ', &
696	& ilansobsmpp
697	IF (ld_nea) THEN
698	WRITE(numout,*) ' Remaining SST data near land points (removed) = ', &
699	& inlasobsmpp
700	ELSE
701	WRITE(numout,*) ' Remaining SST data near land points (kept) = ', &
702	& inlasobsmpp
703	ENDIF
704	WRITE(numout,*) ' SST data accepted = ', &
705	& sstdatqc%nsurfmpp
706
707	WRITE(numout,*)
708	WRITE(numout,*) ' Number of observations per time step :'
709	WRITE(numout,*)
710	WRITE(numout,1997)
711	WRITE(numout,1998)
712	ENDIF
713
714	DO jobs = 1, sstdatqc%nsurf
715	inrc = sstdatqc%mstp(jobs) + 2 - nit000
716	sstdatqc%nsstp(inrc) = sstdatqc%nsstp(inrc) + 1
717	END DO
718
719	CALL obs_mpp_sum_integers( sstdatqc%nsstp, sstdatqc%nsstpmpp, &
720	& nitend - nit000 + 2 )
721
722	IF ( lwp ) THEN
723	DO jstp = nit000 - 1, nitend
724	inrc = jstp - nit000 + 2
725	WRITE(numout,1999) jstp, sstdatqc%nsstpmpp(inrc)
726	END DO
727	ENDIF
728
729	1997 FORMAT(10X,'Time step',5X,'Sea surface temperature')
730	1998 FORMAT(10X,'---------',5X,'-----------------')
731	1999 FORMAT(10X,I9,5X,I17)
732
733	END SUBROUTINE obs_pre_sst
734
735	SUBROUTINE obs_pre_seaice( seaicedata, seaicedatqc, ld_seaice, ld_nea, kqc_cutoff )
736	!!----------------------------------------------------------------------
737	!! * ROUTINE obs_pre_seaice *
738	!!
739	!! ** Purpose : First level check and screening of Sea Ice observations
740	!!
741	!! ** Method : First level check and screening of Sea Ice observations
742	!!
743	!! ** Action :
744	!!
745	!! References :
746	!!
747	!! History :
748	!! ! 2007-11 (D. Lea) based on obs_pre_sst
749	!!----------------------------------------------------------------------
750	!! * Modules used
751	USE domstp ! Domain: set the time-step
752	USE par_oce ! Ocean parameters
753	USE dom_oce, ONLY : & ! Geographical information
754	& glamt, &
755	& gphit, &
756	& tmask, &
757	& nproc
758	!! * Arguments
759	TYPE(obs_surf), INTENT(INOUT) :: seaicedata ! Full set of Sea Ice data
760	TYPE(obs_surf), INTENT(INOUT) :: seaicedatqc ! Subset of sea ice data not failing screening
761	LOGICAL :: ld_seaice ! Switch for sea ice data
762	LOGICAL :: ld_nea ! Switch for rejecting observation near land
763	INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff ! cut off for QC value
764	!! * Local declarations
765	INTEGER :: iqc_cutoff = 255 ! cut off for QC value
766	INTEGER :: iyea0 ! Initial date
767	INTEGER :: imon0 ! - (year, month, day, hour, minute)
768	INTEGER :: iday0
769	INTEGER :: ihou0
770	INTEGER :: imin0
771	INTEGER :: icycle ! Current assimilation cycle
772	! Counters for observations that
773	INTEGER :: iotdobs ! - outside time domain
774	INTEGER :: iosdsobs ! - outside space domain
775	INTEGER :: ilansobs ! - within a model land cell
776	INTEGER :: inlasobs ! - close to land
777	INTEGER :: igrdobs ! - fail the grid search
778	! Global counters for observations that
779	INTEGER :: iotdobsmpp ! - outside time domain
780	INTEGER :: iosdsobsmpp ! - outside space domain
781	INTEGER :: ilansobsmpp ! - within a model land cell
782	INTEGER :: inlasobsmpp ! - close to land
783	INTEGER :: igrdobsmpp ! - fail the grid search
784	LOGICAL, DIMENSION(:), ALLOCATABLE :: &
785	& llvalid ! data selection
786	INTEGER :: jobs ! Obs. loop variable
787	INTEGER :: jstp ! Time loop variable
788	INTEGER :: inrc ! Time index variable
789
790	IF (lwp) WRITE(numout,*)'obs_pre_seaice : Preparing the sea ice observations...'
791
792	! Initial date initialization (year, month, day, hour, minute)
793	iyea0 = ndate0 / 10000
794	imon0 = ( ndate0 - iyea0 * 10000 ) / 100
795	iday0 = ndate0 - iyea0 * 10000 - imon0 * 100
796	ihou0 = 0
797	imin0 = 0
798
799	icycle = no ! Assimilation cycle
800
801	! Diagnotics counters for various failures.
802
803	iotdobs = 0
804	igrdobs = 0
805	iosdsobs = 0
806	ilansobs = 0
807	inlasobs = 0
808
809	! Set QC cutoff to optional value if provided
810
811	IF ( PRESENT(kqc_cutoff) ) iqc_cutoff=kqc_cutoff
812
813	! -----------------------------------------------------------------------
814	! Find time coordinate for sea ice data
815	! -----------------------------------------------------------------------
816
817	CALL obs_coo_tim( icycle, &
818	& iyea0, imon0, iday0, ihou0, imin0, &
819	& seaicedata%nsurf, seaicedata%nyea, seaicedata%nmon, &
820	& seaicedata%nday, seaicedata%nhou, seaicedata%nmin, &
821	& seaicedata%nqc, seaicedata%mstp, iotdobs )
822	CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
823	! -----------------------------------------------------------------------
824	! Check for sea ice data failing the grid search
825	! -----------------------------------------------------------------------
826
827	CALL obs_coo_grd( seaicedata%nsurf, seaicedata%mi, seaicedata%mj, &
828	& seaicedata%nqc, igrdobs )
829	CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
830
831	! -----------------------------------------------------------------------
832	! Check for land points.
833	! -----------------------------------------------------------------------
834
835	CALL obs_coo_spc_2d( seaicedata%nsurf, &
836	& jpi, jpj, &
837	& seaicedata%mi, seaicedata%mj, &
838	& seaicedata%rlam, seaicedata%rphi, &
839	& glamt, gphit, &
840	& tmask(:,:,1), seaicedata%nqc, &
841	& iosdsobs, ilansobs, &
842	& inlasobs, ld_nea, &
843	& iqc_cutoff )
844
845	CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
846	CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
847	CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
848
849	! -----------------------------------------------------------------------
850	! Copy useful data from the seaicedata data structure to
851	! the seaicedatqc data structure
852	! -----------------------------------------------------------------------
853
854	! Allocate the selection arrays
855
856	ALLOCATE( llvalid(seaicedata%nsurf) )
857
858	! We want all data which has qc flags <= iqc_cutoff
859
860	llvalid(:) = ( seaicedata%nqc(:) <= iqc_cutoff )
861
862	! The actual copying
863
864	CALL obs_surf_compress( seaicedata, seaicedatqc, .TRUE., numout, &
865	& lvalid=llvalid )
866
867	! Dellocate the selection arrays
868	DEALLOCATE( llvalid )
869
870	! -----------------------------------------------------------------------
871	! Print information about what observations are left after qc
872	! -----------------------------------------------------------------------
873
874	! Update the total observation counter array
875
876	IF(lwp) THEN
877	WRITE(numout,*)
878	WRITE(numout,*) 'obs_pre_seaice :'
879	WRITE(numout,*) '~~~~~~~~~~~'
880	WRITE(numout,*)
881	WRITE(numout,*) ' Sea ice data outside time domain = ', &
882	& iotdobsmpp
883	WRITE(numout,*) ' Remaining sea ice data that failed grid search = ', &
884	& igrdobsmpp
885	WRITE(numout,*) ' Remaining sea ice data outside space domain = ', &
886	& iosdsobsmpp
887	WRITE(numout,*) ' Remaining sea ice data at land points = ', &
888	& ilansobsmpp
889	IF (ld_nea) THEN
890	WRITE(numout,*) ' Remaining sea ice data near land points (removed) = ', &
891	& inlasobsmpp
892	ELSE
893	WRITE(numout,*) ' Remaining sea ice data near land points (kept) = ', &
894	& inlasobsmpp
895	ENDIF
896	WRITE(numout,*) ' Sea ice data accepted = ', &
897	& seaicedatqc%nsurfmpp
898
899	WRITE(numout,*)
900	WRITE(numout,*) ' Number of observations per time step :'
901	WRITE(numout,*)
902	WRITE(numout,1997)
903	WRITE(numout,1998)
904	ENDIF
905
906	DO jobs = 1, seaicedatqc%nsurf
907	inrc = seaicedatqc%mstp(jobs) + 2 - nit000
908	seaicedatqc%nsstp(inrc) = seaicedatqc%nsstp(inrc) + 1
909	END DO
910
911	CALL obs_mpp_sum_integers( seaicedatqc%nsstp, seaicedatqc%nsstpmpp, &
912	& nitend - nit000 + 2 )
913
914	IF ( lwp ) THEN
915	DO jstp = nit000 - 1, nitend
916	inrc = jstp - nit000 + 2
917	WRITE(numout,1999) jstp, seaicedatqc%nsstpmpp(inrc)
918	END DO
919	ENDIF
920
921	1997 FORMAT(10X,'Time step',5X,'Sea ice data ')
922	1998 FORMAT(10X,'---------',5X,'-----------------')
923	1999 FORMAT(10X,I9,5X,I17)
924
925	END SUBROUTINE obs_pre_seaice
926
927	SUBROUTINE obs_pre_vel( profdata, prodatqc, ld_vel3d, ld_nea, ld_dailyav, kqc_cutoff )
928	!!----------------------------------------------------------------------
929	!! * ROUTINE obs_pre_taovel *
930	!!
931	!! ** Purpose : First level check and screening of U and V profiles
932	!!
933	!! ** Method : First level check and screening of U and V profiles
934	!!
935	!! History :
936	!! ! 2007-06 (K. Mogensen) original : T and S profile data
937	!! ! 2008-09 (M. Valdivieso) : TAO velocity data
938	!! ! 2009-01 (K. Mogensen) : New feedback strictuer
939	!!
940	!!----------------------------------------------------------------------
941	!! * Modules used
942	USE domstp ! Domain: set the time-step
943	USE par_oce ! Ocean parameters
944	USE dom_oce, ONLY : & ! Geographical information
945	& glamt, glamu, glamv, &
946	& gphit, gphiu, gphiv, &
947	& gdept_1d, &
948	& tmask, umask, vmask, &
949	& nproc
950	!! * Arguments
951	TYPE(obs_prof), INTENT(INOUT) :: profdata ! Full set of profile data
952	TYPE(obs_prof), INTENT(INOUT) :: prodatqc ! Subset of profile data not failing screening
953	LOGICAL, INTENT(IN) :: ld_vel3d ! Switch for zonal and meridional velocity components
954	LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land
955	LOGICAL, INTENT(IN) :: ld_dailyav ! Switch for daily average data
956	INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff ! cut off for QC value
957	!! * Local declarations
958	INTEGER :: iqc_cutoff = 255 ! cut off for QC value
959	INTEGER :: iyea0 ! Initial date
960	INTEGER :: imon0 ! - (year, month, day, hour, minute)
961	INTEGER :: iday0
962	INTEGER :: ihou0
963	INTEGER :: imin0
964	INTEGER :: icycle ! Current assimilation cycle
965	! Counters for observations that
966	INTEGER :: iotdobs ! - outside time domain
967	INTEGER :: iosduobs ! - outside space domain (zonal velocity component)
968	INTEGER :: iosdvobs ! - outside space domain (meridional velocity component)
969	INTEGER :: ilanuobs ! - within a model land cell (zonal velocity component)
970	INTEGER :: ilanvobs ! - within a model land cell (meridional velocity component)
971	INTEGER :: inlauobs ! - close to land (zonal velocity component)
972	INTEGER :: inlavobs ! - close to land (meridional velocity component)
973	INTEGER :: igrdobs ! - fail the grid search
974	INTEGER :: iuvchku ! - reject u if v rejected and vice versa
975	INTEGER :: iuvchkv !
976	! Global counters for observations that
977	INTEGER :: iotdobsmpp ! - outside time domain
978	INTEGER :: iosduobsmpp ! - outside space domain (zonal velocity component)
979	INTEGER :: iosdvobsmpp ! - outside space domain (meridional velocity component)
980	INTEGER :: ilanuobsmpp ! - within a model land cell (zonal velocity component)
981	INTEGER :: ilanvobsmpp ! - within a model land cell (meridional velocity component)
982	INTEGER :: inlauobsmpp ! - close to land (zonal velocity component)
983	INTEGER :: inlavobsmpp ! - close to land (meridional velocity component)
984	INTEGER :: igrdobsmpp ! - fail the grid search
985	INTEGER :: iuvchkumpp ! - reject u if v rejected and vice versa
986	INTEGER :: iuvchkvmpp !
987	TYPE(obs_prof_valid) :: llvalid ! Profile selection
988	TYPE(obs_prof_valid), DIMENSION(profdata%nvar) :: &
989	& llvvalid ! U,V selection
990	INTEGER :: jvar ! Variable loop variable
991	INTEGER :: jobs ! Obs. loop variable
992	INTEGER :: jstp ! Time loop variable
993	INTEGER :: inrc ! Time index variable
994
995	IF(lwp) WRITE(numout,*)'obs_pre_vel: Preparing the velocity profile data'
996
997	! Initial date initialization (year, month, day, hour, minute)
998	iyea0 = ndate0 / 10000
999	imon0 = ( ndate0 - iyea0 * 10000 ) / 100
1000	iday0 = ndate0 - iyea0 * 10000 - imon0 * 100
1001	ihou0 = 0
1002	imin0 = 0
1003
1004	icycle = no ! Assimilation cycle
1005
1006	! Diagnotics counters for various failures.
1007
1008	iotdobs = 0
1009	igrdobs = 0
1010	iosduobs = 0
1011	iosdvobs = 0
1012	ilanuobs = 0
1013	ilanvobs = 0
1014	inlauobs = 0
1015	inlavobs = 0
1016	iuvchku = 0
1017	iuvchkv = 0
1018
1019	! Set QC cutoff to optional value if provided
1020
1021	IF ( PRESENT(kqc_cutoff) ) iqc_cutoff=kqc_cutoff
1022
1023	! -----------------------------------------------------------------------
1024	! Find time coordinate for profiles
1025	! -----------------------------------------------------------------------
1026
1027	CALL obs_coo_tim_prof( icycle, &
1028	& iyea0, imon0, iday0, ihou0, imin0, &
1029	& profdata%nprof, profdata%nyea, profdata%nmon, &
1030	& profdata%nday, profdata%nhou, profdata%nmin, &
1031	& profdata%ntyp, profdata%nqc, profdata%mstp, &
1032	& iotdobs, ld_dailyav = ld_dailyav, &
1033	& kqc_cutoff=kqc_cutoff )
1034
1035	CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
1036
1037	! -----------------------------------------------------------------------
1038	! Check for profiles failing the grid search
1039	! -----------------------------------------------------------------------
1040
1041	CALL obs_coo_grd( profdata%nprof, profdata%mi(:,1), profdata%mj(:,1), &
1042	& profdata%nqc, igrdobs )
1043	CALL obs_coo_grd( profdata%nprof, profdata%mi(:,2), profdata%mj(:,2), &
1044	& profdata%nqc, igrdobs )
1045
1046	CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
1047
1048	! -----------------------------------------------------------------------
1049	! Reject all observations for profiles with nqc > 10
1050	! -----------------------------------------------------------------------
1051
1052	CALL obs_pro_rej( profdata, iqc_cutoff )
1053
1054	! -----------------------------------------------------------------------
1055	! Check for land points. This includes points below the model
1056	! bathymetry so this is done for every point in the profile
1057	! -----------------------------------------------------------------------
1058
1059	! Zonal Velocity Component
1060
1061	CALL obs_coo_spc_3d( profdata%nprof, profdata%nvprot(1), &
1062	& profdata%npvsta(:,1), profdata%npvend(:,1), &
1063	& jpi, jpj, &
1064	& jpk, &
1065	& profdata%mi, profdata%mj, &
1066	& profdata%var(1)%mvk, &
1067	& profdata%rlam, profdata%rphi, &
1068	& profdata%var(1)%vdep, &
1069	& glamu, gphiu, &
1070	& gdept_1d, umask, &
1071	& profdata%nqc, profdata%var(1)%nvqc, &
1072	& iosduobs, ilanuobs, &
1073	& inlauobs, ld_nea, &
1074	& iqc_cutoff )
1075
1076	CALL obs_mpp_sum_integer( iosduobs, iosduobsmpp )
1077	CALL obs_mpp_sum_integer( ilanuobs, ilanuobsmpp )
1078	CALL obs_mpp_sum_integer( inlauobs, inlauobsmpp )
1079
1080	! Meridional Velocity Component
1081
1082	CALL obs_coo_spc_3d( profdata%nprof, profdata%nvprot(2), &
1083	& profdata%npvsta(:,2), profdata%npvend(:,2), &
1084	& jpi, jpj, &
1085	& jpk, &
1086	& profdata%mi, profdata%mj, &
1087	& profdata%var(2)%mvk, &
1088	& profdata%rlam, profdata%rphi, &
1089	& profdata%var(2)%vdep, &
1090	& glamv, gphiv, &
1091	& gdept_1d, vmask, &
1092	& profdata%nqc, profdata%var(2)%nvqc, &
1093	& iosdvobs, ilanvobs, &
1094	& inlavobs, ld_nea, &
1095	& iqc_cutoff )
1096
1097	CALL obs_mpp_sum_integer( iosdvobs, iosdvobsmpp )
1098	CALL obs_mpp_sum_integer( ilanvobs, ilanvobsmpp )
1099	CALL obs_mpp_sum_integer( inlavobs, inlavobsmpp )
1100
1101	! -----------------------------------------------------------------------
1102	! Reject u if v is rejected and vice versa
1103	! -----------------------------------------------------------------------
1104
1105	CALL obs_uv_rej( profdata, iuvchku, iuvchkviqc_cutoff )
1106	CALL obs_mpp_sum_integer( iuvchku, iuvchkumpp )
1107	CALL obs_mpp_sum_integer( iuvchkv, iuvchkvmpp )
1108
1109	! -----------------------------------------------------------------------
1110	! Copy useful data from the profdata data structure to
1111	! the prodatqc data structure
1112	! -----------------------------------------------------------------------
1113
1114	! Allocate the selection arrays
1115
1116	ALLOCATE( llvalid%luse(profdata%nprof) )
1117	DO jvar = 1,profdata%nvar
1118	ALLOCATE( llvvalid(jvar)%luse(profdata%nvprot(jvar)) )
1119	END DO
1120
1121	! We want all data which has qc flags <= iqc_cutoff
1122
1123	llvalid%luse(:) = ( profdata%nqc(:) <= iqc_cutoff )
1124	DO jvar = 1,profdata%nvar
1125	llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= iqc_cutoff )
1126	END DO
1127
1128	! The actual copying
1129
1130	CALL obs_prof_compress( profdata, prodatqc, .TRUE., numout, &
1131	& lvalid=llvalid, lvvalid=llvvalid )
1132
1133	! Dellocate the selection arrays
1134	DEALLOCATE( llvalid%luse )
1135	DO jvar = 1,profdata%nvar
1136	DEALLOCATE( llvvalid(jvar)%luse )
1137	END DO
1138
1139	! -----------------------------------------------------------------------
1140	! Print information about what observations are left after qc
1141	! -----------------------------------------------------------------------
1142
1143	! Update the total observation counter array
1144
1145	IF(lwp) THEN
1146	WRITE(numout,*)
1147	WRITE(numout,*) 'obs_pre_vel :'
1148	WRITE(numout,*) '~~~~~~~~~~~'
1149	WRITE(numout,*)
1150	WRITE(numout,*) ' Profiles outside time domain = ', &
1151	& iotdobsmpp
1152	WRITE(numout,*) ' Remaining profiles that failed grid search = ', &
1153	& igrdobsmpp
1154	WRITE(numout,*) ' Remaining U data outside space domain = ', &
1155	& iosduobsmpp
1156	WRITE(numout,*) ' Remaining U data at land points = ', &
1157	& ilanuobsmpp
1158	IF (ld_nea) THEN
1159	WRITE(numout,*) ' Remaining U data near land points (removed) = ',&
1160	& inlauobsmpp
1161	ELSE
1162	WRITE(numout,*) ' Remaining U data near land points (kept) = ',&
1163	& inlauobsmpp
1164	ENDIF
1165	WRITE(numout,*) ' U observation rejected since V rejected = ', &
1166	& iuvchku
1167	WRITE(numout,*) ' U data accepted = ', &
1168	& prodatqc%nvprotmpp(1)
1169	WRITE(numout,*) ' Remaining V data outside space domain = ', &
1170	& iosdvobsmpp
1171	WRITE(numout,*) ' Remaining V data at land points = ', &
1172	& ilanvobsmpp
1173	IF (ld_nea) THEN
1174	WRITE(numout,*) ' Remaining V data near land points (removed) = ',&
1175	& inlavobsmpp
1176	ELSE
1177	WRITE(numout,*) ' Remaining V data near land points (kept) = ',&
1178	& inlavobsmpp
1179	ENDIF
1180	WRITE(numout,*) ' V observation rejected since U rejected = ', &
1181	& iuvchkv
1182	WRITE(numout,*) ' V data accepted = ', &
1183	& prodatqc%nvprotmpp(2)
1184
1185	WRITE(numout,*)
1186	WRITE(numout,*) ' Number of observations per time step :'
1187	WRITE(numout,*)
1188	WRITE(numout,997)
1189	WRITE(numout,998)
1190	ENDIF
1191
1192	DO jobs = 1, prodatqc%nprof
1193	inrc = prodatqc%mstp(jobs) + 2 - nit000
1194	prodatqc%npstp(inrc) = prodatqc%npstp(inrc) + 1
1195	DO jvar = 1, prodatqc%nvar
1196	IF ( prodatqc%npvend(jobs,jvar) > 0 ) THEN
1197	prodatqc%nvstp(inrc,jvar) = prodatqc%nvstp(inrc,jvar) + &
1198	& ( prodatqc%npvend(jobs,jvar) - &
1199	& prodatqc%npvsta(jobs,jvar) + 1 )
1200	ENDIF
1201	END DO
1202	END DO
1203
1204
1205	CALL obs_mpp_sum_integers( prodatqc%npstp, prodatqc%npstpmpp, &
1206	& nitend - nit000 + 2 )
1207	DO jvar = 1, prodatqc%nvar
1208	CALL obs_mpp_sum_integers( prodatqc%nvstp(:,jvar), &
1209	& prodatqc%nvstpmpp(:,jvar), &
1210	& nitend - nit000 + 2 )
1211	END DO
1212
1213	IF ( lwp ) THEN
1214	DO jstp = nit000 - 1, nitend
1215	inrc = jstp - nit000 + 2
1216	WRITE(numout,999) jstp, prodatqc%npstpmpp(inrc), &
1217	& prodatqc%nvstpmpp(inrc,1), &
1218	& prodatqc%nvstpmpp(inrc,2)
1219	END DO
1220	ENDIF
1221
1222	997 FORMAT(10X,'Time step',5X,'Profiles',5X,'Zonal Comp.',5X,'Meridional Comp.')
1223	998 FORMAT(10X,'---------',5X,'--------',5X,'-----------',5X,'----------------')
1224	999 FORMAT(10X,I9,5X,I8,5X,I11,5X,I8)
1225
1226	END SUBROUTINE obs_pre_vel
1227
1228	SUBROUTINE obs_coo_tim( kcycle, &
1229	& kyea0, kmon0, kday0, khou0, kmin0, &
1230	& kobsno, &
1231	& kobsyea, kobsmon, kobsday, kobshou, kobsmin, &
1232	& kobsqc, kobsstp, kotdobs )
1233	!!----------------------------------------------------------------------
1234	!! * ROUTINE obs_coo_tim *
1235	!!
1236	!! ** Purpose : Compute the number of time steps to the observation time.
1237	!!
1238	!! ** Method : For time coordinates ( yea_obs, mon_obs, day_obs,
1239	!! hou_obs, min_obs ), this routine locates the time step
1240	!! that is closest to this time.
1241	!!
1242	!! ** Action :
1243	!!
1244	!! References :
1245	!!
1246	!! History :
1247	!! ! 1997-07 (A. Weaver) Original
1248	!! ! 2006-08 (A. Weaver) NEMOVAR migration
1249	!! ! 2006-10 (A. Weaver) Cleanup
1250	!! ! 2007-01 (K. Mogensen) Rewritten with loop
1251	!! ! 2010-05 (D. Lea) Fix in leap year calculation for NEMO vn3.2
1252	!!----------------------------------------------------------------------
1253	!! * Modules used
1254	USE dom_oce, ONLY : & ! Geographical information
1255	& rdt
1256	USE phycst, ONLY : & ! Physical constants
1257	& rday, &
1258	& rmmss, &
1259	& rhhmm
1260	!! * Arguments
1261	INTEGER, INTENT(IN) :: kcycle ! Current cycle
1262	INTEGER, INTENT(IN) :: kyea0 ! Initial date coordinates
1263	INTEGER, INTENT(IN) :: kmon0
1264	INTEGER, INTENT(IN) :: kday0
1265	INTEGER, INTENT(IN) :: khou0
1266	INTEGER, INTENT(IN) :: kmin0
1267	INTEGER, INTENT(IN) :: kobsno ! Number of observations
1268	INTEGER, INTENT(INOUT) :: kotdobs ! Number of observations failing time check
1269	INTEGER, DIMENSION(kobsno), INTENT(IN ) :: &
1270	& kobsyea, & ! Observation time coordinates
1271	& kobsmon, &
1272	& kobsday, &
1273	& kobshou, &
1274	& kobsmin
1275	INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
1276	& kobsqc ! Quality control flag
1277	INTEGER, DIMENSION(kobsno), INTENT(OUT) :: &
1278	& kobsstp ! Number of time steps up to the
1279	! observation time
1280
1281	!! * Local declarations
1282	INTEGER :: jyea
1283	INTEGER :: jmon
1284	INTEGER :: jday
1285	INTEGER :: jobs
1286	INTEGER :: iyeastr
1287	INTEGER :: iyeaend
1288	INTEGER :: imonstr
1289	INTEGER :: imonend
1290	INTEGER :: idaystr
1291	INTEGER :: idayend
1292	INTEGER :: iskip
1293	INTEGER :: idaystp
1294	REAL(KIND=wp) :: zminstp
1295	REAL(KIND=wp) :: zhoustp
1296	REAL(KIND=wp) :: zobsstp
1297	INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year
1298
1299	!-----------------------------------------------------------------------
1300	! Initialization
1301	!-----------------------------------------------------------------------
1302
1303	! Intialize the number of time steps per day
1304	idaystp = NINT( rday / rdt )
1305
1306	!---------------------------------------------------------------------
1307	! Locate the model time coordinates for interpolation
1308	!---------------------------------------------------------------------
1309
1310	DO jobs = 1, kobsno
1311
1312	! Initialize the time step counter
1313	kobsstp(jobs) = nit000 - 1
1314
1315	! Flag if observation date is less than the initial date
1316
1317	IF ( ( kobsyea(jobs) < kyea0 ) &
1318	& .OR. ( ( kobsyea(jobs) == kyea0 ) &
1319	& .AND. ( kobsmon(jobs) < kmon0 ) ) &
1320	& .OR. ( ( kobsyea(jobs) == kyea0 ) &
1321	& .AND. ( kobsmon(jobs) == kmon0 ) &
1322	& .AND. ( kobsday(jobs) < kday0 ) ) &
1323	& .OR. ( ( kobsyea(jobs) == kyea0 ) &
1324	& .AND. ( kobsmon(jobs) == kmon0 ) &
1325	& .AND. ( kobsday(jobs) == kday0 ) &
1326	& .AND. ( kobshou(jobs) < khou0 ) ) &
1327	& .OR. ( ( kobsyea(jobs) == kyea0 ) &
1328	& .AND. ( kobsmon(jobs) == kmon0 ) &
1329	& .AND. ( kobsday(jobs) == kday0 ) &
1330	& .AND. ( kobshou(jobs) == khou0 ) &
1331	& .AND. ( kobsmin(jobs) <= kmin0 ) ) ) THEN
1332	kobsstp(jobs) = -1
1333	kobsqc(jobs) = kobsqc(jobs) + 11
1334	kotdobs = kotdobs + 1
1335	CYCLE
1336	ENDIF
1337
1338	! Compute the number of time steps to the observation day
1339	iyeastr = kyea0
1340	iyeaend = kobsyea(jobs)
1341
1342	!---------------------------------------------------------------------
1343	! Year loop
1344	!---------------------------------------------------------------------
1345	DO jyea = iyeastr, iyeaend
1346
1347	CALL calc_month_len( jyea, imonth_len )
1348
1349	imonstr = 1
1350	IF ( jyea == kyea0 ) imonstr = kmon0
1351	imonend = 12
1352	IF ( jyea == kobsyea(jobs) ) imonend = kobsmon(jobs)
1353
1354	! Month loop
1355	DO jmon = imonstr, imonend
1356
1357	idaystr = 1
1358	IF ( ( jmon == kmon0 ) &
1359	& .AND. ( jyea == kyea0 ) ) idaystr = kday0
1360	idayend = imonth_len(jmon)
1361	IF ( ( jmon == kobsmon(jobs) ) &
1362	& .AND. ( jyea == kobsyea(jobs) ) ) idayend = kobsday(jobs) - 1
1363
1364	! Day loop
1365	DO jday = idaystr, idayend
1366	kobsstp(jobs) = kobsstp(jobs) + idaystp
1367	END DO
1368
1369	END DO
1370
1371	END DO
1372
1373	! Add in the number of time steps to the observation minute
1374	zminstp = rmmss / rdt
1375	zhoustp = rhhmm * zminstp
1376
1377	zobsstp = REAL( kobsmin(jobs) - kmin0, KIND=wp ) * zminstp &
1378	& + REAL( kobshou(jobs) - khou0, KIND=wp ) * zhoustp
1379	kobsstp(jobs) = kobsstp(jobs) + NINT( zobsstp )
1380
1381	! Flag if observation step outside the time window
1382	IF ( ( kobsstp(jobs) < ( nit000 - 1 ) ) &
1383	& .OR.( kobsstp(jobs) > nitend ) ) THEN
1384	kobsqc(jobs) = kobsqc(jobs) + 12
1385	kotdobs = kotdobs + 1
1386	CYCLE
1387	ENDIF
1388
1389	END DO
1390
1391	END SUBROUTINE obs_coo_tim
1392
1393	SUBROUTINE calc_month_len( iyear, imonth_len )
1394	!!----------------------------------------------------------------------
1395	!! * ROUTINE calc_month_len *
1396	!!
1397	!! ** Purpose : Compute the number of days in a months given a year.
1398	!!
1399	!! ** Method :
1400	!!
1401	!! ** Action :
1402	!!
1403	!! History :
1404	!! ! 10-05 (D. Lea) New routine based on day_init
1405	!!----------------------------------------------------------------------
1406
1407	INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year
1408	INTEGER :: iyear !: year
1409
1410	! length of the month of the current year (from nleapy, read in namelist)
1411	IF ( nleapy < 2 ) THEN
1412	imonth_len(:) = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)
1413	IF ( nleapy == 1 ) THEN ! we are using calendar with leap years
1414	IF ( MOD(iyear, 4) == 0 .AND. ( MOD(iyear, 400) == 0 .OR. MOD(iyear, 100) /= 0 ) ) THEN
1415	imonth_len(2) = 29
1416	ENDIF
1417	ENDIF
1418	ELSE
1419	imonth_len(:) = nleapy ! all months with nleapy days per year
1420	ENDIF
1421
1422	END SUBROUTINE
1423
1424	SUBROUTINE obs_coo_tim_prof( kcycle, &
1425	& kyea0, kmon0, kday0, khou0, kmin0, &
1426	& kobsno, &
1427	& kobsyea, kobsmon, kobsday, kobshou, kobsmin, &
1428	& ktyp, kobsqc, kobsstp, kotdobs, kdailyavtypes, &
1429	& ld_dailyav, kqc_cutoff )
1430	!!----------------------------------------------------------------------
1431	!! * ROUTINE obs_coo_tim *
1432	!!
1433	!! ** Purpose : Compute the number of time steps to the observation time.
1434	!!
1435	!! ** Method : For time coordinates ( yea_obs, mon_obs, day_obs,
1436	!! hou_obs, min_obs ), this routine locates the time step
1437	!! that is closest to this time.
1438	!!
1439	!! ** Action :
1440	!!
1441	!! References :
1442	!!
1443	!! History :
1444	!! ! 1997-07 (A. Weaver) Original
1445	!! ! 2006-08 (A. Weaver) NEMOVAR migration
1446	!! ! 2006-10 (A. Weaver) Cleanup
1447	!! ! 2007-01 (K. Mogensen) Rewritten with loop
1448	!!----------------------------------------------------------------------
1449	!! * Modules used
1450	!! * Arguments
1451	INTEGER, INTENT(IN) :: kcycle ! Current cycle
1452	INTEGER, INTENT(IN) :: kyea0 ! Initial date coordinates
1453	INTEGER, INTENT(IN) :: kmon0
1454	INTEGER, INTENT(IN) :: kday0
1455	INTEGER, INTENT(IN) :: khou0
1456	INTEGER, INTENT(IN) :: kmin0
1457	INTEGER, INTENT(IN) :: kobsno ! Number of observations
1458	INTEGER, INTENT(INOUT) :: kotdobs ! Number of observations failing time check
1459	INTEGER, DIMENSION(kobsno), INTENT(IN ) :: &
1460	& kobsyea, & ! Observation time coordinates
1461	& kobsmon, &
1462	& kobsday, &
1463	& kobshou, &
1464	& kobsmin, &
1465	& ktyp ! Observation type.
1466	INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
1467	& kobsqc ! Quality control flag
1468	INTEGER, DIMENSION(kobsno), INTENT(OUT) :: &
1469	& kobsstp ! Number of time steps up to the
1470	! observation time
1471	INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
1472	& kdailyavtypes ! Types for daily averages
1473	LOGICAL, OPTIONAL :: ld_dailyav ! All types are daily averages
1474	INTEGER, OPTIONAL, INTENT(IN) :: kqc_cutoff ! QC cutoff value
1475	!! * Local declarations
1476	INTEGER :: jobs
1477	INTEGER :: iqc_cutoff=255
1478
1479	!-----------------------------------------------------------------------
1480	! Call standard obs_coo_tim
1481	!-----------------------------------------------------------------------
1482
1483	CALL obs_coo_tim( kcycle, &
1484	& kyea0, kmon0, kday0, khou0, kmin0, &
1485	& kobsno, &
1486	& kobsyea, kobsmon, kobsday, kobshou, kobsmin, &
1487	& kobsqc, kobsstp, kotdobs )
1488
1489	!------------------------------------------------------------------------
1490	! Always reject daily averaged data (e.g. MRB data (820)) at initial time
1491	!------------------------------------------------------------------------
1492
1493	IF ( PRESENT(kdailyavtypes) ) THEN
1494	DO jobs = 1, kobsno
1495
1496	IF ( kobsqc(jobs) <= iqc_cutoff ) THEN
1497
1498	IF ( ( kobsstp(jobs) == (nit000 - 1) ).AND.&
1499	& ( ANY (kdailyavtypes(:) == ktyp(jobs)) ) ) THEN
1500	kobsqc(jobs) = IBSET(kobsqc(jobs),13)
1501	kotdobs = kotdobs + 1
1502	CYCLE
1503	ENDIF
1504
1505	ENDIF
1506	END DO
1507	ENDIF
1508
1509	!------------------------------------------------------------------------
1510	! If ld_dailyav is set then all data assumed to be daily averaged
1511	!------------------------------------------------------------------------
1512
1513	IF ( PRESENT( ld_dailyav) ) THEN
1514	IF (ld_dailyav) THEN
1515	DO jobs = 1, kobsno
1516
1517	IF ( kobsqc(jobs) <= iqc_cutoff ) THEN
1518
1519	IF ( kobsstp(jobs) == (nit000 - 1) ) THEN
1520	kobsqc(jobs) = IBSET(kobsqc(jobs),13)
1521	kotdobs = kotdobs + 1
1522	CYCLE
1523	ENDIF
1524
1525	ENDIF
1526	END DO
1527	ENDIF
1528	ENDIF
1529
1530	END SUBROUTINE obs_coo_tim_prof
1531
1532	SUBROUTINE obs_coo_grd( kobsno, kobsi, kobsj, kobsqc, kgrdobs )
1533	!!----------------------------------------------------------------------
1534	!! * ROUTINE obs_coo_grd *
1535	!!
1536	!! ** Purpose : Verify that the grid search has not failed
1537	!!
1538	!! ** Method : The previously computed i,j indeces are checked
1539	!!
1540	!! ** Action :
1541	!!
1542	!! References :
1543	!!
1544	!! History :
1545	!! ! 2007-01 (K. Mogensen) Original
1546	!!----------------------------------------------------------------------
1547	!! * Arguments
1548	INTEGER, INTENT(IN) :: kobsno ! Number of observations
1549	INTEGER, DIMENSION(kobsno), INTENT(IN ) :: &
1550	& kobsi, & ! i,j indeces previously computed
1551	& kobsj
1552	INTEGER, INTENT(INOUT) :: kgrdobs ! Number of observations failing the check
1553	INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
1554	& kobsqc ! Quality control flag
1555
1556	!! * Local declarations
1557	INTEGER :: jobs ! Loop variable
1558
1559	! Flag if the grid search failed
1560
1561	DO jobs = 1, kobsno
1562	IF ( ( kobsi(jobs) <= 0 ) .AND. ( kobsj(jobs) <= 0 ) ) THEN
1563	kobsqc(jobs) = IBSET(kobsqc(jobs),12)
1564	kgrdobs = kgrdobs + 1
1565	ENDIF
1566	END DO
1567
1568	END SUBROUTINE obs_coo_grd
1569
1570	SUBROUTINE obs_coo_spc_2d( kobsno, kpi, kpj, &
1571	& kobsi, kobsj, pobslam, pobsphi, &
1572	& plam, pphi, pmask, &
1573	& kobsqc, kosdobs, klanobs, &
1574	& knlaobs,ld_nea, kqc_cutoff )
1575	!!----------------------------------------------------------------------
1576	!! * ROUTINE obs_coo_spc_2d *
1577	!!
1578	!! ** Purpose : Check for points outside the domain and land points
1579	!!
1580	!! ** Method : Remove the observations that are outside the model space
1581	!! and time domain or located within model land cells.
1582	!!
1583	!! ** Action :
1584	!!
1585	!! History :
1586	!! ! 2007-03 (A. Weaver, K. Mogensen) Original
1587	!! ! 2007-06 (K. Mogensen et al) Reject obs. near land.
1588	!!----------------------------------------------------------------------
1589	!! * Modules used
1590
1591	!! * Arguments
1592	INTEGER, INTENT(IN) :: kobsno ! Total number of observations
1593	INTEGER, INTENT(IN) :: kpi ! Number of grid points in (i,j)
1594	INTEGER, INTENT(IN) :: kpj
1595	INTEGER, DIMENSION(kobsno), INTENT(IN) :: &
1596	& kobsi, & ! Observation (i,j) coordinates
1597	& kobsj
1598	REAL(KIND=wp), DIMENSION(kobsno), INTENT(IN) :: &
1599	& pobslam, & ! Observation (lon,lat) coordinates
1600	& pobsphi
1601	REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) :: &
1602	& plam, pphi ! Model (lon,lat) coordinates
1603	REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) :: &
1604	& pmask ! Land mask array
1605	INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
1606	& kobsqc ! Observation quality control
1607	INTEGER, INTENT(INOUT) :: kosdobs ! Observations outside space domain
1608	INTEGER, INTENT(INOUT) :: klanobs ! Observations within a model land cell
1609	INTEGER, INTENT(INOUT) :: knlaobs ! Observations near land
1610	LOGICAL, INTENT(IN) :: ld_nea ! Flag observations near land
1611	INTEGER, INTENT(IN) :: kqc_cutoff ! Cutoff QC value
1612	!! * Local declarations
1613	REAL(KIND=wp), DIMENSION(2,2,kobsno) :: &
1614	& zgmsk ! Grid mask
1615	REAL(KIND=wp), DIMENSION(2,2,kobsno) :: &
1616	& zglam, & ! Model longitude at grid points
1617	& zgphi ! Model latitude at grid points
1618	INTEGER, DIMENSION(2,2,kobsno) :: &
1619	& igrdi, & ! Grid i,j
1620	& igrdj
1621	LOGICAL :: lgridobs ! Is observation on a model grid point.
1622	INTEGER :: iig, ijg ! i,j of observation on model grid point.
1623	INTEGER :: jobs, ji, jj
1624
1625	! Get grid point indices
1626
1627	DO jobs = 1, kobsno
1628
1629	! For invalid points use 2,2
1630
1631	IF ( kobsqc(jobs) >= kqc_cutoff ) THEN
1632
1633	igrdi(1,1,jobs) = 1
1634	igrdj(1,1,jobs) = 1
1635	igrdi(1,2,jobs) = 1
1636	igrdj(1,2,jobs) = 2
1637	igrdi(2,1,jobs) = 2
1638	igrdj(2,1,jobs) = 1
1639	igrdi(2,2,jobs) = 2
1640	igrdj(2,2,jobs) = 2
1641
1642	ELSE
1643
1644	igrdi(1,1,jobs) = kobsi(jobs)-1
1645	igrdj(1,1,jobs) = kobsj(jobs)-1
1646	igrdi(1,2,jobs) = kobsi(jobs)-1
1647	igrdj(1,2,jobs) = kobsj(jobs)
1648	igrdi(2,1,jobs) = kobsi(jobs)
1649	igrdj(2,1,jobs) = kobsj(jobs)-1
1650	igrdi(2,2,jobs) = kobsi(jobs)
1651	igrdj(2,2,jobs) = kobsj(jobs)
1652
1653	ENDIF
1654
1655	END DO
1656
1657	CALL obs_int_comm_2d( 2, 2, kobsno, igrdi, igrdj, pmask, zgmsk )
1658	CALL obs_int_comm_2d( 2, 2, kobsno, igrdi, igrdj, plam, zglam )
1659	CALL obs_int_comm_2d( 2, 2, kobsno, igrdi, igrdj, pphi, zgphi )
1660
1661	DO jobs = 1, kobsno
1662
1663	! Skip bad observations
1664	IF ( kobsqc(jobs) >= kqc_cutoff ) CYCLE
1665
1666	! Flag if the observation falls outside the model spatial domain
1667	IF ( ( pobslam(jobs) < -180. ) &
1668	& .OR. ( pobslam(jobs) > 180. ) &
1669	& .OR. ( pobsphi(jobs) < -90. ) &
1670	& .OR. ( pobsphi(jobs) > 90. ) ) THEN
1671	kobsqc(jobs) = IBSET(kobsqc(jobs),11)
1672	kosdobs = kosdobs + 1
1673	CYCLE
1674	ENDIF
1675
1676	! Flag if the observation falls with a model land cell
1677	IF ( SUM( zgmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN
1678	kobsqc(jobs) = IBSET(kobsqc(jobs),10)
1679	klanobs = klanobs + 1
1680	CYCLE
1681	ENDIF
1682
1683	! Check if this observation is on a grid point
1684
1685	lgridobs = .FALSE.
1686	iig = -1
1687	ijg = -1
1688	DO jj = 1, 2
1689	DO ji = 1, 2
1690	IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) &
1691	& .AND. &
1692	& ( ABS( zglam(ji,jj,jobs) - pobslam(jobs) ) < 1.0e-6_wp ) &
1693	& ) THEN
1694	lgridobs = .TRUE.
1695	iig = ji
1696	ijg = jj
1697	ENDIF
1698	END DO
1699	END DO
1700
1701	! For observations on the grid reject them if their are at
1702	! a masked point
1703
1704	IF (lgridobs) THEN
1705	IF (zgmsk(iig,ijg,jobs) == 0.0_wp ) THEN
1706	kobsqc(jobs) = IBSET(kobsqc(jobs),10)
1707	klanobs = klanobs + 1
1708	CYCLE
1709	ENDIF
1710	ENDIF
1711
1712	! Flag if the observation falls is close to land
1713	IF ( MINVAL( zgmsk(1:2,1:2,jobs) ) == 0.0_wp) THEN
1714	IF (ld_nea) kobsqc(jobs) = IBSET(kobsqc(jobs),9)
1715	knlaobs = knlaobs + 1
1716	CYCLE
1717	ENDIF
1718
1719	END DO
1720
1721	END SUBROUTINE obs_coo_spc_2d
1722
1723	SUBROUTINE obs_coo_spc_3d( kprofno, kobsno, kpstart, kpend, &
1724	& kpi, kpj, kpk, &
1725	& kobsi, kobsj, kobsk, &
1726	& pobslam, pobsphi, pobsdep, &
1727	& plam, pphi, pdep, pmask, &
1728	& kpobsqc, kobsqc, kosdobs, &
1729	& klanobs, knlaobs, ld_nea, &
1730	& kqc_cutoff )
1731	!!----------------------------------------------------------------------
1732	!! * ROUTINE obs_coo_spc_3d *
1733	!!
1734	!! ** Purpose : Check for points outside the domain and land points
1735	!! Reset depth of observation above highest model level
1736	!! to the value of highest model level
1737	!!
1738	!! ** Method : Remove the observations that are outside the model space
1739	!! and time domain or located within model land cells.
1740	!!
1741	!! NB. T and S profile observations lying between the ocean
1742	!! surface and the depth of the first model T point are
1743	!! assigned a depth equal to that of the first model T pt.
1744	!!
1745	!! ** Action :
1746	!!
1747	!! History :
1748	!! ! 2007-01 (K. Mogensen) Rewrite of parts of obs_scr
1749	!! ! 2007-06 (K. Mogensen et al) Reject obs. near land.
1750	!!----------------------------------------------------------------------
1751	!! * Modules used
1752	USE dom_oce, ONLY : & ! Geographical information
1753	& gdepw_1d
1754
1755	!! * Arguments
1756	INTEGER, INTENT(IN) :: kprofno ! Number of profiles
1757	INTEGER, INTENT(IN) :: kobsno ! Total number of observations
1758	INTEGER, INTENT(IN) :: kpi ! Number of grid points in (i,j,k)
1759	INTEGER, INTENT(IN) :: kpj
1760	INTEGER, INTENT(IN) :: kpk
1761	INTEGER, DIMENSION(kprofno), INTENT(IN) :: &
1762	& kpstart, & ! Start of individual profiles
1763	& kpend ! End of individual profiles
1764	INTEGER, DIMENSION(kprofno), INTENT(IN) :: &
1765	& kobsi, & ! Observation (i,j) coordinates
1766	& kobsj
1767	INTEGER, DIMENSION(kobsno), INTENT(IN) :: &
1768	& kobsk ! Observation k coordinate
1769	REAL(KIND=wp), DIMENSION(kprofno), INTENT(IN) :: &
1770	& pobslam, & ! Observation (lon,lat) coordinates
1771	& pobsphi
1772	REAL(KIND=wp), DIMENSION(kobsno), INTENT(INOUT) :: &
1773	& pobsdep ! Observation depths
1774	REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) :: &
1775	& plam, pphi ! Model (lon,lat) coordinates
1776	REAL(KIND=wp), DIMENSION(kpk), INTENT(IN) :: &
1777	& pdep ! Model depth coordinates
1778	REAL(KIND=wp), DIMENSION(kpi,kpj,kpk), INTENT(IN) :: &
1779	& pmask ! Land mask array
1780	INTEGER, DIMENSION(kprofno), INTENT(INOUT) :: &
1781	& kpobsqc ! Profile quality control
1782	INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
1783	& kobsqc ! Observation quality control
1784	INTEGER, INTENT(INOUT) :: kosdobs ! Observations outside space domain
1785	INTEGER, INTENT(INOUT) :: klanobs ! Observations within a model land cell
1786	INTEGER, INTENT(INOUT) :: knlaobs ! Observations near land
1787	LOGICAL, INTENT(IN) :: ld_nea ! Flag observations near land
1788	INTEGER, INTENT(IN) :: kqc_cutoff ! Cutoff QC value
1789	!! * Local declarations
1790	REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: &
1791	& zgmsk ! Grid mask
1792	REAL(KIND=wp), DIMENSION(2,2,kprofno) :: &
1793	& zglam, & ! Model longitude at grid points
1794	& zgphi ! Model latitude at grid points
1795	INTEGER, DIMENSION(2,2,kprofno) :: &
1796	& igrdi, & ! Grid i,j
1797	& igrdj
1798	LOGICAL :: lgridobs ! Is observation on a model grid point.
1799	INTEGER :: iig, ijg ! i,j of observation on model grid point.
1800	INTEGER :: jobs, jobsp, jk, ji, jj
1801
1802	! Get grid point indices
1803
1804	DO jobs = 1, kprofno
1805
1806	! For invalid points use 2,2
1807
1808	IF ( kpobsqc(jobs) >= kqc_cutoff ) THEN
1809
1810	igrdi(1,1,jobs) = 1
1811	igrdj(1,1,jobs) = 1
1812	igrdi(1,2,jobs) = 1
1813	igrdj(1,2,jobs) = 2
1814	igrdi(2,1,jobs) = 2
1815	igrdj(2,1,jobs) = 1
1816	igrdi(2,2,jobs) = 2
1817	igrdj(2,2,jobs) = 2
1818
1819	ELSE
1820
1821	igrdi(1,1,jobs) = kobsi(jobs)-1
1822	igrdj(1,1,jobs) = kobsj(jobs)-1
1823	igrdi(1,2,jobs) = kobsi(jobs)-1
1824	igrdj(1,2,jobs) = kobsj(jobs)
1825	igrdi(2,1,jobs) = kobsi(jobs)
1826	igrdj(2,1,jobs) = kobsj(jobs)-1
1827	igrdi(2,2,jobs) = kobsi(jobs)
1828	igrdj(2,2,jobs) = kobsj(jobs)
1829
1830	ENDIF
1831
1832	END DO
1833
1834	CALL obs_int_comm_3d( 2, 2, kprofno, kpk, igrdi, igrdj, pmask, zgmsk )
1835	CALL obs_int_comm_2d( 2, 2, kprofno, igrdi, igrdj, plam, zglam )
1836	CALL obs_int_comm_2d( 2, 2, kprofno, igrdi, igrdj, pphi, zgphi )
1837
1838	DO jobs = 1, kprofno
1839
1840	! Skip bad profiles
1841	IF ( kpobsqc(jobs) >= kqc_cutoff ) CYCLE
1842
1843	! Check if this observation is on a grid point
1844
1845	lgridobs = .FALSE.
1846	iig = -1
1847	ijg = -1
1848	DO jj = 1, 2
1849	DO ji = 1, 2
1850	IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) &
1851	& .AND. &
1852	& ( ABS( zglam(ji,jj,jobs) - pobslam(jobs) ) < 1.0e-6_wp ) &
1853	& ) THEN
1854	lgridobs = .TRUE.
1855	iig = ji
1856	ijg = jj
1857	ENDIF
1858	END DO
1859	END DO
1860
1861	! Reject observations
1862
1863	DO jobsp = kpstart(jobs), kpend(jobs)
1864
1865	! Flag if the observation falls outside the model spatial domain
1866	IF ( ( pobslam(jobs) < -180. ) &
1867	& .OR. ( pobslam(jobs) > 180. ) &
1868	& .OR. ( pobsphi(jobs) < -90. ) &
1869	& .OR. ( pobsphi(jobs) > 90. ) &
1870	& .OR. ( pobsdep(jobsp) < 0.0 ) &
1871	& .OR. ( pobsdep(jobsp) > gdepw_1d(kpk)) ) THEN
1872	kobsqc(jobsp) = IBSET(kobsqc(jobsp),11)
1873	kosdobs = kosdobs + 1
1874	CYCLE
1875	ENDIF
1876
1877	! Flag if the observation falls with a model land cell
1878	IF ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) &
1879	& == 0.0_wp ) THEN
1880	kobsqc(jobsp) = IBSET(kobsqc(jobsp),10)
1881	klanobs = klanobs + 1
1882	CYCLE
1883	ENDIF
1884
1885	! For observations on the grid reject them if their are at
1886	! a masked point
1887
1888	IF (lgridobs) THEN
1889	IF (zgmsk(iig,ijg,kobsk(jobsp)-1,jobs) == 0.0_wp ) THEN
1890	kobsqc(jobsp) = IBSET(kobsqc(jobsp),9)
1891	klanobs = klanobs + 1
1892	CYCLE
1893	ENDIF
1894	ENDIF
1895
1896	! Flag if the observation falls is close to land
1897	IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == &
1898	& 0.0_wp) THEN
1899	IF (ld_nea) kobsqc(jobsp) = IBSET(kobsqc(jobsp),10)
1900	knlaobs = knlaobs + 1
1901	ENDIF
1902
1903	! Set observation depth equal to that of the first model depth
1904	IF ( pobsdep(jobsp) <= pdep(1) ) THEN
1905	pobsdep(jobsp) = pdep(1)
1906	ENDIF
1907
1908	END DO
1909	END DO
1910
1911	END SUBROUTINE obs_coo_spc_3d
1912
1913	SUBROUTINE obs_pro_rej( profdata, kqc_cutoff )
1914	!!----------------------------------------------------------------------
1915	!! * ROUTINE obs_pro_rej *
1916	!!
1917	!! ** Purpose : Reject all data within a rejected profile
1918	!!
1919	!! ** Method :
1920	!!
1921	!! ** Action :
1922	!!
1923	!! References :
1924	!!
1925	!! History :
1926	!! ! 2007-10 (K. Mogensen) Original code
1927	!!----------------------------------------------------------------------
1928	!! * Modules used
1929	!! * Arguments
1930	TYPE(obs_prof), INTENT(INOUT) :: profdata ! Profile data
1931	INTEGER, INTENT(IN) :: kqc_cutoff ! QC cutoff value
1932	!! * Local declarations
1933	INTEGER :: jprof
1934	INTEGER :: jvar
1935	INTEGER :: jobs
1936
1937	! Loop over profiles
1938
1939	DO jprof = 1, profdata%nprof
1940
1941	IF ( profdata%nqc(jprof) > kqc_cutoff ) THEN
1942
1943	DO jvar = 1, profdata%nvar
1944
1945	DO jobs = profdata%npvsta(jprof,jvar), &
1946	& profdata%npvend(jprof,jvar)
1947
1948	profdata%var(jvar)%nvqc(jobs) = &
1949	& IBSET(profdata%var(jvar)%nvqc(jobs),14)
1950
1951	END DO
1952
1953	END DO
1954
1955	ENDIF
1956
1957	END DO
1958
1959	END SUBROUTINE obs_pro_rej
1960
1961	SUBROUTINE obs_uv_rej( profdata, knumu, knumv, kqc_cutoff )
1962	!!----------------------------------------------------------------------
1963	!! * ROUTINE obs_uv_rej *
1964	!!
1965	!! ** Purpose : Reject u if v is rejected and vice versa
1966	!!
1967	!! ** Method :
1968	!!
1969	!! ** Action :
1970	!!
1971	!! References :
1972	!!
1973	!! History :
1974	!! ! 2009-2 (K. Mogensen) Original code
1975	!!----------------------------------------------------------------------
1976	!! * Modules used
1977	!! * Arguments
1978	TYPE(obs_prof), INTENT(INOUT) :: profdata ! Profile data
1979	INTEGER, INTENT(INOUT) :: knumu ! Number of u rejected
1980	INTEGER, INTENT(INOUT) :: knumv ! Number of v rejected
1981	INTEGER, INTENT(IN) :: kqc_cutoff ! QC cutoff value
1982	!! * Local declarations
1983	INTEGER :: jprof
1984	INTEGER :: jvar
1985	INTEGER :: jobs
1986
1987	! Loop over profiles
1988
1989	DO jprof = 1, profdata%nprof
1990
1991	IF ( ( profdata%npvsta(jprof,1) /= profdata%npvsta(jprof,2) ) .OR. &
1992	& ( profdata%npvend(jprof,1) /= profdata%npvend(jprof,2) ) ) THEN
1993
1994	CALL ctl_stop('U,V profiles inconsistent in obs_uv_rej')
1995	RETURN
1996
1997	ENDIF
1998
1999	DO jobs = profdata%npvsta(jprof,1), profdata%npvend(jprof,1)
2000
2001	IF ( ( profdata%var(1)%nvqc(jobs) > kqc_cutoff ) .AND. &
2002	& ( profdata%var(2)%nvqc(jobs) <= kqc_cutoff) ) THEN
2003	profdata%var(2)%nvqc(jobs) = profdata%var(2)%nvqc(jobs) + 42
2004	knumv = knumv + 1
2005	ENDIF
2006	IF ( ( profdata%var(2)%nvqc(jobs) > kqc_cutoff ) .AND. &
2007	& ( profdata%var(1)%nvqc(jobs) <= kqc_cutoff) ) THEN
2008	profdata%var(1)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15)
2009	knumu = knumu + 1
2010	ENDIF
2011
2012	END DO
2013
2014	END DO
2015
2016	END SUBROUTINE obs_uv_rej
2017
2018	END MODULE obs_prep

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

Download in other formats: