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

Context Navigation

source: branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 6073

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

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

Download in other formats: