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

Last change on this file since 11468 was 11468, checked in by mattmartin, 17 months ago

Merged changes to allow writing of climatological information to feedback files.

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