!Subroutine for read input files:
!READ LAKES    - read the input data for lake areas
!READ_MAPS     - read the geomorphology input files
!READ_CLIMATE  - read the climate input files
!GET_RESTART   - read the restart files - it is not working yet

      SUBROUTINE READ_LAKES!(laket, alon, alat, outi, outj, &
!                          sillh, alev, ierr, times)
      
      use comvar, only: outi, &
                        outj, &
                        sillh, &
                        incf, &
                        inrf, &
                        alat, &
                        alon
      implicit none
      include 'netcdf.inc'

      INTEGER ierr
      REAL alev
      INTEGER, DIMENSION(12) :: times
      INTEGER, DIMENSION(4) :: nstart,icount
      CHARACTER*80 :: filen, name3d
      CHARACTER*21 :: aname
! dimensions for files
!
      nstart(1) = 1
      nstart(2) = 1
      nstart(3) = 1
      nstart(4) = 1
      icount(1) = incf
      icount(2) = inrf
      icount(3) = 1
      icount(4) = 1
!
      filen = 'input/outi.nc'
      aname = 'outi'
      name3d = 'level'
!
      CALL readvar(filen,aname,name3d,nstart,icount,outi, &
       alon,alat,alev,times,ierr)
!
      filen = 'input/outj.nc'
      aname = 'outj'
      name3d = 'level'
!
      CALL readvar(filen,aname,name3d,nstart,icount,outj, &
       alon,alat,alev,times,ierr)
!
      filen = 'input/mflac.reg.nc'
      aname = 'mflac'
      name3d = 'level'
!
      CALL readvar(filen,aname,name3d,nstart,icount,sillh, &
       alon,alat,alev,times,ierr)
!
      WRITE(*,*)
      WRITE(*,*)'**********************'
      WRITE(*,*)'ALL LAKES FILES READED'
      WRITE(*,*)'**********************'

      END SUBROUTINE READ_LAKES

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! READ GEOMORPHOLOGY FILES
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      SUBROUTINE READ_MAPS!(grid2, stddev, uparea, &
!                     outdir, sinuos, wetmask)
      
      use comvar, only: grid2, &
                        stddev, &
                        uparea, &
                        outdir, &
                        sinuos, &
                        wetmask, &
                        incf, &
                        inrf, &
                        alat, &
                        alon, &
                        iini, &
                        jini
                        
      
      implicit none
      include 'netcdf.inc'
!      include 'comvar.h'

      REAL, DIMENSION(incf,inrf) :: rivdir, outdir2
      INTEGER ierr, i, j, i2, j2
      REAL alev
      INTEGER, DIMENSION(12) :: times
      INTEGER, DIMENSION(4) :: nstart,icount
      CHARACTER*80 :: filen, name3d
      CHARACTER*21 :: aname
!--------------------------------------------------------------

! dimensions for geomorphology files
!
      nstart(1) = iini
      nstart(2) = jini
      nstart(3) = 1
      nstart(4) = 1
      icount(1) = incf
      icount(2) = inrf
!      icount(1) = 570
!      icount(2) = 498
      icount(3) = 1
      icount(4) = 1
      name3d = 'level'
      
!      print *, incf, inrf, alat, alon
!      print *, nstart, icount
! --- TOPOGRAPHY
      WRITE(*,*)'* ----------------------------*'
      WRITE(*,*)'* reading topography...       *'
      filen = 'input/NOAA3_mdt_global_10km.nc'   !SRTM DEM Hess and Eva
      aname = 'topo'
      CALL readvar(filen,aname,name3d,nstart,icount,grid2, &
       alon,alat,alev,times,ierr)

! --- STANDARD DEVIATION
      WRITE(*,*)'* -------------------------------*'
      WRITE(*,*)'* reading standard deviation...  *'
      filen = 'input/NOAA3_desvio_global_10km.nc'   !SRTM DEM Hess and Eva
      aname = 'desvio'
      CALL readvar(filen,aname,name3d,nstart,icount,stddev, &
       alon,alat,alev,times,ierr)

! --- RIVER DIRECTION
      WRITE(*,*)'* ------------------------------*'
      WRITE(*,*)'* reading river direction...    *'

      filen = 'input/NOAA3_ldd_thmb_global_10km.nc'   !SRTM DEM Hess and Eva
       aname = 'ldd'
      CALL readvar(filen,aname,name3d,nstart,icount,outdir, &
       alon,alat,alev,times,ierr)

