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

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

Last change on this file since 12035 was 12035, checked in by kingr, 5 years ago

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