Path: blob/devel/elmerice/Solvers/GlaDSCoupledSolver.F90
3204 views
!/*****************************************************************************/1! *2! * Elmer/Ice, a glaciological add-on to Elmer3! * http://elmerice.elmerfem.org4! *5! *6! * This program is free software; you can redistribute it and/or7! * modify it under the terms of the GNU General Public License8! * as published by the Free Software Foundation; either version 29! * of the License, or (at your option) any later version.10! *11! * This program is distributed in the hope that it will be useful,12! * but WITHOUT ANY WARRANTY; without even the implied warranty of13! * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the14! * GNU General Public License for more details.15! *16! * You should have received a copy of the GNU General Public License17! * along with this program (in file fem/GPL-2); if not, write to the18! * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,19! * Boston, MA 02110-1301, USA.20! *21! *****************************************************************************/22! ******************************************************************************23! *24! * Authors: Olivier Gagliardini, Mauro Werder25! * Email: [email protected], [email protected]26! * Web: http://www.csc.fi/elmer27! * Address: CSC - Scientific Computing Ltd.28! * Keilaranta 1429! * 02101 Espoo, Finland30! *31! * Original Date: 15 February 201332! *33! *****************************************************************************/34!> Solve for the sheet hydraulic Potential, sheet thickness and channels area35!> similutaneously (GlaDS model) - This solver replace the 3 solvers solving36!> for these 3 variables independently.37!> Equations defined in Werder et al., 201338RECURSIVE SUBROUTINE GlaDSCoupledsolver( Model,Solver,Timestep,TransientSimulation )39!******************************************************************************40!41! ARGUMENTS:42!43! TYPE(Model_t) :: Model,44! INPUT: All model information (mesh,materials,BCs,etc...)45!46! TYPE(Solver_t) :: Solver47! INPUT: Linear equation solver options48!49! REAL(KIND=dp) :: Timestep50! INPUT: Timestep size for time dependent simulations51!52!******************************************************************************53USE Differentials54USE MaterialModels55USE DefUtils56!------------------------------------------------------------------------------57IMPLICIT NONE58!------------------------------------------------------------------------------59!------------------------------------------------------------------------------60! External variables61!------------------------------------------------------------------------------62TYPE(Model_t) :: Model63TYPE(Solver_t), TARGET :: Solver64LOGICAL :: TransientSimulation65REAL(KIND=dp) :: Timestep66!------------------------------------------------------------------------------67! Local variables68!------------------------------------------------------------------------------69TYPE(Solver_t), POINTER :: PointerToSolver70TYPE(Matrix_t), POINTER :: Systemmatrix71TYPE(Nodes_t) :: ElementNodes, EdgeNodes72TYPE(Element_t), POINTER :: Element, Edge, Face, Bulk73TYPE(ValueList_t), POINTER :: Equation, Material, SolverParams, BodyForce, BC, Constants74TYPE(Variable_t), POINTER :: WorkVar, WorkVar275TYPE(Mesh_t), POINTER :: Mesh7677INTEGER :: i, j, k, l, m, n, t, iter, body_id, eq_id, material_id, &78istat, LocalNodes,bf_id, bc_id, DIM, dimSheet, iterC, &79NSDOFs, NonlinearIter, GhostNodes, NonlinearIterMin, Ne, BDForder, &80CoupledIter, Nel, ierror, ChannelSolver, FluxVariable, ThicknessSolver, ierr8182TYPE(Variable_t), POINTER :: HydPotSol83TYPE(Variable_t), POINTER :: ThickSol, AreaSol, VSol, WSol, NSol, &84PwSol, ZbSol, qSol, hstoreSol, QcSol, QmSol, &85MoulinMaskSol, MoulinFluxSol8687INTEGER, POINTER :: NodeIndexes(:), HydPotPerm(:), PwPerm(:), ZbPerm(:), &88ThickPerm(:), VPerm(:), WPerm(:), NPerm(:), AreaPerm(:), &89qPerm(:), hstorePerm(:), QcPerm(:), QmPerm(:),&90CAPerm(:), CFPerm(:), SHPerm(:), &91MoulinMaskPerm(:), MoulinFluxPerm(:)9293REAL(KIND=dp), POINTER :: HydPot(:), HydPotPrev(:,:), ForceVector(:)94REAL(KIND=dp), POINTER :: ThickSolution(:), ThickPrev(:,:), VSolution(:), WSolution(:), &95NSolution(:), PwSolution(:), AreaSolution(:), AreaPrev(:,:), ZbSolution(:), &96qSolution(:), hstoreSolution(:), QcSolution(:), QmSolution(:),&97CAValues(:), CFValues(:), SHValues(:),&98MoulinMask(:), MoulinFluxVal(:)99100CHARACTER(LEN=MAX_NAME_LEN) :: VariableName, SolverName101CHARACTER(LEN=MAX_NAME_LEN) :: SheetThicknessName, ChannelAreaName, ZbName102CHARACTER(LEN=MAX_NAME_LEN) :: methodSheet, methodChannels103104LOGICAL :: Found, FluxBC, Channels, Storage, FirstTime = .TRUE., &105AllocationsDone = .FALSE., SubroutineVisited = .FALSE., &106meltChannels = .TRUE., NeglectH = .TRUE., Calving = .FALSE., &107CycleElement=.FALSE., MABool = .FALSE., MaxHBool = .FALSE., LimitEffPres=.FALSE., &108MinHBool=.FALSE., HaveMoulinMask=.FALSE.109LOGICAL, ALLOCATABLE :: IsGhostNode(:), NoChannel(:), NodalNoChannel(:)110111REAL(KIND=dp) :: NonlinearTol, dt, CumulativeTime, RelativeChange, &112Norm, PrevNorm, S, C, Qc, MaxArea, MaxH, MinH113REAL(KIND=dp), ALLOCATABLE :: MASS(:,:), &114STIFF(:,:), LOAD(:), SheetConductivity(:), ChannelConductivity(:),&115FORCE(:), C1(:), CT(:), OldValues(:), Refq(:)116117REAL(KIND=dp), ALLOCATABLE :: Vvar(:), ublr(:), hr2(:)118119REAL(KIND=dp), ALLOCATABLE :: lr(:), hr(:), Ar(:), Wopen(:), &120ub(:), Snn(:), Ev(:), &121ng(:), alphas(:), betas(:), betac(:), Phi0(:), Phim(:), Afactor(:), Bfactor(:), &122MoulinArea(:), MoulinFlux(:)123124REAL(KIND=dp), ALLOCATABLE :: IceDensity(:), Ac(:), alphac(:), CCt(:), &125CCw(:), lc(:)126REAL(KIND=dp) :: ChannelArea, WaterDensity, gravity, Lw, EdgeTangent(3), &127ALPHA, BETA, CoupledTol, CoupledNorm, PrevCoupledNorm, &128Discharge(3)129130REAL(KIND=dp) :: totat, st, totst, t1131REAL(KIND=dp) :: Np, pw, zb, Wo, he132133REAL(KIND=dp) :: at, at0134135TYPE(ValueHandle_t) :: Load_h136137SAVE &138ElementNodes, EdgeNodes, &139C1, &140CT, &141FirstTime, &142SheetConductivity, &143ChannelConductivity, &144MASS, &145STIFF,LOAD, &146FORCE, &147IsGhostNode, &148AllocationsDone, VariableName, SolverName, NonLinearTol, M, &149lr, hr, Ar, Wopen, Afactor, Bfactor, &150MoulinArea, MoulinFlux, &151ub, Snn, Ev, WaterDensity, gravity, &152ng, alphas, betas, betac, Phi0, Phim, SheetThicknessName, &153ChannelAreaName, ZbName, IceDensity, Ac, alphac, CCt, &154CCw, lc, Lw, NoChannel, NodalNoChannel, &155Channels, meltChannels, NeglectH, BDForder, &156Vvar, ublr, hr2, Refq, Nel,&157Calving, Load_h, LimitEffPres, HaveMoulinMask158159160totst = 0.0_dp161totat = 0.0_dp162163!------------------------------------------------------------------------------164! Get variables needed for solution165!------------------------------------------------------------------------------166VariableName = TRIM(Solver % Variable % Name)167SolverName = 'GlaDSCoupledsolver ('// VariableName // ')'168169CALL ListInitElementKeyword( Load_h, 'Body Force', TRIM(Solver % Variable % Name) // ' Volume Source')170171172IF ( .NOT. ASSOCIATED( Solver % Matrix ) ) RETURN173SystemMatrix => Solver % Matrix174ForceVector => Solver % Matrix % RHS175176PointerToSolver => Solver177178HydPotSol => Solver % Variable179HydPotPerm => HydPotSol % Perm180HydPot => HydPotSol % Values181HydPotPrev => HydPotSol % PrevValues182183LocalNodes = COUNT( HydPotPerm > 0 )184IF ( LocalNodes <= 0 ) RETURN185186!CHANGE (and at all DIMs)187Mesh => Solver % Mesh188DIM = Mesh % MeshDim189M = Mesh % NumberOfNodes190191!------------------------------------------------------------------------------192! Allocate some permanent storage, this is done first time only193!------------------------------------------------------------------------------194IF ( .NOT. AllocationsDone .OR. Mesh % Changed ) THEN195N = Mesh % MaxElementNodes196Ne = Mesh % NumberOfEdges197Nel = Mesh % NumberOfBulkElements + Mesh % NumberOfBoundaryElements198K = SIZE( SystemMatrix % Values )199L = SIZE( SystemMatrix % RHS )200201IF ( AllocationsDone ) THEN202DEALLOCATE( &203ElementNodes % x, &204ElementNodes % y, &205ElementNodes % z, &206EdgeNodes % x, &207EdgeNodes % y, &208EdgeNodes % z, &209C1, &210CT, &211SheetConductivity, &212ChannelConductivity, &213MASS, &214STIFF,LOAD, &215FORCE, &216IsGhostNode, &217lr, hr, Ar, Wopen, Afactor, Bfactor, &218MoulinArea, MoulinFlux, &219ub, Snn, Ev, &220ng, alphas, betas, betac, Phi0, Phim, &221IceDensity, Ac, alphac, CCt, &222CCw, lc, OldValues, NoChannel, NodalNoChannel, &223Vvar, ublr, hr2, &224Refq)225226END IF227228ALLOCATE( &229ElementNodes % x( N ), &230ElementNodes % y( N ), &231ElementNodes % z( N ), &232EdgeNodes % x( N ), &233EdgeNodes % y( N ), &234EdgeNodes % z( N ), &235C1( N ), &236CT( N ), &237SheetConductivity( N ), &238ChannelConductivity( N ), &239MASS( N, N ), &240STIFF( N, N ), LOAD( N ), &241FORCE( N ), &242IsGhostNode( M ), &243lr(N), hr(N), Ar(N), Wopen(N), Afactor(N), Bfactor(N), &244MoulinArea(N), MoulinFlux(N), &245ub(N), Snn(N), Ev(N), &246ng(N), alphas(N), betas(N), betac(N), Phi0(N), Phim(N), &247IceDensity(N), Ac(N), alphac(N), CCt(N), &248CCw(N), lc(N), OldValues(K), NoChannel(M), NodalNoChannel(N), &249Vvar(M), ublr(M), hr2(M), &250refq(dim*M), &251STAT=istat)252253IF ( istat /= 0 ) THEN254CALL FATAL( SolverName, 'Memory allocation error' )255ELSE256CALL INFO(SolverName, 'Memory allocation done', level=1 )257END IF258259! get the ghost nodes of this partition260IF ( ParEnv % PEs > 1 ) THEN !only if we have a parallel run261IsGhostNode( 1:M ) = .FALSE.262GhostNodes = 0;263DO t=1,Solver % NumberOfActiveElements264Element => GetActiveElement(t)265IF (ParEnv % myPe == Element % partIndex) CYCLE266DO i=1,GetElementNOFNodes(Element)267IsGhostNode(Element % NodeIndexes(i)) = .TRUE.268GhostNodes = GhostNodes + 1269END DO270END DO271PRINT *, ParEnv % myPe, ':', GhostNodes, ' ghost nodes'272END IF273274! Find the nodes for which we have no channel (on the boundary)275! Default is False - We allow Channel to growth everywhere276NoChannel = .False.277DO t=1, Mesh % NumberOfBoundaryElements278! get element information279Element => GetBoundaryElement(t)280!IF ( .NOT.ActiveBoundaryElement() ) CYCLE281IF (ParEnv % PEs > 1) THEN282IF (ParEnv % myPe /= Mesh % ParallelInfo % EdgeNeighbourList(t) % Neighbours(1)) CYCLE283END IF284n = GetElementNOFNodes()285IF ( GetElementFamily() == 1 ) CYCLE286287NULLIFY(BC)288BC => GetBC( Element )289bc_id = GetBCId( Element )290NodalNoChannel = .False.291NodalNoChannel(1:n) = GetLogical(BC, 'No Channel BC', Found)292IF (Found) NoChannel(Element % NodeIndexes(1:n)) = NodalNoChannel(1:n)293END DO294295AllocationsDone = .TRUE.296END IF297298299!------------------------------------------------------------------------------300! Read physical and numerical constants and initialize301!------------------------------------------------------------------------------302IF (FirstTime) THEN303FirstTime = .FALSE.304Constants => GetConstants()305306WaterDensity = ListGetConstReal( Constants, 'Fresh Water Density', Found )307IF (.NOT. Found) THEN308WaterDensity = ListGetConstReal( Constants, 'Water Density', Found )309IF (Found) THEN310WRITE(Message,'(A)') 'Constant name >Water Density< has been &311replaced with >Fresh Water Density< due to naming conflict'312CALL WARN(SolverName, Message )313END IF314CALL FATAL(SolverName, 'Constant >Fresh Water Density< not found')315END IF316317gravity = ListGetConstReal( Constants, 'Gravity Norm', UnFoundFatal=.TRUE. )318Lw = ListGetConstReal( Constants, 'Latent Heat', UnFoundFatal=.TRUE. )319320ChannelAreaName = GetString( Constants, &321'Channel Area Variable Name', Found )322IF(.NOT.Found) THEN323CALL WARN(SolverName,'Keyword >Channel Area Variable Name< not found in section Constants')324CALL WARN(SolverName,'Taking default value >Channel Area<')325WRITE(ChannelAreaName,'(A)') 'Channel Area'326END IF327328SheetThicknessName = GetString( Constants, &329'Sheet Thickness Variable Name', Found )330IF(.NOT.Found) THEN331CALL WARN(SolverName,'Keyword >Sheet Thickness Variable Name< not found in section Constants')332CALL WARN(SolverName,'Taking default value >Sheet Thickness<')333WRITE(SheetThicknessName,'(A)') 'Sheet Thickness'334END IF335336ZbName = GetString( Constants, 'Bedrock Variable Name', Found )337IF (Found) THEN338ZbSol => VariableGet( Solver % Mesh % Variables, ZbName, UnFoundFatal=.TRUE. )339ELSE340CALL WARN(SolverName,'Keyword >Bedrock Variable Name< not found in section Constants')341CALL WARN(SolverName,'Taking default value >zb<')342WRITE(ZbName,'(A)') 'Zb'343END IF344345!CHANGE - to get Channel variables added to this solver mesh if346!doing calving and hydrology and consequently having many meshes347Calving = ListGetLogical(Model % Simulation, 'Calving', Found)348IF(.NOT.Found) Calving = .FALSE.349IF(Calving) THEN350DO i=1,Model % NumberOfSolvers351IF(Model % Solvers(i) % Variable % Name == ChannelAreaName) THEN352ChannelSolver = i353EXIT354END IF355END DO356WorkVar => VariableGet(Model % Solvers(ChannelSolver) % Mesh&357% Variables, ChannelAreaName, ThisOnly=.TRUE.)358ALLOCATE(CAPerm(SIZE(WorkVar % Perm)), CAValues(SIZE(WorkVar % Values)))359CAPerm = WorkVar % Perm360CAValues = WorkVar % Values361CALL VariableAdd(Mesh % Variables, Mesh, Solver,&362'Channel Area', 1, CAValues, CAPerm)363WorkVar => VariableGet(Mesh % Variables, 'Channel Area',&364ThisOnly=.TRUE.)365ALLOCATE(WorkVar % PrevValues(SIZE(WorkVar % Values),MAX(Solver&366% Order, Solver % TimeOrder)))367WorkVar % PrevValues(:,1) = WorkVar % Values368369WorkVar => VariableGet(Model % Solvers(ChannelSolver) % Mesh&370% Variables, 'Channel Flux', ThisOnly=.TRUE.)371ALLOCATE(CFPerm(SIZE(WorkVar % Perm)), CFValues(SIZE(WorkVar % Values)))372CFPerm = WorkVar % Perm373CFValues = WorkVar % Values374CALL VariableAdd(Mesh % Variables, Mesh, Solver,&375'Channel Flux', 1, CFValues, CFPerm)376WorkVar => VariableGet(Mesh % Variables, 'Channel Flux',&377ThisOnly=.TRUE.)378ALLOCATE(WorkVar % PrevValues(SIZE(WorkVar % Values),MAX(Solver&379% Order, Solver % TimeOrder)))380WorkVar % PrevValues(:,1) = WorkVar % Values381382!The same for sheet thickness383DO i=1,Model % NumberOfSolvers384IF(Model % Solvers(i) % Variable % Name == SheetThicknessName) THEN385ThicknessSolver = i386EXIT387END IF388END DO389WorkVar => VariableGet(Model % Solvers(ThicknessSolver) % Mesh&390% Variables, SheetThicknessName, ThisOnly=.TRUE.)391ALLOCATE(SHPerm(SIZE(WorkVar % Perm)), SHValues(SIZE(WorkVar % Values)))392SHPerm = WorkVar % Perm393SHValues = WorkVar % Values !Needed to reflect initial condition394CALL VariableAdd(Mesh % Variables, Mesh, Solver,&395'Sheet Thickness', 1, SHValues, SHPerm)396WorkVar => VariableGet(Mesh % Variables, 'Sheet Thickness',&397ThisOnly=.TRUE.)398ALLOCATE(WorkVar % PrevValues(SIZE(WorkVar % Values),MAX(Solver&399% Order, Solver % TimeOrder)))400WorkVar % PrevValues(:,1) = WorkVar % Values401!Necessary to ensure initial condition value reflected in PrevValues402WorkVar % PrevValues(:,1) = WorkVar % Values403NULLIFY(WorkVar)404END IF405406! TODO : implement higher order BDF method407BDForder = GetInteger(GetSimulation(),'BDF Order', Found)408IF (.NOT.Found) BDForder = 1409IF (BDForder /= 1) THEN410WRITE(Message,'(a)') 'Only working for BDF = 1'411CALL FATAL(SolverName, Message)412END IF413END IF ! FirstTime414415SolverParams => GetSolverParams()416417HaveMoulinMask = GetLogical( SolverParams,'Have Moulin Mask', Found)418IF (Found) THEN419WRITE(Message, *) 'Have Moulin Mask set to', HaveMoulinMask420CALL INFO(SolverName,Message,Level=1)421END IF422423NeglectH = GetLogical( SolverParams,'Neglect Sheet Thickness in Potential', Found )424IF ( .NOT.Found ) THEN425CALL FATAL(SolverName, 'No >Neglect Sheet Thickness in Potential< found')426END IF427428methodSheet = GetString( SolverParams, 'Sheet Integration Method', Found )429IF(.NOT.Found) THEN430CALL FATAL(SolverName, 'No >Sheet Integration Methods< found')431ELSE432IF ((methodSheet /='explicit').AND.(methodSheet /='implicit').AND.(methodSheet /= 'crank-nicolson')) THEN433CALL FATAL(SolverName, 'Sheet Integration Method: Implicit, Explicit or Crank-Nicolson')434END IF435END IF436437Channels = GetLogical( SolverParams,'Activate Channels', Found )438IF (.NOT. Found) Channels = .FALSE.439440!CHANGE441!To pick up channel and sheet size limiters, if specified442MaxArea = GetConstReal( SolverParams, &443'Max Channel Area', MABool )444IF ((.NOT. MABool)) CALL WARN(SolverName,'No max channel area specified. &445Channels may grow very large')446447LimitEffPres = GetLogical( SolverParams, &448'Limit Negative Effective Pressure', Found)449IF (.NOT.Found) LimitEffPres= .FALSE.450451MaxH = GetConstReal( SolverParams, &452'Max Sheet Thickness', MaxHBool )453IF ((.NOT. MaxHBool)) CALL WARN(SolverName,'No max sheet thickness specified.&454Sheet may grow very large')455456MinH = GetConstReal( SolverParams, &457'Min Sheet Thickness', MinHBool )458459460IF (Channels) THEN461meltChannels = GetLogical( SolverParams,'Activate Melt From Channels', Found )462IF ( .NOT.Found ) THEN463CALL FATAL(SolverName, 'No >Activate Melt From Channels< found')464END IF465466methodChannels = GetString( SolverParams, 'Channels Integration Method', Found )467IF(.NOT.Found) THEN468CALL FATAL(SolverName, 'No >Channels Integration Methods< found')469ELSE470IF ((methodChannels /= 'explicit').AND.(methodChannels /= 'implicit').AND.(methodChannels /= 'crank-nicolson')) THEN471CALL FATAL(SolverName, 'Channels Integration Method: Implicit, Explicit or Crank-Nicolson')472END IF473END IF474END IF475476NonlinearIter = GetInteger( SolverParams, &477'Nonlinear System Max Iterations', Found )478IF ( .NOT.Found ) NonlinearIter = 1479480NonlinearIterMin = GetInteger( SolverParams, &481'Nonlinear System Min Iterations', Found )482IF ( .NOT.Found ) NonlinearIterMin = 1483484NonlinearTol = GetConstReal( SolverParams, &485'Nonlinear System Convergence Tolerance', Found )486IF ((.Not.Found).AND.(NonlinearIter>1)) CALL FATAL(SolverName,'Need >Nonlinear System Convergence Tolerance<')487488CoupledIter = GetInteger( SolverParams, &489'Coupled Max Iterations', Found)490IF ( .NOT.Found ) CoupledIter = 1491492CoupledTol = GetConstReal( SolverParams, &493'Coupled Convergence Tolerance', Found )494IF ((.Not.Found).AND.(CoupledIter>1)) CALL FATAL(SolverName,'Need >Nonlinear System Convergence Tolerance<')495496ThickSol => VariableGet( Mesh % Variables, SheetThicknessName, UnfoundFatal = .TRUE. )497ThickPerm => ThickSol % Perm498ThickSolution => ThickSol % Values499ThickPrev => ThickSol % PrevValues500IF (HaveMoulinMask) THEN501MoulinMaskSol => VariableGet( Solver % Mesh % Variables, "Moulin Mask", UnFoundFatal=.TRUE. )502IF (.NOT.ASSOCIATED(MoulinMaskSol)) CALL FATAL(SolverName,">Have Moulin Mask< declared but no >Moulin Mask< variable found")503MoulinMaskPerm => MoulinMaskSol % Perm504MoulinMask => MoulinMaskSol % Values505MoulinFluxSol => VariableGet( Solver % Mesh % Variables, "Moulin Flux", UnFoundFatal=.TRUE. )506IF (.NOT.ASSOCIATED(MoulinMaskSol)) CALL FATAL(SolverName,">Have Moulin Mask< declared but no >Moulin Flux< variable found")507MoulinFluxPerm => MoulinFluxSol % Perm508MoulinFluxVal => MoulinFluxSol % Values509END IF510511IF (Channels) THEN512AreaSol => VariableGet( Mesh % Variables, ChannelAreaName, UnfoundFatal = .TRUE. )513AreaPerm => AreaSol % Perm514AreaSolution => AreaSol % Values515AreaPrev => AreaSol % PrevValues516517! flux in the channels (for output only) - edge type variable518QcSol => VariableGet( Mesh % Variables, "Channel Flux" )519IF ( ASSOCIATED( QcSol ) ) THEN520QcPerm => QcSol % Perm521QcSolution => QcSol % Values522END IF523END IF524525! discharge out of the moulins (for output only)526QmSol => VariableGet( Mesh % Variables, "Flux from Moulins" )527IF ( ASSOCIATED( QmSol ) ) THEN528QmPerm => QmSol % Perm529QmSolution => QmSol % Values530END IF531532ZbSol => VariableGet( Mesh % Variables, ZbName )533IF ( ASSOCIATED( ZbSol ) ) THEN534ZbPerm => ZbSol % Perm535ZbSolution => ZbSol % Values536END IF537538dt = Timestep539540!------------------------------------------------------------------------------541! Read the other variables needed542!------------------------------------------------------------------------------543VSol => VariableGet( Mesh % Variables, 'Vclose' )544IF ( ASSOCIATED( VSol ) ) THEN545VPerm => VSol % Perm546VSolution => VSol % Values547END IF548549WSol => VariableGet( Mesh % Variables, 'Wopen' )550IF ( ASSOCIATED( WSol ) ) THEN551WPerm => WSol % Perm552WSolution => WSol % Values553END IF554555NSol => VariableGet( Mesh % Variables, 'Effective Pressure' )556IF ( ASSOCIATED( NSol ) ) THEN557NPerm => NSol % Perm558NSolution => NSol % Values559END IF560561PwSol => VariableGet( Mesh % Variables, 'Water Pressure' )562IF ( ASSOCIATED( PwSol ) ) THEN563PwPerm => PwSol % Perm564PwSolution => PwSol % Values565END IF566567qSol => VariableGet( Mesh % Variables, 'Sheet Discharge' )568IF ( ASSOCIATED( qSol ) ) THEN569qPerm => qSol % Perm570qSolution => qSol % Values571END IF572573hstoreSol => VariableGet( Mesh % Variables, 'Sheet Storage' )574IF ( ASSOCIATED( hstoreSol ) ) THEN575hstorePerm => hstoreSol % Perm576hstoreSolution => hstoreSol % Values577END IF578579!------------------------------------------------------------------------------580! Loop for the coupling of the three equations581!------------------------------------------------------------------------------582! check on the coupled Convergence is done on the potential solution only583PrevCoupledNorm = ComputeNorm( Solver, SIZE(HydPot), HydPot )584585586DO iterC = 1, CoupledIter587588!------------------------------------------------------------------------------589! non-linear system iteration loop590!------------------------------------------------------------------------------591DO iter = 1, NonlinearIter592!------------------------------------------------------------------------------593! print out some information594!------------------------------------------------------------------------------595at = CPUTime()596at0 = RealTime()597598CALL Info( SolverName, ' ', Level=4 )599CALL Info( SolverName, ' ', Level=4 )600CALL Info( SolverName, '-------------------------------------',Level=4 )601WRITE( Message,'(A,A,I3,A,I3)') &602TRIM(Solver % Variable % Name), ' iteration no.', iter,' of ',NonlinearIter603CALL Info( SolverName, Message, Level=4 )604CALL Info( SolverName, '-------------------------------------',Level=4 )605CALL Info( SolverName, ' ', Level=4 )606CALL Info( SolverName, 'Starting Assembly...', Level=4 )607608Ne = Mesh % NumberOfEdges609610!------------------------------------------------------------------------------611! lets start612!------------------------------------------------------------------------------613CALL DefaultInitialize()614body_id = -1615NULLIFY(Material)616!------------------------------------------------------------------------------617! Bulk elements (the sheet layer)618!------------------------------------------------------------------------------619DO t=1,Solver % NumberOfActiveElements620!------------------------------------------------------------------------------621! write some info on status of assembly622!------------------------------------------------------------------------------623IF ( RealTime() - at0 > 1.0 ) THEN624WRITE(Message,'(a,i3,a)' ) ' Assembly: ', INT(100.0 - 100.0 * &625(Solver % NumberOfActiveElements-t) / &626(1.0*Solver % NumberOfActiveElements)), ' % done'627628CALL Info( SolverName, Message, Level=5 )629at0 = RealTime()630END IF631!------------------------------------------------------------------------------632! Check if this element belongs to a body where scalar equation633! should be calculated634!------------------------------------------------------------------------------635Element => GetActiveElement(t,Solver)636637! cycle halo elements638!-------------------639IF (ParEnv % myPe /= Element % partIndex) CYCLE640641IF (.NOT.ASSOCIATED(Element)) CYCLE642IF ( Element % BodyId /= body_id ) THEN643body_id = Element % bodyId644Equation => GetEquation()645IF (.NOT.ASSOCIATED(Equation)) THEN646WRITE (Message,'(A,I3)') 'No Equation found for boundary element no. ', t647CALL FATAL(SolverName,Message)648END IF649650Material => GetMaterial()651IF (.NOT.ASSOCIATED(Material)) THEN652WRITE (Message,'(A,I3)') 'No Material found for boundary element no. ', t653CALL FATAL(SolverName,Message)654ELSE655material_id = GetMaterialId( Element, Found)656IF(.NOT.Found) THEN657WRITE (Message,'(A,I3)') 'No Material ID found for boundary element no. ', t658CALL FATAL(SolverName,Message)659END IF660END IF661END IF662!------------------------------------------------------------------------------663! Check if Element dimension is 2 (in x, y plane)664!------------------------------------------------------------------------------665dimSheet = Element % TYPE % DIMENSION666IF(dimSheet>2) THEN667WRITE (Message,'(A,I3)')' Work only for 1D or 2D elements'668CALL FATAL(SolverName,Message)669END IF670671!------------------------------------------------------------------------------672! Get element material parameters673!------------------------------------------------------------------------------674N = GetElementNOFNodes(Element)675CALL GetElementNodes( ElementNodes )676677!----------------------------------------------------678! Get parameters to compute the Total Conductivity679! K = k . f(h) . |grad HydPot |^(beta-2)680! k = SheetConductivity681! f(h) = h^alphas682!683! And read parameters to evaluate v - w684! w = u_b/l_r (h_r - h)685! N = p_i + rhow . g . b - HydPot686! v = A . h . | N |^{n-1} . N687!----------------------------------------------------688689CALL GetParametersSheet( Element, Material, N, SheetConductivity, alphas, &690betas, Ev, ub, Snn, lr, hr, Ar, ng )691692CT = 0.0_dp693Wopen = 0.0_dp694Phi0 = 0.0_dp695DO i=1, n696j = Element % NodeIndexes(i)697k = ThickPerm(j)698IF ( ASSOCIATED( ZbSol )) THEN699zb = ZbSolution(ZbPerm(j))700ELSE701IF (dimSheet==1) THEN702zb = Mesh % Nodes % y(j)703ELSE704zb = Mesh % Nodes % z(j)705END IF706END IF707CT(i) = Ev(i) /( WaterDensity * gravity)708Wopen(i) = MAX(ub(i) / lr(i) * (hr(i) - ThickSolution(k)), 0.0)709710Phi0(i) = Snn(i) + gravity*WaterDensity*zb711IF (.NOT. NeglectH) THEN712Phi0(i) = Phi0(i) + gravity*WaterDensity*ThickSolution(k)713END IF714END DO715716!------------------------------------------------------------------------------717! Add body forces718!------------------------------------------------------------------------------719LOAD = 0.0_dp720721BodyForce => GetBodyForce()722!IF ( ASSOCIATED( BodyForce ) ) THEN723! bf_id = GetBodyForceId()724! LOAD(1:N) = LOAD(1:N) + &725! GetReal( BodyForce, TRIM(Solver % Variable % Name) // ' Volume Source', Found )726!END IF727! f = m - w + v728! v is not added here as it will be linearized for the assembly729LOAD(1:N) = LOAD(1:N) - Wopen(1:N)730731!------------------------------------------------------------------------------732! Get element local matrices, and RHS vectors733!------------------------------------------------------------------------------734MASS = 0.0_dp735STIFF = 0.0_dp736FORCE = 0.0_dp737! cartesian coords738!----------------739IF ( CurrentCoordinateSystem() == Cartesian ) THEN740CALL SheetCompose( MASS, STIFF, FORCE, LOAD, &741ThickSolution(ThickPerm(Element % NodeIndexes(1:n))), &742HydPot(HydPotPerm(Element % NodeIndexes(1:N))), &743CT(1:N), SheetConductivity(1:n), alphas(1:n), betas(1:n), &744Phi0(1:n), Ar(1:n), ng(1:n), Element, N, ElementNodes )745ELSE746WRITE(Message,'(A)')' Work only for cartesian coordinate'747CALL FATAL( SolverName, Message)748END IF749!------------------------------------------------------------------------------750! If time dependent simulation add mass matrix to stiff matrix751!------------------------------------------------------------------------------752IF ( TransientSimulation ) THEN753CALL Default1stOrderTime( MASS, STIFF, FORCE )754END IF755756CALL DefaultUpdateEquations( STIFF, FORCE )757END DO ! Bulk elements758759! TODO: Is this really needed?760CALL DefaultFinishBulkAssembly()761762!------------------------------------------------------------------------------763! Edge element (the channels)764! Edge element are created in the SIF file using Element = 'n=1 e:1'765!------------------------------------------------------------------------------766! Go only if Activate Channels is True767!------------------------------------------------------------------------------768IF (Channels) THEN769body_id = -1770NULLIFY(Material)771DO t=1, Mesh % NumberOfEdges772Edge => Mesh % Edges(t)773IF (.NOT.ASSOCIATED(Edge)) CYCLE774IF (ParEnv % PEs > 1) THEN775IF (ParEnv % myPe /= Mesh % ParallelInfo % EdgeNeighbourList(t) % Neighbours(1)) CYCLE776END IF777n = Edge % TYPE % NumberOfNodes778779! Work only for 202 elements => n=2780IF (n/=2) CALL FATAL(SolverName, 'Work only for edge element of type 202')781! We keep only the edge which belong in the sheet surface782! i.e. Both nodes have Perm > 0783IF (ANY(HydPotPerm(Edge % NodeIndexes(1:n))==0)) CYCLE784785! We check if we are in a boundary where we want No channel786IF (ALL(NoChannel(Edge % NodeIndexes(1:n)))) CYCLE787788789EdgeNodes % x(1:n) = Mesh % Nodes % x(Edge % NodeIndexes(1:n))790EdgeNodes % y(1:n) = Mesh % Nodes % y(Edge % NodeIndexes(1:n))791EdgeNodes % z(1:n) = Mesh % Nodes % z(Edge % NodeIndexes(1:n))792793! Compute the unit tangent vector of the edge794EdgeTangent(1) = EdgeNodes % x(2) - EdgeNodes % x(1)795EdgeTangent(2) = EdgeNodes % y(2) - EdgeNodes % y(1)796EdgeTangent(3) = EdgeNodes % z(2) - EdgeNodes % z(1)797EdgeTangent = EdgeTangent / SQRT(SUM(EdgeTangent*EdgeTangent))798799! Read in the Material Section the needed Parameters for the800! Channels -801! In case DIM = 2 Faces have a material ASSOCIATED802! In case DIM = 3, Bulk elements have material ASSOCIATED803IF (DIM==2) THEN804Bulk => Edge % BoundaryInfo % Left805IF (.Not.ASSOCIATED(Bulk)) Bulk => Edge % BoundaryInfo % Right806ELSE807Face => Edge % BoundaryInfo % Left808IF (.Not.ASSOCIATED(Face)) Face => Edge % BoundaryInfo % Right809IF (ASSOCIATED(Face)) THEN810Bulk => Face % BoundaryInfo % Left811IF (.Not.ASSOCIATED(Bulk)) Bulk => Face % BoundaryInfo % Right812END IF813END IF814IF (.Not.ASSOCIATED(Bulk)) THEN815WRITE(Message,'(a,i0)')'No face or bulk element associated to edge no:', t816CALL FATAL(SolverName, Message)817END IF818819IF ( Bulk % BodyId /= body_id ) THEN820body_id = Bulk % bodyId821Material => GetMaterial( Bulk )822IF (.NOT.ASSOCIATED(Material)) THEN823WRITE (Message,'(A,I3)') 'No Material found for edge no. ', t824CALL FATAL(SolverName,Message)825ELSE826material_id = GetMaterialId( Bulk, Found)827IF(.NOT.Found) THEN828WRITE (Message,'(A,I3)') 'No Material ID found for edge no. ', t829CALL FATAL(SolverName,Message)830END IF831END IF832END IF833834!----------------------------------------------------835! Get parameters to compute the Total Conductivity836! Kc = kc . fc(h) . |grad HydPot |^(betac-2)837! kc = ChannelConductivity838! fc(h) = S^alphac839! and the closure and opening velocity for Channels840!----------------------------------------------------841842CALL GetParametersChannel( Edge, Material, n, SheetConductivity, &843ChannelConductivity, alphac, betac, alphas, betas, IceDensity, &844Snn, Ac, ng, CCt, CCw, lc )845846Phi0 = 0.0_dp847Afactor = 0.0_dp848Bfactor = 0.0_dp849Phi0 = 0.0_dp850Phim = 0.0_dp851DO i=1,N852IF ( ASSOCIATED( ZbSol )) THEN853zb = ZbSolution(ZbPerm(Edge % NodeIndexes(i)))854ELSE855IF (dimSheet==1) THEN856zb = Mesh % Nodes % y(Edge % NodeIndexes(i))857ELSE858zb = Mesh % Nodes % z(Edge % NodeIndexes(i))859END IF860END If861Phim(i) = gravity*WaterDensity*zb862Phi0(i) = Snn(i) + Phim(i)863IF (.NOT. NeglectH) THEN864k = ThickPerm(Edge % NodeIndexes(i))865Phi0(i) = Phi0(i) + gravity*WaterDensity*ThickSolution(k)866END IF867Afactor(i) = CCt(i) * CCw(i) * WaterDensity868Bfactor(i) = 1.0/(Lw * IceDensity(i))869IF (meltChannels) Bfactor(i) = Bfactor(i) - 1.0/(Lw * WaterDensity)870END DO871872! The variable Channel Area is defined on the edges only873IF (AreaPerm(M+t) <= 0) CYCLE874875!CHANGE876!To stabilise channels877IF(MABool) THEN878IF(AreaSolution(AreaPerm(M+t)) > MaxArea) AreaSolution(AreaPerm(M+t)) = MaxArea879END IF880ChannelArea = AreaSolution(AreaPerm(M+t))881882!------------------------------------------------------------------------------883! Get element local matrices, and RHS vectors884!------------------------------------------------------------------------------885MASS = 0.0_dp886STIFF = 0.0_dp887FORCE = 0.0_dp888! cartesian coords889!----------------890891IF ( CurrentCoordinateSystem() == Cartesian ) THEN892CALL ChannelCompose( MASS, STIFF, FORCE, &893ThickSolution(ThickPerm(Edge % NodeIndexes(1:n))), &894HydPot(HydPotPerm(Edge % NodeIndexes(1:N))), ChannelArea, &895ChannelConductivity, alphac, betac, Phi0, Phim, Ac, lc, ng, &896SheetConductivity, alphas, betas, Afactor, Bfactor, EdgeTangent, &897Edge, n, EdgeNodes )898ELSE899WRITE(Message,'(A)')' Work only for cartesian coordinate'900CALL FATAL( SolverName, Message)901END IF902!CHANGE903!To stop weird channel instability where some channels grow904!exponentially to stupid levels and eventually mess up whole905!mesh. Usually seems to be channels with limited fluxes that906!shouldn't do much, so resetting to 0 seems safest907!TODO Come up with a better way of fixing this908IF(MABool) THEN909k = AreaPerm(M+t)910IF(AreaSolution(k) > MaxArea) AreaSolution(k) = MaxArea911END IF912! This should be not needed as MASS = 0 here913IF ( TransientSimulation ) THEN914CALL Default1stOrderTime( MASS, STIFF, FORCE, Edge )915END IF916917CALL DefaultUpdateEquations( STIFF, FORCE, Edge )918END DO ! Edge elements919END IF920921922!------------------------------------------------------------------------------923! Neumann & Newton boundary conditions924!------------------------------------------------------------------------------925IF (HaveMoulinMask) THEN926DO t=1,Solver % NumberOfActiveElements927Element => GetActiveElement(t,Solver)928BodyForce => GetBodyForce()929Storage = GetLogical(BodyForce,'Moulin Storage', Found)930FORCE = 0.0931932! cycle halo elements933!-------------------934IF (ParEnv % myPe /= Element % partIndex) CYCLE935IF (Storage) THEN936MoulinArea(1:N) = ListGetReal( BodyForce, 'Moulin Area', N, Element % NodeIndexes, Found, &937UnfoundFatal = .TRUE. )938END IF939DO I=1,Element % TYPE % NumberOfNodes940k = Element % NodeIndexes(I)941IF (MoulinMask(MoulinMaskPerm(k)) > 0.0) THEN942FORCE(i) = MoulinFluxVal(MoulinFluxPerm(k))943MASS(i,i) = MoulinArea(i)/(WaterDensity*gravity)944END IF945END DO946947IF ( TransientSimulation ) THEN948CALL Default1stOrderTime( MASS, STIFF, FORCE, Element )949END IF950CALL DefaultUpdateEquations( STIFF, FORCE, Element )951END DO952ELSE953DO t=1, Mesh % NumberOfBoundaryElements954! get element information955Element => GetBoundaryElement(t)956! Don't test this as BC on the side are not active for this solver957!IF ( .NOT.ActiveBoundaryElement() ) CYCLE958959! cycle halo elements960!--------------------961IF (ParEnv % myPe /= Element % partIndex) CYCLE962n = GetElementNOFNodes()963964IF ( GetElementFamily() == 1 ) THEN965! Moulin case (flux at node)966967BC => GetBC( Element )968bc_id = GetBCId( Element )969CALL GetElementNodes( ElementNodes )970IF ( ASSOCIATED( BC ) ) THEN971STIFF=0.0_dp972FORCE=0.0_dp973MASS=0.0_dp974975Storage = GetLogical(BC,'Moulin Storage', Found)976IF (Storage) THEN977MoulinArea(1:N) = ListGetReal( BC, 'Moulin Area', N, Element % NodeIndexes, Found, &978UnfoundFatal = .TRUE. )979! MASS is a scalar here980MASS(1,1) = MoulinArea(1)/(WaterDensity*gravity)981END IF982983! Is there surface input984MoulinFlux(1:N) = ListGetReal( BC, 'Moulin Flux', N, Element % NodeIndexes, Found, &985UnfoundFatal = .TRUE. )986FORCE(1) = MoulinFlux(1)987988! If variable exist, update the Flux from Moulins variable989IF (ASSOCIATED( QmSol )) THEN990j = Element % NodeIndexes(1)991IF ( HydPotPerm(j) <= 0 ) CYCLE992QmSolution(QmPerm(j)) = MoulinFlux(1)-MASS(1,1)*(HydPot(HydPotPerm(j))-HydPotPrev(HydPotPerm(j),1))/dt993END IF994995!------------------------------------------------------------------------------996! Update global matrices from local matrices997!------------------------------------------------------------------------------998! We need the Default1stOrderUpdate for the Moulin999IF ( TransientSimulation ) THEN1000CALL Default1stOrderTime( MASS, STIFF, FORCE, Element )1001END IF1002CALL DefaultUpdateEquations( STIFF, FORCE, Element )1003END IF10041005ELSE1006! flux over the sheet1007BC => GetBC( Element )1008bc_id = GetBCId( Element )1009CALL GetElementNodes( ElementNodes )10101011IF ( ASSOCIATED( BC ) ) THEN1012! Check that we are on the correct boundary part!1013STIFF=0.0_dp1014FORCE=0.0_dp1015MASS=0.0_dp10161017LOAD=0.0_dp1018FluxBC = GetLogical(BC,TRIM(Solver % Variable % Name) // ' Flux BC', Found)10191020IF (FluxBC) THEN1021!---------------1022!BC: -k@T/@n = q1023!---------------1024LOAD(1:N) = LOAD(1:N) + &1025GetReal( BC, 'Sheet Water Flux', Found )1026! -------------------------------------1027! set boundary due to coordinate system1028! -------------------------------------1029IF ( CurrentCoordinateSystem() == Cartesian ) THEN1030CALL SheetBoundary( STIFF, FORCE, &1031LOAD, Element, n, ElementNodes )1032ELSE1033WRITE(Message,'(A)')' Work only for cartesian coordinate'1034CALL FATAL( SolverName, Message)1035END IF1036!!! TODO : do we need that as MASS = 0 here1037IF ( TransientSimulation ) THEN1038CALL Default1stOrderTime( MASS, STIFF, FORCE, Element )1039END IF10401041CALL DefaultUpdateEquations( STIFF, FORCE, Element )1042END IF1043END IF1044END IF ! Element 1O110451046END DO ! Neumann & Newton BCs1047END IF ! END IF(HaveMoulinMask)1048!------------------------------------------------------------------------------10491050CALL DefaultFinishAssembly()10511052CALL DefaultDirichletBCs()10531054CALL Info( SolverName, 'Assembly done', Level=4 )10551056!------------------------------------------------------------------------------1057! Solve the system and check for convergence1058!------------------------------------------------------------------------------1059at = CPUTime() - at1060st = CPUTime()1061Norm = 0.0_dp10621063PrevNorm = Solver % Variable % Norm1064Norm = DefaultSolve()10651066st = CPUTime()-st1067totat = totat + at1068totst = totst + st1069WRITE(Message,'(a,i4,a,F8.2,F8.2)') 'iter: ',iter,' Assembly: (s)', at, totat1070CALL Info( SolverName, Message, Level=4 )1071WRITE(Message,'(a,i4,a,F8.2,F8.2)') 'iter: ',iter,' Solve: (s)', st, totst1072CALL Info( SolverName, Message, Level=4 )107310741075IF ( PrevNorm + Norm /= 0.0d0 ) THEN1076RelativeChange = 2.0d0 * ABS( PrevNorm-Norm ) / (PrevNorm + Norm)1077ELSE1078RelativeChange = 0.0d01079END IF10801081WRITE( Message, * ) 'Result Norm : ',Norm1082CALL Info( SolverName, Message, Level=4 )1083WRITE( Message, * ) 'Relative Change : ',RelativeChange1084CALL Info( SolverName, Message, Level=4 )10851086! !----------------------1087! ! check for convergence1088! !----------------------1089IF ( RelativeChange < NonlinearTol .AND. iter > NonlinearIterMin ) EXIT1090END DO ! of the nonlinear iteration10911092!------------------------------------------------------------------------------1093! Update the Sheet Thickness1094!------------------------------------------------------------------------------1095DO t=1,Solver % NumberOfActiveElements1096Element => GetActiveElement(t,Solver)1097IF (ParEnv % myPe /= Element % partIndex) CYCLE10981099IF (.NOT.ASSOCIATED(Element)) CYCLE1100IF ( Element % BodyId /= body_id ) THEN1101body_id = Element % bodyId1102Equation => GetEquation()1103IF (.NOT.ASSOCIATED(Equation)) THEN1104WRITE (Message,'(A,I3)') 'No Equation found for boundary element no. ', t1105CALL FATAL(SolverName,Message)1106END IF11071108Material => GetMaterial()1109IF (.NOT.ASSOCIATED(Material)) THEN1110WRITE (Message,'(A,I3)') 'No Material found for boundary element no. ', t1111CALL FATAL(SolverName,Message)1112ELSE1113material_id = GetMaterialId( Element, Found)1114IF(.NOT.Found) THEN1115WRITE (Message,'(A,I3)') 'No Material ID found for boundary element no. ', t1116CALL FATAL(SolverName,Message)1117END IF1118END IF1119END IF11201121N = GetElementNOFNodes(Element)1122CALL GetElementNodes( ElementNodes )11231124!CHANGE1125!If calving, cycle elements with ungrounded nodes and zero all1126!hydrology variables1127IF(Calving) THEN1128CycleElement = .FALSE.1129WorkVar => VariableGet(Mesh % Variables, "gmcheck", ThisOnly=.TRUE., UnfoundFatal=.FALSE.)1130WorkVar2 => VariableGet(Mesh % Variables, "groundedmask", ThisOnly=.TRUE., UnfoundFatal=.FALSE.)1131IF(ASSOCIATED(WorkVar)) THEN1132DO i=1, N1133IF(WorkVar % Values(WorkVar % Perm(Element % NodeIndexes(i)))>0.0) THEN1134!IF(WorkVar2 % Values(WorkVar2 % Perm(Element % NodeIndexes(i)))<0.0) THEN1135CycleElement = .TRUE.11361137WSolution(WPerm(Element % NodeIndexes(i))) = 0.01138Vvar(Element % NodeIndexes(i)) = 0.01139NSolution(NPerm(Element % NodeIndexes(i))) = 0.01140!PwSolution(PwPerm(Element % NodeIndexes(i))) = 0.01141hstoreSolution(hstorePerm(Element % NodeIndexes(i))) = 0.01142!END IF1143END IF1144END DO1145END IF1146IF(ASSOCIATED(WorkVar2) .AND. .NOT. ASSOCIATED(WorkVar)) THEN1147DO i=1, N1148IF(WorkVar2 % Values(WorkVar2 % Perm(Element % NodeIndexes(i)))<0.0) THEN1149CycleElement = .TRUE.11501151WSolution(WPerm(Element % NodeIndexes(i))) = 0.01152Vvar(Element % NodeIndexes(i)) = 0.01153NSolution(NPerm(Element % NodeIndexes(i))) = 0.01154!PwSolution(PwPerm(Element % NodeIndexes(i))) = 0.01155hstoreSolution(hstorePerm(Element % NodeIndexes(i))) = 0.01156END IF1157END DO1158END IF1159NULLIFY(WorkVar, WorkVar2)1160IF(CycleElement) CYCLE1161END IF11621163CALL GetParametersSheet( Element, Material, N, SheetConductivity, alphas, &1164betas, Ev, ub, Snn, lr, hr, Ar, ng )11651166CT = 0.0_dp1167Phi0 = 0.0_dp1168DO i=1, n1169j = Element % NodeIndexes(i)1170k = ThickPerm(j)1171IF ( ASSOCIATED( ZbSol )) THEN1172zb = ZbSolution(ZbPerm(j))1173ELSE1174IF (dimSheet==1) THEN1175zb = Mesh % Nodes % y(j)1176ELSE1177zb = Mesh % Nodes % z(j)1178END IF1179END IF1180CT(i) = Ev(i) /( WaterDensity * gravity)1181Phi0(i) = Snn(i) + gravity*WaterDensity*zb1182Np = Phi0(i) - HydPot(HydPotPerm(j))1183IF (.NOT. NeglectH) THEN1184Np = Np + WaterDensity*gravity*ThickSolution(k)1185END IF1186pw = HydPot(HydPotPerm(j)) - gravity*WaterDensity*zb1187Vvar(j) = Ar(i)*ABS(Np)**(ng(i)-1.0)*Np1188Wo = MAX(ub(i) / lr(i) * (hr(i) - ThickSolution(k)), 0.0)1189he = Ev(i)*(HydPot(HydPotPerm(j))/(WaterDensity*gravity)-zb)1190ublr(j) = ub(i)/lr(i)1191hr2(j) = hr(i)11921193!CHANGE1194!To stop it working out values for non-ice covered parts of a1195!hydromesh in a coupled calving-hydro simulation1196IF(Calving) THEN1197IF(Snn(i)==0.0) THEN1198Np = 0.01199!pw = 0.01200he = 0.01201END IF1202END IF12031204! Save this for output only1205IF (ASSOCIATED(WSol)) WSolution(WPerm(j)) = Wo1206IF (ASSOCIATED(NSol)) NSolution(NPerm(j)) = Np1207IF (ASSOCIATED(PwSol)) PwSolution(PwPerm(j)) = pw1208IF (ASSOCIATED(hstoreSol)) hstoreSolution(hstorePerm(j)) = he1209END DO1210END DO ! Bulk elements1211! Loop over all nodes to update ThickSolution1212DO j = 1, Mesh % NumberOfNodes1213k = ThickPerm(j)1214IF (k==0) CYCLE1215!CHANGE1216!If calving, cycle elements with ungrounded nodes and zero all1217!hydrology variables1218IF(Calving) THEN1219CycleElement = .FALSE.1220WorkVar => VariableGet(Mesh % Variables, "gmcheck", ThisOnly=.TRUE., UnfoundFatal=.FALSE.)1221WorkVar2 => VariableGet(Mesh % Variables, "groundedmask", ThisOnly=.TRUE., UnfoundFatal=.FALSE.)1222IF(ASSOCIATED(WorkVar)) THEN1223IF(WorkVar % Values(k)>0.0) THEN !.AND. WorkVar2 % Values(k)<0.0) THEN1224CycleElement = .TRUE.1225ThickSolution(k) = 0.01226ThickPrev(k,1) = 0.01227END IF1228END IF1229IF(ASSOCIATED(WorkVar2) .AND. .NOT. ASSOCIATED(WorkVar)) THEN1230IF(WorkVar2 % Values(k)<0.0) THEN1231CycleElement = .TRUE.1232ThickSolution(k) = 0.01233ThickPrev(k,1) = 0.01234END IF1235END IF1236WorkVar => VariableGet(Mesh % Variables, "hydraulic potential", ThisOnly=.TRUE., UnfoundFatal=.FALSE.)1237IF(WorkVar % Values(k)==0.0) THEN1238ThickSolution(k) = 0.01239ThickPrev(k,1) = 0.01240CycleElement = .TRUE.1241END IF1242NULLIFY(WorkVar, WorkVar2)1243IF(CycleElement) CYCLE1244END IF12451246IF(MaxHBool) THEN1247IF (ThickSolution(k)>MaxH) THEN1248ThickSolution(k) = MaxH1249!ThickPrev(k,1) = 0.01250END IF1251END IF12521253IF(MinHBool) THEN1254IF (ThickSolution(k)<MinH) THEN1255ThickSolution(k) = MinH1256END IF1257END IF12581259SELECT CASE(methodSheet)1260CASE('implicit')1261IF (ThickSolution(k) > hr2(j)) THEN1262ThickSolution(k) = MAX(ThickPrev(k,1)/(1.0_dp + dt*Vvar(j)) , AEPS)1263ELSE1264ThickSolution(k) = MAX((ThickPrev(k,1) + dt*ublr(j)*hr2(j))/(1.0_dp + dt*(Vvar(j)+ublr(j))) , AEPS)1265END IF1266CASE('explicit')1267IF (ThickSolution(k) > hr2(j)) THEN1268ThickSolution(k) = MAX(ThickPrev(k,1)*(1.0_dp - dt*Vvar(j)) , AEPS)1269ELSE1270ThickSolution(k) = MAX(ThickPrev(k,1)*(1.0_dp - dt*(Vvar(j)+ublr(j))) + dt*ublr(j)*hr2(j), AEPS)1271END IF1272CASE('crank-nicolson')1273IF (ThickSolution(k) > hr2(j)) THEN1274ThickSolution(k) = MAX(ThickPrev(k,1)*(1.0_dp - 0.5*dt*Vvar(j))/(1.0_dp + 0.5*dt*Vvar(j)) , AEPS)1275ELSE1276ThickSolution(k) = MAX((ThickPrev(k,1)*(1.0_dp - 0.5*dt*(Vvar(j)+ublr(j))) + dt*ublr(j)*hr2(j)) &1277& /(1.0_dp + 0.5*dt*(Vvar(j)+ublr(j))) , AEPS)1278END IF1279END SELECT12801281! Update Vvar1282Vvar(j) = Vvar(j) * ThickSolution(k)12831284END DO1285!------------------------------------------------------------------------------1286! Update the Channels Area1287!------------------------------------------------------------------------------1288!------------------------------------------------------------------------------1289! non-linear system iteration loop1290!------------------------------------------------------------------------------1291IF (Channels) THEN1292PrevNorm = ChannelAreaNorm()12931294DO iter = 1, NonlinearIter1295DO t=1, Mesh % NumberOfEdges1296Edge => Mesh % Edges(t)1297IF (.NOT.ASSOCIATED(Edge)) CYCLE1298IF (ParEnv % PEs > 1) THEN1299IF (ParEnv % myPe /= Mesh % ParallelInfo % EdgeNeighbourList(t) % Neighbours(1)) CYCLE1300END IF1301n = Edge % TYPE % NumberOfNodes1302IF (ANY(HydPotPerm(Edge % NodeIndexes(1:n))==0)) CYCLE1303IF (ALL(NoChannel(Edge % NodeIndexes(1:n)))) CYCLE13041305!CHANGE1306!If calving, cycle elements with ungrounded nodes and zero all1307!hydrology variables1308IF(Calving) THEN1309CycleElement = .FALSE.1310WorkVar => VariableGet(Mesh % Variables, "gmcheck", ThisOnly=.TRUE., UnfoundFatal=.FALSE.)1311WorkVar2 => VariableGet(Mesh % Variables, "groundedmask", ThisOnly=.TRUE., UnfoundFatal=.FALSE.)1312IF(ASSOCIATED(WorkVar)) THEN1313DO i=1, n1314IF(WorkVar % Values(WorkVar % Perm(Edge % NodeIndexes(i)))>0.0) THEN1315!IF(WorkVar2 % Values(WorkVar2 % Perm(Edge % NodeIndexes(i)))<0.0) THEN1316CycleElement = .TRUE.1317AreaSolution(AreaPerm(M+t)) = 0.01318QcSolution(QcPerm(M+t)) = 0.01319!END IF1320END IF1321END DO1322END IF1323IF(ASSOCIATED(WorkVar2) .AND. .NOT. ASSOCIATED(WorkVar)) THEN1324DO i=1,n1325IF(WorkVar2 % Values(WorkVar2 % Perm(Edge % NodeIndexes(i)))<0.0) THEN1326CycleElement = .TRUE.1327AreaSolution(AreaPerm(M+t)) = 0.01328QcSolution(QcPerm(M+t)) = 0.01329END IF1330END DO1331END IF1332NULLIFY(WorkVar, WorkVar2)1333IF(CycleElement) CYCLE1334END IF13351336EdgeNodes % x(1:n) = Mesh % Nodes % x(Edge % NodeIndexes(1:n))1337EdgeNodes % y(1:n) = Mesh % Nodes % y(Edge % NodeIndexes(1:n))1338EdgeNodes % z(1:n) = Mesh % Nodes % z(Edge % NodeIndexes(1:n))13391340! Compute the unit tangent vector of the edge1341EdgeTangent(1) = EdgeNodes % x(2) - EdgeNodes % x(1)1342EdgeTangent(2) = EdgeNodes % y(2) - EdgeNodes % y(1)1343EdgeTangent(3) = EdgeNodes % z(2) - EdgeNodes % z(1)1344EdgeTangent = EdgeTangent / SQRT(SUM(EdgeTangent*EdgeTangent))13451346IF (DIM==2) THEN1347Bulk => Edge % BoundaryInfo % Left1348IF (.Not.ASSOCIATED(Bulk)) Bulk => Edge % BoundaryInfo % Right1349ELSE1350Face => Edge % BoundaryInfo % Left1351IF (.Not.ASSOCIATED(Face)) Face => Edge % BoundaryInfo % Right1352IF (ASSOCIATED(Face)) THEN1353Bulk => Face % BoundaryInfo % Left1354IF (.Not.ASSOCIATED(Bulk)) Bulk => Face % BoundaryInfo % Right1355END IF1356END IF1357IF (.Not.ASSOCIATED(Bulk)) THEN1358WRITE(Message,'(a,i0)')'No face or bulk element associated to edge no:', t1359CALL FATAL(SolverName, Message)1360END IF13611362IF ( Bulk % BodyId /= body_id ) THEN1363body_id = Bulk % bodyId1364Material => GetMaterial( Bulk )1365IF (.NOT.ASSOCIATED(Material)) THEN1366WRITE (Message,'(A,I3)') 'No Material found for edge no. ', t1367CALL FATAL(SolverName,Message)1368ELSE1369material_id = GetMaterialId( Bulk, Found)1370IF(.NOT.Found) THEN1371WRITE (Message,'(A,I3)') 'No Material ID found for edge no. ', t1372CALL FATAL(SolverName,Message)1373END IF1374END IF1375END IF13761377CALL GetParametersChannel( Edge, Material, n, SheetConductivity, &1378ChannelConductivity, alphac, betac, alphas, betas, IceDensity, &1379Snn, Ac, ng, CCt, CCw, lc )13801381Phi0 = 0.0_dp1382Afactor = 0.0_dp1383Bfactor = 0.0_dp1384Phi0 = 0.0_dp1385Phim = 0.0_dp1386DO i=1,N1387IF ( ASSOCIATED( ZbSol )) THEN1388zb = ZbSolution(ZbPerm(Edge % NodeIndexes(i)))1389ELSE1390IF (dimSheet==1) THEN1391zb = Mesh % Nodes % y(Edge % NodeIndexes(i))1392ELSE1393zb = Mesh % Nodes % z(Edge % NodeIndexes(i))1394END IF1395END If1396Phim(i) = gravity*WaterDensity*zb1397Phi0(i) = Snn(i) + Phim(i)1398IF (.NOT. NeglectH) THEN1399k = ThickPerm(Edge % NodeIndexes(i))1400Phi0(i) = Phi0(i) + gravity*WaterDensity*ThickSolution(k)1401END IF1402Afactor(i) = CCt(i) * CCw(i) * WaterDensity1403Bfactor(i) = 1.0/(Lw * IceDensity(i))1404END DO14051406IF (AreaPerm(M+t) <= 0 ) CYCLE1407ChannelArea = AreaSolution(AreaPerm(M+t))14081409! Compute the force term to evolve the channels area1410! Equation of the form dS/dt = S x ALPHA + BETA1411IF ( CurrentCoordinateSystem() == Cartesian ) THEN1412CALL GetEvolveChannel( ALPHA, BETA, Qc, ChannelArea, &1413HydPot(HydPotPerm(Edge % NodeIndexes(1:n))), &1414ThickSolution(ThickPerm(Edge % NodeIndexes(1:n))), &1415alphac, betac, ChannelConductivity, Phi0, Phim, Ac, lc, ng, &1416SheetConductivity, alphas, betas, Afactor, Bfactor, &1417EdgeTangent, Edge, n, EdgeNodes, LimitEffPres)1418ELSE1419WRITE(Message,'(A)')' Work only for cartesian coordinate'1420CALL FATAL( SolverName, Message)1421END IF14221423k = AreaPerm(M+t)1424IF ( k <= 0 ) CYCLE1425SELECT CASE(methodChannels)1426CASE('implicit')1427AreaSolution(k) = (AreaPrev(k,1) + dt*BETA)/(1.0_dp - dt*ALPHA)1428CASE('explicit')1429AreaSolution(k) = AreaPrev(k,1)*(1.0_dp + dt*ALPHA) + dt*BETA1430CASE('crank-nicolson')1431AreaSolution(k) = (AreaPrev(k,1)*(1.0_dp + 0.5*ALPHA*dt) + dt*BETA)/(1.0_dp - 0.5*dt*ALPHA)1432END SELECT14331434!CHANGE1435!To stop weird channel instability where some channels grow1436!exponentially to stupid levels and eventually mess up whole1437!mesh. Usually seems to be channels with limited fluxes that1438!shouldn't do much, so resetting to 0 seems safest1439!TODO Come up with a better way of fixing this1440IF(MABool) THEN1441IF(AreaPrev(k,1) /= 0.0) THEN1442IF(AreaSolution(k)>1.0 .AND. (AreaSolution(k)/AreaPrev(k,1))>5.0) THEN1443AreaSolution(k) = AreaPrev(k,1)1444END IF1445END IF1446IF(AreaSolution(k) > MaxArea) AreaSolution(k) = MaxArea1447IF(ISNAN(AreaSolution(k))) AreaSolution(k) = 0.01448END IF14491450! Save Qc if variable exists1451IF (ASSOCIATED(QcSol)) THEN1452IF ( QcPerm(M+t) <= 0 ) CYCLE1453QcSolution(QcPerm(M+t)) = Qc1454END IF1455END DO14561457Norm = ChannelAreaNorm()1458t = Mesh % NumberOfEdges14591460IF ( PrevNorm + Norm /= 0.0d0 ) THEN1461RelativeChange = 2.0d0 * ABS( PrevNorm-Norm ) / (PrevNorm + Norm)1462ELSE1463RelativeChange = 0.0d01464END IF14651466PrevNorm = Norm14671468WRITE( Message, * ) 'S (NRM,RELC) : ',iter, Norm, RelativeChange1469CALL Info( SolverName, Message, Level=3 )14701471! !----------------------1472! ! check for convergence1473! !----------------------1474IF ( RelativeChange < NonlinearTol .AND. iter > NonlinearIterMin ) EXIT1475END DO ! of the nonlinear iteration147614771478! Make sure Area >= 01479AreaSolution = MAX(AreaSolution, 0.0_dp )14801481!Stop channels from expanding to eleventy-stupid1482IF(MABool) THEN1483AreaSolution = MIN(AreaSolution,MaxArea)1484END IF14851486END IF ! If Channels14871488! Check for convergence1489CoupledNorm = ComputeNorm( Solver, SIZE(HydPot), HydPot )14901491IF ( PrevCoupledNorm + CoupledNorm /= 0.0d0 ) THEN1492RelativeChange = 2.0d0 * ABS( PrevCoupledNorm-CoupledNorm ) / (PrevCoupledNorm + CoupledNorm)1493ELSE1494RelativeChange = 0.0d01495END IF1496PrevCoupledNorm = CoupledNorm14971498WRITE( Message, * ) 'COUPLING LOOP (NRM,RELC) : ',iterC, CoupledNorm, RelativeChange1499CALL Info( SolverName, Message, Level=3 )15001501IF ((RelativeChange < CoupledTol).AND. (iterC > 1)) EXIT1502END DO ! iterC15031504!--------------------------------------------------------------------------------------------1505! Output some useful variables1506!--------------------------------------------------------------------------------------------1507! Already computed variables1508IF (ASSOCIATED(VSol)) THEN1509DO i=1, M1510IF ( VPerm(i) <= 0 ) CYCLE1511VSolution(VPerm(i)) = Vvar(i)1512ENDDO1513ENDIF15141515! Output the sheet discharge (dimension = dimSheet)1516IF (ASSOCIATED(qSol)) THEN1517Refq = 0.0_dp1518qSolution = 0.0_dp15191520! Loop over all elements are we need to compute grad(Phi)1521DO t=1,Solver % NumberOfActiveElements1522!CHANGE - necessary if using a 2D mesh as is otherwise set to 1 as1523!boundary elements are last in first loop where it's set1524dimSheet = Element % TYPE % DIMENSION1525Element => GetActiveElement(t,Solver)1526IF (ParEnv % myPe /= Element % partIndex) CYCLE15271528IF (.NOT.ASSOCIATED(Element)) CYCLE1529IF ( Element % BodyId /= body_id ) THEN1530body_id = Element % bodyId1531Material => GetMaterial()1532IF (.NOT.ASSOCIATED(Material)) THEN1533WRITE (Message,'(A,I3)') 'No Material found for boundary element no. ', t1534CALL FATAL(SolverName,Message)1535ELSE1536material_id = GetMaterialId( Element, Found)1537IF(.NOT.Found) THEN1538WRITE (Message,'(A,I3)') 'No Material ID found for boundary element no. ', t1539CALL FATAL(SolverName,Message)1540END IF1541END IF1542END IF15431544n = GetElementNOFNodes(Element)1545CALL GetElementNodes( ElementNodes )1546!If calving, cycle elements with ungrounded nodes and zero all1547!hydrology variables1548IF(Calving) THEN1549CycleElement = .FALSE.1550WorkVar => VariableGet(Mesh % Variables, "gmcheck", ThisOnly=.TRUE., UnfoundFatal=.FALSE.)1551WorkVar2 => VariableGet(Mesh % Variables, "groundedmask", ThisOnly=.TRUE., UnfoundFatal=.FALSE.)1552IF(ASSOCIATED(WorkVar)) THEN1553DO i=1, n1554IF(WorkVar % Values(WorkVar % Perm(Element % NodeIndexes(i)))>0.0) THEN1555!IF(WorkVar2 % Values(WorkVar2 % Perm(Element % NodeIndexes(i)))<0.0) THEN1556CycleElement = .TRUE.1557DO j=1,dimSheet1558k = dimSheet*(qPerm(Element % NodeIndexes(i))-1)+j1559qSolution(k) = 0.01560Refq(k) = 0.01561END DO1562EXIT1563!END IF1564END IF1565END DO1566END IF1567IF(ASSOCIATED(WorkVar2) .AND. .NOT. ASSOCIATED(WorkVar)) THEN1568DO i=1,n1569IF(WorkVar2 % Values(WorkVar2 % Perm(Element % NodeIndexes(i)))<0.0) THEN1570CycleElement = .TRUE.1571DO j=1,dimSheet1572k = dimSheet*(qPerm(Element % NodeIndexes(i))-1)+j1573qSolution(k) = 0.01574Refq(k) = 0.01575END DO1576EXIT1577END IF1578END DO1579END IF1580NULLIFY(WorkVar, WorkVar2)1581IF(CycleElement) CYCLE1582END IF15831584! we need the SheetConductivity, alphas, betas1585CALL GetParametersSheet( Element, Material, n, SheetConductivity, alphas, &1586betas, Ev, ub, Snn, lr, hr, Ar, ng )15871588! Go for all nodes of the element1589DO i=1,n1590Discharge = 0.0_dp1591CALL SheetDischargeCompute( &1592HydPot(HydPotPerm(Element % NodeIndexes(1:n))), &1593ThickSolution(ThickPerm(Element % NodeIndexes(i))), &1594SheetConductivity(i), alphas(i), betas(i), &1595Discharge, Element, n, ElementNodes, i )15961597! One more value for that node1598DO j=1,dimSheet1599k = dimSheet*(qPerm(Element % NodeIndexes(i))-1)+j1600Refq(k) = Refq(k) + 1.0_dp1601qSolution(k) = qSolution(k) + Discharge(j)1602END DO1603END DO1604END DO16051606! Mean nodal value1607DO i=1,n1608DO j=1,dimSheet1609k = dimSheet*(qPerm(Element % NodeIndexes(i))-1)+j1610IF ( Refq(k) > 0.0_dp ) THEN1611qSolution(k) = qSolution(k)/Refq(k)1612END IF1613END DO1614END DO16151616END IF16171618SubroutineVisited = .TRUE.16191620!CHANGE - to make sure PrevValues for added variables in calving updated1621IF(Calving) THEN1622WorkVar => VariableGet(Mesh % Variables, 'Sheet Thickness',ThisOnly=.TRUE.)1623WorkVar % PrevValues(:,1) = WorkVar % Values1624WorkVar => VariableGet(Mesh % Variables, 'Channel Area',ThisOnly=.TRUE.)1625WorkVar % PrevValues(:,1) = WorkVar % Values1626NULLIFY(WorkVar)1627END IF16281629CONTAINS16301631! Compute consistent channel norm only considering the edges that also have hydrology defined on the nodes.1632! In parallel only consider the edges in the partition where it is active.1633!----------------------------------------------------------------------------------------------------------1634FUNCTION ChannelAreaNorm() RESULT( AreaNorm )1635REAL(KIND=dp) :: AreaNorm16361637INTEGER :: t,i,n,ncount1638REAL(KIND=dp) :: s, ssum, smin, smax16391640ncount = 01641smin = HUGE(smin)1642smax = 0.0_dp1643ssum = 0.0_dp16441645DO t=1, Mesh % NumberOfEdges1646i = AreaPerm(Mesh % NumberOfNodes + t)1647IF(i==0) CYCLE16481649IF (ParEnv % PEs > 1) THEN1650IF (ParEnv % myPe /= Mesh % ParallelInfo % EdgeNeighbourList(t) % Neighbours(1)) CYCLE1651END IF16521653Edge => Mesh % Edges(t)1654n = Edge % TYPE % NumberOfNodes1655IF (ANY(HydPotPerm(Edge % NodeIndexes(1:n))==0)) CYCLE16561657s = AreaSolution(i)16581659ssum = ssum + s**21660smin = MIN(smin,s)1661smax = MAX(smax,s)16621663ncount = ncount + 11664END DO16651666IF( ParEnv % PEs > 1 ) THEN1667ncount = ParallelReduction(ncount)1668ssum = ParallelReduction(ssum)1669smin = ParallelReduction(smin,1)1670smax = ParallelReduction(smax,2)1671END IF16721673WRITE( Message, * ) 'Range S:', smin, smax1674CALL Info( SolverName, Message )16751676AreaNorm = SQRT(ssum / ncount )16771678END FUNCTION ChannelAreaNorm1679168016811682!------------------------------------------------------------------------------1683SUBROUTINE SheetCompose( MassMatrix, StiffMatrix, ForceVector, &1684LoadVector, NodalH, NodalHydPot, NodalCT, NodalC2, Nodalalphas, &1685Nodalbetas, NodalPhi0, NodalAr, NodalNg, Element, n, Nodes )1686!------------------------------------------------------------------------------1687USE MaterialModels1688USE Integration1689USE Differentials16901691IMPLICIT NONE16921693REAL(KIND=dp), DIMENSION(:) :: ForceVector, LoadVector1694REAL(KIND=dp), DIMENSION(:,:) :: MassMatrix, StiffMatrix1695REAL(KIND=dp) :: NodalHydPot(:), Nodalalphas(:), NodalBetas(:)1696REAL(KIND=dp) :: NodalCT(:), NodalC2(:), NodalH(:)1697REAL(KIND=dp) :: NodalPhi0(:), NodalAr(:), NodalNg(:)16981699INTEGER :: n17001701TYPE(Nodes_t) :: Nodes1702TYPE(Element_t), POINTER :: Element17031704!------------------------------------------------------------------------------1705! Local variables1706!------------------------------------------------------------------------------17071708REAL(KIND=dp) :: dBasisdx(n,3), detJ1709REAL(KIND=dp) :: Basis(n)17101711REAL(KIND=dp) :: Force17121713REAL(KIND=dp) :: A, M1714REAL(KIND=dp) :: Load17151716REAL(KIND=dp) :: x, y, z17171718INTEGER :: i, j, k, c, p, q, t, dim, N_Integ, NBasis17191720TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff1721REAL(KIND=dp) :: s, u, v, w17221723REAL(KIND=dp) :: CT, C2, Vfactor, Phi0, PhiG, LoadAtIP17241725REAL(KIND=dp) :: gradPhi(3), Ngrad, na, nb, ng, hsheet17261727REAL(KIND=dp), DIMENSION(:), POINTER :: U_Integ, V_Integ, W_Integ, S_Integ17281729LOGICAL :: Found, Transient, stat17301731!------------------------------------------------------------------------------17321733Transient = GetString(GetSimulation(),'Simulation type',Found)=='transient'17341735! the dimension for this solver is given by the dim of the elements1736! should work for 2D and 3D general problem1737dim = Element % TYPE % DIMENSION17381739ForceVector = 0.0_dp1740StiffMatrix = 0.0_dp1741MassMatrix = 0.0_dp1742Load = 0.0_dp17431744NBasis = n17451746!------------------------------------------------------------------------------1747! Integration stuff1748!------------------------------------------------------------------------------1749IntegStuff = GaussPoints( element )1750U_Integ => IntegStuff % u1751V_Integ => IntegStuff % v1752W_Integ => IntegStuff % w1753S_Integ => IntegStuff % s1754N_Integ = IntegStuff % n175517561757!------------------------------------------------------------------------------1758! Now we start integrating1759!------------------------------------------------------------------------------1760DO t=1,N_Integ1761u = U_Integ(t)1762v = V_Integ(t)1763w = W_Integ(t)17641765!------------------------------------------------------------------------------1766! Basis function values & derivatives at the integration point1767!------------------------------------------------------------------------------1768stat = ElementInfo( Element, Nodes, u, v, w, detJ, Basis, dBasisdx )17691770s = detJ * S_Integ(t)1771!------------------------------------------------------------------------------1772! Coefficient of derivative terms1773! at the integration point1774!------------------------------------------------------------------------------1775CT = SUM( NodalCT(1:n) * Basis(1:n) )1776!------------------------------------------------------------------------------1777! Compute the effective conductivity K = k h |grad Phi|^{beta-2}1778!------------------------------------------------------------------------------1779gradPhi = 0.0_dp1780Do i=1, dim1781gradPhi(i) = SUM(NodalHydPot(1:n)*dBasisdx(1:n,i))1782END DO1783Ngrad = SQRT(SUM(gradPhi(1:dim)*gradPhi(1:dim)))1784IF (Ngrad < AEPS) Ngrad = AEPS1785nb = SUM(NodalBetas(1:n)*Basis(1:n))1786na = SUM(Nodalalphas(1:n)*Basis(1:n))1787hsheet = SUM(NodalH(1:n)*Basis(1:n))17881789C2 = SUM( NodalC2(1:n) * Basis(1:n))1790C2 = C2 * hsheet**na1791C2 = C2 * Ngrad**(nb-2.0)17921793! Vclose velocity is split into Vfactor . (Phi0 - Phi)^(n-1)1794PhiG = SUM(NodalHydPot(1:n)*Basis(1:n))1795Phi0 = SUM(NodalPhi0(1:n)*Basis(1:n))1796Vfactor = SUM(NodalAr(1:n)*Basis(1:n)) * hsheet1797ng = SUM(NodalNg(1:n)*Basis(1:n))1798Vfactor = Vfactor * ABS(Phi0-PhiG)**(ng-1.0_dp)17991800Force = SUM( LoadVector(1:n)*Basis(1:n) )1801! contribution from volume source (using handle)1802LoadAtIP = ListGetElementReal( Load_h, Basis, Element, Found, GaussPoint=t)1803IF (Found) THEN1804Force = Force + LoadAtIP1805END IF1806Force = Force + Vfactor * (Phi0 + (ng - 1.0_dp) * PhiG)18071808!------------------------------------------------------------------------------1809! Loop over basis functions of both unknowns and weights1810!------------------------------------------------------------------------------1811DO p=1,NBasis1812DO q=1,NBasis1813!------------------------------------------------------------------------------1814! The diffusive-convective equation without stabilization1815!------------------------------------------------------------------------------1816M = CT * Basis(q) * Basis(p)1817!------------------------------------------------------------------------------1818! The diffusion term1819!------------------------------------------------------------------------------1820A = C2 * SUM( dBasisdx(q,1:dim) * dBasisdx(p,1:dim))1821A = A + Vfactor * ng * Basis(q) * Basis(p)18221823StiffMatrix(p,q) = StiffMatrix(p,q) + s * A1824MassMatrix(p,q) = MassMatrix(p,q) + s * M1825END DO1826ForceVector(p) = ForceVector(p) + s * Force * Basis(p)1827END DO1828END DO18291830!------------------------------------------------------------------------------1831END SUBROUTINE SheetCompose1832!------------------------------------------------------------------------------183318341835!------------------------------------------------------------------------------1836!> Return element local matrices and RSH vector for boundary conditions1837!> of diffusion convection equation:1838!------------------------------------------------------------------------------1839SUBROUTINE SheetBoundary( BoundaryMatrix, BoundaryVector, &1840LoadVector, Element, n, Nodes )1841!------------------------------------------------------------------------------1842USE MaterialModels1843USE Integration1844USE Differentials18451846IMPLICIT NONE1847REAL(KIND=dp) :: BoundaryMatrix(:,:), BoundaryVector(:), LoadVector(:)1848TYPE(Nodes_t) :: Nodes1849TYPE(Element_t) :: Element18501851INTEGER :: n18521853REAL(KIND=dp) :: Basis(n)1854REAL(KIND=dp) :: dBasisdx(n,3), detJ18551856REAL(KIND=dp) :: U, V, W, S1857REAL(KIND=dp) :: Force1858REAL(KIND=dp), POINTER :: U_Integ(:), V_Integ(:), W_Integ(:), S_Integ(:)18591860INTEGER :: i, t, q, p, N_Integ18611862TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff18631864LOGICAL :: stat1865!------------------------------------------------------------------------------1866BoundaryVector = 0.0_dp1867BoundaryMatrix = 0.0_dp1868!------------------------------------------------------------------------------1869! Integration stuff1870!------------------------------------------------------------------------------1871IntegStuff = GaussPoints( Element )1872U_Integ => IntegStuff % u1873V_Integ => IntegStuff % v1874W_Integ => IntegStuff % w1875S_Integ => IntegStuff % s1876N_Integ = IntegStuff % n18771878!------------------------------------------------------------------------------1879! Now we start integrating1880!------------------------------------------------------------------------------1881DO t=1,N_Integ1882U = U_Integ(t)1883V = V_Integ(t)1884W = W_Integ(t)1885!------------------------------------------------------------------------------1886! Basis function values & derivatives at the integration point1887!------------------------------------------------------------------------------1888stat = ElementInfo( Element, Nodes, u, v, w, detJ, &1889Basis, dBasisdx )18901891S = detJ * S_Integ(t)1892!------------------------------------------------------------------------------1893Force = SUM( LoadVector(1:n)*Basis(1:n) )18941895DO q=1,n1896BoundaryVector(q) = BoundaryVector(q) + s * Basis(q) * Force1897END DO1898END DO1899END SUBROUTINE SheetBoundary19001901!------------------------------------------------------------------------------1902SUBROUTINE ChannelCompose( MassMatrix, StiffMatrix, ForceVector, &1903NodalH, NodalHydPot, CArea, NodalKc, NodalAlphac, NodalBetac, NodalPhi0, NodalPhim, &1904NodalAc, Nodallc, Nodalng, NodalKs, NodalAlphas, NodalBetas, &1905NodalAfactor, NodalBfactor, EdgeTangent, Element, n, Nodes )1906!------------------------------------------------------------------------------1907USE MaterialModels1908USE Integration1909USE Differentials19101911IMPLICIT NONE19121913REAL(KIND=dp), DIMENSION(:) :: ForceVector1914REAL(KIND=dp), DIMENSION(:,:) :: MassMatrix, StiffMatrix1915REAL(KIND=dp) :: NodalHydPot(:), NodalBetas(:), NodalAlphas(:)1916REAL(KIND=dp) :: NodalH(:), CArea, EdgeTangent(3)1917REAL(KIND=dp) :: NodalPhi0(:), NodalPhim(:), NodalNg(:), NodalAfactor(:)1918REAL(KIND=dp) :: NodalKc(:), NodalAlphac(:), NodalAc(:), Nodallc(:)1919REAL(KIND=dp) :: NodalBetac(:), NodalKs(:), NodalBfactor(:)19201921INTEGER :: n19221923TYPE(Nodes_t) :: Nodes1924TYPE(Element_t), POINTER :: Element19251926!------------------------------------------------------------------------------1927! Local variables1928!------------------------------------------------------------------------------19291930REAL(KIND=dp) :: dBasisdx(n,3), detJ1931REAL(KIND=dp) :: Basis(n), dBasisds(n)19321933REAL(KIND=dp) :: Force19341935REAL(KIND=dp) :: A, M1936REAL(KIND=dp) :: Load19371938REAL(KIND=dp) :: x, y, z19391940INTEGER :: i, j, k, c, p, q, t, N_Integ, NBasis19411942TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff1943REAL(KIND=dp) :: s, u, v, w19441945REAL(KIND=dp) :: Vfactor, Phi0, PhiG, Afactor, Bfactor, GradPhim, dPw, Ffactor19461947REAL(KIND=dp) :: GradPhi, Ngrad, nbc, nbs, nas, nac, hsheet, ng, qc, Kc, Ks, lc1948REAL(KIND=dp), DIMENSION(:), POINTER :: U_Integ, V_Integ, W_Integ, S_Integ19491950LOGICAL :: Vms, Found, Transient, stat1951TYPE(ValueList_t), POINTER :: BodyForce19521953!------------------------------------------------------------------------------19541955Transient = GetString(GetSimulation(),'Simulation type',Found)=='transient'19561957ForceVector = 0.0_dp1958StiffMatrix = 0.0_dp1959MassMatrix = 0.0_dp1960Load = 0.0_dp19611962NBasis = n19631964!------------------------------------------------------------------------------1965! Integration stuff1966!------------------------------------------------------------------------------1967IntegStuff = GaussPoints( element )1968U_Integ => IntegStuff % u1969V_Integ => IntegStuff % v1970W_Integ => IntegStuff % w1971S_Integ => IntegStuff % s1972N_Integ = IntegStuff % n19731974!------------------------------------------------------------------------------1975! Now we start integrating1976!------------------------------------------------------------------------------1977DO t=1,N_Integ19781979u = U_Integ(t)1980v = V_Integ(t)1981w = W_Integ(t)19821983!------------------------------------------------------------------------------1984! Basis function values & derivatives at the integration point1985!------------------------------------------------------------------------------1986stat = ElementInfo( Element, Nodes, u, v, w, detJ, &1987Basis, dBasisdx )19881989s = detJ * S_Integ(t)19901991! derivative along the edge1992DO p=1, n1993dBasisds(p) = SUM(dBasisdx(p,1:3)*EdgeTangent(1:3))1994END DO1995!1996!------------------------------------------------------------------------------1997! Compute the effective conductivity K = k h^alpha |grad Phi|^{beta-2}1998! in the sheet ks and in the Channel Kc1999! gradPhi is here a scalar2000!------------------------------------------------------------------------------2001GradPhi = SUM(NodalHydPot(1:n)*dBasisds(1:n))2002Ngrad = ABS(GradPhi)2003IF (Ngrad < AEPS) Ngrad = AEPS2004nas = SUM(NodalAlphas(1:n)*Basis(1:n))2005nbs = SUM(NodalBetas(1:n)*Basis(1:n))2006nac = SUM(NodalAlphac(1:n)*Basis(1:n))2007nbc = SUM(NodalBetac(1:n)*Basis(1:n))2008hsheet = SUM(NodalH(1:n)*Basis(1:n))20092010Ks = SUM( NodalKs(1:n)*Basis(1:n))2011Ks = Ks * hsheet**nas2012Ks = Ks * Ngrad**(nbs-2.0)20132014Kc = SUM( NodalKc(1:n)*Basis(1:n))2015Kc = Kc * CArea**nac2016Kc = Kc * Ngrad**(nbc-2.0)20172018! Vclose velocity is split into Ac . S . (Phi0 - Phi)^n2019PhiG = SUM(NodalHydPot(1:n)*Basis(1:n))2020Phi0 = SUM(NodalPhi0(1:n)*Basis(1:n))20212022ng = SUM(NodalNg(1:n)*Basis(1:n))2023Vfactor = CArea*SUM(NodalAc(1:n)*Basis(1:n))2024Vfactor = Vfactor*ABS(Phi0-PhiG)**(ng-1.0_dp)20252026Afactor = SUM(NodalAfactor(1:n)*Basis(1:n))2027Bfactor = SUM(NodalBfactor(1:n)*Basis(1:n))2028Afactor = Afactor * Bfactor2029lc = SUM(Nodallc(1:n)*Basis(1:n))20302031qc = -Ks * GradPhi20322033GradPhim = SUM(NodalPhim(1:n)*dBasisds(1:n))2034dPw = GradPhi - GradPhim20352036IF ((CArea>0.0).OR.(qc*dPw>0.0)) THEN2037Ffactor = Afactor * lc * qc2038ELSE2039Ffactor = 0.0_dp2040END IF20412042Bfactor = Bfactor * SIGN((ABS(-Kc*GradPhi)+ABS(lc*qc)),GradPhi)20432044Force = Ffactor * GradPhim2045Force = Force + Vfactor*(Phi0+(ng-1.0_dp)*PhiG)20462047!------------------------------------------------------------------------------2048! Loop over basis functions of both unknowns and weights2049!------------------------------------------------------------------------------2050DO p=1,NBasis20512052Load = Force * Basis(p)2053ForceVector(p) = ForceVector(p) + s * Load20542055DO q=1,NBasis2056!------------------------------------------------------------------------------2057! The diffusive-convective equation without stabilization2058!------------------------------------------------------------------------------2059A = 0.0_dp2060!------------------------------------------------------------------------------2061! The diffusion term2062!------------------------------------------------------------------------------2063A = A + Kc * dBasisds(p) * dBasisds(q)2064A = A - Afactor * Kc * dPw * Basis(p) * dBasisds(q)2065A = A + Ffactor * Basis(p) * dBasisds(q)2066A = A + Vfactor*ng*Basis(p)*Basis(q)2067A = A + Bfactor * Basis(p) * dBasisds(q)20682069StiffMatrix(p,q) = StiffMatrix(p,q) + s * A2070END DO2071END DO20722073END DO2074!------------------------------------------------------------------------------2075!------------------------------------------------------------------------------2076END SUBROUTINE ChannelCompose2077!------------------------------------------------------------------------------20782079!------------------------------------------------------------------------------2080SUBROUTINE GetEvolveChannel(ALPHA, BETA, Qcc, CArea, NodalHydPot, NodalH, &2081NodalAlphac, NodalBetac, NodalKc, NodalPhi0, NodalPhim, NodalAc, Nodallc, Nodalng, &2082NodalKs, NodalAlphas, NodalBetas, NodalAfactor, NodalBfactor, &2083Tangent, Element, n, Nodes, LimitEffPres)2084!------------------------------------------------------------------------------2085USE MaterialModels2086USE Integration2087USE Differentials20882089IMPLICIT NONE209020912092REAL(KIND=dp) :: ALPHA, BETA, Qcc2093REAL(KIND=dp) :: NodalHydPot(:), NodalAlphas(:), NodalBetas(:)2094REAL(KIND=dp) :: NodalH(:), CArea, Tangent(3)2095REAL(KIND=dp) :: NodalPhi0(:), NodalPhim(:), NodalNg(:), NodalAfactor(:)2096REAL(KIND=dp) :: NodalKc(:), Nodalalphac(:), NodalAc(:), Nodallc(:)2097REAL(KIND=dp) :: NodalBetac(:), NodalKs(:), NodalBfactor(:)20982099INTEGER :: n21002101TYPE(Nodes_t) :: Nodes2102TYPE(Element_t), POINTER :: Element21032104LOGICAL :: LimitEffPres21052106!------------------------------------------------------------------------------2107! Local variables2108!------------------------------------------------------------------------------21092110REAL(KIND=dp) :: dBasisdx(n,3), detJ2111REAL(KIND=dp) :: Basis(n), dBasisds(n)211221132114REAL(KIND=dp) :: A, M2115REAL(KIND=dp) :: Load21162117REAL(KIND=dp) :: x, y, z21182119INTEGER :: i, j, k, c, p, q, t, N_Integ,NBasis21202121TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff2122REAL(KIND=dp) :: s, u, v, w21232124REAL(KIND=dp) :: Phi0, PhiG, EffPressatIP, Afactor, Bfactor, GradPhim, dPw, Ffactor21252126REAL(KIND=dp) :: GradPhi, Ngrad, nbc, hsheet, nas, nbs, nac, ng, qc, Kc, Ks, lc21272128REAL(KIND=dp), DIMENSION(:), POINTER :: U_Integ, V_Integ, W_Integ, S_Integ21292130REAL(KIND=dp) :: Xi, Pii, Nef, Vc21312132LOGICAL :: Vms, Found, Transient, stat21332134!------------------------------------------------------------------------------21352136Transient = GetString(GetSimulation(),'Simulation type',Found)=='transient'21372138NBasis = n21392140!------------------------------------------------------------------------------2141! Integration stuff2142!------------------------------------------------------------------------------2143IntegStuff = GaussPoints( element )2144U_Integ => IntegStuff % u2145V_Integ => IntegStuff % v2146W_Integ => IntegStuff % w2147S_Integ => IntegStuff % s2148N_Integ = IntegStuff % n21492150! Compute the force at the middle of the Edge element (u=0)2151u = 0.0_dp2152v = 0.0_dp2153w = 0.0_dp21542155!------------------------------------------------------------------------------2156! Basis function values & derivatives at the integration point2157!------------------------------------------------------------------------------2158stat = ElementInfo( Element, Nodes, u, v, w, detJ, &2159Basis, dBasisdx )216021612162! derivative along the edge2163DO p=1, n2164dBasisds(p) = SUM(dBasisdx(p,1:3)*EdgeTangent(1:3))2165END DO21662167!------------------------------------------------------------------------------2168! Compute the effective conductivity K = k h |grad Phi|^{beta-2}2169! in the sheet ks and in the Channel Kc2170! gradPhi is here a scalar2171!------------------------------------------------------------------------------2172!2173gradPhi = SUM(NodalHydPot(1:n)*dBasisds(1:n))2174Ngrad = ABS(GradPhi)2175IF (Ngrad < AEPS) Ngrad = AEPS21762177nas = SUM(NodalAlphas(1:n)*Basis(1:n))2178nbs = SUM(NodalBetas(1:n)*Basis(1:n))2179nbc = SUM(NodalBetac(1:n)*Basis(1:n))2180nac = SUM(NodalAlphac(1:n)*Basis(1:n))2181hsheet = SUM(NodalH(1:n)*Basis(1:n))21822183Ks = SUM( NodalKs(1:n) * Basis(1:n))2184Ks = Ks * hsheet**nas2185Ks = Ks * Ngrad**(nbs-2.0_dp)21862187Kc = SUM( NodalKc(1:n) * Basis(1:n))2188Kc = Kc * MAX(CArea,0.0)**(nac-1.0_dp)2189Kc = Kc * Ngrad**(nbc-2.0_dp)21902191PhiG = SUM(NodalHydPot(1:n)*Basis(1:n))2192Phi0 = SUM(NodalPhi0(1:n)*Basis(1:n))21932194ng = SUM(NodalNg(1:n)*Basis(1:n))2195Vc = SUM(NodalAc(1:n)*Basis(1:n))21962197IF (LimitEffPres) THEN2198EffPressatIP = MAX(Phi0-PhiG, 0.0_dp)2199ELSE2200EffPressatIP = Phi0-PhiG2201END IF22022203Vc = Vc*ABS(EffPressatIP)**(ng-1.0_dp)2204Vc = Vc*(EffPressatIP)22052206Afactor = SUM(NodalAfactor(1:n)*Basis(1:n))2207Bfactor = SUM(NodalBfactor(1:n)*Basis(1:n))2208lc = SUM(Nodallc(1:n)*Basis(1:n))22092210qc = -Ks * GradPhi22112212GradPhim = SUM(NodalPhim(1:n)*dBasisds(1:n))2213dPw = GradPhi - GradPhim22142215IF ((CArea>0.0).OR.(qc*dPw>0.0)) THEN2216Ffactor = lc * qc2217ELSE2218Ffactor = 0.0_dp2219END IF222022212222! Terms in ALPHA2223Xi = ABS(-Kc*GradPhi*GradPhi)2224Pii = -Afactor*(-Kc*GradPhi)*dPw2225ALPHA = Bfactor*(Xi - Pii) - Vc22262227! Terms in BETA2228Xi = ABS(lc*qc*GradPhi)2229Pii = -Afactor*(Ffactor)*dPw2230BETA = Bfactor*(Xi - Pii)22312232! Channel flux for output2233Qcc = ABS(MAX(CArea,0.0)*Kc*GradPhi)22342235!------------------------------------------------------------------------------2236END SUBROUTINE GetEvolveChannel2237!------------------------------------------------------------------------------22382239!------------------------------------------------------------------------------2240SUBROUTINE GetParametersSheet( Element, Material, N, SheetConductivity, alphas, &2241betas, Ev, ub, Snn, lr, hr, Ar, ng )2242!------------------------------------------------------------------------------2243USE MaterialModels2244USE Integration2245USE Differentials22462247IMPLICIT NONE2248REAL(KIND=dp) :: SheetConductivity(:), alphas(:), betas(:), Ev(:), ub(:), &2249Snn(:), lr(:), hr(:), Ar(:), ng(:)2250INTEGER :: N2251LOGICAL :: Found = .FALSE.2252TYPE(Element_t), POINTER :: Element2253TYPE(ValueList_t), POINTER :: Material22542255!------------------------------------------------------------------------------2256SheetConductivity = 0.0_dp2257SheetConductivity = ListGetReal(Material, 'Sheet Conductivity', n, Element % NodeIndexes, &2258Found, UnfoundFatal = .TRUE. )22592260alphas = 0.0_dp2261alphas(1:N) = ListGetReal( Material, 'Sheet Flow Exponent alpha', n, Element % NodeIndexes, &2262Found, UnfoundFatal = .TRUE. )22632264betas = 0.0_dp2265betas(1:N) = ListGetReal( Material, 'Sheet Flow Exponent beta', n, Element % NodeIndexes, &2266Found, UnfoundFatal = .TRUE. )22672268Ev = 0.0_dp2269Ev(1:N) = ListGetReal( Material, 'Englacial Void Ratio', n, Element % NodeIndexes, &2270Found, UnfoundFatal = .TRUE. )22712272ub = 0.0_dp2273ub(1:N) = ListGetReal( Material, 'Sliding Velocity', n, Element % NodeIndexes, &2274Found, UnfoundFatal = .TRUE. )22752276Snn = 0.0_dp2277Snn(1:N) = ListGetReal( Material, 'Ice Normal Stress', N, Element % NodeIndexes, &2278Found, UnfoundFatal = .TRUE. )2279IF (ANY(Snn(1:N) < 0.0)) THEN2280WRITE(Message,'(A)')'The Ice Normal Stress (overburden pressure) must be positive'2281CALL FATAL(SolverName, Message)2282END IF22832284lr = 0.0_dp2285lr(1:N) = ListGetReal( Material, 'Bedrock Bump Length', N, Element % NodeIndexes, &2286Found, UnfoundFatal = .TRUE. )22872288hr = 0.0_dp2289hr(1:N) = ListGetReal( Material, 'Bedrock Bump High', N, Element % NodeIndexes, &2290Found, UnfoundFatal = .TRUE. )22912292Ar = 0.0_dp2293Ar(1:N) = ListGetReal( Material, 'Sheet Closure Coefficient', N, Element % NodeIndexes, &2294Found, UnfoundFatal = .TRUE. )22952296ng = 0.0_dp2297ng(1:N) = ListGetReal( Material, 'Glen Exponent', N, Element % NodeIndexes, &2298Found, UnfoundFatal = .TRUE. )2299!------------------------------------------------------------------------------2300END SUBROUTINE GetParametersSheet2301!------------------------------------------------------------------------------23022303!------------------------------------------------------------------------------2304SUBROUTINE GetParametersChannel( Edge, Material, N, SheetConductivity, &2305ChannelConductivity, alphac, betac, alphas, betas, IceDensity, &2306Snn, Ac, ng, CCt, CCw, lc )2307!------------------------------------------------------------------------------2308USE MaterialModels2309USE Integration2310USE Differentials23112312IMPLICIT NONE2313REAL(KIND=dp) :: SheetConductivity(:), ChannelConductivity(:), alphac(:), &2314betac(:), alphas(:), betas(:), IceDensity(:), &2315Snn(:), Ac(:), ng(:), CCt(:), CCw(:), lc(:)2316INTEGER :: N2317LOGICAL :: Found = .FALSE.23182319TYPE(Element_t), POINTER :: Edge2320TYPE(ValueList_t), POINTER :: Material23212322!------------------------------------------------------------------------------23232324SheetConductivity = 0.0_dp2325SheetConductivity(1:n) = ListGetReal(Material, 'Sheet Conductivity', n, Edge % NodeIndexes, &2326Found, UnfoundFatal = .TRUE. )23272328ChannelConductivity = 0.0_dp2329ChannelConductivity(1:n) = ListGetReal(Material, 'Channel Conductivity', n, Edge % NodeIndexes, &2330Found, UnfoundFatal = .TRUE. )23312332alphac = 0.0_dp2333alphac(1:n) = ListGetReal( Material, 'Channel Flow Exponent Alpha', n, Edge % NodeIndexes, &2334Found, UnfoundFatal = .TRUE. )23352336betac = 0.0_dp2337betac(1:n) = ListGetReal( Material, 'Channel Flow Exponent Beta', n, Edge % NodeIndexes, &2338Found, UnfoundFatal = .TRUE. )23392340alphas = 0.0_dp2341alphas(1:n) = ListGetReal( Material, 'Sheet Flow Exponent alpha', n, Edge % NodeIndexes, &2342Found, UnfoundFatal = .TRUE. )23432344betas = 0.0_dp2345betas(1:n) = ListGetReal( Material, 'Sheet Flow Exponent beta', n, Edge % NodeIndexes, &2346Found, UnfoundFatal = .TRUE. )23472348IceDensity = 0.0_dp2349IceDensity(1:n) = ListGetReal( Material, 'Density', N, Edge % NodeIndexes, &2350Found, UnfoundFatal = .TRUE. )23512352Snn = 0.0_dp2353Snn(1:n) = ListGetReal( Material, 'Ice Normal Stress', N, Edge % NodeIndexes, &2354Found, UnfoundFatal = .TRUE. )2355IF (ANY(Snn(1:n) < 0.0)) THEN2356WRITE(Message,'(A)')'The Ice Normal Stress (overburden pressure) must be positive'2357CALL FATAL(SolverName, Message)2358END IF23592360Ac = 0.0_dp2361Ac(1:n) = ListGetReal( Material, 'Channel Closure Coefficient', N, Edge % NodeIndexes, &2362Found, UnfoundFatal = .TRUE. )23632364ng = 0.0_dp2365ng(1:n) = ListGetReal( Material, 'Glen Exponent', N, Edge % NodeIndexes, &2366Found, UnfoundFatal = .TRUE. )23672368CCt = 0.0_dp2369CCt(1:n) = ListGetReal( Material, 'Pressure Melting Coefficient', N, Edge % NodeIndexes, &2370Found, UnfoundFatal = .TRUE. )23712372CCw = 0.0_dp2373CCw(1:n) = ListGetReal( Material, 'Water Heat Capacity', N, Edge % NodeIndexes, &2374Found, UnfoundFatal = .TRUE. )23752376lc = 0.0_dp2377lc(1:n) = ListGetReal( Material, 'Sheet Width Over Channel', N, Edge % NodeIndexes, &2378Found, UnfoundFatal = .TRUE. )2379!------------------------------------------------------------------------------2380END SUBROUTINE GetParametersChannel2381!------------------------------------------------------------------------------23822383!------------------------------------------------------------------------------2384SUBROUTINE SheetDischargeCompute( &2385NodalHydPot, hsheet, ks, na, nb, &2386Discharge, Element, n, Nodes, NodeNumber )2387!------------------------------------------------------------------------------2388USE MaterialModels2389USE Integration2390USE Differentials23912392IMPLICIT NONE2393REAL(KIND=dp) :: NodalHydPot(:), Discharge(:)2394REAL(KIND=dp) :: na, nb, ks, hsheet2395INTEGER :: n, NodeNumber2396TYPE(Nodes_t) :: Nodes2397TYPE(Element_t), POINTER :: Element23982399!------------------------------------------------------------------------------2400! Local variables2401!------------------------------------------------------------------------------24022403REAL(KIND=dp) :: dBasisdx(n,3), detJ2404REAL(KIND=dp) :: Basis(n)24052406INTEGER :: i, dim24072408REAL(KIND=dp) :: u, v, w2409REAL(KIND=dp) :: gradPhi(3), Ngrad2410LOGICAL :: stat24112412!------------------------------------------------------------------------------2413! the dimension for this solver is given by the dim of the elements2414! should work for 2D and 3D general problem2415dim = Element % TYPE % DIMENSION24162417u = Element % Type % NodeU(NodeNumber)2418v = Element % Type % NodeV(NodeNumber)2419w = Element % Type % NodeW(NodeNumber)24202421stat = ElementInfo( Element, Nodes, u, v, w, detJ, Basis, dBasisdx )24222423!------------------------------------------------------------------------------2424! Compute Sheet Discharge = k h^alphas |grad Phi|^{beta-2} grad Phi2425!------------------------------------------------------------------------------2426gradPhi = 0.0_dp2427DO i=1, dim2428gradPhi(i) = SUM(NodalHydPot(1:n)*dBasisdx(1:n,i))2429END DO2430Ngrad = SQRT(SUM(gradPhi(1:dim)*gradPhi(1:dim)))2431IF (Ngrad < AEPS) Ngrad = AEPS24322433DO i=1,dim2434Discharge(i) = -ks * (hsheet**na) * (Ngrad**(nb-2.0_dp)) * gradPhi(i)2435END DO24362437!------------------------------------------------------------------------------2438END SUBROUTINE SheetDischargeCompute2439!------------------------------------------------------------------------------2440!------------------------------------------------------------------------------2441END SUBROUTINE GlaDSCoupledsolver2442!------------------------------------------------------------------------------24432444! *****************************************************************************/2445!> Dummy solver just to declare the Sheet Thickness variable2446RECURSIVE SUBROUTINE GlaDSsheetThickDummy( Model,Solver,Timestep,TransientSimulation )2447!******************************************************************************2448!******************************************************************************2449USE Differentials2450USE MaterialModels2451USE DefUtils24522453!------------------------------------------------------------------------------2454IMPLICIT NONE2455!------------------------------------------------------------------------------2456!------------------------------------------------------------------------------2457! External variables2458!------------------------------------------------------------------------------2459TYPE(Model_t) :: Model2460TYPE(Solver_t), TARGET :: Solver2461LOGICAL :: TransientSimulation2462REAL(KIND=dp) :: Timestep24632464RETURN2465END SUBROUTINE GlaDSsheetThickDummy2466!------------------------------------------------------------------------------246724682469