Dev

.gitignore

@@ -9,6 +9,7 @@
 *.user
 *.userosscache
 *.sln.docstates
+ubsubsurface.f90
 
 # User-specific files (MonoDevelop/Xamarin Studio)
 *.userprefs
@@ -28,7 +29,7 @@ bld/
 [Ll]og/
 
 # Visual Studio 2015/2017 cache/options directory
-.vs/
+*.vs/
 # Uncomment if you have tasks that create the project's static files in wwwroot
 #wwwroot/
 

AllocateHydroVars.f90

@@ -31,9 +31,14 @@ Subroutine AllocateHydroVars(HydroParam,MeshParam)
     Allocate(HydroParam%Fvu(MeshParam%Kmax,MeshParam%nEdge))
     Allocate(HydroParam%Fuv(MeshParam%Kmax,MeshParam%nEdge))
     Allocate(HydroParam%u(MeshParam%Kmax,MeshParam%nEdge))
+
+    Allocate(HydroParam%psij(MeshParam%Kmax,MeshParam%nEdge))
+    Allocate(HydroParam%rj(MeshParam%Kmax,MeshParam%nEdge))
+
     Allocate(HydroParam%utang(MeshParam%Kmax,MeshParam%nEdge))
     Allocate(HydroParam%Wu(MeshParam%Kmax,MeshParam%nEdge))
     Allocate(HydroParam%uxyback(MeshParam%Kmax,2,MeshParam%nEdge))
+    Allocate(HydroParam%uArrow(MeshParam%Kmax,3,MeshParam%nEdge))    
     Allocate(HydroParam%ubBack(MeshParam%Kmax,3,MeshParam%nElem))
     Allocate(HydroParam%uxy(MeshParam%Kmax,2,MeshParam%nEdge))
     Allocate(HydroParam%uxyL(MeshParam%Kmax+1,2,MeshParam%nElem))
@@ -49,9 +54,20 @@ Subroutine AllocateHydroVars(HydroParam,MeshParam)
     Allocate(HydroParam%psi_edge(MeshParam%Kmax,MeshParam%nEdge))
     Allocate(HydroParam%psi_cell(MeshParam%Kmax,MeshParam%nElem))
 
+    Allocate (HydroParam%ugt(MeshParam%nEdge,MeshParam%kMax+1))
+    Allocate (HydroParam%vgt(MeshParam%nEdge,MeshParam%kMax+1))
+    Allocate (HydroParam%wgt(MeshParam%nEdge,MeshParam%kMax+1)) 
+    Allocate(HydroParam%ubt(MeshParam%Kmax,3,MeshParam%nElem))
+    Allocate(HydroParam%ubVt(MeshParam%Kmax,3,MeshParam%nNode))
+    Allocate(HydroParam%uxyt(MeshParam%Kmax,2,MeshParam%nEdge))
+    Allocate(HydroParam%uxyLt(MeshParam%Kmax+1,2,MeshParam%nElem))
+    Allocate(HydroParam%uNodet(MeshParam%Kmax+1,3,MeshParam%nNode))
+    Allocate(HydroParam%wfct(MeshParam%Kmax,MeshParam%nEdge))
+
 
     ! 2.2. Others Variables
     Allocate(HydroParam%etaInf(MeshParam%nElem))
+    Allocate(HydroParam%etaInfn(MeshParam%nElem))
     Allocate(HydroParam%etaplus(MeshParam%nElem))
     Allocate(HydroParam%peta(MeshParam%nNode))
     Allocate(HydroParam%petan(MeshParam%nNode))
@@ -116,25 +132,57 @@ Subroutine AllocateHydroVars(HydroParam,MeshParam)
     !Allocate(HydroParam%Fq(MeshParam%Kmax,MeshParam%nElem))
 
     Allocate(HydroParam%Vol(MeshParam%nElem)) !CAYO
-    Allocate(HydroParam%Hs(MeshParam%nEdge)) !CAYO
+
     Allocate(MeshParam%ei(MeshParam%Kmax,MeshParam%nElem))!CAYO