!Testind Xingu correction - AAC - TO THE GLOBAL MAPS
!      DO j=1,inrf
!         DO i=1,incf
!            i2=iini+i-1
!            j2=jini+j-1
!            IF((i2.eq.1533).and.(j2.eq.1190)) print *, i2, j2, outdir(i,j)
!            IF((i2.eq.1533).and.(j2.eq.1191)) print *, i2, j2, outdir(i,j)
!            IF((i2.eq.1533).and.(j2.eq.1192)) print *, i2, j2, outdir(i,j)
!            IF((i2.eq.1533).and.(j2.eq.1193)) print *, i2, j2, outdir(i,j)
!            IF((i2.eq.1533).and.(j2.eq.1194)) print *, i2, j2, outdir(i,j)
!            IF((i2.eq.1533).and.(j2.eq.1195)) print *, i2, j2, outdir(i,j)
!            IF((i2.eq.1533).and.(j2.eq.1196)) print *, i2, j2, outdir(i,j)
!         ENDDO
!      ENDDO
! Xingu correction - AAC
!      DO j=1,inrf
!         DO i=1,incf
!            i2=iini+i-1
!            j2=jini+j-1
!            IF((i2.eq.1533).and.(j2.eq.1190)) outdir(i,j)=1
!            IF((i2.eq.1533).and.(j2.eq.1191)) outdir(i,j)=1
!            IF((i2.eq.1533).and.(j2.eq.1192)) outdir(i,j)=1
!            IF((i2.eq.1533).and.(j2.eq.1193)) outdir(i,j)=1
!            IF((i2.eq.1533).and.(j2.eq.1194)) outdir(i,j)=1
!            IF((i2.eq.1533).and.(j2.eq.1195)) outdir(i,j)=1
!            IF((i2.eq.1533).and.(j2.eq.1196)) outdir(i,j)=1
!         ENDDO
!      ENDDO
!END Xingu correction

!!Testind Xingu correction - AAC - TO THE REGIONAL MAPS (EGU inputs) - AM AND CERRADO
!      DO j=1,inrf
!         DO i=1,incf
!!            i2=iini+i-1
!!            j2=jini+j-1
!            IF((i.eq.357).and.(j.eq.182)) print *, i, j, outdir(i,j)
!            IF((i.eq.357).and.(j.eq.183)) print *, i, j, outdir(i,j)
!            IF((i.eq.357).and.(j.eq.184)) print *, i, j, outdir(i,j)
!            IF((i.eq.357).and.(j.eq.185)) print *, i, j, outdir(i,j)
!            IF((i.eq.357).and.(j.eq.186)) print *, i, j, outdir(i,j)
!            IF((i.eq.357).and.(j.eq.187)) print *, i, j, outdir(i,j)
!            IF((i.eq.357).and.(j.eq.188)) print *, i, j, outdir(i,j)
!         ENDDO
!      ENDDO
!! Xingu correction - AAC
!      DO j=1,inrf
!         DO i=1,incf
!            IF((i.eq.357).and.(j.eq.182)) outdir(i,j)=1
!            IF((i.eq.357).and.(j.eq.183)) outdir(i,j)=1
!            IF((i.eq.357).and.(j.eq.184)) outdir(i,j)=1
!            IF((i.eq.357).and.(j.eq.185)) outdir(i,j)=1
!            IF((i.eq.357).and.(j.eq.186)) outdir(i,j)=1
!            IF((i.eq.357).and.(j.eq.187)) outdir(i,j)=1
!            IF((i.eq.357).and.(j.eq.188)) outdir(i,j)=1
!         ENDDO
!      ENDDO
!!END Xingu correction

! --- SINUOSITY
      WRITE(*,*)'* ----------------------------*'
      WRITE(*,*)'* reading sinuosity...        *'

      filen = 'input/NOAA3_sin_global_10km.nc'
      aname = 'sin'
      CALL readvar(filen,aname,name3d,nstart,icount,sinuos, &
       alon,alat,alev,times,ierr)

