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

source: branches/dev_1784_OBS/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 2001

Last change on this file since 2001 was 2001, checked in by djlea, 14 years ago

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