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

@@ -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

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))

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,17 +128,17 @@ 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
     SELECT CASE (iFriction)
     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,33 +184,33 @@ 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.)
       END IF
 
     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
    END IF
 
    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