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 @ 15228

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

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