+    Allocate(MeshParam%Ki(MeshParam%Kmax,MeshParam%nElem))!CAYO
+    Allocate(MeshParam%Ksat(MeshParam%Kmax,MeshParam%nElem))!CAYO   
+    Allocate(MeshParam%nSoil(MeshParam%Kmax,MeshParam%nElem))!CAYO
+    Allocate(MeshParam%alpha(MeshParam%Kmax,MeshParam%nElem))!CAYO    
     Allocate(MeshParam%Kj(MeshParam%Kmax,MeshParam%nEdge))!CAYO
+
+    If(MeshParam%Kmax > 1) Then
+        Allocate(MeshParam%Si(MeshParam%Kmax,MeshParam%nElem))!CAYO 
+        Allocate(HydroParam%SScalar(MeshParam%nElem*MeshParam%Kmax))
+    Else
+        Allocate(MeshParam%Si(MeshParam%subfactor+1,MeshParam%nElem))!CAYO 
+        Allocate(HydroParam%SScalarSaturation(MeshParam%nElem*(MeshParam%subfactor+1)))
+    EndIf
+
     Allocate(HydroParam%us(MeshParam%Kmax,MeshParam%nEdge)) !CAYO
     Allocate(HydroParam%ust(MeshParam%Kmax,MeshParam%nEdge)) !CAYO
+    Allocate(HydroParam%ustang(MeshParam%Kmax,MeshParam%nEdge)) !CAYO
     Allocate(HydroParam%um(MeshParam%Kmax,MeshParam%nEdge))!CAYO
     Allocate(HydroParam%umt(MeshParam%Kmax,MeshParam%nEdge))!CAYO
+    Allocate(HydroParam%umtang(MeshParam%Kmax,MeshParam%nEdge))!CAYO
+    Allocate(HydroParam%wm(MeshParam%Kmax+1,MeshParam%nElem))
+    Allocate(HydroParam%wmt(MeshParam%Kmax+1,MeshParam%nElem))
+
+    Allocate(HydroParam%uxysub(MeshParam%Kmax,2,MeshParam%nEdge))
+    Allocate(HydroParam%ubsub(MeshParam%Kmax,3,MeshParam%nElem))
 
     Allocate(HydroParam%DZsj(MeshParam%Kmax,MeshParam%nEdge)) !CAYO   
     Allocate(HydroParam%DZsjt(MeshParam%Kmax,MeshParam%nEdge)) !CAYO
     Allocate(HydroParam%DZhj(MeshParam%Kmax,MeshParam%nEdge)) !CAYO   
     Allocate(HydroParam%DZhjt(MeshParam%Kmax,MeshParam%nEdge)) !CAYO
 
-    Allocate(HydroParam%DZsi(MeshParam%Kmax,MeshParam%nEdge)) !CAYO   
-    Allocate(HydroParam%DZsit(MeshParam%Kmax,MeshParam%nEdge)) !CAYO
-    Allocate(HydroParam%DZhi(MeshParam%Kmax,MeshParam%nEdge)) !CAYO   
-    Allocate(HydroParam%DZhit(MeshParam%Kmax,MeshParam%nEdge)) !CAYO 
+    Allocate(HydroParam%DZsi(MeshParam%Kmax,MeshParam%nElem)) !CAYO   
+    Allocate(HydroParam%DZsit(MeshParam%Kmax,MeshParam%nElem)) !CAYO
+    Allocate(HydroParam%DZhi(MeshParam%Kmax,MeshParam%nElem)) !CAYO   
+    Allocate(HydroParam%DZhit(MeshParam%Kmax,MeshParam%nElem)) !CAYO 
+    Allocate(HydroParam%PsiCrit(MeshParam%nElem)) !CAYO 
 
     Allocate(HydroParam%DZK(MeshParam%nEdge)) !Sediment Layer !CAYO 
