! ---------------------------------------------------------------------
SUBROUTINE THMB_loop!(nrf, ncf, incf, inrf, nmons)
!-----------------------------------------------------------------------
!     Calculos a serem realizados a cada timestep
!---------------------------Code history--------------------------------
!     
!
! ---------------------------------------------------------------------
      use comvar, only: jini, &
                        jend, &
                        iini, &
                        iend, &
                        area, &
                        incf, &
                        inrf, &
                        dx, &
                        dy, &
                        scale2, &
                        along, &
                        alatg, &

!Variaveis da configuracao da rodada
                        laket, &
                        wett, &
                        irrwgt, &
                        joff, &
                        ioff, &
                        irrig, &
!Variaveis lidas no arquivo de parametros
                        a, &
                        b, &
                        c, &
                        d, &
                        timer, &
                        delt, &
                        timed, &
                        effref, &
                        refper, &
                        effrefl, &
                        flref, &
!Configuracoes toemporais
                        ktstep, &
                        imon, &
                        ndaypm, &
                        startyear, &
                        tyear, &
                        imon, &
                        nspday, &
                        wday, &
                        icmon, &
                        iwmon2, &
                        iday, &
                        nmons2, &
                        monthini, &
                        monthend, &
                        wrtrestart, &
!Entradas dos mapas de superficie e lagos
                        sinuos, &
                        uparea, &
                        wetmask, &
                        stddev, &
                        grid2, &
                        outdir, &
                        outi, &
                        outj, &
                        sillh, &
!Entradas dos mapas climaticos
                        runin, &
                        drainin, &
                        prcpi, &
                        evapi, &
                        prec, &
!variaveis inicializadas em 0 no programa principal
                        volw, &
                        volr, &
                        volb, &
                        warea, &
                        dvoll, &
                        voll, & 
                        volt, &
                        dvolw, &
!                        larea, &
                        sfluxin, &
                        fluxout, &
                        sfluxinw, &
                        fluxoutw, &
                        outelv, &
                        areat, &
                        basin2, &
                        height, &
                        ioo, &
                        laream, &
                        warea, &
                        cdf, &
!Variaveis para arquivos de saida
                        sflux, &
                        heightm, &
                        velocity, &
                        tmpdir, &
                        waream, &
                        temp, &
                        tempdl, &
                        tempdw, &
                        tempdr, &
                        af, &
                        ddist, &
                        larea,	&
                        vollir,	&	!Height Correction
                        warear		!Height Correction


      implicit none
      include 'netcdf.inc'



      INTEGER :: i,j,k,k2,i1,j1,i2,j2,i3, j3,i4,j4, ix, jx, ii, jj, iii, ijj,tmpdir2, km
      REAL(KIND=8), PARAMETER :: grideps  = 1.e-10
      REAL(KIND=8) :: gridif
      REAL(KIND=8), PARAMETER :: dveps    = 1.e-10
      REAL(KIND=8), PARAMETER :: dvepsmin = -1.e-10 !RAcosta
      REAL(KIND=8), PARAMETER :: dvepsmax = +1.e-10 !RAcosta
      REAL :: bin,bout,rin,rout,irrout,prcpl,evapl, tmph
      REAL :: volw5, slength, upareal, heighti,widthi, volltune, volli, stddev1
      REAL :: rarea,ic,io, effvel,effvelf,wperim, reflngth,fllngth, z0



!       WRITE(*,*)
!       WRITE(*,*) 'FLAG 01!'
!       WRITE(*,*)



!       IF (wrtrestart.eq.0 .AND. imon.eq.monthend .AND. iday.eq.ndaypm(imon) ) THEN
!                    PRINT *, 'sflux: testando!!!!! ', sflux(362,420)
!                    WRITE(*,*) '--- escrevendo RESTART ---- '
!                    CALL w_restart(iwmon2, incf, inrf, along, alatg, sflux, height, &
!                     velocity, tmpdir, voll, warea, dvoll, sfluxin)

!                    WRITE(*,*) '--- saindo RESTART ---- '
!       ENDIF

! This loop calculates the change in volume with time.
!                              print *, 'teste: ', wday, evapi(10,10,wday)*0.864e+8
!WRITE(*,*)' --- virou mes!'
!WRITE(*,*) ': --- INIT SUBROUTINE sFlux: ', sflux(362,420), fluxout(362,420), nspday, ndaypm(imon)

                  DO j = jini+1,jend-1 !120
                     DO i = iini+1,iend-1 !110
!
                        ii = i - (iini-1)
                        jj = j - (jini-1)

                        

!
                        IF ((outdir(ii,jj) .gt. 0.).AND.(outdir(ii,jj).lt.9)) THEN
                          !WRITE(*,*) ': --- 1 FluxOUT: ', fluxout(ii,jj)
!                                 WRITE(*,*)
!                                 WRITE(*,*) 'FLAG 02!'
!                                 WRITE(*,*)
!
!----------------------------------------------------
! Initialize
!
                           bin     = 0.
                           rin     = 0.
                           bout    = 0.
                           rout    = 0.
                           irrout  = 0.
                           prcpl   = 0.
                           evapl   = 0.
!
! Derive 1/2 degree cell which corresponds to current 5' location
!
                           i2 = int((i-1)/scale2)-int((((iini-1)/scale2)-1))!longitude at 1/2 deg
                           j2 = int((j-1)/scale2)-int((((jini-1)/scale2)-1))!latitude at 1/2 deg
!                           print *, area(j), runin(i2,j2,wday)


!                           print *, icmon, iwmon2
!                           stop
!

! Calculate the daily runoff, precip, and evap input for a given cell
! from the monthly means.
!
                           IF (iday .le. 15) THEN
!
                              km = iwmon2-1
!
                              IF (iwmon2 .eq. 1) km = 12
