TITLE Example 12a.--Advective and diffusive transport of heat and solutes.
      Constant boundary condition at one end, closed at other.
      The problem is designed so that temperature should equal Na-conc 
      (in mmol/kgw) after diffusion. Compares with analytical solution
      for 20-cell and 60-cell models.
EXCHANGE_SPECIES
  Na+ + X- = NaX
    log_k       0.0
    -gamma      4.0     0.075
  H+ + X- = HX
    log_k       -99.
    -gamma      9.0     0.0
  K+ + X- = KX
    log_k       0.0
    -gamma       3.5     0.015
#
#   20-cell model, initial conditions
#
SOLUTION 0   Fixed temp 24C, and NaCl conc (first type boundary cond) at inlet
    units  mol/kgw
    temp 24
    pH  7.0
    pe  12.0   O2(g) -0.67
    Na  24.e-3
    Cl  24.e-3
SOLUTION 1-19   24.0 mM KNO3
    units mol/kgw
    temp  0                 # Incoming solution 0C
    pH    7.0
    pe   12.0   O2(g) -0.67
    K    24.e-3
    N(5) 24.e-3
EXCHANGE 1-19
    KX    0.048
SOLUTION 20  Same as soln 0 in cell 20 at closed column end (second type boundary cond)
    units  mol/kgw
    temp 24
    pH  7.0
    pe  12.0   O2(g) -0.67
    Na  24.e-3
    Cl  24.e-3
EXCHANGE 20
    NaX    0.048
PRINT
   -reset      false
   -echo_input true
   -status     false
COPY solution 20 58-60 # for 60-cell model   
COPY exchange 20 58-60 # for 60-cell model   
END
#
#   20-cell model, transport
#
TRANSPORT   # Diffuse 24C, NaCl solution from column ends
   -cells   20
   -shifts  1
   -flow_direction diffusion
   -boundary_conditions  constant  closed
   -lengths           20*1.0
   -thermal_diffusion 3.0        # Heat is retarded equal to Na
   -dispersivities    20*0.0        # No dispersion
   -diffusion_coefficient 0.3e-9 # m^2/s
   -time_step         1.0e+10    # 317 years, 19 substeps will be used
SELECTED_OUTPUT
   -file           ex12a.sel
   -high_precision true
   -reset          false
   -distance       true
   -temperature    true
USER_PUNCH
  -head Na_mmol K_mmol Cl_mmol error_Cl  error_Na
  10  PUNCH TOT("Na")*1000, TOT("K")*1000, TOT("Cl")*1000
#
#   Calculate deviation from analytical solution for Cl and Na
#
  20  x = DIST
  30  if (x > 8.5 OR SIM_TIME <= 0) THEN END
  40  IF (ABS(x MOD 0.5) > 1e-3) OR (TC <= 0) THEN END
  50  DATA 0.254829592, -0.284496736, 1.421413741, -1.453152027, 1.061405429, 0.3275911
  60  READ a1, a2, a3, a4, a5, a6
  70      REM  calculate error in Cl
  80  z = x / (2*SQRT(3e-10 * SIM_TIME / 1.0))
  90  GOSUB 2000 # erfc(z)...
  100 erfc_Cl = erfc
  110     REM  calculate error in Na, 3 times retarded
  120 z = z * SQRT(3.0)
  130 GOSUB 2000 # erfc(z)...
  140 erfc_Na = erfc
  150     REM punch results
  160 error_Cl = 0.024 * erfc_Cl - TOT("Cl")
  170 error_Na = 0.024 * erfc_Na - TOT("Na")
  180 PUNCH error_Cl, error_Na
  190     REM store results
  200 j = x - 0.5
  210 PUT(error_Cl, SIM_NO, j, 1) 
  220 PUT(error_Na, SIM_NO, j, 2) 
  500 END
  2000 REM calculate erfc...
  2050 b = 1 / (1 + a6 * z)
  2060 erfc = b * (a1 + b * (a2 + b * (a3 + b * (a4 + b * a5)))) * EXP(-(z * z))
  2080 RETURN

END
#
#   60-cell model, initial conditions, when not copied from 20-cell model
#
SELECTED_OUTPUT
   -active        false     # See also PRINT; selected_output false
SOLUTION 1-57   24.0 mM KNO3
    units mol/kgw
    temp  0                 # Incoming solution 0C
    pH    7.0
    pe   12.0   O2(g) -0.67
    K    24.e-3
    N(5) 24.e-3
EXCHANGE 1-57
    KX    0.048
END
#
#   60-cell model, transport
#
TRANSPORT  # no need to redefine parameters that don't change from 20 cell model
   -cells        60
   -lengths      60*0.33333333333333333
   -disp         60*0.0
   -punch_cells  1-60
SELECTED_OUTPUT
   -active        true      # See also PRINT; selected_output false
END
#
#   Print comparison with analytical solution for 
#   Cl and Na in 20-cell and 60-cell models
#
SOLUTION                    # Initial solution calculation for pure water
                            # A calculation is needed to invoke USER_PRINT
PRINT
   -reset      false        # Initial solution calculation not printed
   -echo_input true
   -user_print true
   -high_precision false    # formerly USER_PRINT precision controlled by SELECTED_OUTPUT
USER_PRINT
  10  PRINT "             Error in Cl concentration Error in Na concentration"
  20  PRINT "             ------------------------- -------------------------"    
  30  PRINT "    Distance      20-cell      60-cell      20-cell      60-cell"
  40  PRINT " "
  50  FOR j = 0 TO 8
  60    PRINT j + 0.5, GET(2, j, 1), GET(4, j, 1), GET(2, j, 2), GET(4, j, 2)
  70  NEXT j
END