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

source: branches/UKMO/dev_r5518_obs_oper_update_kd490/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90 @ 11259

Last change on this file since 11259 was 11259, checked in by dford, 5 years ago

Check variable names are as expected.
File size: 31.2 KB
Line 
1MODULE obs_read_prof
2   !!======================================================================
3   !!                       ***  MODULE obs_read_prof  ***
4   !! Observation diagnostics: Read the T and S profile observations
5   !!======================================================================
6
7   !!----------------------------------------------------------------------
8   !!   obs_rea_pro_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 lib_mpp                  ! For ctl_warn/stop
28   USE obs_fbm                  ! Feedback routines
29
30   IMPLICIT NONE
31
32   !! * Routine accessibility
33   PRIVATE
34
35   PUBLIC obs_rea_prof  ! 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_prof( profdata, knumfiles, cdfilenames, &
46      &                     kvars, kextr, kstp, ddobsini, ddobsend, &
47      &                     ldvar, ldignmis, ldsatt, &
48      &                     ldmod, cdvars, kdailyavtypes )
49      !!---------------------------------------------------------------------
50      !!
51      !!                   *** ROUTINE obs_rea_prof ***
52      !!
53      !! ** Purpose : Read from file the profile observations
54      !!
55      !! ** Method  : Read feedback data in and transform to NEMO internal
56      !!              profile data structure
57      !!
58      !! ** Action  :
59      !!
60      !! References :
61      !!
62      !! History : 
63      !!      ! :  2009-09 (K. Mogensen) : New merged version of old routines
64      !!      ! :  2015-08 (M. Martin) : Merged profile and velocity routines
65      !!----------------------------------------------------------------------
66
67      !! * Arguments
68      TYPE(obs_prof), INTENT(OUT) :: &
69         & profdata                     ! Profile data to be read
70      INTEGER, INTENT(IN) :: knumfiles  ! Number of files to read
71      CHARACTER(LEN=128), INTENT(IN) ::  &
72         & cdfilenames(knumfiles)        ! File names to read in
73      INTEGER, INTENT(IN) :: kvars      ! Number of variables in profdata
74      INTEGER, INTENT(IN) :: kextr      ! Number of extra fields for each var
75      INTEGER, INTENT(IN) :: kstp       ! Ocean time-step index
76      LOGICAL, DIMENSION(kvars), INTENT(IN) :: ldvar     ! Observed variables switches
77      LOGICAL, INTENT(IN) :: ldignmis   ! Ignore missing files
78      LOGICAL, INTENT(IN) :: ldsatt     ! Compute salinity at all temperature points
79      LOGICAL, INTENT(IN) :: ldmod      ! Initialize model from input data
80      REAL(dp), INTENT(IN) :: ddobsini  ! Obs. ini time in YYYYMMDD.HHMMSS
81      REAL(dp), INTENT(IN) :: ddobsend  ! Obs. end time in YYYYMMDD.HHMMSS
82      CHARACTER(len=8), DIMENSION(kvars), INTENT(IN) :: cdvars
83      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
84         & kdailyavtypes                ! Types of daily average observations
85
86      !! * Local declarations
87      CHARACTER(LEN=15), PARAMETER :: cpname='obs_rea_prof'
88      CHARACTER(len=8) :: clrefdate
89      CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvarsin
90      INTEGER :: jvar
91      INTEGER :: ji
92      INTEGER :: jj
93      INTEGER :: jk
94      INTEGER :: ij
95      INTEGER :: iflag
96      INTEGER :: inobf
97      INTEGER :: i_file_id
98      INTEGER :: inowin
99      INTEGER :: iyea
100      INTEGER :: imon
101      INTEGER :: iday
102      INTEGER :: ihou
103      INTEGER :: imin
104      INTEGER :: isec
105      INTEGER :: iprof
106      INTEGER :: iproftot
107      INTEGER, DIMENSION(kvars) :: ivart0
108      INTEGER, DIMENSION(kvars) :: ivart
109      INTEGER :: ip3dt
110      INTEGER :: ios
111      INTEGER :: ioserrcount
112      INTEGER, DIMENSION(kvars) :: ivartmpp
113      INTEGER :: ip3dtmpp
114      INTEGER :: itype
115      INTEGER, DIMENSION(knumfiles) :: &
116         & irefdate
117      INTEGER, DIMENSION(ntyp1770+1,kvars) :: &
118         & itypvar,    &
119         & itypvarmpp
120      INTEGER, DIMENSION(:,:), ALLOCATABLE :: &
121         & iobsi,    &
122         & iobsj,    &
123         & iproc
124      INTEGER, DIMENSION(:), ALLOCATABLE :: &
125         & iindx,    &
126         & ifileidx, &
127         & iprofidx
128      INTEGER, DIMENSION(imaxavtypes) :: &
129         & idailyavtypes
130      INTEGER, DIMENSION(kvars) :: &
131         & iv3dt
132      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
133         & zphi, &
134         & zlam
135      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
136         & zdat
137      REAL(wp), DIMENSION(knumfiles) :: &
138         & djulini, &
139         & djulend
140      LOGICAL :: llvalprof
141      LOGICAL :: lldavtimset
142      LOGICAL :: llcycle
143      TYPE(obfbdata), POINTER, DIMENSION(:) :: &
144         & inpfiles
145
146      ! Local initialization
147      iprof = 0
148      ivart0(:) = 0
149      ip3dt = 0
150
151      ! Daily average types
152      lldavtimset = .FALSE.
153      IF ( PRESENT(kdailyavtypes) ) THEN
154         idailyavtypes(:) = kdailyavtypes(:)
155         IF ( ANY (idailyavtypes(:) /= -1) ) lldavtimset = .TRUE.
156      ELSE
157         idailyavtypes(:) = -1
158      ENDIF
159
160      !-----------------------------------------------------------------------
161      ! Count the number of files needed and allocate the obfbdata type
162      !-----------------------------------------------------------------------
163
164      inobf = knumfiles
165
166      ALLOCATE( inpfiles(inobf) )
167
168      prof_files : DO jj = 1, inobf
169
170         !---------------------------------------------------------------------
171         ! Prints
172         !---------------------------------------------------------------------
173         IF(lwp) THEN
174            WRITE(numout,*)
175            WRITE(numout,*) ' obs_rea_pro_dri : Reading from file = ', &
176               & TRIM( TRIM( cdfilenames(jj) ) )
177            WRITE(numout,*) ' ~~~~~~~~~~~~~~~'
178            WRITE(numout,*)
179         ENDIF
180
181         !---------------------------------------------------------------------
182         !  Initialization: Open file and get dimensions only
183         !---------------------------------------------------------------------
184
185         iflag = nf90_open( TRIM( cdfilenames(jj) ), nf90_nowrite, &
186            &                      i_file_id )
187
188         IF ( iflag /= nf90_noerr ) THEN
189
190            IF ( ldignmis ) THEN
191               inpfiles(jj)%nobs = 0
192               CALL ctl_warn( 'File ' // TRIM( cdfilenames(jj) ) // &
193                  &           ' not found' )
194            ELSE
195               CALL ctl_stop( 'File ' // TRIM( cdfilenames(jj) ) // &
196                  &           ' not found' )
197            ENDIF
198
199         ELSE 
200
201            !------------------------------------------------------------------
202            !  Close the file since it is opened in read_obfbdata
203            !------------------------------------------------------------------
204
205            iflag = nf90_close( i_file_id )
206
207            !------------------------------------------------------------------
208            !  Read the profile file into inpfiles
209            !------------------------------------------------------------------
210            CALL init_obfbdata( inpfiles(jj) )
211            CALL read_obfbdata( TRIM( cdfilenames(jj) ), inpfiles(jj), &
212               &                ldgrid = .TRUE. )
213
214            IF ( inpfiles(jj)%nvar /= kvars ) THEN
215               CALL ctl_stop( 'Feedback format error: ', &
216                  &           ' unexpected number of vars in profile file' )
217            ENDIF
218
219            IF ( ldmod .AND. ( inpfiles(jj)%nadd == 0 ) ) THEN
220               CALL ctl_stop( 'Model not in input data' )
221            ENDIF
222
223            IF ( jj == 1 ) THEN
224               ALLOCATE( clvarsin( inpfiles(jj)%nvar ) )
225               DO ji = 1, inpfiles(jj)%nvar
226                 clvarsin(ji) = inpfiles(jj)%cname(ji)
227                 IF ( clvarsin(ji) /= cdvars(ji) ) THEN
228                    CALL ctl_stop( 'Feedback file variables do not match', &
229                        &           ' expected variable names for this type' )
230                 ENDIF
231               END DO
232            ELSE
233               DO ji = 1, inpfiles(jj)%nvar
234                  IF ( inpfiles(jj)%cname(ji) /= clvarsin(ji) ) THEN
235                     CALL ctl_stop( 'Feedback file variables not consistent', &
236                        &           ' with previous files for this type' )
237                  ENDIF
238               END DO
239            ENDIF
240
241            !------------------------------------------------------------------
242            !  Change longitude (-180,180)
243            !------------------------------------------------------------------
244
245            DO ji = 1, inpfiles(jj)%nobs 
246
247               IF ( inpfiles(jj)%plam(ji) < -180. ) &
248                  &   inpfiles(jj)%plam(ji) = inpfiles(jj)%plam(ji) + 360.
249
250               IF ( inpfiles(jj)%plam(ji) >  180. ) &
251                  &   inpfiles(jj)%plam(ji) = inpfiles(jj)%plam(ji) - 360.
252
253            END DO
254
255            !------------------------------------------------------------------
256            !  Calculate the date  (change eventually)
257            !------------------------------------------------------------------
258            clrefdate=inpfiles(jj)%cdjuldref(1:8)
259            READ(clrefdate,'(I8)') irefdate(jj)
260
261            CALL ddatetoymdhms( ddobsini, iyea, imon, iday, ihou, imin, isec )
262            CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulini(jj), &
263               &           krefdate = irefdate(jj) )
264            CALL ddatetoymdhms( ddobsend, iyea, imon, iday, ihou, imin, isec )
265            CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulend(jj), &
266               &           krefdate = irefdate(jj) )
267
268            ioserrcount=0
269            IF ( lldavtimset ) THEN
270
271               IF ( ANY ( idailyavtypes(:) /= -1 ) .AND. lwp) THEN
272                  WRITE(numout,*)' Resetting time of daily averaged', &
273                     &           ' observations to the end of the day'
274               ENDIF
275
276               DO ji = 1, inpfiles(jj)%nobs
277                  READ( inpfiles(jj)%cdtyp(ji), '(I4)', IOSTAT = ios, ERR = 900 ) itype
278900               IF ( ios /= 0 ) THEN
279                     ! Set type to zero if there is a problem in the string conversion
280                     itype = 0
281                  ENDIF
282
283                  IF ( ANY ( idailyavtypes(:) == itype ) ) THEN
284                  !  for daily averaged data force the time
285                  !  to be the last time-step of the day, but still within the day.
286                     IF ( inpfiles(jj)%ptim(ji) >= 0. ) THEN
287                        inpfiles(jj)%ptim(ji) = &
288                           & INT(inpfiles(jj)%ptim(ji)) + 0.9999
289                     ELSE
290                        inpfiles(jj)%ptim(ji) = &
291                           & INT(inpfiles(jj)%ptim(ji)) - 0.0001
292                     ENDIF
293                  ENDIF
294
295               END DO
296
297            ENDIF
298
299            IF ( inpfiles(jj)%nobs > 0 ) THEN
300               inpfiles(jj)%iproc(:,:) = -1
301               inpfiles(jj)%iobsi(:,:) = -1
302               inpfiles(jj)%iobsj(:,:) = -1
303            ENDIF
304            inowin = 0
305            DO ji = 1, inpfiles(jj)%nobs
306               IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
307               llcycle = .TRUE.
308               DO jvar = 1, kvars
309                  IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
310                     llcycle = .FALSE.
311                     EXIT
312                  ENDIF
313               END DO
314               IF ( llcycle ) CYCLE
315               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
316                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
317                  inowin = inowin + 1
318               ENDIF
319            END DO
320            ALLOCATE( zlam(inowin)  )
321            ALLOCATE( zphi(inowin)  )
322            ALLOCATE( iobsi(inowin,kvars) )
323            ALLOCATE( iobsj(inowin,kvars) )
324            ALLOCATE( iproc(inowin,kvars) )
325            inowin = 0
326            DO ji = 1, inpfiles(jj)%nobs
327               IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
328               llcycle = .TRUE.
329               DO jvar = 1, kvars
330                  IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
331                     llcycle = .FALSE.
332                     EXIT
333                  ENDIF
334               END DO
335               IF ( llcycle ) CYCLE
336               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
337                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
338                  inowin = inowin + 1
339                  zlam(inowin) = inpfiles(jj)%plam(ji)
340                  zphi(inowin) = inpfiles(jj)%pphi(ji)
341               ENDIF
342            END DO
343
344            ! Assume anything other than velocity is on T grid
345            IF ( TRIM( inpfiles(jj)%cname(1) ) == 'UVEL' ) THEN
346               CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,1), iobsj(:,1), &
347                  &                  iproc(:,1), 'U' )
348               CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,2), iobsj(:,2), &
349                  &                  iproc(:,2), 'V' )
350            ELSE
351               CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,1), iobsj(:,1), &
352                  &                  iproc(:,1), 'T' )
353               IF ( kvars > 1 ) THEN
354                  DO jvar = 2, kvars
355                     iobsi(:,jvar) = iobsi(:,1)
356                     iobsj(:,jvar) = iobsj(:,1)
357                     iproc(:,jvar) = iproc(:,1)
358                  END DO
359               ENDIF
360            ENDIF
361
362            inowin = 0
363            DO ji = 1, inpfiles(jj)%nobs
364               IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
365               llcycle = .TRUE.
366               DO jvar = 1, kvars
367                  IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
368                     llcycle = .FALSE.
369                     EXIT
370                  ENDIF
371               END DO
372               IF ( llcycle ) CYCLE
373               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
374                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
375                  inowin = inowin + 1
376                  DO jvar = 1, kvars
377                     inpfiles(jj)%iproc(ji,jvar) = iproc(inowin,jvar)
378                     inpfiles(jj)%iobsi(ji,jvar) = iobsi(inowin,jvar)
379                     inpfiles(jj)%iobsj(ji,jvar) = iobsj(inowin,jvar)
380                  END DO
381                  IF ( kvars > 1 ) THEN
382                     DO jvar = 2, kvars
383                        IF ( inpfiles(jj)%iproc(ji,jvar) /= &
384                           & inpfiles(jj)%iproc(ji,1) ) THEN
385                           CALL ctl_stop( 'Error in obs_read_prof:', &
386                              & 'observation on different processors for different vars')
387                        ENDIF
388                     END DO
389                  ENDIF
390               ENDIF
391            END DO
392            DEALLOCATE( zlam, zphi, iobsi, iobsj, iproc )
393
394            DO ji = 1, inpfiles(jj)%nobs
395               IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
396               llcycle = .TRUE.
397               DO jvar = 1, kvars
398                  IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
399                     llcycle = .FALSE.
400                     EXIT
401                  ENDIF
402               END DO
403               IF ( llcycle ) CYCLE
404               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
405                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
406                  IF ( nproc == 0 ) THEN
407                     IF ( inpfiles(jj)%iproc(ji,1) >  nproc ) CYCLE
408                  ELSE
409                     IF ( inpfiles(jj)%iproc(ji,1) /= nproc ) CYCLE
410                  ENDIF
411                  llvalprof = .FALSE.
412                  DO jvar = 1, kvars
413                     IF ( ldvar(jvar) ) THEN
414                        DO ij = 1,inpfiles(jj)%nlev
415                           IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
416                              & CYCLE
417                           IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. &
418                              & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN
419                              ivart0(jvar) = ivart0(jvar) + 1
420                           ENDIF
421                        END DO
422                     ENDIF
423                  END DO
424                  DO ij = 1,inpfiles(jj)%nlev
425                     IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
426                        & CYCLE
427                     DO jvar = 1, kvars
428                        IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. &
429                           &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
430                           &    ldvar(jvar) ) ) THEN
431                           ip3dt = ip3dt + 1
432                           llvalprof = .TRUE.
433                           EXIT
434                        ENDIF
435                     END DO
436                  END DO
437
438                  IF ( llvalprof ) iprof = iprof + 1
439
440               ENDIF
441            END DO
442
443         ENDIF
444
445      END DO prof_files
446
447      !-----------------------------------------------------------------------
448      ! Get the time ordered indices of the input data
449      !-----------------------------------------------------------------------
450
451      !---------------------------------------------------------------------
452      !  Loop over input data files to count total number of profiles
453      !---------------------------------------------------------------------
454      iproftot = 0
455      DO jj = 1, inobf
456         DO ji = 1, inpfiles(jj)%nobs
457            IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
458            llcycle = .TRUE.
459            DO jvar = 1, kvars
460               IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
461                  llcycle = .FALSE.
462                  EXIT
463               ENDIF
464            END DO
465            IF ( llcycle ) CYCLE
466            IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
467               & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
468               iproftot = iproftot + 1
469            ENDIF
470         END DO
471      END DO
472
473      ALLOCATE( iindx(iproftot), ifileidx(iproftot), &
474         &      iprofidx(iproftot), zdat(iproftot) )
475      jk = 0
476      DO jj = 1, inobf
477         DO ji = 1, inpfiles(jj)%nobs
478            IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
479            llcycle = .TRUE.
480            DO jvar = 1, kvars
481               IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
482                  llcycle = .FALSE.
483                  EXIT
484               ENDIF
485            END DO
486            IF ( llcycle ) CYCLE
487            IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
488               & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
489               jk = jk + 1
490               ifileidx(jk) = jj
491               iprofidx(jk) = ji
492               zdat(jk)     = inpfiles(jj)%ptim(ji)
493            ENDIF
494         END DO
495      END DO
496      CALL sort_dp_indx( iproftot, &
497         &               zdat,     &
498         &               iindx   )
499
500      iv3dt(:) = -1
501      IF (ldsatt) THEN
502         iv3dt(:) = ip3dt
503      ELSE
504         iv3dt(:) = ivart0(:)
505      ENDIF
506      CALL obs_prof_alloc( profdata, kvars, kextr, iprof, iv3dt, &
507         &                 kstp, jpi, jpj, jpk )
508
509      ! * Read obs/positions, QC, all variable and assign to profdata
510
511      profdata%nprof     = 0
512      profdata%nvprot(:) = 0
513      profdata%cvars(:)  = clvarsin(:)
514      iprof = 0
515
516      ip3dt = 0
517      ivart(:) = 0
518      itypvar   (:,:) = 0
519      itypvarmpp(:,:) = 0
520
521      ioserrcount = 0
522      DO jk = 1, iproftot
523
524         jj = ifileidx(iindx(jk))
525         ji = iprofidx(iindx(jk))
526
527         IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
528         llcycle = .TRUE.
529         DO jvar = 1, kvars
530            IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
531               llcycle = .FALSE.
532               EXIT
533            ENDIF
534         END DO
535         IF ( llcycle ) CYCLE
536
537         IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND.  &
538            & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
539
540            IF ( nproc == 0 ) THEN
541               IF ( inpfiles(jj)%iproc(ji,1) >  nproc ) CYCLE
542            ELSE
543               IF ( inpfiles(jj)%iproc(ji,1) /= nproc ) CYCLE
544            ENDIF
545
546            llvalprof = .FALSE.
547
548            IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
549
550            IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
551            llcycle = .TRUE.
552            DO jvar = 1, kvars
553               IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
554                  llcycle = .FALSE.
555                  EXIT
556               ENDIF
557            END DO
558            IF ( llcycle ) CYCLE
559
560            loop_prof : DO ij = 1, inpfiles(jj)%nlev
561
562               IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
563                  & CYCLE
564
565               DO jvar = 1, kvars
566                  IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. &
567                     & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN
568
569                     llvalprof = .TRUE. 
570                     EXIT loop_prof
571
572                  ENDIF
573               END DO
574
575            END DO loop_prof
576
577            ! Set profile information
578
579            IF ( llvalprof ) THEN
580
581               iprof = iprof + 1
582
583               CALL jul2greg( isec,                   &
584                  &           imin,                   &
585                  &           ihou,                   &
586                  &           iday,                   &
587                  &           imon,                   &
588                  &           iyea,                   &
589                  &           inpfiles(jj)%ptim(ji), &
590                  &           irefdate(jj) )
591
592
593               ! Profile time coordinates
594               profdata%nyea(iprof) = iyea
595               profdata%nmon(iprof) = imon
596               profdata%nday(iprof) = iday
597               profdata%nhou(iprof) = ihou
598               profdata%nmin(iprof) = imin
599
600               ! Profile space coordinates
601               profdata%rlam(iprof) = inpfiles(jj)%plam(ji)
602               profdata%rphi(iprof) = inpfiles(jj)%pphi(ji)
603
604               ! Coordinate search parameters
605               DO jvar = 1, kvars
606                  profdata%mi  (iprof,jvar) = inpfiles(jj)%iobsi(ji,jvar)
607                  profdata%mj  (iprof,jvar) = inpfiles(jj)%iobsj(ji,jvar)
608               END DO
609
610               ! Profile WMO number
611               profdata%cwmo(iprof) = inpfiles(jj)%cdwmo(ji)
612
613               ! Instrument type
614               READ( inpfiles(jj)%cdtyp(ji), '(I4)', IOSTAT = ios, ERR = 901 ) itype
615901            IF ( ios /= 0 ) THEN
616                  IF (ioserrcount == 0) CALL ctl_warn ( 'Problem converting an instrument type to integer. Setting type to zero' )
617                  ioserrcount = ioserrcount + 1
618                  itype = 0
619               ENDIF
620
621               profdata%ntyp(iprof) = itype
622
623               ! QC stuff
624
625               profdata%nqc(iprof)     = inpfiles(jj)%ioqc(ji)
626               profdata%nqcf(:,iprof)  = inpfiles(jj)%ioqcf(:,ji)
627               profdata%ipqc(iprof)    = inpfiles(jj)%ipqc(ji)
628               profdata%ipqcf(:,iprof) = inpfiles(jj)%ipqcf(:,ji)
629               profdata%itqc(iprof)    = inpfiles(jj)%itqc(ji)
630               profdata%itqcf(:,iprof) = inpfiles(jj)%itqcf(:,ji)
631               profdata%ivqc(iprof,:)  = inpfiles(jj)%ivqc(ji,:)
632               profdata%ivqcf(:,iprof,:) = inpfiles(jj)%ivqcf(:,ji,:)
633
634               ! Bookkeeping data to match profiles
635               profdata%npidx(iprof) = iprof
636               profdata%npfil(iprof) = iindx(jk)
637
638               ! Observation QC flag (whole profile)
639               profdata%nqc(iprof)  = 0 !TODO
640
641               loop_p : DO ij = 1, inpfiles(jj)%nlev
642
643                  IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
644                     & CYCLE
645
646                  IF (ldsatt) THEN
647
648                     DO jvar = 1, kvars
649                        IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. &
650                           &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
651                           &    ldvar(jvar) ) ) THEN
652                           ip3dt = ip3dt + 1
653                           EXIT
654                        ELSE IF ( jvar == kvars ) THEN
655                           CYCLE loop_p
656                        ENDIF
657                     END DO
658
659                  ENDIF
660
661                  DO jvar = 1, kvars
662                 
663                     IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. &
664                       &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
665                       &    ldvar(jvar) ) .OR. ldsatt ) THEN
666
667                        IF (ldsatt) THEN
668
669                           ivart(jvar) = ip3dt
670
671                        ELSE
672
673                           ivart(jvar) = ivart(jvar) + 1
674
675                        ENDIF
676
677                        ! Depth of jvar observation
678                        profdata%var(jvar)%vdep(ivart(jvar)) = &
679                           &                inpfiles(jj)%pdep(ij,ji)
680
681                        ! Depth of jvar observation QC
682                        profdata%var(jvar)%idqc(ivart(jvar)) = &
683                           &                inpfiles(jj)%idqc(ij,ji)
684
685                        ! Depth of jvar observation QC flags
686                        profdata%var(jvar)%idqcf(:,ivart(jvar)) = &
687                           &                inpfiles(jj)%idqcf(:,ij,ji)
688
689                        ! Profile index
690                        profdata%var(jvar)%nvpidx(ivart(jvar)) = iprof
691
692                        ! Vertical index in original profile
693                        profdata%var(jvar)%nvlidx(ivart(jvar)) = ij
694
695                        ! Profile jvar value
696                        IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. &
697                           & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN
698                           profdata%var(jvar)%vobs(ivart(jvar)) = &
699                              &                inpfiles(jj)%pob(ij,ji,jvar)
700                           IF ( ldmod ) THEN
701                              profdata%var(jvar)%vmod(ivart(jvar)) = &
702                                 &                inpfiles(jj)%padd(ij,ji,1,jvar)
703                           ENDIF
704                           ! Count number of profile var1 data as function of type
705                           itypvar( profdata%ntyp(iprof) + 1, jvar ) = &
706                              & itypvar( profdata%ntyp(iprof) + 1, jvar ) + 1
707                        ELSE
708                           profdata%var(jvar)%vobs(ivart(jvar)) = fbrmdi
709                        ENDIF
710
711                        ! Profile jvar qc
712                        profdata%var(jvar)%nvqc(ivart(jvar)) = &
713                           & inpfiles(jj)%ivlqc(ij,ji,jvar)
714
715                        ! Profile jvar qc flags
716                        profdata%var(jvar)%nvqcf(:,ivart(jvar)) = &
717                           & inpfiles(jj)%ivlqcf(:,ij,ji,jvar)
718
719                        ! Profile insitu T value
720                        IF ( TRIM( inpfiles(jj)%cname(jvar) ) == 'POTM' ) THEN
721                           profdata%var(jvar)%vext(ivart(jvar),1) = &
722                              &                inpfiles(jj)%pext(ij,ji,1)
723                        ENDIF
724
725                     ENDIF
726                 
727                  END DO
728
729               END DO loop_p
730
731            ENDIF
732
733         ENDIF
734
735      END DO
736
737      !-----------------------------------------------------------------------
738      ! Sum up over processors
739      !-----------------------------------------------------------------------
740
741      DO jvar = 1, kvars
742         CALL obs_mpp_sum_integer ( ivart0(jvar), ivartmpp(jvar) )
743      END DO
744      CALL obs_mpp_sum_integer ( ip3dt,   ip3dtmpp  )
745
746      DO jvar = 1, kvars
747         CALL obs_mpp_sum_integers( itypvar(:,jvar), itypvarmpp(:,jvar), ntyp1770 + 1 )
748      END DO
749
750      !-----------------------------------------------------------------------
751      ! Output number of observations.
752      !-----------------------------------------------------------------------
753      IF(lwp) THEN
754         WRITE(numout,*) 
755         WRITE(numout,'(A)') ' Profile data'
756         WRITE(numout,'(1X,A)') '------------'
757         WRITE(numout,*) 
758         DO jvar = 1, kvars
759            WRITE(numout,'(1X,A)') 'Profile data, '//TRIM( profdata%cvars(jvar) )
760            WRITE(numout,'(1X,A)') '------------------------'
761            DO ji = 0, ntyp1770
762               IF ( itypvarmpp(ji+1,jvar) > 0 ) THEN
763                  WRITE(numout,'(1X,A3,1X,A48,A3,I8)') ctypshort(ji), &
764                     & cwmonam1770(ji)(1:52),' = ', &
765                     & itypvarmpp(ji+1,jvar)
766               ENDIF
767            END DO
768            WRITE(numout,'(1X,A)') &
769               & '---------------------------------------------------------------'
770            WRITE(numout,'(1X,A55,I8)') &
771               & 'Total profile data for variable '//TRIM( profdata%cvars(jvar) )// &
772               & '             = ', ivartmpp(jvar)
773            WRITE(numout,'(1X,A)') &
774               & '---------------------------------------------------------------'
775            WRITE(numout,*) 
776         END DO
777      ENDIF
778
779      IF (ldsatt) THEN
780         profdata%nvprot(:)    = ip3dt
781         profdata%nvprotmpp(:) = ip3dtmpp
782      ELSE
783         DO jvar = 1, kvars
784            profdata%nvprot(jvar)    = ivart(jvar)
785            profdata%nvprotmpp(jvar) = ivartmpp(jvar)
786         END DO
787      ENDIF
788      profdata%nprof        = iprof
789
790      !-----------------------------------------------------------------------
791      ! Model level search
792      !-----------------------------------------------------------------------
793      DO jvar = 1, kvars
794         IF ( ldvar(jvar) ) THEN
795            CALL obs_level_search( jpk, gdept_1d, &
796               & profdata%nvprot(jvar), profdata%var(jvar)%vdep, &
797               & profdata%var(jvar)%mvk )
798         ENDIF
799      END DO
800
801      !-----------------------------------------------------------------------
802      ! Set model equivalent to 99999
803      !-----------------------------------------------------------------------
804      IF ( .NOT. ldmod ) THEN
805         DO jvar = 1, kvars
806            profdata%var(jvar)%vmod(:) = fbrmdi
807         END DO
808      ENDIF
809      !-----------------------------------------------------------------------
810      ! Deallocate temporary data
811      !-----------------------------------------------------------------------
812      DEALLOCATE( ifileidx, iprofidx, zdat, clvarsin )
813
814      !-----------------------------------------------------------------------
815      ! Deallocate input data
816      !-----------------------------------------------------------------------
817      DO jj = 1, inobf
818         CALL dealloc_obfbdata( inpfiles(jj) )
819      END DO
820      DEALLOCATE( inpfiles )
821
822   END SUBROUTINE obs_rea_prof
823
824END MODULE obs_read_prof
Note: See TracBrowser for help on using the repository browser.