! --- FLOODED AREA
      WRITE(*,*)'* ----------------------------*'
      WRITE(*,*)'* reading flooded area ...    *'
      filen = 'input/ALAGAVEL.nc'
      aname = 'ALAGAVEL'
      CALL readvar(filen,aname,name3d,nstart,icount,wetmask, &
       alon,alat,alev,times,ierr)

! --- UPSTREAM AREA
      WRITE(*,*)'* ----------------------------*'
      WRITE(*,*)'* reading upstream area...    *'
      filen = 'input/NOAA3_ca_global_10km.nc'
      aname = 'ca'
      CALL readvar(filen,aname,name3d,nstart,icount,uparea, &
       alon,alat,alev,times,ierr)

      DO i=1, incf
         DO j=1, inrf
            IF(wetmask(i,j).lt.0.) wetmask(i,j)=0.
!              print *, i, j, outdir(i,j), wetmask(i,j)
         ENDDO
      ENDDO

      WRITE(*,*)'* ----------------------------*'
      WRITE(*,*)
      WRITE(*,*)'ALL GEOMORPHOLOGY FILES READED'
      WRITE(*,*)
      WRITE(*,*)'* ----------------------------*'

      
      PRINT *, 'topo    i=100, j=100 ',   grid2(100,100)
      PRINT *, 'Stddev  i=100, j=100 ',  stddev(100,100)
      PRINT *, 'outdir  i=362, j=420 ',  outdir(100,100)
      PRINT *, 'sinuos  i=100, j=100 ',  sinuos(100,100)
      PRINT *, 'wetmask i=100, j=100 ', wetmask(100,100)
      PRINT *, 'uparea  i=100, j=100 ',  uparea(100,100)
      CALL SLEEP(5)
      
      RETURN

      END SUBROUTINE READ_MAPS


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! READ CLIMATE FILES
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      
      SUBROUTINE READ_CLIMATE

      use comvar, only: nmons2, nmons3, & !nmons
                        inc, &
                        inr, &
                        runin, &
                        prcpi, &
                        evapi, &
                        drainin
      implicit none
      include 'netcdf.inc'
      
      INTEGER ntime, nlevel, ierr, k, j, i
      REAL alev
      CHARACTER*80 :: filen,units, title
      CHARACTER*21 :: aname

!--------------------------------------------------------------
! READ netcdf monthly input of runoff and precipitation
!
! GLOBAL INPUT
! get runoff and drainage input (from IBIS), precip and evap 
! (from Coupled Atmospheric-Ocean Model-INPE)
!
      WRITE(*,*)'* ----------------------------*'
      WRITE(*,*)' total months to input climate data: ', nmons2
      WRITE(*,*)
      WRITE(*,*)' start to read climate data'
      WRITE(*,*)' inc:', inc,'inr:', inr,'nmons:', nmons2
      WRITE(*,*)
      WRITE(*,*)'* ----------------------------*'
      CALL SLEEP(5)


       nlevel = 1
! -------------------------------------------------------------
! --- RUNOFF
      WRITE(*,*)'* ----------------------------*'
      WRITE(*,*)'* reading runoff ...          *'
 
      filen = 'input/RSFC_199001_201012_v2_AS.nc' !correction to IBIS output
      aname = 'rsfc$'
      DO ntime  = 1,nmons2
         CALL rd4dnc (runin(1,1,ntime), inc*inr, inc, inr, nlevel, ntime, &
          aname, units, title, filen, ierr )
      ENDDO
      WRITE(*,*)
      WRITE(*,*)'* runoff readed!              *'

! --- DRAINAGE
      WRITE(*,*)'* ----------------------------*'
      WRITE(*,*)'* reading drainage ...        *'
      filen = 'input/RNSG_199001_201012_v2_AS.nc' !correction to IBIS output
      aname = 'rnsg$'
      DO ntime  = 1,nmons2
         CALL rd4dnc (drainin(1,1,ntime), inc*inr, inc, inr, nlevel, ntime, &
          aname, units, title, filen, ierr )
      ENDDO
      WRITE(*,*)
      WRITE(*,*)'* drainage readed!            *'

! --- PRECIPITATION
      WRITE(*,*)'* ----------------------------*'
      WRITE(*,*)'* reading precipitation ...   *'
      filen = 'input/PREC_199001_201012_v2_AS.nc' !correction to IBIS output
      aname = 'prec$'
      DO ntime  = 1,nmons2
         CALL rd4dnc (prcpi(1,1,ntime), inc*inr, inc, inr, nlevel, ntime, &
          aname, units, title, filen, ierr )
      ENDDO
      WRITE(*,*)
      WRITE(*,*)'* precipitation  readed!      *'

