source: NEMO/trunk/src/OCE/OBS/obs_prep.F90 @ 10425

Last change on this file since 10425 was 10068, checked in by nicolasmartin, 3 years ago

First part of modifications to have a common default header : fix typos and SVN keywords properties

  • Property svn:keywords set to Id
File size: 60.2 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 license (see ./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      !! * 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      INTEGER, INTENT(IN) :: kobsno        ! Number of observations
887      INTEGER, DIMENSION(kobsno), INTENT(IN ) :: &
888         & kobsi, &         ! i,j indeces previously computed
889         & kobsj
890      INTEGER, INTENT(INOUT) ::  kgrdobs   ! Number of observations failing the check
891      INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
892         & kobsqc           ! Quality control flag
893
894      !! * Local declarations
895      INTEGER :: jobs       ! Loop variable
896
897      ! Flag if the grid search failed
898
899      DO jobs = 1, kobsno
900         IF ( ( kobsi(jobs) <= 0 ) .AND. ( kobsj(jobs) <= 0 ) ) THEN
901            kobsqc(jobs) = IBSET(kobsqc(jobs),12)
902            kgrdobs = kgrdobs + 1
903         ENDIF
904      END DO
905     
906   END SUBROUTINE obs_coo_grd
907
908   SUBROUTINE obs_coo_spc_2d( kobsno, kpi,     kpj,              &
909      &                       kobsi,  kobsj,   pobslam, pobsphi, & 
910      &                       plam,   pphi,    pmask,            &
911      &                       kobsqc, kosdobs, klanobs,          &
912      &                       knlaobs,ld_nea,                    &
913      &                       kbdyobs,ld_bound_reject,           &
914      &                       kqc_cutoff                         )
915      !!----------------------------------------------------------------------
916      !!                    ***  ROUTINE obs_coo_spc_2d  ***
917      !!
918      !! ** Purpose : Check for points outside the domain and land points
919      !!
920      !! ** Method  : Remove the observations that are outside the model space
921      !!              and time domain or located within model land cells.
922      !!
923      !! ** Action  :
924      !!   
925      !! History :  2007-03  (A. Weaver, K. Mogensen)  Original
926      !!         !  2007-06  (K. Mogensen et al) Reject obs. near land.
927      !!----------------------------------------------------------------------
928      INTEGER , INTENT(in   )                     ::   kobsno            ! Total number of observations
929      INTEGER , INTENT(in   )                     ::   kpi    , kpj      ! Number of grid points in (i,j)
930      INTEGER , INTENT(in   ), DIMENSION(kobsno)  ::   kobsi  , kobsj    ! Observation (i,j) coordinates
931      REAL(wp), INTENT(in   ), DIMENSION(kobsno)  ::   pobslam, pobsphi  ! Observation (lon,lat) coordinates
932      REAL(wp), INTENT(in   ), DIMENSION(kpi,kpj) ::   plam   , pphi     ! Model (lon,lat) coordinates
933      REAL(wp), INTENT(in   ), DIMENSION(kpi,kpj) ::   pmask             ! Land mask array
934      INTEGER , INTENT(inout), DIMENSION(kobsno)  ::   kobsqc            ! Observation quality control
935      INTEGER , INTENT(inout)                     ::   kosdobs           ! Observations outside space domain
936      INTEGER , INTENT(inout)                     ::   klanobs           ! Observations within a model land cell
937      INTEGER , INTENT(inout)                     ::   knlaobs           ! Observations near land
938      INTEGER , INTENT(inout)                     ::   kbdyobs           ! Observations near boundary
939      LOGICAL , INTENT(in   )                     ::   ld_nea            ! Flag observations near land
940      LOGICAL , INTENT(in   )                     ::   ld_bound_reject   ! Flag observations near open boundary
941      INTEGER , INTENT(in   )                     ::   kqc_cutoff        ! Cutoff QC value
942      !
943      REAL(KIND=wp), DIMENSION(2,2,kobsno) ::   zgmsk          ! Grid mask
944      REAL(KIND=wp), DIMENSION(2,2,kobsno) ::   zbmsk          ! Boundary mask
945      REAL(KIND=wp), DIMENSION(jpi,jpj)    ::   zbdymask
946      REAL(KIND=wp), DIMENSION(2,2,kobsno) ::   zglam, zgphi   ! Model Lon/lat at grid points
947      INTEGER      , DIMENSION(2,2,kobsno) ::   igrdi, igrdj   ! Grid i,j
948      LOGICAL ::   lgridobs           ! Is observation on a model grid point.
949      INTEGER ::   iig, ijg           ! i,j of observation on model grid point.
950      INTEGER ::   jobs, ji, jj
951      !!----------------------------------------------------------------------
952     
953      ! Get grid point indices
954
955      DO jobs = 1, kobsno
956         
957         ! For invalid points use 2,2
958
959         IF ( kobsqc(jobs) >= kqc_cutoff ) THEN
960
961            igrdi(1,1,jobs) = 1
962            igrdj(1,1,jobs) = 1
963            igrdi(1,2,jobs) = 1
964            igrdj(1,2,jobs) = 2
965            igrdi(2,1,jobs) = 2
966            igrdj(2,1,jobs) = 1
967            igrdi(2,2,jobs) = 2
968            igrdj(2,2,jobs) = 2
969
970         ELSE
971
972            igrdi(1,1,jobs) = kobsi(jobs)-1
973            igrdj(1,1,jobs) = kobsj(jobs)-1
974            igrdi(1,2,jobs) = kobsi(jobs)-1
975            igrdj(1,2,jobs) = kobsj(jobs)
976            igrdi(2,1,jobs) = kobsi(jobs)
977            igrdj(2,1,jobs) = kobsj(jobs)-1
978            igrdi(2,2,jobs) = kobsi(jobs)
979            igrdj(2,2,jobs) = kobsj(jobs)
980
981         ENDIF
982
983      END DO
984
985      IF (ln_bdy) THEN
986        ! Create a mask grid points in boundary rim
987        IF (ld_bound_reject) THEN
988           zbdymask(:,:) = 1.0_wp
989           DO ji = 1, nb_bdy
990              DO jj = 1, idx_bdy(ji)%nblen(1)
991                 zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp
992              ENDDO
993           ENDDO
994
995           CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, zbdymask, zbmsk )
996        ENDIF
997      ENDIF
998
999     
1000      CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, pmask, zgmsk )
1001      CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, plam, zglam )
1002      CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, pphi, zgphi )
1003
1004      DO jobs = 1, kobsno
1005
1006         ! Skip bad observations
1007         IF ( kobsqc(jobs) >= kqc_cutoff ) CYCLE
1008
1009         ! Flag if the observation falls outside the model spatial domain
1010         IF (       ( pobslam(jobs) < -180. ) &
1011            &  .OR. ( pobslam(jobs) >  180. ) &
1012            &  .OR. ( pobsphi(jobs) <  -90. ) &
1013            &  .OR. ( pobsphi(jobs) >   90. ) ) THEN
1014            kobsqc(jobs) = IBSET(kobsqc(jobs),11)
1015            kosdobs = kosdobs + 1
1016            CYCLE
1017         ENDIF
1018
1019         ! Flag if the observation falls with a model land cell
1020         IF ( SUM( zgmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN
1021            kobsqc(jobs) = IBSET(kobsqc(jobs),10)
1022            klanobs = klanobs + 1
1023            CYCLE
1024         ENDIF
1025
1026         ! Check if this observation is on a grid point
1027
1028         lgridobs = .FALSE.
1029         iig = -1
1030         ijg = -1
1031         DO jj = 1, 2
1032            DO ji = 1, 2
1033               IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) &
1034                  & .AND. &
1035                  & ( ABS( MOD( zglam(ji,jj,jobs) - pobslam(jobs),360.0) )  &
1036                  & < 1.0e-6_wp ) ) THEN
1037                  lgridobs = .TRUE.
1038                  iig = ji
1039                  ijg = jj
1040               ENDIF
1041            END DO
1042         END DO
1043 
1044         ! For observations on the grid reject them if their are at
1045         ! a masked point
1046         
1047         IF (lgridobs) THEN
1048            IF (zgmsk(iig,ijg,jobs) == 0.0_wp ) THEN
1049               kobsqc(jobs) = IBSET(kobsqc(jobs),10)
1050               klanobs = klanobs + 1
1051               CYCLE
1052            ENDIF
1053         ENDIF
1054                     
1055         ! Flag if the observation falls is close to land
1056         IF ( MINVAL( zgmsk(1:2,1:2,jobs) ) == 0.0_wp) THEN
1057            knlaobs = knlaobs + 1
1058            IF (ld_nea) THEN
1059               kobsqc(jobs) = IBSET(kobsqc(jobs),9)
1060               CYCLE
1061            ENDIF
1062         ENDIF
1063
1064         IF (ln_bdy) THEN
1065         ! Flag if the observation falls close to the boundary rim
1066           IF (ld_bound_reject) THEN
1067              IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN
1068                 kobsqc(jobs) = IBSET(kobsqc(jobs),8)
1069                 kbdyobs = kbdyobs + 1
1070                 CYCLE
1071              ENDIF
1072              ! for observations on the grid...
1073              IF (lgridobs) THEN
1074                 IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN
1075                    kobsqc(jobs) = IBSET(kobsqc(jobs),8)
1076                    kbdyobs = kbdyobs + 1
1077                    CYCLE
1078                 ENDIF
1079              ENDIF
1080            ENDIF
1081         ENDIF
1082         !
1083      END DO
1084      !
1085   END SUBROUTINE obs_coo_spc_2d
1086
1087
1088   SUBROUTINE obs_coo_spc_3d( kprofno, kobsno,  kpstart, kpend, &
1089      &                       kpi,     kpj,     kpk,            &
1090      &                       kobsi,   kobsj,   kobsk,          &
1091      &                       pobslam, pobsphi, pobsdep,        &
1092      &                       plam,    pphi,    pdep,    pmask, &
1093      &                       kpobsqc, kobsqc,  kosdobs,        &
1094      &                       klanobs, knlaobs, ld_nea,         &
1095      &                       kbdyobs, ld_bound_reject,         &
1096      &                       kqc_cutoff                        )
1097      !!----------------------------------------------------------------------
1098      !!                    ***  ROUTINE obs_coo_spc_3d  ***
1099      !!
1100      !! ** Purpose : Check for points outside the domain and land points
1101      !!              Reset depth of observation above highest model level
1102      !!              to the value of highest model level
1103      !!
1104      !! ** Method  : Remove the observations that are outside the model space
1105      !!              and time domain or located within model land cells.
1106      !!
1107      !!              NB. T and S profile observations lying between the ocean
1108      !!              surface and the depth of the first model T point are
1109      !!              assigned a depth equal to that of the first model T pt.
1110      !!
1111      !! ** Action  :
1112      !!   
1113      !! History :
1114      !!        !  2007-01  (K. Mogensen) Rewrite of parts of obs_scr
1115      !!        !  2007-06  (K. Mogensen et al) Reject obs. near land.
1116      !!----------------------------------------------------------------------
1117      !! * Modules used
1118      USE dom_oce, ONLY : &       ! Geographical information
1119         & gdepw_1d,      &
1120         & gdepw_0,       &                       
1121         & gdepw_n,       &
1122         & gdept_n,       &
1123         & ln_zco,        &
1124         & ln_zps             
1125
1126      !! * Arguments
1127      INTEGER, INTENT(IN) :: kprofno      ! Number of profiles
1128      INTEGER, INTENT(IN) :: kobsno       ! Total number of observations
1129      INTEGER, INTENT(IN) :: kpi          ! Number of grid points in (i,j,k)
1130      INTEGER, INTENT(IN) :: kpj
1131      INTEGER, INTENT(IN) :: kpk   
1132      INTEGER, DIMENSION(kprofno), INTENT(IN) :: &
1133         & kpstart, &         ! Start of individual profiles
1134         & kpend              ! End of individual profiles
1135      INTEGER, DIMENSION(kprofno), INTENT(IN) :: &
1136         & kobsi, &           ! Observation (i,j) coordinates
1137         & kobsj
1138      INTEGER, DIMENSION(kobsno), INTENT(IN) :: &
1139         & kobsk              ! Observation k coordinate
1140      REAL(KIND=wp), DIMENSION(kprofno), INTENT(IN) :: &
1141         & pobslam, &         ! Observation (lon,lat) coordinates
1142         & pobsphi
1143      REAL(KIND=wp), DIMENSION(kobsno), INTENT(INOUT) :: &
1144         & pobsdep            ! Observation depths 
1145      REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) :: &
1146         & plam, pphi         ! Model (lon,lat) coordinates
1147      REAL(KIND=wp), DIMENSION(kpk), INTENT(IN) :: &
1148         & pdep               ! Model depth coordinates
1149      REAL(KIND=wp), DIMENSION(kpi,kpj,kpk), INTENT(IN) :: &
1150         & pmask              ! Land mask array
1151      INTEGER, DIMENSION(kprofno), INTENT(INOUT) :: &
1152         & kpobsqc             ! Profile quality control
1153      INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
1154         & kobsqc             ! Observation quality control
1155      INTEGER, INTENT(INOUT) :: kosdobs     ! Observations outside space domain
1156      INTEGER, INTENT(INOUT) :: klanobs     ! Observations within a model land cell
1157      INTEGER, INTENT(INOUT) :: knlaobs     ! Observations near land
1158      INTEGER, INTENT(INOUT) :: kbdyobs     ! Observations near boundary
1159      LOGICAL, INTENT(IN) :: ld_nea         ! Flag observations near land
1160      LOGICAL, INTENT(IN) :: ld_bound_reject  ! Flag observations near open boundary
1161      INTEGER, INTENT(IN) :: kqc_cutoff     ! Cutoff QC value
1162
1163      !! * Local declarations
1164      REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: &
1165         & zgmsk              ! Grid mask
1166      REAL(KIND=wp), DIMENSION(2,2,kprofno) :: &
1167         & zbmsk              ! Boundary mask
1168      REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask
1169      REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: &
1170         & zgdepw
1171      REAL(KIND=wp), DIMENSION(2,2,kprofno) :: &
1172         & zglam, &           ! Model longitude at grid points
1173         & zgphi              ! Model latitude at grid points
1174      INTEGER, DIMENSION(2,2,kprofno) :: &
1175         & igrdi, &           ! Grid i,j
1176         & igrdj
1177      LOGICAL :: lgridobs           ! Is observation on a model grid point.
1178      LOGICAL :: ll_next_to_land    ! Is a profile next to land
1179      INTEGER :: iig, ijg           ! i,j of observation on model grid point.
1180      INTEGER :: jobs, jobsp, jk, ji, jj
1181      !!----------------------------------------------------------------------
1182
1183      ! Get grid point indices
1184     
1185      DO jobs = 1, kprofno
1186
1187         ! For invalid points use 2,2
1188
1189         IF ( kpobsqc(jobs) >= kqc_cutoff ) THEN
1190           
1191            igrdi(1,1,jobs) = 1
1192            igrdj(1,1,jobs) = 1
1193            igrdi(1,2,jobs) = 1
1194            igrdj(1,2,jobs) = 2
1195            igrdi(2,1,jobs) = 2
1196            igrdj(2,1,jobs) = 1
1197            igrdi(2,2,jobs) = 2
1198            igrdj(2,2,jobs) = 2
1199           
1200         ELSE
1201           
1202            igrdi(1,1,jobs) = kobsi(jobs)-1
1203            igrdj(1,1,jobs) = kobsj(jobs)-1
1204            igrdi(1,2,jobs) = kobsi(jobs)-1
1205            igrdj(1,2,jobs) = kobsj(jobs)
1206            igrdi(2,1,jobs) = kobsi(jobs)
1207            igrdj(2,1,jobs) = kobsj(jobs)-1
1208            igrdi(2,2,jobs) = kobsi(jobs)
1209            igrdj(2,2,jobs) = kobsj(jobs)
1210           
1211         ENDIF
1212         
1213      END DO
1214
1215      IF (ln_bdy) THEN
1216        ! Create a mask grid points in boundary rim
1217        IF (ld_bound_reject) THEN           
1218           zbdymask(:,:) = 1.0_wp
1219           DO ji = 1, nb_bdy
1220              DO jj = 1, idx_bdy(ji)%nblen(1)
1221                 zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp
1222              ENDDO
1223           ENDDO
1224        ENDIF
1225   
1226        CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, zbdymask, zbmsk )
1227      ENDIF
1228     
1229      CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, pmask, zgmsk )
1230      CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, plam, zglam )
1231      CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, pphi, zgphi )
1232      CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, gdepw_n(:,:,:), &
1233        &                     zgdepw )
1234
1235      DO jobs = 1, kprofno
1236
1237         ! Skip bad profiles
1238         IF ( kpobsqc(jobs) >= kqc_cutoff ) CYCLE
1239
1240         ! Check if this observation is on a grid point
1241
1242         lgridobs = .FALSE.
1243         iig = -1
1244         ijg = -1
1245         DO jj = 1, 2
1246            DO ji = 1, 2
1247               IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) &
1248                  & .AND. &
1249                  & ( ABS( MOD( zglam(ji,jj,jobs) - pobslam(jobs),360.0) ) < 1.0e-6_wp ) &
1250                  & ) THEN
1251                  lgridobs = .TRUE.
1252                  iig = ji
1253                  ijg = jj
1254               ENDIF
1255            END DO
1256         END DO
1257
1258         ! Check if next to land
1259         IF (  ANY( zgmsk(1:2,1:2,1,jobs) == 0.0_wp ) ) THEN
1260            ll_next_to_land=.TRUE.
1261         ELSE
1262            ll_next_to_land=.FALSE.
1263         ENDIF 
1264
1265         ! Reject observations
1266
1267         DO jobsp = kpstart(jobs), kpend(jobs)
1268
1269            ! Flag if the observation falls outside the model spatial domain
1270            IF (       ( pobslam(jobs) < -180.         )       &
1271               &  .OR. ( pobslam(jobs) >  180.         )       &
1272               &  .OR. ( pobsphi(jobs) <  -90.         )       &
1273               &  .OR. ( pobsphi(jobs) >   90.         )       &
1274               &  .OR. ( pobsdep(jobsp) < 0.0          )       &
1275               &  .OR. ( pobsdep(jobsp) > gdepw_1d(kpk)) ) THEN
1276               kobsqc(jobsp) = IBSET(kobsqc(jobsp),11)
1277               kosdobs = kosdobs + 1
1278               CYCLE
1279            ENDIF
1280
1281            ! To check if an observations falls within land:
1282             
1283            ! Flag if the observation is deeper than the bathymetry
1284            ! Or if it is within the mask
1285            IF ( ALL( zgdepw(1:2,1:2,kpk,jobs) < pobsdep(jobsp) ) &
1286               &     .OR. &
1287               &  ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) &
1288               &  == 0.0_wp) ) THEN
1289               kobsqc(jobsp) = IBSET(kobsqc(jobsp),10)
1290               klanobs = klanobs + 1
1291               CYCLE
1292            ENDIF
1293               
1294            ! Flag if the observation is close to land
1295            IF ( ll_next_to_land ) THEN
1296               knlaobs = knlaobs + 1
1297               IF (ld_nea) THEN   
1298                  kobsqc(jobsp) = IBSET(kobsqc(jobsp),10)
1299               ENDIF
1300            ENDIF
1301           
1302            ! For observations on the grid reject them if their are at
1303            ! a masked point
1304           
1305            IF (lgridobs) THEN
1306               IF (zgmsk(iig,ijg,kobsk(jobsp)-1,jobs) == 0.0_wp ) THEN
1307                  kobsqc(jobsp) = IBSET(kobsqc(jobs),10)
1308                  klanobs = klanobs + 1
1309                  CYCLE
1310               ENDIF
1311            ENDIF
1312           
1313            ! Flag if the observation falls is close to land
1314            IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == &
1315               &  0.0_wp) THEN
1316               IF (ld_nea) kobsqc(jobsp) = kobsqc(jobsp) + 14
1317               knlaobs = knlaobs + 1
1318            ENDIF
1319
1320            ! Set observation depth equal to that of the first model depth
1321            IF ( pobsdep(jobsp) <= pdep(1) ) THEN
1322               pobsdep(jobsp) = pdep(1) 
1323            ENDIF
1324           
1325            IF (ln_bdy) THEN
1326               ! Flag if the observation falls close to the boundary rim
1327               IF (ld_bound_reject) THEN
1328                  IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN
1329                     kobsqc(jobsp) = IBSET(kobsqc(jobs),8)
1330                     kbdyobs = kbdyobs + 1
1331                     CYCLE
1332                  ENDIF
1333                  ! for observations on the grid...
1334                  IF (lgridobs) THEN
1335                     IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN
1336                        kobsqc(jobsp) = IBSET(kobsqc(jobs),8)
1337                        kbdyobs = kbdyobs + 1
1338                        CYCLE
1339                     ENDIF
1340                  ENDIF
1341               ENDIF
1342            ENDIF
1343            !
1344         END DO
1345      END DO
1346      !
1347   END SUBROUTINE obs_coo_spc_3d
1348
1349
1350   SUBROUTINE obs_pro_rej( profdata, kqc_cutoff )
1351      !!----------------------------------------------------------------------
1352      !!                    ***  ROUTINE obs_pro_rej ***
1353      !!
1354      !! ** Purpose : Reject all data within a rejected profile
1355      !!
1356      !! ** Method  :
1357      !!
1358      !! ** Action  :
1359      !!
1360      !! References :
1361      !!   
1362      !! History :   2007-10  (K. Mogensen) Original code
1363      !!----------------------------------------------------------------------
1364      TYPE(obs_prof), INTENT(inout) ::   profdata     ! Profile data
1365      INTEGER       , INTENT(in   ) ::   kqc_cutoff   ! QC cutoff value
1366      !
1367      INTEGER :: jprof
1368      INTEGER :: jvar
1369      INTEGER :: jobs
1370      !!----------------------------------------------------------------------
1371     
1372      ! Loop over profiles
1373
1374      DO jprof = 1, profdata%nprof
1375
1376         IF ( profdata%nqc(jprof) > kqc_cutoff ) THEN
1377           
1378            DO jvar = 1, profdata%nvar
1379
1380               DO jobs = profdata%npvsta(jprof,jvar),  &
1381                  &      profdata%npvend(jprof,jvar)
1382                 
1383                  profdata%var(jvar)%nvqc(jobs) = &
1384                     & IBSET(profdata%var(jvar)%nvqc(jobs),14)
1385
1386               END DO
1387
1388            END DO
1389
1390         ENDIF
1391
1392      END DO
1393      !
1394   END SUBROUTINE obs_pro_rej
1395
1396
1397   SUBROUTINE obs_uv_rej( profdata, knumu, knumv, kqc_cutoff )
1398      !!----------------------------------------------------------------------
1399      !!                    ***  ROUTINE obs_uv_rej ***
1400      !!
1401      !! ** Purpose : Reject u if v is rejected and vice versa
1402      !!
1403      !! ** Method  :
1404      !!
1405      !! ** Action  :
1406      !!
1407      !! References :
1408      !!   
1409      !! History :   2009-2  (K. Mogensen) Original code
1410      !!----------------------------------------------------------------------
1411      TYPE(obs_prof), INTENT(INOUT) :: profdata   ! Profile data
1412      INTEGER, INTENT(INOUT) :: knumu             ! Number of u rejected
1413      INTEGER, INTENT(INOUT) :: knumv             ! Number of v rejected
1414      INTEGER, INTENT(IN) :: kqc_cutoff           ! QC cutoff value
1415      !
1416      INTEGER :: jprof
1417      INTEGER :: jvar
1418      INTEGER :: jobs
1419      !!----------------------------------------------------------------------
1420
1421      DO jprof = 1, profdata%nprof      !==  Loop over profiles  ==!
1422         !
1423         IF ( ( profdata%npvsta(jprof,1) /= profdata%npvsta(jprof,2) ) .OR. &
1424            & ( profdata%npvend(jprof,1) /= profdata%npvend(jprof,2) ) ) THEN
1425            !
1426            CALL ctl_stop('U,V profiles inconsistent in obs_uv_rej')
1427            RETURN
1428            !
1429         ENDIF
1430         !
1431         DO jobs = profdata%npvsta(jprof,1), profdata%npvend(jprof,1)
1432           
1433            IF ( ( profdata%var(1)%nvqc(jobs) >  kqc_cutoff ) .AND. &
1434               & ( profdata%var(2)%nvqc(jobs) <=  kqc_cutoff) ) THEN
1435               profdata%var(2)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15)
1436               knumv = knumv + 1
1437            ENDIF
1438            IF ( ( profdata%var(2)%nvqc(jobs) >  kqc_cutoff ) .AND. &
1439               & ( profdata%var(1)%nvqc(jobs) <=  kqc_cutoff) ) THEN
1440               profdata%var(1)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15)
1441               knumu = knumu + 1
1442            ENDIF
1443            !
1444         END DO
1445         !
1446      END DO
1447      !
1448   END SUBROUTINE obs_uv_rej
1449
1450   !!=====================================================================
1451END MODULE obs_prep
Note: See TracBrowser for help on using the repository browser.