source: NEMO/trunk/src/OCE/OBS/obs_read_surf.F90

Last change on this file was 14056, checked in by ayoung, 41 hours ago

Adding branch for ticket #2567 to trunk.

  • Property svn:keywords set to Id
File size: 19.2 KB
Line 
1MODULE obs_read_surf
2   !!======================================================================
3   !!                       ***  MODULE obs_read_surf  ***
4   !! Observation diagnostics: Read the surface data from feedback files
5   !!======================================================================
6
7   !!----------------------------------------------------------------------
8   !!   obs_rea_surf : Driver for reading surface data from feedback files
9   !!----------------------------------------------------------------------
10
11   !! * Modules used
12   USE par_kind                 ! Precision variables
13   USE in_out_manager           ! I/O manager
14   USE dom_oce                  ! Ocean space and time domain variables
15   USE obs_mpp                  ! MPP support routines for observation diagnostics
16   USE julian                   ! Julian date routines
17   USE obs_utils                ! Observation operator utility functions
18   USE obs_grid                 ! Grid search
19   USE obs_sort                 ! Sorting observation arrays
20   USE obs_surf_def             ! Surface observation definitions
21   USE obs_types                ! Observation type definitions
22   USE obs_fbm                  ! Feedback routines
23   USE netcdf                   ! NetCDF library
24
25   IMPLICIT NONE
26
27   !! * Routine accessibility
28   PRIVATE
29
30   PUBLIC obs_rea_surf      ! Read the surface observations from the point data
31
32   !!----------------------------------------------------------------------
33   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
34   !! $Id$
35   !! Software governed by the CeCILL license (see ./LICENSE)
36   !!----------------------------------------------------------------------
37
38CONTAINS
39
40   SUBROUTINE obs_rea_surf( surfdata, knumfiles, cdfilenames, &
41      &                     kvars, kextr, kstp, ddobsini, ddobsend, &
42      &                     ldignmis, ldmod, ldnightav, cdvars )
43      !!---------------------------------------------------------------------
44      !!
45      !!                   *** ROUTINE obs_rea_surf ***
46      !!
47      !! ** Purpose : Read from file the surface data
48      !!
49      !! ** Method  : Read in the data from feedback format files and
50      !!              put into the NEMO internal surface data structure
51      !!
52      !! ** Action  :
53      !!
54      !!
55      !! History : 
56      !!      ! :  2009-01 (K. Mogensen) Initial version based on old versions
57      !!      ! :  2015-02 (M. Martin)   Unify the different surface data type reading.
58      !!----------------------------------------------------------------------
59      !! * Modules used
60
61      !! * Arguments
62      TYPE(obs_surf), INTENT(INOUT) :: &
63         & surfdata                     ! Surface data to be read
64      INTEGER, INTENT(IN) :: knumfiles  ! Number of corio format files to read
65      CHARACTER(LEN=128), INTENT(IN) :: &
66         & cdfilenames(knumfiles)       ! File names to read in
67      INTEGER, INTENT(IN) :: kvars      ! Number of variables in surfdata
68      INTEGER, INTENT(IN) :: kextr      ! Number of extra fields for each var
69      INTEGER, INTENT(IN) :: kstp       ! Ocean time-step index
70      LOGICAL, INTENT(IN) :: ldignmis   ! Ignore missing files
71      LOGICAL, INTENT(IN) :: ldmod      ! Initialize model from input data
72      LOGICAL, INTENT(IN) :: ldnightav  ! Observations represent a night-time average
73      REAL(dp), INTENT(IN) :: ddobsini   ! Obs. ini time in YYYYMMDD.HHMMSS
74      REAL(dp), INTENT(IN) :: ddobsend   ! Obs. end time in YYYYMMDD.HHMMSS
75      CHARACTER(len=8), DIMENSION(kvars), INTENT(IN) :: cdvars
76
77      !! * Local declarations
78      CHARACTER(LEN=11), PARAMETER :: cpname='obs_rea_surf'
79      CHARACTER(len=8) :: clrefdate
80      CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvarsin
81      INTEGER :: ji
82      INTEGER :: jj
83      INTEGER :: jk
84      INTEGER :: iflag
85      INTEGER :: inobf
86      INTEGER :: i_file_id
87      INTEGER :: inowin
88      INTEGER :: iyea
89      INTEGER :: imon
90      INTEGER :: iday
91      INTEGER :: ihou
92      INTEGER :: imin
93      INTEGER :: isec
94      INTEGER :: itype
95      INTEGER :: iobsmpp
96      INTEGER :: iobs
97      INTEGER :: iobstot
98      INTEGER :: ios
99      INTEGER :: ioserrcount
100      INTEGER, PARAMETER :: jpsurfmaxtype = 1024
101      INTEGER, DIMENSION(knumfiles) :: irefdate
102      INTEGER, DIMENSION(jpsurfmaxtype+1) :: &
103         & ityp, &
104         & itypmpp
105      INTEGER, DIMENSION(:), ALLOCATABLE :: &
106         & iobsi,    &
107         & iobsj,    &
108         & iproc,    &
109         & iindx,    &
110         & ifileidx, &
111         & isurfidx
112      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
113         & zphi, &
114         & zlam
115      REAL(dp), DIMENSION(:), ALLOCATABLE :: &
116         & zdat
117      REAL(dp), DIMENSION(knumfiles) :: &
118         & djulini, &
119         & djulend
120      LOGICAL :: llvalprof
121      TYPE(obfbdata), POINTER, DIMENSION(:) :: &
122         & inpfiles
123
124      ! Local initialization
125      iobs = 0
126
127      !-----------------------------------------------------------------------
128      ! Count the number of files needed and allocate the obfbdata type
129      !-----------------------------------------------------------------------
130
131      inobf = knumfiles
132
133      ALLOCATE( inpfiles(inobf) )
134
135      surf_files : DO jj = 1, inobf
136
137         !---------------------------------------------------------------------
138         ! Prints
139         !---------------------------------------------------------------------
140         IF(lwp) THEN
141            WRITE(numout,*)
142            WRITE(numout,*) ' obs_rea_surf : Reading from file = ', &
143               & TRIM( TRIM( cdfilenames(jj) ) )
144            WRITE(numout,*) ' ~~~~~~~~~~~'
145            WRITE(numout,*)
146         ENDIF
147
148         !---------------------------------------------------------------------
149         !  Initialization: Open file and get dimensions only
150         !---------------------------------------------------------------------
151
152         iflag = nf90_open( TRIM( TRIM( cdfilenames(jj) ) ), nf90_nowrite, &
153            &                      i_file_id )
154
155         IF ( iflag /= nf90_noerr ) THEN
156
157            IF ( ldignmis ) THEN
158               inpfiles(jj)%nobs = 0
159               CALL ctl_warn( 'File ' // TRIM( TRIM( cdfilenames(jj) ) ) // &
160                  &           ' not found' )
161            ELSE
162               CALL ctl_stop( 'File ' // TRIM( TRIM( cdfilenames(jj) ) ) // &
163                  &           ' not found' )
164            ENDIF
165
166         ELSE 
167
168            !------------------------------------------------------------------
169            !  Close the file since it is opened in read_obfbdata
170            !------------------------------------------------------------------
171
172            iflag = nf90_close( i_file_id )
173
174            !------------------------------------------------------------------
175            !  Read the surface file into inpfiles
176            !------------------------------------------------------------------
177            CALL init_obfbdata( inpfiles(jj) )
178            CALL read_obfbdata( TRIM( cdfilenames(jj) ), inpfiles(jj), &
179               &                ldgrid = .TRUE. )
180
181            IF ( inpfiles(jj)%nvar /= kvars ) THEN
182               CALL ctl_stop( 'Feedback format error: ', &
183                  &           ' unexpected number of vars in feedback file' )
184            ENDIF
185
186            IF ( ldmod .AND. ( inpfiles(jj)%nadd == 0 ) ) THEN
187               CALL ctl_stop( 'Model not in input data' )
188               RETURN
189            ENDIF
190
191            IF ( jj == 1 ) THEN
192               ALLOCATE( clvarsin( inpfiles(jj)%nvar ) )
193               DO ji = 1, inpfiles(jj)%nvar
194                 clvarsin(ji) = inpfiles(jj)%cname(ji)
195                 IF ( clvarsin(ji) /= cdvars(ji) ) THEN
196                    CALL ctl_stop( 'Feedback file variables do not match', &
197                        &           ' expected variable names for this type' )
198                 ENDIF
199               END DO
200            ELSE
201               DO ji = 1, inpfiles(jj)%nvar
202                  IF ( inpfiles(jj)%cname(ji) /= clvarsin(ji) ) THEN
203                     CALL ctl_stop( 'Feedback file variables not consistent', &
204                        &           ' with previous files for this type' )
205                  ENDIF
206               END DO
207            ENDIF
208
209            IF (lwp) WRITE(numout,*)'Observation file contains ',inpfiles(jj)%nobs,' observations'
210
211            !------------------------------------------------------------------
212            !  Change longitude (-180,180)
213            !------------------------------------------------------------------
214
215            DO ji = 1, inpfiles(jj)%nobs
216
217               IF ( inpfiles(jj)%plam(ji) < -180. ) &
218                  &   inpfiles(jj)%plam(ji) = inpfiles(jj)%plam(ji) + 360.
219
220               IF ( inpfiles(jj)%plam(ji) >  180. ) &
221                  &   inpfiles(jj)%plam(ji) = inpfiles(jj)%plam(ji) - 360.
222
223            END DO
224
225            !------------------------------------------------------------------
226            !  Calculate the date  (change eventually)
227            !------------------------------------------------------------------
228            clrefdate=inpfiles(jj)%cdjuldref(1:8)
229            READ(clrefdate,'(I8)') irefdate(jj)
230
231            CALL ddatetoymdhms( ddobsini, iyea, imon, iday, ihou, imin, isec )
232            CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulini(jj), &
233               &           krefdate = irefdate(jj) )
234            CALL ddatetoymdhms( ddobsend, iyea, imon, iday, ihou, imin, isec )
235            CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulend(jj), &
236               &           krefdate = irefdate(jj) )
237
238            IF ( ldnightav ) THEN
239
240               IF ( lwp ) THEN
241                  WRITE(numout,*)'Resetting time of night-time averaged observations', &
242                     &             ' to the end of the day'
243               ENDIF
244
245               DO ji = 1, inpfiles(jj)%nobs
246                  !  for night-time averaged data force the time
247                  !  to be the last time-step of the day, but still within the day.
248                  IF ( inpfiles(jj)%ptim(ji) >= 0. ) THEN
249                     inpfiles(jj)%ptim(ji) = &
250                        & INT(inpfiles(jj)%ptim(ji)) + 0.9999
251                  ELSE
252                     inpfiles(jj)%ptim(ji) = &
253                        & INT(inpfiles(jj)%ptim(ji)) - 0.0001
254                  ENDIF
255               END DO
256            ENDIF
257
258            IF ( inpfiles(jj)%nobs > 0 ) THEN
259               inpfiles(jj)%iproc = -1
260               inpfiles(jj)%iobsi = -1
261               inpfiles(jj)%iobsj = -1
262            ENDIF
263            inowin = 0
264            DO ji = 1, inpfiles(jj)%nobs
265               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
266                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
267                  inowin = inowin + 1
268               ENDIF
269            END DO
270            ALLOCATE( zlam(inowin)  )
271            ALLOCATE( zphi(inowin)  )
272            ALLOCATE( iobsi(inowin) )
273            ALLOCATE( iobsj(inowin) )
274            ALLOCATE( iproc(inowin) )
275            inowin = 0
276            DO ji = 1, inpfiles(jj)%nobs
277               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
278                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
279                  inowin = inowin + 1
280                  zlam(inowin) = inpfiles(jj)%plam(ji)
281                  zphi(inowin) = inpfiles(jj)%pphi(ji)
282               ENDIF
283            END DO
284
285            CALL obs_grid_search( inowin, zlam, zphi, iobsi, iobsj, iproc, 'T' )
286
287            inowin = 0
288            DO ji = 1, inpfiles(jj)%nobs
289               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
290                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
291                  inowin = inowin + 1
292                  inpfiles(jj)%iproc(ji,1) = iproc(inowin)
293                  inpfiles(jj)%iobsi(ji,1) = iobsi(inowin)
294                  inpfiles(jj)%iobsj(ji,1) = iobsj(inowin)
295               ENDIF
296            END DO
297            DEALLOCATE( zlam, zphi, iobsi, iobsj, iproc )
298
299            DO ji = 1, inpfiles(jj)%nobs
300               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
301                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
302                  IF ( nproc == 0 ) THEN
303                     IF ( inpfiles(jj)%iproc(ji,1) >  nproc ) CYCLE
304                  ELSE
305                     IF ( inpfiles(jj)%iproc(ji,1) /= nproc ) CYCLE
306                  ENDIF
307                  llvalprof = .FALSE.
308                  IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(1,ji,1),2) ) THEN
309                     iobs = iobs + 1
310                  ENDIF
311               ENDIF
312            END DO
313
314         ENDIF
315
316      END DO surf_files
317
318      !-----------------------------------------------------------------------
319      ! Get the time ordered indices of the input data
320      !-----------------------------------------------------------------------
321
322      !---------------------------------------------------------------------
323      !  Loop over input data files to count total number of profiles
324      !---------------------------------------------------------------------
325      iobstot = 0
326      DO jj = 1, inobf
327         DO ji = 1, inpfiles(jj)%nobs
328            IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
329               & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
330               iobstot = iobstot + 1
331            ENDIF
332         END DO
333      END DO
334
335      ALLOCATE( iindx(iobstot), ifileidx(iobstot), &
336         &      isurfidx(iobstot), zdat(iobstot) )
337      jk = 0
338      DO jj = 1, inobf
339         DO ji = 1, inpfiles(jj)%nobs
340            IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
341               & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
342               jk = jk + 1
343               ifileidx(jk) = jj
344               isurfidx(jk) = ji
345               zdat(jk)     = inpfiles(jj)%ptim(ji)
346            ENDIF
347         END DO
348      END DO
349      CALL sort_dp_indx( iobstot, &
350         &               zdat,     &
351         &               iindx   )
352
353      CALL obs_surf_alloc( surfdata, iobs, kvars, kextr, kstp, jpi, jpj )
354
355      ! Read obs/positions, QC, all variable and assign to surfdata
356
357      iobs = 0
358
359      surfdata%cvars(:)  = clvarsin(:)
360
361      ityp   (:) = 0
362      itypmpp(:) = 0
363
364      ioserrcount = 0
365
366      DO jk = 1, iobstot
367
368         jj = ifileidx(iindx(jk))
369         ji = isurfidx(iindx(jk))
370         IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND.  &
371            & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
372
373            IF ( nproc == 0 ) THEN
374               IF ( inpfiles(jj)%iproc(ji,1) >  nproc ) CYCLE
375            ELSE
376               IF ( inpfiles(jj)%iproc(ji,1) /= nproc ) CYCLE
377            ENDIF
378
379            ! Set observation information
380
381            IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(1,ji,1),2) ) THEN
382
383               iobs = iobs + 1
384
385               CALL jul2greg( isec,                   &
386                  &           imin,                   &
387                  &           ihou,                   &
388                  &           iday,                   &
389                  &           imon,                   &
390                  &           iyea,                   &
391                  &           inpfiles(jj)%ptim(ji), &
392                  &           irefdate(jj) )
393
394
395               ! Surface time coordinates
396               surfdata%nyea(iobs) = iyea
397               surfdata%nmon(iobs) = imon
398               surfdata%nday(iobs) = iday
399               surfdata%nhou(iobs) = ihou
400               surfdata%nmin(iobs) = imin
401
402               ! Surface space coordinates
403               surfdata%rlam(iobs) = inpfiles(jj)%plam(ji)
404               surfdata%rphi(iobs) = inpfiles(jj)%pphi(ji)
405
406               ! Coordinate search parameters
407               surfdata%mi  (iobs) = inpfiles(jj)%iobsi(ji,1)
408               surfdata%mj  (iobs) = inpfiles(jj)%iobsj(ji,1)
409
410               ! WMO number
411               surfdata%cwmo(iobs) = inpfiles(jj)%cdwmo(ji)
412
413               ! Instrument type
414               READ( inpfiles(jj)%cdtyp(ji), '(I4)', IOSTAT = ios, ERR = 901 ) itype
415901            IF ( ios /= 0 ) THEN
416                  IF (ioserrcount == 0) THEN
417                     CALL ctl_warn ( 'Problem converting an instrument type ', &
418                        &            'to integer. Setting type to zero' )
419                  ENDIF
420                  ioserrcount = ioserrcount + 1
421                  itype = 0
422               ENDIF
423               surfdata%ntyp(iobs) = itype
424               IF ( itype < jpsurfmaxtype + 1 ) THEN
425                  ityp(itype+1) = ityp(itype+1) + 1
426               ELSE
427                  IF(lwp)WRITE(numout,*)'WARNING:Increase jpsurfmaxtype in ',&
428                     &                  cpname
429               ENDIF
430
431               ! Bookkeeping data to match observations
432               surfdata%nsidx(iobs) = iobs
433               surfdata%nsfil(iobs) = iindx(jk)
434
435               ! QC flags
436               surfdata%nqc(iobs) = inpfiles(jj)%ivqc(ji,1)
437
438               ! Observed value
439               surfdata%robs(iobs,1) = inpfiles(jj)%pob(1,ji,1)
440
441
442               ! Model and MDT is set to fbrmdi unless read from file
443               IF ( ldmod ) THEN
444                  surfdata%rmod(iobs,1) = inpfiles(jj)%padd(1,ji,1,1)
445                  IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) THEN
446                     surfdata%rext(iobs,1) = inpfiles(jj)%padd(1,ji,2,1)
447                     surfdata%rext(iobs,2) = inpfiles(jj)%pext(1,ji,1)
448                  ENDIF
449                ELSE
450                  surfdata%rmod(iobs,1) = fbrmdi
451                  IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) surfdata%rext(iobs,:) = fbrmdi
452               ENDIF
453            ENDIF
454         ENDIF
455
456      END DO
457
458      !-----------------------------------------------------------------------
459      ! Sum up over processors
460      !-----------------------------------------------------------------------
461
462      CALL obs_mpp_sum_integer( iobs, iobsmpp )
463      CALL obs_mpp_sum_integers( ityp, itypmpp, jpsurfmaxtype + 1 )
464
465      !-----------------------------------------------------------------------
466      ! Output number of observations.
467      !-----------------------------------------------------------------------
468      IF (lwp) THEN
469
470         WRITE(numout,*)
471         WRITE(numout,'(1X,A)')TRIM( surfdata%cvars(1) )//' data'
472         WRITE(numout,'(1X,A)')'--------------'
473         DO jj = 1,8
474            IF ( itypmpp(jj) > 0 ) THEN
475               WRITE(numout,'(1X,A4,I4,A3,I10)')'Type ', jj,' = ',itypmpp(jj)
476            ENDIF
477         END DO
478         WRITE(numout,'(1X,A)') &
479            & '---------------------------------------------------------------'
480         WRITE(numout,'(1X,A,I8)') &
481            & 'Total data for variable '//TRIM( surfdata%cvars(1) )// &
482            & '           = ', iobsmpp
483         WRITE(numout,'(1X,A)') &
484            & '---------------------------------------------------------------'
485         WRITE(numout,*)
486
487      ENDIF
488
489      !-----------------------------------------------------------------------
490      ! Deallocate temporary data
491      !-----------------------------------------------------------------------
492      DEALLOCATE( ifileidx, isurfidx, zdat, clvarsin )
493
494      !-----------------------------------------------------------------------
495      ! Deallocate input data
496      !-----------------------------------------------------------------------
497      DO jj = 1, inobf
498         IF ( inpfiles(jj)%lalloc ) THEN
499            CALL dealloc_obfbdata( inpfiles(jj) )
500         ENDIF
501      END DO
502      DEALLOCATE( inpfiles )
503
504   END SUBROUTINE obs_rea_surf
505
506END MODULE obs_read_surf
Note: See TracBrowser for help on using the repository browser.