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 @ 15487

Last change on this file since 15487 was 15487, checked in by dford, 8 months ago

Remove some comments and add some spaces.

File size: 61.9 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         & mbkt
1124
1125      !! * Arguments
1126      INTEGER, INTENT(IN) :: kprofno      ! Number of profiles
1127      INTEGER, INTENT(IN) :: kobsno       ! Total number of observations
1128      INTEGER, INTENT(IN) :: kpi          ! Number of grid points in (i,j,k)
1129      INTEGER, INTENT(IN) :: kpj
1130      INTEGER, INTENT(IN) :: kpk   
1131      INTEGER, DIMENSION(kprofno), INTENT(IN) :: &
1132         & kpstart, &         ! Start of individual profiles
1133         & kpend              ! End of individual profiles
1134      INTEGER, DIMENSION(kprofno), INTENT(IN) :: &
1135         & kobsi, &           ! Observation (i,j) coordinates
1136         & kobsj
1137      INTEGER, DIMENSION(kobsno), INTENT(IN) :: &
1138         & kobsk              ! Observation k coordinate
1139      REAL(KIND=wp), DIMENSION(kprofno), INTENT(IN) :: &
1140         & pobslam, &         ! Observation (lon,lat) coordinates
1141         & pobsphi
1142      REAL(KIND=wp), DIMENSION(kobsno), INTENT(INOUT) :: &
1143         & pobsdep            ! Observation depths 
1144      REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) :: &
1145         & plam, pphi         ! Model (lon,lat) coordinates
1146      REAL(KIND=wp), DIMENSION(kpk), INTENT(IN) :: &
1147         & pdep               ! Model depth coordinates
1148      REAL(KIND=wp), DIMENSION(kpi,kpj,kpk), INTENT(IN) :: &
1149         & pmask              ! Land mask array
1150      INTEGER, DIMENSION(kprofno), INTENT(INOUT) :: &
1151         & kpobsqc             ! Profile quality control
1152      INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
1153         & kobsqc             ! Observation quality control
1154      INTEGER, INTENT(INOUT) :: kosdobs     ! Observations outside space domain
1155      INTEGER, INTENT(INOUT) :: klanobs     ! Observations within a model land cell
1156      INTEGER, INTENT(INOUT) :: knlaobs     ! Observations near land
1157      INTEGER, INTENT(INOUT) :: kbdyobs     ! Observations near boundary
1158      LOGICAL, INTENT(IN) :: ld_nea         ! Flag observations near land
1159      LOGICAL, INTENT(IN) :: ld_bound_reject  ! Flag observations near open boundary
1160      INTEGER, INTENT(IN) :: kqc_cutoff     ! Cutoff QC value
1161
1162      !! * Local declarations
1163      REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: &
1164         & zgmsk              ! Grid mask
1165      REAL(KIND=wp), DIMENSION(2,2,kprofno) :: &
1166         & zbmsk              ! Boundary mask
1167      REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask
1168      REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: &
1169         & zgdept, &
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         & zbathy             ! Index of deepest wet level at grid points
1175      INTEGER, DIMENSION(2,2,kprofno) :: &
1176         & igrdi, &           ! Grid i,j
1177         & igrdj
1178      LOGICAL :: lgridobs           ! Is observation on a model grid point.
1179      LOGICAL :: ll_next_to_land    ! Is a profile next to land
1180      INTEGER :: iig, ijg           ! i,j of observation on model grid point.
1181      INTEGER :: jobs, jobsp, jk, ji, jj
1182      REAL(KIND=wp) :: maxdept, maxdepw
1183      !!----------------------------------------------------------------------
1184
1185      ! Get grid point indices
1186     
1187      DO jobs = 1, kprofno
1188
1189         ! For invalid points use 2,2
1190
1191         IF ( kpobsqc(jobs) >= kqc_cutoff ) THEN
1192           
1193            igrdi(1,1,jobs) = 1
1194            igrdj(1,1,jobs) = 1
1195            igrdi(1,2,jobs) = 1
1196            igrdj(1,2,jobs) = 2
1197            igrdi(2,1,jobs) = 2
1198            igrdj(2,1,jobs) = 1
1199            igrdi(2,2,jobs) = 2
1200            igrdj(2,2,jobs) = 2
1201           
1202         ELSE
1203           
1204            igrdi(1,1,jobs) = kobsi(jobs)-1
1205            igrdj(1,1,jobs) = kobsj(jobs)-1
1206            igrdi(1,2,jobs) = kobsi(jobs)-1
1207            igrdj(1,2,jobs) = kobsj(jobs)
1208            igrdi(2,1,jobs) = kobsi(jobs)
1209            igrdj(2,1,jobs) = kobsj(jobs)-1
1210            igrdi(2,2,jobs) = kobsi(jobs)
1211            igrdj(2,2,jobs) = kobsj(jobs)
1212           
1213         ENDIF
1214         
1215      END DO
1216
1217      IF (ln_bdy) THEN
1218        ! Create a mask grid points in boundary rim
1219        IF (ld_bound_reject) THEN           
1220           zbdymask(:,:) = 1.0_wp
1221           DO ji = 1, nb_bdy
1222              DO jj = 1, idx_bdy(ji)%nblen(1)
1223                 zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp
1224              ENDDO
1225           ENDDO
1226        ENDIF
1227   
1228        CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, zbdymask, zbmsk )
1229      ENDIF
1230     
1231      CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, pmask, zgmsk )
1232      CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, plam, zglam )
1233      CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, pphi, zgphi )
1234      CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, REAL(mbkt), zbathy )
1235      ! Need to know the bathy depth for each observation for sco
1236      CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, gdepw_n(:,:,:), &
1237        &                   zgdepw )
1238      CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, gdept_n(:,:,:), &
1239         &                  zgdept )
1240
1241      DO jobs = 1, kprofno
1242
1243         ! Skip bad profiles
1244         IF ( kpobsqc(jobs) >= kqc_cutoff ) CYCLE
1245
1246         ! Check if this observation is on a grid point
1247
1248         lgridobs = .FALSE.
1249         iig = -1
1250         ijg = -1
1251         DO jj = 1, 2
1252            DO ji = 1, 2
1253               IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) &
1254                  & .AND. &
1255                  & ( ABS( MOD( zglam(ji,jj,jobs) - pobslam(jobs),360.0) ) < 1.0e-6_wp ) &
1256                  & ) THEN
1257                  lgridobs = .TRUE.
1258                  iig = ji
1259                  ijg = jj
1260               ENDIF
1261            END DO
1262         END DO
1263
1264         ! Check if next to land
1265         IF (  ANY( zgmsk(1:2,1:2,1,jobs) == 0.0_wp ) ) THEN
1266            ll_next_to_land=.TRUE.
1267         ELSE
1268            ll_next_to_land=.FALSE.
1269         ENDIF 
1270
1271         ! Reject observations
1272
1273         DO jobsp = kpstart(jobs), kpend(jobs)
1274
1275            ! Calculate max T and W depths of 2x2 grid
1276            maxdept = zgdept(1,1,NINT(zbathy(1,1,jobs)),jobs)
1277            maxdepw = zgdepw(1,1,NINT(zbathy(1,1,jobs))+1,jobs)
1278            DO jj = 1, 2
1279               DO ji = 1, 2
1280                  IF ( zgdept(ji,jj,NINT(zbathy(ji,jj,jobs)),jobs) > maxdept ) THEN
1281                     maxdept = zgdept(ji,jj,NINT(zbathy(ji,jj,jobs)),jobs)
1282                  END IF
1283                  IF ( zgdepw(ji,jj,NINT(zbathy(ji,jj,jobs))+1,jobs) > maxdepw ) THEN
1284                     maxdepw = zgdepw(ji,jj,NINT(zbathy(ji,jj,jobs))+1,jobs)
1285                  END IF
1286               END DO
1287            END DO
1288
1289            ! Flag if the observation falls outside the model spatial domain
1290            IF (       ( pobslam(jobs) < -180.    )       &
1291               &  .OR. ( pobslam(jobs) >  180.    )       &
1292               &  .OR. ( pobsphi(jobs) <  -90.    )       &
1293               &  .OR. ( pobsphi(jobs) >   90.    )       &
1294               &  .OR. ( pobsdep(jobsp) < 0.0     )       &
1295               &  .OR. ( pobsdep(jobsp) > maxdepw ) ) THEN
1296               kobsqc(jobsp) = IBSET(kobsqc(jobsp),11)
1297               kosdobs = kosdobs + 1
1298               CYCLE
1299            ENDIF
1300
1301            ! To check if an observations falls within land:
1302             
1303            ! Flag if the observation is deeper than the bathymetry
1304            ! Or if it is within the mask
1305            IF ( ALL( zgdepw(1:2,1:2,kpk,jobs) < pobsdep(jobsp) ) &
1306               &     .OR. &
1307               &  ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) &
1308               &  == 0.0_wp) ) THEN
1309               kobsqc(jobsp) = IBSET(kobsqc(jobsp),10)
1310               klanobs = klanobs + 1
1311               CYCLE
1312            ENDIF
1313               
1314            ! Flag if the observation is close to land
1315            IF ( ll_next_to_land ) THEN
1316               knlaobs = knlaobs + 1
1317               IF (ld_nea) THEN   
1318                  kobsqc(jobsp) = IBSET(kobsqc(jobsp),10)
1319               ENDIF
1320            ENDIF
1321           
1322            ! For observations on the grid reject them if their are at
1323            ! a masked point
1324           
1325            IF (lgridobs) THEN
1326               IF (zgmsk(iig,ijg,kobsk(jobsp)-1,jobs) == 0.0_wp ) THEN
1327                  kobsqc(jobsp) = IBSET(kobsqc(jobs),10)
1328                  klanobs = klanobs + 1
1329                  CYCLE
1330               ENDIF
1331            ENDIF
1332           
1333            ! Flag if the observation falls is close to land
1334            IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == &
1335               &  0.0_wp) THEN
1336               IF (ld_nea) kobsqc(jobsp) = kobsqc(jobsp) + 14
1337               knlaobs = knlaobs + 1
1338            ENDIF
1339
1340            ! Set observation depth equal to that of the first model depth
1341            IF ( pobsdep(jobsp) < MINVAL(zgdept(1:2,1:2,1,jobs) ) ) THEN
1342               pobsdep(jobsp) = MINVAL(zgdept(1:2,1:2,1,jobs))
1343            ENDIF
1344
1345            ! Set observation depth equal to that of the last wet T-point
1346            IF ( ( pobsdep(jobsp) > maxdept ) .AND. &
1347               & ( pobsdep(jobsp) < maxdepw ) ) THEN
1348               pobsdep(jobsp) = maxdept
1349            END IF
1350           
1351            IF (ln_bdy) THEN
1352               ! Flag if the observation falls close to the boundary rim
1353               IF (ld_bound_reject) THEN
1354                  IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN
1355                     kobsqc(jobsp) = IBSET(kobsqc(jobs),8)
1356                     kbdyobs = kbdyobs + 1
1357                     CYCLE
1358                  ENDIF
1359                  ! for observations on the grid...
1360                  IF (lgridobs) THEN
1361                     IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN
1362                        kobsqc(jobsp) = IBSET(kobsqc(jobs),8)
1363                        kbdyobs = kbdyobs + 1
1364                        CYCLE
1365                     ENDIF
1366                  ENDIF
1367               ENDIF
1368            ENDIF
1369            !
1370         END DO
1371      END DO
1372      !
1373   END SUBROUTINE obs_coo_spc_3d
1374
1375
1376   SUBROUTINE obs_pro_rej( profdata, kqc_cutoff )
1377      !!----------------------------------------------------------------------
1378      !!                    ***  ROUTINE obs_pro_rej ***
1379      !!
1380      !! ** Purpose : Reject all data within a rejected profile
1381      !!
1382      !! ** Method  :
1383      !!
1384      !! ** Action  :
1385      !!
1386      !! References :
1387      !!   
1388      !! History :   2007-10  (K. Mogensen) Original code
1389      !!----------------------------------------------------------------------
1390      TYPE(obs_prof), INTENT(inout) ::   profdata     ! Profile data
1391      INTEGER       , INTENT(in   ) ::   kqc_cutoff   ! QC cutoff value
1392      !
1393      INTEGER :: jprof
1394      INTEGER :: jvar
1395      INTEGER :: jobs
1396      !!----------------------------------------------------------------------
1397     
1398      ! Loop over profiles
1399
1400      DO jprof = 1, profdata%nprof
1401
1402         IF ( profdata%nqc(jprof) > kqc_cutoff ) THEN
1403           
1404            DO jvar = 1, profdata%nvar
1405
1406               DO jobs = profdata%npvsta(jprof,jvar),  &
1407                  &      profdata%npvend(jprof,jvar)
1408                 
1409                  profdata%var(jvar)%nvqc(jobs) = &
1410                     & IBSET(profdata%var(jvar)%nvqc(jobs),14)
1411
1412               END DO
1413
1414            END DO
1415
1416         ENDIF
1417
1418      END DO
1419      !
1420   END SUBROUTINE obs_pro_rej
1421
1422
1423   SUBROUTINE obs_uv_rej( profdata, knumu, knumv, kqc_cutoff, kuvar, kvvar )
1424      !!----------------------------------------------------------------------
1425      !!                    ***  ROUTINE obs_uv_rej ***
1426      !!
1427      !! ** Purpose : Reject u if v is rejected and vice versa
1428      !!
1429      !! ** Method  :
1430      !!
1431      !! ** Action  :
1432      !!
1433      !! References :
1434      !!   
1435      !! History :   2009-2  (K. Mogensen) Original code
1436      !!----------------------------------------------------------------------
1437      TYPE(obs_prof), INTENT(INOUT) :: profdata   ! Profile data
1438      INTEGER, INTENT(INOUT) :: knumu             ! Number of u rejected
1439      INTEGER, INTENT(INOUT) :: knumv             ! Number of v rejected
1440      INTEGER, INTENT(IN)    :: kqc_cutoff        ! QC cutoff value
1441      INTEGER, INTENT(IN)    :: kuvar             ! Index of u
1442      INTEGER, INTENT(IN)    :: kvvar             ! Index of v
1443      !
1444      INTEGER :: jprof
1445      INTEGER :: jvar
1446      INTEGER :: jobs
1447      !!----------------------------------------------------------------------
1448
1449      DO jprof = 1, profdata%nprof      !==  Loop over profiles  ==!
1450         !
1451         IF ( ( profdata%npvsta(jprof,kuvar) /= profdata%npvsta(jprof,kvvar) ) .OR. &
1452            & ( profdata%npvend(jprof,kuvar) /= profdata%npvend(jprof,kvvar) ) ) THEN
1453            !
1454            CALL ctl_stop('U,V profiles inconsistent in obs_uv_rej')
1455            RETURN
1456            !
1457         ENDIF
1458         !
1459         DO jobs = profdata%npvsta(jprof,kuvar), profdata%npvend(jprof,kuvar)
1460           
1461            IF ( ( profdata%var(kuvar)%nvqc(jobs) >  kqc_cutoff ) .AND. &
1462               & ( profdata%var(kvvar)%nvqc(jobs) <=  kqc_cutoff) ) THEN
1463               profdata%var(kvvar)%nvqc(jobs) = IBSET(profdata%var(kuvar)%nvqc(jobs),15)
1464               knumv = knumv + 1
1465            ENDIF
1466            IF ( ( profdata%var(kvvar)%nvqc(jobs) >  kqc_cutoff ) .AND. &
1467               & ( profdata%var(kuvar)%nvqc(jobs) <=  kqc_cutoff) ) THEN
1468               profdata%var(kuvar)%nvqc(jobs) = IBSET(profdata%var(kuvar)%nvqc(jobs),15)
1469               knumu = knumu + 1
1470            ENDIF
1471            !
1472         END DO
1473         !
1474      END DO
1475      !
1476   END SUBROUTINE obs_uv_rej
1477
1478   !!=====================================================================
1479END MODULE obs_prep
Note: See TracBrowser for help on using the repository browser.