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/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 9023

Last change on this file since 9023 was 9023, checked in by timgraham, 7 years ago

Merged METO_MERCATOR branch and resolved all conflicts in OPA_SRC
  • 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/OPA 3.3 , NEMO Consortium (2010)
37   !! $Id$
38   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
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 = 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 = 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      !! * Modules used
613      USE dom_oce, ONLY : &  ! Geographical information
614         & rdt
615      USE phycst, ONLY : &   ! Physical constants
616         & rday,  &             
617         & rmmss, &             
618         & rhhmm                       
619      !! * Arguments
620      INTEGER, INTENT(IN) :: kcycle     ! Current cycle
621      INTEGER, INTENT(IN) :: kyea0      ! Initial date coordinates
622      INTEGER, INTENT(IN) :: kmon0
623      INTEGER, INTENT(IN) :: kday0 
624      INTEGER, INTENT(IN) :: khou0
625      INTEGER, INTENT(IN) :: kmin0
626      INTEGER, INTENT(IN) :: kobsno     ! Number of observations
627      INTEGER, INTENT(INOUT) :: kotdobs   ! Number of observations failing time check
628      INTEGER, DIMENSION(kobsno), INTENT(IN ) :: &
629         & kobsyea,  &      ! Observation time coordinates
630         & kobsmon,  &
631         & kobsday,  & 
632         & kobshou,  &
633         & kobsmin
634      INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
635         & kobsqc           ! Quality control flag
636      INTEGER, DIMENSION(kobsno), INTENT(OUT) :: &
637         & kobsstp          ! Number of time steps up to the
638                            ! observation time
639
640      !! * Local declarations
641      INTEGER :: jyea
642      INTEGER :: jmon
643      INTEGER :: jday
644      INTEGER :: jobs
645      INTEGER :: iyeastr
646      INTEGER :: iyeaend
647      INTEGER :: imonstr
648      INTEGER :: imonend
649      INTEGER :: idaystr
650      INTEGER :: idayend 
651      INTEGER :: iskip
652      INTEGER :: idaystp
653      REAL(KIND=wp) :: zminstp
654      REAL(KIND=wp) :: zhoustp
655      REAL(KIND=wp) :: zobsstp 
656      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
657 
658      !-----------------------------------------------------------------------
659      ! Initialization
660      !-----------------------------------------------------------------------
661
662      ! Intialize the number of time steps per day
663      idaystp = NINT( rday / rdt )
664
665      !---------------------------------------------------------------------
666      ! Locate the model time coordinates for interpolation
667      !---------------------------------------------------------------------
668
669      DO jobs = 1, kobsno
670
671         ! Initialize the time step counter
672         kobsstp(jobs) = nit000 - 1
673
674         ! Flag if observation date is less than the initial date
675
676         IF ( ( kobsyea(jobs) < kyea0 )                   &
677            & .OR. ( ( kobsyea(jobs) == kyea0 )           &
678            &        .AND. ( kobsmon(jobs) <  kmon0 ) )   &
679            & .OR. ( ( kobsyea(jobs) == kyea0 )           &
680            &        .AND. ( kobsmon(jobs) == kmon0 )     &
681            &        .AND. ( kobsday(jobs) <  kday0 ) )   &
682            & .OR. ( ( kobsyea(jobs) == kyea0 )           &
683            &        .AND. ( kobsmon(jobs) == kmon0 )     &
684            &        .AND. ( kobsday(jobs) == kday0 )     &
685            &        .AND. ( kobshou(jobs) <  khou0 ) )   &
686            &  .OR. ( ( kobsyea(jobs) == kyea0 )          &
687            &        .AND. ( kobsmon(jobs) == kmon0 )     &
688            &        .AND. ( kobsday(jobs) == kday0 )          &
689            &        .AND. ( kobshou(jobs) == khou0 )          &
690            &        .AND. ( kobsmin(jobs) <= kmin0 ) ) ) THEN
691            kobsstp(jobs) = -1
692            kobsqc(jobs)  = IBSET(kobsqc(jobs),13)
693            kotdobs       = kotdobs + 1
694            CYCLE
695         ENDIF
696
697         ! Compute the number of time steps to the observation day
698         iyeastr = kyea0
699         iyeaend = kobsyea(jobs)
700         
701         !---------------------------------------------------------------------
702         ! Year loop
703         !---------------------------------------------------------------------
704         DO jyea = iyeastr, iyeaend
705
706            CALL calc_month_len( jyea, imonth_len )
707           
708            imonstr = 1
709            IF ( jyea == kyea0         ) imonstr = kmon0
710            imonend = 12
711            IF ( jyea == kobsyea(jobs) ) imonend = kobsmon(jobs)
712           
713            ! Month loop
714            DO jmon = imonstr, imonend
715               
716               idaystr = 1
717               IF (       ( jmon == kmon0   ) &
718                  & .AND. ( jyea == kyea0   ) ) idaystr = kday0
719               idayend = imonth_len(jmon)
720               IF (       ( jmon == kobsmon(jobs) ) &
721                  & .AND. ( jyea == kobsyea(jobs) ) ) idayend = kobsday(jobs) - 1
722               
723               ! Day loop
724               DO jday = idaystr, idayend
725                  kobsstp(jobs) = kobsstp(jobs) + idaystp
726               END DO
727               
728            END DO
729
730         END DO
731
732         ! Add in the number of time steps to the observation minute
733         zminstp = rmmss / rdt
734         zhoustp = rhhmm * zminstp
735
736         zobsstp =   REAL( kobsmin(jobs) - kmin0, KIND=wp ) * zminstp &
737            &      + REAL( kobshou(jobs) - khou0, KIND=wp ) * zhoustp
738         kobsstp(jobs) = kobsstp(jobs) + NINT( zobsstp )
739
740         ! Flag if observation step outside the time window
741         IF ( ( kobsstp(jobs) < ( nit000 - 1 ) ) &
742            & .OR.( kobsstp(jobs) > nitend ) ) THEN
743            kobsqc(jobs) = IBSET(kobsqc(jobs),13)
744            kotdobs = kotdobs + 1
745            CYCLE
746         ENDIF
747
748      END DO
749
750   END SUBROUTINE obs_coo_tim
751
752   SUBROUTINE calc_month_len( iyear, imonth_len )
753      !!----------------------------------------------------------------------
754      !!                    ***  ROUTINE calc_month_len  ***
755      !!         
756      !! ** Purpose : Compute the number of days in a months given a year.
757      !!
758      !! ** Method  :
759      !!
760      !! ** Action  :
761      !!
762      !! History :
763      !!        !  10-05  (D. Lea)   New routine based on day_init
764      !!----------------------------------------------------------------------
765
766      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
767      INTEGER :: iyear         !: year
768     
769      ! length of the month of the current year (from nleapy, read in namelist)
770      IF ( nleapy < 2 ) THEN
771         imonth_len(:) = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)
772         IF ( nleapy == 1 ) THEN   ! we are using calendar with leap years
773            IF ( MOD(iyear, 4) == 0 .AND. ( MOD(iyear, 400) == 0 .OR. MOD(iyear, 100) /= 0 ) ) THEN
774               imonth_len(2) = 29
775            ENDIF
776         ENDIF
777      ELSE
778         imonth_len(:) = nleapy   ! all months with nleapy days per year
779      ENDIF
780
781   END SUBROUTINE
782
783   SUBROUTINE obs_coo_tim_prof( kcycle,                                   &
784      &                    kyea0,   kmon0,   kday0,   khou0,   kmin0,     &
785      &                    kobsno,                                        &
786      &                    kobsyea, kobsmon, kobsday, kobshou, kobsmin,   &
787      &                    ktyp,    kobsqc,  kobsstp, kotdobs, kdailyavtypes, &
788      &                    kqc_cutoff )
789      !!----------------------------------------------------------------------
790      !!                    ***  ROUTINE obs_coo_tim ***
791      !!
792      !! ** Purpose : Compute the number of time steps to the observation time.
793      !!
794      !! ** Method  : For time coordinates ( yea_obs, mon_obs, day_obs,
795      !!              hou_obs, min_obs ), this routine locates the time step
796      !!              that is closest to this time.
797      !!
798      !! ** Action  :
799      !!
800      !! References :
801      !!   
802      !! History :
803      !!        !  1997-07  (A. Weaver) Original
804      !!        !  2006-08  (A. Weaver) NEMOVAR migration
805      !!        !  2006-10  (A. Weaver) Cleanup
806      !!        !  2007-01  (K. Mogensen) Rewritten with loop
807      !!----------------------------------------------------------------------
808      !! * Modules used
809      !! * Arguments
810      INTEGER, INTENT(IN) :: kcycle      ! Current cycle
811      INTEGER, INTENT(IN) :: kyea0       ! Initial date coordinates
812      INTEGER, INTENT(IN) :: kmon0
813      INTEGER, INTENT(IN) :: kday0
814      INTEGER, INTENT(IN) :: khou0
815      INTEGER, INTENT(IN) :: kmin0
816      INTEGER, INTENT(IN) :: kobsno      ! Number of observations
817      INTEGER, INTENT(INOUT) ::  kotdobs    ! Number of observations failing time check
818      INTEGER, DIMENSION(kobsno), INTENT(IN ) :: &
819         & kobsyea,  &      ! Observation time coordinates
820         & kobsmon,  &
821         & kobsday,  & 
822         & kobshou,  &
823         & kobsmin,  &
824         & ktyp             ! Observation type.
825      INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
826         & kobsqc           ! Quality control flag
827      INTEGER, DIMENSION(kobsno), INTENT(OUT) :: &
828         & kobsstp          ! Number of time steps up to the
829                            ! observation time
830      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
831         & kdailyavtypes    ! Types for daily averages
832      INTEGER, OPTIONAL, INTENT(IN) :: kqc_cutoff           ! QC cutoff value
833
834      !! * Local declarations
835      INTEGER :: jobs
836      INTEGER :: iqc_cutoff=255
837
838      !-----------------------------------------------------------------------
839      ! Call standard obs_coo_tim
840      !-----------------------------------------------------------------------
841
842      CALL obs_coo_tim( kcycle, &
843         &              kyea0,   kmon0,   kday0,   khou0,   kmin0,     &
844         &              kobsno,                                        &
845         &              kobsyea, kobsmon, kobsday, kobshou, kobsmin,   &
846         &              kobsqc,  kobsstp, kotdobs                      )
847
848      !------------------------------------------------------------------------
849      ! Always reject daily averaged data (e.g. MRB data (820)) at initial time
850      !------------------------------------------------------------------------
851
852      IF ( PRESENT(kdailyavtypes) ) THEN
853         DO jobs = 1, kobsno
854           
855            IF ( kobsqc(jobs) <= iqc_cutoff ) THEN
856               
857               IF ( ( kobsstp(jobs) == (nit000 - 1) ).AND.&
858                  & ( ANY (kdailyavtypes(:) == ktyp(jobs)) ) ) THEN
859                  kobsqc(jobs) = IBSET(kobsqc(jobs),13)
860                  kotdobs      = kotdobs + 1
861                  CYCLE
862               ENDIF
863               
864            ENDIF
865         END DO
866      ENDIF
867
868
869   END SUBROUTINE obs_coo_tim_prof
870
871   SUBROUTINE obs_coo_grd( kobsno, kobsi, kobsj, kobsqc, kgrdobs )
872      !!----------------------------------------------------------------------
873      !!                    ***  ROUTINE obs_coo_grd ***
874      !!
875      !! ** Purpose : Verify that the grid search has not failed
876      !!
877      !! ** Method  : The previously computed i,j indeces are checked 
878      !!
879      !! ** Action  :
880      !!
881      !! References :
882      !!   
883      !! History :
884      !!        !  2007-01  (K. Mogensen)  Original
885      !!----------------------------------------------------------------------
886      !! * Arguments
887      INTEGER, INTENT(IN) :: kobsno        ! Number of observations
888      INTEGER, DIMENSION(kobsno), INTENT(IN ) :: &
889         & kobsi, &         ! i,j indeces previously computed
890         & kobsj
891      INTEGER, INTENT(INOUT) ::  kgrdobs   ! Number of observations failing the check
892      INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
893         & kobsqc           ! Quality control flag
894
895      !! * Local declarations
896      INTEGER :: jobs       ! Loop variable
897
898      ! Flag if the grid search failed
899
900      DO jobs = 1, kobsno
901         IF ( ( kobsi(jobs) <= 0 ) .AND. ( kobsj(jobs) <= 0 ) ) THEN
902            kobsqc(jobs) = IBSET(kobsqc(jobs),12)
903            kgrdobs = kgrdobs + 1
904         ENDIF
905      END DO
906     
907   END SUBROUTINE obs_coo_grd
908
909   SUBROUTINE obs_coo_spc_2d( kobsno, kpi,     kpj,              &
910      &                       kobsi,  kobsj,   pobslam, pobsphi, & 
911      &                       plam,   pphi,    pmask,            &
912      &                       kobsqc, kosdobs, klanobs,          &
913      &                       knlaobs,ld_nea,                    &
914      &                       kbdyobs,ld_bound_reject,           &
915      &                       kqc_cutoff                         )
916      !!----------------------------------------------------------------------
917      !!                    ***  ROUTINE obs_coo_spc_2d  ***
918      !!
919      !! ** Purpose : Check for points outside the domain and land points
920      !!
921      !! ** Method  : Remove the observations that are outside the model space
922      !!              and time domain or located within model land cells.
923      !!
924      !! ** Action  :
925      !!   
926      !! History :
927      !!        !  2007-03  (A. Weaver, K. Mogensen)  Original
928      !!        !  2007-06  (K. Mogensen et al) Reject obs. near land.
929      !!----------------------------------------------------------------------
930      !! * Modules used
931
932      !! * Arguments
933      INTEGER, INTENT(IN) :: kobsno    ! Total number of observations
934      INTEGER, INTENT(IN) :: kpi       ! Number of grid points in (i,j)
935      INTEGER, INTENT(IN) :: kpj
936      INTEGER, DIMENSION(kobsno), INTENT(IN) :: &
937         & kobsi, &           ! Observation (i,j) coordinates
938         & kobsj
939      REAL(KIND=wp), DIMENSION(kobsno), INTENT(IN) :: &
940         & pobslam, &         ! Observation (lon,lat) coordinates
941         & pobsphi
942      REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) :: &
943         & plam, pphi         ! Model (lon,lat) coordinates
944      REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) :: &
945         & pmask              ! Land mask array
946      INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
947         & kobsqc             ! Observation quality control
948      INTEGER, INTENT(INOUT) :: kosdobs          ! Observations outside space domain
949      INTEGER, INTENT(INOUT) :: klanobs          ! Observations within a model land cell
950      INTEGER, INTENT(INOUT) :: knlaobs          ! Observations near land
951      INTEGER, INTENT(INOUT) :: kbdyobs          ! Observations near boundary
952      LOGICAL, INTENT(IN)    :: ld_nea           ! Flag observations near land
953      LOGICAL, INTENT(IN)    :: ld_bound_reject  ! Flag observations near open boundary
954      INTEGER, INTENT(IN)    :: kqc_cutoff       ! Cutoff QC value
955
956      !! * Local declarations
957      REAL(KIND=wp), DIMENSION(2,2,kobsno) :: &
958         & zgmsk              ! Grid mask
959
960      REAL(KIND=wp), DIMENSION(2,2,kobsno) :: &
961         & zbmsk              ! Boundary mask
962      REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask
963      REAL(KIND=wp), DIMENSION(2,2,kobsno) :: &
964         & zglam, &           ! Model longitude at grid points
965         & zgphi              ! Model latitude at grid points
966      INTEGER, DIMENSION(2,2,kobsno) :: &
967         & igrdi, &           ! Grid i,j
968         & igrdj
969      LOGICAL :: lgridobs           ! Is observation on a model grid point.
970      INTEGER :: iig, ijg           ! i,j of observation on model grid point.
971      INTEGER :: jobs, ji, jj
972     
973      ! Get grid point indices
974
975      DO jobs = 1, kobsno
976         
977         ! For invalid points use 2,2
978
979         IF ( kobsqc(jobs) >= kqc_cutoff ) THEN
980
981            igrdi(1,1,jobs) = 1
982            igrdj(1,1,jobs) = 1
983            igrdi(1,2,jobs) = 1
984            igrdj(1,2,jobs) = 2
985            igrdi(2,1,jobs) = 2
986            igrdj(2,1,jobs) = 1
987            igrdi(2,2,jobs) = 2
988            igrdj(2,2,jobs) = 2
989
990         ELSE
991
992            igrdi(1,1,jobs) = kobsi(jobs)-1
993            igrdj(1,1,jobs) = kobsj(jobs)-1
994            igrdi(1,2,jobs) = kobsi(jobs)-1
995            igrdj(1,2,jobs) = kobsj(jobs)
996            igrdi(2,1,jobs) = kobsi(jobs)
997            igrdj(2,1,jobs) = kobsj(jobs)-1
998            igrdi(2,2,jobs) = kobsi(jobs)
999            igrdj(2,2,jobs) = kobsj(jobs)
1000
1001         ENDIF
1002
1003      END DO
1004
1005      IF (ln_bdy) THEN
1006        ! Create a mask grid points in boundary rim
1007        IF (ld_bound_reject) THEN
1008           zbdymask(:,:) = 1.0_wp
1009           DO ji = 1, nb_bdy
1010              DO jj = 1, idx_bdy(ji)%nblen(1)
1011                 zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp
1012              ENDDO
1013           ENDDO
1014
1015           CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, zbdymask, zbmsk )
1016        ENDIF
1017      ENDIF
1018
1019     
1020      CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, pmask, zgmsk )
1021      CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, plam, zglam )
1022      CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, pphi, zgphi )
1023
1024      DO jobs = 1, kobsno
1025
1026         ! Skip bad observations
1027         IF ( kobsqc(jobs) >= kqc_cutoff ) CYCLE
1028
1029         ! Flag if the observation falls outside the model spatial domain
1030         IF (       ( pobslam(jobs) < -180. ) &
1031            &  .OR. ( pobslam(jobs) >  180. ) &
1032            &  .OR. ( pobsphi(jobs) <  -90. ) &
1033            &  .OR. ( pobsphi(jobs) >   90. ) ) THEN
1034            kobsqc(jobs) = IBSET(kobsqc(jobs),11)
1035            kosdobs = kosdobs + 1
1036            CYCLE
1037         ENDIF
1038
1039         ! Flag if the observation falls with a model land cell
1040         IF ( SUM( zgmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN
1041            kobsqc(jobs) = IBSET(kobsqc(jobs),10)
1042            klanobs = klanobs + 1
1043            CYCLE
1044         ENDIF
1045
1046         ! Check if this observation is on a grid point
1047
1048         lgridobs = .FALSE.
1049         iig = -1
1050         ijg = -1
1051         DO jj = 1, 2
1052            DO ji = 1, 2
1053               IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) &
1054                  & .AND. &
1055                  & ( ABS( MOD( zglam(ji,jj,jobs) - pobslam(jobs),360.0) )  &
1056                  & < 1.0e-6_wp ) ) THEN
1057                  lgridobs = .TRUE.
1058                  iig = ji
1059                  ijg = jj
1060               ENDIF
1061            END DO
1062         END DO
1063 
1064         ! For observations on the grid reject them if their are at
1065         ! a masked point
1066         
1067         IF (lgridobs) THEN
1068            IF (zgmsk(iig,ijg,jobs) == 0.0_wp ) THEN
1069               kobsqc(jobs) = IBSET(kobsqc(jobs),10)
1070               klanobs = klanobs + 1
1071               CYCLE
1072            ENDIF
1073         ENDIF
1074                     
1075         ! Flag if the observation falls is close to land
1076         IF ( MINVAL( zgmsk(1:2,1:2,jobs) ) == 0.0_wp) THEN
1077            knlaobs = knlaobs + 1
1078            IF (ld_nea) THEN
1079               kobsqc(jobs) = IBSET(kobsqc(jobs),9)
1080               CYCLE
1081            ENDIF
1082         ENDIF
1083
1084         IF (ln_bdy) THEN
1085         ! Flag if the observation falls close to the boundary rim
1086           IF (ld_bound_reject) THEN
1087              IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN
1088                 kobsqc(jobs) = IBSET(kobsqc(jobs),8)
1089                 kbdyobs = kbdyobs + 1
1090                 CYCLE
1091              ENDIF
1092              ! for observations on the grid...
1093              IF (lgridobs) THEN
1094                 IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN
1095                    kobsqc(jobs) = IBSET(kobsqc(jobs),8)
1096                    kbdyobs = kbdyobs + 1
1097                    CYCLE
1098                 ENDIF
1099              ENDIF
1100            ENDIF
1101         ENDIF
1102           
1103      END DO
1104
1105   END SUBROUTINE obs_coo_spc_2d
1106
1107   SUBROUTINE obs_coo_spc_3d( kprofno, kobsno,  kpstart, kpend, &
1108      &                       kpi,     kpj,     kpk,            &
1109      &                       kobsi,   kobsj,   kobsk,          &
1110      &                       pobslam, pobsphi, pobsdep,        &
1111      &                       plam,    pphi,    pdep,    pmask, &
1112      &                       kpobsqc, kobsqc,  kosdobs,        &
1113      &                       klanobs, knlaobs, ld_nea,         &
1114      &                       kbdyobs, ld_bound_reject,         &
1115      &                       kqc_cutoff                        )
1116      !!----------------------------------------------------------------------
1117      !!                    ***  ROUTINE obs_coo_spc_3d  ***
1118      !!
1119      !! ** Purpose : Check for points outside the domain and land points
1120      !!              Reset depth of observation above highest model level
1121      !!              to the value of highest model level
1122      !!
1123      !! ** Method  : Remove the observations that are outside the model space
1124      !!              and time domain or located within model land cells.
1125      !!
1126      !!              NB. T and S profile observations lying between the ocean
1127      !!              surface and the depth of the first model T point are
1128      !!              assigned a depth equal to that of the first model T pt.
1129      !!
1130      !! ** Action  :
1131      !!   
1132      !! History :
1133      !!        !  2007-01  (K. Mogensen) Rewrite of parts of obs_scr
1134      !!        !  2007-06  (K. Mogensen et al) Reject obs. near land.
1135      !!----------------------------------------------------------------------
1136      !! * Modules used
1137      USE dom_oce, ONLY : &       ! Geographical information
1138         & gdepw_1d,      &
1139         & gdepw_0,       &                       
1140         & gdepw_n,       &
1141         & gdept_n,       &
1142         & ln_zco,        &
1143         & ln_zps             
1144
1145      !! * Arguments
1146      INTEGER, INTENT(IN) :: kprofno      ! Number of profiles
1147      INTEGER, INTENT(IN) :: kobsno       ! Total number of observations
1148      INTEGER, INTENT(IN) :: kpi          ! Number of grid points in (i,j,k)
1149      INTEGER, INTENT(IN) :: kpj
1150      INTEGER, INTENT(IN) :: kpk   
1151      INTEGER, DIMENSION(kprofno), INTENT(IN) :: &
1152         & kpstart, &         ! Start of individual profiles
1153         & kpend              ! End of individual profiles
1154      INTEGER, DIMENSION(kprofno), INTENT(IN) :: &
1155         & kobsi, &           ! Observation (i,j) coordinates
1156         & kobsj
1157      INTEGER, DIMENSION(kobsno), INTENT(IN) :: &
1158         & kobsk              ! Observation k coordinate
1159      REAL(KIND=wp), DIMENSION(kprofno), INTENT(IN) :: &
1160         & pobslam, &         ! Observation (lon,lat) coordinates
1161         & pobsphi
1162      REAL(KIND=wp), DIMENSION(kobsno), INTENT(INOUT) :: &
1163         & pobsdep            ! Observation depths 
1164      REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) :: &
1165         & plam, pphi         ! Model (lon,lat) coordinates
1166      REAL(KIND=wp), DIMENSION(kpk), INTENT(IN) :: &
1167         & pdep               ! Model depth coordinates
1168      REAL(KIND=wp), DIMENSION(kpi,kpj,kpk), INTENT(IN) :: &
1169         & pmask              ! Land mask array
1170      INTEGER, DIMENSION(kprofno), INTENT(INOUT) :: &
1171         & kpobsqc             ! Profile quality control
1172      INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
1173         & kobsqc             ! Observation quality control
1174      INTEGER, INTENT(INOUT) :: kosdobs     ! Observations outside space domain
1175      INTEGER, INTENT(INOUT) :: klanobs     ! Observations within a model land cell
1176      INTEGER, INTENT(INOUT) :: knlaobs     ! Observations near land
1177      INTEGER, INTENT(INOUT) :: kbdyobs     ! Observations near boundary
1178      LOGICAL, INTENT(IN) :: ld_nea         ! Flag observations near land
1179      LOGICAL, INTENT(IN) :: ld_bound_reject  ! Flag observations near open boundary
1180      INTEGER, INTENT(IN) :: kqc_cutoff     ! Cutoff QC value
1181
1182      !! * Local declarations
1183      REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: &
1184         & zgmsk              ! Grid mask
1185      REAL(KIND=wp), DIMENSION(2,2,kprofno) :: &
1186         & zbmsk              ! Boundary mask
1187      REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask
1188      REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: &
1189         & zgdepw
1190      REAL(KIND=wp), DIMENSION(2,2,kprofno) :: &
1191         & zglam, &           ! Model longitude at grid points
1192         & zgphi              ! Model latitude at grid points
1193      INTEGER, DIMENSION(2,2,kprofno) :: &
1194         & igrdi, &           ! Grid i,j
1195         & igrdj
1196      LOGICAL :: lgridobs           ! Is observation on a model grid point.
1197      LOGICAL :: ll_next_to_land    ! Is a profile next to land
1198      INTEGER :: iig, ijg           ! i,j of observation on model grid point.
1199      INTEGER :: jobs, jobsp, jk, ji, jj
1200
1201      ! Get grid point indices
1202     
1203      DO jobs = 1, kprofno
1204
1205         ! For invalid points use 2,2
1206
1207         IF ( kpobsqc(jobs) >= kqc_cutoff ) THEN
1208           
1209            igrdi(1,1,jobs) = 1
1210            igrdj(1,1,jobs) = 1
1211            igrdi(1,2,jobs) = 1
1212            igrdj(1,2,jobs) = 2
1213            igrdi(2,1,jobs) = 2
1214            igrdj(2,1,jobs) = 1
1215            igrdi(2,2,jobs) = 2
1216            igrdj(2,2,jobs) = 2
1217           
1218         ELSE
1219           
1220            igrdi(1,1,jobs) = kobsi(jobs)-1
1221            igrdj(1,1,jobs) = kobsj(jobs)-1
1222            igrdi(1,2,jobs) = kobsi(jobs)-1
1223            igrdj(1,2,jobs) = kobsj(jobs)
1224            igrdi(2,1,jobs) = kobsi(jobs)
1225            igrdj(2,1,jobs) = kobsj(jobs)-1
1226            igrdi(2,2,jobs) = kobsi(jobs)
1227            igrdj(2,2,jobs) = kobsj(jobs)
1228           
1229         ENDIF
1230         
1231      END DO
1232
1233      IF (ln_bdy) THEN
1234        ! Create a mask grid points in boundary rim
1235        IF (ld_bound_reject) THEN           
1236           zbdymask(:,:) = 1.0_wp
1237           DO ji = 1, nb_bdy
1238              DO jj = 1, idx_bdy(ji)%nblen(1)
1239                 zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp
1240              ENDDO
1241           ENDDO
1242        ENDIF
1243   
1244        CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, zbdymask, zbmsk )
1245      ENDIF
1246     
1247      CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, pmask, zgmsk )
1248      CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, plam, zglam )
1249      CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, pphi, zgphi )
1250      CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, gdepw_n(:,:,:), &
1251        &                     zgdepw )
1252
1253      DO jobs = 1, kprofno
1254
1255         ! Skip bad profiles
1256         IF ( kpobsqc(jobs) >= kqc_cutoff ) CYCLE
1257
1258         ! Check if this observation is on a grid point
1259
1260         lgridobs = .FALSE.
1261         iig = -1
1262         ijg = -1
1263         DO jj = 1, 2
1264            DO ji = 1, 2
1265               IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) &
1266                  & .AND. &
1267                  & ( ABS( MOD( zglam(ji,jj,jobs) - pobslam(jobs),360.0) ) < 1.0e-6_wp ) &
1268                  & ) THEN
1269                  lgridobs = .TRUE.
1270                  iig = ji
1271                  ijg = jj
1272               ENDIF
1273            END DO
1274         END DO
1275
1276         ! Check if next to land
1277         IF (  ANY( zgmsk(1:2,1:2,1,jobs) == 0.0_wp ) ) THEN
1278            ll_next_to_land=.TRUE.
1279         ELSE
1280            ll_next_to_land=.FALSE.
1281         ENDIF 
1282
1283         ! Reject observations
1284
1285         DO jobsp = kpstart(jobs), kpend(jobs)
1286
1287            ! Flag if the observation falls outside the model spatial domain
1288            IF (       ( pobslam(jobs) < -180.         )       &
1289               &  .OR. ( pobslam(jobs) >  180.         )       &
1290               &  .OR. ( pobsphi(jobs) <  -90.         )       &
1291               &  .OR. ( pobsphi(jobs) >   90.         )       &
1292               &  .OR. ( pobsdep(jobsp) < 0.0          )       &
1293               &  .OR. ( pobsdep(jobsp) > gdepw_1d(kpk)) ) THEN
1294               kobsqc(jobsp) = IBSET(kobsqc(jobsp),11)
1295               kosdobs = kosdobs + 1
1296               CYCLE
1297            ENDIF
1298
1299            ! To check if an observations falls within land:
1300             
1301            ! Flag if the observation is deeper than the bathymetry
1302            ! Or if it is within the mask
1303            IF ( ALL( zgdepw(1:2,1:2,kpk,jobs) < pobsdep(jobsp) ) &
1304               &     .OR. &
1305               &  ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) &
1306               &  == 0.0_wp) ) THEN
1307               kobsqc(jobsp) = IBSET(kobsqc(jobsp),10)
1308               klanobs = klanobs + 1
1309               CYCLE
1310            ENDIF
1311               
1312            ! Flag if the observation is close to land
1313            IF ( ll_next_to_land ) THEN
1314               knlaobs = knlaobs + 1
1315               IF (ld_nea) THEN   
1316                  kobsqc(jobsp) = IBSET(kobsqc(jobsp),10)
1317               ENDIF
1318            ENDIF
1319           
1320            ! For observations on the grid reject them if their are at
1321            ! a masked point
1322           
1323            IF (lgridobs) THEN
1324               IF (zgmsk(iig,ijg,kobsk(jobsp)-1,jobs) == 0.0_wp ) THEN
1325                  kobsqc(jobsp) = IBSET(kobsqc(jobs),10)
1326                  klanobs = klanobs + 1
1327                  CYCLE
1328               ENDIF
1329            ENDIF
1330           
1331            ! Flag if the observation falls is close to land
1332            IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == &
1333               &  0.0_wp) THEN
1334               IF (ld_nea) kobsqc(jobsp) = kobsqc(jobsp) + 14
1335               knlaobs = knlaobs + 1
1336            ENDIF
1337
1338            ! Set observation depth equal to that of the first model depth
1339            IF ( pobsdep(jobsp) <= pdep(1) ) THEN
1340               pobsdep(jobsp) = pdep(1) 
1341            ENDIF
1342           
1343            IF (ln_bdy) THEN
1344               ! Flag if the observation falls close to the boundary rim
1345               IF (ld_bound_reject) THEN
1346                  IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN
1347                     kobsqc(jobsp) = IBSET(kobsqc(jobs),8)
1348                     kbdyobs = kbdyobs + 1
1349                     CYCLE
1350                  ENDIF
1351                  ! for observations on the grid...
1352                  IF (lgridobs) THEN
1353                     IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN
1354                        kobsqc(jobsp) = IBSET(kobsqc(jobs),8)
1355                        kbdyobs = kbdyobs + 1
1356                        CYCLE
1357                     ENDIF
1358                  ENDIF
1359               ENDIF
1360            ENDIF
1361           
1362         END DO
1363      END DO
1364
1365   END SUBROUTINE obs_coo_spc_3d
1366
1367   SUBROUTINE obs_pro_rej( profdata, kqc_cutoff )
1368      !!----------------------------------------------------------------------
1369      !!                    ***  ROUTINE obs_pro_rej ***
1370      !!
1371      !! ** Purpose : Reject all data within a rejected profile
1372      !!
1373      !! ** Method  :
1374      !!
1375      !! ** Action  :
1376      !!
1377      !! References :
1378      !!   
1379      !! History :
1380      !!        !  2007-10  (K. Mogensen) Original code
1381      !!----------------------------------------------------------------------
1382      !! * Modules used
1383      !! * Arguments
1384      TYPE(obs_prof), INTENT(INOUT) :: profdata     ! Profile data
1385      INTEGER, INTENT(IN) :: kqc_cutoff             ! QC cutoff value
1386
1387      !! * Local declarations
1388      INTEGER :: jprof
1389      INTEGER :: jvar
1390      INTEGER :: jobs
1391     
1392      ! Loop over profiles
1393
1394      DO jprof = 1, profdata%nprof
1395
1396         IF ( profdata%nqc(jprof) > kqc_cutoff ) THEN
1397           
1398            DO jvar = 1, profdata%nvar
1399
1400               DO jobs = profdata%npvsta(jprof,jvar),  &
1401                  &      profdata%npvend(jprof,jvar)
1402                 
1403                  profdata%var(jvar)%nvqc(jobs) = &
1404                     & IBSET(profdata%var(jvar)%nvqc(jobs),14)
1405
1406               END DO
1407
1408            END DO
1409
1410         ENDIF
1411
1412      END DO
1413
1414   END SUBROUTINE obs_pro_rej
1415
1416   SUBROUTINE obs_uv_rej( profdata, knumu, knumv, kqc_cutoff )
1417      !!----------------------------------------------------------------------
1418      !!                    ***  ROUTINE obs_uv_rej ***
1419      !!
1420      !! ** Purpose : Reject u if v is rejected and vice versa
1421      !!
1422      !! ** Method  :
1423      !!
1424      !! ** Action  :
1425      !!
1426      !! References :
1427      !!   
1428      !! History :
1429      !!        !  2009-2  (K. Mogensen) Original code
1430      !!----------------------------------------------------------------------
1431      !! * Modules used
1432      !! * Arguments
1433      TYPE(obs_prof), INTENT(INOUT) :: profdata   ! Profile data
1434      INTEGER, INTENT(INOUT) :: knumu             ! Number of u rejected
1435      INTEGER, INTENT(INOUT) :: knumv             ! Number of v rejected
1436      INTEGER, INTENT(IN) :: kqc_cutoff           ! QC cutoff value
1437
1438      !! * Local declarations
1439      INTEGER :: jprof
1440      INTEGER :: jvar
1441      INTEGER :: jobs
1442     
1443      ! Loop over profiles
1444
1445      DO jprof = 1, profdata%nprof
1446
1447         IF ( ( profdata%npvsta(jprof,1) /= profdata%npvsta(jprof,2) ) .OR. &
1448            & ( profdata%npvend(jprof,1) /= profdata%npvend(jprof,2) ) ) THEN
1449
1450            CALL ctl_stop('U,V profiles inconsistent in obs_uv_rej')
1451            RETURN
1452
1453         ENDIF
1454
1455         DO jobs = profdata%npvsta(jprof,1), profdata%npvend(jprof,1)
1456           
1457            IF ( ( profdata%var(1)%nvqc(jobs) >  kqc_cutoff ) .AND. &
1458               & ( profdata%var(2)%nvqc(jobs) <=  kqc_cutoff) ) THEN
1459               profdata%var(2)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15)
1460               knumv = knumv + 1
1461            ENDIF
1462            IF ( ( profdata%var(2)%nvqc(jobs) >  kqc_cutoff ) .AND. &
1463               & ( profdata%var(1)%nvqc(jobs) <=  kqc_cutoff) ) THEN
1464               profdata%var(1)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15)
1465               knumu = knumu + 1
1466            ENDIF
1467           
1468         END DO
1469           
1470      END DO
1471
1472   END SUBROUTINE obs_uv_rej
1473
1474END MODULE obs_prep
Note: See TracBrowser for help on using the repository browser.