subroutine read_goesglm(nread,ndata,nodata,infile,obstype,lunout,twindin,sis)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:  read_goesglm                reads lightning obs from a BUFR file
!   prgmmr:      k apodaca <karina.apodaca@colostate.edu>
!      org:      CSU/CIRA, Data Assimilation Group
!     date:      2015-03-12
!
! abstract:  This routine reads lightning data (Earth-relative location, frequency) from 
!            file and prepares it for assimiliation in the form of lightning flash rate 
!            (#hits km-2 hr-1).  
!
!            Note: when running GSI in regional mode, the code only
!            retains those observations that fall within the regional
!            domain
!
! program history log:
!   2015-06-12  zupanski  - include the convert_to_flash_rate subroutine to convert 
!                           lightning strike observations into lightning flash rate 
!                           (#hits km-2 hr-1). 
!   2015-11-18  apodaca   - include the convert_time subroutine to deal with computer 
!                           precision dependencies.
!   2018-02-07  apodaca   - add further documentation
!
!   input argument list:
!     infile   - unit from which to read BUFR data
!     obstype  - observation type to process
!     lunout   - unit to which to write data for further processing
!
!   output argument list:
!     nread    - number of type "obstype" observations read
!     nodata   - number of individual "obstype" observations read
!     ndata    - number of type "obstype" observations retained for further processing
!     twindin  - input group time window (hours)
!     sis      - satellite/instrument/sensor indicator
!
! attributes:
!   language: Fortran 90 and/or above
!   machine:  
!
!$$$
  use kinds, only: r_single,r_kind,r_double,i_kind
  use constants, only: zero,one_tenth,one,deg2rad,&
      three,rad2deg,&
      r60inv,ten
  use gridmod, only: diagnostic_reg,wrf_mass_regional,regional,nlon,nlat,&
      tll2xy,txy2ll,&
      rlats,rlons
  use lightinfo, only: iuse_light,nlighttype
  use obsmod, only: iadate
  use obsmod, only: offtime_data
  use gsi_4dvar, only: l4dvar,l4densvar,time_4dvar,winlen

  implicit none

! Declare passed variables
  character(len=*)                      ,intent(in   ) :: infile,obstype
  character(len=*)                      ,intent(in   ) :: sis
  integer(i_kind)                       ,intent(in   ) :: lunout
  integer(i_kind)                       ,intent(inout) :: nread,ndata
  real(r_kind)                          ,intent(in   ) :: twindin
  integer(i_kind)                                      :: nodata

! Declare local parameters
  real(r_kind),parameter:: r6   = 6.0_r_kind
  real(r_kind),parameter:: r90  = 90.0_r_kind
  real(r_kind),parameter:: r180 = 180.0_r_kind
  real(r_kind),parameter:: r360 = 360.0_r_kind

!--- Declare local variables
  logical lob
  logical outside

  character(40) hdstr,oestr,qcstr
  character(10) date
  character(8)  subset
  character(1)  sidchr(8)

  integer(i_kind) ireadmg,ireadsb,icntpnt,icount
  integer(i_kind) lunin,i
  integer(i_kind) itx
  integer(i_kind) ihh,idd,idate,iret,im,iy,k
  integer(i_kind) nchanl,nreal,ilat,ilon
  integer(i_kind) lqm
  integer(i_kind) iout
  integer(i_kind) ntest,nvtest
  integer(i_kind) minobs,minan
  integer(i_kind) ntb
  integer(i_kind) nmsg                ! message index
  integer(i_kind),parameter :: maxobs=2000000
  integer(i_kind),dimension(5):: idate5
  integer(i_kind),allocatable,dimension(:):: isort,iloc

  real(r_kind) time
  real(r_kind) usage
  real(r_kind) loe,lmerr
  real(r_kind) time_correction
  real(r_kind) dlat,dlon,dlat_earth,dlon_earth
  real(r_kind) cdist,disterr,disterrmax,rlon00,rlat00
  real(r_kind) vdisterrmax
  real(r_kind) timex,timeobs,toff,t4dv,zeps
  real(r_kind),allocatable,dimension(:,:):: cdata_all,cdata_out
!--- flash rate
  real(r_kind),allocatable,dimension(:,:):: cdata_flash,cdata_flash_h
  integer(i_kind)                        :: ndata_flash,ndata_flash_h

  real(r_double) rstation_id

  real(r_double),dimension(3):: hdr
  real(r_double),dimension(1,1):: qcmark,obserr

!  equivalence to handle character names
  equivalence(rstation_id,sidchr)

!--- data statements
  data hdstr  /'XOB YOB DHR'/
  data oestr  /'LOE'/ 
  data qcstr  /'LQM'/   

  data lunin / 13 /


  nreal=13
  lob = obstype == 'goes_glm'

!                .      .    .                                       .

! Open, then read date from BUFR file


  call closbf(lunin) 
  open(lunin,file=infile,form='unformatted') 
  call openbf(lunin,'IN',lunin) 
  call datelen(10) 


! Initialization

  ntb = 0 
  nmsg = 0 
  disterrmax=-9999.0_r_kind 

  allocate(cdata_all(nreal,maxobs),isort(maxobs))
  isort = 0
  cdata_all=zero
  nread=0
  ntest=0
  nvtest=0
  nchanl=0
  ilon=2
  ilat=3
  icntpnt=0

!                .      .    .                                       .
! Big loop over glmbufr file : READING THE BUFR FILE
     
  loop_msg: do while (ireadmg(lunin,subset,idate)== 0) 
     nmsg = nmsg+1 
     write(*,'(3a,i10)') 'subset=',subset,' cycle time =',idate 


     loop_readsb: do while(ireadsb(lunin) == 0)  
!       use msg lookup table to decide which messages to skip
!       use report id lookup table to only process matching reports
        ntb = ntb+1 

          
!       Extract location and date information
        call ufbint(lunin,hdr,3,1,iret,hdstr)

        if(abs(hdr(2))>r90 .or. abs(hdr(1))>r360) cycle loop_readsb
        if(hdr(1) > r180)hdr(1)=hdr(1)-r360
        if(hdr(1) < zero)hdr(1)=hdr(1)+r360
        dlon_earth=hdr(1)*deg2rad
        dlat_earth=hdr(2)*deg2rad
                
        if (regional) then

!-- WRF-ARW

           if (wrf_mass_regional) then

              call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside)    ! convert to rotated coordinates
              if(diagnostic_reg) then
                 call txy2ll(dlon,dlat,rlon00,rlat00)
                 ntest=ntest+1
                 cdist=sin(dlat_earth)*sin(rlat00)+cos(dlat_earth)*cos(rlat00)* &
                       (sin(dlon_earth)*sin(rlon00)+cos(dlon_earth)*cos(rlon00))
                 cdist=max(-one,min(cdist,one))
                 disterr=acos(cdist)*rad2deg
                 disterrmax=max(disterrmax,disterr)
              end if
              if(outside) cycle loop_readsb   ! check to see if outside regional domain

           endif ! wrf_mass_regional

        endif !if (regional) then

! Global

        if (.not. regional) then
            dlat = dlat_earth
            dlon = dlon_earth
            call grdcrd1(dlat,rlats,nlat,1)
            call grdcrd1(dlon,rlons,nlon,1)
        endif !  end global block 

        if (offtime_data) then
 
!           in time correction for observations to account for analysis
!           time being different from obs file time.
            write(date,'( i10)') idate
            read (date,'(i4,3i2)') iy,im,idd,ihh
            idate5(1)=iy
            idate5(2)=im
            idate5(3)=idd
            idate5(4)=ihh
            idate5(5)=0
            call w3fs21(idate5,minobs)    !  obs ref time in seconds relative to historic date
            idate5(1)=iadate(1)
            idate5(2)=iadate(2)
            idate5(3)=iadate(3)
            idate5(4)=iadate(4)
            idate5(5)=0
            call w3fs21(idate5,minan)    !  analysis ref time in seconds relative to historic date
           
!           Add obs reference time, then subtract analysis time to get obs time relative to analysis
            time_correction=float(minobs-minan)*r60inv
        else
            time_correction=zero
        end if

        timeobs=real(real(hdr(3),r_single),r_double)
        t4dv=timeobs + toff
        zeps=1.0e-8_r_kind
        if (t4dv<zero  .and.t4dv>      -zeps) t4dv=zero
        if (t4dv>winlen.and.t4dv<winlen+zeps) t4dv=winlen
        t4dv=t4dv + time_correction
        time=timeobs + time_correction


        if (l4dvar.or.l4densvar) then
           if (t4dv<zero.or.t4dv>winlen) cycle loop_readsb ! outside time window
        else
           if((real(abs(time)) > real(twindin)))cycle loop_readsb ! outside time window
        endif

        timex=time


!       Note: An assessment of the GLM detection error is still a work in
!       progress and at present, there are no clear metrics to assign a 
!       "quality mark" for an effective quality control procedure for
!       lightning observations from this sensor. Nonetheless, the infrastructure
!       for passing observations errors and quality control information to other 
!       routines within the GSI source code has been left in place. In addition, a 
!       temporarty undefined value (-.9999) has been assigned to the mnemonics
!       corresponding to these variables in the lightning BUFR file. Developments 
!       for the calculation of observation errors and quality control (sanity checks) are 
!       expected in future upgrades to the "GOES/GLM variational lightning 
!       assimilation package."

!       Extract observation error informatiom
        call ufbint(lunin,obserr,1,1,iret,oestr)

        loe=obserr(1,1)
        lmerr=loe
 
!       Extract quality control information
        call ufbint(lunin,qcmark,1,1,iret,qcstr)
 
        lqm=qcmark(1,1)

!       Data counter

        nread=nread+1
        icntpnt=icntpnt+1

        ndata=ndata+1
        nodata=nodata+1
        iout=ndata
        isort(icntpnt)=iout

        if (ndata > maxobs) then
           write(6,*)'READ_GOESGLM:  ***WARNING*** ndata > maxobs for ',obstype
           ndata = maxobs
        end if

             
!       Set usage variable              
        usage = zero

        if (iuse_light(nlighttype) <= 0)usage=100._r_kind
        if (lob) then
           cdata_all(1,iout) =loe                  ! lightning observation error
           cdata_all(2,iout) =dlon                 ! grid relative longitude
           cdata_all(3,iout) =dlat                 ! grid relative latitude
           cdata_all(4,iout) =iout                 ! lightning obs
           cdata_all(5,iout) =rstation_id          ! station id
           cdata_all(6,iout) =t4dv                 ! analysis time
           cdata_all(7,iout) =nlighttype           ! type
           cdata_all(8,iout) =lmerr                ! lightning max error
           cdata_all(9,iout) =lqm                  ! quality mark
           cdata_all(10,iout)=loe                  ! original lightning obs error loe 
           cdata_all(11,iout)=usage                ! usage parameter
           cdata_all(12,iout)=dlon_earth*rad2deg   ! earth relative lon (degrees)
           cdata_all(13,iout)=dlat_earth*rad2deg   ! earth relative lat (degrees)
        end if


! end loop on read line BUFR

     end do loop_readsb

! end of BUFR read loop

  enddo loop_msg   !Uncomment if reading messeges in a loop

!                .      .    .                                       .

! Close unit to bufr file

  call closbf(lunin)

! Write header record and data to output file for further processing
  allocate(iloc(ndata))
  icount=0.
  do i=1,maxobs
     if(isort(i) > 0)then
       icount=icount+1
       iloc(icount)=isort(i)
     end if
  end do
  if(ndata /= icount)then
     write(6,*) ' READ_GOESGLM: mix up in read_goesglm ,ndata,icount ',ndata,icount
     call stop2(50)
  end if
 
  allocate(cdata_out(nreal,ndata))
  do i=1,ndata
     itx=iloc(i)
     do k=1,nreal
        cdata_out(k,i)=cdata_all(k,itx)
     end do
  end do

  deallocate(iloc,isort,cdata_all)

!                .      .    .                                       .

! Call to the subroutine that transforms lightning strikes into lightning flash rate


  if(ndata /= 0) then

      !! count flash rate data and alocate temporary domain 
      !! begin with the current number of strikes as the theoretical upper limit
   
     ndata_flash_h=ndata
  
     allocate(cdata_flash_h(nreal,ndata_flash_h))

     call convert_to_flash_rate   &
              (nreal,ndata,cdata_out,ndata_flash_h,cdata_flash_h,ndata_flash)
 
     deallocate(cdata_out)
     ndata=ndata_flash
     allocate(cdata_flash(nreal,ndata))

     do i=1,ndata
        do k=1,nreal
           cdata_flash(k,i)=cdata_flash_h(k,i)
        end do
     end do

     deallocate(cdata_flash_h)

     write(lunout) obstype,sis,nreal,nchanl,ilat,ilon
     write(lunout) cdata_flash

!!!!!!!!!
!!! Write lat, lon, time, lightning flash rate "superobs" into a file: 
!!! cdata_flash_h(2,iout),cdata_flash_h(3,iout),cdata_flash_h(4,iout), &
!!! cdata_flash_h(6,iout)cdata_flash(4,:)

     deallocate(cdata_flash)

  else  ! ndata=0

     write(lunout) obstype,sis,nreal,nchanl,ilat,ilon
     write(lunout) cdata_out
     deallocate(cdata_out)

  end if  !!  if(ndata =/ 0) then


900 continue
  if(diagnostic_reg .and. ntest>0) write(6,*)'READ_GOESGLM:  ',&
     'ntest,disterrmax=',ntest,disterrmax
  if(diagnostic_reg .and. nvtest>0) write(6,*)'READ_GOESGLM:  ',&
     'nvtest,vdisterrmax=',ntest,vdisterrmax

  if (ndata == 0) then 
     call closbf(lunin)
     write(6,*)'READ_GOESGLM:  closbf(',lunin,')'
  endif


  close(lunin)

  close(55)

! End of routine
  return

end subroutine read_goesglm


!!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

subroutine convert_to_flash_rate   &
   (nreal,ndata_strike,cdata_strike,ndata_flash_h,cdata_flash_h,ndata_flash)

!$$$  documentation block
!                .      .    .                                       .
! subroutine:  convert_to_flash_rate      converts geo-located lightning strikes into
!                                         lightning flash rate (#hits km-2 hr-1)

!   prgmmr: zupanski <milija.zupanski@colostate.edu>
!      org: CSU/CIRA, Data Assimilation group
!     date: 2015-06-12
!
! abstract:  This subroutine does the following: 
!----
!- 1- counts the number of hits surrounding a GSI analysis grid point
!- 2- calculates the flash rate averaged over time and area (hr*km**2)
!- 3- finds the center of mass in terms of lon,lat and glon,glat
!- 4- assigns a center of mass to be the flash rate observation point
!---- 
!
! program history log:
!   2015-07-01  apodaca  - several updates in the calculation of lightning flashrate

  use kinds, only: r_single,r_kind,r_double,i_kind
  use constants, only: zero,one_tenth,one,deg2rad,&
      three,half,zero,&
      r10,r100,ten, r1000, rearth
  use gridmod, only: wrf_mass_regional,regional,nlon,nlat,&
      tll2xy,txy2ll
  use gridmod, only: lat2, lon2
  use wrf_mass_guess_mod, only: ges_xlon, ges_xlat
  use gsi_4dvar, only: nhr_assimilation

  implicit none

  integer(i_kind),intent(in   )               :: nreal,ndata_strike,ndata_flash_h
  integer(i_kind),intent(inout)               :: ndata_flash
  real(r_kind),intent(inout),dimension(nreal,ndata_strike)  :: cdata_strike
  real(r_kind),intent(inout),dimension(nreal,ndata_flash_h) :: cdata_flash_h

  real(r_kind),allocatable,dimension(:) :: gtim_central
  real(r_kind),allocatable,dimension(:) :: glon_central
  real(r_kind),allocatable,dimension(:) :: glat_central
  real(r_kind),allocatable,dimension(:) ::  lon_central
  real(r_kind),allocatable,dimension(:) ::  lat_central
  integer(i_kind),allocatable,dimension(:) ::  ind_min
  integer(i_kind),allocatable,dimension(:) ::  lcount
  
  real(r_kind)      :: rearth2
  real(r_kind)      :: dtime,darea,cosine
  real(r_kind)      :: darea_sum
  real(r_kind)      :: delta_lon,delta_lat
  real(r_kind)      :: lat_ref
  real(r_kind)      :: xx,yy
  real(r_kind)      :: dist2,dist_min
  integer(i_kind)   :: ii0,jj0
  integer(i_kind)   :: ngridh
  integer(i_kind)   :: index
  integer(i_kind)   :: iobs,usage 
  logical           :: xflag,yflag

  real(r_kind)      :: xbound,ybound
  integer(i_kind)   :: nxdim,nydim
  integer(i_kind)   :: icount
  
  !  Output files
!----


!- unit grid lat-lon distances (Earth difference divided by grid difference)
!- calculated from all averaged grid areas
!!- note: this can be relaxed if the unit grid area (km) is known (i.e. darea)
!! darea = (r*cos(lat)*dlon)*(dlat)


  rearth2=(rearth/r1000)**2  !! squared earth radius in km (need for darea calculation)

  if (ndata_strike>0) then
 
     darea_sum=0._r_kind
     do iobs=1,ndata_strike

        ii0=INT(cdata_strike(2,iobs))
        jj0=INT(cdata_strike(3,iobs))

        delta_lon=ges_xlon(jj0,ii0+1,1)-ges_xlon(jj0,ii0,1)
        delta_lat=ges_xlat(jj0+1,ii0,1)-ges_xlat(jj0,ii0,1)
  
        lat_ref  =half*(cdata_strike(13,iobs)+cdata_strike(13,iobs-1))

        cosine=cos(lat_ref*deg2rad)

        darea_sum=darea_sum+rearth2*cos(lat_ref*deg2rad)*(abs(delta_lon)*deg2rad)*&
                  (abs(delta_lat)*deg2rad)

     end do  !! do iobs=2,ndata_strike

     darea=darea_sum/float(ndata_strike)
  else   !! ndata_strike=0
     darea=zero

  end if  !! if(ndata_strike>0) then

  dtime=float(nhr_assimilation)

  ! Regional

  if (regional) then

!-- WRF-ARW

     if (wrf_mass_regional) then
      
        nxdim=lon2
        nydim=lat2

     endif ! wrf_mass_regional

  endif !if (regional) then

! Global

  if (.not. regional) then

     nxdim=nlon
     nydim=nlat

  endif !  end global block 

!!! Allocate new var for flash rate here (nxdim,nydim)
! asign zero to all points
! update the relevant points in the loop

  ngridh=nxdim*nydim


  allocate(gtim_central(1:ngridh))
  allocate(glon_central(1:ngridh))
  allocate(glat_central(1:ngridh))
  allocate( lon_central(1:ngridh))
  allocate( lat_central(1:ngridh))
  allocate(      lcount(1:ngridh))

  lcount(:)=0
  glon_central(:)=zero
  glat_central(:)=zero
  lon_central(:) =zero
  lat_central(:) =zero
  gtim_central(:)=zero
     
  do iobs=1,ndata_strike

     xx=cdata_strike(2,iobs)   !! glon
     yy=cdata_strike(3,iobs)   !! glat
     ii0=INT(cdata_strike(2,iobs))
     jj0=INT(cdata_strike(3,iobs))

!!   find lightning strikes near the (ii0,jj0) point

     xbound=float(ii0)
     ybound=float(jj0)

     xflag=(xx>xbound) .AND. (xx<xbound+1.)
     yflag=(yy>ybound) .AND. (yy<ybound+1.)

     if (xflag .AND. yflag ) then

        index=(jj0-1)*nxdim+ii0
        lcount(index)=lcount(index)+1
       
        glon_central(index)=glon_central(index)+cdata_strike(2,iobs)
        glat_central(index)=glat_central(index)+cdata_strike(3,iobs)
        lon_central(index)= lon_central(index)+cdata_strike(12,iobs) 
        lat_central(index)= lat_central(index)+cdata_strike(13,iobs)

        if (lcount(index)<2) then
            gtim_central(index)=cdata_strike(6,iobs)
        else
            call convert_time (gtim_central(index),cdata_strike(6,iobs),lcount(index))
        end if 

     end if  !! if(xflag .AND. yflag ) 
      
  enddo !! do iobs=1,ndata_strike



!-- find the center of mass

  do index=1,ngridh
     if (lcount(index)>0) then
        glon_central(index)=glon_central(index)/float(lcount(index))
        glat_central(index)=glat_central(index)/float(lcount(index))
        lon_central(index)= lon_central(index)/float(lcount(index))
        lat_central(index)= lat_central(index)/float(lcount(index))
     endif  !! if(lcount(index)>0) then
  enddo  !! do index=1,ngridh

!-- find the original index of the nearest strike (need for transfer of input obs)

  allocate(ind_min(1:ngridh))
  ind_min(:)=-99

  dist_min=1.e10_r_kind
  do iobs=1,ndata_strike

     xx=cdata_strike(2,iobs)   !! glon
     yy=cdata_strike(3,iobs)   !! glat
     ii0=INT(cdata_strike(2,iobs))
     jj0=INT(cdata_strike(3,iobs))
     index=(jj0-1)*nxdim+ii0

     if (lcount(index)>0) then
        dist2=(xx-glon_central(index))**2+(yy-glat_central(index))**2
        if (dist2<dist_min) then
           dist_min=dist2
           ind_min(index)=iobs
        end if
     end if  !! if(lcount(index)>0) then

  enddo  !! do iobs=1,ndata_strike

!----
!---- Output
!----

!-- count the non-zero flash rates and assign a temporary domain cdata_flash_h
!-- Note: it is assumed that only non-zero flash rates are true observations

  icount=0
  do index=1,ngridh
     if (lcount(index)>0) then
        icount=icount+1
        cdata_flash_h( 1,icount)=cdata_strike( 1,ind_min(index))
        cdata_flash_h( 2,icount)=glon_central(index)
        cdata_flash_h( 3,icount)=glat_central(index)

        if (darea>0._r_kind) then
           cdata_flash_h( 4,icount)=float(lcount(index))/(darea*dtime)
        else 
           cdata_flash_h( 4,icount)=0. 
        end if

        cdata_flash_h( 5,icount)=cdata_strike( 5,ind_min(index))
        cdata_flash_h( 6,icount)=gtim_central(index)
        cdata_flash_h( 7,icount)=cdata_strike( 7,ind_min(index))
        cdata_flash_h( 8,icount)=cdata_strike( 8,ind_min(index))
        cdata_flash_h( 9,icount)=cdata_strike( 9,ind_min(index))
        cdata_flash_h(10,icount)=cdata_strike(10,ind_min(index))
        cdata_flash_h(11,icount)=usage                 
        cdata_flash_h(12,icount)=lon_central(index)
        cdata_flash_h(13,icount)=lat_central(index)

     endif  !! if(lcount(index)>0) then

  enddo  !! do index=1,ngridh

  ndata_flash=icount

  deallocate(ind_min)
  deallocate(lcount)
  deallocate(gtim_central)
  deallocate(glon_central) 
  deallocate(glat_central)
  deallocate(lon_central) 
  deallocate(lat_central)

!-----
! End of routine
  return

end subroutine convert_to_flash_rate

!-----
subroutine convert_time (date_old,date_new,nmax)

!$$$  documentation block
!                .      .    .                                       .
! subroutine:  convert_time      
! prgmmr: k apodaca               date: 2015-11-18 
!
! abstract:  This subroutine performs a date conversion to deal with
!            computer precission dependencies associated with a 10-digit 
!            float for the analysis date/time.
!-- 
  use kinds, only: r_kind,i_kind

  implicit none

  integer(i_kind), intent(in) :: nmax
  real(r_kind), intent(inout) :: date_old
  real(r_kind), intent(in)  :: date_new
  integer(i_kind) :: i,sumidd
  integer(i_kind) :: idd,jdd,kdd
  real(r_kind), allocatable :: xdate(:) 
  real(r_kind) :: dd,hh,ysumidd,xsumidd
  real(r_kind) :: xdd,xhh,ydate
  real(r_kind) :: xccyy

  allocate(xdate(1:nmax))

  xdate(1:nmax-1) = date_old
  xdate(nmax) = date_new
    
  sumidd=0._r_kind
  do i=1,nmax
     xccyy = INT(1.0e-8_r_kind*xdate(i))*1.0e8_r_kind 
     xdate(i) = INT(xdate(i))-xccyy   

     jdd=INT(0.0001_r_kind*xdate(i))
     idd=INT(xdate(i))-jdd*10000

     ysumidd=float(idd)
     dd=float(INT(0.01_r_kind*ysumidd))
     hh=ysumidd-dd*100._r_kind

     sumidd=sumidd+dd*24._r_kind+hh

  enddo  !! do i=1,nmax

  xsumidd=float(sumidd)/nmax
  ysumidd=float(INT(xsumidd))

  kdd=INT(xsumidd/24._r_kind)
  xdd=float(kdd)
  xhh=ysumidd-float(kdd)*24._r_kind

  ydate=float(jdd)*10000._r_kind+xdd*100._r_kind+xhh+xccyy
      
  date_old=ydate 
      
  deallocate(xdate)

end subroutine convert_time

