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

source: branches/2017/dev_r8657_UKMO_OBSoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 8667

Last change on this file since 8667 was 8667, checked in by timgraham, 6 years ago

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