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

source: branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 6990

Last change on this file since 6990 was 6990, checked in by kingr, 7 years ago

Added code from nemo3.4 OBS branch to allow rejection of observations near open boundaries.

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