From ee057a9c6c85e1200ea0a2a964769c431cb285c6 Mon Sep 17 00:00:00 2001 From: RupertGladstone Date: Fri, 4 Oct 2024 17:46:45 +0300 Subject: [PATCH] Allow SSA code to use a different grounded mask by setting "Grounded Mask Variable Name" in constants section. (implemented in SSAmaterial properties, SSA solver and thickness solver; default is GroundedMask for backward compatibility) --- elmerice/Solvers/SSASolver.F90 | 37 ++++++++--- elmerice/Solvers/ThicknessSolver.F90 | 10 ++- elmerice/Utils/SSAMaterialModels.F90 | 96 ++++++++++++++++------------ 3 files changed, 89 insertions(+), 54 deletions(-) diff --git a/elmerice/Solvers/SSASolver.F90 b/elmerice/Solvers/SSASolver.F90 index 862bba0053..a9396477c3 100644 --- a/elmerice/Solvers/SSASolver.F90 +++ b/elmerice/Solvers/SSASolver.F90 @@ -70,8 +70,8 @@ SUBROUTINE SSABasalSolver( Model,Solver,dt,TransientSimulation ) TYPE(Nodes_t) :: ElementNodes TYPE(Element_t),POINTER :: CurrentElement, Element, ParentElement, BoundaryElement TYPE(Matrix_t),POINTER :: StiffMatrix - TYPE(ValueList_t), POINTER :: SolverParams, BodyForce, Material, BC - TYPE(Variable_t), POINTER :: PointerToVariable, ZsSol, ZbSol, VeloSol,strbasemag,Ceff + TYPE(ValueList_t), POINTER :: SolverParams, BodyForce, Material, BC, Constants + TYPE(Variable_t), POINTER :: PointerToVariable, ZsSol, ZbSol, VeloSol, strbasemag, Ceff, GMSol LOGICAL :: AllocationsDone = .FALSE., Found, GotIt, CalvingFront, UnFoundFatal=.TRUE. LOGICAL :: stat @@ -91,15 +91,16 @@ SUBROUTINE SSABasalSolver( Model,Solver,dt,TransientSimulation ) REAL(KIND=dp), ALLOCATABLE :: STIFF(:,:), LOAD(:), FORCE(:), & NodalGravity(:), NodalViscosity(:), NodalDensity(:), & - NodalZs(:), NodalZb(:), NodalU(:), NodalV(:),Basis(:) + NodalZs(:), NodalZb(:), NodalU(:), NodalV(:),Basis(:), NodalGM(:) REAL(KIND=dp) :: DetJ,UnLimit,un,un_max,FillValue - CHARACTER(LEN=MAX_NAME_LEN) :: SolverName, ZsName, ZbName + CHARACTER(LEN=MAX_NAME_LEN) :: SolverName, ZsName, ZbName, MaskName #ifdef USE_ISO_C_BINDINGS REAL(KIND=dp) :: at, at0 #else REAL(KIND=dp) :: at, at0, CPUTime, RealTime #endif + LOGICAL :: PartlyGroundedElement LOGICAL :: SEP ! Sub-element parametrization for Grounding line INTEGER :: GLnIP ! number of Integ. Points for GL Sub-element parametrization @@ -108,7 +109,7 @@ SUBROUTINE SSABasalSolver( Model,Solver,dt,TransientSimulation ) SAVE NodalGravity, NodalViscosity, NodalDensity, & NodalZs, NodalZb, & NodalU, NodalV, & - Basis + Basis, nodalGM !------------------------------------------------------------------------------ PointerToVariable => Solver % Variable @@ -185,13 +186,13 @@ SUBROUTINE SSABasalSolver( Model,Solver,dt,TransientSimulation ) M = Model % Mesh % NumberOfNodes IF (AllocationsDone) DEALLOCATE(FORCE, LOAD, STIFF, NodalGravity, & NodalViscosity, NodalDensity, & - NodalZb, NodalZs, NodalU, NodalV, & + NodalZb, NodalZs, NodalU, NodalV, NodalGM, & ElementNodes % x, & ElementNodes % y, ElementNodes % z ,Basis) ALLOCATE( FORCE(STDOFs*N), LOAD(N), STIFF(STDOFs*N,STDOFs*N), & NodalGravity(N), NodalDensity(N), NodalViscosity(N), & - NodalZb(N), NodalZs(N), NodalU(N), NodalV(N), & + NodalZb(N), NodalZs(N), NodalU(N), NodalV(N), NodalGM(N), & ElementNodes % x(N), ElementNodes % y(N), ElementNodes % z(N), & Basis(N), & STAT=istat ) @@ -534,10 +535,20 @@ SUBROUTINE SSABasalSolver( Model,Solver,dt,TransientSimulation ) End do un=sqrt(un) + PartlyGroundedElement = .FALSE. + IF (SEP) THEN + Constants => GetConstants() + MaskName = ListGetString(Constants,'Grounded Mask Variable Name',UnFoundFatal=.FALSE.,DefValue='GroundedMask') + GMSol => VariableGet( CurrentModel % Variables, MaskName, UnFoundFatal=.TRUE. ) + CALL GetLocalSolution( NodalGM,UElement=Element,UVariable=GMSol) + PartlyGroundedElement=(ANY(NodalGM(1:n).GE.0._dp).AND.ANY(NodalGM(1:n).LT.0._dp)) + END IF + h=MAX(SUM(Basis(1:n)*(NodalZs(1:n)-NodalZb(1:n))),MinH) Ceff%Values(Ceff%Perm(NodeIndexes(i)))= & - SSAEffectiveFriction(Element,n,Basis,un,SEP,.TRUE.,h,SUM(NodalDensity(1:n)*Basis(1:n)),rhow,sealevel) - End do + SSAEffectiveFriction(Element,n,Basis,un,SEP,PartlyGroundedElement,h,SUM(NodalDensity(1:n)*Basis(1:n)),rhow,sealevel) +! SSAEffectiveFriction(Element,n,Basis,un,SEP,.TRUE.,h,SUM(NodalDensity(1:n)*Basis(1:n)),rhow,sealevel) + END DO End IF END DO @@ -596,6 +607,8 @@ SUBROUTINE LocalMatrixUVSSA( STIFF, FORCE, Element, n, Nodes, gravity, & LOGICAL :: Stat, NewtonLin INTEGER :: i, j, t, p, q , dim TYPE(GaussIntegrationPoints_t) :: IP + CHARACTER(LEN=MAX_NAME_LEN) :: MaskName + TYPE(ValueList_t), POINTER :: Constants TYPE(Nodes_t) :: Nodes !------------------------------------------------------------------------------ @@ -608,8 +621,12 @@ SUBROUTINE LocalMatrixUVSSA( STIFF, FORCE, Element, n, Nodes, gravity, & ! Use Newton Linearisation NewtonLin = (Newton.AND.(cm.NE.1.0_dp)) + PartlyGroundedElement = .FALSE. IF (SEP) THEN - GMSol => VariableGet( CurrentModel % Variables, 'GroundedMask',UnFoundFatal=.TRUE. ) + Constants => GetConstants() + MaskName = ListGetString(Constants,'Grounded Mask Variable Name',UnFoundFatal=.FALSE.,DefValue='GroundedMask') + GMSol => VariableGet( CurrentModel % Variables, MaskName, UnFoundFatal=.TRUE. ) +! GMSol => VariableGet( CurrentModel % Variables, 'GroundedMask',UnFoundFatal=.TRUE. ) CALL GetLocalSolution( NodalGM,UElement=Element,UVariable=GMSol) PartlyGroundedElement=(ANY(NodalGM(1:n).GE.0._dp).AND.ANY(NodalGM(1:n).LT.0._dp)) IF (PartlyGroundedElement) THEN diff --git a/elmerice/Solvers/ThicknessSolver.F90 b/elmerice/Solvers/ThicknessSolver.F90 index b4ae53684e..c3ee78ffdf 100644 --- a/elmerice/Solvers/ThicknessSolver.F90 +++ b/elmerice/Solvers/ThicknessSolver.F90 @@ -706,10 +706,15 @@ SUBROUTINE LocalMatrix( STIFF, MASS, FORCE,& INTEGER :: i,j,t,p,q, n,FIPcount REAL(KIND=dp) :: smbE, bmbE, area, MinH REAL(KIND=dp) :: smbAtIP, bmbAtIP, GMatIP, rho, rhow, hh, sealevel,FFI + TYPE(ValueList_t), POINTER :: Constants + CHARACTER(LEN=MAX_NAME_LEN) :: MaskName !------------------------------------------------------------------------------ IF (SEM) THEN - GMSol => VariableGet( CurrentModel % Variables, 'GroundedMask',UnFoundFatal=.TRUE. ) + Constants => GetConstants() + MaskName = ListGetString(Constants,'Grounded Mask Variable Name',UnFoundFatal=.FALSE.,DefValue='GroundedMask') + GMSol => VariableGet( CurrentModel % Variables,MaskName,UnFoundFatal=.TRUE. ) +! GMSol => VariableGet( CurrentModel % Variables, 'GroundedMask',UnFoundFatal=.TRUE. ) CALL GetLocalSolution( NodalGM,UElement=Element,UVariable=GMSol) PartlyGroundedElement=(ANY(NodalGM(1:nCoord).GE.0._dp).AND.ANY(NodalGM(1:nCoord).LT.0._dp)) IF (PartlyGroundedElement) THEN @@ -791,7 +796,8 @@ SUBROUTINE LocalMatrix( STIFF, MASS, FORCE,& ! Numerical integration: ! ---------------------- IF (SEM) THEN - GMSol => VariableGet( CurrentModel % Variables, 'GroundedMask',UnFoundFatal=.TRUE. ) + GMSol => VariableGet( CurrentModel % Variables,MaskName,UnFoundFatal=.TRUE. ) +! GMSol => VariableGet( CurrentModel % Variables, 'GroundedMask',UnFoundFatal=.TRUE. ) CALL GetLocalSolution( NodalGM,UElement=Element,UVariable=GMSol) CALL GetLocalSolution( NodalThick,UElement=Element,UVariable=Solver % Variable) PartlyGroundedElement=(ANY(NodalGM(1:nCoord).GE.0._dp).AND.ANY(NodalGM(1:nCoord).LT.0._dp)) diff --git a/elmerice/Utils/SSAMaterialModels.F90 b/elmerice/Utils/SSAMaterialModels.F90 index 382d3c2415..b2d42f2871 100644 --- a/elmerice/Utils/SSAMaterialModels.F90 +++ b/elmerice/Utils/SSAMaterialModels.F90 @@ -45,12 +45,12 @@ MODULE SSAMaterialModels !-------------------------------------------------------------------------------- !> Return the effective friction coefficient !-------------------------------------------------------------------------------- - FUNCTION SSAEffectiveFriction(Element,n,Basis,ub,SEP,PartlyGrounded,h,rho,rhow,sealevel,SlipDer) RESULT(Slip) + FUNCTION SSAEffectiveFriction(Element,nn,Basis,ub,SEP,PartlyGrounded,h,rho,rhow,sealevel,SlipDer) RESULT(Slip) IMPLICIT NONE REAL(KIND=dp) :: Slip ! the effective friction coefficient TYPE(Element_t), POINTER :: Element ! the current element - INTEGER :: n ! number of nodes + INTEGER :: nn ! number of nodes REAL(KIND=dp) :: Basis(:) ! basis functions REAL(KIND=dp) :: ub ! the velocity for non-linear friction laws LOGICAL :: SEP ! Sub-Element Parametrisation of the friction @@ -75,25 +75,26 @@ FUNCTION SSAEffectiveFriction(Element,n,Basis,ub,SEP,PartlyGrounded,h,rho,rhow,s REAL(KIND=dp) :: Slip2, gravity, qq, hafq REAL(KIND=dp) :: fm,fq,MinN,MaxN,U0 REAL(KIND=dp) :: alpha,beta,fB - INTEGER :: GLnIP + INTEGER :: GLnIP,ii - REAL(KIND=dp),DIMENSION(n) :: NodalBeta, NodalGM, NodalBed, NodalLinVelo,NodalC,NodalN + REAL(KIND=dp),DIMENSION(nn) :: NodalBeta, NodalGM, NodalBed, NodalLinVelo,NodalC,NodalN REAL(KIND=dp) :: bedrock,Hf,fC,fN,LinVelo LOGICAL :: Found, NeedN + SAVE FirstTime Material => GetMaterial(Element) ! Allow user-named grounded mask - MaskName = ListGetString(Material, 'SSA Friction mask name',Found, UnFoundFatal=.FALSE.) - IF (.NOT.Found) THEN - MaskName = 'GroundedMask' + Constants => GetConstants() + MaskName = ListGetString(Constants,'Grounded Mask Variable Name',UnFoundFatal=.FALSE.,DefValue='GroundedMask') + IF (FirstTime) THEN + WRITE( Message, * ) 'Grounded mask name for SSA friction is:', MaskName + CALL INFO("SSAEffectiveFriction", Message, level=5) END IF - WRITE( Message, * ) 'Grounded mask name for SSA friction is:', MaskName - CALL INFO("SSAEffectiveFriction", Message, level=5) - + ! Sub - element GL parameterisation IF (SEP) THEN GMSol => VariableGet( CurrentModel % Variables, MaskName,UnFoundFatal=.TRUE. ) @@ -103,7 +104,7 @@ FUNCTION SSAEffectiveFriction(Element,n,Basis,ub,SEP,PartlyGrounded,h,rho,rhow,s CALL GetLocalSolution( NodalBed,UElement=Element,UVariable= BedrockSol) END IF -! Friction law + ! Friction law NodeIndexes => Element % NodeIndexes Friction = ListGetString(Material, 'SSA Friction Law',Found, UnFoundFatal=.TRUE.) @@ -127,8 +128,8 @@ FUNCTION SSAEffectiveFriction(Element,n,Basis,ub,SEP,PartlyGrounded,h,rho,rhow,s ! coefficient for all friction parameterisations NodalBeta = 0.0_dp - NodalBeta(1:n) = ListGetReal( & - Material, 'SSA Friction Parameter', n, NodeIndexes(1:n), Found,& + NodalBeta(1:nn) = ListGetReal( & + Material, 'SSA Friction Parameter', nn, NodeIndexes(1:nn), Found,& UnFoundFatal=.TRUE.) ! for nonlinear powers of sliding velocity @@ -136,8 +137,8 @@ FUNCTION SSAEffectiveFriction(Element,n,Basis,ub,SEP,PartlyGrounded,h,rho,rhow,s CASE(REG_COULOMB_JOU,REG_COULOMB_GAG,WEERTMAN,BUDD) fm = ListGetConstReal( Material, 'SSA Friction Exponent', Found , UnFoundFatal=.TRUE.) NodalLinVelo = 0.0_dp - NodalLinVelo(1:n) = ListGetReal( & - Material, 'SSA Friction Linear Velocity', n, NodeIndexes(1:n), Found,& + NodalLinVelo(1:nn) = ListGetReal( & + Material, 'SSA Friction Linear Velocity', nn, NodeIndexes(1:nn), Found,& UnFoundFatal=.TRUE.) CASE DEFAULT END SELECT @@ -153,7 +154,6 @@ FUNCTION SSAEffectiveFriction(Element,n,Basis,ub,SEP,PartlyGrounded,h,rho,rhow,s IF (.NOT. Found) THEN IF (FirstTime) THEN CALL INFO("SSAEffectiveFriction","> SSA Friction need N < not found, assuming false",level=3) - FirstTime = .FALSE. END IF NeedN = .FALSE. END IF @@ -165,7 +165,7 @@ FUNCTION SSAEffectiveFriction(Element,n,Basis,ub,SEP,PartlyGrounded,h,rho,rhow,s NSol => VariableGet( CurrentModel % Variables, 'Effective Pressure', UnFoundFatal=.TRUE. ) CALL GetLocalSolution( NodalN,UElement=Element, UVariable=NSol) MinN = ListGetConstReal( Material, 'SSA Min Effective Pressure', Found, UnFoundFatal=.TRUE.) - fN = SUM( NodalN(1:n) * Basis(1:n) ) + fN = SUM( NodalN(1:nn) * Basis(1:nn) ) fN = MAX(fN, MinN) ! Effective pressure should be >0 (for the friction law) MaxN = ListGetConstReal( Material, 'SSA Max Effective Pressure', Found, UnFoundFatal=.FALSE.) IF (Found) fN = MIN(fN, MaxN) @@ -175,7 +175,6 @@ FUNCTION SSAEffectiveFriction(Element,n,Basis,ub,SEP,PartlyGrounded,h,rho,rhow,s SELECT CASE (iFriction) CASE(BUDD) - Constants => GetConstants() gravity = ListGetConstReal( Constants, 'Gravity Norm', UnFoundFatal=.TRUE. ) ! calculate haf from N = rho_i g z* qq = ListGetConstReal( Material, 'SSA Haf Exponent', Found, UnFoundFatal=.TRUE.) @@ -185,10 +184,10 @@ FUNCTION SSAEffectiveFriction(Element,n,Basis,ub,SEP,PartlyGrounded,h,rho,rhow,s IF (iFriction .NE. REG_COULOMB_JOU) THEN fq = ListGetConstReal( Material, 'SSA Friction Post-Peak', Found, UnFoundFatal=.TRUE. ) NodalC = 0.0_dp - NodalC(1:n) = ListGetReal( & - Material, 'SSA Friction Maximum Value', n, NodeIndexes(1:n), Found,& + NodalC(1:nn) = ListGetReal( & + Material, 'SSA Friction Maximum Value', nn, NodeIndexes(1:nn), Found,& UnFoundFatal=.TRUE.) - fC = SUM( NodalC(1:n) * Basis(1:n) ) + fC = SUM( NodalC(1:nn) * Basis(1:nn) ) END IF IF (iFriction .NE. REG_COULOMB_GAG) THEN U0 = ListGetConstReal( Material, 'SSA Friction Threshold Velocity', Found, UnFoundFatal=.TRUE.) @@ -196,14 +195,14 @@ FUNCTION SSAEffectiveFriction(Element,n,Basis,ub,SEP,PartlyGrounded,h,rho,rhow,s END SELECT - Beta=SUM(Basis(1:n)*NodalBeta(1:n)) + Beta=SUM(Basis(1:nn)*NodalBeta(1:nn)) IF (SEP) THEN ! Floating - IF (ALL(NodalGM(1:n).LT.0._dp)) THEN + IF (ALL(NodalGM(1:nn).LT.0._dp)) THEN beta=0._dp ELSE IF (PartlyGrounded) THEN - bedrock = SUM( NodalBed(1:n) * Basis(1:n) ) + bedrock = SUM( NodalBed(1:nn) * Basis(1:nn) ) Hf= rhow * (sealevel-bedrock) / rho if (h.lt.Hf) beta=0._dp END IF @@ -211,7 +210,7 @@ FUNCTION SSAEffectiveFriction(Element,n,Basis,ub,SEP,PartlyGrounded,h,rho,rhow,s Slip2=0.0_dp IF (iFriction .NE. LINEAR) THEN - LinVelo = SUM( NodalLinVelo(1:n) * Basis(1:n) ) + LinVelo = SUM( NodalLinVelo(1:nn) * Basis(1:nn) ) IF ((iFriction == WEERTMAN).AND.(fm==1.0_dp)) iFriction=LINEAR Slip2=1.0_dp IF (ub < LinVelo) then @@ -246,10 +245,10 @@ FUNCTION SSAEffectiveFriction(Element,n,Basis,ub,SEP,PartlyGrounded,h,rho,rhow,s IF (PRESENT(SlipDer)) SlipDer = Slip2 * Slip * ((fm-1.0_dp) / (ub*ub) - & fm*fq*fB*ub**(fq-2.0_dp)/(1.0_dp+fB*ub**fq)) + CASE(REG_COULOMB_HYB) ! The sandard "SSA friction parameter" is taken as the effective pressure threshold. ! Max val is same as REG_COULMB_GAG ! Threshold vel is same as REG_COULOMB_JOU - CASE(REG_COULOMB_HYB) IF (fq.NE.1.0_dp) THEN CALL Fatal('SSAEffectiveFriction','Expecting unity post peak exponent') END IF @@ -262,9 +261,11 @@ FUNCTION SSAEffectiveFriction(Element,n,Basis,ub,SEP,PartlyGrounded,h,rho,rhow,s IF (NeedN) Slip = Slip * fN IF (PRESENT(SlipDer)) SlipDer = Slip2 * Slip * ((fm-1.0_dp) / (ub*ub) - & fm*ub**(-1.0_dp)/(ub+U0)) - - END SELECT + END SELECT + + FirstTime = .FALSE. + END FUNCTION SSAEffectiveFriction !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -283,7 +284,7 @@ FUNCTION ComputeMeanFriction(Element,n,ElementNodes,STDOFs,NodalU,NodalV,NodalZs INTEGER :: GLnIP REAL(KIND=dp) :: sealevel,rhow - TYPE(ValueList_t), POINTER :: Material + TYPE(ValueList_t), POINTER :: Material, Constants CHARACTER(LEN=MAX_NAME_LEN) :: MaskName LOGICAL :: PartlyGroundedElement TYPE(Variable_t),POINTER :: GMSol @@ -295,16 +296,20 @@ FUNCTION ComputeMeanFriction(Element,n,ElementNodes,STDOFs,NodalU,NodalV,NodalZs REAL(KIND=dp) :: Ceff LOGICAL :: stat, Found INTEGER :: t + LOGICAL :: FirstTime = .TRUE. + + SAVE FirstTime ! Allow user-named grounded mask Material => GetMaterial(Element) - MaskName = ListGetString(Material, 'SSA Friction mask name',Found, UnFoundFatal=.FALSE.) - IF (.NOT.Found) THEN - MaskName = 'GroundedMask' - END IF - WRITE( Message, * ) 'Grounded mask name for SSA friction is:', MaskName - CALL INFO("ComputeMeanFriction", Message, level=5) + Constants => GetConstants() + MaskName = ListGetString(Constants,'Grounded Mask Variable Name',UnFoundFatal=.FALSE.,DefValue='GroundedMask') + IF (FirstTime) THEN + WRITE( Message, * ) 'Grounded mask name for SSA friction is:', MaskName + CALL INFO("ComputeMeanFriction", Message, level=5) + END IF + strbasemag=0._dp IF (SEP) THEN GMSol => VariableGet( CurrentModel % Variables, MaskName,UnFoundFatal=.TRUE. ) @@ -343,6 +348,8 @@ FUNCTION ComputeMeanFriction(Element,n,ElementNodes,STDOFs,NodalU,NodalV,NodalZs strbasemag=tb/area + FirstTime = .FALSE. + END FUNCTION ComputeMeanFriction !-------------------------------------------------------------------------------- @@ -367,25 +374,28 @@ FUNCTION SSAEffectiveBMB(Element,nn,Basis,SEM,BMB,hh,FIPcount,rho,rhow,sealevel, REAL(KIND=dp),INTENT(IN),OPTIONAL :: rho,rhow,sealevel ! to calculate floatation for SEM3 REAL(KIND=dp),INTENT(IN),OPTIONAL :: FAF ! Floating area fraction for SEM1 - TYPE(ValueList_t), POINTER :: Material + TYPE(ValueList_t), POINTER :: Material, Constants TYPE(Variable_t), POINTER :: GMSol,BedrockSol CHARACTER(LEN=MAX_NAME_LEN) :: MeltParam, MaskName REAL(KIND=dp),DIMENSION(nn) :: NodalBeta, NodalGM, NodalBed, NodalLinVelo,NodalC REAL(KIND=dp) :: bedrock,Hf + LOGICAL :: FirstTime = .TRUE. LOGICAL :: Found + SAVE FirstTime + Material => GetMaterial(Element) ! Allow user-named grounded mask - MaskName = ListGetString(Material, 'SSA BMB mask name',Found, UnFoundFatal=.FALSE.) - IF (.NOT.Found) THEN - MaskName = 'GroundedMask' + Constants => GetConstants() + MaskName = ListGetString(Constants,'Grounded Mask Variable Name',UnFoundFatal=.FALSE.,DefValue='GroundedMask') + IF (FirstTime) THEN + WRITE( Message, * ) 'Grounded mask name for SSA BMB is:', MaskName + CALL INFO("SSAEffectiveBMB", Message, level=5) END IF - WRITE( Message, * ) 'Grounded mask name for SSA BMB is:', MaskName - CALL INFO("SSAEffectiveBMB", Message, level=5) - + ! Sub - element GL parameterisation IF (SEM) THEN GMSol => VariableGet( CurrentModel % Variables, MaskName,UnFoundFatal=.TRUE. ) @@ -432,6 +442,8 @@ FUNCTION SSAEffectiveBMB(Element,nn,Basis,SEM,BMB,hh,FIPcount,rho,rhow,sealevel, CALL FATAL("SSAEffectiveBMB",Message) END SELECT + + FirstTime = .FALSE. END FUNCTION SSAEffectiveBMB