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/r8395_coupling_sequence/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: NEMO/branches/UKMO/r8395_coupling_sequence/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90 @ 10763

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

Remove svn keywords properly

File size: 33.0 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      &                     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=6), 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 ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
310               IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. &
311                  & ( 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 ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
328               IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. &
329                  & ( 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 ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
354               IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. &
355                  & ( 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 ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
376               IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. &
377                  & ( 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 ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
391                           & ( 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 ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
401                           & ( 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 ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
410                        &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. &
411                        &     ldvar1 ) .OR. &
412                        & ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
413                        &     ( 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 ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
440            IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. &
441               & ( 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 ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
455            IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. &
456               & ( 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 ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
504         IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. &
505            & ( 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 ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. &
521               & ( inpfiles(jj)%ivqc(ji,2) > 2 ) ) CYCLE
522
523            loop_prof : DO ij = 1, inpfiles(jj)%nlev
524
525               IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
526                  & CYCLE
527
528               IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
529                  & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN
530
531                  llvalprof = .TRUE. 
532                  EXIT loop_prof
533
534               ENDIF
535
536               IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
537                  & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN
538
539                  llvalprof = .TRUE. 
540                  EXIT loop_prof
541
542               ENDIF
543
544            END DO loop_prof
545
546            ! Set profile information
547
548            IF ( llvalprof ) THEN
549
550               iprof = iprof + 1
551
552               CALL jul2greg( isec,                   &
553                  &           imin,                   &
554                  &           ihou,                   &
555                  &           iday,                   &
556                  &           imon,                   &
557                  &           iyea,                   &
558                  &           inpfiles(jj)%ptim(ji), &
559                  &           irefdate(jj) )
560
561
562               ! Profile time coordinates
563               profdata%nyea(iprof) = iyea
564               profdata%nmon(iprof) = imon
565               profdata%nday(iprof) = iday
566               profdata%nhou(iprof) = ihou
567               profdata%nmin(iprof) = imin
568
569               ! Profile space coordinates
570               profdata%rlam(iprof) = inpfiles(jj)%plam(ji)
571               profdata%rphi(iprof) = inpfiles(jj)%pphi(ji)
572
573               ! Coordinate search parameters
574               profdata%mi  (iprof,1) = inpfiles(jj)%iobsi(ji,1)
575               profdata%mj  (iprof,1) = inpfiles(jj)%iobsj(ji,1)
576               profdata%mi  (iprof,2) = inpfiles(jj)%iobsi(ji,2)
577               profdata%mj  (iprof,2) = inpfiles(jj)%iobsj(ji,2)
578
579               ! Profile WMO number
580               profdata%cwmo(iprof) = inpfiles(jj)%cdwmo(ji)
581
582               ! Instrument type
583               READ( inpfiles(jj)%cdtyp(ji), '(I4)', IOSTAT = ios, ERR = 901 ) itype
584901            IF ( ios /= 0 ) THEN
585                  IF (ioserrcount == 0) CALL ctl_warn ( 'Problem converting an instrument type to integer. Setting type to zero' )
586                  ioserrcount = ioserrcount + 1
587                  itype = 0
588               ENDIF
589
590               profdata%ntyp(iprof) = itype
591
592               ! QC stuff
593
594               profdata%nqc(iprof)     = inpfiles(jj)%ioqc(ji)
595               profdata%nqcf(:,iprof)  = inpfiles(jj)%ioqcf(:,ji)
596               profdata%ipqc(iprof)    = inpfiles(jj)%ipqc(ji)
597               profdata%ipqcf(:,iprof) = inpfiles(jj)%ipqcf(:,ji)
598               profdata%itqc(iprof)    = inpfiles(jj)%itqc(ji)
599               profdata%itqcf(:,iprof) = inpfiles(jj)%itqcf(:,ji)
600               profdata%ivqc(iprof,:)  = inpfiles(jj)%ivqc(ji,:)
601               profdata%ivqcf(:,iprof,:) = inpfiles(jj)%ivqcf(:,ji,:)
602
603               ! Bookkeeping data to match profiles
604               profdata%npidx(iprof) = iprof
605               profdata%npfil(iprof) = iindx(jk)
606
607               ! Observation QC flag (whole profile)
608               profdata%nqc(iprof)  = 0 !TODO
609
610               loop_p : DO ij = 1, inpfiles(jj)%nlev
611
612                  IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
613                     & CYCLE
614
615                  IF (ldsatt) THEN
616
617                     IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
618                        &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. &
619                        &     ldvar1 ) .OR. &
620                        & ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
621                        &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. &
622                        &     ldvar2 ) ) THEN
623                        ip3dt = ip3dt + 1
624                     ELSE
625                        CYCLE
626                     ENDIF
627
628                  ENDIF
629
630                  IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
631                     &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. &
632                     &       ldvar1 ) .OR. ldsatt ) THEN
633
634                     IF (ldsatt) THEN
635
636                        ivar1t = ip3dt
637
638                     ELSE
639
640                        ivar1t = ivar1t + 1
641
642                     ENDIF
643
644                     ! Depth of var1 observation
645                     profdata%var(1)%vdep(ivar1t) = &
646                        &                inpfiles(jj)%pdep(ij,ji)
647
648                     ! Depth of var1 observation QC
649                     profdata%var(1)%idqc(ivar1t) = &
650                        &                inpfiles(jj)%idqc(ij,ji)
651
652                     ! Depth of var1 observation QC flags
653                     profdata%var(1)%idqcf(:,ivar1t) = &
654                        &                inpfiles(jj)%idqcf(:,ij,ji)
655
656                     ! Profile index
657                     profdata%var(1)%nvpidx(ivar1t) = iprof
658
659                     ! Vertical index in original profile
660                     profdata%var(1)%nvlidx(ivar1t) = ij
661
662                     ! Profile var1 value
663                     IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
664                        & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN
665                        profdata%var(1)%vobs(ivar1t) = &
666                           &                inpfiles(jj)%pob(ij,ji,1)
667                        IF ( ldmod ) THEN
668                           profdata%var(1)%vmod(ivar1t) = &
669                              &                inpfiles(jj)%padd(ij,ji,1,1)
670                        ENDIF
671                        ! Count number of profile var1 data as function of type
672                        itypvar1( profdata%ntyp(iprof) + 1 ) = &
673                           & itypvar1( profdata%ntyp(iprof) + 1 ) + 1
674                     ELSE
675                        profdata%var(1)%vobs(ivar1t) = fbrmdi
676                     ENDIF
677
678                     ! Profile var1 qc
679                     profdata%var(1)%nvqc(ivar1t) = &
680                        & inpfiles(jj)%ivlqc(ij,ji,1)
681
682                     ! Profile var1 qc flags
683                     profdata%var(1)%nvqcf(:,ivar1t) = &
684                        & inpfiles(jj)%ivlqcf(:,ij,ji,1)
685
686                     ! Profile insitu T value
687                     IF ( TRIM( inpfiles(jj)%cname(1) ) == 'POTM' ) THEN
688                        profdata%var(1)%vext(ivar1t,1) = &
689                           &                inpfiles(jj)%pext(ij,ji,1)
690                     ENDIF
691
692                  ENDIF
693
694                  IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
695                     &   ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. &
696                     &   ldvar2 ) .OR. ldsatt ) THEN
697
698                     IF (ldsatt) THEN
699
700                        ivar2t = ip3dt
701
702                     ELSE
703
704                        ivar2t = ivar2t + 1
705
706                     ENDIF
707
708                     ! Depth of var2 observation
709                     profdata%var(2)%vdep(ivar2t) = &
710                        &                inpfiles(jj)%pdep(ij,ji)
711
712                     ! Depth of var2 observation QC
713                     profdata%var(2)%idqc(ivar2t) = &
714                        &                inpfiles(jj)%idqc(ij,ji)
715
716                     ! Depth of var2 observation QC flags
717                     profdata%var(2)%idqcf(:,ivar2t) = &
718                        &                inpfiles(jj)%idqcf(:,ij,ji)
719
720                     ! Profile index
721                     profdata%var(2)%nvpidx(ivar2t) = iprof
722
723                     ! Vertical index in original profile
724                     profdata%var(2)%nvlidx(ivar2t) = ij
725
726                     ! Profile var2 value
727                     IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
728                        & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN
729                        profdata%var(2)%vobs(ivar2t) = &
730                           &                inpfiles(jj)%pob(ij,ji,2)
731                        IF ( ldmod ) THEN
732                           profdata%var(2)%vmod(ivar2t) = &
733                              &                inpfiles(jj)%padd(ij,ji,1,2)
734                        ENDIF
735                        ! Count number of profile var2 data as function of type
736                        itypvar2( profdata%ntyp(iprof) + 1 ) = &
737                           & itypvar2( profdata%ntyp(iprof) + 1 ) + 1
738                     ELSE
739                        profdata%var(2)%vobs(ivar2t) = fbrmdi
740                     ENDIF
741
742                     ! Profile var2 qc
743                     profdata%var(2)%nvqc(ivar2t) = &
744                        & inpfiles(jj)%ivlqc(ij,ji,2)
745
746                     ! Profile var2 qc flags
747                     profdata%var(2)%nvqcf(:,ivar2t) = &
748                        & inpfiles(jj)%ivlqcf(:,ij,ji,2)
749
750                  ENDIF
751
752               END DO loop_p
753
754            ENDIF
755
756         ENDIF
757
758      END DO
759
760      !-----------------------------------------------------------------------
761      ! Sum up over processors
762      !-----------------------------------------------------------------------
763
764      CALL obs_mpp_sum_integer ( ivar1t0, ivar1tmpp )
765      CALL obs_mpp_sum_integer ( ivar2t0, ivar2tmpp )
766      CALL obs_mpp_sum_integer ( ip3dt,   ip3dtmpp  )
767
768      CALL obs_mpp_sum_integers( itypvar1, itypvar1mpp, ntyp1770 + 1 )
769      CALL obs_mpp_sum_integers( itypvar2, itypvar2mpp, ntyp1770 + 1 )
770
771      !-----------------------------------------------------------------------
772      ! Output number of observations.
773      !-----------------------------------------------------------------------
774      IF(lwp) THEN
775         WRITE(numout,*) 
776         WRITE(numout,'(A)') ' Profile data'
777         WRITE(numout,'(1X,A)') '------------'
778         WRITE(numout,*) 
779         WRITE(numout,'(1X,A)') 'Profile data, '//TRIM( profdata%cvars(1) )
780         WRITE(numout,'(1X,A)') '------------------------'
781         DO ji = 0, ntyp1770
782            IF ( itypvar1mpp(ji+1) > 0 ) THEN
783               WRITE(numout,'(1X,A3,1X,A48,A3,I8)') ctypshort(ji), &
784                  & cwmonam1770(ji)(1:52),' = ', &
785                  & itypvar1mpp(ji+1)
786            ENDIF
787         END DO
788         WRITE(numout,'(1X,A)') &
789            & '---------------------------------------------------------------'
790         WRITE(numout,'(1X,A55,I8)') &
791            & 'Total profile data for variable '//TRIM( profdata%cvars(1) )// &
792            & '             = ', ivar1tmpp
793         WRITE(numout,'(1X,A)') &
794            & '---------------------------------------------------------------'
795         WRITE(numout,*) 
796         WRITE(numout,'(1X,A)') 'Profile data, '//TRIM( profdata%cvars(2) )
797         WRITE(numout,'(1X,A)') '------------------------'
798         DO ji = 0, ntyp1770
799            IF ( itypvar2mpp(ji+1) > 0 ) THEN
800               WRITE(numout,'(1X,A3,1X,A48,A3,I8)') ctypshort(ji), &
801                  & cwmonam1770(ji)(1:52),' = ', &
802                  & itypvar2mpp(ji+1)
803            ENDIF
804         END DO
805         WRITE(numout,'(1X,A)') &
806            & '---------------------------------------------------------------'
807         WRITE(numout,'(1X,A55,I8)') &
808            & 'Total profile data for variable '//TRIM( profdata%cvars(2) )// &
809            & '             = ', ivar2tmpp
810         WRITE(numout,'(1X,A)') &
811            & '---------------------------------------------------------------'
812         WRITE(numout,*) 
813      ENDIF
814
815      IF (ldsatt) THEN
816         profdata%nvprot(1)    = ip3dt
817         profdata%nvprot(2)    = ip3dt
818         profdata%nvprotmpp(1) = ip3dtmpp
819         profdata%nvprotmpp(2) = ip3dtmpp
820      ELSE
821         profdata%nvprot(1)    = ivar1t
822         profdata%nvprot(2)    = ivar2t
823         profdata%nvprotmpp(1) = ivar1tmpp
824         profdata%nvprotmpp(2) = ivar2tmpp
825      ENDIF
826      profdata%nprof        = iprof
827
828      !-----------------------------------------------------------------------
829      ! Model level search
830      !-----------------------------------------------------------------------
831      IF ( ldvar1 ) THEN
832         CALL obs_level_search( jpk, gdept_1d, &
833            & profdata%nvprot(1), profdata%var(1)%vdep, &
834            & profdata%var(1)%mvk )
835      ENDIF
836      IF ( ldvar2 ) THEN
837         CALL obs_level_search( jpk, gdept_1d, &
838            & profdata%nvprot(2), profdata%var(2)%vdep, &
839            & profdata%var(2)%mvk )
840      ENDIF
841
842      !-----------------------------------------------------------------------
843      ! Set model equivalent to 99999
844      !-----------------------------------------------------------------------
845      IF ( .NOT. ldmod ) THEN
846         DO jvar = 1, kvars
847            profdata%var(jvar)%vmod(:) = fbrmdi
848         END DO
849      ENDIF
850      !-----------------------------------------------------------------------
851      ! Deallocate temporary data
852      !-----------------------------------------------------------------------
853      DEALLOCATE( ifileidx, iprofidx, zdat, clvars )
854
855      !-----------------------------------------------------------------------
856      ! Deallocate input data
857      !-----------------------------------------------------------------------
858      DO jj = 1, inobf
859         CALL dealloc_obfbdata( inpfiles(jj) )
860      END DO
861      DEALLOCATE( inpfiles )
862
863   END SUBROUTINE obs_rea_prof
864
865END MODULE obs_read_prof
Note: See TracBrowser for help on using the repository browser.