source: branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 11200

Last change on this file since 11200 was 11200, checked in by emmafiedler, 15 months ago

Avoid multiple lines of write statement and change of variable name following CICE review

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