ecearth3/sources/nemo-3.6/NEMO/OPA_SRC/OBS/obssla_io.h90

   !!----------------------------------------------------------------------
   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
   !! $Id: obssla_io.h90 2287 2010-10-18 07:53:52Z smasson $
   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------

   SUBROUTINE read_avisofile( cdfilename, inpfile, kunit, ldwp, ldgrid )
      !!---------------------------------------------------------------------
      !!
      !!                     ** ROUTINE read_avisofile **
      !!
      !! ** Purpose : Read from file the AVISO SLA observations.
      !!
      !! ** Method  : The data file is a NetCDF file. 
      !!
      !! ** Action  :
      !!
      !! References : http://www.aviso.oceanobs.com
      !!
      !! History : 
      !!          ! 09-01 (K. Mogensen) Original based on old versions
      !!----------------------------------------------------------------------
      !! * Arguments
      CHARACTER(LEN=*) :: cdfilename ! Input filename
      TYPE(obfbdata)   :: inpfile    ! Output obfbdata structure
      INTEGER          :: kunit      ! Unit for output
      LOGICAL          :: ldwp       ! Print info
      LOGICAL          :: ldgrid     ! Save grid info in data structure
      !! * Local declarations
      CHARACTER(LEN=14),PARAMETER :: cpname = 'read_avisofile'
      INTEGER :: i_file_id     ! netcdf IDS
      INTEGER :: i_tracks_id
      INTEGER :: i_cycles_id
      INTEGER :: i_data_id
      INTEGER :: i_var_id
      INTEGER, PARAMETER :: imaxdim = 2    ! Assumed maximum for no. dims. in file
      INTEGER, DIMENSION(2) :: idims       ! Dimensions in file
      INTEGER :: iilen         ! Length of netCDF attributes
      INTEGER :: itype         ! Typeof netCDF attributes
      REAL(fbdp) :: zsca       ! Scale factor
      REAL(fbdp) :: zfill      ! Fill value
      CHARACTER(len=3) ::  cmission      ! Mission global attribute
      INTEGER :: itracks       ! Maximum number of passes in file
      INTEGER :: icycles       ! Maximum number of cycles for each pass
      INTEGER :: idata         ! Number of data per parameter in current file
      REAL(fbdp) :: zdeltat    ! Time gap getween two measurements in seconds
      INTEGER, DIMENSION(:), POINTER :: & 
         & iptracks,     &  ! List of passes contained in current file
         & ipnbpoints,   &  ! Number of points per pass
         & ipdataindexes    ! Index of data in theoretical profile
      INTEGER, DIMENSION(:,:), POINTER :: & 
         & ipcycles         ! List of cycles per pass
      REAL(fbdp), DIMENSION(:), POINTER :: &
         & zphi,   &        ! Latitudes
         & zlam             ! Longitudes
      REAL(fbdp), DIMENSION(:,:), POINTER :: &
         & zbegindates      ! Date of point with index 0
      REAL(fbdp) :: zbeginmiss    ! Missing data for dates
      REAL(fbsp), DIMENSION(:,:), POINTER :: &
         & zsla             ! SLA data
      REAL(fbdp), DIMENSION(:), POINTER :: &
         & zjuld            ! Julian date
      LOGICAL, DIMENSION(:), POINTER :: &
         & llskip           ! Skip observation
      CHARACTER(len=14) :: cdjuldref    ! Julian data reference
      INTEGER :: imission   ! Mission number converted from Mission global 
                            ! netCDF atttribute.
      CHARACTER(len=255) :: ctmp
      INTEGER :: iobs
      INTEGER :: jl
      INTEGER :: jm
      INTEGER :: jj
      INTEGER :: ji
      INTEGER :: jk
      INTEGER :: jobs
      INTEGER :: jcycle

      ! Open the file

      CALL chkerr( nf90_open( TRIM( cdfilename ), nf90_nowrite, i_file_id ), &
         &         cpname, __LINE__ )

      ! Get the netCDF dimensions
      
      CALL chkerr( nf90_inq_dimid( i_file_id, 'Tracks', i_tracks_id ),  &
         &         cpname, __LINE__ )
      CALL chkerr( nf90_inquire_dimension( i_file_id, i_tracks_id, &
         &                                 len = itracks ),  &
         &         cpname, __LINE__ )
      
      CALL chkerr( nf90_inq_dimid( i_file_id, 'Cycles', i_cycles_id ),  &
         &         cpname, __LINE__ )
      CALL chkerr( nf90_inquire_dimension( i_file_id, i_cycles_id, &
         &                                 len = icycles ),  &
         &      cpname, __LINE__ )
      
      CALL chkerr( nf90_inq_dimid( i_file_id, 'Data', i_data_id ),  &
         &         cpname, __LINE__ )
      CALL chkerr( nf90_inquire_dimension( i_file_id, i_data_id, &
         &                                 len = idata ),  &
         &         cpname, __LINE__ )

      ! Allocate memory for input data

      ALLOCATE( &
         & iptracks     ( itracks          ), &  
         & ipnbpoints   ( itracks          ), &  
         & ipcycles     ( icycles, itracks ), &
         & zphi         ( idata            ), &   
         & zlam         ( idata            ), &  
         & zbegindates  ( icycles, itracks ), &
         & ipdataindexes( idata            ), & 
         & zsla         ( icycles, idata   ), &
         & zjuld        ( idata*icycles    ), &
         & llskip       ( idata*icycles    )  &
         & )

      ! Get time gap getween two measurements in seconds
      
      CALL chkerr( nf90_inq_varid( i_file_id, 'DeltaT', i_var_id ), & 
         &         cpname, __LINE__ )
      CALL chkdim( i_file_id, i_var_id, 0, idims, cpname, __LINE__ )
      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, zdeltat), &
         &         cpname, __LINE__ )
      zsca = 1.0
      IF (nf90_inquire_attribute( i_file_id, i_var_id, 'scale_factor') &
         &                      == nf90_noerr) THEN
         CALL chkerr( nf90_get_att( i_file_id, i_var_id, &
            &                     'scale_factor',zsca), cpname,  __LINE__)
      ENDIF
      zdeltat = zsca * zdeltat
         
      ! Get List of passes contained in current file      
      
      CALL chkerr( nf90_inq_varid( i_file_id, 'Tracks', i_var_id ), & 
         &         cpname, __LINE__ )
      idims(1) = itracks
      CALL chkdim( i_file_id, i_var_id, 1, idims, cpname, __LINE__ )
      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, iptracks), &
         &         cpname, __LINE__ )
      
      ! Get Number of points per pass
      
      CALL chkerr( nf90_inq_varid( i_file_id, 'NbPoints', i_var_id ), & 
         &         cpname, __LINE__ )
      idims(1) = itracks
      CALL chkdim( i_file_id, i_var_id, 1, idims, cpname, __LINE__ )
      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, ipnbpoints),&
         &         cpname, __LINE__ )
      
      ! Get list of cycles per pass
      
      CALL chkerr( nf90_inq_varid( i_file_id, 'Cycles', i_var_id ), & 
         &         cpname, __LINE__ )
      idims(1) = icycles
      idims(2) = itracks
      CALL chkdim( i_file_id, i_var_id, 2, idims, cpname, __LINE__ )
      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, ipcycles),   &
         &         cpname, __LINE__ )

      ! Get longitudes

      CALL chkerr( nf90_inq_varid( i_file_id, 'Longitudes', i_var_id ), & 
         &         cpname, __LINE__ )
      idims(1) = idata
      CALL chkdim( i_file_id, i_var_id, 1, idims, cpname, __LINE__ )
      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, zlam),   &
         &         cpname, __LINE__ )
      zsca = 1.0
      IF (nf90_inquire_attribute( i_file_id, i_var_id, 'scale_factor') &
         &                         == nf90_noerr) THEN
         CALL chkerr( nf90_get_att( i_file_id, i_var_id, &
            &                     'scale_factor',zsca), cpname,  __LINE__)
      ENDIF
      zlam(:) = zsca * zlam(:)
      
      ! Get latitudes
      
      CALL chkerr( nf90_inq_varid( i_file_id, 'Latitudes', i_var_id ), & 
         &         cpname, __LINE__ )
      idims(1) = idata
      CALL chkdim( i_file_id, i_var_id, 1, idims, cpname, __LINE__ )
      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, zphi),   &
         &         cpname, __LINE__ )
      zsca = 1.0
      IF (nf90_inquire_attribute( i_file_id, i_var_id, 'scale_factor') &
         &                       == nf90_noerr) THEN
         CALL chkerr( nf90_get_att( i_file_id, i_var_id, &
            &                     'scale_factor',zsca), cpname,  __LINE__)
      ENDIF
      zphi(:) = zsca * zphi(:)
      
      ! Get date of point with index 0
      
      CALL chkerr( nf90_inq_varid( i_file_id, 'BeginDates', i_var_id ), & 
         &        cpname, __LINE__ )
      idims(1) = icycles
      idims(2) = itracks
      CALL chkdim( i_file_id, i_var_id, 2, idims, cpname, __LINE__ )
      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, zbegindates), &
         &         cpname, __LINE__ )
      CALL chkerr( nf90_inquire_attribute( i_file_id, i_var_id, &
         &                                 'units', len = iilen,  &
         &                                  xtype = itype), cpname, __LINE__ )
      IF (( itype /= NF90_CHAR ) .OR. ( iilen > 255 )) THEN
         CALL fatal_error('Error decoding BeginDates unit', __LINE__ )
      ENDIF
      CALL chkerr( nf90_get_att( i_file_id, i_var_id, 'units', &
         &                       ctmp ), cpname, __LINE__ )
      jl=1
      DO jm = 1, LEN(TRIM(ctmp))
         IF ((ctmp(jm:jm)>='0').AND.(ctmp(jm:jm)<='9')) THEN
            cdjuldref(jl:jl) = ctmp(jm:jm)
            jl=jl+1
         ENDIF
         IF (jl>14) EXIT
      END DO
      CALL chkerr( nf90_inquire_attribute( i_file_id, i_var_id, '_FillValue', &
         &                                  xtype = itype), cpname, __LINE__ )
      IF ( itype /= NF90_DOUBLE ) THEN
         CALL fatal_error('Error decoding BeginDates missing data', __LINE__ )
      ENDIF
      CALL chkerr( nf90_get_att( i_file_id, i_var_id, '_FillValue', &
         &                       zbeginmiss ), cpname, __LINE__ )

      ! Get indices of data in theoretical profile
      
      CALL chkerr( nf90_inq_varid( i_file_id, 'DataIndexes', i_var_id ), & 
         &         cpname, __LINE__ )
      idims(1) = idata
      CALL chkdim( i_file_id, i_var_id, 1, idims, cpname, __LINE__ )
      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, ipdataindexes),   &
         &       cpname, __LINE__ )
      
      ! Get SLA data
      
      CALL chkerr( nf90_inq_varid( i_file_id, 'SLA', i_var_id ), & 
         &         cpname, __LINE__ )
      idims(1) = icycles
      idims(2) = idata
      CALL chkdim( i_file_id, i_var_id, 2, idims, cpname, __LINE__ )
      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, zsla),   &
         &         cpname, __LINE__ )
      zsca = 1.0
      IF (nf90_inquire_attribute( i_file_id, i_var_id, 'scale_factor') &
         &                        == nf90_noerr) THEN
         CALL chkerr( nf90_get_att( i_file_id, i_var_id, &
            &                       'scale_factor',zsca), cpname, __LINE__ )
      ENDIF
      zfill = 0.0
      IF (nf90_inquire_attribute( i_file_id, i_var_id, '_FillValue') &
         &                      == nf90_noerr) THEN
         CALL chkerr( nf90_get_att( i_file_id, i_var_id, &
            &                     '_FillValue',zfill), cpname,  __LINE__ )
      ENDIF
      WHERE(zsla(:,:) /=  zfill)
         zsla(:,:) = zsca * zsla(:,:)
      ELSEWHERE
         zsla(:,:) = fbrmdi
      END WHERE
      
      ! Get the global Mission netCDF attribute
      
      cmission='   '
      CALL chkerr( nf90_inquire_attribute( i_file_id, nf90_global, &
         &                                 'Mission', len = iilen,  &
         &                                 xtype = itype), cpname, __LINE__ )
      IF (( itype /= NF90_CHAR ) .OR. ( iilen > 3 )) THEN
         CALL fatal_error('Error decoding Mission global attribute', __LINE__ )
      ENDIF
      CALL chkerr( nf90_get_att( i_file_id, nf90_global, 'Mission', &
         &                       cmission ), cpname, __LINE__ )
      
      ! Convert it to an integer
      imission = 0
      DO jm = 1, imaxmissions
         IF (cmission==cmissions(jm)) THEN
            imission = jm
            EXIT
         ENDIF
      END DO
      
      ! Close the file
      
      CALL chkerr( nf90_close( i_file_id ), cpname, __LINE__ )

      ! Compute Julian dates for all observations
      
      jm = 0
      jl = 0
      DO jj = 1, itracks
         DO ji = 1, ipnbpoints(jj)
            jl = jl + 1
            DO jk = 1, icycles
               jm = jm + 1
               IF (zbegindates(jk,jj)==zbeginmiss) THEN
                  llskip(jm) = .TRUE.
                  zjuld(jm)  = fbrmdi
               ELSE
                  llskip(jm) = .FALSE.
                  zjuld(jm)  = zbegindates(jk,jj) + &
                     &         (ipdataindexes(jl) * zdeltat / 86400._wp )
               ENDIF
            END DO
         END DO
      END DO

      ! Get rid of missing data

      jm = 0
      DO jobs = 1, idata
         DO jcycle = 1, icycles
            jm = jm + 1
            IF (zsla(jcycle,jobs) == fbrmdi) llskip(jm) = .TRUE.
         END DO
      END DO
      
      ! Allocate obfbdata

      iobs = COUNT( .NOT.llskip(:) )
      CALL init_obfbdata( inpfile )
      CALL alloc_obfbdata( inpfile, 1, iobs, 1, 0, 0, ldgrid )
      inpfile%cname(1) = 'SLA'

      ! Fill the obfbdata structure from input data

      inpfile%cdjuldref = cdjuldref
      iobs = 0
      jm = 0
      DO jobs = 1, idata
         DO jcycle = 1, icycles
            jm = jm + 1
            IF (llskip(jm)) CYCLE
            iobs = iobs + 1
            ! Characters
            WRITE(inpfile%cdwmo(iobs),'(A3,A5)') cmissions(imission), '     '
            WRITE(inpfile%cdtyp(iobs),'(I4)') imission
            ! Real values
            inpfile%plam(iobs)         = zlam(jobs)
            inpfile%pphi(iobs)         = zphi(jobs)
            inpfile%pob(1,iobs,1)      = zsla(jcycle,jobs)
            inpfile%ptim(iobs)         = zjuld(jm)
            inpfile%pdep(1,iobs)       = 0.0
            ! Integers
            inpfile%kindex(iobs)       = iobs
            inpfile%ioqc(iobs)      = 1
            inpfile%ivqc(iobs,1)    = 1
            inpfile%ivlqc(1,iobs,1) = 1
            inpfile%ipqc(iobs)         = 0 
            inpfile%ipqcf(:,iobs)      = 0
            inpfile%itqc(iobs)         = 0
            inpfile%itqcf(:,iobs)      = 0
            inpfile%ivqcf(:,iobs,1)    = 0
            inpfile%ioqcf(:,iobs)      = 0
            inpfile%idqc(1,iobs)       = 0
            inpfile%idqcf(1,1,iobs)    = 0
            inpfile%ivlqcf(:,1,iobs,1) = 0
         END DO
      END DO


      ! Deallocate memory for input data

      DEALLOCATE( &
         & iptracks,      &
         & ipnbpoints,    &
         & ipcycles,      &
         & zphi,          &
         & zlam,          &  
         & zbegindates,   &
         & ipdataindexes, &
         & zsla,          &
         & zjuld,         &
         & llskip         &
         & )

   END SUBROUTINE read_avisofile