      PROGRAM THMB
!
! ------------------------  06/2003 ---------------------------
!
! This is the beginning of the LBA-ECO phase II model (funding period 
! 2003-2006). It starts with hydra.v2.f that was used in the original
! LBA-Ecology Phase I projects (Coe et al., Foley et al., 2002).
!
! v2 HYDRA differs from hydra.f in one way.
! The area and height of the floodwaters are calculated. They are
! calculated by applying any water above a pre-determined
! bankfull volume to the gridcell. Then shuttling the floodwaters
! across the landscape in a manner identical to the river transport 
! except that the flow direction is dynamic rather than prescribed.
! The floodwaters are calculated entirely seperately from the river
! waters. This guarantees that the river calculation remains stable
! but does not allow for the calculation of the impact of flood
! waters on the river discharge timing. This very simple method
! is comparable to a combination of Bates and de Roo, J Hydrol. 236, 54-77 
! and Vorosmarty 1989.
!
! ************************************************************
! hydra.v2.1 -- includes river sinuosity in distance calculation
!
! Simply READ in sinusoity data (from Marcos) and multiply
! the geometric cell length by this dimensionless value to get
! the stream length in each grid cell.

! hydra.v2.2 -- has all i/o in netcdf no hdf at all. 11/2003
!
! hydra.v2.3 -- has P-E over flooded areas           11/2003
!
!	For P-E over flooded areas, I substituted warea for larea
! in the calculation of the local water budget. Therefore, the 
! the local water budget is warea*(P-E) + 1-warea*(trunoff)
!
! hydra.v2.4 -- has flood initiation height changes  11/2003
!
!	The flood initiation has been altered to be a function
! of the water height and river width. We used an empirical
! function (derived by Marcos from the river data) to indicate
! at what river height the river leaves its banks. Additionally,
! we can derive the river width that cooresponds to that height.
! Then the program simply asks if the local volume is > than
! the product of the width, height, and river length. If so,
! initiate flooding as in previous version.
!
! hydra.v2.5 -- has alteration to the flood plain algorithm 03/2004
!
!	The floodplain is allowed  to drain after the river has gone
! below flood stage. This allows the floodplain to totally drain.
!
! hydra.v2.6 -- has alteration to the velocity and initiation algorithm 03/2004
!
!	The velocity is developed from an algorithm that looks more
! like the Chezy streamflow equation and produces discharge timing
! in better agreement with the observations at Obidos. 
!	The height at flood initiation is a power function of 
! upstream area rather than log.
!
! hydra.v2.7 -- has alteration to the velocity 03/2004
!
! changed velocity calculation so that now is a function of the 
! dependence on gradient and scale of the river (represented by the 
! discharge).
!
! hydra.v2.8 -- has alteration to the velocity 03/2004
!
! Altered dependence of river velocity on scale of the river to 
! be the wetted perimeter. This is done to conform more closely to
! the empirical Chezy streamflow equation.
!
! hydra.v2.9.f -- alteration to flood calculation 09/2004
!
! Changed the flooding velocity algorithm to be the same structure
! as the river velocity (head and perimeter dependence) but to be 
! slower. I also put in a parameter for the calculation of the flood
! initiation volume. It is a simple multiplier. I realized that the
! equation by itself created too much flooded area on the mainstem.
! This is not surprising since the equation was developed from data
! of streams much smaller than the mainstem.  
!
! THMB.v1.f  -- new name, SRTM data 5/23/2005
!
! This version is the same as v2.9 except that it now uses the SRTM 
! DEM averaged to 5-minute resolution and pit-filled with ArcView as 
! the base topography. The name has also been changed to THMB 
! Terrestrial Hydrology Model with Biogeochemistry because of the 
! copyright infringement issue.
!
! THMB.v1.2.f  -- new name, SRTM data 5/23/2005
!
! THis is the same as THMB.v1.1.f except that the height of the 
! water in the cell is calulated as a function of the standard
! deviation of the cell elevation. In other words, the water height
! can be less than the mean topography of the cell based on a
! standard normal cumulative distribution function of the elevation.
!
! THMB.v1f.f  -- application version 8/23/2006
!
! This is the final version used for publication. It is the same as v1.2.f
! but is no longer experimental. It has been tuned (effref and volltune)
! so that the r2 of the river discharge at 9 stations and the flooded
! area on the 10 Sippel reaches is maximized.
! The results of this version are SRTM_109
!
! ************************************************************
!
      use comvar
      implicit none
      include 'netcdf.inc'!'/usr/include/netcdf.inc'
