source: branches/UKMO/dev_r5518_obs_oper_update_NearSurfDepths/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 11451

Last change on this file since 11451 was 11451, checked in by kingr, 14 months ago

Modified sheck on min obs deth to use 3D depth array.

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