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 NEMO/branches/UKMO/NEMO_4.0.4_generic_obs/src/OCE/OBS – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.4_generic_obs/src/OCE/OBS/obs_read_prof.F90 @ 15187

Last change on this file since 15187 was 15187, checked in by dford, 3 years ago

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