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

source: branches/dev_1784_OBS/NEMO/OPA_SRC/OBS/obs_read_prof.F90 @ 2001

Last change on this file since 2001 was 2001, checked in by djlea, 14 years ago

Adding observation operator code
File size: 29.8 KB
Line 
1MODULE obs_read_prof
2   !!======================================================================
3   !!                       ***  MODULE obs_read_prof  ***
4   !! Observation diagnostics: Read the T and S profile observations
5   !!======================================================================
6
7   !!----------------------------------------------------------------------
8   !!   obs_rea_pro_dri : Driver for reading profile obs
9   !!----------------------------------------------------------------------
10
11   !! * Modules used   
12   USE par_kind                 ! Precision variables
13   USE par_oce                  ! Ocean parameters
14   USE in_out_manager           ! I/O manager
15   USE dom_oce                  ! Ocean space and time domain variables
16   USE obs_mpp                  ! MPP support routines for observation diagnostics
17   USE julian                   ! Julian date routines
18   USE obs_utils                ! Observation operator utility functions
19   USE obs_prep                 ! Prepare observation arrays
20   USE obs_grid                 ! Grid search
21   USE obs_sort                 ! Sorting observation arrays
22   USE obs_profiles_def         ! Profile definitions
23   USE obs_conv                 ! Various conversion routines
24   USE obs_types                ! Observation type definitions
25   USE netcdf                   ! NetCDF library
26   USE obs_oper                 ! Observation operators
27   USE 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 :: &
63         & kformat        ! Format of input data
64      !                   ! 1: ENACT
65      !                   ! 2: Coriolis
66      TYPE(obs_prof), INTENT(OUT) :: &
67         & profdata       ! Profile data to be read
68      INTEGER, INTENT(IN) :: &
69         & knumfiles      ! Number of files to read in
70      CHARACTER(LEN=128), INTENT(IN) ::  &
71         & cfilenames(knumfiles) ! File names to read in
72      INTEGER, INTENT(IN) ::  &
73         & kvars,  &      ! Number of variables in profdata
74         & kextr,  &      ! Number of extra fields for each var in profdata
75         & kstp           ! Ocean time-step index
76      LOGICAL, INTENT(IN) :: &
77         & ldt3d, &       ! Observed variables switches
78         & lds3d, &
79         & ldignmis, &    ! Ignore missing files
80         & ldsatt, &      ! Compute salinity at all temperature points
81         & ldavtimset, &  ! Correct time for daily averaged data
82         & ldmod          ! Initialize model from input data
83      REAL(KIND=dp), INTENT(IN) ::  &
84         & ddobsini, &    ! Obs. ini time in YYYYMMDD.HHMMSS
85         & ddobsend       ! Obs. end time in YYYYMMDD.HHMMSS
86      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
87         & kdailyavtypes
88
89      !! * Local declarations
90      CHARACTER(LEN=15), PARAMETER :: cpname='obs_rea_pro_dri'
91      INTEGER :: &
92         & jvar, &
93         & ji,   &
94         & jj,   & 
95         & jk
96      INTEGER :: &
97         & ij,        &
98         & iflag,     &
99         & inobf,     &
100         & i_file_id, &
101         & inowin,    &
102         & iyea,      &
103         & imon,      &
104         & iday,      &
105         & ihou,      &
106         & imin,      &
107         & isec
108      INTEGER, DIMENSION(knumfiles) :: &
109         & irefdate
110      INTEGER, DIMENSION(ntyp1770+1) :: &
111         & itypt,    &
112         & ityptmpp, &
113         & ityps,    &
114         & itypsmpp 
115      INTEGER :: &
116         & it3dtmpp, &
117         & is3dtmpp, &
118         & ip3dtmpp
119      INTEGER, DIMENSION(:), ALLOCATABLE :: &
120         & iobsi,    &
121         & iobsj,    &
122         & iproc,    &
123         & iindx,    &
124         & ifileidx, &
125         & iprofidx
126      INTEGER :: itype
127      INTEGER, DIMENSION(imaxavtypes) :: &
128         & idailyavtypes
129      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
130         & zphi, &
131         & zlam
132      REAL(dp), DIMENSION(:), ALLOCATABLE :: &
133         & zdat
134      LOGICAL :: &
135         & llvalprof
136      TYPE(obfbdata), POINTER, DIMENSION(:) :: &
137         & inpfiles
138      REAL(dp), DIMENSION(knumfiles) :: &
139         & djulini, &
140         & djulend
141      INTEGER :: &
142         & iprof,    &
143         & iproftot, &
144         & it3dt0,    &
145         & is3dt0,    &
146         & it3dt,    &
147         & is3dt,    &
148         & ip3dt
149      INTEGER, DIMENSION(kvars) :: &
150         & iv3dt
151      CHARACTER(len=8) :: &
152         & cl_refdate
153   
154      ! Local initialization
155      iprof = 0
156      it3dt0 = 0
157      is3dt0 = 0
158      ip3dt = 0
159
160      ! Daily average types
161      IF ( PRESENT(kdailyavtypes) ) THEN
162         idailyavtypes(:) = kdailyavtypes(:)
163      ELSE
164         idailyavtypes(:) = -1
165      ENDIF
166
167      !-----------------------------------------------------------------------
168      ! Check data the model part is just with feedback data files
169      !-----------------------------------------------------------------------
170      IF ( ldmod .AND. ( kformat /= 0 ) ) THEN
171         CALL ctl_stop( 'Model can only be read from feedback data' )
172         RETURN
173      ENDIF
174
175      !-----------------------------------------------------------------------
176      ! Count the number of files needed and allocate the obfbdata type
177      !-----------------------------------------------------------------------
178     
179      inobf = knumfiles
180     
181      ALLOCATE( inpfiles(inobf) )
182
183      prof_files : DO jj = 1, inobf
184         
185         !---------------------------------------------------------------------
186         ! Prints
187         !---------------------------------------------------------------------
188         IF(lwp) THEN
189            WRITE(numout,*)
190            WRITE(numout,*) ' obs_rea_pro_dri : Reading from file = ', &
191               & TRIM( TRIM( cfilenames(jj) ) )
192            WRITE(numout,*) ' ~~~~~~~~~~~~~~~'
193            WRITE(numout,*)
194         ENDIF
195
196         !---------------------------------------------------------------------
197         !  Initialization: Open file and get dimensions only
198         !---------------------------------------------------------------------
199         
200         iflag = nf90_open( TRIM( TRIM( cfilenames(jj) ) ), nf90_nowrite, &
201            &                      i_file_id )
202         
203         IF ( iflag /= nf90_noerr ) THEN
204
205            IF ( ldignmis ) THEN
206               inpfiles(jj)%nobs = 0
207               CALL ctl_warn( 'File ' // TRIM( TRIM( cfilenames(jj) ) ) // &
208                  &           ' not found' )
209            ELSE
210               CALL ctl_stop( 'File ' // TRIM( TRIM( cfilenames(jj) ) ) // &
211                  &           ' not found' )
212            ENDIF
213
214         ELSE 
215           
216            !------------------------------------------------------------------
217            !  Close the file since it is opened in read_proffile
218            !------------------------------------------------------------------
219           
220            iflag = nf90_close( i_file_id )
221
222            !------------------------------------------------------------------
223            !  Read the profile file into inpfiles
224            !------------------------------------------------------------------
225            IF ( kformat == 0 ) THEN
226               CALL init_obfbdata( inpfiles(jj) )
227               IF(lwp) THEN
228                  WRITE(numout,*)
229                  WRITE(numout,*)'Reading from feedback file :', &
230                     &           TRIM( cfilenames(jj) )
231               ENDIF
232               CALL read_obfbdata( TRIM( cfilenames(jj) ), inpfiles(jj), &
233                  &                ldgrid = .TRUE. )
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)%ptim(ji) >  djulini(jj) ) .AND. &
298                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
299                  inowin = inowin + 1
300               ENDIF
301            END DO
302            ALLOCATE( zlam(inowin)  )
303            ALLOCATE( zphi(inowin)  )
304            ALLOCATE( iobsi(inowin) )
305            ALLOCATE( iobsj(inowin) )
306            ALLOCATE( iproc(inowin) )
307            inowin = 0
308            DO ji = 1, inpfiles(jj)%nobs
309               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
310                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
311                  inowin = inowin + 1
312                  zlam(inowin) = inpfiles(jj)%plam(ji)
313                  zphi(inowin) = inpfiles(jj)%pphi(ji)
314               ENDIF
315            END DO
316
317            CALL obs_grid_search( inowin, zlam, zphi, iobsi, iobsj, iproc, 'T' )
318
319            inowin = 0
320            DO ji = 1, inpfiles(jj)%nobs
321               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
322                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
323                  inowin = inowin + 1
324                  inpfiles(jj)%iproc(ji,1) = iproc(inowin)
325                  inpfiles(jj)%iobsi(ji,1) = iobsi(inowin)
326                  inpfiles(jj)%iobsj(ji,1) = iobsj(inowin)
327               ENDIF
328            END DO
329            DEALLOCATE( zlam, zphi, iobsi, iobsj, iproc )
330
331            DO ji = 1, inpfiles(jj)%nobs
332               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
333                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
334                  IF ( nproc == 0 ) THEN
335                     IF ( inpfiles(jj)%iproc(ji,1) >  nproc ) CYCLE
336                  ELSE
337                     IF ( inpfiles(jj)%iproc(ji,1) /= nproc ) CYCLE
338                  ENDIF
339                  llvalprof = .FALSE.
340                  IF ( ldt3d ) THEN
341                     loop_t_count : DO ij = 1,inpfiles(jj)%nlev
342                        IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
343                           & CYCLE
344                        IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
345                           & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN
346                           it3dt0 = it3dt0 + 1
347                        ENDIF
348                     END DO loop_t_count
349                  ENDIF
350                  IF ( lds3d ) THEN
351                     loop_s_count : DO ij = 1,inpfiles(jj)%nlev
352                        IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
353                           & CYCLE
354                        IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
355                           & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN
356                           is3dt0 = is3dt0 + 1
357                        ENDIF
358                     END DO loop_s_count
359                  ENDIF
360                  loop_p_count : DO ij = 1,inpfiles(jj)%nlev
361                     IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
362                        & CYCLE
363                     IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
364                        &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. &
365                        &     ldt3d ) .OR. &
366                        & ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
367                        &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. &
368                        &     lds3d ) ) THEN
369                        ip3dt = ip3dt + 1
370                        llvalprof = .TRUE.
371                     ENDIF
372                  END DO loop_p_count
373                  IF ( llvalprof ) iprof = iprof + 1
374               ENDIF
375            END DO
376
377         ENDIF
378
379      END DO prof_files
380     
381      !-----------------------------------------------------------------------
382      ! Get the time ordered indices of the input data
383      !-----------------------------------------------------------------------
384
385      !---------------------------------------------------------------------
386      !  Loop over input data files to count total number of profiles
387      !---------------------------------------------------------------------
388      iproftot = 0
389      DO jj = 1, inobf
390         DO ji = 1, inpfiles(jj)%nobs
391            IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
392               & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
393               iproftot = iproftot + 1
394            ENDIF
395         END DO
396      END DO
397
398      ALLOCATE( iindx(iproftot), ifileidx(iproftot), &
399         &      iprofidx(iproftot), zdat(iproftot) )
400      jk = 0
401      DO jj = 1, inobf
402         DO ji = 1, inpfiles(jj)%nobs
403            IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
404               & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
405               jk = jk + 1
406               ifileidx(jk) = jj
407               iprofidx(jk) = ji
408               zdat(jk)     = inpfiles(jj)%ptim(ji)
409            ENDIF
410         ENDDO
411      ENDDO
412      CALL sort_dp_indx( iproftot, &
413         &               zdat,     &
414         &               iindx   )
415     
416      iv3dt(:) = -1
417      IF (ldsatt) THEN
418         iv3dt(1) = ip3dt
419         iv3dt(2) = ip3dt
420      ELSE
421         iv3dt(1) = it3dt0
422         iv3dt(2) = is3dt0
423      ENDIF
424      CALL obs_prof_alloc( profdata, kvars, kextr, iprof, iv3dt, &
425         &                 kstp, jpi, jpj, jpk )
426     
427      ! * Read obs/positions, QC, all variable and assign to profdata
428
429      profdata%nprof     = 0
430      profdata%nvprot(:) = 0
431
432      iprof = 0
433
434      ip3dt = 0
435      it3dt = 0
436      is3dt = 0
437      itypt   (:) = 0
438      ityptmpp(:) = 0
439     
440      ityps   (:) = 0
441      itypsmpp(:) = 0
442     
443     
444      DO jk = 1, iproftot
445         
446         jj = ifileidx(iindx(jk))
447         ji = iprofidx(iindx(jk))
448         IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND.  &
449            & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
450           
451            IF ( nproc == 0 ) THEN
452               IF ( inpfiles(jj)%iproc(ji,1) >  nproc ) CYCLE
453            ELSE
454               IF ( inpfiles(jj)%iproc(ji,1) /= nproc ) CYCLE
455            ENDIF
456           
457            llvalprof = .FALSE.
458
459            loop_prof : DO ij = 1, inpfiles(jj)%nlev
460               
461               IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
462                  & CYCLE
463               
464               IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
465                  & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN
466                 
467                  llvalprof = .TRUE. 
468                  EXIT loop_prof
469                 
470               ENDIF
471               
472               IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
473                  & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN
474                 
475                  llvalprof = .TRUE. 
476                  EXIT loop_prof
477                 
478               ENDIF
479               
480            END DO loop_prof
481           
482            ! Set profile information
483           
484            IF ( llvalprof ) THEN
485               
486               iprof = iprof + 1
487
488               CALL jul2greg( isec,                   &
489                  &           imin,                   &
490                  &           ihou,                   &
491                  &           iday,                   &
492                  &           imon,                   &
493                  &           iyea,                   &
494                  &           inpfiles(jj)%ptim(ji), &
495                  &           irefdate(jj) )
496
497
498               ! Profile time coordinates
499               profdata%nyea(iprof) = iyea
500               profdata%nmon(iprof) = imon
501               profdata%nday(iprof) = iday
502               profdata%nhou(iprof) = ihou
503               profdata%nmin(iprof) = imin
504               
505               ! Profile space coordinates
506               profdata%rlam(iprof) = inpfiles(jj)%plam(ji)
507               profdata%rphi(iprof) = inpfiles(jj)%pphi(ji)
508
509               ! Coordinate search parameters
510               profdata%mi  (iprof,:) = inpfiles(jj)%iobsi(ji,1)
511               profdata%mj  (iprof,:) = inpfiles(jj)%iobsj(ji,1)
512               
513               ! Profile WMO number
514               profdata%cwmo(iprof) = inpfiles(jj)%cdwmo(ji)
515               
516               ! Instrument type
517               READ( inpfiles(jj)%cdtyp(ji), '(I4)' ) itype
518               profdata%ntyp(iprof) = itype
519               
520               ! QC stuff
521
522               profdata%nqc(iprof)     = inpfiles(jj)%ioqc(ji)
523               profdata%nqcf(:,iprof)  = inpfiles(jj)%ioqcf(:,ji)
524               profdata%ipqc(iprof)    = inpfiles(jj)%ipqc(ji)
525               profdata%ipqcf(:,iprof) = inpfiles(jj)%ipqcf(:,ji)
526               profdata%itqc(iprof)    = inpfiles(jj)%itqc(ji)
527               profdata%itqcf(:,iprof) = inpfiles(jj)%itqcf(:,ji)
528               profdata%ivqc(iprof,:)  = inpfiles(jj)%ivqc(ji,:)
529               profdata%ivqcf(:,iprof,:) = inpfiles(jj)%ivqcf(:,ji,:)
530
531               ! Bookkeeping data to match profiles
532               profdata%npidx(iprof) = iprof
533               profdata%npfil(iprof) = iindx(jk)
534
535               ! Observation QC flag (whole profile)
536               profdata%nqc(iprof)  = 0 !TODO
537
538               loop_p : DO ij = 1, inpfiles(jj)%nlev           
539                 
540                  IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
541                     & CYCLE
542
543                  IF (ldsatt) THEN
544
545                     IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
546                        &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. &
547                        &     ldt3d ) .OR. &
548                        & ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
549                        &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. &
550                        &     lds3d ) ) THEN
551                        ip3dt = ip3dt + 1
552                     ELSE
553                        CYCLE
554                     ENDIF
555                     
556                  ENDIF
557
558                  IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
559                     &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. &
560                     &       ldt3d ) .OR. ldsatt ) THEN
561                     
562                     IF (ldsatt) THEN
563
564                        it3dt = ip3dt
565
566                     ELSE
567
568                        it3dt = it3dt + 1
569                       
570                     ENDIF
571
572                     ! Depth of T observation
573                     profdata%var(1)%vdep(it3dt) = &
574                        &                inpfiles(jj)%pdep(ij,ji)
575                     
576                     ! Depth of T observation QC
577                     profdata%var(1)%idqc(it3dt) = &
578                        &                inpfiles(jj)%idqc(ij,ji)
579                     
580                     ! Depth of T observation QC flags
581                     profdata%var(1)%idqcf(:,it3dt) = &
582                        &                inpfiles(jj)%idqcf(:,ij,ji)
583                     
584                     ! Profile index
585                     profdata%var(1)%nvpidx(it3dt) = iprof
586                     
587                     ! Vertical index in original profile
588                     profdata%var(1)%nvlidx(it3dt) = ij
589
590                     ! Profile potential T value
591                     IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
592                        & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN
593                        profdata%var(1)%vobs(it3dt) = &
594                           &                inpfiles(jj)%pob(ij,ji,1)
595                        IF ( ldmod ) THEN
596                           profdata%var(1)%vmod(it3dt) = &
597                              &                inpfiles(jj)%padd(ij,ji,1,1)
598                        ENDIF
599                        ! Count number of profile T data as function of type
600                        itypt( profdata%ntyp(iprof) + 1 ) = &
601                           & itypt( profdata%ntyp(iprof) + 1 ) + 1
602                     ELSE
603                        profdata%var(1)%vobs(it3dt) = fbrmdi
604                     ENDIF
605
606                     ! Profile T qc
607                     profdata%var(1)%nvqc(it3dt) = &
608                        & inpfiles(jj)%ivlqc(ij,ji,1)
609
610                     ! Profile T qc flags
611                     profdata%var(1)%nvqcf(:,it3dt) = &
612                        & inpfiles(jj)%ivlqcf(:,ij,ji,1)
613
614                     ! Profile insitu T value
615                     profdata%var(1)%vext(it3dt,1) = &
616                        &                inpfiles(jj)%pext(ij,ji,1)
617                     
618                  ENDIF
619                 
620                  IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
621                     &   ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. &
622                     &   lds3d ) .OR. ldsatt ) THEN
623                     
624                     IF (ldsatt) THEN
625
626                        is3dt = ip3dt
627
628                     ELSE
629
630                        is3dt = is3dt + 1
631                       
632                     ENDIF
633
634                     ! Depth of S observation
635                     profdata%var(2)%vdep(is3dt) = &
636                        &                inpfiles(jj)%pdep(ij,ji)
637                     
638                     ! Depth of S observation QC
639                     profdata%var(2)%idqc(is3dt) = &
640                        &                inpfiles(jj)%idqc(ij,ji)
641                     
642                     ! Depth of S observation QC flags
643                     profdata%var(2)%idqcf(:,is3dt) = &
644                        &                inpfiles(jj)%idqcf(:,ij,ji)
645                     
646                     ! Profile index
647                     profdata%var(2)%nvpidx(is3dt) = iprof
648                     
649                     ! Vertical index in original profile
650                     profdata%var(2)%nvlidx(is3dt) = ij
651
652                     ! Profile S value
653                     IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
654                        & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN
655                        profdata%var(2)%vobs(is3dt) = &
656                           &                inpfiles(jj)%pob(ij,ji,2)
657                        IF ( ldmod ) THEN
658                           profdata%var(2)%vmod(is3dt) = &
659                              &                inpfiles(jj)%padd(ij,ji,1,2)
660                        ENDIF
661                        ! Count number of profile S data as function of type
662                        ityps( profdata%ntyp(iprof) + 1 ) = &
663                           & ityps( profdata%ntyp(iprof) + 1 ) + 1
664                     ELSE
665                        profdata%var(2)%vobs(is3dt) = fbrmdi
666                     ENDIF
667                     
668                     ! Profile S qc
669                     profdata%var(2)%nvqc(is3dt) = &
670                        & inpfiles(jj)%ivlqc(ij,ji,2)
671
672                     ! Profile S qc flags
673                     profdata%var(2)%nvqcf(:,is3dt) = &
674                        & inpfiles(jj)%ivlqcf(:,ij,ji,2)
675
676                  ENDIF
677           
678               END DO loop_p
679
680            ENDIF
681
682         ENDIF
683
684      END DO
685
686      !-----------------------------------------------------------------------
687      ! Sum up over processors
688      !-----------------------------------------------------------------------
689     
690      CALL obs_mpp_sum_integer ( it3dt0, it3dtmpp )
691      CALL obs_mpp_sum_integer ( is3dt0, is3dtmpp )
692      CALL obs_mpp_sum_integer ( ip3dt, ip3dtmpp )
693     
694      CALL obs_mpp_sum_integers( itypt, ityptmpp, ntyp1770 + 1 )
695      CALL obs_mpp_sum_integers( ityps, itypsmpp, ntyp1770 + 1 )
696     
697      !-----------------------------------------------------------------------
698      ! Output number of observations.
699      !-----------------------------------------------------------------------
700      IF(lwp) THEN
701         WRITE(numout,*) 
702         WRITE(numout,'(1X,A)') 'Profile data'
703         WRITE(numout,'(1X,A)') '------------'
704         WRITE(numout,*) 
705         WRITE(numout,'(1X,A)') 'Profile T data'
706         WRITE(numout,'(1X,A)') '--------------'
707         DO ji = 0, ntyp1770
708            IF ( ityptmpp(ji+1) > 0 ) THEN
709               WRITE(numout,'(1X,A3,1X,A48,A3,I8)') ctypshort(ji), &
710                  & cwmonam1770(ji)(1:52),' = ', &
711                  & ityptmpp(ji+1)
712            ENDIF
713         ENDDO
714         WRITE(numout,'(1X,A)') &
715            & '---------------------------------------------------------------'
716         WRITE(numout,'(1X,A55,I8)') &
717            & 'Total profile T data                                 = ',&
718            & it3dtmpp
719         WRITE(numout,'(1X,A)') &
720            & '---------------------------------------------------------------'
721         WRITE(numout,*) 
722         WRITE(numout,'(1X,A)') 'Profile S data'
723         WRITE(numout,'(1X,A)') '--------------'
724         DO ji = 0, ntyp1770
725            IF ( itypsmpp(ji+1) > 0 ) THEN
726               WRITE(numout,'(1X,A3,1X,A48,A3,I8)') ctypshort(ji), &
727                  & cwmonam1770(ji)(1:52),' = ', &
728                  & itypsmpp(ji+1)
729            ENDIF
730         ENDDO
731         WRITE(numout,'(1X,A)') &
732            & '---------------------------------------------------------------'
733         WRITE(numout,'(1X,A55,I8)') &
734            & 'Total profile S data                                 = ',&
735            & is3dtmpp
736         WRITE(numout,'(1X,A)') &
737            & '---------------------------------------------------------------'
738         WRITE(numout,*) 
739      ENDIF
740     
741      IF (ldsatt) THEN
742         profdata%nvprot(1)    = ip3dt
743         profdata%nvprot(2)    = ip3dt
744         profdata%nvprotmpp(1) = ip3dtmpp
745         profdata%nvprotmpp(2) = ip3dtmpp
746      ELSE
747         profdata%nvprot(1)    = it3dt
748         profdata%nvprot(2)    = is3dt
749         profdata%nvprotmpp(1) = it3dtmpp
750         profdata%nvprotmpp(2) = is3dtmpp
751      ENDIF
752      profdata%nprof        = iprof
753
754      !-----------------------------------------------------------------------
755      ! Model level search
756      !-----------------------------------------------------------------------
757      IF ( ldt3d ) THEN
758         CALL obs_level_search( jpk, gdept_0, &
759            & profdata%nvprot(1), profdata%var(1)%vdep, &
760            & profdata%var(1)%mvk )
761      ENDIF
762      IF ( lds3d ) THEN
763         CALL obs_level_search( jpk, gdept_0, &
764            & profdata%nvprot(2), profdata%var(2)%vdep, &
765            & profdata%var(2)%mvk )
766      ENDIF
767     
768      !-----------------------------------------------------------------------
769      ! Set model equivalent to 99999
770      !-----------------------------------------------------------------------
771      IF ( .NOT. ldmod ) THEN
772         DO jvar = 1, kvars
773            profdata%var(jvar)%vmod(:) = fbrmdi
774         END DO
775      ENDIF
776      !-----------------------------------------------------------------------
777      ! Deallocate temporary data
778      !-----------------------------------------------------------------------
779      DEALLOCATE( ifileidx, iprofidx, zdat )
780
781      !-----------------------------------------------------------------------
782      ! Deallocate input data
783      !-----------------------------------------------------------------------
784      DO jj = 1, inobf
785         CALL dealloc_obfbdata( inpfiles(jj) )
786      END DO
787      DEALLOCATE( inpfiles )
788
789   END SUBROUTINE obs_rea_pro_dri
790
791END MODULE obs_read_prof
Note: See TracBrowser for help on using the repository browser.