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 NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 11361

Last change on this file since 11361 was 11361, checked in by jcastill, 5 years ago

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