Merge pull request #50 from gemshub/CCM-model-devel

added option for CCM model where capacitance is sqrt(IS)/alpha

Author: gdmiron

GEMS3K/ipm_chemical.cpp

@@ -523,7 +523,7 @@ case DC_SCM_SPECIES:
                   }
                }
             }
-            else if( pm.SCM[k][ist] == SC_BSM || pm.SCM[k][ist] == SC_CCM )
+            else if( pm.SCM[k][ist] == SC_BSM || pm.SCM[k][ist] == SC_CCM || pm.SCM[k][ist] == SC_ISCCM )
             { // Basic Stern model, Christl & Kretzschmar, 1999
             // New CD version of TLM  added 25.10.2004
                if( fabs( CD0 ) > 1e-20 && fabs( CDb ) > 1e-20 )

GEMS3K/ipm_chemical2.cpp

@@ -595,9 +595,10 @@ void TMultiBase::IS_EtaCalc()
                         pm.XetaB[k][ist] += pm.X[j]*Ez;
                         pm.XetaA[k][ist] += pm.X[j]*CD0;  // added for testing
                     }
-                    else if( pm.SCM[k][ist] == SC_CCM )
+                    else if( pm.SCM[k][ist] == SC_CCM || pm.SCM[k][ist] == SC_ISCCM )
                     { // Added 25.07.03 to implement the extended CCM Nilsson ea 1996
 // New CD version of BSM  added 25.10.2004
+                        // 2026 SC_ISCCM uses dependcy of capacity with IS C = sqrt(IS)/alpha, alpha is given instead of capacity
                            pm.XetaB[k][ist] += pm.X[j] * CDb;
                            pm.XetaA[k][ist] += pm.X[j] * CD0;
                     }
@@ -697,6 +698,7 @@ TMultiBase::GouyChapman(  long int, long int, long int k )
         switch( pm.SCM[k][ist] )
         {
         case SC_CCM:  // Constant-Capacitance Model Schindler, extension Nilsson
+        case SC_ISCCM: // 2026 SC_ISCCM uses dependcy of capacity with IS C = sqrt(IS)/alpha, alpha is given instead of capacity
             SigA = pm.Xetaf[k][ist] + XetaA[ist];
             SigDDL = -SigA - XetaB[ist];
             SigB = XetaB[ist];
@@ -796,9 +798,13 @@ TMultiBase::GouyChapman(  long int, long int, long int k )
             pm.XpsiB[k][ist] = PsiD;
             break;
         case SC_CCM:  // Constant-Capacitance Model Schindler, ext. Nilsson
+        case SC_ISCCM: // 2026 SC_ISCCM uses dependcy of capacity with IS C = sqrt(IS)/alpha, alpha is given instead of capacity
             if( pm.XcapB[k][ist] > 0.001 )
             {  // Classic CCM Schindler with inner-sphere species only
-               PsiA = SigA / pm.XcapA[k][ist];
+                if (pm.SCM[k][ist] == SC_ISCCM)
+                    PsiA = SigA / (sqrt(I)/pm.XcapA[k][ist]);
+                else
+                    PsiA = SigA / pm.XcapA[k][ist];
                if( fabs( PsiA ) > 0.7 ) // truncated 0-plane potential
                {
                    PsiA = PsiA<0? -0.7: 0.7;

GEMS3K/ipm_simplex.cpp

@@ -56,7 +56,7 @@ void TMultiBase::AutoInitialApproximation( )
         for( i=0; i<pm.N; i++)
              DU[i+Q] = 0.;
         EPS = base_param()->EPS; //  13.10.00  KC  DK
-        GZ = 1./EPS;    
+        GZ = 1./EPS;
 
         T=0; // Calcuation of all non-zero values in A and G arrays
         for(i=0;i<pm.L;i++)

GEMS3K/m_const_base.h

@@ -253,7 +253,7 @@ enum sorption_control {
     // EDL interface models - separate for site types in v. 3.1
     SC_DDLM = 'D',  SC_CCM = 'C',     SC_TLM = 'T',   SC_MTL = 'M',
     SC_MXC = 'E',   SC_NNE = 'X',     SC_IEV  = 'I',  SC_BSM = 'S',
-    SC_3LM = '3', SC_NOT_USED = 'N',
+    SC_3LM = '3', SC_NOT_USED = 'N',  SC_ISCCM = 'P',
 
     // Methods of Surface Activity Terms calculation
     SAT_COMP = 'C', SAT_NCOMP = 'N', SAT_SOLV = 'S', SAT_INDEF = 'I',