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

source: branches/UKMO/dev_r5518_obs_oper_update_sit/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 12557

Last change on this file since 12557 was 12557, checked in by dcarneir, 4 years ago

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