source: branches/UKMO/r6232_tracer_advection/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_sla.F90

Last change on this file was 9295, checked in by jcastill, 3 years ago

Remove svn keywords

File size: 20.9 KB
Line 
1MODULE obs_read_sla
2   !!======================================================================
3   !!                       ***  MODULE obs_read_sla  ***
4   !! Observation diagnostics: Read the along track SLA data from
5   !!                          the AVISO/CLS database or any SLA data
6   !!                          from feedback files
7   !!======================================================================
8
9   !!----------------------------------------------------------------------
10   !!   obs_rea_sla : Driver for reading SLA data from the AVISO/CLS/feedback
11   !!----------------------------------------------------------------------
12
13   !! * Modules used   
14   USE par_kind                 ! Precision variables
15   USE in_out_manager           ! I/O manager
16   USE dom_oce                  ! Ocean space and time domain variables
17   USE obs_mpp                  ! MPP support routines for observation diagnostics
18   USE julian                   ! Julian date routines
19   USE obs_utils                ! Observation operator utility functions
20   USE obs_grid                 ! Grid search
21   USE obs_sort                 ! Sorting observation arrays
22   USE obs_surf_def             ! Surface observation definitions
23   USE obs_types                ! Observation type definitions
24   USE obs_sla_io               ! I/O for sla files
25   USE netcdf                   ! NetCDF library
26
27   IMPLICIT NONE
28
29   !! * Routine accessibility
30   PRIVATE
31
32   PUBLIC obs_rea_sla  ! Read the SLA observations from the AVISO/SLA database
33
34   !!----------------------------------------------------------------------
35   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
36   !! $Id$
37   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
38   !!----------------------------------------------------------------------
39
40CONTAINS
41
42   SUBROUTINE obs_rea_sla( kformat, &
43      &                    sladata, knumfiles, cfilenames, &
44      &                    kvars, kextr, kstp, ddobsini, ddobsend, &
45      &                    ldignmis, ldmod, ldobstd )
46      !!---------------------------------------------------------------------
47      !!
48      !!                   *** ROUTINE obs_rea_sla ***
49      !!
50      !! ** Purpose : Read from file the SLA data
51      !!
52      !! ** Method  : Depending on kformat either AVISO or
53      !!              feedback data files are read
54      !!
55      !! ** Action  :
56      !!
57      !!
58      !! History : 
59      !!      ! :  2009-01 (K. Mogensen) Initial version based on old versions
60      !!----------------------------------------------------------------------
61      !! * Modules used
62
63      !! * Arguments
64      INTEGER :: kformat    ! Format of input data
65      !                     ! 0: Feedback
66      !                     ! 1: AVISO
67      TYPE(obs_surf), INTENT(INOUT) :: sladata    ! SLA data to be read
68      INTEGER, INTENT(IN) :: knumfiles      ! Number of files to read in
69      CHARACTER(LEN=128), INTENT(IN) ::  &
70         & cfilenames(knumfiles) ! File names to read in
71      INTEGER, INTENT(IN) :: kvars     ! Number of variables in sladata
72      INTEGER, INTENT(IN) :: kextr     ! Number of extra fields for each var in sladata
73      INTEGER, INTENT(IN) :: kstp      ! Ocean time-step index
74      LOGICAL, INTENT(IN) :: ldignmis    ! Ignore missing files
75      LOGICAL, INTENT(IN) :: ldmod       ! Initialize model from input data
76      LOGICAL, INTENT(INOUT), optional :: &
77         & ldobstd        ! Read observation standard deviation from fb. file
78      REAL(KIND=dp), INTENT(IN) :: ddobsini  ! Obs. ini time in YYYYMMDD.HHMMSS
79      REAL(KIND=dp), INTENT(IN) :: ddobsend  ! Obs. end time in YYYYMMDD.HHMMSS
80         
81      !! * Local declarations
82      CHARACTER(LEN=11), PARAMETER :: cpname='obs_rea_sla'
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, DIMENSION(knumfiles) :: &
97         & irefdate
98      INTEGER :: iobsmpp
99      INTEGER, DIMENSION(imaxmissions+1) :: &
100         & ityp, &
101         & itypmpp
102      INTEGER, DIMENSION(:), ALLOCATABLE :: &
103         & iobsi,    &
104         & iobsj,    &
105         & iproc,    &
106         & iindx,    &
107         & ifileidx, &
108         & islaidx
109      INTEGER :: itype 
110      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
111         & zphi, &
112         & zlam
113      real(wp), DIMENSION(:), ALLOCATABLE :: &
114         & zdat
115      LOGICAL :: llvalprof
116      LOGICAL :: llobstd
117      TYPE(obfbdata), POINTER, DIMENSION(:) :: &
118         & inpfiles
119      real(wp), DIMENSION(knumfiles) :: &
120         & djulini, &
121         & djulend
122      INTEGER, DIMENSION(knumfiles) :: &
123         & iobspos, &
124         & iobcpos
125      INTEGER :: iobs
126      INTEGER :: iobstot
127      INTEGER :: ios
128      INTEGER :: ioserrcount
129      CHARACTER(len=8) :: cl_refdate
130   
131      ! Local initialization
132      iobs = 0
133      IF ( PRESENT(ldobstd) ) THEN
134         IF (.NOT.ldmod) THEN
135            llobstd = .false.
136         ELSE
137            llobstd = ldobstd
138         ENDIF
139      ELSE
140         llobstd = .FALSE.
141      ENDIF
142 
143      !-----------------------------------------------------------------------
144      ! Check that the model part is just with feedback data files
145      !-----------------------------------------------------------------------
146      IF ( ldmod .AND. ( kformat /= 0 ) ) THEN
147         CALL ctl_stop( 'Model can only be read from feedback data' )
148         RETURN
149      ENDIF
150
151      !-----------------------------------------------------------------------
152      ! Check that the prescribed obs err is just with feedback data files
153      !-----------------------------------------------------------------------
154      IF ( llobstd .AND. ( kformat /= 0 ) ) THEN
155         CALL ctl_stop( 'Observation error can only be read from feedback files' )
156         RETURN
157      ENDIF
158
159      !-----------------------------------------------------------------------
160      ! Count the number of files needed and allocate the obfbdata type
161      !-----------------------------------------------------------------------
162     
163      inobf = knumfiles
164     
165      ALLOCATE( inpfiles(inobf) )
166
167      sla_files : DO jj = 1, inobf
168         
169         !---------------------------------------------------------------------
170         ! Prints
171         !---------------------------------------------------------------------
172         IF(lwp) THEN
173            WRITE(numout,*)
174            WRITE(numout,*) ' obs_rea_sla : Reading from file = ', &
175               & TRIM( TRIM( cfilenames(jj) ) )
176            WRITE(numout,*) ' ~~~~~~~~~~~'
177            WRITE(numout,*)
178         ENDIF
179
180         !---------------------------------------------------------------------
181         !  Initialization: Open file and get dimensions only
182         !---------------------------------------------------------------------
183         
184         iflag = nf90_open( TRIM( TRIM( cfilenames(jj) ) ), nf90_nowrite, &
185            &                      i_file_id )
186         
187         IF ( iflag /= nf90_noerr ) THEN
188
189            IF ( ldignmis ) THEN
190               inpfiles(jj)%nobs = 0
191               CALL ctl_warn( 'File ' // TRIM( TRIM( cfilenames(jj) ) ) // &
192                  &           ' not found' )
193            ELSE
194               CALL ctl_stop( 'File ' // TRIM( TRIM( cfilenames(jj) ) ) // &
195                  &           ' not found' )
196            ENDIF
197
198         ELSE 
199           
200            !------------------------------------------------------------------
201            !  Close the file since it is opened in read_proffile
202            !------------------------------------------------------------------
203           
204            iflag = nf90_close( i_file_id )
205
206            !------------------------------------------------------------------
207            !  Read the profile file into inpfiles
208            !------------------------------------------------------------------
209            IF ( kformat == 0 ) THEN
210               CALL init_obfbdata( inpfiles(jj) )
211               IF(lwp) THEN
212                  WRITE(numout,*)
213                  WRITE(numout,*)'Reading from feedback file :', &
214                     &           TRIM( cfilenames(jj) )
215               ENDIF
216               CALL read_obfbdata( TRIM( cfilenames(jj) ), inpfiles(jj), &
217                  &                ldgrid = .TRUE. )
218               IF ( inpfiles(jj)%nvar < 1 ) THEN
219                  CALL ctl_stop( 'Feedback format error' )
220                  RETURN
221               ENDIF
222               IF ( TRIM(inpfiles(jj)%cname(1)) /= 'SLA' ) THEN
223                  CALL ctl_stop( 'Feedback format error' )
224                  RETURN
225               ENDIF
226               IF ( ldmod .AND. ( ( inpfiles(jj)%nadd == 0 ) .OR.&
227                  &               ( inpfiles(jj)%next == 0 ) ) ) THEN
228                  CALL ctl_stop( 'Model not in input data' )
229                  RETURN
230               ENDIF
231            ELSEIF ( kformat == 1) THEN
232               CALL read_avisofile( TRIM( cfilenames(jj) ), inpfiles(jj), &
233               &                    numout, lwp, .TRUE. )
234            ELSE
235               CALL ctl_stop( 'File format unknown' )
236            ENDIF
237
238            !------------------------------------------------------------------
239            !  Change longitude (-180,180)
240            !------------------------------------------------------------------
241
242            DO ji = 1, inpfiles(jj)%nobs 
243
244               IF ( inpfiles(jj)%plam(ji) < -180. ) &
245                  &   inpfiles(jj)%plam(ji) = inpfiles(jj)%plam(ji) + 360.
246
247               IF ( inpfiles(jj)%plam(ji) >  180. ) &
248                  &   inpfiles(jj)%plam(ji) = inpfiles(jj)%plam(ji) - 360.
249
250            END DO
251
252            !------------------------------------------------------------------
253            !  Calculate the date  (change eventually)
254            !------------------------------------------------------------------
255            cl_refdate=inpfiles(jj)%cdjuldref(1:8)
256            READ(cl_refdate,'(I8)') irefdate(jj)
257           
258            CALL ddatetoymdhms( ddobsini, iyea, imon, iday, ihou, imin, isec )
259            CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulini(jj), &
260               &           krefdate = irefdate(jj) )
261            CALL ddatetoymdhms( ddobsend, iyea, imon, iday, ihou, imin, isec )
262            CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulend(jj), &
263               &           krefdate = irefdate(jj) )
264            IF ( inpfiles(jj)%nobs > 0 ) THEN
265               inpfiles(jj)%iproc = -1
266               inpfiles(jj)%iobsi = -1
267               inpfiles(jj)%iobsj = -1
268            ENDIF
269            inowin = 0
270            DO ji = 1, inpfiles(jj)%nobs
271               IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
272               IF ( inpfiles(jj)%ivqc(ji,1) > 2 ) CYCLE
273               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
274                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
275                  inowin = inowin + 1
276               ENDIF
277            END DO
278            ALLOCATE( zlam(inowin)  )
279            ALLOCATE( zphi(inowin)  )
280            ALLOCATE( iobsi(inowin) )
281            ALLOCATE( iobsj(inowin) )
282            ALLOCATE( iproc(inowin) )
283            inowin = 0
284            DO ji = 1, inpfiles(jj)%nobs
285               IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
286               IF ( inpfiles(jj)%ivqc(ji,1) > 2 ) CYCLE
287               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
288                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
289                  inowin = inowin + 1
290                  zlam(inowin) = inpfiles(jj)%plam(ji)
291                  zphi(inowin) = inpfiles(jj)%pphi(ji)
292               ENDIF
293            END DO
294
295            CALL obs_grid_search( inowin, zlam, zphi, iobsi, iobsj, iproc, 'T' )
296
297            inowin = 0
298            DO ji = 1, inpfiles(jj)%nobs
299               IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
300               IF ( inpfiles(jj)%ivqc(ji,1) > 2 ) CYCLE
301               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
302                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
303                  inowin = inowin + 1
304                  inpfiles(jj)%iproc(ji,1) = iproc(inowin)
305                  inpfiles(jj)%iobsi(ji,1) = iobsi(inowin)
306                  inpfiles(jj)%iobsj(ji,1) = iobsj(inowin)
307               ENDIF
308            END DO
309            DEALLOCATE( zlam, zphi, iobsi, iobsj, iproc )
310
311            DO ji = 1, inpfiles(jj)%nobs
312               IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
313               IF ( inpfiles(jj)%ivqc(ji,1) > 2 ) CYCLE
314               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
315                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
316                  IF ( nproc == 0 ) THEN
317                     IF ( inpfiles(jj)%iproc(ji,1) >  nproc ) CYCLE
318                  ELSE
319                     IF ( inpfiles(jj)%iproc(ji,1) /= nproc ) CYCLE
320                  ENDIF
321                  llvalprof = .FALSE.
322                  IF ( ( inpfiles(jj)%ivlqc(1,ji,1) == 1 ) .OR. &
323                     & ( inpfiles(jj)%ivlqc(1,ji,1) == 2 ) ) THEN
324                     iobs = iobs + 1
325                  ENDIF
326               ENDIF
327            END DO
328
329         ENDIF
330     
331      END DO sla_files
332
333      IF (llobstd) THEN
334         
335         DO jj = 1, inobf
336            iobspos(jj) = -1
337            iobcpos(jj) = -1
338            DO ji = 1,inpfiles(jj)%nadd
339               IF ( TRIM(inpfiles(jj)%caddname(ji)) == 'OSTD' ) THEN
340                  iobspos(jj)=ji
341               ENDIF
342               IF ( TRIM(inpfiles(jj)%caddname(ji)) == 'OCNT' ) THEN
343                  iobcpos(jj)=ji
344               ENDIF
345            END DO
346         END DO
347         llobstd = ( ( MINVAL(iobspos) > 0 ) .AND. ( MINVAL(iobcpos) > 0 ) )
348         IF (llobstd) THEN
349            IF (lwp) WRITE(numout,*)'SLA superobs information present.'
350         ELSE
351            IF (lwp) WRITE(numout,*)'SLA superobs information not present.'
352         ENDIF
353
354      ENDIF
355     
356      !-----------------------------------------------------------------------
357      ! Get the time ordered indices of the input data
358      !-----------------------------------------------------------------------
359
360      !---------------------------------------------------------------------
361      !  Loop over input data files to count total number of profiles
362      !---------------------------------------------------------------------
363      iobstot = 0
364      DO jj = 1, inobf
365         DO ji = 1, inpfiles(jj)%nobs
366            IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
367            IF ( inpfiles(jj)%ivqc(ji,1) > 2 ) CYCLE
368            IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
369               & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
370               iobstot = iobstot + 1
371            ENDIF
372         END DO
373      END DO
374
375      ALLOCATE( iindx(iobstot), ifileidx(iobstot), &
376         &      islaidx(iobstot), zdat(iobstot) )
377      jk = 0
378      DO jj = 1, inobf
379         DO ji = 1, inpfiles(jj)%nobs
380            IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
381            IF ( inpfiles(jj)%ivqc(ji,1) > 2 ) CYCLE
382            IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
383               & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
384               jk = jk + 1
385               ifileidx(jk) = jj
386               islaidx(jk) = ji
387               zdat(jk)     = inpfiles(jj)%ptim(ji)
388            ENDIF
389         END DO
390      END DO
391      CALL sort_dp_indx( iobstot, &
392         &               zdat,     &
393         &               iindx   )
394     
395      CALL obs_surf_alloc( sladata, iobs, kvars, kextr, kstp, jpi, jpj )
396     
397      ! * Read obs/positions, QC, all variable and assign to sladata
398 
399      iobs = 0
400
401      ityp   (:) = 0
402      itypmpp(:) = 0
403
404      ioserrcount = 0
405     
406      DO jk = 1, iobstot
407         
408         jj = ifileidx(iindx(jk))
409         ji = islaidx(iindx(jk))
410
411         IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
412         IF ( inpfiles(jj)%ivqc(ji,1) > 2 ) CYCLE
413
414         IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND.  &
415            & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
416           
417            IF ( nproc == 0 ) THEN
418               IF ( inpfiles(jj)%iproc(ji,1) >  nproc ) CYCLE
419            ELSE
420               IF ( inpfiles(jj)%iproc(ji,1) /= nproc ) CYCLE
421            ENDIF
422           
423            ! Set observation information
424           
425            IF ( ( inpfiles(jj)%ivlqc(1,ji,1) == 1 ) .OR. &
426               & ( inpfiles(jj)%ivlqc(1,ji,1) == 2 ) ) THEN
427
428               iobs = iobs + 1
429
430               CALL jul2greg( isec,                   &
431                  &           imin,                   &
432                  &           ihou,                   &
433                  &           iday,                   &
434                  &           imon,                   &
435                  &           iyea,                   &
436                  &           inpfiles(jj)%ptim(ji), &
437                  &           irefdate(jj) )
438
439
440               ! SLA time coordinates
441               sladata%nyea(iobs) = iyea
442               sladata%nmon(iobs) = imon
443               sladata%nday(iobs) = iday
444               sladata%nhou(iobs) = ihou
445               sladata%nmin(iobs) = imin
446               
447               ! SLA space coordinates
448               sladata%rlam(iobs) = inpfiles(jj)%plam(ji)
449               sladata%rphi(iobs) = inpfiles(jj)%pphi(ji)
450
451               ! Coordinate search parameters
452               sladata%mi  (iobs) = inpfiles(jj)%iobsi(ji,1)
453               sladata%mj  (iobs) = inpfiles(jj)%iobsj(ji,1)
454               
455               ! Instrument type
456               READ( inpfiles(jj)%cdtyp(ji), '(I4)', IOSTAT = ios, ERR = 901 ) itype
457901            IF ( ios /= 0 ) THEN
458                  IF (ioserrcount == 0) CALL ctl_warn ( 'Problem converting an instrument type to integer. Setting type to zero' ) 
459                  ioserrcount = ioserrcount + 1
460                  itype = 0
461               ENDIF
462               sladata%ntyp(iobs) = itype
463               ityp(itype+1) = ityp(itype+1) + 1
464
465               ! Identifier
466               sladata%cwmo(iobs) = inpfiles(jj)%cdwmo(ji)
467
468               ! Bookkeeping data to match observations
469               sladata%nsidx(iobs) = iobs
470               sladata%nsfil(iobs) = iindx(jk)
471
472               ! QC flags
473               sladata%nqc(iobs) = inpfiles(jj)%ivqc(ji,1)
474
475               ! Observed value
476               sladata%robs(iobs,1) = inpfiles(jj)%pob(1,ji,1)
477
478
479               ! Model and MDT is set to fbrmdi unless read from file
480               IF ( ldmod ) THEN
481                  sladata%rmod(iobs,1) = inpfiles(jj)%padd(1,ji,1,1)
482                  sladata%rext(iobs,1) = inpfiles(jj)%padd(1,ji,2,1)
483                  sladata%rext(iobs,2) = inpfiles(jj)%pext(1,ji,1)
484                  IF (llobstd) THEN
485                     sladata%rext(iobs,3) = &
486                        & inpfiles(jj)%padd(1,ji,iobspos(jj),1)
487                     sladata%rext(iobs,4) = &
488                        & inpfiles(jj)%padd(1,ji,iobcpos(jj),1)
489                  ENDIF
490               ELSE
491                  sladata%rmod(iobs,1) = fbrmdi
492                  sladata%rext(iobs,:) = fbrmdi
493               ENDIF
494
495            ENDIF
496         ENDIF
497
498      END DO
499
500      !-----------------------------------------------------------------------
501      ! Sum up over processors
502      !-----------------------------------------------------------------------
503     
504      CALL obs_mpp_sum_integer( iobs, iobsmpp )
505      CALL obs_mpp_sum_integers( ityp, itypmpp, imaxmissions + 1 )
506     
507      !-----------------------------------------------------------------------
508      ! Output number of observations.
509      !-----------------------------------------------------------------------
510      IF (lwp) THEN
511
512         WRITE(numout,*)
513         WRITE(numout,'(1X,A)')'Altimeter satellites'
514         WRITE(numout,'(1X,A)')'--------------------'
515         DO jj = 1,8
516            IF ( itypmpp(jj) > 0 ) THEN
517               WRITE(numout,'(1X,A38,A2,I10)')calttyp(jj-1),'= ',itypmpp(jj)
518            ENDIF
519         END DO
520         WRITE(numout,'(1X,A50)')'--------------------------------------------------'
521         WRITE(numout,'(1X,A40,I10)')'Total                                 = ',iobsmpp
522         WRITE(numout,*)
523
524      ENDIF
525
526      !-----------------------------------------------------------------------
527      ! Deallocate temporary data
528      !-----------------------------------------------------------------------
529      DEALLOCATE( ifileidx, islaidx, zdat )
530
531      !-----------------------------------------------------------------------
532      ! Deallocate input data
533      !-----------------------------------------------------------------------
534      DO jj = 1, inobf
535         CALL dealloc_obfbdata( inpfiles(jj) )
536      END DO
537      DEALLOCATE( inpfiles )
538
539      !-----------------------------------------------------------------------
540      ! Reset ldobstd if the data is present
541      !-----------------------------------------------------------------------
542      IF ( PRESENT(ldobstd) ) THEN
543         ldobstd = llobstd
544      ENDIF
545
546   END SUBROUTINE obs_rea_sla
547
548END MODULE obs_read_sla
Note: See TracBrowser for help on using the repository browser.