source: branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 11202

Last change on this file since 11202 was 11202, checked in by jcastill, 15 months ago

Copy of branch branches/UKMO/dev_r5518_obs_oper_update@11130 without namelist_ref changes to allow merging with coupling and biogeochemistry branches

  • Property svn:keywords set to Id
File size: 58.9 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         & zgdepw         
1162      REAL(KIND=wp), DIMENSION(2,2,kprofno) :: &
1163         & zglam, &           ! Model longitude at grid points
1164         & zgphi              ! Model latitude at grid points
1165      INTEGER, DIMENSION(2,2,kprofno) :: &
1166         & igrdi, &           ! Grid i,j
1167         & igrdj
1168      LOGICAL :: lgridobs           ! Is observation on a model grid point.
1169      LOGICAL :: ll_next_to_land    ! Is a profile next to land
1170      INTEGER :: iig, ijg           ! i,j of observation on model grid point.
1171      INTEGER :: jobs, jobsp, jk, ji, jj
1172
1173      ! Get grid point indices
1174     
1175      DO jobs = 1, kprofno
1176
1177         ! For invalid points use 2,2
1178
1179         IF ( kpobsqc(jobs) >= kqc_cutoff ) THEN
1180           
1181            igrdi(1,1,jobs) = 1
1182            igrdj(1,1,jobs) = 1
1183            igrdi(1,2,jobs) = 1
1184            igrdj(1,2,jobs) = 2
1185            igrdi(2,1,jobs) = 2
1186            igrdj(2,1,jobs) = 1
1187            igrdi(2,2,jobs) = 2
1188            igrdj(2,2,jobs) = 2
1189           
1190         ELSE
1191           
1192            igrdi(1,1,jobs) = kobsi(jobs)-1
1193            igrdj(1,1,jobs) = kobsj(jobs)-1
1194            igrdi(1,2,jobs) = kobsi(jobs)-1
1195            igrdj(1,2,jobs) = kobsj(jobs)
1196            igrdi(2,1,jobs) = kobsi(jobs)
1197            igrdj(2,1,jobs) = kobsj(jobs)-1
1198            igrdi(2,2,jobs) = kobsi(jobs)
1199            igrdj(2,2,jobs) = kobsj(jobs)
1200           
1201         ENDIF
1202         
1203      END DO
1204
1205#if defined key_bdy 
1206      ! Create a mask grid points in boundary rim
1207      IF (ld_bound_reject) THEN           
1208         zbdymask(:,:) = 1.0_wp
1209         DO ji = 1, nb_bdy
1210            DO jj = 1, idx_bdy(ji)%nblen(1)
1211               zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp
1212            ENDDO
1213         ENDDO
1214      ENDIF
1215 
1216      CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, zbdymask, zbmsk )
1217#endif
1218     
1219      CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, pmask, zgmsk )
1220      CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, plam, zglam )
1221      CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, pphi, zgphi )
1222      ! Need to know the bathy depth for each observation for sco
1223      CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, fsdepw(:,:,:), &
1224         &                  zgdepw )
1225
1226      DO jobs = 1, kprofno
1227
1228         ! Skip bad profiles
1229         IF ( kpobsqc(jobs) >= kqc_cutoff ) CYCLE
1230
1231         ! Check if this observation is on a grid point
1232
1233         lgridobs = .FALSE.
1234         iig = -1
1235         ijg = -1
1236         DO jj = 1, 2
1237            DO ji = 1, 2
1238               IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) &
1239                  & .AND. &
1240                  & ( ABS( MOD( zglam(ji,jj,jobs) - pobslam(jobs),360.0) ) < 1.0e-6_wp ) &
1241                  & ) THEN
1242                  lgridobs = .TRUE.
1243                  iig = ji
1244                  ijg = jj
1245               ENDIF
1246            END DO
1247         END DO
1248
1249         ! Check if next to land
1250         IF (  ANY( zgmsk(1:2,1:2,1,jobs) == 0.0_wp ) ) THEN
1251            ll_next_to_land=.TRUE. 
1252         ELSE
1253            ll_next_to_land=.FALSE. 
1254         ENDIF 
1255         
1256         ! Reject observations
1257
1258         DO jobsp = kpstart(jobs), kpend(jobs)
1259
1260            ! Flag if the observation falls outside the model spatial domain
1261            IF (       ( pobslam(jobs) < -180.         )       &
1262               &  .OR. ( pobslam(jobs) >  180.         )       &
1263               &  .OR. ( pobsphi(jobs) <  -90.         )       &
1264               &  .OR. ( pobsphi(jobs) >   90.         )       &
1265               &  .OR. ( pobsdep(jobsp) < 0.0          )       &
1266               &  .OR. ( pobsdep(jobsp) > gdepw_1d(kpk)) ) THEN
1267               kobsqc(jobsp) = IBSET(kobsqc(jobsp),11)
1268               kosdobs = kosdobs + 1
1269               CYCLE
1270            ENDIF
1271
1272            ! To check if an observations falls within land there are two cases:
1273            ! 1: z-coordibnates, where the check uses the mask
1274            ! 2: terrain following (eg s-coordinates), 
1275            !    where we use the depth of the bottom cell to mask observations
1276             
1277            IF( (.NOT. lk_vvl) .AND. ( ln_zps .OR. ln_zco )  ) THEN !(CASE 1)
1278               
1279               ! Flag if the observation falls with a model land cell
1280               IF ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 
1281                  &  == 0.0_wp ) THEN
1282                  kobsqc(jobsp) = IBSET(kobsqc(jobsp),10)
1283                  klanobs = klanobs + 1 
1284                  CYCLE
1285               ENDIF 
1286             
1287               ! Flag if the observation is close to land
1288               IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == & 
1289                  &  0.0_wp) THEN
1290                  knlaobs = knlaobs + 1 
1291                  IF (ld_nea) THEN   
1292                     kobsqc(jobsp) = IBSET(kobsqc(jobsp),10)
1293                  ENDIF 
1294               ENDIF
1295             
1296            ELSE ! Case 2
1297               ! Flag if the observation is deeper than the bathymetry
1298               ! Or if it is within the mask
1299               IF ( ANY( zgdepw(1:2,1:2,kpk,jobs) < pobsdep(jobsp) ) &
1300                  &     .OR. & 
1301                  &  ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) &
1302                  &  == 0.0_wp) ) THEN
1303                  kobsqc(jobsp) = IBSET(kobsqc(jobsp),10)
1304                  klanobs = klanobs + 1 
1305                  CYCLE
1306               ENDIF 
1307               
1308               ! Flag if the observation is close to land
1309               IF ( ll_next_to_land ) THEN
1310                  knlaobs = knlaobs + 1 
1311                  IF (ld_nea) THEN   
1312                     kobsqc(jobsp) = IBSET(kobsqc(jobsp),10)
1313                  ENDIF 
1314               ENDIF
1315             
1316            ENDIF
1317
1318            ! For observations on the grid reject them if their are at
1319            ! a masked point
1320           
1321            IF (lgridobs) THEN
1322               IF (zgmsk(iig,ijg,kobsk(jobsp)-1,jobs) == 0.0_wp ) THEN
1323                  kobsqc(jobsp) = IBSET(kobsqc(jobs),10)
1324                  klanobs = klanobs + 1
1325                  CYCLE
1326               ENDIF
1327            ENDIF
1328           
1329            ! Set observation depth equal to that of the first model depth
1330            IF ( pobsdep(jobsp) <= pdep(1) ) THEN
1331               pobsdep(jobsp) = pdep(1) 
1332            ENDIF
1333           
1334#if defined key_bdy
1335            ! Flag if the observation falls close to the boundary rim
1336            IF (ld_bound_reject) THEN
1337               IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN
1338                  kobsqc(jobsp) = IBSET(kobsqc(jobs),8)
1339                  kbdyobs = kbdyobs + 1
1340                  CYCLE
1341               ENDIF
1342               ! for observations on the grid...
1343               IF (lgridobs) THEN
1344                  IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN
1345                     kobsqc(jobsp) = IBSET(kobsqc(jobs),8)
1346                     kbdyobs = kbdyobs + 1
1347                     CYCLE
1348                  ENDIF
1349               ENDIF
1350            ENDIF
1351#endif
1352           
1353         END DO
1354      END DO
1355
1356   END SUBROUTINE obs_coo_spc_3d
1357
1358   SUBROUTINE obs_pro_rej( profdata, kqc_cutoff )
1359      !!----------------------------------------------------------------------
1360      !!                    ***  ROUTINE obs_pro_rej ***
1361      !!
1362      !! ** Purpose : Reject all data within a rejected profile
1363      !!
1364      !! ** Method  :
1365      !!
1366      !! ** Action  :
1367      !!
1368      !! References :
1369      !!   
1370      !! History :
1371      !!        !  2007-10  (K. Mogensen) Original code
1372      !!----------------------------------------------------------------------
1373      !! * Modules used
1374      !! * Arguments
1375      TYPE(obs_prof), INTENT(INOUT) :: profdata     ! Profile data
1376      INTEGER, INTENT(IN) :: kqc_cutoff             ! QC cutoff value
1377
1378      !! * Local declarations
1379      INTEGER :: jprof
1380      INTEGER :: jvar
1381      INTEGER :: jobs
1382     
1383      ! Loop over profiles
1384
1385      DO jprof = 1, profdata%nprof
1386
1387         IF ( profdata%nqc(jprof) > kqc_cutoff ) THEN
1388           
1389            DO jvar = 1, profdata%nvar
1390
1391               DO jobs = profdata%npvsta(jprof,jvar),  &
1392                  &      profdata%npvend(jprof,jvar)
1393                 
1394                  profdata%var(jvar)%nvqc(jobs) = &
1395                     & IBSET(profdata%var(jvar)%nvqc(jobs),14)
1396
1397               END DO
1398
1399            END DO
1400
1401         ENDIF
1402
1403      END DO
1404
1405   END SUBROUTINE obs_pro_rej
1406
1407   SUBROUTINE obs_uv_rej( profdata, knumu, knumv, kqc_cutoff )
1408      !!----------------------------------------------------------------------
1409      !!                    ***  ROUTINE obs_uv_rej ***
1410      !!
1411      !! ** Purpose : Reject u if v is rejected and vice versa
1412      !!
1413      !! ** Method  :
1414      !!
1415      !! ** Action  :
1416      !!
1417      !! References :
1418      !!   
1419      !! History :
1420      !!        !  2009-2  (K. Mogensen) Original code
1421      !!----------------------------------------------------------------------
1422      !! * Modules used
1423      !! * Arguments
1424      TYPE(obs_prof), INTENT(INOUT) :: profdata   ! Profile data
1425      INTEGER, INTENT(INOUT) :: knumu             ! Number of u rejected
1426      INTEGER, INTENT(INOUT) :: knumv             ! Number of v rejected
1427      INTEGER, INTENT(IN) :: kqc_cutoff           ! QC cutoff value
1428
1429      !! * Local declarations
1430      INTEGER :: jprof
1431      INTEGER :: jvar
1432      INTEGER :: jobs
1433     
1434      ! Loop over profiles
1435
1436      DO jprof = 1, profdata%nprof
1437
1438         IF ( ( profdata%npvsta(jprof,1) /= profdata%npvsta(jprof,2) ) .OR. &
1439            & ( profdata%npvend(jprof,1) /= profdata%npvend(jprof,2) ) ) THEN
1440
1441            CALL ctl_stop('U,V profiles inconsistent in obs_uv_rej')
1442            RETURN
1443
1444         ENDIF
1445
1446         DO jobs = profdata%npvsta(jprof,1), profdata%npvend(jprof,1)
1447           
1448            IF ( ( profdata%var(1)%nvqc(jobs) >  kqc_cutoff ) .AND. &
1449               & ( profdata%var(2)%nvqc(jobs) <=  kqc_cutoff) ) THEN
1450               profdata%var(2)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15)
1451               knumv = knumv + 1
1452            ENDIF
1453            IF ( ( profdata%var(2)%nvqc(jobs) >  kqc_cutoff ) .AND. &
1454               & ( profdata%var(1)%nvqc(jobs) <=  kqc_cutoff) ) THEN
1455               profdata%var(1)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15)
1456               knumu = knumu + 1
1457            ENDIF
1458           
1459         END DO
1460           
1461      END DO
1462
1463   END SUBROUTINE obs_uv_rej
1464
1465END MODULE obs_prep
Note: See TracBrowser for help on using the repository browser.