source: NEMO/trunk/src/OCE/OBS/obs_prep.F90

Last change on this file was 14275, checked in by smasson, 3 months ago

trunk: suppress nproc ( = mpprank = narea-1)

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