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

source: branches/UKMO/dev_r5107_iceshelf_fw_input_coupled_model/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 5511

Last change on this file since 5511 was 5511, checked in by davestorkey, 9 years ago

UKMO/dev_r5107_iceshelf_fw_input_coupled_model branch: clear SVN keywords

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