!
      INTEGER i,j,i2,j2, ii, jj, k, kcount, kt
      REAL sum
      INTEGER iyyyy,imm,idd, fyyyy,fmm,fdd !RAcosta 20200825
      INTEGER iyyyyrun,immrun,iddrun, fyyyyrun,fmmrun,fddrun !RAcosta 20200825
      !INTEGER nmons2,nmons3 !RAcosta 20200825 
!--------------------------------------------------------------
! PRINT LOGO THMB/BESM DEVELOPMENT - RAcosta 20200818
!--------------------------------------------------------------
!
      CALL LOGO
!    
!--------------------------------------------------------------
! READ information file to determine run characteristics
!--------------------------------------------------------------
!
      OPEN(unit = 23, file = 'THMB.v1.inf',status = 'old')
      READ(23,*)latii
      READ(23,*)latif
      READ(23,*)longi
      READ(23,*)longf
      READ(23,*)dgper
      READ(23,*)nmin
      READ(23,*)startyear
      READ(23,*)nyrs
      READ(23,*)spin
!      READ(23,*)restart
      READ(23,*)wrtrestart  !RAcosta 20200818
      READ(23,*)getrestart  !RAcosta 20200818
      READ(23,*)laket
      READ(23,*)wett
      READ(23,*)irrig
      READ(23,*)outfile
      READ(23,*)monthini    !RAcosta 20200818
      READ(23,*)monthend    !RAcosta 20200818
      READ(23,*)startinput  !RAcosta 20200825
      READ(23,*)finalinput    !RAcosta 20200825
      READ(23,*)startrun    !RAcosta 20200825
      READ(23,*)finalrun    !RAcosta 20200825

      CLOSE(23)
!
!--------------------------------------------------------------
! transfor date information to integer
!--------------------------------------------------------------
!
      READ(startinput(1:4),'(i)') iyyyy
      READ(startinput(5:6),'(i)') imm
      READ(startinput(7:8),'(i)') idd
      READ(finalinput(1:4),'(i)') fyyyy
      READ(finalinput(5:6),'(i)') fmm
      READ(finalinput(7:8),'(i)') fdd
      WRITE(*,*)
      WRITE(*,*) 'initial readclimmate yyyy:mm:dd', iyyyy, imm, idd
      WRITE(*,*)
      WRITE(*,*) ' finish  readclimate yyyy:mm:dd', fyyyy, fmm, fdd
      WRITE(*,*)
!
! transfor date information to integer
!
      READ(startrun(1:4),'(i)') iyyyyrun
      READ(startrun(5:6),'(i)') immrun
      READ(startrun(7:8),'(i)') iddrun
      READ(finalrun(1:4),'(i)') fyyyyrun
      READ(finalrun(5:6),'(i)') fmmrun
      READ(finalrun(7:8),'(i)') fddrun
      WRITE(*,*)
      WRITE(*,*) 'initial run yyyy:mm:dd', iyyyyrun, immrun, iddrun
      WRITE(*,*)
      WRITE(*,*) ' finish run yyyy:mm:dd', fyyyyrun, fmmrun, fddrun
      WRITE(*,*)

      nmons2 = (fyyyy - iyyyy + 1 )*12
      WRITE(*,*) ':-- nmons2 to READ CLIMATE: ', nmons2 

      iwmon2 = ((iyyyyrun - iyyyy )*12)+1
      WRITE(*,*) ':-- iwmon2:startyear to INIT INDEX MATRIX: ', iwmon2

      startyear = iyyyyrun
      nyrs      = fyyyyrun


!       DO iyear = startyear,nyrs+spin !130
! !
!          IF ((iyear-startyear)+1 .le. spin)then
!             tyear = 1            !reset for spin-up
!             WRITE(*,*) '((iyear-startyear)+1 .le. spin)'
!             WRITE(*,*) iyear,  startyear, spin
!             WRITE(*,*)
!          ELSE
!             tyear = (iyear-spin-startyear)+1
!             WRITE(*,*) 'tyear=(iyear-spin-startyear)+1'
!             WRITE(*,*) tyear, iyear, spin, startyear
!             WRITE(*,*)
!          ENDIF
! !
!          wyear = (tyear-1) + (startyear -1)
! !
!          WRITE(*,*)'BEGIN YEAR:year/nyrs/tyear/wyear '
!          WRITE(*,*) iyear, nyrs,tyear,wyear
! !
! ! Leap year
!          IF (mod(wyear,4) .eq. 0)then
!             ndaypm(2) = 29
!          ELSE
!             ndaypm(2) = 28
!          ENDIF

