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.
obssla_io.h90 in NEMO/branches/UKMO/r8395_restart_datestamp/NEMOGCM/TOOLS/OBSTOOLS/src – NEMO

source: NEMO/branches/UKMO/r8395_restart_datestamp/NEMOGCM/TOOLS/OBSTOOLS/src/obssla_io.h90 @ 10758

Last change on this file since 10758 was 10758, checked in by jcastill, 5 years ago

Remove svn keywords

File size: 14.4 KB
Line 
1   !!----------------------------------------------------------------------
2   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
3   !! $Id$
4   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
5   !!----------------------------------------------------------------------
6
7   SUBROUTINE read_avisofile( cdfilename, inpfile, kunit, ldwp, ldgrid )
8      !!---------------------------------------------------------------------
9      !!
10      !!                     ** ROUTINE read_avisofile **
11      !!
12      !! ** Purpose : Read from file the AVISO SLA observations.
13      !!
14      !! ** Method  : The data file is a NetCDF file.
15      !!
16      !! ** Action  :
17      !!
18      !! References : http://www.aviso.oceanobs.com
19      !!
20      !! History :
21      !!          ! 09-01 (K. Mogensen) Original based on old versions
22      !!----------------------------------------------------------------------
23      !! * Arguments
24      CHARACTER(LEN=*) :: cdfilename ! Input filename
25      TYPE(obfbdata)   :: inpfile    ! Output obfbdata structure
26      INTEGER          :: kunit      ! Unit for output
27      LOGICAL          :: ldwp       ! Print info
28      LOGICAL          :: ldgrid     ! Save grid info in data structure
29      !! * Local declarations
30      CHARACTER(LEN=14),PARAMETER :: cpname = 'read_avisofile'
31      INTEGER :: i_file_id     ! netcdf IDS
32      INTEGER :: i_tracks_id
33      INTEGER :: i_cycles_id
34      INTEGER :: i_data_id
35      INTEGER :: i_var_id
36      INTEGER, PARAMETER :: imaxdim = 2    ! Assumed maximum for no. dims. in file
37      INTEGER, DIMENSION(2) :: idims       ! Dimensions in file
38      INTEGER :: iilen         ! Length of netCDF attributes
39      INTEGER :: itype         ! Typeof netCDF attributes
40      REAL(fbdp) :: zsca       ! Scale factor
41      REAL(fbdp) :: zfill      ! Fill value
42      CHARACTER(len=3) ::  cmission      ! Mission global attribute
43      INTEGER :: itracks       ! Maximum number of passes in file
44      INTEGER :: icycles       ! Maximum number of cycles for each pass
45      INTEGER :: idata         ! Number of data per parameter in current file
46      REAL(fbdp) :: zdeltat    ! Time gap getween two measurements in seconds
47      INTEGER, DIMENSION(:), POINTER :: &
48         & iptracks,     &  ! List of passes contained in current file
49         & ipnbpoints,   &  ! Number of points per pass
50         & ipdataindexes    ! Index of data in theoretical profile
51      INTEGER, DIMENSION(:,:), POINTER :: &
52         & ipcycles         ! List of cycles per pass
53      REAL(fbdp), DIMENSION(:), POINTER :: &
54         & zphi,   &        ! Latitudes
55         & zlam             ! Longitudes
56      REAL(fbdp), DIMENSION(:,:), POINTER :: &
57         & zbegindates      ! Date of point with index 0
58      REAL(fbdp) :: zbeginmiss    ! Missing data for dates
59      REAL(fbsp), DIMENSION(:,:), POINTER :: &
60         & zsla             ! SLA data
61      REAL(fbdp), DIMENSION(:), POINTER :: &
62         & zjuld            ! Julian date
63      LOGICAL, DIMENSION(:), POINTER :: &
64         & llskip           ! Skip observation
65      CHARACTER(len=14) :: cdjuldref    ! Julian data reference
66      INTEGER :: imission   ! Mission number converted from Mission global
67                            ! netCDF atttribute.
68      CHARACTER(len=255) :: ctmp
69      INTEGER :: iobs
70      INTEGER :: jl
71      INTEGER :: jm
72      INTEGER :: jj
73      INTEGER :: ji
74      INTEGER :: jk
75      INTEGER :: jobs
76      INTEGER :: jcycle
77
78      ! Open the file
79
80      CALL chkerr( nf90_open( TRIM( cdfilename ), nf90_nowrite, i_file_id ), &
81         &         cpname, __LINE__ )
82
83      ! Get the netCDF dimensions
84     
85      CALL chkerr( nf90_inq_dimid( i_file_id, 'Tracks', i_tracks_id ),  &
86         &         cpname, __LINE__ )
87      CALL chkerr( nf90_inquire_dimension( i_file_id, i_tracks_id, &
88         &                                 len = itracks ),  &
89         &         cpname, __LINE__ )
90     
91      CALL chkerr( nf90_inq_dimid( i_file_id, 'Cycles', i_cycles_id ),  &
92         &         cpname, __LINE__ )
93      CALL chkerr( nf90_inquire_dimension( i_file_id, i_cycles_id, &
94         &                                 len = icycles ),  &
95         &      cpname, __LINE__ )
96     
97      CALL chkerr( nf90_inq_dimid( i_file_id, 'Data', i_data_id ),  &
98         &         cpname, __LINE__ )
99      CALL chkerr( nf90_inquire_dimension( i_file_id, i_data_id, &
100         &                                 len = idata ),  &
101         &         cpname, __LINE__ )
102
103      ! Allocate memory for input data
104
105      ALLOCATE( &
106         & iptracks     ( itracks          ), & 
107         & ipnbpoints   ( itracks          ), & 
108         & ipcycles     ( icycles, itracks ), &
109         & zphi         ( idata            ), &   
110         & zlam         ( idata            ), & 
111         & zbegindates  ( icycles, itracks ), &
112         & ipdataindexes( idata            ), &
113         & zsla         ( icycles, idata   ), &
114         & zjuld        ( idata*icycles    ), &
115         & llskip       ( idata*icycles    )  &
116         & )
117
118      ! Get time gap getween two measurements in seconds
119     
120      CALL chkerr( nf90_inq_varid( i_file_id, 'DeltaT', i_var_id ), &
121         &         cpname, __LINE__ )
122      CALL chkdim( i_file_id, i_var_id, 0, idims, cpname, __LINE__ )
123      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, zdeltat), &
124         &         cpname, __LINE__ )
125      zsca = 1.0
126      IF (nf90_inquire_attribute( i_file_id, i_var_id, 'scale_factor') &
127         &                      == nf90_noerr) THEN
128         CALL chkerr( nf90_get_att( i_file_id, i_var_id, &
129            &                     'scale_factor',zsca), cpname,  __LINE__)
130      ENDIF
131      zdeltat = zsca * zdeltat
132         
133      ! Get List of passes contained in current file     
134     
135      CALL chkerr( nf90_inq_varid( i_file_id, 'Tracks', i_var_id ), &
136         &         cpname, __LINE__ )
137      idims(1) = itracks
138      CALL chkdim( i_file_id, i_var_id, 1, idims, cpname, __LINE__ )
139      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, iptracks), &
140         &         cpname, __LINE__ )
141     
142      ! Get Number of points per pass
143     
144      CALL chkerr( nf90_inq_varid( i_file_id, 'NbPoints', i_var_id ), &
145         &         cpname, __LINE__ )
146      idims(1) = itracks
147      CALL chkdim( i_file_id, i_var_id, 1, idims, cpname, __LINE__ )
148      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, ipnbpoints),&
149         &         cpname, __LINE__ )
150     
151      ! Get list of cycles per pass
152     
153      CALL chkerr( nf90_inq_varid( i_file_id, 'Cycles', i_var_id ), &
154         &         cpname, __LINE__ )
155      idims(1) = icycles
156      idims(2) = itracks
157      CALL chkdim( i_file_id, i_var_id, 2, idims, cpname, __LINE__ )
158      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, ipcycles),   &
159         &         cpname, __LINE__ )
160
161      ! Get longitudes
162
163      CALL chkerr( nf90_inq_varid( i_file_id, 'Longitudes', i_var_id ), &
164         &         cpname, __LINE__ )
165      idims(1) = idata
166      CALL chkdim( i_file_id, i_var_id, 1, idims, cpname, __LINE__ )
167      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, zlam),   &
168         &         cpname, __LINE__ )
169      zsca = 1.0
170      IF (nf90_inquire_attribute( i_file_id, i_var_id, 'scale_factor') &
171         &                         == nf90_noerr) THEN
172         CALL chkerr( nf90_get_att( i_file_id, i_var_id, &
173            &                     'scale_factor',zsca), cpname,  __LINE__)
174      ENDIF
175      zlam(:) = zsca * zlam(:)
176     
177      ! Get latitudes
178     
179      CALL chkerr( nf90_inq_varid( i_file_id, 'Latitudes', i_var_id ), &
180         &         cpname, __LINE__ )
181      idims(1) = idata
182      CALL chkdim( i_file_id, i_var_id, 1, idims, cpname, __LINE__ )
183      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, zphi),   &
184         &         cpname, __LINE__ )
185      zsca = 1.0
186      IF (nf90_inquire_attribute( i_file_id, i_var_id, 'scale_factor') &
187         &                       == nf90_noerr) THEN
188         CALL chkerr( nf90_get_att( i_file_id, i_var_id, &
189            &                     'scale_factor',zsca), cpname,  __LINE__)
190      ENDIF
191      zphi(:) = zsca * zphi(:)
192     
193      ! Get date of point with index 0
194     
195      CALL chkerr( nf90_inq_varid( i_file_id, 'BeginDates', i_var_id ), &
196         &        cpname, __LINE__ )
197      idims(1) = icycles
198      idims(2) = itracks
199      CALL chkdim( i_file_id, i_var_id, 2, idims, cpname, __LINE__ )
200      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, zbegindates), &
201         &         cpname, __LINE__ )
202      CALL chkerr( nf90_inquire_attribute( i_file_id, i_var_id, &
203         &                                 'units', len = iilen,  &
204         &                                  xtype = itype), cpname, __LINE__ )
205      IF (( itype /= NF90_CHAR ) .OR. ( iilen > 255 )) THEN
206         CALL fatal_error('Error decoding BeginDates unit', __LINE__ )
207      ENDIF
208      CALL chkerr( nf90_get_att( i_file_id, i_var_id, 'units', &
209         &                       ctmp ), cpname, __LINE__ )
210      jl=1
211      DO jm = 1, LEN(TRIM(ctmp))
212         IF ((ctmp(jm:jm)>='0').AND.(ctmp(jm:jm)<='9')) THEN
213            cdjuldref(jl:jl) = ctmp(jm:jm)
214            jl=jl+1
215         ENDIF
216         IF (jl>14) EXIT
217      END DO
218      CALL chkerr( nf90_inquire_attribute( i_file_id, i_var_id, '_FillValue', &
219         &                                  xtype = itype), cpname, __LINE__ )
220      IF ( itype /= NF90_DOUBLE ) THEN
221         CALL fatal_error('Error decoding BeginDates missing data', __LINE__ )
222      ENDIF
223      CALL chkerr( nf90_get_att( i_file_id, i_var_id, '_FillValue', &
224         &                       zbeginmiss ), cpname, __LINE__ )
225
226      ! Get indices of data in theoretical profile
227     
228      CALL chkerr( nf90_inq_varid( i_file_id, 'DataIndexes', i_var_id ), &
229         &         cpname, __LINE__ )
230      idims(1) = idata
231      CALL chkdim( i_file_id, i_var_id, 1, idims, cpname, __LINE__ )
232      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, ipdataindexes),   &
233         &       cpname, __LINE__ )
234     
235      ! Get SLA data
236     
237      CALL chkerr( nf90_inq_varid( i_file_id, 'SLA', i_var_id ), &
238         &         cpname, __LINE__ )
239      idims(1) = icycles
240      idims(2) = idata
241      CALL chkdim( i_file_id, i_var_id, 2, idims, cpname, __LINE__ )
242      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, zsla),   &
243         &         cpname, __LINE__ )
244      zsca = 1.0
245      IF (nf90_inquire_attribute( i_file_id, i_var_id, 'scale_factor') &
246         &                        == nf90_noerr) THEN
247         CALL chkerr( nf90_get_att( i_file_id, i_var_id, &
248            &                       'scale_factor',zsca), cpname, __LINE__ )
249      ENDIF
250      zfill = 0.0
251      IF (nf90_inquire_attribute( i_file_id, i_var_id, '_FillValue') &
252         &                      == nf90_noerr) THEN
253         CALL chkerr( nf90_get_att( i_file_id, i_var_id, &
254            &                     '_FillValue',zfill), cpname,  __LINE__ )
255      ENDIF
256      WHERE(zsla(:,:) /=  zfill)
257         zsla(:,:) = zsca * zsla(:,:)
258      ELSEWHERE
259         zsla(:,:) = fbrmdi
260      END WHERE
261     
262      ! Get the global Mission netCDF attribute
263     
264      cmission='   '
265      CALL chkerr( nf90_inquire_attribute( i_file_id, nf90_global, &
266         &                                 'Mission', len = iilen,  &
267         &                                 xtype = itype), cpname, __LINE__ )
268      IF (( itype /= NF90_CHAR ) .OR. ( iilen > 3 )) THEN
269         CALL fatal_error('Error decoding Mission global attribute', __LINE__ )
270      ENDIF
271      CALL chkerr( nf90_get_att( i_file_id, nf90_global, 'Mission', &
272         &                       cmission ), cpname, __LINE__ )
273     
274      ! Convert it to an integer
275      imission = 0
276      DO jm = 1, imaxmissions
277         IF (cmission==cmissions(jm)) THEN
278            imission = jm
279            EXIT
280         ENDIF
281      END DO
282     
283      ! Close the file
284     
285      CALL chkerr( nf90_close( i_file_id ), cpname, __LINE__ )
286
287      ! Compute Julian dates for all observations
288     
289      jm = 0
290      jl = 0
291      DO jj = 1, itracks
292         DO ji = 1, ipnbpoints(jj)
293            jl = jl + 1
294            DO jk = 1, icycles
295               jm = jm + 1
296               IF (zbegindates(jk,jj)==zbeginmiss) THEN
297                  llskip(jm) = .TRUE.
298                  zjuld(jm)  = fbrmdi
299               ELSE
300                  llskip(jm) = .FALSE.
301                  zjuld(jm)  = zbegindates(jk,jj) + &
302                     &         (ipdataindexes(jl) * zdeltat / 86400._wp )
303               ENDIF
304            END DO
305         END DO
306      END DO
307
308      ! Get rid of missing data
309
310      jm = 0
311      DO jobs = 1, idata
312         DO jcycle = 1, icycles
313            jm = jm + 1
314            IF (zsla(jcycle,jobs) == fbrmdi) llskip(jm) = .TRUE.
315         END DO
316      END DO
317     
318      ! Allocate obfbdata
319
320      iobs = COUNT( .NOT.llskip(:) )
321      CALL init_obfbdata( inpfile )
322      CALL alloc_obfbdata( inpfile, 1, iobs, 1, 0, 0, ldgrid )
323      inpfile%cname(1) = 'SLA'
324
325      ! Fill the obfbdata structure from input data
326
327      inpfile%cdjuldref = cdjuldref
328      iobs = 0
329      jm = 0
330      DO jobs = 1, idata
331         DO jcycle = 1, icycles
332            jm = jm + 1
333            IF (llskip(jm)) CYCLE
334            iobs = iobs + 1
335            ! Characters
336            WRITE(inpfile%cdwmo(iobs),'(A3,A5)') cmissions(imission), '     '
337            WRITE(inpfile%cdtyp(iobs),'(I4)') imission
338            ! Real values
339            inpfile%plam(iobs)         = zlam(jobs)
340            inpfile%pphi(iobs)         = zphi(jobs)
341            inpfile%pob(1,iobs,1)      = zsla(jcycle,jobs)
342            inpfile%ptim(iobs)         = zjuld(jm)
343            inpfile%pdep(1,iobs)       = 0.0
344            ! Integers
345            inpfile%kindex(iobs)       = iobs
346            inpfile%ioqc(iobs)      = 1
347            inpfile%ivqc(iobs,1)    = 1
348            inpfile%ivlqc(1,iobs,1) = 1
349            inpfile%ipqc(iobs)         = 0
350            inpfile%ipqcf(:,iobs)      = 0
351            inpfile%itqc(iobs)         = 0
352            inpfile%itqcf(:,iobs)      = 0
353            inpfile%ivqcf(:,iobs,1)    = 0
354            inpfile%ioqcf(:,iobs)      = 0
355            inpfile%idqc(1,iobs)       = 0
356            inpfile%idqcf(1,1,iobs)    = 0
357            inpfile%ivlqcf(:,1,iobs,1) = 0
358         END DO
359      END DO
360
361
362      ! Deallocate memory for input data
363
364      DEALLOCATE( &
365         & iptracks,      &
366         & ipnbpoints,    &
367         & ipcycles,      &
368         & zphi,          &
369         & zlam,          & 
370         & zbegindates,   &
371         & ipdataindexes, &
372         & zsla,          &
373         & zjuld,         &
374         & llskip         &
375         & )
376
377   END SUBROUTINE read_avisofile
378
Note: See TracBrowser for help on using the repository browser.