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_FOAM_package/src/OCE/OBS – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.4_FOAM_package/src/OCE/OBS/obs_read_prof.F90 @ 15799

Last change on this file since 15799 was 15799, checked in by dford, 2 years ago

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