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_surf.F90 in branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

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, 5 years 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.