! --- Evaporation
      WRITE(*,*)'* ----------------------------*'
      WRITE(*,*)'* reading evaporation ...     *'

       filen = 'input/EVTP_199001_201012_v2_AS.nc' !correction to IBIS output
       aname = 'evtp$'
       DO ntime  = 1,nmons2
!          ntimee = ntime + ((1915)-1901)*12 !adjust for evap file starting in 1901
          CALL rd4dnc (evapi(1,1,ntime), inc*inr, inc, inr, nlevel, ntime, &
          aname, units, title, filen, ierr )
       ENDDO
      WRITE(*,*)
      WRITE(*,*)'* evaporation  readed!        *'

!
!--------------------------------------------------------------
! convert input from mm/day to m/s
!

      DO k = 1,nmons2
      !DO k = 1,24
         DO 335 j = 1,inr
            DO 336 i = 1,inc
!
               IF (  runin(i,j,k) .ge. 1.e+20)  runin(i,j,k) = 0.
               IF (drainin(i,j,k) .ge. 1.e+20)drainin(i,j,k) = 0.
               IF (  prcpi(i,j,k) .ge. 1.e+20)  prcpi(i,j,k) = 0.
               IF (  evapi(i,j,k) .ge. 1.e+20)  evapi(i,j,k) = 0.
!
! Daily input
               runin(i,j,k)   = max(   (runin(i,j,k))/0.864e+8,0.)      !IBIS data
               drainin(i,j,k) = max( (drainin(i,j,k))/0.864e+8,0.)      !IBIS data
               prcpi(i,j,k)   = max(   (prcpi(i,j,k))/0.864e+8,0.)      !IBIS data
               evapi(i,j,k)   = max((-1*evapi(i,j,k))/0.864e+8,0.)      !RAcosta ERA5 data
!              print *, i, j, evapi(31,37,k)
! Hourly input
!               runin(i,j,k)   = max((runin(i,j,k))/0.36e+7,0.)      !IBIS data
!               drainin(i,j,k) = max((drainin(i,j,k))/0.36e+7,0.)    !IBIS data
!               prcpi(i,j,k)   = max((prcpi(i,j,k))/0.36e+7,0.)      !IBIS data
!               evapi(i,j,k)   = max((evapi(i,j,k))/0.36e+7,0.)      !IBIS data
!
 336        CONTINUE
 335     CONTINUE
      ENDDO

      PRINT *,'  runoff i=100, j=100, t = 4 ',   runin(100,100,4)*0.864e+8
      PRINT *,'drainage i=100, j=100, t = 4 ', drainin(100,100,4)*0.864e+8
      PRINT *,'  precip i=100, j=100, t = 4 ',   prcpi(100,100,4)*0.864e+8
      PRINT *,'   evapi i=100, j=100, t = 4 ',   evapi(100,100,4)*0.864e+8
      CALL SLEEP(5)

      WRITE(*,*) 
      WRITE(*,*)'-------------------------------'
      WRITE(*,*)      
      WRITE(*,*)'  ALL CLIMATE FILES READED     '
      WRITE(*,*)      
      WRITE(*,*)'-------------------------------'

      


      END SUBROUTINE READ_CLIMATE

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! GETTING RESTART
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!      SUBROUTINE GET_RESTART(iwmon, incf, inrf, along, alatg, dvoll, height, sflux, nrf, ncf, &
!                          alat, alon, alev, times, sfluxin, voll, warea, jstart, istart)

      SUBROUTINE GET_RESTART(iwmon, incf, inrf, along, alatg, height, sflux, sfluxin, fluxout, &
                            nrf, ncf, alat, alon, alev, times, &
                            voll, dvoll, warea, jini, iini)
      
      include 'netcdf.inc'

      INTEGER iwmon, incf, inrf, nlevel, ntime, jj, ii, i, j, istart
      INTEGER iend, jstart, jend, ierr, ncf, nrf
      INTEGER, DIMENSION(4) :: nstart,icount
      INTEGER, DIMENSION(12) :: times
      REAL alev, atim
      REAL, DIMENSION(nrf) :: alat
      REAL, DIMENSION(ncf) :: alon
      CHARACTER*80 :: filen,units, title,uidim,ujdim,ukdim,uldim, name3d
      CHARACTER*21 :: aname
      REAL, DIMENSION(inrf) :: alatg 
      REAL, DIMENSION(incf) :: along 
      REAL, DIMENSION(incf,inrf) :: height, sflux, sfluxin, fluxout 
      REAL, DIMENSION(incf,inrf) :: dummy, warea
      REAL(KIND=8), DIMENSION(incf,inrf) :: voll, dvoll


      WRITE(*,*)
           WRITE(*,*) ': --- RESTART INDEX iwmon/incf/inrf/nrf/ncf/jini/iini: '
           WRITE(*,*) ': --- ', iwmon,incf,inrf,nrf,ncf,jini,iini

            WRITE(*,*) ': --- RESTART INDEX JSTART/JEND & ISTART/IEND: '
           WRITE(*,*) ': --- ', jstart, jend, istart, iend
           WRITE(*,*)



