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

source: NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/OBS/obs_read_prof.F90 @ 9957

Last change on this file since 9957 was 9957, checked in by acc, 6 years ago

New development branch for the adaptive-implicit vertical advection (05_Gurvan-Vertical_advection)

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