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

source: branches/UKMO/dev_CO6_obs_bound_reject/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 6980

Last change on this file since 6980 was 6980, checked in by kingr, 8 years ago

Update to add ability to reject observations near open boundaries.

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