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/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90 @ 5838

Last change on this file since 5838 was 5838, checked in by timgraham, 8 years ago

Cleared svn keywords to allow use with fcm make at vn3.6

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