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/trunk/src/OCE/OBS – NEMO

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

Last change on this file was 15062, checked in by jchanut, 3 years ago

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