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/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/OBS – NEMO

source: NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/OBS/obs_prep.F90 @ 14652

Last change on this file since 14652 was 14652, checked in by sparonuz, 3 years ago

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