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

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

source: NEMO/branches/UKMO/NEMO_4.0.4_generic_obs/src/OCE/OBS/obs_prep.F90 @ 15228

Last change on this file since 15228 was 15228, checked in by dford, 3 years ago

Rename module.

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