!
                              rin   = max((runin(i2,j2,km)*area(j)) + &
                                      ((real(iday)+15.)/30.)*         &
                                      ((runin(i2,j2,iwmon2)*area(j)) - &
                                      (runin(i2,j2,km)*area(j))),0.)
                              bin   = max((drainin(i2,j2,km)*area(j)) + &
                                     ((real(iday)+15.)/30.)*            &
                                     ((drainin(i2,j2,iwmon2)*area(j)) -  &
                                     (drainin(i2,j2,km)*area(j))),0.)
                              prcpl =  (prcpi(i2,j2,km)*area(j)) +   &
                                     ((real(iday)+15.)/30.)*         &
                                     ((prcpi(i2,j2,iwmon2)*area(j)) - &
                                     (prcpi(i2,j2,km)*area(j)))
                              evapl = (evapi(i2,j2,km)*area(j)) +      &
                                      ((real(iday)+15.)/30.)*          &
                                      ((evapi(i2,j2,iwmon2)*area(j)) -  &
                                      (evapi(i2,j2,km)*area(j)))
!
                           ELSE
!
                              km = iwmon2+1
!
                              IF (iwmon2 .eq. nmons2) km = nmons2-11

                              rin   = max((runin(i2,j2,iwmon2)*area(j)) + &
                                      ((real(iday)-15.)/30.)*            &
                                      ((runin(i2,j2,km)*area(j)) -       &
                                      (runin(i2,j2,iwmon2)*area(j))),0.)  
                              bin   = max((drainin(i2,j2,iwmon2)*area(j)) + &
                                     ((real(iday)-15.)/30.)*               &
                                     ((drainin(i2,j2,km)*area(j)) -        &
                                     (drainin(i2,j2,iwmon2)*area(j))),0.)   
                              prcpl   = (prcpi(i2,j2,iwmon2)*area(j)) + &
                                        ((real(iday)-15.)/30.)*        &
                                        ((prcpi(i2,j2,km)*area(j)) -   &
                                        (prcpi(i2,j2,iwmon2)*area(j)))
                              evapl   = (evapi(i2,j2,iwmon2)*area(j)) + &
                                        ((real(iday)-15.)/30.)*        &
                                        ((evapi(i2,j2,km)*area(j)) -   &
                                        (evapi(i2,j2,iwmon2)*area(j)))

                           ENDIF
!! Variaveis para verificacao da conservacao de massa
      !WRITE(*,*)
      !WRITE(*,*) 'FLAG 03!'
      !WRITE(*,*)
                          !prec(ii,jj)=prcpl
                          !WRITE(*,*)'prec: ', prec(ii,jj)
!                          etp(ii,jj)=evapl
!                          run(ii,jj)=rin
!                          dren(ii,jj)=bin

!Mudei para a dorma de mensal.... acima, entao comento aqui

!**********************************************
! Aqui para dados horarios e diarios - muda no read_data
!**********************************************


! Calculate the daily runoff, precip, and evap input for a given cell
! from the monthly means.
!
!                              rin   = max((runin(i2,j2,wday)*area(j)),0.)
!                              bin   = max((drainin(i2,j2,wday)*area(j)),0.)
!                              prcpl =  (prcpi(i2,j2,wday)*area(j))
!                              evapl = (evapi(i2,j2,wday)*area(j))

!**********************************************
! Termina a parte que para dados horarios
!**********************************************



! --------------------------------------------------------------------
! IRRIGATION
!
                           IF (irrig .ne. 0)then
!
! This is a model for accounting for irrigation withdrawals from the stream.
! It was used for the Lake Chad study (Coe and Foley, 2001) and no where ELSE.
! Therefore, to use this irrigation values per month at specific locations
! would need to be collected and put into simple time and space if
! statements below.
! 
!                              data irrwgt /.092,.07,.038,.054,.075,.145,.104,.081, &
!                                           .103,.083,.081,.074/
                              IF (tyear+startyear-1 .gt. 1965) THEN
                                 IF (tyear+startyear-1 .gt. 1980)then
                                    IF ((i .eq. 467) .and. (j .eq. 338))then      !Chari/Lagone
                                       irrout = (11.22*irrwgt(imon))*(1.e+09/(ndaypm(imon)*86400.))
                                    ELSE
                                       irrout = 0.
                                    ENDIF
                                 ELSE
                                    IF ((i .eq. 467) .and. (j .eq. 338))then      !Chari/Lagone
                                       irrout = (2.5*irrwgt(imon))*(1.e+09/(ndaypm(imon)*86400.))
                                    ELSE
                                       irrout = 0.
                                    ENDIF
                                 ENDIF
                              ENDIF
                           ENDIF
! --------------------------------------------------------------------
!
! TEST -- This section is used to test the model, in terms of mass 
! convergence etc. I leave it here becuase it is always a very good
! idea to test that mass is conserved after any changes are made to 
! the model.
!
!        IF (iwmon2 .gt. 8)rin = min(rin,150.)
!        IF (imon .ge. 2)then
!         rin = 0.
!         bin = 0.
!         prcpl = 0.
!         evapl = ((13./0.864e+8)*area(j))
!        ENDIF
!
!         bin = (1000./0.864e+8)*area(j)
!         prcpl = (1000./0.864e+8)*area(j)
!         evapl = 0.
!
! I have shut the lake calculation off in the Amazon version of the model
! They are too small to worry about. Therefore, if there is no evap or 
! precip on the lakes then they are not needed. However, if we start 
! calculating the water budget of the wetlands and rivers then we will
! need these terms again.
!
!        IF ((ii .eq. X).and.(jj .eq. Y))then
!         prcpl = 0.
!         evapl = 0.
!         bin   = 0.
!         rin = (1000./0.864e+8)*area(j)
!        ELSE
!         prcpl = 0.
!         evapl = 0.
!         bin   = 0.
!         rin   = 0.
!        ENDIF
!
! end TEST
! --------------------------------------------------------------------
! set i and j locations at smaller grid over which model is being run (ii,jj)
!
! calculate volume in runoff reservoir for land
!
                           rout = volr(ii,jj)/timer
                           volr(ii,jj) = max(volr(ii,jj) + (rin-rout)*delt,0.)
