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

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

source: NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/OBS/obs_prep.F90 @ 9923

Last change on this file since 9923 was 9923, checked in by gm, 6 years ago

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