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

source: branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 8035

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