!                           print *, volr(200,200)
!
! calculate volume in baseflow reservoir, land only
!
                           bout = volb(ii,jj)/timed
                           volb(ii,jj) = max(volb(ii,jj) + (bin-bout)*delt,0.)
!----------------------------------------------------
! calculate volume in transport reservoir
! The idea is to calculate the volume
! of each non-grid lake cell as usual. For the lakes
! cells add any water volume to the volume of the
! sill grid cell not to the local cell and don't 
! include the (fluxin-fluxout).
!

                           tempdl(ii,jj) = 0.  !reset to 0. 
                           tempdr(ii,jj) = 0.  !reset to 0.
                           tempdw(ii,jj) = 0.  !reset to 0. 

                           af(ii,jj)=warea(ii,jj)+larea(i,j) !Lake and flood area 

                           temp(ii,jj) = (rout+bout)*(1.-af(ii,jj))

                           tempdr(ii,jj) = (temp(ii,jj)-irrout + (sfluxin(ii,jj) - fluxout(ii,jj)))*delt
!               print *, temp(ii,jj), tempdr(ii,jj)

                           IF (warea(ii,jj) .gt. 0.) THEN
                              tempdw(ii,jj) = tempdw(ii,jj)  +  & 
                                              ((prcpl-evapl)*warea(ii,jj))*delt !have flooded area impact H2O budget
                           ENDIF

                           IF (laket.ne.0) THEN !lake area impact H2O budget
                              i2=min((max(int(outi(i,j)-iini-1),1)),incf)
                              j2=min((max(int(outj(i,j)-jini-1),1)),inrf)
!
                              IF (larea(i,j) .gt. 0.) THEN
                                 tempdl(i2,j2) = tempdl(i2,j2) +   & 
                                                 ((prcpl-evapl)*larea(i,j))*delt !have lake area impact H2O budget
                              ENDIF
                           ENDIF
                           

! calculate the changes to and the total volume in the river reservoir.
! The change in volume (dvoll) is equal the sum of the P-E over the lake 
! and flooded area (tempdl and tempdw) and the river flow (tempdr). 
! Set the minimum value of dvoll.
! Add dvoll to existing river reservoir volume (voll). 

                           dvoll(ii,jj)  = tempdr(ii,jj) + tempdl(ii,jj) + tempdw(ii,jj) 
                           IF ( ABS(dvoll(ii,jj))/delt/area(j) .lt. dveps) dvoll(ii,jj) = 0.
                           !R ACOSTA -- TESTE SEM RESULTADOS
                           !IF ( dvoll(ii,jj) .lt. dvepsmin .OR. dvoll(ii,jj) .gt. dvepsmax) dvoll(ii,jj) = 0.

                           voll(ii,jj)  = MAX(voll(ii,jj) + dvoll(ii,jj), DBLE(0.))
                           IF (abs(voll(ii,jj))/delt/area(j) .lt. dveps) voll(ii,jj) = 0.
                           !R ACOSTA -- TESTE SEM RESULTADOS
                           !IF ( voll(ii,jj) .lt. dvepsmin .OR. voll(ii,jj) .gt. dvepsmax) voll(ii,jj) = 0.
                            !print *, ii, jj

                        ELSE                 !basin .ne. requested number
                           dvoll(ii,jj)  = 9.e+20 
                           voll(ii,jj)   = 9.e+20
                           tempdr(ii,jj) = 9.e+20
                           tempdl(ii,jj) = 9.e+20
                           tempdw(ii,jj) = 9.e+20
                        ENDIF
!
!        print *, dvoll(200,200), tempdl(200,200)

                     ENDDO  !110
                  ENDDO !120
!
!----------------------------------------------------------------
!
                  DO     j = jini+1,jend-1 !123
                     DO     i = iini+1,iend-1 !113
!
                        ii = i - (iini-1)  !small array local cell
                        jj = j - (jini-1)
!
                        IF ((outdir(ii,jj) .gt. 0.).and. (outdir(ii,jj) .lt.9))then
!
!----------------------------------------------------------------
! Calculate changes to and the total floodplain reservoir. It depends on
! the flux in and out of the floodplain and the local addition of any volume
! in the river in excess of the bankfull volume (voll-height*width*length).
!
! Calculate stream length
!
!               print *, outdir(ii,jj)
                           j2 = j + joff(int(outdir(ii,jj)))
                           i2 = i + ioff(int(outdir(ii,jj)))
                           dx = (area(j)/dy+area(j2)/dy)/2.
                           slength = sqrt((dx*dx*abs(i2-i)*abs(i2-i)) + &
                                     (dy*dy*abs(j2-j)*abs(j2-j)))
                           slength = max(dx,slength)*max(sinuos(ii,jj)/100.,0.5)
!
                           upareal = uparea(ii,jj)/1.e+06             !upstream area in km2
!
! May 31, 2005 derived from combination of ANA and HIBAM radar data
!
                           heighti = max((a*(upareal**b)),1.)! Coe et al, 2007 - eq. 1
                           widthi  = max((c*(upareal**d)),1.)! Coe et al, 2007 - eq.2
!
				                   warear(ii,jj) = widthi*slength      !surface area of river at flood initiation !Height Correction