+    Allocate(HydroParam%Gusub(MeshParam%Kmax,MeshParam%nEdge)) 
+
+    Allocate(HydroParam%etak(MeshParam%nElem))
+    Allocate(HydroParam%etam(MeshParam%nElem))
+    Allocate(HydroParam%Qk(MeshParam%nElem))
+    Allocate(HydroParam%Ci(MeshParam%nElem))
+    Allocate(HydroParam%d(MeshParam%nElem))    
+    Allocate(HydroParam%Vol2(MeshParam%nElem))
+    Allocate(HydroParam%Vol1(MeshParam%nElem))
+
+    Allocate(HydroParam%utangNodes(2,MeshParam%nEdge))
 
     ! 3. Hydrodynamic output variables (VTK)
     Allocate(MeshParam%xPoint(MeshParam%nPoint*(MeshParam%Kmax+1)))
@@ -152,5 +200,4 @@ Subroutine AllocateHydroVars(HydroParam,MeshParam)
     Allocate(HydroParam%IndexWaterLevelEdge(MeshParam%nEdge))
 
 
-End Subroutine AllocateHydroVars
-
+End Subroutine AllocateHydroVars

CGOp.f90

@@ -30,6 +30,7 @@ Subroutine CGOp(a,b,dt,HydroParam,MeshParam)
     N = MeshParam%nElem
     b = a                           ! Initial guess
     Call MatOp(b,Ab,dt,HydroParam,MeshParam)
+    ! b - Ab = %Deta - (P+T).eta(k,0)
     r = b - Ab
     pCG = r                         ! Steepest Descent Direction
     alpha = Dot_Product(r,r)        ! Square of the norm of r
@@ -41,6 +42,7 @@ Subroutine CGOp(a,b,dt,HydroParam,MeshParam)
            !Print*, 'The system has been solved with: ',k, 'CGOp iterations'
            Return
        EndIf
+
        Call MatOp(pCG,v,dt,HydroParam,MeshParam)
 
        lambda = alpha/Dot_Product(pCG,v)
@@ -49,6 +51,7 @@ Subroutine CGOp(a,b,dt,HydroParam,MeshParam)
        alphanew = Dot_Product(r,r)
        pCG = r + alphanew/alpha*pCG
        alpha = alphanew
+
     EndDo
 
     Return

CGOpSub.f90

@@ -0,0 +1,60 @@
+Subroutine CGOpSub(a,b,dt,HydroParam,MeshParam)
+
+    ! Matrix Free Conjugate Gradient Method
+    ! Input:
+    ! a -> Input Matrix
+    ! Output:
+    ! b -> Solution
+
+    ! List of Modifications: 
+    !   -> 10.03.2014: Routine Implementation (Rafael Cavalcanti)
+    !   -> 10.03.2014: Fortran Sintax         (Rafael Cavalcanti)
+    ! Programmer: Michael Dumbser
+
+    Use MeshVars !, Only:  nElem, Edge, Neighbor, EdgeLength, Cirdistance
+    Use Hydrodynamic
+    !Use SimulationModel, Only: dt,NearZero
+    Implicit None
+    Real:: alpha, tol, lambda, alphanew 
+    Integer:: k, N
+    Integer:: iElem, iEdge, Pij, Face
+    Real:: Soma, Coef
+    Real:: NearZero = 1e-10
+    Real:: dt
+    type(MeshGridParam) :: MeshParam
+    type(HydrodynamicParam) :: HydroParam
+    Real, intent(in) :: a(MeshParam%nElem)
+    Real, intent(out) :: b(MeshParam%nElem)
+    Real:: r(MeshParam%nElem), Ab(MeshParam%nElem), pCG(MeshParam%nElem), v(MeshParam%nElem)
+
+    N = MeshParam%nElem
+    !%Deta = %F == residual k == V1(eta) + [T - Q(k-1)]eta - d(k-1) = Error in Water Volume
+    b = a                           ! Initial guess
+    !Ab == (P(k,l-1) + T - Q(k-1,l-1)).eta(k,l)
+    Call MatOp(HydroParam%eta,Ab,dt,HydroParam,MeshParam)
+    ! b - Ab = V1(eta) + [T - Q(k-1)]eta - d(k-1) + P(k,l-1).eta(k,l-1)
+    r = b - Ab + HydroParam%P*HydroParam%eta
+    pCG = r                         ! Steepest Descent Direction
+    alpha = Dot_Product(r,r)        ! Square of the norm of r
+    tol = 1e-14                     ! Tolerance
+
+    Do k = 1,1000*N
+       If ( sqrt(alpha) < tol ) Then
+           ! System has been solved
+           !Print*, 'The system has been solved with: ',k, 'CGOp iterations'
+           Return
+       EndIf
+
+       Call MatOp(pCG,v,dt,HydroParam,MeshParam)
+
+       lambda = alpha/Dot_Product(pCG,v)
+       b = b + lambda*pCG
+       r = r - lambda*v
+       alphanew = Dot_Product(r,r)
+       pCG = r + alphanew/alpha*pCG
+       alpha = alphanew
+
+    EndDo
+
+    Return
+End Subroutine CGOpSub