!--------------------------------------------------------------
! Get volume and lake area files if this is a Restart
!
      nstart(1) = 1
      nstart(2) = 1
      nstart(3) = 1
      nstart(4) = 1
      icount(1) = incf
      icount(2) = inrf
      icount(3) = 1
      icount(4) = 1

  
! --- HEIGHT
      WRITE(*,*)
      WRITE(*,*)
      WRITE(*,*) '---   RESTART HEIGHT OK!  ---- '

      filen =  'restart/height.nc'
      aname = 'height'
      name3d = 'level'
      CALL readvar(filen,aname,name3d,nstart,icount,height, &
       alon,alat,alev,times,ierr)

! --- SFLUX
      WRITE(*,*)
      WRITE(*,*)
      WRITE(*,*) '---   RESTART SFLUX OK!  ---- '

      filen =  'restart/sflux.nc'
      aname = 'sflux'
      name3d = 'level'
      CALL readvar(filen,aname,name3d,nstart,icount,sflux, &
       alon,alat,alev,times,ierr)

! --- SFLUXIN
      WRITE(*,*)
      WRITE(*,*)
      WRITE(*,*) '---   RESTART SFLUXIN OK!  ---- '

      filen =  'restart/sfluxin.nc'
      aname = 'sfluxin'
      name3d = 'level'
      CALL readvar(filen,aname,name3d,nstart,icount,sfluxin, &
       alon,alat,alev,times,ierr)

! --- FLUXOUT
      WRITE(*,*)
      WRITE(*,*)
      WRITE(*,*) '---   RESTART FLUXOUT OK!  ---- '

      filen =  'restart/fluxout.nc'
      aname = 'fluxout'
      name3d = 'level'
      CALL readvar(filen,aname,name3d,nstart,icount,fluxout, &
       alon,alat,alev,times,ierr)      

! --- VOLL
      WRITE(*,*)
      WRITE(*,*)
      WRITE(*,*) '---   RESTART VOLL OK!  ---- '

      filen =  'restart/voll.nc'
      aname = 'voll'
      name3d = 'level'
      CALL readvar(filen,aname,name3d,nstart,icount,dummy, &
       alon,alat,alev,times,ierr)

      DO j = 1,inrf
        DO i = 1,incf
           voll(i,j) = dble(dummy(i,j))
        ENDDO
      ENDDO

! --- DVOLL
      WRITE(*,*)
      WRITE(*,*)
      WRITE(*,*) '---   RESTART DVOLL OK!  ---- '

      filen =  'restart/dvoll.nc'
      aname = 'dvoll'
      name3d = 'level'
      CALL readvar(filen,aname,name3d,nstart,icount,dummy, &
       alon,alat,alev,times,ierr)
      DO j = 1,inrf
        DO i = 1,incf
           dvoll(i,j) = dble(dummy(i,j))
        ENDDO
      ENDDO

! --- WAREA
      WRITE(*,*)
      WRITE(*,*)
      WRITE(*,*) '---   RESTART WAREA OK!  ---- '

      filen =  'restart/warea.nc'
      aname = 'warea'
      name3d = 'level'
      CALL readvar(filen,aname,name3d,nstart,icount,warea, &
       alon,alat,alev,times,ierr)

      END SUBROUTINE GET_RESTART

