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

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

source: NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 11350

Last change on this file since 11350 was 11350, checked in by jcastill, 5 years ago

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