New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
obs_read_prof.F90 in branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90 @ 11202

Last change on this file since 11202 was 11202, checked in by jcastill, 5 years ago

Copy of branch branches/UKMO/dev_r5518_obs_oper_update@11130 without namelist_ref changes to allow merging with coupling and biogeochemistry branches

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