!           !DO imon = 1,12 !131
!           IF(iyear.eq.fyyyyrun) monthend = fmmrun
!           DO imon = monthini, monthend
!             WRITE(*,*)' --- virou mes! : ', imon
!             WRITE(*,*)
!           ENDDO
!       ENDDO

!       STOP

!
      dy = circ*((nmin/60.)/360.)
      scale2=(60./nmin)*dgper
!
      j2start = 1
      j2end = 180./dgper
      i2start = 1
      i2end = 360/dgper
!
      nr = (j2end+1)-j2start
      nc = (i2end+1)-i2start
!
      j3start = 1
      j3end = (180*60)/nmin
      i3start = 1
      i3end = (360*60)/nmin
!
      nrf = (j3end+1)-j3start
      ncf = (i3end+1)-i3start
!
      lati=nint(latii)
      latf=nint(latif)
      loni=nint(longi)
      lonf=nint(longf)
      WRITE(*,*)
      WRITE(*,*) "ilatINT:", lati, "flatINT:",latf
      WRITE(*,*)
      WRITE(*,*)
      WRITE(*,*) "ilonINT:", loni, "flonINT:",lonf
      WRITE(*,*)
      WRITE(*,*) "ncf: ", ncf, "nrf : ", nrf
      WRITE(*,*)
!     
      lonfst = (-180.+((i2start-1)*dgper))+(dgper*(1/(60/nmin))) 		!HighRes
      latfst = (90.-((j2start-1)*dgper))-(dgper*(1/(60/nmin)))			!HighRes
      lonend = (-180.+((i2end-1)*dgper))+(dgper*(1/(60/nmin)))			!HighRes
      latend = (90.-((j2end-1)*dgper))-(dgper*(1/(60/nmin)))			!HighRes
      WRITE(*,*)
      WRITE(*,*) "latfst:", latfst, "latend:",latend
      WRITE(*,*)
      WRITE(*,*)
      WRITE(*,*) "lonfst:", lonfst, "lonend:",lonend
      WRITE(*,*)
!
! Number of  points at  0,5 graus - climate data
!       inr = (abs(latf)+abs(lati))/dgper						!HighRes TODO urgente!!!!!!
!       inc = int(lonf-loni)/dgper							!HighRes
      
!       WRITE(*,*)
!       WRITE(*,*) "inc: ", inc, "inr : ", inr
!       WRITE(*,*)


! Number of  points at  0,1 graus - climate data
      inr = (abs(latif)+abs(latii))/dgper            !HighRes TODO urgente!!!!!!
      inc = (longf-longi)/dgper                   !HighRes
      
      WRITE(*,*)
      WRITE(*,*) "inc: ", inc, "inr : ", inr
      WRITE(*,*)
      

! Start and end poits for the run
      jini = ((nint(latfst)-lati)*(60/nmin))+1					!HighRes
      jend = (jini+(dgper*(60/nmin)*inr))-1					   !HighRes
      iini = ((abs(nint(lonfst)-loni))*(60/nmin))+1			!HighRes
      iend = (iini+(dgper*(60/nmin)*inc))-1					   !HighRes

! Number of points at  5' - morphologic data and output
      incf = iend-iini+1.
      inrf = jend-jini+1.

      !print *, jini, jend, iini, iend, incf, inrf
      WRITE(*,*) "jini :", jini, "jend :", jend
      WRITE(*,*) "iini :", iini, "iend :", iend
      WRITE(*,*) "incf :", incf, "inrf :", inrf



!     stop
      CALL aloca!(nrf, ncf, incf, inrf, nmons)
!
! Define lat and lon grids for writing out 5x5' data.
!
      DO j = 1,nrf
        latscale(j) = (latfst+(1/(60/nmin)))-j*(1/(60/nmin))			!HighRes
      ENDDO

      DO i = 1,ncf
        lonscale(i) = (lonfst-(1/(60/nmin)))+i*(1/(60/nmin))			!HighRes
      ENDDO

      DO i = iini,iend
         along(i-(iini-1)) =  lonscale(i)
      ENDDO

      DO j = jini,jend
         alatg(j-(jini-1)) =  latscale(j)
      ENDDO

      WRITE(*,*)
      WRITE(*,*)'* ------------------------------------------------------ *'
      WRITE(*,*)
      WRITE(*,*)' HIGH-RESOLUTION INDEX '
      WRITE(*,*)' i init= ',iini,'j init = ',jini
      WRITE(*,*)' i  end= ',iend,'j  end = ',jend
      WRITE(*,*)
      WRITE(*,*)' LOW-RESOLUTION INDEX '     
      WRITE(*,*)' inc = ',inc,'inr = ',inr
