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

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90 @ 2281

Last change on this file since 2281 was 2281, checked in by smasson, 14 years ago

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