source: branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_surf.F90 @ 11546

Last change on this file since 11546 was 11546, checked in by dford, 23 months ago

Add an observation operator for Kd490, and add some error checking so that the observation variable name read in from a feedback file matches a defined name for each observation type. Previously it could read in an incorrect variable without failure. See https://code.metoffice.gov.uk/trac/utils/ticket/247.

File size: 20.5 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/OPA 3.3 , NEMO Consortium (2010)
34   !! $Id$
35   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
36   !!----------------------------------------------------------------------
37
38CONTAINS
39
40   SUBROUTINE obs_rea_surf( surfdata, knumfiles, cdfilenames, &
41      &                     kvars, kextr, kstp, ddobsini, ddobsend, &
42      &                     ldignmis, ldmod, ldnightav, ldclim, 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      LOGICAL, INTENT(IN) :: ldclim     ! Will include climatology at obs points.
74      REAL(dp), INTENT(IN) :: ddobsini   ! Obs. ini time in YYYYMMDD.HHMMSS
75      REAL(dp), INTENT(IN) :: ddobsend   ! Obs. end time in YYYYMMDD.HHMMSS
76      CHARACTER(len=8), DIMENSION(kvars), INTENT(IN) :: cdvars
77
78      !! * Local declarations
79      CHARACTER(LEN=11), PARAMETER :: cpname='obs_rea_surf'
80      CHARACTER(len=8) :: clrefdate
81      CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvarsin
82      INTEGER :: ji
83      INTEGER :: jj
84      INTEGER :: jk
85      INTEGER :: iflag
86      INTEGER :: inobf
87      INTEGER :: i_file_id
88      INTEGER :: inowin
89      INTEGER :: iyea
90      INTEGER :: imon
91      INTEGER :: iday
92      INTEGER :: ihou
93      INTEGER :: imin
94      INTEGER :: isec
95      INTEGER :: itype
96      INTEGER :: iobsmpp
97      INTEGER :: iobs
98      INTEGER :: iobstot
99      INTEGER :: ios
100      INTEGER :: ioserrcount
101      INTEGER :: iextr
102      INTEGER, PARAMETER :: jpsurfmaxtype = 1024
103      INTEGER, DIMENSION(knumfiles) :: irefdate
104      INTEGER, DIMENSION(jpsurfmaxtype+1) :: &
105         & ityp, &
106         & itypmpp
107      INTEGER, DIMENSION(:), ALLOCATABLE :: &
108         & iobsi,    &
109         & iobsj,    &
110         & iproc,    &
111         & iindx,    &
112         & ifileidx, &
113         & isurfidx, &
114         & iadd_std
115      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
116         & zphi, &
117         & zlam
118      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
119         & zdat
120      REAL(wp), DIMENSION(knumfiles) :: &
121         & djulini, &
122         & djulend
123      LOGICAL :: llvalprof
124      TYPE(obfbdata), POINTER, DIMENSION(:) :: &
125         & inpfiles
126
127      ! Local initialization
128      iobs  = 0
129      iextr = kextr
130      !-----------------------------------------------------------------------
131      ! Count the number of files needed and allocate the obfbdata type
132      !-----------------------------------------------------------------------
133
134      inobf = knumfiles
135
136      ALLOCATE( inpfiles(inobf) )
137      ALLOCATE( iadd_std(inobf) )
138
139      surf_files : DO jj = 1, inobf
140
141         !---------------------------------------------------------------------
142         ! Prints
143         !---------------------------------------------------------------------
144         IF(lwp) THEN
145            WRITE(numout,*)
146            WRITE(numout,*) ' obs_rea_surf : Reading from file = ', &
147               & TRIM( TRIM( cdfilenames(jj) ) )
148            WRITE(numout,*) ' ~~~~~~~~~~~'
149            WRITE(numout,*)
150         ENDIF
151
152         !---------------------------------------------------------------------
153         !  Initialization: Open file and get dimensions only
154         !---------------------------------------------------------------------
155
156         iflag = nf90_open( TRIM( TRIM( cdfilenames(jj) ) ), nf90_nowrite, &
157            &                      i_file_id )
158
159         IF ( iflag /= nf90_noerr ) THEN
160
161            IF ( ldignmis ) THEN
162               inpfiles(jj)%nobs = 0
163               CALL ctl_warn( 'File ' // TRIM( TRIM( cdfilenames(jj) ) ) // &
164                  &           ' not found' )
165            ELSE
166               CALL ctl_stop( 'File ' // TRIM( TRIM( cdfilenames(jj) ) ) // &
167                  &           ' not found' )
168            ENDIF
169
170         ELSE 
171
172            !------------------------------------------------------------------
173            !  Close the file since it is opened in read_obfbdata
174            !------------------------------------------------------------------
175
176            iflag = nf90_close( i_file_id )
177
178            !------------------------------------------------------------------
179            !  Read the surface file into inpfiles
180            !------------------------------------------------------------------
181            CALL init_obfbdata( inpfiles(jj) )
182            CALL read_obfbdata( TRIM( cdfilenames(jj) ), inpfiles(jj), &
183               &                ldgrid = .TRUE. )
184
185            IF ( inpfiles(jj)%nvar /= kvars ) THEN
186               CALL ctl_stop( 'Feedback format error: ', &
187                  &           ' unexpected number of vars in feedback file' )
188            ENDIF
189
190            IF ( ldmod .AND. ( inpfiles(jj)%nadd == 0 ) ) THEN
191               CALL ctl_stop( 'Model not in input data' )
192               RETURN
193            ENDIF
194
195            IF ( jj == 1 ) THEN
196               ALLOCATE( clvarsin( inpfiles(jj)%nvar ) )
197               DO ji = 1, inpfiles(jj)%nvar
198                 clvarsin(ji) = inpfiles(jj)%cname(ji)
199                 IF ( clvarsin(ji) /= cdvars(ji) ) THEN
200                    CALL ctl_stop( 'Feedback file variables do not match', &
201                        &           ' expected variable names for this type' )
202                 ENDIF
203               END DO
204            ELSE
205               DO ji = 1, inpfiles(jj)%nvar
206                  IF ( inpfiles(jj)%cname(ji) /= clvarsin(ji) ) THEN
207                     CALL ctl_stop( 'Feedback file variables not consistent', &
208                        &           ' with previous files for this type' )
209                  ENDIF
210               END DO
211            ENDIF
212
213            iadd_std(jj) = -1
214            IF ( inpfiles(jj)%nadd > 0 ) THEN
215               DO ji = 1, inpfiles(jj)%nadd
216                  IF ( TRIM( inpfiles(jj)%caddname(ji) ) == 'STD' ) THEN
217                     iextr = kextr + 1
218                     iadd_std(jj) = ji
219                  ENDIF
220               END DO
221            ENDIF
222
223            IF(lwp) THEN
224               IF ( iadd_std(jj) /= -1 ) THEN
225                  WRITE(numout,*) ' STD variable available in input file so passing it through the obs oper'
226                  WRITE(numout,*)
227               ENDIF
228            ENDIF
229
230            !------------------------------------------------------------------
231            !  Change longitude (-180,180)
232            !------------------------------------------------------------------
233
234            DO ji = 1, inpfiles(jj)%nobs
235
236               IF ( inpfiles(jj)%plam(ji) < -180. ) &
237                  &   inpfiles(jj)%plam(ji) = inpfiles(jj)%plam(ji) + 360.
238
239               IF ( inpfiles(jj)%plam(ji) >  180. ) &
240                  &   inpfiles(jj)%plam(ji) = inpfiles(jj)%plam(ji) - 360.
241
242            END DO
243
244            !------------------------------------------------------------------
245            !  Calculate the date  (change eventually)
246            !------------------------------------------------------------------
247            clrefdate=inpfiles(jj)%cdjuldref(1:8)
248            READ(clrefdate,'(I8)') irefdate(jj)
249
250            CALL ddatetoymdhms( ddobsini, iyea, imon, iday, ihou, imin, isec )
251            CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulini(jj), &
252               &           krefdate = irefdate(jj) )
253            CALL ddatetoymdhms( ddobsend, iyea, imon, iday, ihou, imin, isec )
254            CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulend(jj), &
255               &           krefdate = irefdate(jj) )
256
257            IF ( ldnightav ) THEN
258
259               IF ( lwp ) THEN
260                  WRITE(numout,*)'Resetting time of night-time averaged observations', &
261                     &             ' to the end of the day'
262               ENDIF
263
264               DO ji = 1, inpfiles(jj)%nobs
265                  !  for night-time averaged data force the time
266                  !  to be the last time-step of the day, but still within the day.
267                  IF ( inpfiles(jj)%ptim(ji) >= 0. ) THEN
268                     inpfiles(jj)%ptim(ji) = &
269                        & INT(inpfiles(jj)%ptim(ji)) + 0.9999
270                  ELSE
271                     inpfiles(jj)%ptim(ji) = &
272                        & INT(inpfiles(jj)%ptim(ji)) - 0.0001
273                  ENDIF
274               END DO
275            ENDIF
276
277            IF ( inpfiles(jj)%nobs > 0 ) THEN
278               inpfiles(jj)%iproc = -1
279               inpfiles(jj)%iobsi = -1
280               inpfiles(jj)%iobsj = -1
281            ENDIF
282            inowin = 0
283            DO ji = 1, inpfiles(jj)%nobs
284               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
285                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
286                  inowin = inowin + 1
287               ENDIF
288            END DO
289            ALLOCATE( zlam(inowin)  )
290            ALLOCATE( zphi(inowin)  )
291            ALLOCATE( iobsi(inowin) )
292            ALLOCATE( iobsj(inowin) )
293            ALLOCATE( iproc(inowin) )
294            inowin = 0
295            DO ji = 1, inpfiles(jj)%nobs
296               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
297                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
298                  inowin = inowin + 1
299                  zlam(inowin) = inpfiles(jj)%plam(ji)
300                  zphi(inowin) = inpfiles(jj)%pphi(ji)
301               ENDIF
302            END DO
303
304            CALL obs_grid_search( inowin, zlam, zphi, iobsi, iobsj, iproc, 'T' )
305
306            inowin = 0
307            DO ji = 1, inpfiles(jj)%nobs
308               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
309                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
310                  inowin = inowin + 1
311                  inpfiles(jj)%iproc(ji,1) = iproc(inowin)
312                  inpfiles(jj)%iobsi(ji,1) = iobsi(inowin)
313                  inpfiles(jj)%iobsj(ji,1) = iobsj(inowin)
314               ENDIF
315            END DO
316            DEALLOCATE( zlam, zphi, iobsi, iobsj, iproc )
317
318            DO ji = 1, inpfiles(jj)%nobs
319               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
320                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
321                  IF ( nproc == 0 ) THEN
322                     IF ( inpfiles(jj)%iproc(ji,1) >  nproc ) CYCLE
323                  ELSE
324                     IF ( inpfiles(jj)%iproc(ji,1) /= nproc ) CYCLE
325                  ENDIF
326                  llvalprof = .FALSE.
327                  IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(1,ji,1),2) ) THEN
328                     iobs = iobs + 1
329                  ENDIF
330               ENDIF
331            END DO
332
333         ENDIF
334
335      END DO surf_files
336
337      !-----------------------------------------------------------------------
338      ! Get the time ordered indices of the input data
339      !-----------------------------------------------------------------------
340
341      !---------------------------------------------------------------------
342      !  Loop over input data files to count total number of profiles
343      !---------------------------------------------------------------------
344      iobstot = 0
345      DO jj = 1, inobf
346         DO ji = 1, inpfiles(jj)%nobs
347            IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
348               & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
349               iobstot = iobstot + 1
350            ENDIF
351         END DO
352      END DO
353
354      ALLOCATE( iindx(iobstot), ifileidx(iobstot), &
355         &      isurfidx(iobstot), zdat(iobstot) )
356      jk = 0
357      DO jj = 1, inobf
358         DO ji = 1, inpfiles(jj)%nobs
359            IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
360               & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
361               jk = jk + 1
362               ifileidx(jk) = jj
363               isurfidx(jk) = ji
364               zdat(jk)     = inpfiles(jj)%ptim(ji)
365            ENDIF
366         END DO
367      END DO
368      CALL sort_dp_indx( iobstot, &
369         &               zdat,     &
370         &               iindx   )
371
372      CALL obs_surf_alloc( surfdata, iobs, kvars, iextr, kstp, jpi, jpj, ldclim )
373
374      ! Read obs/positions, QC, all variable and assign to surfdata
375
376      iobs = 0
377      surfdata%cvars(:)  = clvarsin(:)
378      IF ( ldmod .AND. ( TRIM( surfdata%cvars(1) ) == 'SLA' ) ) THEN
379         surfdata%cext(1) = 'SSH'
380         surfdata%cext(2) = 'MDT'
381      ENDIF
382      IF ( iextr > kextr ) surfdata%cext(iextr) = 'STD'
383
384      ityp   (:) = 0
385      itypmpp(:) = 0
386
387      ioserrcount = 0
388
389      DO jk = 1, iobstot
390
391         jj = ifileidx(iindx(jk))
392         ji = isurfidx(iindx(jk))
393         IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND.  &
394            & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
395
396            IF ( nproc == 0 ) THEN
397               IF ( inpfiles(jj)%iproc(ji,1) >  nproc ) CYCLE
398            ELSE
399               IF ( inpfiles(jj)%iproc(ji,1) /= nproc ) CYCLE
400            ENDIF
401
402            ! Set observation information
403
404            IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(1,ji,1),2) ) THEN
405
406               iobs = iobs + 1
407
408               CALL jul2greg( isec,                   &
409                  &           imin,                   &
410                  &           ihou,                   &
411                  &           iday,                   &
412                  &           imon,                   &
413                  &           iyea,                   &
414                  &           inpfiles(jj)%ptim(ji), &
415                  &           irefdate(jj) )
416
417
418               ! Surface time coordinates
419               surfdata%nyea(iobs) = iyea
420               surfdata%nmon(iobs) = imon
421               surfdata%nday(iobs) = iday
422               surfdata%nhou(iobs) = ihou
423               surfdata%nmin(iobs) = imin
424
425               ! Surface space coordinates
426               surfdata%rlam(iobs) = inpfiles(jj)%plam(ji)
427               surfdata%rphi(iobs) = inpfiles(jj)%pphi(ji)
428
429               ! Coordinate search parameters
430               surfdata%mi  (iobs) = inpfiles(jj)%iobsi(ji,1)
431               surfdata%mj  (iobs) = inpfiles(jj)%iobsj(ji,1)
432
433               ! WMO number
434               surfdata%cwmo(iobs) = inpfiles(jj)%cdwmo(ji)
435
436               ! Instrument type
437               READ( inpfiles(jj)%cdtyp(ji), '(I4)', IOSTAT = ios, ERR = 901 ) itype
438901            IF ( ios /= 0 ) THEN
439                  IF (ioserrcount == 0) THEN
440                     CALL ctl_warn ( 'Problem converting an instrument type ', &
441                        &            'to integer. Setting type to zero' )
442                  ENDIF
443                  ioserrcount = ioserrcount + 1
444                  itype = 0
445               ENDIF
446               surfdata%ntyp(iobs) = itype
447               IF ( itype < jpsurfmaxtype + 1 ) THEN
448                  ityp(itype+1) = ityp(itype+1) + 1
449               ELSE
450                  IF(lwp)WRITE(numout,*)'WARNING:Increase jpsurfmaxtype in ',&
451                     &                  cpname
452               ENDIF
453
454               ! Bookkeeping data to match observations
455               surfdata%nsidx(iobs) = iobs
456               surfdata%nsfil(iobs) = iindx(jk)
457
458               ! QC flags
459               surfdata%nqc(iobs) = inpfiles(jj)%ivqc(ji,1)
460
461               ! Observed value
462               surfdata%robs(iobs,1) = inpfiles(jj)%pob(1,ji,1)
463
464               ! Model and MDT is set to fbrmdi unless read from file
465               IF ( ldmod ) THEN
466                  surfdata%rmod(iobs,1) = inpfiles(jj)%padd(1,ji,1,1)
467                  IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) THEN
468                     surfdata%rext(iobs,1) = inpfiles(jj)%padd(1,ji,2,1)
469                     surfdata%rext(iobs,2) = inpfiles(jj)%pext(1,ji,1)
470                  ENDIF
471                ELSE
472                  surfdata%rmod(iobs,1) = fbrmdi
473                  IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) surfdata%rext(iobs,:) = fbrmdi
474               ENDIF
475               
476               ! Initialise climatology if set
477               IF ( surfdata%lclim ) surfdata%rclm(iobs,1) = fbrmdi
478               
479               ! STD (obs error standard deviation) read from file and passed through obs operator
480               IF ( iadd_std(jj) /= -1 ) THEN
481                  surfdata%rext(iobs,iextr) = inpfiles(jj)%padd(1,ji,iadd_std(jj),1)
482               ENDIF
483            ENDIF
484         ENDIF
485
486      END DO
487
488      !-----------------------------------------------------------------------
489      ! Sum up over processors
490      !-----------------------------------------------------------------------
491
492      CALL obs_mpp_sum_integer( iobs, iobsmpp )
493      CALL obs_mpp_sum_integers( ityp, itypmpp, jpsurfmaxtype + 1 )
494
495      !-----------------------------------------------------------------------
496      ! Output number of observations.
497      !-----------------------------------------------------------------------
498      IF (lwp) THEN
499
500         WRITE(numout,*)
501         WRITE(numout,'(1X,A)')TRIM( surfdata%cvars(1) )//' data'
502         WRITE(numout,'(1X,A)')'--------------'
503         DO jj = 1,8
504            IF ( itypmpp(jj) > 0 ) THEN
505               WRITE(numout,'(1X,A4,I4,A3,I10)')'Type ', jj,' = ',itypmpp(jj)
506            ENDIF
507         END DO
508         WRITE(numout,'(1X,A)') &
509            & '---------------------------------------------------------------'
510         WRITE(numout,'(1X,A,I8)') &
511            & 'Total data for variable '//TRIM( surfdata%cvars(1) )// &
512            & '           = ', iobsmpp
513         WRITE(numout,'(1X,A)') &
514            & '---------------------------------------------------------------'
515         WRITE(numout,*)
516
517      ENDIF
518
519      !-----------------------------------------------------------------------
520      ! Deallocate temporary data
521      !-----------------------------------------------------------------------
522      DEALLOCATE( ifileidx, isurfidx, zdat, clvarsin )
523
524      !-----------------------------------------------------------------------
525      ! Deallocate input data
526      !-----------------------------------------------------------------------
527      DO jj = 1, inobf
528         IF ( inpfiles(jj)%lalloc ) THEN
529            CALL dealloc_obfbdata( inpfiles(jj) )
530         ENDIF
531      END DO
532      DEALLOCATE( inpfiles )
533      DEALLOCATE( iadd_std )
534
535   END SUBROUTINE obs_rea_surf
536
537END MODULE obs_read_surf
Note: See TracBrowser for help on using the repository browser.