!      WRITE(*,*)'ncf = ',ncf,'nrf = ',nrf
!      WRITE(*,*)'along = ',along(1),'alat = ',alatg(1)
!      WRITE(*,*)'lonscale = ', lonscale(1),'latscale = ', latscale(1)
      WRITE(*,*)'incf = ',incf,'inrf = ',inrf
      WRITE(*,*)
      WRITE(*,*)'* ------------------------------------------------------ *'
      WRITE(*,*)
      WRITE(*,*)
      WRITE(*,*)'* ------------------------------------------------------ *'
      WRITE(*,*)
      WRITE(*,*)'j2start = ',j2start,'j2end = ',j2end
      WRITE(*,*)'latfst = ',latfst,'lonfst = ',lonfst
      WRITE(*,*)'latend = ',latend,'lonend = ',lonend
      WRITE(*,*)
      WRITE(*,*)'* ------------------------------------------------------ *'
      WRITE(*,*)
      WRITE(*,*)
      WRITE(*,*)'* -------------  RUN - CHARACTERISTICS   --------------- *'
      WRITE(*,*)      
      WRITE(*,*)'1 = ',nint(lonfst),'lon init: ',loni, 'lon end: ',lonf
      WRITE(*,*)'1 = ',nint(latfst),'lat init: ',lati, 'lat end: ',latf
      WRITE(*,*)
      WRITE(*,*)'nyrs to run: ', nyrs
      WRITE(*,*)'  startyear: ', startyear
      WRITE(*,*)'   yrs spin: ', spin
      WRITE(*,*)'wrt_restart: ', wrtrestart
      WRITE(*,*)'get_restart: ', getrestart
      WRITE(*,*)'      irrig: ', irrig
      WRITE(*,*)'     output: ', outfile
      WRITE(*,*)
      WRITE(*,*)
      WRITE(*,*)'* ------------------------------------------------------ *'
      WRITE(*,*)

      IF (wett.ne. 0)then
      WRITE(*,*)
      WRITE(*,*)'* ------------------------------------------------------ *'
      WRITE(*,*)'                Running with flooded areas                '
      WRITE(*,*)'                   FLOODED AREAS ON                       '
      WRITE(*,*)'* ------------------------------------------------------ *'
      ENDIF

      IF (irrig .ne. 0)then
      WRITE(*,*)
      WRITE(*,*)'* ------------------------------------------------------ *'
      WRITE(*,*)'                Running with irrigation                   '
      WRITE(*,*)'                     IRRIGATION ON                        '
      WRITE(*,*)'* ------------------------------------------------------ *'
      ENDIF
      
      CALL READPAR!(timer, timed, effref, refper, effrefl, &
!                         flref, a, b, c, d)


!--------------------------------------------------------------
!--------------------------------------------------------------
!READing Maps - Geomorphology and surface datas
!--------------------------------------------------------------
!--------------------------------------------------------------
! READ input geomorphology files
!
      CALL READ_MAPS!(grid2, stddev, uparea, &
!                     outdir, sinuos, wetmask)
!      WRITE(*,*)'grid2(13,29) = ',grid2(220,50)
!      WRITE(*,*)'outdir(13,29) = ',outdir(220,50)
!      WRITE(*,*)'uparea(13,29) = ',uparea(220,50)

!--------------------------------------------------------------
!--------------------------------------------------------------
!Initialize variables and constants
!--------------------------------------------------------------
! initialize constants
! nyr is the number of years to run the model, delt is the timestep.
!
      ioo    = 1.0e-04   !increase river velocity with smaller value 
      delt   = 3600.
      iday   = 0
      imon   = 1
      iyear  = 1
      itime  = 0
      ktstep = 0
      nspday = 24./(delt/3600.)
      volchk = 0.

!--------------------------------------------------------------
! initialize variables
!
      DO j = 1,inrf           !small array
        DO i = 1,incf
             outelv(i,j) = grid2(i,j)
             height(i,j) = grid2(i,j)
            wetmask(i,j) = wetmask(i,j)/10000.
               voll(i,j) = 0.
               volw(i,j) = 0.
               volb(i,j) = 0.
               volr(i,j) = 0.
               volt(i,j) = 0.
              dvoll(i,j) = 0.
              dvolw(i,j) = 0.
            fluxout(i,j) = 0.
            sfluxin(i,j) = 0.
           fluxoutw(i,j) = 0.
           sfluxinw(i,j) = 0.
              areat(i,j) = 0.
             basin2(i,j) = 0.
              warea(i,j) = 0.
        ENDDO
      ENDDO

!--------------------------------------------------------------
!--------------------------------------------------------------
!Calculations for model initialization
!--------------------------------------------------------------