! River volume at flood initiation
!
! I have created a tuning parameter for the initiation volume. 
! It is necessasry because the linear reservoir model does not
! simultaneously calculate discharge and volume. Volume is a function 
! of the velocity (it increases until the discharge reaches equilibrium
! with the input terms). Therefore, the parameter needs to be tuned for 
! each change in the effvel. It also appears to be important to scale
! volltune based on the upstream area of the location.
!
                           IF (upareal .le. 0.8e+06) volltune  = 1.00 !>SRTM_116
                           IF (upareal .gt. 0.8e+06) volltune  = 1.8 
                           IF (upareal .gt. 1.0e+06) volltune  = 1.7
                           IF (upareal .gt. 1.5e+06) volltune  = 1.6
                           IF (upareal .gt. 2.2e+06) volltune  = 1.4   !SRTM_117
                           IF (upareal .gt. 4.2e+06) volltune  = 1.1   !SRTM_117! (>SRTM_116)
                           IF (upareal .gt. 4.4e+06) volltune  = 1.0
!

                           volli = (heighti*widthi*slength)*volltune
				                   vollir(ii,jj) = volli   !volume at flood initiation 				!Height Correction
!                           print *, volli
!
!----------------------------------------------------------------
! Calculate the volume of water that may move from the river to the floodplain
! If there is more water in the river volume than the initiation value
! then allow water to move into floodplain reservoir (volw)
!
                          IF (wett.ne.0) then
                           !
                              IF (voll(ii,jj)-volt(ii,jj)-volli .gt. 0.) THEN !volume for flood initiation
                                 dvolw(ii,jj)  = (sfluxinw(ii,jj) - fluxoutw(ii,jj))* delt+dvoll(ii,jj) 
                              ELSE                                        !if below flood stage  
                                 IF (dvoll(ii,jj) .gt. 0.)then                !drain water out of floodplain
                                    dvolw(ii,jj)  = (sfluxinw(ii,jj) - fluxoutw(ii,jj))*delt
                                 ELSE
                                    dvolw(ii,jj)  = (sfluxinw(ii,jj) - fluxoutw(ii,jj))*delt + dvoll(ii,jj) 
                                 ENDIF
                              ENDIF
!
                              IF (abs(dvolw(ii,jj))/delt/area(j) .lt. dveps) dvolw(ii,jj) = 0.
                           
                              volw(ii,jj)  = max(volw(ii,jj) + dvolw(ii,jj),dble(0.))
                              WRITE(*,*)
                              WRITE(*,*) 'volw(362,420)', volw(362,420)
                              WRITE(*,*)
                              WRITE(*,*) 'dvolw(362,420)', dvolw(362,420)
                              WRITE(*,*)

                              WRITE(*,*)
                              WRITE(*,*) 'voll(362,420)', voll(362,420)
                              WRITE(*,*)
                              WRITE(*,*) 'dvoll(362,420)', dvoll(362,420)
                              WRITE(*,*)
                           !
                           ENDIF
                        ELSE
!
                           !dvoll(ii,jj)  = 9.e+20 !RAcosta 20200813
                           voll(ii,jj)   = 9.e+20
                           tempdr(ii,jj) = 9.e+20 
                           tempdl(ii,jj) = 9.e+20 
                           tempdw(ii,jj) = 9.e+20 
!   
                        ENDIF
! 
                     ENDDO !113
                  ENDDO !123

!----------------------------------------------------------------
! Distribute the dvoll into the existing lake  or flooded area
! It will be evened out in loop 121
!
                  DO  j = jini+1,jend-1 !122
                     DO  i = iini+1,iend-1 !112
!
                        ii = i - (iini-1)
                        jj = j - (jini-1)

                        i2 = min(max(outi(i,j)-(iini-1),0.),real(incf))
                        j2 = min(max(outj(i,j)-(jini-1),0.),real(inrf))

!
                        IF ((outdir(ii,jj) .gt. 0.).and. (outdir(ii,jj) .lt.9))then
!
!  -----------------------------------------------------------------
!  ------------------ FLOODED AREA AND HEIGHT ----------------------
!  -----------------------------------------------------------------

                           IF (wett.ne.0) then
                           !
                              upareal = uparea(ii,jj)/1.e+06
! 
! Calculate the area and height of water on the inundated floodplain
! based on cumulative probability function assuming standard normal distribution of
! the elevation within each 5-minute gridcell -- v1.2, v1f
!
                              IF (volw(ii,jj) .gt. grideps)then
                                 stddev1 = max(stddev(ii,jj),1.)
!
!  Calculate 50% of cell volume, where +/- 4*stddev encomposes 100%
!  of the cell volume. The volume is not a rectangle it is a non-
!  linear function based on pdf. Therefore, scale to the mean of the 
!  standard normal distribution (0.39894228). I believe that means
!  that the area under the curve is 40% of the area of a rectangle.
! 
                                 volw5 = 0.
                                 volw5   = (1*stddev1*area(j)*0.39894228)*0.5
!
! Calculate the % of cell that is water covered for a given volume.
! By using volw-volw5 the lefthand tail of the distribution is accounted
! for. In other words because z0 can be negative and positive the portions 
! of the cell that are less than the mean elevation can be flooded
! as can those that are greater. The  height of the water can be less than
! or greater than the mean cell elevation.
!
! Calculate critical value (z0). It will be negative if the volume on the
! floodplain (volw) is less than the 50% volume (volw5). Any value beyond
! +/-4 is truncated to +/-4 since it makes no statistical difference.
!           
                                 z0 = 0.
                                 IF (volw(ii,jj).gt. 0.)then
                                    z0 = log10(volw(ii,jj)/volw5)
                                 ELSE
                                    z0 = -4
                                 ENDIF
!
! Convert z0 to an integer counter of value 1-801 to access the CDF 
! previously calculated. +/- 400 is the value associated with 0%
! and 100% of a cell being flooded.
!
                                z0 = z0*100
                                z0 = max(min(z0,400.),-400.)+401.
