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

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

source: branches/UKMO/dev_r5518_obs_oper_update_amm7ps45/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 15430

Last change on this file since 15430 was 15430, checked in by petesykes, 12 months ago

merging vertical interp bug fix from parent branch

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