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_read_vel.F90 in branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_vel.F90

Last change on this file was 7367, checked in by deazer, 8 years ago

Contains merged code for CO5 reference.

File size: 27.5 KB
Line 
1MODULE obs_read_vel
2   !!======================================================================
3   !!                       ***  MODULE obs_read_vel  ***
4   !! Observation diagnostics: Read the velocity profile observations
5   !!======================================================================
6
7   !!----------------------------------------------------------------------
8   !!   obs_rea_vel_dri : Driver for reading profile obs
9   !!----------------------------------------------------------------------
10
11   !! * Modules used   
12   USE par_kind                 ! Precision variables
13   USE par_oce                  ! Ocean parameters
14   USE in_out_manager           ! I/O manager
15   USE dom_oce                  ! Ocean space and time domain variables
16   USE obs_mpp                  ! MPP support routines for observation diagnostics
17   USE julian                   ! Julian date routines
18   USE obs_utils                ! Observation operator utility functions
19   USE obs_prep                 ! Prepare observation arrays
20   USE obs_grid                 ! Grid search
21   USE obs_sort                 ! Sorting observation arrays
22   USE obs_profiles_def         ! Profile definitions
23   USE obs_conv                 ! Various conversion routines
24   USE obs_types                ! Observation type definitions
25   USE netcdf                   ! NetCDF library
26   USE obs_oper                 ! Observation operators
27   USE obs_vel_io               ! Velocity profile files I/O (non-FB files)
28   USE lib_mpp                  ! For ctl_warn/stop
29
30   IMPLICIT NONE
31
32   !! * Routine accessibility
33   PRIVATE
34
35   PUBLIC obs_rea_vel_dri  ! Read the profile observations
36
37   !!----------------------------------------------------------------------
38   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
39   !! $Id$
40   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
41   !!----------------------------------------------------------------------
42
43CONTAINS
44 
45   SUBROUTINE obs_rea_vel_dri( kformat, &
46      &                        profdata, knumfiles, cfilenames, &
47      &                        kvars, kextr, kstp, ddobsini, ddobsend, &
48      &                        ldignmis, ldavtimset, ldmod )
49      !!---------------------------------------------------------------------
50      !!
51      !!                   *** ROUTINE obs_rea_pro_dri ***
52      !!
53      !! ** Purpose : Read from file the profile observations
54      !!
55      !! ** Method  : Depending on kformat either ENACT, CORIOLIS or
56      !!              feedback data files are read
57      !!
58      !! ** Action  :
59      !!
60      !! References :
61      !!
62      !! History : 
63      !!      ! :  2009-01 (K. Mogensen) : New merged version of old routines
64      !!----------------------------------------------------------------------
65      !! * Modules used
66   
67      !! * Arguments
68      INTEGER :: kformat   ! Format of input data
69      !                    ! 1: ENACT
70      !                    ! 2: Coriolis
71      TYPE(obs_prof), INTENT(OUT) :: profdata    ! Profile data to be read
72      INTEGER, INTENT(IN) :: knumfiles           ! Number of files to read in
73      CHARACTER(LEN=128), INTENT(IN) ::  cfilenames(knumfiles) ! File names to read in
74      INTEGER, INTENT(IN) :: kvars       ! Number of variables in profdata
75      INTEGER, INTENT(IN) :: kextr       ! Number of extra fields for each var in profdata
76      INTEGER, INTENT(IN) :: kstp        !  Ocean time-step index
77      LOGICAL, INTENT(IN) :: ldignmis    ! Ignore missing files
78      LOGICAL, INTENT(IN) :: ldavtimset  ! Set time to be equal to the end of the day
79      LOGICAL, INTENT(IN) :: ldmod       ! Initialize model from input data
80      REAL(KIND=dp), INTENT(IN) :: ddobsini   ! Obs. ini time in YYYYMMDD.HHMMSS
81      REAL(KIND=dp), INTENT(IN) :: ddobsend   ! Obs. end time in YYYYMMDD.HHMMSS
82
83      !! * Local declarations
84      CHARACTER(LEN=15), PARAMETER :: cpname='obs_rea_vel_dri'
85      INTEGER :: jvar
86      INTEGER :: ji
87      INTEGER :: jj
88      INTEGER :: jk
89      INTEGER :: ij
90      INTEGER :: iflag
91      INTEGER :: inobf
92      INTEGER :: i_file_id
93      INTEGER :: inowin
94      INTEGER :: iyea
95      INTEGER :: imon
96      INTEGER :: iday
97      INTEGER :: ihou
98      INTEGER :: imin
99      INTEGER :: isec
100      INTEGER, DIMENSION(knumfiles) :: &
101         & irefdate
102      INTEGER, DIMENSION(ntyp1770+1) :: &
103         & itypuv,    &
104         & itypuvmpp 
105      INTEGER :: iuv3dtmpp
106      INTEGER, DIMENSION(:), ALLOCATABLE :: &
107         & iobsiu,   &
108         & iobsju,   &
109         & iprocu,   &
110         & iobsiv,   &
111         & iobsjv,   &
112         & iprocv,   &
113         & iindx,    &
114         & ifileidx, &
115         & iprofidx
116      INTEGER :: itype
117      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
118         & zphi, &
119         & zlam
120      REAL(dp), DIMENSION(:), ALLOCATABLE :: &
121         & zdat
122      LOGICAL :: &
123         & llvalprof
124      TYPE(obfbdata), POINTER, DIMENSION(:) :: &
125         & inpfiles
126      REAL(dp), DIMENSION(knumfiles) :: &
127         & djulini, &
128         & djulend
129      INTEGER :: iprof
130      INTEGER :: iproftot
131      INTEGER :: iuv3dt
132      INTEGER, DIMENSION(kvars) :: iv3dt
133      CHARACTER(len=8) :: cl_refdate
134   
135      ! Local initialization
136      iprof = 0
137      iuv3dt = 0
138
139      !-----------------------------------------------------------------------
140      ! Check data the model part is just with feedback data files
141      !-----------------------------------------------------------------------
142      IF ( ldmod .AND. ( kformat /= 0 ) ) THEN
143         CALL ctl_stop( 'Model can only be read from feedback data' )
144         RETURN
145      ENDIF
146
147      !-----------------------------------------------------------------------
148      ! Count the number of files needed and allocate the obfbdata type
149      !-----------------------------------------------------------------------
150     
151      inobf = knumfiles
152
153      ALLOCATE( inpfiles(inobf) )
154
155      prof_files : DO jj = 1, inobf
156         
157         !---------------------------------------------------------------------
158         ! Prints
159         !---------------------------------------------------------------------
160         IF(lwp) THEN
161            WRITE(numout,*)
162            WRITE(numout,*) ' obs_rea_vel_dri : Reading from file = ', &
163               & TRIM( TRIM( cfilenames(jj) ) )
164            WRITE(numout,*) ' ~~~~~~~~~~~~~~~'
165            WRITE(numout,*)
166         ENDIF
167
168         !---------------------------------------------------------------------
169         !  Initialization: Open file and get dimensions only
170         !---------------------------------------------------------------------
171         
172         iflag = nf90_open( TRIM( TRIM( cfilenames(jj) ) ), nf90_nowrite, &
173            &                      i_file_id )
174         
175         IF ( iflag /= nf90_noerr ) THEN
176
177            IF ( ldignmis ) THEN
178               inpfiles(jj)%nobs = 0
179               CALL ctl_warn( 'File ' // TRIM( TRIM( cfilenames(jj) ) ) // &
180                  &           ' not found' )
181            ELSE
182               CALL ctl_stop( 'File ' // TRIM( TRIM( cfilenames(jj) ) ) // &
183                  &           ' not found' )
184            ENDIF
185
186         ELSE 
187           
188            !------------------------------------------------------------------
189            !  Close the file since it is opened in read_proffile
190            !------------------------------------------------------------------
191           
192            iflag = nf90_close( i_file_id )
193
194            !------------------------------------------------------------------
195            !  Read the profile file into inpfiles
196            !------------------------------------------------------------------
197            IF ( kformat == 0 ) THEN
198               CALL init_obfbdata( inpfiles(jj) )
199               IF(lwp) THEN
200                  WRITE(numout,*)
201                  WRITE(numout,*)'Reading from feedback file :', &
202                     &           TRIM( cfilenames(jj) )
203               ENDIF
204               CALL read_obfbdata( TRIM( cfilenames(jj) ), inpfiles(jj), &
205                  &                ldgrid = .TRUE. )
206               IF ( ldmod .AND. ( inpfiles(jj)%nadd == 0 ) ) THEN
207                  CALL ctl_stop( 'Model not in input data' )
208                  RETURN
209               ENDIF
210            ELSEIF ( kformat == 1 ) THEN
211               CALL read_taondbc( TRIM( cfilenames(jj) ), inpfiles(jj), &
212                  &               numout, lwp, .TRUE. )
213            ELSE
214               CALL ctl_stop( 'File format unknown' )
215            ENDIF
216           
217            !------------------------------------------------------------------
218            !  Change longitude (-180,180)
219            !------------------------------------------------------------------
220
221            DO ji = 1, inpfiles(jj)%nobs 
222
223               IF ( inpfiles(jj)%plam(ji) < -180. ) &
224                  &   inpfiles(jj)%plam(ji) = inpfiles(jj)%plam(ji) + 360.
225
226               IF ( inpfiles(jj)%plam(ji) >  180. ) &
227                  &   inpfiles(jj)%plam(ji) = inpfiles(jj)%plam(ji) - 360.
228
229            END DO
230
231            !------------------------------------------------------------------
232            !  Calculate the date  (change eventually)
233            !------------------------------------------------------------------
234            cl_refdate=inpfiles(jj)%cdjuldref(1:8)
235            READ(cl_refdate,'(I8)') irefdate(jj)
236           
237            CALL ddatetoymdhms( ddobsini, iyea, imon, iday, ihou, imin, isec )
238            CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulini(jj), &
239               &           krefdate = irefdate(jj) )
240            CALL ddatetoymdhms( ddobsend, iyea, imon, iday, ihou, imin, isec )
241            CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulend(jj), &
242               &           krefdate = irefdate(jj) )
243
244            IF ( ldavtimset ) THEN
245               DO ji = 1, inpfiles(jj)%nobs
246                  !
247                  !  for daily averaged data force the time
248                  !  to be the  end of the day
249                  !
250                  inpfiles(jj)%ptim(ji) = &
251                     & INT(inpfiles(jj)%ptim(ji)) + 1
252               END DO
253            ENDIF
254           
255            IF ( inpfiles(jj)%nobs > 0 ) THEN
256               inpfiles(jj)%iproc = -1
257               inpfiles(jj)%iobsi = -1
258               inpfiles(jj)%iobsj = -1
259            ENDIF
260            inowin = 0
261            DO ji = 1, inpfiles(jj)%nobs
262               IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
263               IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .OR. &
264                  & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE
265               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
266                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
267                  inowin = inowin + 1
268               ENDIF
269            END DO
270            ALLOCATE( zlam(inowin)  )
271            ALLOCATE( zphi(inowin)  )
272            ALLOCATE( iobsiu(inowin) )
273            ALLOCATE( iobsju(inowin) )
274            ALLOCATE( iprocu(inowin) )
275            ALLOCATE( iobsiv(inowin) )
276            ALLOCATE( iobsjv(inowin) )
277            ALLOCATE( iprocv(inowin) )
278            inowin = 0
279            DO ji = 1, inpfiles(jj)%nobs
280               IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
281               IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .OR. &
282                  & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE
283               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
284                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
285                  inowin = inowin + 1
286                  zlam(inowin) = inpfiles(jj)%plam(ji)
287                  zphi(inowin) = inpfiles(jj)%pphi(ji)
288               ENDIF
289            END DO
290
291            CALL obs_grid_search( inowin, zlam, zphi, iobsiu, iobsju, iprocu, &
292               & 'U' )
293            CALL obs_grid_search( inowin, zlam, zphi, iobsiv, iobsjv, iprocv, &
294               & 'V' )
295
296            ! Check that grid search has not failed for one component !
297            ! and not the other
298
299            DO ji = 1, inowin
300               IF ( ( iprocu(ji) < 0 ) .AND. ( iprocv(ji) >= 0 ) ) THEN
301                  IF (lwp) THEN
302                     WRITE(numout,*)
303                     WRITE(numout,'(1X,A,2F14.4)') &
304                        & 'Grid search for u failed at ', &
305                        & zphi(ji),zlam(ji)
306                     WRITE(numout,*)'Changing v grid search status to failed.'
307                  ENDIF
308                  iprocv(ji) = -1
309                  iobsiv(ji) = -1
310                  iobsjv(ji) = -1
311               ENDIF
312               IF ( ( iprocv(ji) < 0 ) .AND. ( iprocu(ji) >= 0 ) ) THEN
313                  IF (lwp) THEN
314                     WRITE(numout,*)
315                     WRITE(numout,'(1X,A,2F14.4)') &
316                        & 'Grid search for v failed at ', &
317                        & zphi(ji),zlam(ji)
318                     WRITE(numout,*)'Changing u grid search status to failed.'
319                  ENDIF
320                  iprocu(ji) = -1
321                  iobsiu(ji) = -1
322                  iobsju(ji) = -1
323               ENDIF
324            ENDDO
325
326            inowin = 0
327            DO ji = 1, inpfiles(jj)%nobs
328               IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
329               IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .OR. &
330                  & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE
331               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
332                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
333                  inowin = inowin + 1
334                  inpfiles(jj)%iproc(ji,1) = iprocu(inowin)
335                  inpfiles(jj)%iobsi(ji,1) = iobsiu(inowin)
336                  inpfiles(jj)%iobsj(ji,1) = iobsju(inowin)
337                  inpfiles(jj)%iproc(ji,2) = iprocv(inowin)
338                  inpfiles(jj)%iobsi(ji,2) = iobsiv(inowin)
339                  inpfiles(jj)%iobsj(ji,2) = iobsjv(inowin)
340                  IF ( inpfiles(jj)%iproc(ji,1) /= &
341                     & inpfiles(jj)%iproc(ji,2) ) THEN
342                     CALL ctl_stop( 'Error in obs_read_vel:', &
343                        & 'U and V observation on different processors')
344                  ENDIF
345               ENDIF
346            END DO
347            DEALLOCATE( zlam, zphi, iobsiu, iobsju, iprocu,  &
348               & iobsiv, iobsjv, iprocv )
349
350            DO ji = 1, inpfiles(jj)%nobs
351               IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
352               IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .OR. &
353                  & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE
354               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
355                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
356                  IF ( nproc == 0 ) THEN
357                     IF ( inpfiles(jj)%iproc(ji,1) >  nproc ) CYCLE
358                  ELSE
359                     IF ( inpfiles(jj)%iproc(ji,1) /=  nproc ) CYCLE
360                  ENDIF
361                  llvalprof = .FALSE.
362                  IF ( ( ( nproc == 0 ) .AND. &
363                     & ( inpfiles(jj)%iproc(ji,1) <=  nproc ) ) .OR. &
364                     & ( inpfiles(jj)%iproc(ji,1) ==  nproc ) ) THEN
365                     loop_uv_count : DO ij = 1,inpfiles(jj)%nlev
366                        IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
367                           & CYCLE
368                        IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
369                           & ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
370                           & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN
371                           iuv3dt = iuv3dt + 1
372                           llvalprof = .TRUE.
373                        ENDIF
374                     END DO loop_uv_count
375                  ENDIF
376                  IF ( llvalprof ) iprof = iprof + 1
377               ENDIF
378            END DO
379
380         ENDIF
381         
382      END DO prof_files
383
384      !-----------------------------------------------------------------------
385      ! Get the time ordered indices of the input data
386      !-----------------------------------------------------------------------
387
388      !---------------------------------------------------------------------
389      !  Loop over input data files to count total number of profiles
390      !---------------------------------------------------------------------
391      iproftot = 0
392      DO jj = 1, inobf
393         DO ji = 1, inpfiles(jj)%nobs
394            IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
395            IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .OR. &
396               & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE
397            IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
398               & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
399               iproftot = iproftot + 1
400            ENDIF
401         END DO
402      END DO
403
404      ALLOCATE( iindx(iproftot), ifileidx(iproftot), &
405         &      iprofidx(iproftot), zdat(iproftot) )
406      jk = 0
407      DO jj = 1, inobf
408         DO ji = 1, inpfiles(jj)%nobs
409            IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
410            IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .OR. &
411               & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE
412            IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
413               & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
414               jk = jk + 1
415               ifileidx(jk) = jj
416               iprofidx(jk) = ji
417               zdat(jk)     = inpfiles(jj)%ptim(ji)
418            ENDIF
419         END DO
420      END DO
421      CALL sort_dp_indx( iproftot, &
422         &               zdat,     &
423         &               iindx   )
424     
425      iv3dt(:) = -1
426      iv3dt(1:2) = iuv3dt
427      CALL obs_prof_alloc( profdata, kvars, kextr, iprof, iv3dt, &
428         &                 kstp, jpi, jpj, jpk )
429     
430      ! * Read obs/positions, QC, all variable and assign to profdata
431
432      profdata%nprof     = 0
433      profdata%nvprot(:) = 0
434
435      iprof = 0
436
437      iuv3dt = 0
438      itypuv   (:) = 0
439      itypuvmpp(:) = 0
440     
441      DO jk = 1, iproftot
442         
443         jj = ifileidx(iindx(jk))
444         ji = iprofidx(iindx(jk))
445
446         IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
447         IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .OR. &
448            & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE
449
450         IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND.  &
451            & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
452           
453            IF ( nproc == 0 ) THEN
454               IF ( inpfiles(jj)%iproc(ji,1) >  nproc ) CYCLE
455            ELSE
456               IF ( inpfiles(jj)%iproc(ji,1) /= nproc ) CYCLE
457            ENDIF
458           
459            llvalprof = .FALSE.
460
461            loop_prof : DO ij = 1, inpfiles(jj)%nlev
462               
463               IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
464                  & CYCLE
465               
466               IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
467                  & ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
468                  & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN
469
470                  llvalprof = .TRUE.
471                  EXIT loop_prof
472
473               ENDIF
474
475            END DO loop_prof
476           
477            ! Set profile information
478           
479            IF ( llvalprof ) THEN
480               
481               iprof = iprof + 1
482
483               CALL jul2greg( isec,                   &
484                  &           imin,                   &
485                  &           ihou,                   &
486                  &           iday,                   &
487                  &           imon,                   &
488                  &           iyea,                   &
489                  &           inpfiles(jj)%ptim(ji), &
490                  &           irefdate(jj) )
491
492
493               ! Profile time coordinates
494               profdata%nyea(iprof) = iyea
495               profdata%nmon(iprof) = imon
496               profdata%nday(iprof) = iday
497               profdata%nhou(iprof) = ihou
498               profdata%nmin(iprof) = imin
499               
500               ! Profile space coordinates
501               profdata%rlam(iprof) = inpfiles(jj)%plam(ji)
502               profdata%rphi(iprof) = inpfiles(jj)%pphi(ji)
503
504               ! Coordinate search parameters
505               profdata%mi  (iprof,1) = inpfiles(jj)%iobsi(ji,1)
506               profdata%mj  (iprof,1) = inpfiles(jj)%iobsj(ji,1)
507               profdata%mi  (iprof,2) = inpfiles(jj)%iobsi(ji,2)
508               profdata%mj  (iprof,2) = inpfiles(jj)%iobsj(ji,2)
509               
510               ! Profile WMO number
511               profdata%cwmo(iprof) = inpfiles(jj)%cdwmo(ji)
512               
513               ! Instrument type
514               READ( inpfiles(jj)%cdtyp(ji), '(I4)' ) itype
515               profdata%ntyp(iprof) = itype
516               
517               ! QC stuff
518
519               profdata%nqc(iprof)     = inpfiles(jj)%ioqc(ji)
520               profdata%nqcf(:,iprof)  = inpfiles(jj)%ioqcf(:,ji)
521               profdata%ipqc(iprof)    = inpfiles(jj)%ipqc(ji)
522               profdata%ipqcf(:,iprof) = inpfiles(jj)%ipqcf(:,ji)
523               profdata%itqc(iprof)    = inpfiles(jj)%itqc(ji)
524               profdata%itqcf(:,iprof) = inpfiles(jj)%itqcf(:,ji)
525               profdata%ivqc(iprof,:)  = inpfiles(jj)%ivqc(ji,:)
526               profdata%ivqcf(:,iprof,:) = inpfiles(jj)%ivqcf(:,ji,:)
527
528               ! Bookkeeping data to match profiles
529               profdata%npidx(iprof) = iprof
530               profdata%npfil(iprof) = iindx(jk)
531
532               ! Observation QC flag (whole profile)
533               profdata%nqc(iprof)  = 0 !TODO
534
535               loop_uv : DO ij = 1, inpfiles(jj)%nlev           
536                   
537                  IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
538                     & CYCLE
539
540                  IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
541                     & ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
542                     & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN
543                     iuv3dt = iuv3dt + 1
544                  ELSE
545                     CYCLE
546                  ENDIF
547
548                  ! Depth of U observation
549                  profdata%var(1)%vdep(iuv3dt) = &
550                     &                inpfiles(jj)%pdep(ij,ji)
551                 
552                  ! Depth of U observation QC
553                  profdata%var(1)%idqc(iuv3dt) = &
554                     &                inpfiles(jj)%idqc(ij,ji)
555                 
556                  ! Depth of U observation QC flags
557                  profdata%var(1)%idqcf(:,iuv3dt) = &
558                     &                inpfiles(jj)%idqcf(:,ij,ji)
559                 
560                  ! Profile index
561                  profdata%var(1)%nvpidx(iuv3dt) = iprof
562                 
563                  ! Vertical index in original profile
564                  profdata%var(1)%nvlidx(iuv3dt) = ij
565
566                  ! Profile U value
567
568                  profdata%var(1)%vobs(iuv3dt) = &
569                     &                inpfiles(jj)%pob(ij,ji,1)
570                  IF ( ldmod ) THEN
571                     profdata%var(1)%vmod(iuv3dt) = &
572                        &                inpfiles(jj)%padd(ij,ji,1,1)
573                  ENDIF
574                 
575                  ! Profile U qc
576                  profdata%var(1)%nvqc(iuv3dt) = &
577                     & inpfiles(jj)%ivlqc(ij,ji,1)
578                 
579                  ! Profile U qc flags
580                  profdata%var(1)%nvqcf(:,iuv3dt) = &
581                     & inpfiles(jj)%ivlqcf(:,ij,ji,1)
582                 
583                 
584                  ! Depth of V observation
585                  profdata%var(2)%vdep(iuv3dt) = &
586                     &                inpfiles(jj)%pdep(ij,ji)
587                 
588                  ! Depth of V observation QC
589                  profdata%var(2)%idqc(iuv3dt) = &
590                     &                inpfiles(jj)%idqc(ij,ji)
591                 
592                  ! Depth of V observation QC flags
593                  profdata%var(2)%idqcf(:,iuv3dt) = &
594                     &                inpfiles(jj)%idqcf(:,ij,ji)
595                 
596                  ! Profile index
597                  profdata%var(2)%nvpidx(iuv3dt) = iprof
598                 
599                  ! Vertical index in original profile
600                  profdata%var(2)%nvlidx(iuv3dt) = ij
601                 
602                  ! Profile V value
603                  profdata%var(2)%vobs(iuv3dt) = &
604                     &                inpfiles(jj)%pob(ij,ji,2)
605                  IF ( ldmod ) THEN
606                     profdata%var(2)%vmod(iuv3dt) = &
607                        &                inpfiles(jj)%padd(ij,ji,1,2)
608                  ENDIF
609                     
610                  ! Profile V qc
611                  profdata%var(2)%nvqc(iuv3dt) = &
612                     & inpfiles(jj)%ivlqc(ij,ji,2)
613                 
614                  ! Profile V qc flags
615                  profdata%var(2)%nvqcf(:,iuv3dt) = &
616                        & inpfiles(jj)%ivlqcf(:,ij,ji,2)
617                 
618                  ! Observation type
619                  itypuv( profdata%ntyp(iprof) + 1 ) = &
620                     & itypuv( profdata%ntyp(iprof) + 1 ) + 1
621
622
623               END DO loop_uv
624
625            ENDIF
626
627         ENDIF
628
629      END DO
630
631      !-----------------------------------------------------------------------
632      ! Sum up over processors
633      !-----------------------------------------------------------------------
634     
635      CALL obs_mpp_sum_integer ( iuv3dt, iuv3dtmpp )
636     
637      CALL obs_mpp_sum_integers( itypuv, itypuvmpp, ntyp1770 + 1 )
638     
639      !-----------------------------------------------------------------------
640      ! Output number of observations.
641      !-----------------------------------------------------------------------
642      IF(lwp) THEN
643         WRITE(numout,*) 
644         WRITE(numout,'(1X,A)') 'Profile U,V velocity data'
645         WRITE(numout,'(1X,A)') '-------------------------'
646         WRITE(numout,*) 
647         DO ji = 0, ntyp1770
648            IF ( itypuvmpp(ji+1) > 0 ) THEN
649               WRITE(numout,'(1X,A3,A3,I8)') ctypshort(ji), ' = ', &
650                  & itypuvmpp(ji+1)
651            ENDIF
652         END DO
653         WRITE(numout,'(1X,A)') '--------------'
654         WRITE(numout,'(1X,A6,I8)') &
655            & 'Total profile UV data                                 = ',&
656            & iuv3dtmpp
657         WRITE(numout,'(1X,A)') '--------------'
658      ENDIF
659     
660      profdata%nvprot(1)    = iuv3dt
661      profdata%nvprot(2)    = iuv3dt
662      profdata%nvprotmpp(1) = iuv3dtmpp
663      profdata%nvprotmpp(2) = iuv3dtmpp
664      profdata%nprof        = iprof
665
666      !-----------------------------------------------------------------------
667      ! Model level search
668      ! Only calculated here for z-levels and partial steps.
669      ! Otherwise calculated in obs_oper
670      !-----------------------------------------------------------------------
671      IF ( ln_zco .OR. ln_zps ) THEN
672         CALL obs_level_search( jpk, gdept_0, &
673            & profdata%nvprot(1), profdata%var(1)%vdep, &
674            & profdata%var(1)%mvk )
675         CALL obs_level_search( jpk, gdept_0, &
676            & profdata%nvprot(2), profdata%var(2)%vdep, &
677            & profdata%var(2)%mvk )
678      ELSE
679         profdata%var(1)%mvk = 0
680         profdata%var(2)%mvk = 0
681      ENDIF   
682     
683      !-----------------------------------------------------------------------
684      ! Set model equivalent to 99999
685      !-----------------------------------------------------------------------
686      IF ( .NOT. ldmod ) THEN
687         DO jvar = 1, kvars
688            profdata%var(jvar)%vmod(:) = fbrmdi
689         END DO
690      ENDIF
691      !-----------------------------------------------------------------------
692      ! Deallocate temporary data
693      !-----------------------------------------------------------------------
694      DEALLOCATE( ifileidx, iprofidx, zdat )
695
696      !-----------------------------------------------------------------------
697      ! Deallocate input data
698      !-----------------------------------------------------------------------
699      DO jj = 1, inobf
700         CALL dealloc_obfbdata( inpfiles(jj) )
701      END DO
702      DEALLOCATE( inpfiles )
703
704   END SUBROUTINE obs_rea_vel_dri
705   
706END MODULE obs_read_vel
Note: See TracBrowser for help on using the repository browser.