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

source: branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 5785

Last change on this file since 5785 was 5785, checked in by mattmartin, 8 years ago

Updated obs_prep.F90 based on reviewers comment about the call to obs_coo_grd. This change was tested with SETTE.

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