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

source: branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90 @ 5659

Last change on this file since 5659 was 5659, checked in by mattmartin, 9 years ago

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