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

source: trunk/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90 @ 2733

Last change on this file since 2733 was 2715, checked in by rblod, 13 years ago

First attempt to put dynamic allocation on the trunk

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