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

source: branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 7837

Last change on this file since 7837 was 7837, checked in by mattmartin, 7 years ago

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