source: branches/UKMO/dev_r5518_obs_oper_update_sit/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 11932

Last change on this file since 11932 was 11932, checked in by dcarneir, 2 years ago

Changes in OBS and SBC routines for sea ice thickness data assimilation

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