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

source: branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 6855

Last change on this file since 6855 was 6855, checked in by dford, 8 years ago

Initial implementation of observation operator for SPM.

File size: 97.6 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_pro  : First level check and screening of T/S profiles
10   !!   obs_pre_sla  : First level check and screening of SLA observations
11   !!   obs_pre_sst  : First level check and screening of SLA observations
12   !!   obs_pre_seaice : First level check and screening of sea ice observations
13   !!   obs_pre_vel  : First level check and screening of velocity obs.
14   !!   obs_pre_logchl : First level check and screening of logchl obs.
15   !!   obs_pre_spm  : First level check and screening of spm obs.
16   !!   obs_scr      : Basic screening of the observations
17   !!   obs_coo_tim  : Compute number of time steps to the observation time
18   !!   obs_sor      : Sort the observation arrays
19   !!---------------------------------------------------------------------
20   !! * Modules used
21   USE par_kind, ONLY : & ! Precision variables
22      & wp   
23   USE in_out_manager     ! I/O manager
24   USE obs_profiles_def   ! Definitions for storage arrays for profiles
25   USE obs_surf_def       ! Definitions for storage arrays for surface data
26   USE obs_mpp, ONLY : &  ! MPP support routines for observation diagnostics
27      & obs_mpp_sum_integer, &
28      & obs_mpp_sum_integers
29   USE obs_inter_sup      ! Interpolation support
30   USE obs_oper           ! Observation operators
31   USE lib_mpp, ONLY : &
32      & ctl_warn, ctl_stop
33
34   IMPLICIT NONE
35
36   !! * Routine accessibility
37   PRIVATE
38
39   PUBLIC &
40      & obs_pre_pro, &    ! First level check and screening of profiles
41      & obs_pre_sla, &    ! First level check and screening of SLA data
42      & obs_pre_sst, &    ! First level check and screening of SLA data
43      & obs_pre_seaice, & ! First level check and screening of sea ice data
44      & obs_pre_vel, &     ! First level check and screening of velocity profiles
45      & obs_pre_logchl, & ! First level check and screening of logchl data
46      & obs_pre_spm, &    ! First level check and screening of spm data
47      & calc_month_len     ! Calculate the number of days in the months of a year 
48
49   !!----------------------------------------------------------------------
50   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
51   !! $Id$
52   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
53   !!----------------------------------------------------------------------
54
55!! * Substitutions
56#  include "domzgr_substitute.h90" 
57
58CONTAINS
59
60   SUBROUTINE obs_pre_pro( profdata, prodatqc, ld_t3d, ld_s3d, ld_nea, &
61      &                    kdailyavtypes )
62      !!----------------------------------------------------------------------
63      !!                    ***  ROUTINE obs_pre_pro  ***
64      !!
65      !! ** Purpose : First level check and screening of T and S profiles
66      !!
67      !! ** Method  : First level check and screening of T and S profiles
68      !!
69      !! ** Action  :
70      !!
71      !! References :
72      !!   
73      !! History :
74      !!        !  2007-01  (K. Mogensen) Merge of obs_pre_t3d and obs_pre_s3d
75      !!        !  2007-03  (K. Mogensen) General handling of profiles
76      !!        !  2007-06  (K. Mogensen et al) Reject obs. near land.
77      !!----------------------------------------------------------------------
78      !! * Modules used
79      USE domstp              ! Domain: set the time-step
80      USE par_oce             ! Ocean parameters
81      USE dom_oce, ONLY : &   ! Geographical information
82         & glamt,   &
83         & gphit,   &
84         & gdept_1d,&
85         & tmask     
86      !! * Arguments
87      TYPE(obs_prof), INTENT(INOUT) :: profdata     ! Full set of profile data
88      TYPE(obs_prof), INTENT(INOUT) :: prodatqc     ! Subset of profile data not failing screening
89      LOGICAL, INTENT(IN) :: ld_t3d         ! Switch for temperature
90      LOGICAL, INTENT(IN) :: ld_s3d         ! Switch for salinity
91      LOGICAL, INTENT(IN) :: ld_nea         ! Switch for rejecting observation near land
92      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
93         & kdailyavtypes! Types for daily averages
94      !! * Local declarations   
95      INTEGER :: iyea0         ! Initial date
96      INTEGER :: imon0         !  - (year, month, day, hour, minute)
97      INTEGER :: iday0   
98      INTEGER :: ihou0
99      INTEGER :: imin0
100      INTEGER :: icycle        ! Current assimilation cycle
101                               ! Counters for observations that
102      INTEGER :: iotdobs       !  - outside time domain
103      INTEGER :: iosdtobs      !  - outside space domain (temperature)
104      INTEGER :: iosdsobs      !  - outside space domain (salinity)
105      INTEGER :: ilantobs      !  - within a model land cell (temperature)
106      INTEGER :: ilansobs      !  - within a model land cell (salinity)
107      INTEGER :: inlatobs      !  - close to land (temperature)
108      INTEGER :: inlasobs      !  - close to land (salinity)
109      INTEGER :: igrdobs       !  - fail the grid search
110                               ! Global counters for observations that
111      INTEGER :: iotdobsmpp    !  - outside time domain
112      INTEGER :: iosdtobsmpp   !  - outside space domain (temperature)
113      INTEGER :: iosdsobsmpp   !  - outside space domain (salinity)
114      INTEGER :: ilantobsmpp   !  - within a model land cell (temperature)
115      INTEGER :: ilansobsmpp   !  - within a model land cell (salinity)
116      INTEGER :: inlatobsmpp   !  - close to land (temperature)
117      INTEGER :: inlasobsmpp   !  - close to land (salinity)
118      INTEGER :: igrdobsmpp    !  - fail the grid search
119      TYPE(obs_prof_valid) ::  llvalid     ! Profile selection
120      TYPE(obs_prof_valid), DIMENSION(profdata%nvar) :: &
121         & llvvalid            ! T,S selection
122      INTEGER :: jvar          ! Variable loop variable
123      INTEGER :: jobs          ! Obs. loop variable
124      INTEGER :: jstp          ! Time loop variable
125      INTEGER :: inrc          ! Time index variable
126     
127      IF(lwp) WRITE(numout,*)'obs_pre_pro : Preparing the profile observations...'
128
129      ! Initial date initialization (year, month, day, hour, minute)
130      iyea0 =   ndate0 / 10000
131      imon0 = ( ndate0 - iyea0 * 10000 ) / 100
132      iday0 =   ndate0 - iyea0 * 10000 - imon0 * 100
133      ihou0 = 0
134      imin0 = 0
135
136      icycle = no     ! Assimilation cycle
137
138      ! Diagnotics counters for various failures.
139
140      iotdobs  = 0
141      igrdobs  = 0
142      iosdtobs = 0
143      iosdsobs = 0
144      ilantobs = 0
145      ilansobs = 0
146      inlatobs = 0
147      inlasobs = 0
148
149      ! -----------------------------------------------------------------------
150      ! Find time coordinate for profiles
151      ! -----------------------------------------------------------------------
152
153      IF ( PRESENT(kdailyavtypes) ) THEN
154         CALL obs_coo_tim_prof( icycle, &
155            &                iyea0,   imon0,   iday0,   ihou0,   imin0,      &
156            &                profdata%nprof,   profdata%nyea, profdata%nmon, &
157            &                profdata%nday,    profdata%nhou, profdata%nmin, &
158            &                profdata%ntyp,    profdata%nqc,  profdata%mstp, &
159            &                iotdobs, kdailyavtypes = kdailyavtypes        )
160      ELSE
161         CALL obs_coo_tim_prof( icycle, &
162            &                iyea0,   imon0,   iday0,   ihou0,   imin0,      &
163            &                profdata%nprof,   profdata%nyea, profdata%nmon, &
164            &                profdata%nday,    profdata%nhou, profdata%nmin, &
165            &                profdata%ntyp,    profdata%nqc,  profdata%mstp, &
166            &                iotdobs )
167      ENDIF
168      CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
169     
170      ! -----------------------------------------------------------------------
171      ! Check for profiles failing the grid search
172      ! -----------------------------------------------------------------------
173
174      CALL obs_coo_grd( profdata%nprof,   profdata%mi, profdata%mj, &
175         &              profdata%nqc,     igrdobs                         )
176
177      CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
178
179      ! -----------------------------------------------------------------------
180      ! Reject all observations for profiles with nqc > 10
181      ! -----------------------------------------------------------------------
182
183      CALL obs_pro_rej( profdata )
184
185      ! -----------------------------------------------------------------------
186      ! Check for land points. This includes points below the model
187      ! bathymetry so this is done for every point in the profile
188      ! -----------------------------------------------------------------------
189
190      ! Temperature
191
192      CALL obs_coo_spc_3d( profdata%nprof,        profdata%nvprot(1),   &
193         &                 profdata%npvsta(:,1),  profdata%npvend(:,1), &
194         &                 jpi,                   jpj,                  &
195         &                 jpk,                                         &
196         &                 profdata%mi,           profdata%mj,          & 
197         &                 profdata%var(1)%mvk,                         &
198         &                 profdata%rlam,         profdata%rphi,        &
199         &                 profdata%var(1)%vdep,                        &
200         &                 glamt,                 gphit,                &
201         &                 gdept_1d,              tmask,                &
202         &                 profdata%nqc,          profdata%var(1)%nvqc, &
203         &                 iosdtobs,              ilantobs,             &
204         &                 inlatobs,              ld_nea                )
205
206      CALL obs_mpp_sum_integer( iosdtobs, iosdtobsmpp )
207      CALL obs_mpp_sum_integer( ilantobs, ilantobsmpp )
208      CALL obs_mpp_sum_integer( inlatobs, inlatobsmpp )
209
210      ! Salinity
211
212      CALL obs_coo_spc_3d( profdata%nprof,        profdata%nvprot(2),   &
213         &                 profdata%npvsta(:,2),  profdata%npvend(:,2), &
214         &                 jpi,                   jpj,                  &
215         &                 jpk,                                         &
216         &                 profdata%mi,           profdata%mj,          & 
217         &                 profdata%var(2)%mvk,                         &
218         &                 profdata%rlam,         profdata%rphi,        &
219         &                 profdata%var(2)%vdep,                        &
220         &                 glamt,                 gphit,                &
221         &                 gdept_1d,              tmask,                &
222         &                 profdata%nqc,          profdata%var(2)%nvqc, &
223         &                 iosdsobs,              ilansobs,             &
224         &                 inlasobs,              ld_nea                )
225
226      CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
227      CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
228      CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
229
230      ! -----------------------------------------------------------------------
231      ! Copy useful data from the profdata data structure to
232      ! the prodatqc data structure
233      ! -----------------------------------------------------------------------
234
235      ! Allocate the selection arrays
236
237      ALLOCATE( llvalid%luse(profdata%nprof) )
238      DO jvar = 1,profdata%nvar
239         ALLOCATE( llvvalid(jvar)%luse(profdata%nvprot(jvar)) )
240      END DO
241
242      ! We want all data which has qc flags <= 10
243
244      llvalid%luse(:) = ( profdata%nqc(:)  <= 10 )
245      DO jvar = 1,profdata%nvar
246         llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= 10 )
247      END DO
248
249      ! The actual copying
250
251      CALL obs_prof_compress( profdata,     prodatqc,       .TRUE.,  numout, &
252         &                    lvalid=llvalid, lvvalid=llvvalid )
253
254      ! Dellocate the selection arrays
255      DEALLOCATE( llvalid%luse )
256      DO jvar = 1,profdata%nvar
257         DEALLOCATE( llvvalid(jvar)%luse )
258      END DO
259
260      ! -----------------------------------------------------------------------
261      ! Print information about what observations are left after qc
262      ! -----------------------------------------------------------------------
263
264      ! Update the total observation counter array
265     
266      IF(lwp) THEN
267         WRITE(numout,*)
268         WRITE(numout,*) 'obs_pre_pro :'
269         WRITE(numout,*) '~~~~~~~~~~~'
270         WRITE(numout,*)
271         WRITE(numout,*) ' Profiles outside time domain                = ', &
272            &            iotdobsmpp
273         WRITE(numout,*) ' Remaining profiles that failed grid search  = ', &
274            &            igrdobsmpp
275         WRITE(numout,*) ' Remaining T data outside space domain       = ', &
276            &            iosdtobsmpp
277         WRITE(numout,*) ' Remaining T data at land points             = ', &
278            &            ilantobsmpp
279         IF (ld_nea) THEN
280            WRITE(numout,*) ' Remaining T data near land points (removed) = ',&
281               &            inlatobsmpp
282         ELSE
283            WRITE(numout,*) ' Remaining T data near land points (kept)    = ',&
284               &            inlatobsmpp
285         ENDIF
286         WRITE(numout,*) ' T data accepted                             = ', &
287            &            prodatqc%nvprotmpp(1)
288         WRITE(numout,*) ' Remaining S data outside space domain       = ', &
289            &            iosdsobsmpp
290         WRITE(numout,*) ' Remaining S data at land points             = ', &
291            &            ilansobsmpp
292         IF (ld_nea) THEN
293            WRITE(numout,*) ' Remaining S data near land points (removed) = ',&
294               &            inlasobsmpp
295         ELSE
296            WRITE(numout,*) ' Remaining S data near land points (kept)    = ',&
297               &            inlasobsmpp
298         ENDIF
299         WRITE(numout,*) ' S data accepted                             = ', &
300            &            prodatqc%nvprotmpp(2)
301
302         WRITE(numout,*)
303         WRITE(numout,*) ' Number of observations per time step :'
304         WRITE(numout,*)
305         WRITE(numout,997)
306         WRITE(numout,998)
307      ENDIF
308     
309      DO jobs = 1, prodatqc%nprof
310         inrc = prodatqc%mstp(jobs) + 2 - nit000
311         prodatqc%npstp(inrc)  = prodatqc%npstp(inrc) + 1
312         DO jvar = 1, prodatqc%nvar
313            IF ( prodatqc%npvend(jobs,jvar) > 0 ) THEN
314               prodatqc%nvstp(inrc,jvar) = prodatqc%nvstp(inrc,jvar) + &
315                  &                      ( prodatqc%npvend(jobs,jvar) - &
316                  &                        prodatqc%npvsta(jobs,jvar) + 1 )
317            ENDIF
318         END DO
319      END DO
320     
321     
322      CALL obs_mpp_sum_integers( prodatqc%npstp, prodatqc%npstpmpp, &
323         &                       nitend - nit000 + 2 )
324      DO jvar = 1, prodatqc%nvar
325         CALL obs_mpp_sum_integers( prodatqc%nvstp(:,jvar), &
326            &                       prodatqc%nvstpmpp(:,jvar), &
327            &                       nitend - nit000 + 2 )
328      END DO
329
330      IF ( lwp ) THEN
331         DO jstp = nit000 - 1, nitend
332            inrc = jstp - nit000 + 2
333            WRITE(numout,999) jstp, prodatqc%npstpmpp(inrc), &
334               &                    prodatqc%nvstpmpp(inrc,1), &
335               &                    prodatqc%nvstpmpp(inrc,2)
336         END DO
337      ENDIF
338
339997   FORMAT(10X,'Time step',5X,'Profiles',5X,'Temperature',5X,'Salinity')
340998   FORMAT(10X,'---------',5X,'--------',5X,'-----------',5X,'--------')
341999   FORMAT(10X,I9,5X,I8,5X,I11,5X,I8)
342     
343   END SUBROUTINE obs_pre_pro
344
345   SUBROUTINE obs_pre_sla( sladata, sladatqc, ld_sla, ld_nea )
346      !!----------------------------------------------------------------------
347      !!                    ***  ROUTINE obs_pre_sla  ***
348      !!
349      !! ** Purpose : First level check and screening of SLA observations
350      !!
351      !! ** Method  : First level check and screening of SLA observations
352      !!
353      !! ** Action  :
354      !!
355      !! References :
356      !!   
357      !! History :
358      !!        !  2007-03  (A. Weaver, K. Mogensen) Original
359      !!        !  2007-06  (K. Mogensen et al) Reject obs. near land.
360      !!----------------------------------------------------------------------
361      !! * Modules used
362      USE domstp              ! Domain: set the time-step
363      USE par_oce             ! Ocean parameters
364      USE dom_oce, ONLY : &   ! Geographical information
365         & glamt,   &
366         & gphit,   &
367         & tmask         
368      !! * Arguments
369      TYPE(obs_surf), INTENT(INOUT) :: sladata    ! Full set of SLA data
370      TYPE(obs_surf), INTENT(INOUT) :: sladatqc   ! Subset of SLA data not failing screening
371      LOGICAL, INTENT(IN) :: ld_sla         ! Switch for SLA data
372      LOGICAL, INTENT(IN) :: ld_nea         ! Switch for rejecting observation near land
373      !! * Local declarations
374      INTEGER :: iyea0        ! Initial date
375      INTEGER :: imon0        !  - (year, month, day, hour, minute)
376      INTEGER :: iday0   
377      INTEGER :: ihou0   
378      INTEGER :: imin0
379      INTEGER :: icycle       ! Current assimilation cycle
380                              ! Counters for observations that
381      INTEGER :: iotdobs      !  - outside time domain
382      INTEGER :: iosdsobs     !  - outside space domain
383      INTEGER :: ilansobs     !  - within a model land cell
384      INTEGER :: inlasobs     !  - close to land
385      INTEGER :: igrdobs      !  - fail the grid search
386                              ! Global counters for observations that
387      INTEGER :: iotdobsmpp     !  - outside time domain
388      INTEGER :: iosdsobsmpp    !  - outside space domain
389      INTEGER :: ilansobsmpp    !  - within a model land cell
390      INTEGER :: inlasobsmpp    !  - close to land
391      INTEGER :: igrdobsmpp     !  - fail the grid search
392      LOGICAL, DIMENSION(:), ALLOCATABLE :: & 
393         & llvalid            ! SLA data selection
394      INTEGER :: jobs         ! Obs. loop variable
395      INTEGER :: jstp         ! Time loop variable
396      INTEGER :: inrc         ! Time index variable
397
398      IF(lwp) WRITE(numout,*)'obs_pre_sla : Preparing the SLA observations...'
399
400      ! Initial date initialization (year, month, day, hour, minute)
401      iyea0 =   ndate0 / 10000
402      imon0 = ( ndate0 - iyea0 * 10000 ) / 100
403      iday0 =   ndate0 - iyea0 * 10000 - imon0 * 100
404      ihou0 = 0
405      imin0 = 0
406
407      icycle = no     ! Assimilation cycle
408
409      ! Diagnotics counters for various failures.
410
411      iotdobs  = 0
412      igrdobs  = 0
413      iosdsobs = 0
414      ilansobs = 0
415      inlasobs = 0
416
417      ! -----------------------------------------------------------------------
418      ! Find time coordinate for SLA data
419      ! -----------------------------------------------------------------------
420
421      CALL obs_coo_tim( icycle, &
422         &              iyea0,   imon0,   iday0,   ihou0,   imin0,      &
423         &              sladata%nsurf,   sladata%nyea, sladata%nmon, &
424         &              sladata%nday,    sladata%nhou, sladata%nmin, &
425         &              sladata%nqc,     sladata%mstp, iotdobs        )
426
427      CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
428     
429      ! -----------------------------------------------------------------------
430      ! Check for SLA data failing the grid search
431      ! -----------------------------------------------------------------------
432
433      CALL obs_coo_grd( sladata%nsurf,   sladata%mi, sladata%mj, &
434         &              sladata%nqc,     igrdobs                         )
435
436      CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
437
438      ! -----------------------------------------------------------------------
439      ! Check for land points.
440      ! -----------------------------------------------------------------------
441
442      CALL obs_coo_spc_2d( sladata%nsurf,              &
443         &                 jpi,          jpj,          &
444         &                 sladata%mi,   sladata%mj,   & 
445         &                 sladata%rlam, sladata%rphi, &
446         &                 glamt,        gphit,        &
447         &                 tmask(:,:,1), sladata%nqc,  &
448         &                 iosdsobs,     ilansobs,     &
449         &                 inlasobs,     ld_nea        )
450
451      CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
452      CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
453      CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
454
455      ! -----------------------------------------------------------------------
456      ! Copy useful data from the sladata data structure to
457      ! the sladatqc data structure
458      ! -----------------------------------------------------------------------
459
460      ! Allocate the selection arrays
461
462      ALLOCATE( llvalid(sladata%nsurf) )
463     
464      ! We want all data which has qc flags <= 10
465
466      llvalid(:)  = ( sladata%nqc(:)  <= 10 )
467
468      ! The actual copying
469
470      CALL obs_surf_compress( sladata,     sladatqc,       .TRUE.,  numout, &
471         &                    lvalid=llvalid )
472
473      ! Dellocate the selection arrays
474      DEALLOCATE( llvalid )
475
476      ! -----------------------------------------------------------------------
477      ! Print information about what observations are left after qc
478      ! -----------------------------------------------------------------------
479
480      ! Update the total observation counter array
481     
482      IF(lwp) THEN
483         WRITE(numout,*)
484         WRITE(numout,*) 'obs_pre_sla :'
485         WRITE(numout,*) '~~~~~~~~~~~'
486         WRITE(numout,*)
487         WRITE(numout,*) ' SLA data outside time domain                  = ', &
488            &            iotdobsmpp
489         WRITE(numout,*) ' Remaining SLA data that failed grid search    = ', &
490            &            igrdobsmpp
491         WRITE(numout,*) ' Remaining SLA data outside space domain       = ', &
492            &            iosdsobsmpp
493         WRITE(numout,*) ' Remaining SLA data at land points             = ', &
494            &            ilansobsmpp
495         IF (ld_nea) THEN
496            WRITE(numout,*) ' Remaining SLA data near land points (removed) = ', &
497               &            inlasobsmpp
498         ELSE
499            WRITE(numout,*) ' Remaining SLA data near land points (kept)    = ', &
500               &            inlasobsmpp
501         ENDIF
502         WRITE(numout,*) ' SLA data accepted                             = ', &
503            &            sladatqc%nsurfmpp
504
505         WRITE(numout,*)
506         WRITE(numout,*) ' Number of observations per time step :'
507         WRITE(numout,*)
508         WRITE(numout,1997)
509         WRITE(numout,1998)
510      ENDIF
511     
512      DO jobs = 1, sladatqc%nsurf
513         inrc = sladatqc%mstp(jobs) + 2 - nit000
514         sladatqc%nsstp(inrc)  = sladatqc%nsstp(inrc) + 1
515      END DO
516     
517      CALL obs_mpp_sum_integers( sladatqc%nsstp, sladatqc%nsstpmpp, &
518         &                       nitend - nit000 + 2 )
519
520      IF ( lwp ) THEN
521         DO jstp = nit000 - 1, nitend
522            inrc = jstp - nit000 + 2
523            WRITE(numout,1999) jstp, sladatqc%nsstpmpp(inrc)
524         END DO
525      ENDIF
526
5271997  FORMAT(10X,'Time step',5X,'Sea level anomaly')
5281998  FORMAT(10X,'---------',5X,'-----------------')
5291999  FORMAT(10X,I9,5X,I17)
530
531   END SUBROUTINE obs_pre_sla
532
533   SUBROUTINE obs_pre_sst( sstdata, sstdatqc, ld_sst, ld_nea )
534      !!----------------------------------------------------------------------
535      !!                    ***  ROUTINE obs_pre_sst  ***
536      !!
537      !! ** Purpose : First level check and screening of SST observations
538      !!
539      !! ** Method  : First level check and screening of SST observations
540      !!
541      !! ** Action  :
542      !!
543      !! References :
544      !!   
545      !! History :
546      !!        !  2007-03  (S. Ricci) SST data preparation
547      !!----------------------------------------------------------------------
548      !! * Modules used
549      USE domstp              ! Domain: set the time-step
550      USE par_oce             ! Ocean parameters
551      USE dom_oce, ONLY : &   ! Geographical information
552         & glamt,   &
553         & gphit,   &
554         & tmask
555      !! * Arguments
556      TYPE(obs_surf), INTENT(INOUT) :: sstdata     ! Full set of SST data
557      TYPE(obs_surf), INTENT(INOUT) :: sstdatqc    ! Subset of SST data not failing screening
558      LOGICAL :: ld_sst             ! Switch for SST data
559      LOGICAL :: ld_nea             ! Switch for rejecting observation near land
560      !! * Local declarations
561      INTEGER :: iyea0        ! Initial date
562      INTEGER :: imon0        !  - (year, month, day, hour, minute)
563      INTEGER :: iday0   
564      INTEGER :: ihou0   
565      INTEGER :: imin0
566      INTEGER :: icycle       ! Current assimilation cycle
567                              ! Counters for observations that
568      INTEGER :: iotdobs      !  - outside time domain
569      INTEGER :: iosdsobs     !  - outside space domain
570      INTEGER :: ilansobs     !  - within a model land cell
571      INTEGER :: inlasobs     !  - close to land
572      INTEGER :: igrdobs      !  - fail the grid search
573                              ! Global counters for observations that
574      INTEGER :: iotdobsmpp   !  - outside time domain
575      INTEGER :: iosdsobsmpp  !  - outside space domain
576      INTEGER :: ilansobsmpp  !  - within a model land cell
577      INTEGER :: inlasobsmpp  !  - close to land
578      INTEGER :: igrdobsmpp   !  - fail the grid search
579      LOGICAL, DIMENSION(:), ALLOCATABLE :: & 
580         & llvalid            ! SST data selection
581      INTEGER :: jobs         ! Obs. loop variable
582      INTEGER :: jstp         ! Time loop variable
583      INTEGER :: inrc         ! Time index variable
584
585      IF(lwp) WRITE(numout,*)'obs_pre_sst : Preparing the SST observations...'
586
587      ! Initial date initialization (year, month, day, hour, minute)
588      iyea0 =   ndate0 / 10000
589      imon0 = ( ndate0 - iyea0 * 10000 ) / 100
590      iday0 =   ndate0 - iyea0 * 10000 - imon0 * 100
591      ihou0 = 0
592      imin0 = 0
593
594      icycle = no     ! Assimilation cycle
595
596      ! Diagnotics counters for various failures.
597
598      iotdobs  = 0
599      igrdobs  = 0
600      iosdsobs = 0
601      ilansobs = 0
602      inlasobs = 0
603
604      ! -----------------------------------------------------------------------
605      ! Find time coordinate for SST data
606      ! -----------------------------------------------------------------------
607
608      CALL obs_coo_tim( icycle, &
609         &              iyea0,   imon0,   iday0,   ihou0,   imin0,      &
610         &              sstdata%nsurf,   sstdata%nyea, sstdata%nmon, &
611         &              sstdata%nday,    sstdata%nhou, sstdata%nmin, &
612         &              sstdata%nqc,     sstdata%mstp, iotdobs        )
613      CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
614      ! -----------------------------------------------------------------------
615      ! Check for SST data failing the grid search
616      ! -----------------------------------------------------------------------
617
618      CALL obs_coo_grd( sstdata%nsurf,   sstdata%mi, sstdata%mj, &
619         &              sstdata%nqc,     igrdobs                         )
620      CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
621
622      ! -----------------------------------------------------------------------
623      ! Check for land points.
624      ! -----------------------------------------------------------------------
625
626      CALL obs_coo_spc_2d( sstdata%nsurf,              &
627         &                 jpi,          jpj,          &
628         &                 sstdata%mi,   sstdata%mj,   & 
629         &                 sstdata%rlam, sstdata%rphi, &
630         &                 glamt,        gphit,        &
631         &                 tmask(:,:,1), sstdata%nqc,  &
632         &                 iosdsobs,     ilansobs,     &
633         &                 inlasobs,     ld_nea        )
634
635      CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
636      CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
637      CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
638
639      ! -----------------------------------------------------------------------
640      ! Copy useful data from the sstdata data structure to
641      ! the sstdatqc data structure
642      ! -----------------------------------------------------------------------
643
644      ! Allocate the selection arrays
645
646      ALLOCATE( llvalid(sstdata%nsurf) )
647     
648      ! We want all data which has qc flags <= 0
649
650      llvalid(:)  = ( sstdata%nqc(:)  <= 10 )
651
652      ! The actual copying
653
654      CALL obs_surf_compress( sstdata,     sstdatqc,       .TRUE.,  numout, &
655         &                    lvalid=llvalid )
656
657      ! Dellocate the selection arrays
658      DEALLOCATE( llvalid )
659
660      ! -----------------------------------------------------------------------
661      ! Print information about what observations are left after qc
662      ! -----------------------------------------------------------------------
663
664      ! Update the total observation counter array
665     
666      IF(lwp) THEN
667         WRITE(numout,*)
668         WRITE(numout,*) 'obs_pre_sst :'
669         WRITE(numout,*) '~~~~~~~~~~~'
670         WRITE(numout,*)
671         WRITE(numout,*) ' SST data outside time domain                  = ', &
672            &            iotdobsmpp
673         WRITE(numout,*) ' Remaining SST data that failed grid search    = ', &
674            &            igrdobsmpp
675         WRITE(numout,*) ' Remaining SST data outside space domain       = ', &
676            &            iosdsobsmpp
677         WRITE(numout,*) ' Remaining SST data at land points             = ', &
678            &            ilansobsmpp
679         IF (ld_nea) THEN
680            WRITE(numout,*) ' Remaining SST data near land points (removed) = ', &
681               &            inlasobsmpp
682         ELSE
683            WRITE(numout,*) ' Remaining SST data near land points (kept)    = ', &
684               &            inlasobsmpp
685         ENDIF
686         WRITE(numout,*) ' SST data accepted                             = ', &
687            &            sstdatqc%nsurfmpp
688
689         WRITE(numout,*)
690         WRITE(numout,*) ' Number of observations per time step :'
691         WRITE(numout,*)
692         WRITE(numout,1997)
693         WRITE(numout,1998)
694      ENDIF
695     
696      DO jobs = 1, sstdatqc%nsurf
697         inrc = sstdatqc%mstp(jobs) + 2 - nit000
698         sstdatqc%nsstp(inrc)  = sstdatqc%nsstp(inrc) + 1
699      END DO
700     
701      CALL obs_mpp_sum_integers( sstdatqc%nsstp, sstdatqc%nsstpmpp, &
702         &                       nitend - nit000 + 2 )
703
704      IF ( lwp ) THEN
705         DO jstp = nit000 - 1, nitend
706            inrc = jstp - nit000 + 2
707            WRITE(numout,1999) jstp, sstdatqc%nsstpmpp(inrc)
708         END DO
709      ENDIF
710
7111997  FORMAT(10X,'Time step',5X,'Sea surface temperature')
7121998  FORMAT(10X,'---------',5X,'-----------------')
7131999  FORMAT(10X,I9,5X,I17)
714     
715   END SUBROUTINE obs_pre_sst
716
717   SUBROUTINE obs_pre_seaice( seaicedata, seaicedatqc, ld_seaice, ld_nea )
718      !!----------------------------------------------------------------------
719      !!                    ***  ROUTINE obs_pre_seaice  ***
720      !!
721      !! ** Purpose : First level check and screening of Sea Ice observations
722      !!
723      !! ** Method  : First level check and screening of Sea Ice observations
724      !!
725      !! ** Action  :
726      !!
727      !! References :
728      !!   
729      !! History :
730      !!        !  2007-11 (D. Lea) based on obs_pre_sst
731      !!----------------------------------------------------------------------
732      !! * Modules used
733      USE domstp              ! Domain: set the time-step
734      USE par_oce             ! Ocean parameters
735      USE dom_oce, ONLY : &   ! Geographical information
736         & glamt,   &
737         & gphit,   &
738         & tmask
739      !! * Arguments
740      TYPE(obs_surf), INTENT(INOUT) :: seaicedata     ! Full set of Sea Ice data
741      TYPE(obs_surf), INTENT(INOUT) :: seaicedatqc    ! Subset of sea ice data not failing screening
742      LOGICAL :: ld_seaice     ! Switch for sea ice data
743      LOGICAL :: ld_nea        ! Switch for rejecting observation near land
744      !! * Local declarations
745      INTEGER :: iyea0         ! Initial date
746      INTEGER :: imon0         !  - (year, month, day, hour, minute)
747      INTEGER :: iday0   
748      INTEGER :: ihou0   
749      INTEGER :: imin0
750      INTEGER :: icycle       ! Current assimilation cycle
751                              ! Counters for observations that
752      INTEGER :: iotdobs      !  - outside time domain
753      INTEGER :: iosdsobs     !  - outside space domain
754      INTEGER :: ilansobs     !  - within a model land cell
755      INTEGER :: inlasobs     !  - close to land
756      INTEGER :: igrdobs      !  - fail the grid search
757                              ! Global counters for observations that
758      INTEGER :: iotdobsmpp   !  - outside time domain
759      INTEGER :: iosdsobsmpp  !  - outside space domain
760      INTEGER :: ilansobsmpp  !  - within a model land cell
761      INTEGER :: inlasobsmpp  !  - close to land
762      INTEGER :: igrdobsmpp   !  - fail the grid search
763      LOGICAL, DIMENSION(:), ALLOCATABLE :: & 
764         & llvalid            ! data selection
765      INTEGER :: jobs         ! Obs. loop variable
766      INTEGER :: jstp         ! Time loop variable
767      INTEGER :: inrc         ! Time index variable
768
769      IF (lwp) WRITE(numout,*)'obs_pre_seaice : Preparing the sea ice observations...'
770
771      ! Initial date initialization (year, month, day, hour, minute)
772      iyea0 =   ndate0 / 10000
773      imon0 = ( ndate0 - iyea0 * 10000 ) / 100
774      iday0 =   ndate0 - iyea0 * 10000 - imon0 * 100
775      ihou0 = 0
776      imin0 = 0
777
778      icycle = no     ! Assimilation cycle
779
780      ! Diagnotics counters for various failures.
781
782      iotdobs  = 0
783      igrdobs  = 0
784      iosdsobs = 0
785      ilansobs = 0
786      inlasobs = 0
787
788      ! -----------------------------------------------------------------------
789      ! Find time coordinate for sea ice data
790      ! -----------------------------------------------------------------------
791
792      CALL obs_coo_tim( icycle, &
793         &              iyea0,   imon0,   iday0,   ihou0,   imin0,      &
794         &              seaicedata%nsurf,   seaicedata%nyea, seaicedata%nmon, &
795         &              seaicedata%nday,    seaicedata%nhou, seaicedata%nmin, &
796         &              seaicedata%nqc,     seaicedata%mstp, iotdobs        )
797      CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
798      ! -----------------------------------------------------------------------
799      ! Check for sea ice data failing the grid search
800      ! -----------------------------------------------------------------------
801
802      CALL obs_coo_grd( seaicedata%nsurf,   seaicedata%mi, seaicedata%mj, &
803         &              seaicedata%nqc,     igrdobs                         )
804      CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
805
806      ! -----------------------------------------------------------------------
807      ! Check for land points.
808      ! -----------------------------------------------------------------------
809
810      CALL obs_coo_spc_2d( seaicedata%nsurf,                 &
811         &                 jpi,             jpj,             &
812         &                 seaicedata%mi,   seaicedata%mj,   & 
813         &                 seaicedata%rlam, seaicedata%rphi, &
814         &                 glamt,           gphit,           &
815         &                 tmask(:,:,1),    seaicedata%nqc,  &
816         &                 iosdsobs,        ilansobs,        &
817         &                 inlasobs,        ld_nea           )
818
819      CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
820      CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
821      CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
822
823      ! -----------------------------------------------------------------------
824      ! Copy useful data from the seaicedata data structure to
825      ! the seaicedatqc data structure
826      ! -----------------------------------------------------------------------
827
828      ! Allocate the selection arrays
829
830      ALLOCATE( llvalid(seaicedata%nsurf) )
831     
832      ! We want all data which has qc flags <= 0
833
834      llvalid(:)  = ( seaicedata%nqc(:)  <= 10 )
835
836      ! The actual copying
837
838      CALL obs_surf_compress( seaicedata,     seaicedatqc,       .TRUE.,  numout, &
839         &                    lvalid=llvalid )
840
841      ! Dellocate the selection arrays
842      DEALLOCATE( llvalid )
843
844      ! -----------------------------------------------------------------------
845      ! Print information about what observations are left after qc
846      ! -----------------------------------------------------------------------
847
848      ! Update the total observation counter array
849     
850      IF(lwp) THEN
851         WRITE(numout,*)
852         WRITE(numout,*) 'obs_pre_seaice :'
853         WRITE(numout,*) '~~~~~~~~~~~'
854         WRITE(numout,*)
855         WRITE(numout,*) ' Sea ice data outside time domain                  = ', &
856            &            iotdobsmpp
857         WRITE(numout,*) ' Remaining sea ice data that failed grid search    = ', &
858            &            igrdobsmpp
859         WRITE(numout,*) ' Remaining sea ice data outside space domain       = ', &
860            &            iosdsobsmpp
861         WRITE(numout,*) ' Remaining sea ice data at land points             = ', &
862            &            ilansobsmpp
863         IF (ld_nea) THEN
864            WRITE(numout,*) ' Remaining sea ice data near land points (removed) = ', &
865               &            inlasobsmpp
866         ELSE
867            WRITE(numout,*) ' Remaining sea ice data near land points (kept)    = ', &
868               &            inlasobsmpp
869         ENDIF
870         WRITE(numout,*) ' Sea ice data accepted                             = ', &
871            &            seaicedatqc%nsurfmpp
872
873         WRITE(numout,*)
874         WRITE(numout,*) ' Number of observations per time step :'
875         WRITE(numout,*)
876         WRITE(numout,1997)
877         WRITE(numout,1998)
878      ENDIF
879     
880      DO jobs = 1, seaicedatqc%nsurf
881         inrc = seaicedatqc%mstp(jobs) + 2 - nit000
882         seaicedatqc%nsstp(inrc)  = seaicedatqc%nsstp(inrc) + 1
883      END DO
884     
885      CALL obs_mpp_sum_integers( seaicedatqc%nsstp, seaicedatqc%nsstpmpp, &
886         &                       nitend - nit000 + 2 )
887
888      IF ( lwp ) THEN
889         DO jstp = nit000 - 1, nitend
890            inrc = jstp - nit000 + 2
891            WRITE(numout,1999) jstp, seaicedatqc%nsstpmpp(inrc)
892         END DO
893      ENDIF
894
8951997  FORMAT(10X,'Time step',5X,'Sea ice data           ')
8961998  FORMAT(10X,'---------',5X,'-----------------')
8971999  FORMAT(10X,I9,5X,I17)
898     
899   END SUBROUTINE obs_pre_seaice
900
901   SUBROUTINE obs_pre_vel( profdata, prodatqc, ld_vel3d, ld_nea, ld_dailyav )
902      !!----------------------------------------------------------------------
903      !!                    ***  ROUTINE obs_pre_taovel  ***
904      !!
905      !! ** Purpose : First level check and screening of U and V profiles
906      !!
907      !! ** Method  : First level check and screening of U and V profiles
908      !!
909      !! History :
910      !!        !  2007-06  (K. Mogensen) original : T and S profile data
911      !!        !  2008-09  (M. Valdivieso) : TAO velocity data
912      !!        !  2009-01  (K. Mogensen) : New feedback strictuer
913      !!
914      !!----------------------------------------------------------------------
915      !! * Modules used
916      USE domstp              ! Domain: set the time-step
917      USE par_oce             ! Ocean parameters
918      USE dom_oce, ONLY : &   ! Geographical information
919         & glamt, glamu, glamv,    &
920         & gphit, gphiu, gphiv,    &
921         & gdept_1d,             &
922         & tmask, umask, vmask
923      !! * Arguments
924      TYPE(obs_prof), INTENT(INOUT) :: profdata   ! Full set of profile data
925      TYPE(obs_prof), INTENT(INOUT) :: prodatqc   ! Subset of profile data not failing screening
926      LOGICAL, INTENT(IN) :: ld_vel3d      ! Switch for zonal and meridional velocity components
927      LOGICAL, INTENT(IN) :: ld_nea        ! Switch for rejecting observation near land
928      LOGICAL, INTENT(IN) :: ld_dailyav    ! Switch for daily average data
929      !! * Local declarations
930      INTEGER :: iyea0        ! Initial date
931      INTEGER :: imon0        !  - (year, month, day, hour, minute)
932      INTEGER :: iday0   
933      INTEGER :: ihou0   
934      INTEGER :: imin0
935      INTEGER :: icycle       ! Current assimilation cycle
936                              ! Counters for observations that
937      INTEGER :: iotdobs      !  - outside time domain
938      INTEGER :: iosduobs     !  - outside space domain (zonal velocity component)
939      INTEGER :: iosdvobs     !  - outside space domain (meridional velocity component)
940      INTEGER :: ilanuobs     !  - within a model land cell (zonal velocity component)
941      INTEGER :: ilanvobs     !  - within a model land cell (meridional velocity component)
942      INTEGER :: inlauobs     !  - close to land (zonal velocity component)
943      INTEGER :: inlavobs     !  - close to land (meridional velocity component)
944      INTEGER :: igrdobs      !  - fail the grid search
945      INTEGER :: iuvchku      !  - reject u if v rejected and vice versa
946      INTEGER :: iuvchkv      !
947                              ! Global counters for observations that
948      INTEGER :: iotdobsmpp   !  - outside time domain
949      INTEGER :: iosduobsmpp  !  - outside space domain (zonal velocity component)
950      INTEGER :: iosdvobsmpp  !  - outside space domain (meridional velocity component)
951      INTEGER :: ilanuobsmpp  !  - within a model land cell (zonal velocity component)
952      INTEGER :: ilanvobsmpp  !  - within a model land cell (meridional velocity component)
953      INTEGER :: inlauobsmpp  !  - close to land (zonal velocity component)
954      INTEGER :: inlavobsmpp  !  - close to land (meridional velocity component)
955      INTEGER :: igrdobsmpp   !  - fail the grid search
956      INTEGER :: iuvchkumpp   !  - reject u if v rejected and vice versa
957      INTEGER :: iuvchkvmpp   !
958      TYPE(obs_prof_valid) ::  llvalid      ! Profile selection
959      TYPE(obs_prof_valid), DIMENSION(profdata%nvar) :: &
960         & llvvalid           ! U,V selection
961      INTEGER :: jvar         ! Variable loop variable
962      INTEGER :: jobs         ! Obs. loop variable
963      INTEGER :: jstp         ! Time loop variable
964      INTEGER :: inrc         ! Time index variable
965
966      IF(lwp) WRITE(numout,*)'obs_pre_vel: Preparing the velocity profile data'
967
968      ! Initial date initialization (year, month, day, hour, minute)
969      iyea0 =   ndate0 / 10000
970      imon0 = ( ndate0 - iyea0 * 10000 ) / 100
971      iday0 =   ndate0 - iyea0 * 10000 - imon0 * 100
972      ihou0 = 0
973      imin0 = 0
974
975      icycle = no     ! Assimilation cycle
976
977      ! Diagnotics counters for various failures.
978
979      iotdobs  = 0
980      igrdobs  = 0
981      iosduobs = 0
982      iosdvobs = 0
983      ilanuobs = 0
984      ilanvobs = 0
985      inlauobs = 0
986      inlavobs = 0
987      iuvchku  = 0
988      iuvchkv = 0
989
990      ! -----------------------------------------------------------------------
991      ! Find time coordinate for profiles
992      ! -----------------------------------------------------------------------
993
994      CALL obs_coo_tim_prof( icycle, &
995         &              iyea0,   imon0,   iday0,   ihou0,   imin0,      &
996         &              profdata%nprof,   profdata%nyea, profdata%nmon, &
997         &              profdata%nday,    profdata%nhou, profdata%nmin, &
998         &              profdata%ntyp,    profdata%nqc,  profdata%mstp, &
999         &              iotdobs, ld_dailyav = ld_dailyav        )
1000   
1001      CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
1002     
1003      ! -----------------------------------------------------------------------
1004      ! Check for profiles failing the grid search
1005      ! -----------------------------------------------------------------------
1006
1007      CALL obs_coo_grd( profdata%nprof,   profdata%mi(:,1), profdata%mj(:,1), &
1008         &              profdata%nqc,     igrdobs                         )
1009      CALL obs_coo_grd( profdata%nprof,   profdata%mi(:,2), profdata%mj(:,2), &
1010         &              profdata%nqc,     igrdobs                         )
1011
1012      CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
1013
1014      ! -----------------------------------------------------------------------
1015      ! Reject all observations for profiles with nqc > 10
1016      ! -----------------------------------------------------------------------
1017
1018      CALL obs_pro_rej( profdata )
1019
1020      ! -----------------------------------------------------------------------
1021      ! Check for land points. This includes points below the model
1022      ! bathymetry so this is done for every point in the profile
1023      ! -----------------------------------------------------------------------
1024
1025      ! Zonal Velocity Component
1026
1027      CALL obs_coo_spc_3d( profdata%nprof,        profdata%nvprot(1),   &
1028         &                 profdata%npvsta(:,1),  profdata%npvend(:,1), &
1029         &                 jpi,                   jpj,                  &
1030         &                 jpk,                                         &
1031         &                 profdata%mi,           profdata%mj,          & 
1032         &                 profdata%var(1)%mvk,                         &
1033         &                 profdata%rlam,         profdata%rphi,        &
1034         &                 profdata%var(1)%vdep,                        &
1035         &                 glamu,                 gphiu,                &
1036         &                 gdept_1d,              umask,                &
1037         &                 profdata%nqc,          profdata%var(1)%nvqc, &
1038         &                 iosduobs,              ilanuobs,             &
1039         &                 inlauobs,              ld_nea                )
1040
1041      CALL obs_mpp_sum_integer( iosduobs, iosduobsmpp )
1042      CALL obs_mpp_sum_integer( ilanuobs, ilanuobsmpp )
1043      CALL obs_mpp_sum_integer( inlauobs, inlauobsmpp )
1044
1045      ! Meridional Velocity Component
1046
1047      CALL obs_coo_spc_3d( profdata%nprof,        profdata%nvprot(2),   &
1048         &                 profdata%npvsta(:,2),  profdata%npvend(:,2), &
1049         &                 jpi,                   jpj,                  &
1050         &                 jpk,                                         &
1051         &                 profdata%mi,           profdata%mj,          & 
1052         &                 profdata%var(2)%mvk,                         &
1053         &                 profdata%rlam,         profdata%rphi,        &
1054         &                 profdata%var(2)%vdep,                        &
1055         &                 glamv,                 gphiv,                &
1056         &                 gdept_1d,              vmask,                &
1057         &                 profdata%nqc,          profdata%var(2)%nvqc, &
1058         &                 iosdvobs,              ilanvobs,             &
1059         &                 inlavobs,              ld_nea                )
1060
1061      CALL obs_mpp_sum_integer( iosdvobs, iosdvobsmpp )
1062      CALL obs_mpp_sum_integer( ilanvobs, ilanvobsmpp )
1063      CALL obs_mpp_sum_integer( inlavobs, inlavobsmpp )
1064
1065      ! -----------------------------------------------------------------------
1066      ! Reject u if v is rejected and vice versa
1067      ! -----------------------------------------------------------------------
1068
1069      CALL obs_uv_rej( profdata, iuvchku, iuvchkv )
1070      CALL obs_mpp_sum_integer( iuvchku, iuvchkumpp )
1071      CALL obs_mpp_sum_integer( iuvchkv, iuvchkvmpp )
1072
1073      ! -----------------------------------------------------------------------
1074      ! Copy useful data from the profdata data structure to
1075      ! the prodatqc data structure
1076      ! -----------------------------------------------------------------------
1077
1078      ! Allocate the selection arrays
1079
1080      ALLOCATE( llvalid%luse(profdata%nprof) )
1081      DO jvar = 1,profdata%nvar
1082         ALLOCATE( llvvalid(jvar)%luse(profdata%nvprot(jvar)) )
1083      END DO
1084
1085      ! We want all data which has qc flags = 0
1086
1087      llvalid%luse(:) = ( profdata%nqc(:)  <= 10 )
1088      DO jvar = 1,profdata%nvar
1089         llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= 10 )
1090      END DO
1091
1092      ! The actual copying
1093
1094      CALL obs_prof_compress( profdata,     prodatqc,       .TRUE.,  numout, &
1095         &                    lvalid=llvalid, lvvalid=llvvalid )
1096
1097      ! Dellocate the selection arrays
1098      DEALLOCATE( llvalid%luse )
1099      DO jvar = 1,profdata%nvar
1100         DEALLOCATE( llvvalid(jvar)%luse )
1101      END DO
1102
1103      ! -----------------------------------------------------------------------
1104      ! Print information about what observations are left after qc
1105      ! -----------------------------------------------------------------------
1106
1107      ! Update the total observation counter array
1108     
1109      IF(lwp) THEN
1110         WRITE(numout,*)
1111         WRITE(numout,*) 'obs_pre_vel :'
1112         WRITE(numout,*) '~~~~~~~~~~~'
1113         WRITE(numout,*)
1114         WRITE(numout,*) ' Profiles outside time domain                = ', &
1115            &            iotdobsmpp
1116         WRITE(numout,*) ' Remaining profiles that failed grid search  = ', &
1117            &            igrdobsmpp
1118         WRITE(numout,*) ' Remaining U data outside space domain       = ', &
1119            &            iosduobsmpp
1120         WRITE(numout,*) ' Remaining U data at land points             = ', &
1121            &            ilanuobsmpp
1122         IF (ld_nea) THEN
1123            WRITE(numout,*) ' Remaining U data near land points (removed) = ',&
1124               &            inlauobsmpp
1125         ELSE
1126            WRITE(numout,*) ' Remaining U data near land points (kept)    = ',&
1127               &            inlauobsmpp
1128         ENDIF
1129         WRITE(numout,*) ' U observation rejected since V rejected     = ', &
1130            &            iuvchku     
1131         WRITE(numout,*) ' U data accepted                             = ', &
1132            &            prodatqc%nvprotmpp(1)
1133         WRITE(numout,*) ' Remaining V data outside space domain       = ', &
1134            &            iosdvobsmpp
1135         WRITE(numout,*) ' Remaining V data at land points             = ', &
1136            &            ilanvobsmpp
1137         IF (ld_nea) THEN
1138            WRITE(numout,*) ' Remaining V data near land points (removed) = ',&
1139               &            inlavobsmpp
1140         ELSE
1141            WRITE(numout,*) ' Remaining V data near land points (kept)    = ',&
1142               &            inlavobsmpp
1143         ENDIF
1144         WRITE(numout,*) ' V observation rejected since U rejected     = ', &
1145            &            iuvchkv     
1146         WRITE(numout,*) ' V data accepted                             = ', &
1147            &            prodatqc%nvprotmpp(2)
1148
1149         WRITE(numout,*)
1150         WRITE(numout,*) ' Number of observations per time step :'
1151         WRITE(numout,*)
1152         WRITE(numout,997)
1153         WRITE(numout,998)
1154      ENDIF
1155     
1156      DO jobs = 1, prodatqc%nprof
1157         inrc = prodatqc%mstp(jobs) + 2 - nit000
1158         prodatqc%npstp(inrc)  = prodatqc%npstp(inrc) + 1
1159         DO jvar = 1, prodatqc%nvar
1160            IF ( prodatqc%npvend(jobs,jvar) > 0 ) THEN
1161               prodatqc%nvstp(inrc,jvar) = prodatqc%nvstp(inrc,jvar) + &
1162                  &                      ( prodatqc%npvend(jobs,jvar) - &
1163                  &                        prodatqc%npvsta(jobs,jvar) + 1 )
1164            ENDIF
1165         END DO
1166      END DO
1167     
1168     
1169      CALL obs_mpp_sum_integers( prodatqc%npstp, prodatqc%npstpmpp, &
1170         &                       nitend - nit000 + 2 )
1171      DO jvar = 1, prodatqc%nvar
1172         CALL obs_mpp_sum_integers( prodatqc%nvstp(:,jvar), &
1173            &                       prodatqc%nvstpmpp(:,jvar), &
1174            &                       nitend - nit000 + 2 )
1175      END DO
1176
1177      IF ( lwp ) THEN
1178         DO jstp = nit000 - 1, nitend
1179            inrc = jstp - nit000 + 2
1180            WRITE(numout,999) jstp, prodatqc%npstpmpp(inrc), &
1181               &                    prodatqc%nvstpmpp(inrc,1), &
1182               &                    prodatqc%nvstpmpp(inrc,2)
1183         END DO
1184      ENDIF
1185
1186997   FORMAT(10X,'Time step',5X,'Profiles',5X,'Zonal Comp.',5X,'Meridional Comp.')
1187998   FORMAT(10X,'---------',5X,'--------',5X,'-----------',5X,'----------------')
1188999   FORMAT(10X,I9,5X,I8,5X,I11,5X,I8)
1189
1190   END SUBROUTINE obs_pre_vel
1191
1192   SUBROUTINE obs_pre_logchl( logchldata, logchldatqc, ld_logchl, ld_nea )
1193      !!----------------------------------------------------------------------
1194      !!                    ***  ROUTINE obs_pre_logchl  ***
1195      !!
1196      !! ** Purpose : First level check and screening of logchl observations
1197      !!
1198      !! ** Method  : First level check and screening of logchl observations
1199      !!
1200      !! ** Action  :
1201      !!
1202      !! References :
1203      !!   
1204      !! History :
1205      !!----------------------------------------------------------------------
1206      !! * Modules used
1207      USE domstp              ! Domain: set the time-step
1208      USE par_oce             ! Ocean parameters
1209      USE dom_oce, ONLY : &   ! Geographical information
1210         & glamt,   &
1211         & gphit,   &
1212         & tmask
1213      !! * Arguments
1214      TYPE(obs_surf), INTENT(INOUT) :: logchldata     ! Full set of logchl data
1215      TYPE(obs_surf), INTENT(INOUT) :: logchldatqc    ! Subset of logchl data not failing screening
1216      LOGICAL :: ld_logchl     ! Switch for logchl data
1217      LOGICAL :: ld_nea        ! Switch for rejecting observation near land
1218      !! * Local declarations
1219      INTEGER :: iyea0         ! Initial date
1220      INTEGER :: imon0         !  - (year, month, day, hour, minute)
1221      INTEGER :: iday0   
1222      INTEGER :: ihou0   
1223      INTEGER :: imin0
1224      INTEGER :: icycle       ! Current assimilation cycle
1225                              ! Counters for observations that
1226      INTEGER :: iotdobs      !  - outside time domain
1227      INTEGER :: iosdsobs     !  - outside space domain
1228      INTEGER :: ilansobs     !  - within a model land cell
1229      INTEGER :: inlasobs     !  - close to land
1230      INTEGER :: igrdobs      !  - fail the grid search
1231                              ! Global counters for observations that
1232      INTEGER :: iotdobsmpp   !  - outside time domain
1233      INTEGER :: iosdsobsmpp  !  - outside space domain
1234      INTEGER :: ilansobsmpp  !  - within a model land cell
1235      INTEGER :: inlasobsmpp  !  - close to land
1236      INTEGER :: igrdobsmpp   !  - fail the grid search
1237      LOGICAL, DIMENSION(:), ALLOCATABLE :: & 
1238         & llvalid            ! data selection
1239      INTEGER :: jobs         ! Obs. loop variable
1240      INTEGER :: jstp         ! Time loop variable
1241      INTEGER :: inrc         ! Time index variable
1242
1243      IF (lwp) WRITE(numout,*)'obs_pre_logchl : Preparing the logchl observations...'
1244
1245      ! Initial date initialization (year, month, day, hour, minute)
1246      iyea0 =   ndate0 / 10000
1247      imon0 = ( ndate0 - iyea0 * 10000 ) / 100
1248      iday0 =   ndate0 - iyea0 * 10000 - imon0 * 100
1249      ihou0 = 0
1250      imin0 = 0
1251
1252      icycle = no     ! Assimilation cycle
1253
1254      ! Diagnostics counters for various failures.
1255
1256      iotdobs  = 0
1257      igrdobs  = 0
1258      iosdsobs = 0
1259      ilansobs = 0
1260      inlasobs = 0
1261
1262      ! -----------------------------------------------------------------------
1263      ! Find time coordinate for logchl data
1264      ! -----------------------------------------------------------------------
1265
1266      CALL obs_coo_tim( icycle, &
1267         &              iyea0,   imon0,   iday0,   ihou0,   imin0,      &
1268         &              logchldata%nsurf,   logchldata%nyea, logchldata%nmon, &
1269         &              logchldata%nday,    logchldata%nhou, logchldata%nmin, &
1270         &              logchldata%nqc,     logchldata%mstp, iotdobs        )
1271      CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
1272      ! -----------------------------------------------------------------------
1273      ! Check for logchl data failing the grid search
1274      ! -----------------------------------------------------------------------
1275
1276      CALL obs_coo_grd( logchldata%nsurf,   logchldata%mi, logchldata%mj, &
1277         &              logchldata%nqc,     igrdobs                         )
1278      CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
1279
1280      ! -----------------------------------------------------------------------
1281      ! Check for land points.
1282      ! -----------------------------------------------------------------------
1283
1284      CALL obs_coo_spc_2d( logchldata%nsurf,                 &
1285         &                 jpi,             jpj,             &
1286         &                 logchldata%mi,   logchldata%mj,   & 
1287         &                 logchldata%rlam, logchldata%rphi, &
1288         &                 glamt,           gphit,           &
1289         &                 tmask(:,:,1),    logchldata%nqc,  &
1290         &                 iosdsobs,        ilansobs,        &
1291         &                 inlasobs,        ld_nea           ) 
1292         
1293      CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
1294      CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
1295      CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
1296
1297      ! -----------------------------------------------------------------------
1298      ! Copy useful data from the logchldata data structure to
1299      ! the logchldatqc data structure
1300      ! -----------------------------------------------------------------------
1301
1302      ! Allocate the selection arrays
1303
1304      ALLOCATE( llvalid(logchldata%nsurf) )
1305     
1306      ! We want all data which has qc flags <= 0
1307
1308      llvalid(:)  = ( logchldata%nqc(:)  <= 10 )
1309
1310      ! The actual copying
1311
1312      CALL obs_surf_compress( logchldata,     logchldatqc,       .TRUE.,  numout, &
1313         &                    lvalid=llvalid )
1314
1315      ! Dellocate the selection arrays
1316      DEALLOCATE( llvalid )
1317
1318      ! -----------------------------------------------------------------------
1319      ! Print information about what observations are left after qc
1320      ! -----------------------------------------------------------------------
1321
1322      ! Update the total observation counter array
1323     
1324      IF(lwp) THEN
1325         WRITE(numout,*)
1326         WRITE(numout,*) 'obs_pre_logchl :'
1327         WRITE(numout,*) '~~~~~~~~~~~'
1328         WRITE(numout,*)
1329         WRITE(numout,*) ' logchl data outside time domain                  = ', &
1330            &            iotdobsmpp
1331         WRITE(numout,*) ' Remaining logchl data that failed grid search    = ', &
1332            &            igrdobsmpp
1333         WRITE(numout,*) ' Remaining logchl data outside space domain       = ', &
1334            &            iosdsobsmpp
1335         WRITE(numout,*) ' Remaining logchl data at land points             = ', &
1336            &            ilansobsmpp
1337         IF (ld_nea) THEN
1338            WRITE(numout,*) ' Remaining logchl data near land points (removed) = ', &
1339               &            inlasobsmpp
1340         ELSE
1341            WRITE(numout,*) ' Remaining logchl data near land points (kept)    = ', &
1342               &            inlasobsmpp
1343         ENDIF
1344         WRITE(numout,*) ' logchl data accepted                             = ', &
1345            &            logchldatqc%nsurfmpp
1346
1347         WRITE(numout,*)
1348         WRITE(numout,*) ' Number of observations per time step :'
1349         WRITE(numout,*)
1350         WRITE(numout,1997)
1351         WRITE(numout,1998)
1352      ENDIF
1353     
1354      DO jobs = 1, logchldatqc%nsurf
1355         inrc = logchldatqc%mstp(jobs) + 2 - nit000
1356         logchldatqc%nsstp(inrc)  = logchldatqc%nsstp(inrc) + 1
1357      END DO
1358     
1359      CALL obs_mpp_sum_integers( logchldatqc%nsstp, logchldatqc%nsstpmpp, &
1360         &                       nitend - nit000 + 2 )
1361
1362      IF ( lwp ) THEN
1363         DO jstp = nit000 - 1, nitend
1364            inrc = jstp - nit000 + 2
1365            WRITE(numout,1999) jstp, logchldatqc%nsstpmpp(inrc)
1366         END DO
1367      ENDIF
1368
13691997  FORMAT(10X,'Time step',5X,'logchl data')
13701998  FORMAT(10X,'---------',5X,'------------')
13711999  FORMAT(10X,I9,5X,I17)
1372     
1373   END SUBROUTINE obs_pre_logchl
1374
1375   SUBROUTINE obs_pre_spm( spmdata, spmdatqc, ld_spm, ld_nea )
1376      !!----------------------------------------------------------------------
1377      !!                    ***  ROUTINE obs_pre_spm  ***
1378      !!
1379      !! ** Purpose : First level check and screening of spm observations
1380      !!
1381      !! ** Method  : First level check and screening of spm observations
1382      !!
1383      !! ** Action  :
1384      !!
1385      !! References :
1386      !!   
1387      !! History :
1388      !!----------------------------------------------------------------------
1389      !! * Modules used
1390      USE domstp              ! Domain: set the time-step
1391      USE par_oce             ! Ocean parameters
1392      USE dom_oce, ONLY : &   ! Geographical information
1393         & glamt,   &
1394         & gphit,   &
1395         & tmask
1396      !! * Arguments
1397      TYPE(obs_surf), INTENT(INOUT) :: spmdata     ! Full set of spm data
1398      TYPE(obs_surf), INTENT(INOUT) :: spmdatqc    ! Subset of spm data not failing screening
1399      LOGICAL :: ld_spm     ! Switch for spm data
1400      LOGICAL :: ld_nea        ! Switch for rejecting observation near land
1401      !! * Local declarations
1402      INTEGER :: iyea0         ! Initial date
1403      INTEGER :: imon0         !  - (year, month, day, hour, minute)
1404      INTEGER :: iday0   
1405      INTEGER :: ihou0   
1406      INTEGER :: imin0
1407      INTEGER :: icycle       ! Current assimilation cycle
1408                              ! Counters for observations that
1409      INTEGER :: iotdobs      !  - outside time domain
1410      INTEGER :: iosdsobs     !  - outside space domain
1411      INTEGER :: ilansobs     !  - within a model land cell
1412      INTEGER :: inlasobs     !  - close to land
1413      INTEGER :: igrdobs      !  - fail the grid search
1414                              ! Global counters for observations that
1415      INTEGER :: iotdobsmpp   !  - outside time domain
1416      INTEGER :: iosdsobsmpp  !  - outside space domain
1417      INTEGER :: ilansobsmpp  !  - within a model land cell
1418      INTEGER :: inlasobsmpp  !  - close to land
1419      INTEGER :: igrdobsmpp   !  - fail the grid search
1420      LOGICAL, DIMENSION(:), ALLOCATABLE :: & 
1421         & llvalid            ! data selection
1422      INTEGER :: jobs         ! Obs. loop variable
1423      INTEGER :: jstp         ! Time loop variable
1424      INTEGER :: inrc         ! Time index variable
1425
1426      IF (lwp) WRITE(numout,*)'obs_pre_spm : Preparing the spm observations...'
1427
1428      ! Initial date initialization (year, month, day, hour, minute)
1429      iyea0 =   ndate0 / 10000
1430      imon0 = ( ndate0 - iyea0 * 10000 ) / 100
1431      iday0 =   ndate0 - iyea0 * 10000 - imon0 * 100
1432      ihou0 = 0
1433      imin0 = 0
1434
1435      icycle = no     ! Assimilation cycle
1436
1437      ! Diagnostics counters for various failures.
1438
1439      iotdobs  = 0
1440      igrdobs  = 0
1441      iosdsobs = 0
1442      ilansobs = 0
1443      inlasobs = 0
1444
1445      ! -----------------------------------------------------------------------
1446      ! Find time coordinate for spm data
1447      ! -----------------------------------------------------------------------
1448
1449      CALL obs_coo_tim( icycle, &
1450         &              iyea0,   imon0,   iday0,   ihou0,   imin0,      &
1451         &              spmdata%nsurf,   spmdata%nyea, spmdata%nmon, &
1452         &              spmdata%nday,    spmdata%nhou, spmdata%nmin, &
1453         &              spmdata%nqc,     spmdata%mstp, iotdobs        )
1454      CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
1455      ! -----------------------------------------------------------------------
1456      ! Check for spm data failing the grid search
1457      ! -----------------------------------------------------------------------
1458
1459      CALL obs_coo_grd( spmdata%nsurf,   spmdata%mi, spmdata%mj, &
1460         &              spmdata%nqc,     igrdobs                         )
1461      CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
1462
1463      ! -----------------------------------------------------------------------
1464      ! Check for land points.
1465      ! -----------------------------------------------------------------------
1466
1467      CALL obs_coo_spc_2d( spmdata%nsurf,                 &
1468         &                 jpi,             jpj,             &
1469         &                 spmdata%mi,   spmdata%mj,   & 
1470         &                 spmdata%rlam, spmdata%rphi, &
1471         &                 glamt,           gphit,           &
1472         &                 tmask(:,:,1),    spmdata%nqc,  &
1473         &                 iosdsobs,        ilansobs,        &
1474         &                 inlasobs,        ld_nea           ) 
1475         
1476      CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
1477      CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
1478      CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
1479
1480      ! -----------------------------------------------------------------------
1481      ! Copy useful data from the spmdata data structure to
1482      ! the spmdatqc data structure
1483      ! -----------------------------------------------------------------------
1484
1485      ! Allocate the selection arrays
1486
1487      ALLOCATE( llvalid(spmdata%nsurf) )
1488     
1489      ! We want all data which has qc flags <= 0
1490
1491      llvalid(:)  = ( spmdata%nqc(:)  <= 10 )
1492
1493      ! The actual copying
1494
1495      CALL obs_surf_compress( spmdata,     spmdatqc,       .TRUE.,  numout, &
1496         &                    lvalid=llvalid )
1497
1498      ! Dellocate the selection arrays
1499      DEALLOCATE( llvalid )
1500
1501      ! -----------------------------------------------------------------------
1502      ! Print information about what observations are left after qc
1503      ! -----------------------------------------------------------------------
1504
1505      ! Update the total observation counter array
1506     
1507      IF(lwp) THEN
1508         WRITE(numout,*)
1509         WRITE(numout,*) 'obs_pre_spm :'
1510         WRITE(numout,*) '~~~~~~~~~~~'
1511         WRITE(numout,*)
1512         WRITE(numout,*) ' spm data outside time domain                  = ', &
1513            &            iotdobsmpp
1514         WRITE(numout,*) ' Remaining spm data that failed grid search    = ', &
1515            &            igrdobsmpp
1516         WRITE(numout,*) ' Remaining spm data outside space domain       = ', &
1517            &            iosdsobsmpp
1518         WRITE(numout,*) ' Remaining spm data at land points             = ', &
1519            &            ilansobsmpp
1520         IF (ld_nea) THEN
1521            WRITE(numout,*) ' Remaining spm data near land points (removed) = ', &
1522               &            inlasobsmpp
1523         ELSE
1524            WRITE(numout,*) ' Remaining spm data near land points (kept)    = ', &
1525               &            inlasobsmpp
1526         ENDIF
1527         WRITE(numout,*) ' spm data accepted                             = ', &
1528            &            spmdatqc%nsurfmpp
1529
1530         WRITE(numout,*)
1531         WRITE(numout,*) ' Number of observations per time step :'
1532         WRITE(numout,*)
1533         WRITE(numout,1997)
1534         WRITE(numout,1998)
1535      ENDIF
1536     
1537      DO jobs = 1, spmdatqc%nsurf
1538         inrc = spmdatqc%mstp(jobs) + 2 - nit000
1539         spmdatqc%nsstp(inrc)  = spmdatqc%nsstp(inrc) + 1
1540      END DO
1541     
1542      CALL obs_mpp_sum_integers( spmdatqc%nsstp, spmdatqc%nsstpmpp, &
1543         &                       nitend - nit000 + 2 )
1544
1545      IF ( lwp ) THEN
1546         DO jstp = nit000 - 1, nitend
1547            inrc = jstp - nit000 + 2
1548            WRITE(numout,1999) jstp, spmdatqc%nsstpmpp(inrc)
1549         END DO
1550      ENDIF
1551
15521997  FORMAT(10X,'Time step',5X,'spm data')
15531998  FORMAT(10X,'---------',5X,'------------')
15541999  FORMAT(10X,I9,5X,I17)
1555     
1556   END SUBROUTINE obs_pre_spm
1557
1558   SUBROUTINE obs_coo_tim( kcycle, &
1559      &                    kyea0,   kmon0,   kday0,   khou0,   kmin0,     &
1560      &                    kobsno,                                        &
1561      &                    kobsyea, kobsmon, kobsday, kobshou, kobsmin,   &
1562      &                    kobsqc,  kobsstp, kotdobs                      )
1563      !!----------------------------------------------------------------------
1564      !!                    ***  ROUTINE obs_coo_tim ***
1565      !!
1566      !! ** Purpose : Compute the number of time steps to the observation time.
1567      !!
1568      !! ** Method  : For time coordinates ( yea_obs, mon_obs, day_obs,
1569      !!              hou_obs, min_obs ), this routine locates the time step
1570      !!              that is closest to this time.
1571      !!
1572      !! ** Action  :
1573      !!
1574      !! References :
1575      !!   
1576      !! History :
1577      !!        !  1997-07  (A. Weaver) Original
1578      !!        !  2006-08  (A. Weaver) NEMOVAR migration
1579      !!        !  2006-10  (A. Weaver) Cleanup
1580      !!        !  2007-01  (K. Mogensen) Rewritten with loop
1581      !!        !  2010-05  (D. Lea) Fix in leap year calculation for NEMO vn3.2
1582      !!----------------------------------------------------------------------
1583      !! * Modules used
1584      USE dom_oce, ONLY : &  ! Geographical information
1585         & rdt
1586      USE phycst, ONLY : &   ! Physical constants
1587         & rday,  &             
1588         & rmmss, &             
1589         & rhhmm                       
1590      !! * Arguments
1591      INTEGER, INTENT(IN) :: kcycle     ! Current cycle
1592      INTEGER, INTENT(IN) :: kyea0      ! Initial date coordinates
1593      INTEGER, INTENT(IN) :: kmon0
1594      INTEGER, INTENT(IN) :: kday0 
1595      INTEGER, INTENT(IN) :: khou0
1596      INTEGER, INTENT(IN) :: kmin0
1597      INTEGER, INTENT(IN) :: kobsno     ! Number of observations
1598      INTEGER, INTENT(INOUT) :: kotdobs   ! Number of observations failing time check
1599      INTEGER, DIMENSION(kobsno), INTENT(IN ) :: &
1600         & kobsyea,  &      ! Observation time coordinates
1601         & kobsmon,  &
1602         & kobsday,  & 
1603         & kobshou,  &
1604         & kobsmin
1605      INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
1606         & kobsqc           ! Quality control flag
1607      INTEGER, DIMENSION(kobsno), INTENT(OUT) :: &
1608         & kobsstp          ! Number of time steps up to the
1609                            ! observation time
1610
1611      !! * Local declarations
1612      INTEGER :: jyea
1613      INTEGER :: jmon
1614      INTEGER :: jday
1615      INTEGER :: jobs
1616      INTEGER :: iyeastr
1617      INTEGER :: iyeaend
1618      INTEGER :: imonstr
1619      INTEGER :: imonend
1620      INTEGER :: idaystr
1621      INTEGER :: idayend 
1622      INTEGER :: iskip
1623      INTEGER :: idaystp
1624      REAL(KIND=wp) :: zminstp
1625      REAL(KIND=wp) :: zhoustp
1626      REAL(KIND=wp) :: zobsstp 
1627      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
1628 
1629      !-----------------------------------------------------------------------
1630      ! Initialization
1631      !-----------------------------------------------------------------------
1632
1633      ! Intialize the number of time steps per day
1634      idaystp = NINT( rday / rdt )
1635
1636      !---------------------------------------------------------------------
1637      ! Locate the model time coordinates for interpolation
1638      !---------------------------------------------------------------------
1639
1640      DO jobs = 1, kobsno
1641
1642         ! Initialize the time step counter
1643         kobsstp(jobs) = nit000 - 1
1644
1645         ! Flag if observation date is less than the initial date
1646
1647         IF ( ( kobsyea(jobs) < kyea0 )                   &
1648            & .OR. ( ( kobsyea(jobs) == kyea0 )           &
1649            &        .AND. ( kobsmon(jobs) <  kmon0 ) )   &
1650            & .OR. ( ( kobsyea(jobs) == kyea0 )           &
1651            &        .AND. ( kobsmon(jobs) == kmon0 )     &
1652            &        .AND. ( kobsday(jobs) <  kday0 ) )   &
1653            & .OR. ( ( kobsyea(jobs) == kyea0 )           &
1654            &        .AND. ( kobsmon(jobs) == kmon0 )     &
1655            &        .AND. ( kobsday(jobs) == kday0 )     &
1656            &        .AND. ( kobshou(jobs) <  khou0 ) )   &
1657            &  .OR. ( ( kobsyea(jobs) == kyea0 )          &
1658            &        .AND. ( kobsmon(jobs) == kmon0 )     &
1659            &        .AND. ( kobsday(jobs) == kday0 )          &
1660            &        .AND. ( kobshou(jobs) == khou0 )          &
1661            &        .AND. ( kobsmin(jobs) <= kmin0 ) ) ) THEN
1662            kobsstp(jobs) = -1
1663            kobsqc(jobs)  = kobsqc(jobs) + 11
1664            kotdobs       = kotdobs + 1
1665            CYCLE
1666         ENDIF
1667
1668         ! Compute the number of time steps to the observation day
1669         iyeastr = kyea0
1670         iyeaend = kobsyea(jobs)
1671         
1672         !---------------------------------------------------------------------
1673         ! Year loop
1674         !---------------------------------------------------------------------
1675         DO jyea = iyeastr, iyeaend
1676
1677            CALL calc_month_len( jyea, imonth_len )
1678           
1679            imonstr = 1
1680            IF ( jyea == kyea0         ) imonstr = kmon0
1681            imonend = 12
1682            IF ( jyea == kobsyea(jobs) ) imonend = kobsmon(jobs)
1683           
1684            ! Month loop
1685            DO jmon = imonstr, imonend
1686               
1687               idaystr = 1
1688               IF (       ( jmon == kmon0   ) &
1689                  & .AND. ( jyea == kyea0   ) ) idaystr = kday0
1690               idayend = imonth_len(jmon)
1691               IF (       ( jmon == kobsmon(jobs) ) &
1692                  & .AND. ( jyea == kobsyea(jobs) ) ) idayend = kobsday(jobs) - 1
1693               
1694               ! Day loop
1695               DO jday = idaystr, idayend
1696                  kobsstp(jobs) = kobsstp(jobs) + idaystp
1697               END DO
1698               
1699            END DO
1700
1701         END DO
1702
1703         ! Add in the number of time steps to the observation minute
1704         zminstp = rmmss / rdt
1705         zhoustp = rhhmm * zminstp
1706
1707         zobsstp =   REAL( kobsmin(jobs) - kmin0, KIND=wp ) * zminstp &
1708            &      + REAL( kobshou(jobs) - khou0, KIND=wp ) * zhoustp
1709         kobsstp(jobs) = kobsstp(jobs) + NINT( zobsstp )
1710
1711         ! Flag if observation step outside the time window
1712         IF ( ( kobsstp(jobs) < ( nit000 - 1 ) ) &
1713            & .OR.( kobsstp(jobs) > nitend ) ) THEN
1714            kobsqc(jobs) = kobsqc(jobs) + 12
1715            kotdobs = kotdobs + 1
1716            CYCLE
1717         ENDIF
1718
1719      END DO
1720
1721   END SUBROUTINE obs_coo_tim
1722
1723   SUBROUTINE calc_month_len( iyear, imonth_len )
1724      !!----------------------------------------------------------------------
1725      !!                    ***  ROUTINE calc_month_len  ***
1726      !!         
1727      !! ** Purpose : Compute the number of days in a months given a year.
1728      !!
1729      !! ** Method  :
1730      !!
1731      !! ** Action  :
1732      !!
1733      !! History :
1734      !!        !  10-05  (D. Lea)   New routine based on day_init
1735      !!----------------------------------------------------------------------
1736
1737      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
1738      INTEGER :: iyear         !: year
1739     
1740      ! length of the month of the current year (from nleapy, read in namelist)
1741      IF ( nleapy < 2 ) THEN
1742         imonth_len(:) = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)
1743         IF ( nleapy == 1 ) THEN   ! we are using calendar with leap years
1744            IF ( MOD(iyear, 4) == 0 .AND. ( MOD(iyear, 400) == 0 .OR. MOD(iyear, 100) /= 0 ) ) THEN
1745               imonth_len(2) = 29
1746            ENDIF
1747         ENDIF
1748      ELSE
1749         imonth_len(:) = nleapy   ! all months with nleapy days per year
1750      ENDIF
1751
1752   END SUBROUTINE
1753
1754   SUBROUTINE obs_coo_tim_prof( kcycle,                                   &
1755      &                    kyea0,   kmon0,   kday0,   khou0,   kmin0,     &
1756      &                    kobsno,                                        &
1757      &                    kobsyea, kobsmon, kobsday, kobshou, kobsmin,   &
1758      &                    ktyp,    kobsqc,  kobsstp, kotdobs, kdailyavtypes, &
1759      &                    ld_dailyav )
1760      !!----------------------------------------------------------------------
1761      !!                    ***  ROUTINE obs_coo_tim ***
1762      !!
1763      !! ** Purpose : Compute the number of time steps to the observation time.
1764      !!
1765      !! ** Method  : For time coordinates ( yea_obs, mon_obs, day_obs,
1766      !!              hou_obs, min_obs ), this routine locates the time step
1767      !!              that is closest to this time.
1768      !!
1769      !! ** Action  :
1770      !!
1771      !! References :
1772      !!   
1773      !! History :
1774      !!        !  1997-07  (A. Weaver) Original
1775      !!        !  2006-08  (A. Weaver) NEMOVAR migration
1776      !!        !  2006-10  (A. Weaver) Cleanup
1777      !!        !  2007-01  (K. Mogensen) Rewritten with loop
1778      !!----------------------------------------------------------------------
1779      !! * Modules used
1780      !! * Arguments
1781      INTEGER, INTENT(IN) :: kcycle      ! Current cycle
1782      INTEGER, INTENT(IN) :: kyea0       ! Initial date coordinates
1783      INTEGER, INTENT(IN) :: kmon0
1784      INTEGER, INTENT(IN) :: kday0
1785      INTEGER, INTENT(IN) :: khou0
1786      INTEGER, INTENT(IN) :: kmin0
1787      INTEGER, INTENT(IN) :: kobsno      ! Number of observations
1788      INTEGER, INTENT(INOUT) ::  kotdobs    ! Number of observations failing time check
1789      INTEGER, DIMENSION(kobsno), INTENT(IN ) :: &
1790         & kobsyea,  &      ! Observation time coordinates
1791         & kobsmon,  &
1792         & kobsday,  & 
1793         & kobshou,  &
1794         & kobsmin,  &
1795         & ktyp             ! Observation type.
1796      INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
1797         & kobsqc           ! Quality control flag
1798      INTEGER, DIMENSION(kobsno), INTENT(OUT) :: &
1799         & kobsstp          ! Number of time steps up to the
1800                            ! observation time
1801      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
1802         & kdailyavtypes    ! Types for daily averages
1803      LOGICAL, OPTIONAL :: ld_dailyav    ! All types are daily averages
1804      !! * Local declarations
1805      INTEGER :: jobs
1806
1807      !-----------------------------------------------------------------------
1808      ! Call standard obs_coo_tim
1809      !-----------------------------------------------------------------------
1810
1811      CALL obs_coo_tim( kcycle, &
1812         &              kyea0,   kmon0,   kday0,   khou0,   kmin0,     &
1813         &              kobsno,                                        &
1814         &              kobsyea, kobsmon, kobsday, kobshou, kobsmin,   &
1815         &              kobsqc,  kobsstp, kotdobs                      )
1816
1817      !------------------------------------------------------------------------
1818      ! Always reject daily averaged data (e.g. MRB data (820)) at initial time
1819      !------------------------------------------------------------------------
1820
1821      IF ( PRESENT(kdailyavtypes) ) THEN
1822         DO jobs = 1, kobsno
1823           
1824            IF ( kobsqc(jobs) <= 10 ) THEN
1825               
1826               IF ( ( kobsstp(jobs) == (nit000 - 1) ).AND.&
1827                  & ( ANY (kdailyavtypes(:) == ktyp(jobs)) ) ) THEN
1828                  kobsqc(jobs) = kobsqc(jobs) + 14
1829                  kotdobs      = kotdobs + 1
1830                  CYCLE
1831               ENDIF
1832               
1833            ENDIF
1834         END DO
1835      ENDIF
1836
1837      !------------------------------------------------------------------------
1838      ! If ld_dailyav is set then all data assumed to be daily averaged
1839      !------------------------------------------------------------------------
1840     
1841      IF ( PRESENT( ld_dailyav) ) THEN
1842         IF (ld_dailyav) THEN
1843            DO jobs = 1, kobsno
1844               
1845               IF ( kobsqc(jobs) <= 10 ) THEN
1846                 
1847                  IF ( kobsstp(jobs) == (nit000 - 1) ) THEN
1848                     kobsqc(jobs) = kobsqc(jobs) + 14
1849                     kotdobs      = kotdobs + 1
1850                     CYCLE
1851                  ENDIF
1852                 
1853               ENDIF
1854            END DO
1855         ENDIF
1856      ENDIF
1857
1858   END SUBROUTINE obs_coo_tim_prof
1859
1860   SUBROUTINE obs_coo_grd( kobsno, kobsi, kobsj, kobsqc, kgrdobs )
1861      !!----------------------------------------------------------------------
1862      !!                    ***  ROUTINE obs_coo_grd ***
1863      !!
1864      !! ** Purpose : Verify that the grid search has not failed
1865      !!
1866      !! ** Method  : The previously computed i,j indeces are checked 
1867      !!
1868      !! ** Action  :
1869      !!
1870      !! References :
1871      !!   
1872      !! History :
1873      !!        !  2007-01  (K. Mogensen)  Original
1874      !!----------------------------------------------------------------------
1875      !! * Arguments
1876      INTEGER, INTENT(IN) :: kobsno        ! Number of observations
1877      INTEGER, DIMENSION(kobsno), INTENT(IN ) :: &
1878         & kobsi, &         ! i,j indeces previously computed
1879         & kobsj
1880      INTEGER, INTENT(INOUT) ::  kgrdobs   ! Number of observations failing the check
1881      INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
1882         & kobsqc           ! Quality control flag
1883
1884      !! * Local declarations
1885      INTEGER :: jobs       ! Loop variable
1886
1887      ! Flag if the grid search failed
1888
1889      DO jobs = 1, kobsno
1890         IF ( ( kobsi(jobs) <= 0 ) .AND. ( kobsj(jobs) <= 0 ) ) THEN
1891            kobsqc(jobs) = kobsqc(jobs) + 18
1892            kgrdobs = kgrdobs + 1
1893         ENDIF
1894      END DO
1895     
1896   END SUBROUTINE obs_coo_grd
1897
1898   SUBROUTINE obs_coo_spc_2d( kobsno, kpi,     kpj,              &
1899      &                       kobsi,  kobsj,   pobslam, pobsphi, & 
1900      &                       plam,   pphi,    pmask,            &
1901      &                       kobsqc, kosdobs, klanobs,          &
1902      &                       knlaobs,ld_nea                     )
1903      !!----------------------------------------------------------------------
1904      !!                    ***  ROUTINE obs_coo_spc_2d  ***
1905      !!
1906      !! ** Purpose : Check for points outside the domain and land points
1907      !!
1908      !! ** Method  : Remove the observations that are outside the model space
1909      !!              and time domain or located within model land cells.
1910      !!
1911      !! ** Action  :
1912      !!   
1913      !! History :
1914      !!        !  2007-03  (A. Weaver, K. Mogensen)  Original
1915      !!        !  2007-06  (K. Mogensen et al) Reject obs. near land.
1916      !!----------------------------------------------------------------------
1917      !! * Modules used
1918
1919      !! * Arguments
1920      INTEGER, INTENT(IN) :: kobsno    ! Total number of observations
1921      INTEGER, INTENT(IN) :: kpi       ! Number of grid points in (i,j)
1922      INTEGER, INTENT(IN) :: kpj
1923      INTEGER, DIMENSION(kobsno), INTENT(IN) :: &
1924         & kobsi, &           ! Observation (i,j) coordinates
1925         & kobsj
1926      REAL(KIND=wp), DIMENSION(kobsno), INTENT(IN) :: &
1927         & pobslam, &         ! Observation (lon,lat) coordinates
1928         & pobsphi
1929      REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) :: &
1930         & plam, pphi         ! Model (lon,lat) coordinates
1931      REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) :: &
1932         & pmask              ! Land mask array
1933      INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
1934         & kobsqc             ! Observation quality control
1935      INTEGER, INTENT(INOUT) :: kosdobs   ! Observations outside space domain
1936      INTEGER, INTENT(INOUT) :: klanobs   ! Observations within a model land cell
1937      INTEGER, INTENT(INOUT) :: knlaobs   ! Observations near land
1938      LOGICAL, INTENT(IN) :: ld_nea       ! Flag observations near land
1939      !! * Local declarations
1940      REAL(KIND=wp), DIMENSION(2,2,kobsno) :: &
1941         & zgmsk              ! Grid mask
1942      REAL(KIND=wp), DIMENSION(2,2,kobsno) :: &
1943         & zglam, &           ! Model longitude at grid points
1944         & zgphi              ! Model latitude at grid points
1945      INTEGER, DIMENSION(2,2,kobsno) :: &
1946         & igrdi, &           ! Grid i,j
1947         & igrdj
1948      LOGICAL :: lgridobs           ! Is observation on a model grid point.
1949      INTEGER :: iig, ijg           ! i,j of observation on model grid point.
1950      INTEGER :: jobs, ji, jj
1951     
1952      ! Get grid point indices
1953
1954      DO jobs = 1, kobsno
1955         
1956         ! For invalid points use 2,2
1957
1958         IF ( kobsqc(jobs) >= 10 ) THEN
1959
1960            igrdi(1,1,jobs) = 1
1961            igrdj(1,1,jobs) = 1
1962            igrdi(1,2,jobs) = 1
1963            igrdj(1,2,jobs) = 2
1964            igrdi(2,1,jobs) = 2
1965            igrdj(2,1,jobs) = 1
1966            igrdi(2,2,jobs) = 2
1967            igrdj(2,2,jobs) = 2
1968
1969         ELSE
1970
1971            igrdi(1,1,jobs) = kobsi(jobs)-1
1972            igrdj(1,1,jobs) = kobsj(jobs)-1
1973            igrdi(1,2,jobs) = kobsi(jobs)-1
1974            igrdj(1,2,jobs) = kobsj(jobs)
1975            igrdi(2,1,jobs) = kobsi(jobs)
1976            igrdj(2,1,jobs) = kobsj(jobs)-1
1977            igrdi(2,2,jobs) = kobsi(jobs)
1978            igrdj(2,2,jobs) = kobsj(jobs)
1979
1980         ENDIF
1981
1982      END DO
1983     
1984      CALL obs_int_comm_2d( 2, 2, kobsno, igrdi, igrdj, pmask, zgmsk )
1985      CALL obs_int_comm_2d( 2, 2, kobsno, igrdi, igrdj, plam, zglam )
1986      CALL obs_int_comm_2d( 2, 2, kobsno, igrdi, igrdj, pphi, zgphi )
1987
1988      DO jobs = 1, kobsno
1989
1990         ! Skip bad observations
1991         IF ( kobsqc(jobs) >= 10 ) CYCLE
1992
1993         ! Flag if the observation falls outside the model spatial domain
1994         IF (       ( pobslam(jobs) < -180. ) &
1995            &  .OR. ( pobslam(jobs) >  180. ) &
1996            &  .OR. ( pobsphi(jobs) <  -90. ) &
1997            &  .OR. ( pobsphi(jobs) >   90. ) ) THEN
1998            kobsqc(jobs) = kobsqc(jobs) + 11
1999            kosdobs = kosdobs + 1
2000            CYCLE
2001         ENDIF
2002
2003         ! Flag if the observation falls with a model land cell
2004         IF ( SUM( zgmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN
2005            kobsqc(jobs) = kobsqc(jobs)  + 12
2006            klanobs = klanobs + 1
2007            CYCLE
2008         ENDIF
2009
2010         ! Check if this observation is on a grid point
2011
2012         lgridobs = .FALSE.
2013         iig = -1
2014         ijg = -1
2015         DO jj = 1, 2
2016            DO ji = 1, 2
2017               IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) &
2018                  & .AND. &
2019                  & ( ABS( zglam(ji,jj,jobs) - pobslam(jobs) ) < 1.0e-6_wp ) &
2020                  & ) THEN
2021                  lgridobs = .TRUE.
2022                  iig = ji
2023                  ijg = jj
2024               ENDIF
2025            END DO
2026         END DO
2027 
2028         ! For observations on the grid reject them if their are at
2029         ! a masked point
2030         
2031         IF (lgridobs) THEN
2032            IF (zgmsk(iig,ijg,jobs) == 0.0_wp ) THEN
2033               kobsqc(jobs) = kobsqc(jobs) + 12
2034               klanobs = klanobs + 1
2035               CYCLE
2036            ENDIF
2037         ENDIF
2038                     
2039         ! Flag if the observation falls is close to land
2040         IF ( MINVAL( zgmsk(1:2,1:2,jobs) ) == 0.0_wp) THEN
2041            IF (ld_nea) kobsqc(jobs) = kobsqc(jobs) + 14
2042            knlaobs = knlaobs + 1
2043            CYCLE
2044         ENDIF
2045           
2046      END DO
2047
2048   END SUBROUTINE obs_coo_spc_2d
2049
2050   SUBROUTINE obs_coo_spc_3d( kprofno, kobsno,  kpstart, kpend, &
2051      &                       kpi,     kpj,     kpk,            &
2052      &                       kobsi,   kobsj,   kobsk,          &
2053      &                       pobslam, pobsphi, pobsdep,        &
2054      &                       plam,    pphi,    pdep,    pmask, &
2055      &                       kpobsqc, kobsqc,  kosdobs,        &
2056      &                       klanobs, knlaobs, ld_nea          )
2057      !!----------------------------------------------------------------------
2058      !!                    ***  ROUTINE obs_coo_spc_3d  ***
2059      !!
2060      !! ** Purpose : Check for points outside the domain and land points
2061      !!              Reset depth of observation above highest model level
2062      !!              to the value of highest model level
2063      !!
2064      !! ** Method  : Remove the observations that are outside the model space
2065      !!              and time domain or located within model land cells.
2066      !!
2067      !!              NB. T and S profile observations lying between the ocean
2068      !!              surface and the depth of the first model T point are
2069      !!              assigned a depth equal to that of the first model T pt.
2070      !!
2071      !! ** Action  :
2072      !!   
2073      !! History :
2074      !!        !  2007-01  (K. Mogensen) Rewrite of parts of obs_scr
2075      !!        !  2007-06  (K. Mogensen et al) Reject obs. near land.
2076      !!----------------------------------------------------------------------
2077      !! * Modules used
2078      USE dom_oce, ONLY : &       ! Geographical information
2079         & gdepw_1d,      &
2080         & gdepw_0,       &                       
2081#if defined key_vvl
2082         & gdepw_n,       & 
2083         & gdept_n,       &
2084#endif
2085         & ln_zco,        &
2086         & ln_zps,        &
2087         & lk_vvl                       
2088
2089      !! * Arguments
2090      INTEGER, INTENT(IN) :: kprofno      ! Number of profiles
2091      INTEGER, INTENT(IN) :: kobsno       ! Total number of observations
2092      INTEGER, INTENT(IN) :: kpi          ! Number of grid points in (i,j,k)
2093      INTEGER, INTENT(IN) :: kpj
2094      INTEGER, INTENT(IN) :: kpk   
2095      INTEGER, DIMENSION(kprofno), INTENT(IN) :: &
2096         & kpstart, &         ! Start of individual profiles
2097         & kpend              ! End of individual profiles
2098      INTEGER, DIMENSION(kprofno), INTENT(IN) :: &
2099         & kobsi, &           ! Observation (i,j) coordinates
2100         & kobsj
2101      INTEGER, DIMENSION(kobsno), INTENT(IN) :: &
2102         & kobsk              ! Observation k coordinate
2103      REAL(KIND=wp), DIMENSION(kprofno), INTENT(IN) :: &
2104         & pobslam, &         ! Observation (lon,lat) coordinates
2105         & pobsphi
2106      REAL(KIND=wp), DIMENSION(kobsno), INTENT(INOUT) :: &
2107         & pobsdep            ! Observation depths 
2108      REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) :: &
2109         & plam, pphi         ! Model (lon,lat) coordinates
2110      REAL(KIND=wp), DIMENSION(kpk), INTENT(IN) :: &
2111         & pdep               ! Model depth coordinates
2112      REAL(KIND=wp), DIMENSION(kpi,kpj,kpk), INTENT(IN) :: &
2113         & pmask              ! Land mask array
2114      INTEGER, DIMENSION(kprofno), INTENT(INOUT) :: &
2115         & kpobsqc             ! Profile quality control
2116      INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
2117         & kobsqc             ! Observation quality control
2118      INTEGER, INTENT(INOUT) :: kosdobs     ! Observations outside space domain
2119      INTEGER, INTENT(INOUT) :: klanobs     ! Observations within a model land cell
2120      INTEGER, INTENT(INOUT) :: knlaobs     ! Observations near land
2121      LOGICAL, INTENT(IN) :: ld_nea         ! Flag observations near land
2122      !! * Local declarations
2123      REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: &
2124         & zgmsk              ! Grid mask
2125      REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: &
2126         & zgdepw         
2127      REAL(KIND=wp), DIMENSION(2,2,kprofno) :: &
2128         & zglam, &           ! Model longitude at grid points
2129         & zgphi              ! Model latitude at grid points
2130      INTEGER, DIMENSION(2,2,kprofno) :: &
2131         & igrdi, &           ! Grid i,j
2132         & igrdj
2133      LOGICAL :: lgridobs           ! Is observation on a model grid point.
2134      LOGICAL :: ll_next_to_land    ! Is a profile next to land
2135      INTEGER :: iig, ijg           ! i,j of observation on model grid point.
2136      INTEGER :: jobs, jobsp, jk, ji, jj
2137
2138      ! Get grid point indices
2139     
2140      DO jobs = 1, kprofno
2141
2142         ! For invalid points use 2,2
2143
2144         IF ( kpobsqc(jobs) >= 10 ) THEN
2145           
2146            igrdi(1,1,jobs) = 1
2147            igrdj(1,1,jobs) = 1
2148            igrdi(1,2,jobs) = 1
2149            igrdj(1,2,jobs) = 2
2150            igrdi(2,1,jobs) = 2
2151            igrdj(2,1,jobs) = 1
2152            igrdi(2,2,jobs) = 2
2153            igrdj(2,2,jobs) = 2
2154           
2155         ELSE
2156           
2157            igrdi(1,1,jobs) = kobsi(jobs)-1
2158            igrdj(1,1,jobs) = kobsj(jobs)-1
2159            igrdi(1,2,jobs) = kobsi(jobs)-1
2160            igrdj(1,2,jobs) = kobsj(jobs)
2161            igrdi(2,1,jobs) = kobsi(jobs)
2162            igrdj(2,1,jobs) = kobsj(jobs)-1
2163            igrdi(2,2,jobs) = kobsi(jobs)
2164            igrdj(2,2,jobs) = kobsj(jobs)
2165           
2166         ENDIF
2167         
2168      END DO
2169     
2170      CALL obs_int_comm_3d( 2, 2, kprofno, kpk, igrdi, igrdj, pmask, zgmsk )
2171      CALL obs_int_comm_2d( 2, 2, kprofno, igrdi, igrdj, plam, zglam )
2172      CALL obs_int_comm_2d( 2, 2, kprofno, igrdi, igrdj, pphi, zgphi )
2173      ! Need to know the bathy depth for each observation for sco
2174      CALL obs_int_comm_3d( 2, 2, kprofno, kpk, igrdi, igrdj, fsdepw(:,:,:), &
2175      &                     zgdepw )
2176
2177      DO jobs = 1, kprofno
2178
2179         ! Skip bad profiles
2180         IF ( kpobsqc(jobs) >= 10 ) CYCLE
2181
2182         ! Check if this observation is on a grid point
2183
2184         lgridobs = .FALSE.
2185         iig = -1
2186         ijg = -1
2187         DO jj = 1, 2
2188            DO ji = 1, 2
2189               IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) &
2190                  & .AND. &
2191                  & ( ABS( zglam(ji,jj,jobs) - pobslam(jobs) ) < 1.0e-6_wp ) &
2192                  & ) THEN
2193                  lgridobs = .TRUE.
2194                  iig = ji
2195                  ijg = jj
2196               ENDIF
2197            END DO
2198         END DO
2199
2200         ! Check if next to land
2201         IF (  ANY( zgmsk(1:2,1:2,1,jobs) == 0.0_wp ) ) THEN
2202            ll_next_to_land=.TRUE. 
2203         ELSE
2204            ll_next_to_land=.FALSE. 
2205         ENDIF 
2206         
2207         ! Reject observations
2208
2209         DO jobsp = kpstart(jobs), kpend(jobs)
2210
2211            ! Flag if the observation falls outside the model spatial domain
2212            IF (       ( pobslam(jobs) < -180.         )       &
2213               &  .OR. ( pobslam(jobs) >  180.         )       &
2214               &  .OR. ( pobsphi(jobs) <  -90.         )       &
2215               &  .OR. ( pobsphi(jobs) >   90.         )       &
2216               &  .OR. ( pobsdep(jobsp) < 0.0          )       &
2217               &  .OR. ( pobsdep(jobsp) > gdepw_1d(kpk)) ) THEN
2218               kobsqc(jobsp) = kobsqc(jobsp) + 11
2219               kosdobs = kosdobs + 1
2220               CYCLE
2221            ENDIF
2222
2223            ! To check if an observations falls within land there are two cases:
2224            ! 1: z-coordibnates, where the check uses the mask
2225            ! 2: terrain following (eg s-coordinates), 
2226            !    where we use the depth of the bottom cell to mask observations
2227             
2228            IF( (.NOT. lk_vvl) .AND. ( ln_zps .OR. ln_zco )  ) THEN !(CASE 1)
2229               
2230               ! Flag if the observation falls with a model land cell
2231               IF ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 
2232                  &  == 0.0_wp ) THEN
2233                  kobsqc(jobsp) = kobsqc(jobsp) + 12 
2234                  klanobs = klanobs + 1 
2235                  CYCLE
2236               ENDIF 
2237             
2238               ! Flag if the observation is close to land
2239               IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == & 
2240                  &  0.0_wp) THEN
2241                  knlaobs = knlaobs + 1 
2242                  IF (ld_nea) THEN   
2243                     kobsqc(jobsp) = kobsqc(jobsp) + 14 
2244                  ENDIF 
2245               ENDIF
2246             
2247            ELSE ! Case 2
2248 
2249               ! Flag if the observation is deeper than the bathymetry
2250               ! Or if it is within the mask
2251               IF ( ANY( zgdepw(1:2,1:2,kpk,jobs) < pobsdep(jobsp) ) &
2252                  &     .OR. & 
2253                  &  ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) &
2254                  &  == 0.0_wp) ) THEN
2255                  kobsqc(jobsp) = kobsqc(jobsp) + 12 
2256                  klanobs = klanobs + 1 
2257                  CYCLE
2258               ENDIF 
2259               
2260               ! Flag if the observation is close to land
2261               IF ( ll_next_to_land ) THEN
2262                  knlaobs = knlaobs + 1 
2263                  IF (ld_nea) THEN   
2264                     kobsqc(jobsp) = kobsqc(jobsp) + 14 
2265                  ENDIF 
2266               ENDIF
2267             
2268            ENDIF
2269
2270            ! For observations on the grid reject them if their are at
2271            ! a masked point
2272           
2273            IF (lgridobs) THEN
2274               IF (zgmsk(iig,ijg,kobsk(jobsp)-1,jobs) == 0.0_wp ) THEN
2275                  kobsqc(jobsp) = kobsqc(jobsp) + 12
2276                  klanobs = klanobs + 1
2277                  CYCLE
2278               ENDIF
2279            ENDIF
2280           
2281            ! Flag if the observation falls is close to land
2282            IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == &
2283               &  0.0_wp) THEN
2284               IF (ld_nea) kobsqc(jobsp) = kobsqc(jobsp) + 14
2285               knlaobs = knlaobs + 1
2286            ENDIF
2287
2288            ! Set observation depth equal to that of the first model depth
2289            IF ( pobsdep(jobsp) <= pdep(1) ) THEN
2290               pobsdep(jobsp) = pdep(1) 
2291            ENDIF
2292           
2293         END DO
2294      END DO
2295
2296   END SUBROUTINE obs_coo_spc_3d
2297
2298   SUBROUTINE obs_pro_rej( profdata )
2299      !!----------------------------------------------------------------------
2300      !!                    ***  ROUTINE obs_pro_rej ***
2301      !!
2302      !! ** Purpose : Reject all data within a rejected profile
2303      !!
2304      !! ** Method  :
2305      !!
2306      !! ** Action  :
2307      !!
2308      !! References :
2309      !!   
2310      !! History :
2311      !!        !  2007-10  (K. Mogensen) Original code
2312      !!----------------------------------------------------------------------
2313      !! * Modules used
2314      !! * Arguments
2315      TYPE(obs_prof), INTENT(INOUT) :: profdata     ! Profile data
2316      !! * Local declarations
2317      INTEGER :: jprof
2318      INTEGER :: jvar
2319      INTEGER :: jobs
2320     
2321      ! Loop over profiles
2322
2323      DO jprof = 1, profdata%nprof
2324
2325         IF ( profdata%nqc(jprof) > 10 ) THEN
2326           
2327            DO jvar = 1, profdata%nvar
2328
2329               DO jobs = profdata%npvsta(jprof,jvar),  &
2330                  &      profdata%npvend(jprof,jvar)
2331                 
2332                  profdata%var(jvar)%nvqc(jobs) = &
2333                     & profdata%var(jvar)%nvqc(jobs) + 26
2334
2335               END DO
2336
2337            END DO
2338
2339         ENDIF
2340
2341      END DO
2342
2343   END SUBROUTINE obs_pro_rej
2344
2345   SUBROUTINE obs_uv_rej( profdata, knumu, knumv )
2346      !!----------------------------------------------------------------------
2347      !!                    ***  ROUTINE obs_uv_rej ***
2348      !!
2349      !! ** Purpose : Reject u if v is rejected and vice versa
2350      !!
2351      !! ** Method  :
2352      !!
2353      !! ** Action  :
2354      !!
2355      !! References :
2356      !!   
2357      !! History :
2358      !!        !  2009-2  (K. Mogensen) Original code
2359      !!----------------------------------------------------------------------
2360      !! * Modules used
2361      !! * Arguments
2362      TYPE(obs_prof), INTENT(INOUT) :: profdata   ! Profile data
2363      INTEGER, INTENT(INOUT) :: knumu             ! Number of u rejected
2364      INTEGER, INTENT(INOUT) :: knumv             ! Number of v rejected
2365      !! * Local declarations
2366      INTEGER :: jprof
2367      INTEGER :: jvar
2368      INTEGER :: jobs
2369     
2370      ! Loop over profiles
2371
2372      DO jprof = 1, profdata%nprof
2373
2374         IF ( ( profdata%npvsta(jprof,1) /= profdata%npvsta(jprof,2) ) .OR. &
2375            & ( profdata%npvend(jprof,1) /= profdata%npvend(jprof,2) ) ) THEN
2376
2377            CALL ctl_stop('U,V profiles inconsistent in obs_uv_rej')
2378            RETURN
2379
2380         ENDIF
2381
2382         DO jobs = profdata%npvsta(jprof,1), profdata%npvend(jprof,1)
2383           
2384            IF ( ( profdata%var(1)%nvqc(jobs) > 10 ) .AND. &
2385               & ( profdata%var(2)%nvqc(jobs) <= 10) ) THEN
2386               profdata%var(2)%nvqc(jobs) = profdata%var(2)%nvqc(jobs) + 42
2387               knumv = knumv + 1
2388            ENDIF
2389            IF ( ( profdata%var(2)%nvqc(jobs) > 10 ) .AND. &
2390               & ( profdata%var(1)%nvqc(jobs) <= 10) ) THEN
2391               profdata%var(1)%nvqc(jobs) = profdata%var(1)%nvqc(jobs) + 42
2392               knumu = knumu + 1
2393            ENDIF
2394           
2395         END DO
2396           
2397      END DO
2398
2399   END SUBROUTINE obs_uv_rej
2400
2401END MODULE obs_prep
Note: See TracBrowser for help on using the repository browser.