!--------------------------------------------------------------
! calculate the area of each gridcell phi starts at latitude 
! corresponding to area chosen. Therefore, 90-number of gridcells
! one is starting from top of array. dgper is the number of degrees
! per unit of coarse resolution array.
!
!Victor      phi = ((90.-((j2start-1)*dgper))-(2.5/60.))*rad
!      phi = (90.-(2.5/60.))*rad
!
      DO j = 1,jend
!Victor        dx = cos(phi-((j-1.)*(5./60.)*rad))*dy
        dx=cos(latscale(j)*rad)*dy								!Victor
        area(j) = dx*dy
!        print *, j, dx, dy, area(j)
      ENDDO
!
      basarea = 0.      !total area of the basin of interest
! 
      DO j = 1,inrf
        DO i = 1,incf
          IF (outdir(i,j) .lt. 9.)basarea = basarea + area(j+jini-1)
        ENDDO
      ENDDO
!
      WRITE(*,*) sum
      WRITE(*,*)
      WRITE(*,*)'* ------------------------------------------------------ *'
      WRITE(*,*)
      WRITE(*,*)'total area of Basin (km2) = ',basarea/1.e+06 
      WRITE(*,*)
      WRITE(*,*)'* ------------------------------------------------------ *'
      
!--------------------------------------------------------------
!Calculations for the Cumulative Distribuction Function (CDF)
! This will be used on the loop 122-112
!For calculate the distribuction of water inside the grid with different resolution

      DO k = -400,400
        pdf = (exp(-0.5*((float(k)/100.)**2)))/sqrt(2*pi)
        kcount = k+401
        IF (kcount .gt.1)then
           cdf(kcount) = cdf(kcount-1)+pdf/100.
        ELSE
           cdf(kcount) = pdf/100.
        ENDIF
      ENDDO
!
!--------------------------------------------------------------
!--------------------------------------------------------------
!Lakes Inicialization - Will just go inside here if laket is different of 0
!--------------------------------------------------------------
! starting lakes - 0 if without lakes areas,
!                  1 if predicted lakes

      IF (laket.ne.0)then
      WRITE(*,*)
      WRITE(*,*)'* ------------------------------------------------------ *'
      WRITE(*,*)'                Running with lakes                        '
      WRITE(*,*)'                    LAKES ON                              '
      WRITE(*,*)'* ------------------------------------------------------ *'
      CALL READ_LAKES!(laket, ncf, nrf, alon, alat, outi, outj, &
!                     sillh, alev, ierr, times)
!--------------------------------------------------------------
! Calculate the total volume within each potential lake area.
! first set the volume of the outlet location to initial value
!
         DO j = jini,jend
           DO i = iini,iend
              i2 = i -(iini-1)
              j2 = j -(jini-1)
!
              IF (sillh(i,j) .gt. 0.)then
              ii = min(max(outi(i,j)-(iini-1),1.),real(incf))
              jj = min(max(outj(i,j)-(jini-1),1.),real(inrf))
!              WRITE(*,*)'through lake volume calc', ii, jj, volt(ii,jj)
              IF ((ii .gt. 0).and. (ii .le. iend))then
                 IF ((jj .gt. 0).and. (jj .le. jend))then
                    volt(ii,jj) = volt(ii,jj) + &
                        max(sillh(i,j)-grid2(i2,j2),0.1)*area(j)
                 ENDIF
              ENDIF
              ENDIF
!
           ENDDO
         ENDDO
!
      WRITE(*,*)'through lake volume calc'
      WRITE(*,*)
!--------------------------------------------------------------
! Find the lake kernal location, the location that a lake will
! start to fill from. For now the kernal will be the first pit
! encountered. Later it can be modified to provide a more 
! physiCALLy realistic kernal location.
!
! set basin2 as a mask for lake kernal location. Currently am
! setting as the outlet location. Corrections to that will
! be made outside this loop
!
         DO j = jini,jend
           DO i = iini,iend
              ii = i-(iini-1)
              jj = j-(jini-1)
              IF (outi(i,j) .gt. 0.) THEN
                 IF ((i .eq. int(outi(i,j))) .AND. &
                    (j .eq. int(outj(i,j)))) THEN
                        basin2(ii,jj) = 1.
                 ENDIF
              ENDIF
           ENDDO
         ENDDO
!
! correct lakes where want kernal other than outlet
!
!--------------------------------------------------------------
! Get volume and lake area files if this is a Restart
!
      !   IF (get_restart .eq. 0)then 

       !    CALL get_restart(iwmon, incf, inrf, along, alatg, dvoll, height, sflux, nrf, ncf, &
        !                    alat, alon, alev, times, sfluxin, voll, warea, jini, iini)

         !ENDIF
      ENDIF
      
