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 @ 14219

Last change on this file since 14219 was 14219, checked in by mcastril, 4 years ago

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