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

source: branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90 @ 11546

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

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