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

Last change on this file since 12140 was 12140, checked in by kingr, 20 months ago

Added option to remove tides from SLA bkg by taking average over 24h50m.

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