SolveSoilBalance Subroutine

SUBROUTINE SolveSoilBalance   & 
  !
  (time, rain, irrigation, flowdir, vf, vadose)       

IMPLICIT NONE

!Arguments with intent(in):
!INTEGER (KIND = short), INTENT(IN) :: dt !!computation time step [s]
TYPE (DateTime), INTENT(IN) :: time !current date and time
TYPE (grid_real), INTENT(IN) :: rain !!rainfall rate [m/s]
TYPE (grid_real), INTENT(IN) :: irrigation !!irrigation rate [m/s]
TYPE (grid_integer), INTENT(IN) :: flowdir
TYPE (grid_real), INTENT(IN) :: vf !!vegetation fraction [0-1]
TYPE (grid_real), INTENT(IN) :: vadose !vadose zone depth [m]


!local declarations:
INTEGER (KIND = short) :: k, iin, jin, is, js, id, i, j
REAL (KIND = double) :: runoffcell !runoff of single cell [m/s]
REAL (KIND = float) :: percolationcellRZ !root zone percolation of single cell [m/s]
REAL (KIND = float) :: percolationcellTZ !transmission zone percolation of single cell [m/s]
REAL (KIND = float) :: Tdeep !deep surface temperature in lake
REAL (KIND = float) :: alpha, beta !!used to compute actual evapotranspiration [0-1]
REAL (KIND = float) :: sm, wp, fc !!soil moisture, wilting point and field capacity of current cell
REAL (KIND = float) :: actualE, actualT !!actual evaporation and transpiration of current cell
REAL (KIND = double) :: soilDepthCell !!soil depth of current cell [m]
REAL (KIND = float) :: prunoffcum !!cumulated runoff computed by Ross at previous timestep [cm]
REAL (KIND = float) :: pinfilcum !!cumulated infiltration computed by Ross at previous timestep [cm]
REAL (KIND = float) :: pdrncum !!cumulated drainage computed by Ross at previous timestep [cm]
REAL (KIND = float) :: pevapcum !!cumulated evapotranspiration used and possibly modified by Ross at previous timestep [cm]
REAL (KIND = float) :: water !!water in soil [m]
REAL (KIND = float) :: ddx	!actual cell length (m)
REAL (KIND = float) :: meanHydCond  !!mean hydraulic conductivity 
                             !for computing lateral subsurface flow (m/s)
REAL (KIND = float) :: ds !!thickness of the saturated depth (m)
REAL (KIND = float) :: cellWidth !!cell width (m)
REAL (KIND = float) :: ijSlope
LOGICAL             :: snowCovered, iceCovered


!------------end of declaration------------------------------------------------ 

!reset variables for the time step 
QinSoilSub = 0.
QoutSoilSub = 0.
runoff = 0.
rainBalance = 0.


!update PET
!----------
CALL PETupdate (time)


!compute lateral subsurface intercell fluxes
!------------------------------------------
DO i = 1, mask % idim
    DO j = 1, mask % jdim
        IF ( mask % mat (i,j) == mask % nodata ) CYCLE
        
         !check balance id
         IF ( balanceId % mat (i,j) == LANDPLAIN .OR. &
              balanceId % mat (i,j) == LAKE ) THEN
             CYCLE
         END IF
         
        !set downstream cell  (is,js)
        CALL DownstreamCell(i, j, flowdir % mat(i,j), is, js, ddx, flowdir ) 
        
        IF ( CheckOutlet (i, j, is, js, flowdir) ) CYCLE
      
        !harmonic mean saturated conductivity
        meanHydCond = 2. * ksat_sub % mat (i,j) * ksat_sub % mat (is,js) / &
                    ( ksat_sub % mat (i,j) + ksat_sub % mat (is,js) )
      
        !thickness of the saturated depth
        ds = soilDepthTZ % mat(i,j) * soilMoistureTZ % mat (i,j) / thetas % mat (i,j)
        
        !cell width
        cellWidth = ( CellArea (soilMoisture, i, j) ) ** 0.5
        
        !local slope
        ijSlope = ( dem % mat (i,j) - dem % mat (is,js) ) / ddx
        
        IF ( ijSlope <= 0. ) THEN
            ijSlope = 0.0001
        END IF
      
        QoutSoilSub % mat (i,j) = cellWidth * ds * meanHydCond *  ijSlope
      
        !output becomes input of downstream cell
        QinSoilSub % mat (is,js) = QinSoilSub % mat (is,js) + QoutSoilSub % mat (i,j)
      
    END DO