!
! Calculate the area and height of the water
!
                                warea(ii,jj) = cdf(nint(z0))
                                height(ii,jj) = grid2(ii,jj) + max(volw(ii,jj)/ &
                                             (max(warea(ii,jj),0.01)*area(j)),dble(0.))
!
                              ELSE
                                 warea(ii,jj)  = 0.
                                 height(ii,jj) = 0.
                              ENDIF
                           !
                           ENDIF
!
! ------------------------------------------------------------------
! -------------------- LAKE AREA AND HEIGHT ------------------------
!
                           IF (((i2 .gt. 0).and. (j2 .gt. 0)) .and. &
                                (laket .eq. 0))then
!
! if volume in basin > lake volume larea = 1. everywhere
! outelv = sill height
!
                              IF (voll(i2,j2) .ge. volt(i2,j2))then
                                 outelv(ii,jj) = max(outelv(ii,jj) + dble(larea(i,j)* &
                                                 dvoll(i2,j2)/areat(i2,j2)),dble(grid2(ii,jj)))  !0.0 if larea = 0.
                                 larea(i,j) = max(min(sngl(outelv(ii,jj))-grid2(ii,jj),1.),0.)
!
! if there is no volume then larea = 0.
!
                              ELSEIF (int(voll(i2,j2)) .eq. 0)then 
                                 larea(i,j)  = 0.
                                 outelv(ii,jj) = grid2(ii,jj)
                                 areat(i2,j2) = 0.   !probably not necessary
!
                              ELSEIF ((voll(i2,j2).gt.0.).and. &
                                      (voll(i2,j2).lt.volt(i2,j2)))then
!
! if some lake alREADy exists in closed basin distribute dvoll
! evenly to those existing cells 
!
                                 IF (areat(i2,j2) .gt. 0.)then
!
! set outelv if larea > 0. Add depth of water if positive
! or negative. 
!
                                    outelv(ii,jj) = max(outelv(ii,jj) + dble(larea(i,j)* &
                                                    dvoll(i2,j2)/areat(i2,j2)),dble(grid2(ii,jj)))  !0.0 if larea = 0.
!
! If no existing lake; set larea of outlet location = 1.
! Ideally would like to choose a kernal location in a 
! realistic location within a lake. This would be stored
! in array basin2.
!
                                 ELSE   
                                    IF (int(basin2(ii,jj)) .eq. 1.)then
                                       outelv(ii,jj) = max(outelv(ii,jj) + &
                                                           dvoll(ii,jj)/area(j2+(jini-1)),dble(grid2(ii,jj)))
                                       larea(i,j) = max(min(sngl(outelv(ii,jj))- &
                                                            grid2(ii,jj),1.),0.)
                                       areat(i2,j2) = max(area(j)*larea(i,j),0.)
!
                                    ENDIF
                                 ENDIF
!
                              ENDIF
                           ENDIF
!
                        ELSE                 !basin .ne. requested number
                           dvoll(ii,jj) = 9.e+20
                           voll(ii,jj)  = 9.e+20
                        ENDIF
!
                        !fluxout(ii,jj) = 0.
                        sfluxin(ii,jj)  = 0.
                        fluxoutw(ii,jj) = 0.
                        sfluxinw(ii,jj) = 0.
!
                     ENDDO !loop 112
                  ENDDO !loop 122 
 
!--------------------------------------------------------------
! Calculate the flux out direction for the floodplain based on the greatest
! difference in the water height between neighboring cells.
!
                  IF(wett.ne.0) then
                  !
                  DO j = 2,inrf-1 !Loop 124 - small array
                     DO i = 2,incf-1 !Loop 114
!
                        IF ((outdir(i,j) .gt. 0.).and. (outdir(i,j) .lt.9))then
!
                           tmpdir(i,j) = outdir(i,j)
                           ddist(i,j) = (dy**2+(area(j+jini-1)/dy)**2)**.5
!                           print *, outdir
!                           stop
!
                           tmph = height(i,j)
                           DO k2 = 1,8
!      data ioff /0,1,1,1,0,-1,-1,-1/ ! traditional way
!      data joff /-1,-1,0,1,1,1,0,-1/
                              j1 = j + joff(k2)            !neighbor cell
                              i1 = i + ioff(k2)
!                              print *, k2, joff(k2), j, j1, int(tmpdir(i,j))
!                              print *, k2, ioff(k2), i, i1, int(tmpdir(i,j))
                              j4 = j1 + joff(int(tmpdir(i,j)))  !downstream cell of neighbor cell
                              i4 = i1 + ioff(int(tmpdir(i,j)))  !tells if is going to conflicting cell
!
! ------------------------ wetmask constraint --------------------------
                              IF ((((i.ge. 24).and.(i.le.340)).and. &
                                   ((j.ge. 24).and.(j.le.264))))then  !TEST
                                 IF ((wetmask(i1,j1) .gt. 0.).and. &
                                       (warea(i1,j1) .lt. wetmask(i1,j1)))then  ! TEST limiting flooded area to Hess floodplain
                                    IF ((tmph .gt. height(i1,j1)+0.05).and. ((i4.ne.i) .and. &
                                        (j4.ne.j)))then
                                       tmpdir(i,j) = k2
                                       tmph = height(i1,j1)
                                       ddist(i,j) = (dy**2+(area(j1+jini-1)/dy)**2)**.5
                                    ENDIF
                                 ENDIF
                              ELSE
!
! ---------------- Standard function for water head ------------------
!
!                              print *, j1, i1, j4, i4, j, i
                                 IF ((tmph .gt. height(i1,j1)+0.05).and. ((i4.ne.i) .and. &
                                     (j4.ne.j)))then