!
!--------------------------------------------------------------
!--------------------------------------------------------------
! READ netcdf monthly input of runoff and precipitation

      !PRINT *, 'n months defined in read_par.f90: ', nmons
      PRINT *, 'n months defined in comvar.f90: ', nmons2
      
      CALL READ_CLIMATE!(nmons, inc, inr)!, runin, prcpi, evapi, drainin)

      !WRITE(*,*) 'past nc READ'
      WRITE(*,*) ''
!      WRITE(*,*)'runin(13,29,1) = ',runin(30,30,10)+drainin(30,30,10)
!      WRITE(*,*)'prcpi(13,29,1) = ',prcpi(30,30,10)
!      WRITE(*,*)'evapi(13,29,1) = ',evapi(30,30,10)

!--------------------------------------------------------------
!--------------------------------------------------------------
!Startin the time loop - just will be endded at the end of run
!--------------------------------------------------------------
!
!
      WRITE(*,*)
      WRITE(*,*)'* ------------------------------------------------------ *'
      WRITE(*,*)
      WRITE(*,*)'*           -----      MAIN LOOP      -----              *'
      WRITE(*,*)
      WRITE(*,*)'* ------------------------------------------------------ *'
!
      cday = 0
!
      DO iyear = startyear,nyrs+spin !130
!
         IF ((iyear-startyear)+1 .le. spin)then
            tyear = 1            !reset for spin-up
            WRITE(*,*) '((iyear-startyear)+1 .le. spin)'
            WRITE(*,*) iyear,  startyear, spin
            WRITE(*,*)
         ELSE
            tyear = (iyear-spin-startyear)+1
            WRITE(*,*) 'tyear=(iyear-spin-startyear)+1'
            WRITE(*,*) tyear, iyear, spin, startyear
            WRITE(*,*)
         ENDIF
!
         wyear = (tyear-1) + (startyear -1)
!
         WRITE(*,*)'BEGIN YEAR:year/nyrs/tyear/wyear '
         WRITE(*,*) iyear, nyrs,tyear,wyear
!
! Leap year
         IF (mod(wyear,4) .eq. 0)then
            ndaypm(2) = 29
         ELSE
            ndaypm(2) = 28
         ENDIF
!
        PRINT *, 'iyear//iyyyyrun'
        PRINT *, iyear, iyyyyrun

        !DO imon = 1,12 !131
        IF(iyear.eq.iyyyyrun) THEN
          monthini = immrun
        ELSE
          monthini = 1
        ENDIF

        
      PRINT *, 'iyear//fyyyyrun'
      PRINT *, iyear, fyyyyrun
  

        IF(iyear.eq.fyyyyrun) THEN
          monthend = fmmrun
        ELSE
          monthend = 12
        ENDIF

        PRINT *, 'monthini//monthend'
        PRINT *, monthini, monthend

        DO imon = monthini, monthend
          WRITE(*,*)
          WRITE(*,*)': START yyyy/mm:', iyear,'/',imon
          WRITE(*,*)
                    
!
! number of months from start of run
!
            IF (tyear .gt. 1)then
               iwmon  = int((tyear-1)*12 + imon)

               iwmon2 = ((iyyyyrun - iyyyy )*12)+ iwmon
               WRITE(*,*) ':-- iwmon2: LOOP MONTHLY INDEX MATRIX: ', iwmon2
               !iwmon2 = int((tyear-1)*12 + imon) 


               icmon = iwmon + ((startyear-1)*12)
               WRITE(*,*) 'iwmon=(tyear-1)*12 + imon'
               WRITE(*,*) iwmon, tyear, imon
               WRITE(*,*)
               WRITE(*,*) 'icmon=iwmon+((startyear-1)*12)::tyear'
               WRITE(*,*) icmon, iwmon, startyear, tyear
               WRITE(*,*)
            ELSE
               iwmon = imon
               iwmon2 = ((iyyyyrun - iyyyy )*12)+ iwmon
               WRITE(*,*) ':-- iwmon2: LOOP MONTHLY IF tyear < 1: ', iwmon2


               icmon = int(((startyear-1)*12)+imon)
               WRITE(*,*) 'iwmon=imon'
               WRITE(*,*) iwmon, imon
               WRITE(*,*)
               WRITE(*,*) 'icmon=(((startyear-1)*12)+imon)::tyear'
               WRITE(*,*) icmon, startyear, imon, tyear
               WRITE(*,*)

            ENDIF
            
!


