source: branches/UKMO/dev_r5518_obs_oper_update_SlaAssim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_surf.F90 @ 12044

Last change on this file since 12044 was 12044, checked in by kingr, 10 months ago

Bugfix to only alter times of SLA when bkg is time-averaged.

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.