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

source: branches/UKMO/dev_r5518_obs_oper_update_kd490/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_surf.F90 @ 11255

Last change on this file since 11255 was 11255, checked in by dford, 5 years ago

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