! Get volume and lake area files if this is a Restart


          IF(iyear.eq.iyyyyrun) THEN
            mstep2=(imon-monthini+1)*tyear !RAcosta 20200813
            iwmon=mstep2
            
            WRITE(*,*) 'iyear .EQ. iyyyyrun:- mstep2 // iwmon // iwmon2: '
            WRITE(*,*) mstep2 , iwmon, iwmon2
         ELSE
            !mstep2=iwmon
!            mstep2=(imon-monthini+1) + (tyear*12)
            mstep2 = mstep2 + 1
!            mstep2=int((tyear-1)*12 + imon)
            iwmon=mstep2

            WRITE(*,*) 'iyear .NE. iyyyyrun:- mstep2 // iwmon // iwmon2: '
            WRITE(*,*) mstep2 , iwmon, iwmon2
         ENDIF

          WRITE(*,*) '-- MSTEP2:  ', mstep2
          WRITE(*,*)
!
         
        
         IF (getrestart .eq. 0 .and. mstep2 .eq. 1 .and. iyear.eq.iyyyyrun ) THEN 
          
          WRITE(*,*) ' ----------------------------------------'
          WRITE(*,*)
          WRITE(*,*) '        --- ENTREI NO RESTART ---        '
          WRITE(*,*)
          WRITE(*,*) ' ----------------------------------------'

          WRITE(*,*) '                                         '
          WRITE(*,*) 'iwmon: ', iwmon
          WRITE(*,*) '                                         '

           CALL get_restart(iwmon, incf, inrf, along, alatg, height, sflux, sfluxin, fluxout, &
                            nrf, ncf, alat, alon, alev, times, &
                            voll, dvoll, warea, jini, iini)


           WRITE(*,*) ': - WRITE RESTART VARS sflux/sFluxIn/fluxout/voll/dvoll: '
           WRITE(*,*) ': - ', sflux(362,420), sfluxin(362,420), fluxout(362,420)
           WRITE(*,*)
           !WRITE(*,*) ': - WRITE RESTART VARS voll/dvoll/warea/height: '
           !WRITE(*,*) ': - ',voll(362,420),dvoll(362,420),warea(362,420),height(362,420)

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


           !CALL get_restart(iwmon, incf, inrf, along, alatg, height, sflux, nrf, ncf, &
           !                 alat, alon, alev, times, voll, warea, jini, iini)
           
           DO j = 1,inrf
               DO i = 1,incf
                  sflux(i,j)   = 0. !sflux(i,j) RAcosta quando ZEREI ESTAS VARIAVEIS O APORTE DE AGUA DIMINUI SIGNIFICAMENTE
                  sfluxin(i,j) = sfluxin(i,j)
                  fluxout(i,j) = fluxout(i,j)
                  larea(i,j)   = 0.
                  !vollm(i,j)   = voll(i,j)
                  voll(i,j)    = voll(i,j)
                  dvoll(i,j)   = dvoll(i,j)
                  waream(i,j)  = warea(i,j)
                  heightm(i,j) = height(i,j)
               ENDDO
           ENDDO


         ELSE
          WRITE(*,*) '--- NAO ENTREI NO RESTART --'


            DO j = 1,inrf
               DO i = 1,incf
                  sflux(i,j)   = 0.
                  laream(i,j)  = 0.
                  vollm(i,j)   = 0.
                  waream(i,j)  = 0.
                  heightm(i,j) = 0.
               ENDDO
            ENDDO
        ENDIF




!
            WRITE(*,*)'BEGIN MONTH:',imon, iwmon, icmon
            WRITE(*,*)' ----------------------------------------'
            WRITE(*,*)
!
            DO iday = 1,ndaypm(imon) !132
               
               WRITE(*,*) 'Running dd/mm/yyyy: ', iday,imon,iyear
!
! number of days from start of run
!
               cday = cday + 1
               IF ((iyear-startyear)+1 .gt. spin)then
                  wday = wday + 1
               ENDIF
               !WRITE(*,*)'Running hh/dd/mm/yyyy: ', kt,iday,imon,iyear
!
               DO kt = 1,int(nspday) !133
                  ktstep = ktstep + 1
                  

!                   IF( (kt.eq.1 .and. iday.eq.1).OR.(kt.eq.24 .and. iday.eq.ndaypm(imon)) ) THEN
!                      WRITE(*,*)
!                      WRITE(*,*) ':hh/dd/mm/yyyy: ', kt,iday,imon,iyear
!                      WRITE(*,*) ':INIT SUBROUTINE sflux/sFluxIn/fluxout/voll/dvoll: '
!                      WRITE(*,*) ':', sflux(362,420), sfluxin(362,420), fluxout(362,420),voll(362,420),dvoll(362,420)
!                      WRITE(*,*)
!                   ENDIF

                  CALL THMB_loop