DestroyHydro.f90

@@ -22,7 +22,22 @@ Subroutine DestroyHydro(HydroParam)
     Deallocate(HydroParam%ElSmallm)
     Deallocate(HydroParam%ElCapitalM)
 
+    Deallocate(HydroParam%psij)
+    Deallocate(HydroParam%rj)
+
+    Deallocate(HydroParam%ubt)
+    Deallocate(HydroParam%uxyt)
+    Deallocate(HydroParam%uNodet)  
+    Deallocate(HydroParam%ugt)
+    Deallocate(HydroParam%vgt)
+    Deallocate(HydroParam%wgt)
+    Deallocate(HydroParam%ubVt)
+    Deallocate(HydroParam%uxyLt)
+    Deallocate(HydroParam%wfct)
+
+
     Deallocate(HydroParam%etaInf)
+    Deallocate(HydroParam%etaInfn)
     Deallocate(HydroParam%etaplus)
     Deallocate(HydroParam%peta)
     Deallocate(HydroParam%petan)
@@ -34,6 +49,7 @@ Subroutine DestroyHydro(HydroParam)
     Deallocate(HydroParam%u)
     Deallocate(HydroParam%ut)
     Deallocate(HydroParam%uxyback)
+    Deallocate(HydroParam%uArrow)
     Deallocate(HydroParam%uNode)
     Deallocate(HydroParam%uxy)
     Deallocate(HydroParam%Wu)
@@ -93,6 +109,8 @@ Subroutine DestroyHydro(HydroParam)
     Deallocate(HydroParam%ElSmallms)
     Deallocate(HydroParam%ElCapitalMs)    
 
+    Deallocate(HydroParam%etak)
+    Deallocate(HydroParam%etam)
 
     Deallocate(HydroParam%DZsj) !CAYO
     Deallocate(HydroParam%DZsjt)!CAYO
@@ -104,11 +122,28 @@ Subroutine DestroyHydro(HydroParam)
     Deallocate(HydroParam%DZhit)!CAYO
     Deallocate(HydroParam%Vol)!CAYO
     Deallocate(HydroParam%DZK) !Sediment Layer
+    Deallocate(HydroParam%Gusub) !Sediment Layer
+    Deallocate(HydroParam%PsiCrit) !Sediment Layer
+
+    Deallocate(HydroParam%etak)
+    Deallocate(HydroParam%Qk)
+    Deallocate(HydroParam%Ci)
+    Deallocate(HydroParam%d)    
+    Deallocate(HydroParam%Vol2)
+    Deallocate(HydroParam%Vol1)
 
     Deallocate(HydroParam%us)!CAYO
     Deallocate(HydroParam%ust)!CAYO
+    Deallocate(HydroParam%ustang)!CAYO
     Deallocate(HydroParam%um)!CAYO
     Deallocate(HydroParam%umt)!CAYO
+    Deallocate(HydroParam%umtang)!CAYO
+    Deallocate(HydroParam%wm)
+    Deallocate(HydroParam%wmt)
+    Deallocate(HydroParam%uxysub)
+    Deallocate(HydroParam%ubsub)
+
+    Deallocate(HydroParam%utangNodes)
 
     Return
 End Subroutine DestroyHydro