New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
obs_prep.F90 in branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 15400

Last change on this file since 15400 was 15400, checked in by kingr, 12 months ago

Bug-fix to obsoper to avoid division by zero in rare cases.

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