#ifdef IBM
@PROCESS NOOPT
#endif
MODULE bufrdata

#ifdef HAS_BUFR

   USE toolspar_kind
   USE obs_utils
   USE obs_fbm

   IMPLICIT NONE

CONTAINS

   SUBROUTINE read_bufrfile( cdfilename, inptsm, inpts1, inpuvm, inpuv1 )
      !!---------------------------------------------------------------------
      !!
      !!                     ** ROUTINE read_bufrfile **
      !!
      !!  ** Purpose : Read from file the profile BUFR observations.
      !!
      !!  ** Method  : The data file is a BUFR file. 
      !!
      !!  ** Action  :
      !!
      !!   History :
      !!----------------------------------------------------------------------
      !! * Arguments
      CHARACTER(LEN=*) :: cdfilename  ! Input file.
      TYPE(obfbdata) :: inptsm        ! Output TS data structure (prof).
      TYPE(obfbdata) :: inpts1        ! Output TS data structure (single).
      TYPE(obfbdata) :: inpuvm        ! Output UV data structure (prof).
      TYPE(obfbdata) :: inpuv1        ! Output UV data structure (single).
      !! * Local variables
      INTEGER :: iunit, iret
      INTEGER :: ibytesize, isubtype, itype
      INTEGER :: iyea, imon, iday, ihou, imin
      INTEGER :: jk
      REAL :: zlat,zlon
      CHARACTER(len=8) :: csign
      INTEGER, PARAMETER :: imaxlev = 300
      INTEGER :: ilevts, iobstsm, iobsts1, ilevtsmax
      INTEGER :: ilevuv, iobsuvm, iobsuv1, ilevuvmax
      INTEGER :: iobs
      REAL, DIMENSION(imaxlev) :: zdts, ztem, zsal, zduv, zuv, zvv
      real :: zmdi

      ! Initialize data
      
      ibytesize = BIT_SIZE(1)/8
      zmdi = REAL(fbrmdi)
      
      ! Open file
      CALL pbopen( iunit, cdfilename, 'r', iret )
      IF ( iret /= 0 ) THEN
         WRITE(*,*)'Error opening file : ', TRIM(cdfilename)
         WRITE(*,*)'Error code : ', iret
         CALL abort
      ENDIF

      ! Count number of observations

      ilevtsmax = -1
      ilevuvmax = -1
      iobstsm = 0
      iobsts1 = 0
      iobsuvm = 0
      iobsuv1 = 0

      DO

         CALL decode_bufr( &
            & iunit, ibytesize,             &
            & iyea, imon, iday, ihou, imin, &
            & zlat, zlon,                   &
            & csign, isubtype,              &
            & imaxlev, ilevts, ilevuv,      & 
            & zdts, ztem, zsal,             &
            & zduv, zuv, zvv,               &
            & zmdi, iret )
         
         ! No more observations if iret==-1

         IF ( iret == -1 ) EXIT

         ! Ignore observations without levels

         IF ( ilevts > 0 ) THEN 

            IF ( ilevts == 1 ) THEN
               iobsts1 = iobsts1 + 1
            ELSE
               iobstsm = iobstsm + 1
               IF ( ilevts > ilevtsmax ) ilevtsmax = ilevts
            ENDIF

         ENDIF

         IF ( ilevuv > 0 ) THEN 

            IF ( ilevuv == 1 ) THEN
               iobsuv1 = iobsuv1 + 1
            ELSE
               iobsuvm = iobsuvm + 1
               IF ( ilevuv > ilevuvmax ) ilevuvmax = ilevuv
            ENDIF

         ENDIF

      ENDDO

      ! Setup data structures

      CALL init_obfbdata( inptsm )
      CALL alloc_obfbdata( inptsm, 2, iobstsm, ilevtsmax, 0, 1, .FALSE. )
      inptsm%cname(1) = 'POTM'
      inptsm%cname(2) = 'PSAL'
      inptsm%coblong(1) = 'Potential temperature'
      inptsm%coblong(2) = 'Practical salinity'
      inptsm%cobunit(1) = 'Degrees Celsius'
      inptsm%cobunit(2) = 'PSU'
      inptsm%cextname(1) = 'TEMP'
      inptsm%cextlong(1) = 'Insitu temperature'
      inptsm%cextunit(1) = 'Degrees Celsius'
      CALL init_obfbdata( inpts1 )
      CALL alloc_obfbdata( inpts1, 2, iobsts1, 1, 0, 1, .FALSE. )
      inpts1%cname(1) = 'POTM'
      inpts1%cname(2) = 'PSAL'
      inpts1%coblong(1) = 'Potential temperature'
      inpts1%coblong(2) = 'Practical salinity'
      inpts1%cobunit(1) = 'Degrees Celsius'
      inpts1%cobunit(2) = 'PSU'
      inpts1%cextname(1) = 'TEMP'
      inpts1%cextlong(1) = 'Insitu temperature'
      inpts1%cextunit(1) = 'Degrees Celsius'
      CALL init_obfbdata( inpuvm )
      CALL alloc_obfbdata( inpuvm, 2, iobsuvm, ilevuvmax, 0, 0, .FALSE. )
      inpuvm%cname(1) = 'UVEL'
      inpuvm%cname(2) = 'VVEL'
      inpuvm%coblong(1) = 'Zonal current'
      inpuvm%coblong(2) = 'Meridional current'
      inpuvm%cobunit(1) = 'Meters per second'
      inpuvm%cobunit(2) = 'Meters per second'
      CALL init_obfbdata( inpuv1 )
      CALL alloc_obfbdata( inpuv1, 2, iobsuv1, 1, 0, 0, .FALSE. )
      inpuv1%cname(1) = 'UVEL'
      inpuv1%cname(2) = 'VVEL'
      inpuv1%coblong(1) = 'Zonal current'
      inpuv1%coblong(2) = 'Meridional current'
      inpuv1%cobunit(1) = 'Meters per second'
      inpuv1%cobunit(2) = 'Meters per second'

      ! Rewind the file

      CALL pbseek( iunit, 0, 0, iret )
      IF ( iret /= 0 ) THEN
         WRITE(*,*)'Error rewinding file : ',TRIM(cdfilename)
         WRITE(*,*)'Error code : ',iret
         CALL abort
      ENDIF
      
      ! Retrieve the observations and put them into the data structure

      iobstsm = 0
      iobsts1 = 0
      iobsuvm = 0
      iobsuv1 = 0
      iret = 0
      inptsm%cdjuldref = '19500101000000'
      inpts1%cdjuldref = '19500101000000'
      inpuvm%cdjuldref = '19500101000000'
      inpuv1%cdjuldref = '19500101000000'
      iobs = 0
      DO

         CALL decode_bufr( &
            & iunit, ibytesize,             &
            & iyea, imon, iday, ihou, imin, &
            & zlat, zlon,                   &
            & csign, isubtype,              &
            & imaxlev, ilevts, ilevuv,      & 
            & zdts, ztem, zsal,             &
            & zduv, zuv, zvv,               &
            & zmdi, iret )

         ! No more observations if iret==-1

         IF ( iret == -1 ) EXIT
         iobs = iobs + 1
         
         ! Convert isubtype to something EN3 like.

         IF ( isubtype == 131 ) THEN
            itype = 820
         ELSEIF ( isubtype == 132 ) THEN
            itype = 401
         ELSEIF ( isubtype == 133 ) THEN
            IF (LEN_TRIM(CSIGN)>=7) THEN
               itype = 831
            ELSE
               itype = 741
            ENDIF
         ELSE
            WRITE(*,*)'Unknown subtype = ',isubtype
            itype = isubtype
         ENDIF

         ! TS data

         ! Ignore observations without levels

         IF ( ilevts > 0 ) THEN

            IF ( ilevts==1 ) THEN

               iobsts1 = iobsts1 + 1

               ! Position and time
               
               inpts1%kindex(iobsts1) = iobs
               inpts1%pphi(iobsts1) = zlat
               inpts1%plam(iobsts1) = zlon
               CALL greg2jul( 0, imin, ihou, iday, imon, iyea, &
                  &           inpts1%ptim(iobsts1) )

               ! Call sign and type
               
               inpts1%cdwmo(iobsts1) = csign
               WRITE(inpts1%cdtyp(iobsts1),'(I4.4)') itype

               ! Depth, Salinity and Insitu Temperature
               
               DO jk = 1, ilevts
                  
                  inpts1%pdep(jk,iobsts1) = zdts(jk)
                  inpts1%pext(jk,iobsts1,1) = ztem(jk)
                  inpts1%pob(jk,iobsts1,2) = zsal(jk)
                  
                  ! If no insitu temperature set the QC flag to 4
                  
                  IF ( ztem(jk) == zmdi ) THEN
                     inpts1%ivlqc(jk,iobsts1,1) = 4
                  ELSE
                     inpts1%ivlqc(jk,iobsts1,1) = 1
                  ENDIF

                  ! If no salinity set the QC flag to 4
                  
                  IF ( zsal(jk) == zmdi ) THEN
                     inpts1%ivlqc(jk,iobsts1,2) = 4
                  ELSE
                     inpts1%ivlqc(jk,iobsts1,2) = 1
                  ENDIF
                  
               ENDDO

            ELSE

               iobstsm = iobstsm + 1
               
               ! Position and time
               
               inptsm%kindex(iobstsm) = iobs
               inptsm%pphi(iobstsm) = zlat
               inptsm%plam(iobstsm) = zlon
               CALL greg2jul( 0, imin, ihou, iday, imon, iyea, &
                  &           inptsm%ptim(iobstsm) )

               ! Call sign and type
               
               inptsm%cdwmo(iobstsm) = csign
               WRITE(inptsm%cdtyp(iobstsm),'(I4.4)') itype

               ! Depth, Salinity and Insitu Temperature
               
               DO jk = 1, ilevts
               
                  inptsm%pdep(jk,iobstsm) = zdts(jk)
                  inptsm%pext(jk,iobstsm,1) = ztem(jk)
                  inptsm%pob(jk,iobstsm,2) = zsal(jk)

                  ! If no insitu temperature set the QC flag to 4
                  
                  IF ( ztem(jk) == zmdi ) THEN
                     inptsm%ivlqc(jk,iobstsm,1) = 4
                  ELSE
                     inptsm%ivlqc(jk,iobstsm,1) = 1
                  ENDIF

                  ! If no salinity set the QC flag to 4
                  
                  IF ( zsal(jk) == zmdi ) THEN
                     inptsm%ivlqc(jk,iobstsm,2) = 4
                  ELSE
                     inptsm%ivlqc(jk,iobstsm,2) = 1
                  ENDIF
                  
               ENDDO
               
            ENDIF

         ENDIF

         ! UV data

         ! Ignore observations without levels

         IF ( ilevuv > 0 ) THEN

            ! Skip surface only observations (TBC).

            IF ( ilevuv == 1 ) THEN

               iobsuv1 = iobsuv1 + 1
               
               ! Position and time
               
               inpuv1%kindex(iobsuv1) = iobs
               inpuv1%pphi(iobsuv1) = zlat
               inpuv1%plam(iobsuv1) = zlon
               CALL greg2jul( 0, imin, ihou, iday, imon, iyea, &
                  &           inpuv1%ptim(iobsuv1) )
               
               ! Call sign and type
               
               inpuv1%cdwmo(iobsuv1) = csign
               WRITE(inpuv1%cdtyp(iobsuv1),'(I4.4)') itype

               ! Depth, Salinity and Insitu Temperature
               
               DO jk = 1, ilevuv
                  
                  inpuv1%pdep(jk,iobsuv1)  =  zduv(jk)
                  inpuv1%pob(jk,iobsuv1,1)  =  zuv(jk)
                  inpuv1%pob(jk,iobsuv1,2)  =  zvv(jk)
                  
                  ! If no insitu temperature set the QC flag to 4
                  
                  IF ( zuv(jk) == zmdi ) THEN
                     inpuv1%ivlqc(jk,iobsuv1,1) = 4
                  ELSE
                     inpuv1%ivlqc(jk,iobsuv1,1) = 1
                  ENDIF
                  
                  ! If no v velocity set the QC flag to 4
                  
                  IF ( zvv(jk) == zmdi ) THEN
                     inpuv1%ivlqc(jk,iobsuv1,2) = 4
                  ELSE
                     inpuv1%ivlqc(jk,iobsuv1,2) = 1
                  ENDIF
                  
               ENDDO
               
            ELSE

               iobsuvm = iobsuvm + 1
               
               ! Position and time
               
               inpuvm%kindex(iobsuvm) = iobs
               inpuvm%pphi(iobsuvm) = zlat
               inpuvm%plam(iobsuvm) = zlon
               CALL greg2jul( 0, imin, ihou, iday, imon, iyea, &
                  &           inpuvm%ptim(iobsuvm) )
               
               ! Call sign and type
               
               inpuvm%cdwmo(iobsuvm) = csign
               WRITE(inpuvm%cdtyp(iobsuvm),'(I4.4)') itype

               ! Depth, Salinity and Insitu Temperature
               
               DO jk = 1, ilevuv
                  
                  inpuvm%pdep(jk,iobsuvm)  =  zduv(jk)
                  inpuvm%pob(jk,iobsuvm,1)  =  zuv(jk)
                  inpuvm%pob(jk,iobsuvm,2)  =  zvv(jk)
                  
                  ! If no insitu temperature set the QC flag to 4
                  
                  IF ( zuv(jk) == zmdi ) THEN
                     inpuvm%ivlqc(jk,iobsuvm,1) = 4
                  ELSE
                     inpuvm%ivlqc(jk,iobsuvm,1) = 1
                  ENDIF
                  
                  ! If no v velocity set the QC flag to 4
                  
                  IF ( zvv(jk) == zmdi ) THEN
                     inpuvm%ivlqc(jk,iobsuvm,2) = 4
                  ELSE
                     inpuvm%ivlqc(jk,iobsuvm,2) = 1
                  ENDIF
                  
               ENDDO

            ENDIF


         ENDIF


      ENDDO
      
      ! Close the file

      CALL pbclose( iunit, iret )
      IF ( iret /= 0) THEN
         WRITE(*,*)'Error closing file : ',TRIM(cdfilename)
         WRITE(*,*)'Error code : ',iret
         CALL abort
      ENDIF
      
      ! Set the QC flags which can not be read from BUFR to 0

      inptsm%ipqc(:)         = 0
      inptsm%ipqcf(:,:)      = 0
      inptsm%ivqc(:,:)       = 0
      inptsm%ivqcf(:,:,:)    = 0
      inptsm%ivqc(:,:)       = 0
      inptsm%ivqcf(:,:,: )   = 0
      inptsm%ioqc(:)         = 0
      inptsm%ioqcf(:,:)      = 0
      inptsm%idqc(:,:)       = 0
      inptsm%idqcf(:,:,:)    = 0
      inptsm%ivlqcf(:,:,:,:) = 0
      inpts1%ipqc(:)         = 0
      inpts1%ipqcf(:,:)      = 0
      inpts1%ivqc(:,:)       = 0
      inpts1%ivqcf(:,:,:)    = 0
      inpts1%ivqc(:,:)       = 0
      inpts1%ivqcf(:,:,: )   = 0
      inpts1%ioqc(:)         = 0
      inpts1%ioqcf(:,:)      = 0
      inpts1%idqc(:,:)       = 0
      inpts1%idqcf(:,:,:)    = 0
      inpts1%ivlqcf(:,:,:,:) = 0
      inpuvm%ipqc(:)         = 0
      inpuvm%ipqcf(:,:)      = 0
      inpuvm%ivqc(:,:)       = 0
      inpuvm%ivqcf(:,:,:)    = 0
      inpuvm%ivqc(:,:)       = 0
      inpuvm%ivqcf(:,:,: )   = 0
      inpuvm%ioqc(:)         = 0
      inpuvm%ioqcf(:,:)      = 0
      inpuvm%idqc(:,:)       = 0
      inpuvm%idqcf(:,:,:)    = 0
      inpuvm%ivlqcf(:,:,:,:) = 0
      inpuv1%ipqc(:)         = 0
      inpuv1%ipqcf(:,:)      = 0
      inpuv1%ivqc(:,:)       = 0
      inpuv1%ivqcf(:,:,:)    = 0
      inpuv1%ivqc(:,:)       = 0
      inpuv1%ivqcf(:,:,: )   = 0
      inpuv1%ioqc(:)         = 0
      inpuv1%ioqcf(:,:)      = 0
      inpuv1%idqc(:,:)       = 0
      inpuv1%idqcf(:,:,:)    = 0
      inpuv1%ivlqcf(:,:,:,:) = 0

   END SUBROUTINE read_bufrfile

   SUBROUTINE decode_bufr( kunit, kwsize, &
      & kyea, kmon, kday, khou, kmin,     &
      & plat,plon, csign, ksty,           &
      & kmlevs, klevts, klevve,           &
      & pdts, ptem, psal,                 &
      & pdve, puve, pvve,                 &
      & pvalnull, ierr)
      !!---------------------------------------------------------------------
      !!
      !!                     ** ROUTINE decode_bufr **
      !!
      !!  ** Purpose : Decode a singe BUFR record for ocean data
      !!
      !!  ** Method  : Call to BUFREX/BUSEL in EMOS lib.
      !!
      !!  ** Action  :
      !!
      !!   History :   (??-??) A. Vidard  . ODASYS version. 
      !!               (08-12) K. Mogensen. NEMOVAR version
      !!----------------------------------------------------------------------
      !! * Arguments
      INTEGER, INTENT(IN) :: kunit      ! Bufr file unit number
      INTEGER, INTENT(IN) :: kwsize     ! No of bytes in a word
      INTEGER, INTENT(OUT) :: &         ! Date
         & kyea, kmon, kday, khou, kmin
      REAL, INTENT(OUT) :: &            ! Position
         & plat, plon
      CHARACTER(LEN=*), INTENT(OUT) ::  &
         & csign                        ! Call sign
      INTEGER, INTENT(OUT) :: ksty      ! Subtype
      INTEGER, INTENT(IN) :: kmlevs     ! Maximum number of levels
      INTEGER, INTENT(OUT) :: klevts    ! Actual number of levels
      INTEGER, INTENT(OUT) :: klevve    ! Actual number of levels
      REAL, DIMENSION(kmlevs), INTENT(OUT) :: &
         & pdts, &                      ! Depths for t,s
         & ptem, &                      ! Insitu temperature
         & psal, &                      ! Salinity
         & pdve, &                      ! Depths for currents
         & puve, &                      ! U velocity
         & pvve                         ! V velocity
      REAL, INTENT(IN) :: pvalnull      ! Missing value
      INTEGER, INTENT(INOUT) :: ierr    ! Error code
      !! * Define BUFR parameters
      INTEGER, PARAMETER :: imax_bufr_bytes=80000 ! Max bufr length in bytes
      INTEGER, PARAMETER :: imax_elements=20000   ! Kax no of elements
      INTEGER, PARAMETER :: imax_values=20000     ! Max no of values 
      INTEGER, PARAMETER :: imax_sec0=3           ! Length of section 0
      INTEGER, PARAMETER :: imax_sec1=40          ! Length of section 1
      INTEGER, PARAMETER :: imax_sec2=64          ! Length of section 2
      INTEGER, PARAMETER :: imax_sec3=4           ! Length of section 3
      INTEGER, PARAMETER :: imax_sec4=2           ! Length of section 4
      INTEGER, PARAMETER :: idim_ksup=9           ! Length of array ksup
      !! * Define local variables
      INTEGER :: ksup(idim_ksup)        ! Expanded bufr additional information
      INTEGER :: ksec0(imax_sec0)       ! Expanded bufr section 0
      INTEGER :: ksec1(imax_sec1)       ! Expanded bufr section 1
      INTEGER :: ksec2(imax_sec2)       ! Expanded bufr section 2
      INTEGER :: ksec3(imax_sec3)       ! Expanded bufr section 3
      INTEGER :: ksec4(imax_sec4)       ! Expanded bufr section 4
      INTEGER :: kdtlen,ktdexl          ! Bufr lengths
      INTEGER :: &                      ! BUSEL variables
         & ktdlst(imax_elements), ktdexp(imax_elements)
      INTEGER :: kbuff(imax_bufr_bytes) ! Bufr message array 
      REAL :: values(imax_values)       ! Bufr message values
      CHARACTER(LEN=64) :: &
         & cnames(imax_elements), &     ! Element names
         & cunits(imax_elements), &     ! Element units
         & cvals(imax_values)           ! Character values
      REAL, DIMENSION((kmlevs)) :: &
         & zspd, &                      ! Current speed
         & zdir                         ! Current direction
      INTEGER :: ilenbuf                ! Word len. of compressed bufr message
      INTEGER :: kerr                   ! Error return code
      INTEGER :: ibytes                 ! Actual no. of bytes in bufr message
      INTEGER :: ji                     ! Loop counter
      INTEGER :: icur                   ! Pointer to current data
      ! Define physical constants
      REAL, PARAMETER :: zero_celcius = 273.15

      ! Initialize variables

      cvals(1) = '         '

      DO ji = 1, kmlevs
         pdts(ji) = pvalnull
         ptem(ji) = pvalnull
         psal(ji) = pvalnull
         puve(ji) = pvalnull
         pvve(ji) = pvalnull
      END DO
      
      DO ji = 1, imax_elements
         ktdlst(ji) = 0
         ktdexp(ji) = 0
      END DO
      values(:) = pvalnull

      ! Read bufr message

      kerr = 0
      CALL pbbufr( kunit, kbuff, imax_bufr_bytes, ibytes, kerr) 
      IF ( kerr == -1 ) THEN
         ierr=-1
         RETURN
      ELSE IF ( kerr < -1 ) THEN
         WRITE(*,*)'ERROR in PBBUFR with code ',kerr
         CALL abort
      ENDIF
      
      ! Decode and expand bufr message

      ! Calculate length of bufr message in words
      ilenbuf = ibytes/kwsize+1
      
      ! Expand bufr message

      kerr=0
      CALL bufrex( ilenbuf, kbuff, ksup, ksec0, ksec1, ksec2, ksec3, ksec4,   &
         &         imax_elements, cnames, cunits, imax_values, values, cvals, &
         &         kerr)
      IF ( kerr /= 0 ) THEN
         WRITE(*,*)'ERROR in BUFREX with code ',kerr
         klevts = 0
         klevve = 0
         RETURN
      ENDIF
      CALL busel( kdtlen, ktdlst, ktdexl, ktdexp, kerr )
      IF ( kerr /= 0 ) THEN
         WRITE(*,*)'ERROR in BUSEL with code ',kerr
         klevts = 0
         klevve = 0
         RETURN
      ENDIF

      ! Set oceanographic variables

      ksty=ksec1(7)
      kyea=ksec1(9)

      ! Patch for 2 digit only date (only matter if > 2000)

      IF ( kyea < 40 ) kyea = 2000 + kyea

      SELECT CASE(ksty) 
         
      CASE(131)
         
         ! DRIBU

         kyea = values(5)
         kmon = values(6)
         kday = values(7)
         khou = values(8)
         kmin = values(9)
         plat = values(10)
         plon = values(11)

         ! Check that this dribu contains salinity

         IF ( ktdexp(16) == 2033 ) THEN

            ! Get number of levels
            IF ( values(17) < 100000. ) THEN
               klevts = INT( values(17) )
            ELSE
               ! Patch for nasty data.. without data (AV 200307)
               klevts = 0
            ENDIF

         ELSE
            
            klevts = 0

         ENDIF

         ! Patch for missing WMO number (AV 200307)
         IF ( values(1) >=  9999999. ) THEN
            csign = 'XXX'
         ELSE
            IF ( INT( values(1) ) .LT. 100000 ) THEN
               WRITE(csign(1:5),'(i5)') INT( values(1) )
            ELSE
               WRITE(csign,     '(i8)') INT( values(1) )
            ENDIF
         ENDIF

         CALL bufr_get_data( 7062, imax_elements, ktdexp, &
            &                imax_values, values, klevts, 18, 3, pdts )
         
         CALL bufr_get_data( 22043, imax_elements, ktdexp, &
            &                imax_values, values, klevts, 19, 3, ptem )
         
         CALL bufr_get_data( 22062, imax_elements, ktdexp, &
            &                imax_values, values, klevts, 20, 3, psal )

         IF ( ktdexp(15) == 2031 ) THEN

            klevve = INT( values(16) )

            CALL bufr_get_data( 7062, imax_elements, ktdexp, &
               &                imax_values, values, klevve, &
               &                17, 3, pdve )
            
            CALL bufr_get_data( 22004, imax_elements, ktdexp, &
               &                imax_values, values, klevve, &
               &                18, 3, zdir )
            
            CALL bufr_get_data( 22031, imax_elements, ktdexp, &
               &                imax_values, values, klevve, &
               &                19, 3, zspd )

         ELSE
            
            klevve = 0

         ENDIF


      CASE(132)
         
         ! BATHY

         kyea = values(2)
         kmon = values(3)
         kday = values(4)
         khou = values(5)
         kmin = values(6)
         plat = values(7)
         plon = values(8)
         
         csign = cvals(1)
         
         klevts = INT( values(10) )

         CALL bufr_get_data( 7062, imax_elements, ktdexp, &
            &                imax_values, values, klevts, 11, 2, pdts )
         
         CALL bufr_get_data( 22042, imax_elements, ktdexp, &
            &                imax_values, values, klevts, 12, 2, ptem )
         
         psal=pvalnull

         klevve = 0

      CASE(133)
         
         ! TESAC

         kyea = values(2)
         kmon = values(3)
         kday = values(4)
         khou = values(5)
         kmin = values(6)
         plat = values(7)
         plon = values(8)
         
         csign = cvals(1)
         
         klevts = INT( values(11) )
         
         CALL bufr_get_data( 7062, imax_elements, ktdexp, &
            &                imax_values, values, klevts, 12, 3, pdts )
         
         CALL bufr_get_data( 22043, imax_elements, ktdexp, &
            &                imax_values, values, klevts, 13, 3, ptem )
         
         CALL bufr_get_data( 22062, imax_elements, ktdexp, &
            &                imax_values, values, klevts, 14, 3, psal )

         icur = 14 + 3 * (klevts - 1 ) + 1
         
         IF ( ktdexp(icur) == 2031 ) THEN

            klevve = INT( values(icur+1) )

            CALL bufr_get_data( 7062, imax_elements, ktdexp, &
               &                imax_values, values, klevve, &
               &                icur + 2, 3, pdve )
            
            CALL bufr_get_data( 22004, imax_elements, ktdexp, &
               &                imax_values, values, klevve, &
               &                icur + 3, 3, zdir )
            
            CALL bufr_get_data( 22031, imax_elements, ktdexp, &
               &                imax_values, values, klevve, &
               &                icur + 4, 3, zspd )

         ELSE

            klevve = 0

         ENDIF
         
      CASE DEFAULT
         
         WRITE(*,*)'ERROR: SUBTYPE ',ksty, &
            &           ' is not recognised as an Ocean type'
         CALL abort
         
      END SELECT
      
      ! Convert data units

      DO ji = 1, klevts
         ! Convert temperature from Kelvin to Celcius
         IF (( ptem(ji) > 250. ) .AND. ( ptem(ji) < 323. )) THEN
            ptem(ji) = ptem(ji) - zero_celcius
         ELSE
            ptem(ji) = pvalnull
         ENDIF
         ! Check for physical salinities
         IF (( psal(ji) < 0.0 ) .OR. ( psal(ji) > 50.0 )) THEN
            psal(ji) = pvalnull
         ENDIF
      END DO

      DO ji = 1, klevve
         ! Convert speed and direction to u,v
         IF (( zspd(ji) < 100 ).AND.( zdir(ji) < 361. )) THEN
            CALL spd_to_uv( zdir(ji), zspd(ji), puve(ji), pvve(ji) )
         ELSE
            puve(ji) = pvalnull
            pvve(ji) = pvalnull
         ENDIF
      END DO

   END SUBROUTINE decode_bufr


   SUBROUTINE bufr_get_data( ktdval, ktdlen, ktdexp, &
      &                      kmax_values, zvalues, klevs, &
      &                      istart, istep, zovalues )
      !!---------------------------------------------------------------------
      !!
      !!                     ** ROUTINE bufr_get_data **
      !!
      !!  ** Purpose : Extract data from decode buffer message
      !!
      !!  ** Method  : 
      !!
      !!  ** Action  :
      !!
      !!   History :   (??-??) A. Vidard  . ODASYS version. 
      !!               (08-12) K. Mogensen. NEMOVAR version
      !!----------------------------------------------------------------------
      !! * Arguments
      INTEGER, INTENT(IN) :: ktdval                ! Expected descriptor
      INTEGER, INTENT(IN) :: ktdlen                ! Size of descriptors
      INTEGER, INTENT(IN), DIMENSION(ktdlen) :: &
         & ktdexp                                  ! Descriptors
      INTEGER, INTENT(IN) :: kmax_values           ! Size of input array
      REAL, INTENT(IN) :: zvalues(kmax_values)     ! Decoded values
      INTEGER, INTENT(IN) :: klevs                 ! Number of levels 
      INTEGER, INTENT(IN) :: istart                ! Start in input array
      INTEGER, INTENT(IN) :: istep                 ! Step between values
      REAL, INTENT(OUT) :: zovalues(klevs)         ! Output values
      !! * Local variables
      INTEGER :: ji, ipos

      DO ji = 1, klevs
         ipos = istart + (ji-1) * istep
         IF ( ktdexp(ipos) /= ktdval ) THEN
            WRITE(*,*)'Problem decoding bufr data in bufr_get_data'
            WRITE(*,*)'Expected descriptor : ', ktdval
            WRITE(*,*)'Found descriptor : ', ktdexp(ipos)
            CALL abort
         ENDIF
         zovalues(ji) = zvalues( ipos )
      END DO

   END SUBROUTINE bufr_get_data

   SUBROUTINE spd_to_uv( pdir, pspd, pu, pv )
      !!----------------------------------------------------------------------
      !!   
      !!                    *** ROUTINE spd_to_uv ***
      !! 
      !! ** Purpose : Convert ocean current direction and speed to uv.
      !! 
      !! ** Method  : According to WMO current directors are coded according
      !!              to normal oceanographic convention.
      !!
      !! ** Action  :
      !!
      !! ** References:
      !!     http://www.wmo.int/pages/prog/www/WMOCodes/OperationalCodes.html
      !!
      !! History
      !!----------------------------------------------------------------------
      !! * Modules used
      !! * Arguments
      REAL, INTENT(IN) :: &
         & pdir, &       ! Current direction
         & pspd          ! Current speed
      REAL, INTENT(OUT) :: &
         & pu, &         ! Zonal corrent
         & pv            ! Meridonal current
      !! * Local variables
      REAL(KIND=dp), PARAMETER ::  rad = 3.141592653589793_dp / 180.0_dp
      REAL(KIND=dp) :: zdir
      
      zdir = pdir * rad

      pu = pspd * SIN( zdir )
      pv = pspd * COS( zdir )

   END SUBROUTINE spd_to_uv

#include "ctl_stop.h90"

#include "greg2jul.h90"

#endif

END MODULE bufrdata