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 NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/OCE/OBS – NEMO

source: NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/OCE/OBS/obs_prep.F90 @ 13463

Last change on this file since 13463 was 13463, checked in by andmirek, 4 years ago

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