!
                                    tmpdir(i,j) = k2
                                    tmph = height(i1,j1)
                                    ddist(i,j) = (dy**2+(area(j1+jini-1)/dy)**2)**.5
                                 ENDIF
! ---------------- Standard function for water head ------------------
!
                              ENDIF  !belongs to wetmask
                           ENDDO  !k2 loop
!
                           IF (tmpdir(i,j) .eq. 0.)then
                              tmpdir(i,j) = outdir(i,j)
                              ddist(i,j) = (dy**2+(area(j+jini-1)/dy)**2)**.5
                           ENDIF
!
                        ENDIF   !basin mask
!
                     ENDDO ! 114
                  ENDDO ! 124
                  !
                  ENDIF

!
!
!--------------------------------------------------------------
! Calculate the flux out of the cell, to which direction,
! and the sum of the fluxes into cells. Also calculate the
! water depth and lake area (larea) for each cell. 
!
                  DO j = jini+1,jend-1   !Loop 121
                     DO i = iini+1,iend-1 !111
!
                        ii  = i - (iini-1)
                        jj  = j - (jini-1)
                        iii = min(max(outi(i,j)-(iini-1),1.),real(incf))
                        ijj = min(max(outj(i,j)-(jini-1),1.),real(inrf))
!
                        IF ((outdir(ii,jj) .gt. 0.).AND.(outdir(ii,jj) .lt.9)) THEN
!
! -----------------------------------------------------------------
! calculate the river flow length between cells. It is a function
! of the flow direction and the river sinuosity (both prescribed)
! and of course the latitude and geometry of the relationship
! between cells.
!
                           j2 = j + joff(int(outdir(ii,jj)))
                           i2 = i + ioff(int(outdir(ii,jj)))
                           dx = (area(j)/dy+area(j2)/dy)/2.
                           slength = sqrt((dx*dx*abs(i2-i)*abs(i2-i)) + &
                                          (dy*dy*abs(j2-j)*abs(j2-j)))
                           slength = max(dx,slength)*max(sinuos(ii,jj)/100.,0.5)
!
! -----------------------------------------------------------------
! Calculate the effective velocity of each cell. If it is not 
! a lake cell then velocity in the river is a function of two
! factors, the gradient and the total discharge. This follows
! the basic form of the Chezy equation in that large rivers
! increase in velocity despite the gradient becoming very small.
!
!###################### RIVER ##############################
!
! Chezy formula u = C(R*S)**0.5; velocity is product of a constant (C, effref),
! and the square root of the gradient (S) and the radius (wetted perimeter, R).
!
! ---Gradient dependence---
! As gradient increases velocity increases.
!
                           i3 = i2-(iini-1)
                           j3 = j2-(jini-1)
                           ic = max(grid2(ii,jj)-grid2(i3,j3),1.)/slength
                           io = ioo/(dx/dy)   !scale reference gradient to latitude
                           effvel = effref*sqrt(ic/io) 
                           effvel = max(effvel,effref)  !make sure isn't too slow

!
! ---Wetted Perimeter dependence---
! Scale velocity with river wetted perimeter, as perimeter increases 
! velocity increases. 
! The idea is to mimic the wetted perimeter function in the Chezy
! stream-flow equation.
!
! 1st calculate the river geometry
!
! flood initiation height -- function of uparea in km2
!
                           upareal = uparea(ii,jj)/1.e+06             !upstream area in km2
!
! May 31, 2005 derived from combination of ANEEL and HIBAM radar data
!
                           heighti = max((a*((upareal)**b)),1.)  
!
! flood initiation river width -- uparea is in km2
! Width at flood initiation. Per Marcos, widthi is a constant i.e. the river
! is a rectangle in shape. Therefore, this width will serve as river width at
! all discharge rates.
!
! May 31, 2005 derived from combination of ANEEL and HIBAM radar data
!
                           widthi  = max((c*(upareal**d)),1.)
!
! river area is a function of the current volume
!
                           rarea  = slength*widthi               !river area in that gridcell
!
!Height Correction Include river water area in the total cell warea
		                        IF(warea(ii,jj) .lt. MIN(rarea/area(j),1.)) THEN
		                        	warea(ii,jj) = MIN(rarea/area(j),1.)
	                         	ENDIF

