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_sla.F90 in branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_sla.F90 @ 2281

Last change on this file since 2281 was 2281, checked in by smasson, 14 years ago

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