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

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

source: branches/UKMO/dev_r6393_CO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 7019

Last change on this file since 7019 was 7019, checked in by deazer, 7 years ago

Cleared svn keywords

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