! wetted perimeter is 2*height + width of river.
! use the width weighted by the ratio of the current volume to the volume at
! initiation or the widthi itself, whichever is smaller.
! note that voll/(h*d) = w*(voll/(w*h*d)
!
                           wperim = (voll(ii,jj)/rarea)*2 + min(sngl(voll(ii,jj))/(heighti*slength),widthi)
!
! ---Product---
! Make velocity a product of the gradient and discharge dependence.
! Instead of using discharge to calculate a minimum velocity value
! as done by Donner, and Costa here the velocity is the product of both.
! This is consistent with Chezy equation and ensures that both 
! gradient and momentum determine velocity, not one or other.
!
                           effvel = effvel*(effref*((max(wperim,20.)/refper)**0.5))
                           effvel = max(min(effvel,effref*10.),effref)  !constrain upper and lower
!
!###################### END RIVER ##########################
!
!###################### FLOODPLAIN #########################
! Similar procedure as with the river but make the water velocity 
! slower

                           IF (wett.ne.0) then

! ---Gradient dependence---
! flref is a tunable parameter for the speed of the water on the
! floodplain. The speed impacts the total area flooded by changing 
! how quickly it may move from one cell to the next.
!
                           effvelf = (effrefl/flref)*sqrt(ic/io) 
                           effvelf = max(effvelf,effrefl/flref)  !make sure isn't too slow
!
! ---Wetted Perimeter dependence---
! For flood waters wetted perimeter dependence now becomes length of
! flood water boundary between cells. I am basing this on the % area
! of the cell flooded and the distance between center points of cells
! (slength).
!
                           reflngth = 0.1*slength   !10% of cell flooded
                           fllngth = slength*warea(ii,jj)
!
! ---Product---
                           effvelf = effvelf*((effrefl/flref)*((max(fllngth,1.)/reflngth)**0.5))
                           effvelf = min(effvelf,effrefl) !constrain upper bound
                           
                           ENDIF
!
!#################### END FLOODPLAIN #######################
!
                           velocity(ii,jj) = effvel                     !WRITE out to check
!
! -----------------------------------------------------------------
! calculate fluxout of each cell and send it as fluxin to
! either the cell downstream if it is not a lake downstream
! or to the outlet location. The fluxout of the cell which
! corresponds to the sill is calculated for only that water 
! volume in excess of the volume required to fill the lake.
!
                           !WRITE(*,*) ': --- 1 FluxOUT: ', fluxout(ii,jj)
                            fluxout(ii,jj) = 0.  !initialize to 0. for new calculation
                            !WRITE(*,*) ': --- reinit FluxOUT: ', fluxout(ii,jj)
!
                           fluxout(ii,jj) = max((voll(ii,jj)-volt(ii,jj))* &
                                                (effvel/slength),dble(0.))

                           !WRITE(*,*) ': --- 2 FluxOUT: ', fluxout(ii,jj)
                           

! make sure fluxout isn't larger than 1/2 of the volume in the cell

                           fluxout(ii,jj) = max(min(dble(fluxout(ii,jj)), &
                                                sfluxin(ii,jj) + temp(ii,jj) + &
                                                ((voll(ii,jj)-volt(ii,jj))/(delt*2.))),dble(0.))
                           !WRITE(*,*) ': --- 3 FluxOUT: ', fluxout(ii,jj)

! Truncate fluxout if too small for computation

                           IF (fluxout(ii,jj)/area(j) .lt. dveps) fluxout(ii,jj) = 0. 
                           i3 = i2-(iini-1)
                           j3 = j2-(jini-1)
                           IF ((i3.gt.0).AND.(i3.le.incf) .AND. (j3.gt.0).AND.(j3.le.inrf)) THEN
                              sfluxin(i3,j3) = sfluxin(i3,j3) + fluxout(ii,jj)
                           ENDIF
! calculate the flux into and out of the floodplain
!
                           fluxoutw(ii,jj) = max(volw(ii,jj)* &
                                                   (effvelf/ddist(ii,jj)),dble(0.))
!
! make sure fluxout isn't larger than 1/2 of the volume in the cell
!
                           fluxoutw(ii,jj) = max(min(dble(fluxoutw(ii,jj)), &
                                                          sfluxinw(ii,jj) + &
                                                         (volw(ii,jj)/(delt*2.))),dble(0.))
!
!
! Truncate fluxoutw if too small for computation
!
                           IF (fluxoutw(ii,jj)/area(j) .lt. dveps) fluxoutw(ii,jj) = 0. 
                           j4 = jj + joff(int(tmpdir(ii,jj)))
                           i4 = ii + ioff(int(tmpdir(ii,jj)))
                           IF ((i4.gt.0).and.(i4.le.incf).and. &
                               (j4.gt.0).and.(j4.le.inrf)) then
                              sfluxinw(i4,j4) = sfluxinw(i4,j4) + fluxoutw(ii,jj)
                           ENDIF
!
!--------------------------------------------------------------
!
! Adjust height of floodplain water column in each cell using the cellular
! automata. 
!
!         IF (warea(ii,jj).gt.0)then  !start with cells alREADy inundated
!         tmpdir2 = 0
!         tmph = height(ii,jj)
!          DO k = 1,8
!           i2 = ii+ioff(k)   !neighbor directions
!           j2 = jj+joff(k)
!           IF ((j2 .gt. 0.).and.(j2.le.inrf))then
!           IF ((i2 .gt. 0.).and.(i2.le.incf))then
!
! look for neighbors that have no water and are lower in elevation
!
!            IF ((height(i2,j2)+1. .lt. height(ii,jj)).and.
!    *         (height(i2,j2)+1 .lt. tmph))then
!             tmpdir2 = k
!             tmph = height(i2,j2)
!            ENDIF
!
!           ENDIF
!           ENDIF
!          ENDDO
!         ENDIF   !warea > 0
!
!          IF (tmpdir2 .ne. 0)then
!
!           j2 = jj + joff(tmpdir2)
!           i2 = ii + ioff(tmpdir2)
!
!           gridif = min(max((height(ii,jj) -
!    *               height(i2,j2)+1.),0.),
!    *               max((height(ii,jj) - grid2(ii,jj)),0.))
!           IF (abs(gridif) .lt. grideps)gridif = 0.
!           gridif = gridif*area(j)*0.5  !move only 0.5 of difference
!
!           height(i2,j2) = height(i2,j2) +
!    *           gridif/area(j)
!           height(ii,jj) = height(ii,jj) -
!    *           gridif/area(j)
!           warea(i2,j2) = min(max(height(i2,j2)-
!    *           grid2(i2,j2),0.),1.)
!
!         ENDIF   !tmpdir ne 1
!         gridif = 0.
!
!--------------------------------------------------------------
! Adjust height of lake water column in each cell using the cellular
! automata. Distribute water height within a lake basin only
!
                           IF ((laket .eq. 0) .and. ((iii .ne. 0).and. (ijj.ne.0)))then
!                              if (mod(iday,2) .eq. 0)then
!
                              IF ((outelv(ii,jj) .gt. grid2(ii,jj)) .and. &
                                    (voll(iii,ijj) .lt. volt(iii,ijj)))then
!
                                 IF ((voll(iii,ijj) .ge. 0.).and. &
                                    (sillh(i,j) .gt. 0.))then
!
                                    tmpdir2 = 0
                                    tmph = outelv(ii,jj)
                                    DO k2 = 1,8
                                       j3 = (j + joff(k2))-(jini-1)
                                       i3 = (i + ioff(k2))-(iini-1)
                                       IF ((outelv(i3,j3) .lt. outelv(ii,jj)) .and. &
                                           (sillh(i3+(iini-1),j3+(jini-1)) .eq. int(sillh(i,j))))then
                                          IF (outelv(i3,j3) .lt. tmph)then
                                             tmpdir2 = k2
                                             tmph = outelv(i3,j3)
                                          ENDIF
                                       ENDIF
                                    ENDDO
!
                                    IF (tmpdir2 .ne. 0.)then
!
                                       j3 = (j + joff(tmpdir2))-(jini-1)
                                       i3 = (i + ioff(tmpdir2))-(iini-1)
!
                                       gridif = min(max((outelv(ii,jj) - outelv(i3,j3)),dble(0.)), &
                                                    max((outelv(ii,jj) - grid2(ii,jj)),dble(0.)))
                                       IF (abs(gridif) .lt. grideps) gridif = 0.
                                       gridif = gridif*area(j)*0.5  !move only 0.5 of difference
!
                                       outelv(i3,j3) = outelv(i3,j3) + &
                                                          gridif/area(j3+(jini-1))
                                       outelv(ii,jj) = outelv(ii,jj) - &
                                                          gridif/area(j)
!
                                    ENDIF
!
                                 ELSE
                                    larea(i,j) = 0.
                                    outelv(ii,jj) = grid2(ii,jj)
                                    areat(iii,ijj) = 0.

                                 ENDIF  !voll(i2,j2) .ge. 0.
                              ENDIF  !outelv .gt. grid2
              
!
! Set lake area = either; 1 if depth is greater than 1 meter, to
! a % of the lake cell equal to the depth of water, or to 0. if
! depth = 0.  depth = outelv-grid
! Adjust the total lake area (areat) to represent current larea.
! Do this by first subtracting current larea from total then
! adding new larea to total. If larea = 0. area added = 0.
!
                              IF ((outi(i,j) .gt. 0.) .and. &
                                  (voll(iii,ijj) .lt. volt(iii,ijj)))then
!
                                 IF ((voll(iii,ijj) .gt. 0.) .and. &
                                           (sillh(i,j) .gt. 0.))then
!
                                    ix = min(max(int(outi(i,j)-(iini-1)),1),incf)
                                    jx = min(max(int(outj(i,j)-(jini-1)),1),inrf)
                                    areat(ix,jx) = max(areat(ix,jx) - area(j)* &
                                                       larea(i,j),0.)
                                    larea(i,j) = max(min(outelv(ii,jj)- &
                                                     grid2(ii,jj),dble(1.)),dble(0.))
                                    areat(ix,jx) = max(areat(ix,jx) + area(j)* &
                                                       larea(i,j),0.) !previous
                                 ELSE
                                    larea(i,j) = 0.
                                    outelv(ii,jj) = grid2(ii,jj)
                                    areat(iii,ijj) = 0.
                                 ENDIF
                              ENDIF
!
                              laream(ii,jj) = (laream(ii,jj)+((outelv(ii,jj)-grid2(ii,jj)) / &
                                                      (ndaypm(imon)*nspday)))
                              IF (outelv(ii,jj)-grid2(ii,jj) .lt. 0.1) laream(ii,jj) = 0.
                           ENDIF        !laket = 0
!
!WRITE(*,*) ': --- 4 FluxOUT: ', fluxout(ii,jj)
!PRINT *
!WRITE(*,*) ': --- 1 sFlux: ', sflux(ii,jj)
                           sflux(ii,jj)   = sflux(ii,jj)+fluxout(ii,jj) / (ndaypm(imon)*nspday)
!!WRITE(*,*) ': --- 2 sFlux: ', sflux(ii,jj), fluxout(ii,jj), nspday, ndaypm(imon)


!                           print *, fluxout(ii,jj), sflux(ii,jj), nspday
                           heightm(ii,jj) = heightm(ii,jj)+max(((height(ii,jj)-grid2(ii,jj)) / &
                                            (ndaypm(imon)*nspday)),0.)

!Height Correction calculate negative height values -- river height when volume is < bankfull volume
!
                           IF((volw(ii,jj) .lt. grideps) .AND. (voll(ii,jj) .gt. grideps)) THEN
                             IF(warea(ii,jj) .gt. 0.05) THEN
 	                             heightm(ii,jj) = heightm(ii,jj) - (log(max(abs((voll(ii,jj)-vollir(ii,jj))/warear(ii,jj)), 1. ))/(ndaypm(imon)*nspday))   ! values less than bankfull
                             ENDIF
                           ENDIF
!
!Height Correction end
                           waream(ii,jj)  = waream(ii,jj)+warea(ii,jj) / (ndaypm(imon)*nspday)
!
                        ELSE
!
                           sflux(ii,jj)    = 9.e+20
                           fluxout(ii,jj)  = 9.e+20
                           sfluxin(ii,jj)  = 9.e+20
                           heightm(ii,jj)  = 9.e+20
                           waream(ii,jj)   = 9.e+20
                           dvoll(ii,jj)    = 9.e+20
                           voll(ii,jj)     = 9.e+20
!
                        ENDIF        !basin = whatever
!
                     ENDDO !111
                  ENDDO !121
!                  print *, ii, jj, sflux(180,216), fluxout(180,216)
!WRITE(*,*) ': --- END SUBROTINE sflux/sFluxIn/fluxout/voll/dvoll: ', 
!WRITE(*,*) ': --- ' sflux(362,420), sfluxin(362,420), fluxout(362,420),voll(362,420),dvoll(362,420), nspday
!
! end of flux calculations
      END SUBROUTINE THMB_loop