!                   WRITE(*,*)
!                   WRITE(*,*) 'FLAG ERROR!'
!                   WRITE(*,*)

!                   IF( (kt.eq.1 .and. iday.eq.1).OR.(kt.eq.24 .and. iday.eq.ndaypm(imon)) ) THEN
!                     WRITE(*,*)
!                      WRITE(*,*) ':hh/dd/mm/yyyy: ', kt,iday,imon,iyear
!                      WRITE(*,*) ':END SUBROTINE sflux/sFluxIn/fluxout/voll/dvoll: '
!                      WRITE(*,*) ':', sflux(362,420), sfluxin(362,420), fluxout(362,420),voll(362,420),dvoll(362,420)
!                      WRITE(*,*)
!                   ENDIF
!
!--------------------------------------
!       WRITE RESTART
!--------------------------------------
!
!            IF( imon.eq.monthend) THEN
!               WRITE(*,*)
!               WRITE(*,*) ': wrtrestart: imon/monthend/iday/ndaypm(imon): '
!               WRITE(*,*) ':', imon,monthend, iday, ndaypm(imon)
!               WRITE(*,*)
!            ENDIF


            IF (wrtrestart.eq.0 .AND. kt.eq.int(nspday) .AND. (iday).eq.ndaypm(imon) .AND. imon.eq.monthend .AND. iyear.eq.nyrs  ) THEN
            !IF (wrtrestart.eq.0 .AND. imon.eq.monthend .AND. (iday-1).eq.ndaypm(imon) ) THEN
              WRITE(*,*)
              WRITE(*,*) ': wrtrestart: imon/monthend/iday/ndaypm(imon): ', imon,monthend, iday, ndaypm(imon)
              WRITE(*,*)
              WRITE(*,*) ': IN RESTART sflux/sFluxIn/fluxout/voll/dvoll: '
              WRITE(*,*) ': ', sflux(362,420), sfluxin(362,420), fluxout(362,420),voll(362,420),dvoll(362,420)
              WRITE(*,*)
              WRITE(*,*)
              WRITE(*,*) '---   escrevendo RESTART  ---- '
              
              CALL w_restart(iwmon, incf, inrf, along, alatg, sflux, height, &
                    velocity, tmpdir, voll, warea, sfluxin, fluxout, dvoll)
              
              WRITE(*,*)
              WRITE(*,*)
              WRITE(*,*) '---    saindo RESTART     ---- '
            ENDIF
!
!---------------------------------------
!       WRITE RESTART DONE
!---------------------------------------
!
               ENDDO !133   kstep loop
            ENDDO !132  daily loop
            
            WRITE(*,*) ': --- OUTPUT MENSAL --- :'
            WRITE(*,*)
            WRITE(*,*)
            WRITE(*,*) 'Running hh/dd/mm/yyyy: ', kt,iday,imon,iyear
            WRITE(*,*)
            WRITE(*,*) ':END SUBROTINE sflux/sFluxIn/fluxout/voll/dvoll: '
            WRITE(*,*) ':', sflux(362,420), sfluxin(362,420), fluxout(362,420),voll(362,420),dvoll(362,420)
            WRITE(*,*)
            WRITE(*,*)


!
!
!---------------------------
!       WRITE OUTPUT MONTHLY
!---------------------------
!
            IF (outfile.eq.1) then
!                  print *, 'sflux: testando!!!!! ', sflux(180,216)
               CALL w_output1!(iwmon, incf, inrf, along, alatg, sflux, heightm, &
!                    velocity, tmpdir, vollm, waream)
            ELSE
               CALL w_output0!(imon,wday,startyear,wyear,iwmon, iday,ndaypm, & 
!                    incf, inrf, along, alatg, iini, jini, &
!                    sflux, heightm, velocity, tmpdir, vollm, waream)
            ENDIF
!    
         IF ( iyear .NE. fyyyyrun .AND. imon .ne. fmmrun) THEN
            WRITE(*,*)': DONE yyyy/mm:', iyear,'/',imon
         ENDIF

         ENDDO   ! end imonth loop 131
!
         IF ( iyear .NE. fyyyyrun) THEN
            WRITE(*,*)': NEXT YEAR: ', iyear+1
         ENDIF      
      ENDDO    ! end iyear loop 130

      WRITE(*,*)'* ------------------------------------------------------ *'
      WRITE(*,*)
      WRITE(*,*)'                 RUN COMPLETED                            '
      WRITE(*,*)
      WRITE(*,*)'* ------------------------------------------------------ *'
!
      STOP 0
      END