END DO

!loop through mask
DO i = 1, mask % idim 
  DO j = 1, mask % jdim
     IF ( mask % mat (i,j) == mask % nodata ) CYCLE
        
     !reset fluxes of current cell
     runoffcell = 0.
     percolationcellRZ = 0.
     percolationcellTZ = 0.
     
     !set balance id
     id = balanceId % mat (i,j)
     
     !rain contributes to soil balance
     rainBalance % mat (i,j) = rain % mat (i,j)
     
     !add irrigation
     IF ( ALLOCATED ( irrigation % mat ) ) THEN
         rainBalance % mat (i,j) = rainBalance % mat (i,j) + irrigation % mat (i,j)
     END IF
    
     snowCovered = .FALSE.
     iceCovered  = .FALSE.
     !check if snow and ice cover the current cell
     IF ( dtIce > 0 ) THEN  !glaciers are simulated
        IF ( icewe % mat (i,j) <= 0. ) THEN !liquid water in ice 
                                            !contributes to soil balance
            rainBalance % mat (i,j) =  rainBalance % mat (i,j) + &
                                       waterInIce % mat (i,j) / dtSoilBalance
            waterInIce % mat (i,j) = 0.
            IF ( dtSnow > 0 ) THEN !snow is simulated
               IF ( swe % mat (i,j) <= 0. ) THEN !liquid water in snow
                                                 ! contributes to soil balance
                  rainBalance % mat (i,j) =  rainBalance % mat (i,j) + &
                                  waterInSnow % mat (i,j) / dtSoilBalance
                  waterInSnow % mat (i,j) = 0.
               ELSE
                  rainBalance % mat (i,j) = 0.
                  snowCovered = .TRUE.
               END IF
            END IF
        ELSE  !soil is frozen, no vertical fluxes in it 
            !runoff % mat(i,j) = 0.
            !infilt % mat(i,j) = 0.
            !percolation % mat(i,j) = 0.
            !et % mat(i,j) = 0.
            rainBalance % mat (i,j) = 0.
            iceCovered = .TRUE.
            
        END IF
     END IF
    
    IF ( dtSnow > 0 .AND. dtIce <= 0 ) THEN
       IF ( swe % mat (i,j) > 0.) THEN
          !soil is frozen, no vertical fluxes in it 
          !runoff % mat(i,j) = 0.
          !infilt % mat(i,j) = 0.
          !percolation % mat(i,j) = 0.
          !et % mat(i,j) = 0.
          rainBalance % mat (i,j) = 0.
          snowCovered = .TRUE.
       ELSE
          rainBalance % mat (i,j) = rainBalance % mat (i,j) + &
                                    waterInSnow % mat (i,j) / dtSoilBalance
          waterInSnow % mat (i,j) = 0.
       END IF
    END IF
      
   !decide if storm or interstorm period
    IF ( rainFlag % mat(i,j) == 0 .AND. rainBalance % mat (i, j) >= thresholdStartEvent ) THEN
			!storm starts
			rainFlag % mat(i,j) = 1
			interstormDuration % mat(i,j) = 0 !reset interstorm duration
            
			IF (infiltrationModel == SCS_CN) THEN
				CALL Sat2scn (curveNumber % mat(i,j), storativity % mat(i,j), &
					                soilSatRZ % mat(i,j), sEff % mat(i,j) )
			END IF
      
            IF (infiltrationModel == PHILIPEQ .OR. &
                  infiltrationModel == GREEN_AMPT) THEN
				        ism % mat(i,j) = soilMoistureRZ % mat(i,j)
            END IF
  
            IF (infiltrationModel == ROSS_BC .OR. &
                  infiltrationModel == ROSS_VG) THEN
				        runoffcum % mat(i,j) = 0.
                drncum % mat(i,j) = 0.
                infilcum % mat(i,j) = 0.
                evapcum % mat(i,j) = 0.
            END IF
     

    ELSE IF ( rainFlag % mat(i,j) == 1 .AND. rainBalance % mat (i, j) >= thresholdStartEvent ) THEN
			!storm continues
			rainFlag % mat(i,j) = 2
    ELSE IF ( rainBalance % mat (i, j) < thresholdStartEvent ) THEN
			!update interstorm duration
			interstormDuration % mat(i,j) = interstormDuration % mat(i,j) + dtSoilBalance
    END IF
	!check if storm event must be terminated
	IF (interstormDuration % mat(i,j) >= interstorm) THEN
          rainFlag % mat(i,j) = 0
            
          IF (infiltrationModel == SCS_CN) THEN
			      !reset cumulated rainfall
			      raincum % mat(i,j) = 0.
          END IF
          
          IF (infiltrationModel == PHILIPEQ .OR. &
              infiltrationModel == GREEN_AMPT ) THEN
			      !reset cumulated infiltration
			      cuminf % mat(i,j) = 0.
          END IF
    END IF
    
    !compute vertical fluxes according to soil balance computation scheme 
    SELECT CASE ( id )
    CASE (LAKE)
       !Lake drains subsurface flux
       !runoff can become < 0 when pet > rain + QinSoilSub
       runoffcell = rainBalance % mat (i, j) - pet % mat(i,j) + &
                  QinSoilSub % mat(i,j) / CellArea (QinSoilSub, i, j) 
       et % mat(i,j) = pet % mat(i,j) !et is potential 
       percolationcellRZ = 0. ! percolation is zero
       percolationcellTZ = 0.
       
       
    !CASE (LAKE_EWB)
    !    Tdeep = 25. ![°C]
        ! GR TODO
        !CALL ET_energy_nr_lago(SM%mat(iin,jin),albedo%mat(iin,jin),bpress%mat(iin,jin),&
		      !                              dtm%mat(iin,jin),temperatura_oss%mat(iin,jin),&
		      !                              radiazione_oss%mat(iin,jin),vento_oss%mat(iin,jin),umidita_oss%mat(iin,jin),&
		      !                              indice_nuv%mat(iin,jin),ET%mat(iin,jin),&
		      !                              Ts%mat(iin,jin),Rnetta%mat(iin,jin),Xle%mat(iin,jin),Hse%mat(iin,jin),&
		      !                              Ge%mat(iin,jin),Tdeep)                    	   
			!runoffcell = rainBalance % mat (i, j) - et % mat(i,j) + &
   !               QinSoilSub % mat(i,j) / CellArea (QinSoilSub, i, j)
   !        percolationcell = 0. !deep percolation is zero
           
    CASE DEFAULT ! all other cells(SLOPE, SLOPE_EWB, SLOPE_EWB_WET, CHANNEL, CHANNEL_EWB, &
         !CHANNEL_EWB_WET, PLAIN, PLAIN_EWB, PLAIN_EWB_WET)
      
         !compute infiltration and runoff during storm period
         IF ( rainFlag % mat(i,j) /= 0 ) THEN 
        
              SELECT CASE (infiltrationModel)  
          
              CASE (SCS_CN) !use soil conservation services curve number method
                infilt % mat (i,j) = SCScurveNumber (rain = rainBalance % mat (i, j), &
                                         raincum = raincum % mat (i,j), &
                                         c = abstractionRatio % mat(i,j), &
                                         s = sEff % mat(i,j), &
                                         dt = dtSoilBalance, runoff = runoffcell) 

          
              CASE (PHILIPEQ) !use Philips infiltration model  
                infilt % mat (i,j) = Philip (rain = rainBalance % mat (i, j), &
                                         psic = psic % mat(i,j), &
                                         ksat = ksat % mat(i,j), &
                                         theta = soilMoistureRZ % mat(i,j), &
                                         thetas = thetas % mat(i,j), &
                                         thetar = thetar % mat(i,j), &
                                         psdi = psdi % mat (i,j), &
                                         itheta = ism % mat (i,j), dt = dtSoilBalance, &
                                         cuminf = cuminf % mat (i,j) ) 
         
                runoffcell = rainBalance % mat (i, j) - infilt % mat (i,j)
         
          
              CASE (GREEN_AMPT) !use Green-Ampt infiltration model  
                infilt % mat (i,j) = GreenAmpt (rainFlag % mat(i,j), &
                                         rain = rainBalance % mat (i, j),  &
                                         ksat = ksat % mat(i,j), &
                                         theta = soilMoistureRZ % mat(i,j), &
                                         thetas = thetas % mat(i,j), &
                                         thetar = thetar % mat(i,j), &
                                         phy = phy % mat(i,j), &
                                         itheta = ism % mat (i,j), dt = dtSoilBalance, &
                                         cuminf = cuminf % mat (i,j) )    
         
                runoffcell = rainBalance % mat (i, j) - infilt % mat (i,j) 
         
             !CASE (ROSS_BC)
                 !the subroutine SolveRichardsBC is always called even during interstorm period  
         
              END SELECT
      ELSE !interstorm period, runoff is zero. Low intensity rainfall infiltrates
          runoffcell = 0.
          infilt % mat(i,j) = rainBalance % mat (i, j)
      END IF
      
      !compute actual evapotranspiration
      !---------------------------------
      !IF (id == SLOPE_EWB .OR. id == SLOPE_EWB_WET .OR. id == CHANNEL_EWB .OR. &
      !    id == CHANNEL_EWB_WET .OR. id == PLAIN_EWB .OR. id == PLAIN_EWB_WET) THEN
             
           !solve energy balance to compute actual evapotranspiration  
           !IF ( WetlandIsWet (wetland, iin, jin, time) ) THEN
       
               ! TODO
               !  Tsprof=26.5
				           !
		             !CALL ET_energy_nr_wetland(SM%mat(iin,jin),albedo%mat(iin,jin),bpress%mat(iin,jin),&
		             !                       dtm%mat(iin,jin),temperatura_oss%mat(iin,jin),&
		             !                       radiazione_oss%mat(iin,jin),vento_oss%mat(iin,jin),umidita_oss%mat(iin,jin),&
		             !                       indice_nuv%mat(iin,jin),ET%mat(iin,jin),&
		             !                       Ts%mat(iin,jin),Rnetta%mat(iin,jin),Xle%mat(iin,jin),Hse%mat(iin,jin),&
		             !                       Ge%mat(iin,jin),Tsprof)                    
		             !   
	              !  if(ET%mat(iin,jin)<0)then
		             !   ET%mat(iin,jin)=0
		             ! end if 
           !ELSE
         !      CALL ET_ENERGETICO_suolo(SM%mat(iin,jin),SMsat%mat(iin,jin),SMres%mat(iin,jin),albedo%mat(iin,jin),&
		       !                             res_stom_min%mat(iin,jin),LAI%mat(iin,jin),h_veg%mat(iin,jin),bpress%mat(iin,jin),&
		       !                             dtm%mat(iin,jin),B%mat(iin,jin),fraz_veg%mat(iin,jin),temperatura_oss%mat(iin,jin),&
		       !                             radiazione_oss%mat(iin,jin),vento_oss%mat(iin,jin),umidita_oss%mat(iin,jin),&
		       !                             indice_nuv%mat(iin,jin),WP%mat(iin,jin),FC%mat(iin,jin),ET%mat(iin,jin),&
		       !                             Ts%mat(iin,jin),Rnetta%mat(iin,jin),Xle%mat(iin,jin),Hse%mat(iin,jin),&
		       !                             Ge%mat(iin,jin),Ta_prec%mat(iin,jin),Ts_prec%mat(iin,jin))
		       !     Ta_prec%mat(iin,jin)=temperatura_oss%mat(iin,jin)
         !           Ts_prec%mat(iin,jin)=Ts%mat(iin,jin)
	        !        IF(ET%mat(iin,jin)<0)THEN
		       !        ET%mat(iin,jin)=0
		       !     END IF
		       !     IF (PRESENT(stato_neve)) THEN
					    !IF (stato_neve%mat(iin,jin) /= 0.0) THEN
					    !    ET % mat(iin,jin) = 0.0
					    !    IF(temperatura_oss%mat(iin,jin)<0.)THEN
		       !               Ts%mat(iin,jin)=temperatura_oss%mat(iin,jin)
	        !                ELSE
		       !               Ts%mat(iin,jin)=0.
	        !                END IF
         !                   Rnetta%mat(iin,jin) = 0.0
         !                   Xle%mat(iin,jin) = 0.0
		       !             Ge%mat(iin,jin) = 0.0
		       !             Hse%mat(iin,jin) = 0.0
   		    !            END IF    
           !END IF
      !ELSE 
              
           IF ( snowCovered .OR. iceCovered ) THEN
               et % mat(i,j) = 0.0
           ELSE
              !compute actual evaporation as a fraction alpha of potential
               sm = soilMoistureRZ % mat(i,j)
               wp = wiltingPoint % mat(i,j)
               fc = fieldCapacity % mat(i,j)
               alpha = 0.082 * sm + 9.173 * sm**2. - 9.815 * sm**3.
               IF (alpha > 1.) alpha = 1.
               actualE = pet % mat(i,j) * alpha
			    
               !compute actual transpiration as a fraction beta of potential		
               beta =  (sm - wp) / (fc - wp)
               IF (beta > 1.) beta = 1.
               IF (beta < 0.) beta = 0.
               actualT = pet % mat(i,j) * beta
				   
		      !compute actual evapotranspiration
		      et % mat(i,j) = vf % mat(i,j) * actualT + (1. - vf % mat(i,j)) * actualE
		      IF ( et % mat(i,j) < 0.) THEN
			      et % mat(i,j) = 0.0
		      END IF
         END IF
      !END IF
      
      !compute percolation and capilary rise
      !------------------------------------------------------------------------
      CALL PercolationAndCaprise (id, i, j, rainBalance % mat (i, j), vadose, pet, &
                         soilDepthCell, percolationcellRZ, &
                         percolationcellTZ, runoffcell)
   

      !update soil moisture
      !--------------------
      SELECT CASE (infiltrationModel)
      CASE (SCS_CN, PHILIPEQ, GREEN_AMPT)
          CALL UpdateSoilMoisture (id, i, j, dtSoilBalance, &
                                   soilDepthCell, &
                                   percolationcellRZ, &
                                   percolationcellTZ, &
                                   runoffcell,  &
                                   QinSoilSub % mat (i,j), &
                                   QoutSoilSub % mat (i,j) )   
          
         

      CASE (ROSS_BC)
          !store values of cumulated variables at previous time step
          pevapcum = evapcum % mat(i,j)
          prunoffcum = runoffcum % mat(i,j)
          pinfilcum = infilcum % mat(i,j)
          pdrncum = drncum % mat(i,j)
 
          !update soil moisture profile and cumulated fluxes in cm
          CALL SolveRichardsBC (idt = dtSoilBalance, &
                                rain = rainBalance % mat (i, j) , &
                                pet = et % mat (i,j), &
                                p = parameters (cell % mat(i,j) ), &
                                Scell = S (cell % mat(i,j),:), &
                                dxcell = dx (cell % mat(i,j),:), &
                                nsol = 0, &
                                h0 = h0grid % mat(i,j), &
                                nsteps = nsteps % mat(i,j) , &
                                evap = evapcum % mat(i,j), &
                                runoff = runoffcum % mat(i,j), &
                                infil = infilcum % mat(i,j), &
                                drn = drncum % mat(i,j), &
                                watercontent = water )
          !compute fluxes of this time step in m/s
          runoffcell = ( runoffcum % mat(i,j) - prunoffcum ) / 100. / dtSoilBalance
          IF (runoffcell < 0.) runoffcell = 0.

          infilt % mat (i,j)  = ( infilcum % mat(i,j) - pinfilcum ) / 100. / dtSoilBalance
          IF (infilt % mat (i,j)  < 0.) infilt % mat (i,j)  = 0.

          et % mat (i,j)  = ( evapcum % mat(i,j) - pevapcum ) / 100. / dtSoilBalance
          IF (et % mat (i,j)  < 0.) et % mat (i,j)  = 0.

          percolationcellTZ  = ( drncum % mat(i,j) - pdrncum ) / 100. / dtSoilBalance
          IF (percolationcellTZ  < 0.) percolationcellTZ  = 0.

          !compute average soil moisture
          soilMoisture % mat (i,j) = water / soilDepthCell

       END SELECT
              
    END SELECT

    !update percolation grid
    percolation % mat (i,j) = percolationcellTZ
    
    !update runoff grid
    runoff % mat (i,j) = runoffcell
     
  END DO
END DO
 
RETURN
END SUBROUTINE SolveSoilBalance