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

source: branches/UKMO/dev_rev5518_OBS_DoNotAssim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90 @ 7841

Last change on this file since 7841 was 7841, checked in by jwhile, 7 years ago

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