Hemt ph.tcl

   #simulation parameters
   pdbSetDouble Math iterLimit 1000              ;#sets the max number of Newton iterations to 1000
   math device dim=2 col umf none scale          
   #-----------------------------------------------------------------------------------------------------------------------------------------
   #----------------------------------Include Source Files-----------------------------------------------------------------------------------
   # include files    
   source GaN_modelfile_mo.tcl              ;# this file add the materials in the structure and specifies their material properties and specifies Ohmic contacts for S D B and Schottky for G
   source Poisson.tcl
   source Continuity.tcl
   #-----------------------------------------------------------------------------------------------------------------------------------------
   #----------------------------------Create Structure---------------------------------------------------------------------------------------
   #proc called Struct2D that defines grid (default unit is microns)
   #gate Length increment
   set gateLinc 0.0 ; #a 0.15 increment is a +0.30um growth with gate length of 1.3um. 
  proc Struct2D {} {                            
   if {1} {                                  ;# the [if {1}] tcl statement allows one to comment out the text between the brackets if {1} is set to {0}. This is nice to use if you    want to make multiple structures and easily change between them.
   line x loc=-0.3 spac=0.05 tag=NTop
   line x loc=-0.0005 spac=0.0001 tag=Oxtop
   line x loc=0.0 spac=0.0001 tag=AlGaNTop
   line x loc=0.015 spac=0.0008 tag=AlGaNBottom
   line x loc=1.8 spac=0.2 tag=GaNBottom
   line x loc=2.0 spac=0.02 tag=AlNBottom
   line x loc=3.0 spac=1.5 tag=BBottom
   }
   if {1} {        ;# 1 um gate (no T gate)
   #gate Length increment
   set gateLinc 0.0
   line y loc=(-2.0-$gateLinc) spac=0.1 tag=Left 
   line y loc=(-1.25-$gateLinc) spac=0.5 
   line y loc=(-0.50-$gateLinc) spac=0.01 tag=Oxleft
   line y loc=(-0.15-$gateLinc) spac=0.1
   line y loc=0.0  spac=0.10
   line y loc=(0.15+$gateLinc) spac=0.1
   line y loc=(0.50+$gateLinc) spac=0.01 tag=Oxright
   line y loc=(1.25+$gateLinc) spac=0.5
   line y loc=(2.0+$gateLinc) spac=0.1 tag=Right
   }
    #use the tags above to create the material regions
	#Bulk
   region SiC xlo=AlNBottom xhi=BBottom ylo=Left yhi=Right
   #thin AlN layer
   region AlN xlo=GaNBottom xhi=AlNBottom ylo=Left yhi=Right
   #Buffer
   region GaN xlo=AlGaNBottom xhi=GaNBottom ylo=Left yhi=Right
   #AlGaN under gate
   region AlGaN xlo=AlGaNTop xhi=AlGaNBottom ylo=Left yhi=Right
   #DL
   region GaO xlo=Oxtop xhi=AlGaNTop ylo=Oxleft yhi=Oxright
   #gate metal/electrolyte
   region Electrolyte xlo=NTop xhi=Oxtop ylo=Oxleft yhi=Oxright
   #passivation layers
   region Nitride xlo=NTop xhi=AlGaNTop ylo=Left yhi=Oxleft
   region Nitride xlo=NTop xhi=AlGaNTop ylo=Oxright yhi=Right
   #initialize the grid (this command creates the grid with the specified regions)
   init                 
   #Create the contacts 
   contact name=G Electrolyte xlo=-0.40 xhi=-0.28 ylo=(-0.75-$gateLinc) yhi=(0.75+$gateLinc) add depth=1.0 
   contact name=B SiC xlo=2.9 xhi=7.0 add depth=1.0
   contact name=S AlGaN ylo=(-3.4-$gateLinc) yhi=(-1.99-$gateLinc) xlo=-1.5 xhi=0.0149 add depth=1.0 
   contact name=D AlGaN ylo=(1.99+$gateLinc) yhi=(3.4+$gateLinc) xlo=-1.5 xhi=0.0149 add depth=1.0        
   # the x and y pts need to be a little larger than the grid specifications to capture all of the grid points
   # you can specify the depth of the contact into the page by depth=X. With depth=1.0 the current units are A/um. With the depth specified, the current unit will be in A.
   } 
   #end Struct 2D procedure (remember you still need to call the procedure)
   Struct2D   
   #calls the procedure
   #plot the grid and contacts to make sure everything is correct
   plot.2d bound grid
   plot.2d contact=G !cle
   plot.2d contact=B !cle
   plot.2d contact=S !cle
   plot.2d contact=D !cle
   
   #----------------------------------------------------------------------------------------------------------------------------------------
   #--------------------------------------------Specify Solution Variables------------------------------------------------------------------
   # we are solving for Qfn , Qfp, and DevPsi (device potential), temp is a constant and set to 300 K
   solution add name=DevPsi solve negative damp continuous pde   	;# potential in semiconductor
   solution add name=Qfp solve negative damp continuous pde
   solution add name=Qfn solve negative damp continuous pde

   solution add name=Cl solve !negative damp pde         			;# Cl ion concentration
   solution add name=H solve !negative damp pde         			;# hydrogen ion concentration
   solution add name=Na solve !negative damp pde         			;# Na ion concentration
   solution add name=A solve !negative damp pde         			;# acid species concentration
   solution add name=SOH2 solve !negative damp pde					;#positive interface charge

   solution add name=Temp const val=300.0
   #solution add name=Temp solve pde !negative continuous damp

   # these are needed if donor or acceptor traps are included
   solution add name=ETemp solve const val = 300.0 continuous
   solution add name=HTemp solve const val = 300.0 continuous    
   solution add name=Donor solve const val = 10.0
   solution add name=Acceptor solve const val = 10.0

   #----------------------------------------------------------------------------------------------------------------------------------------
   #-------------------------------------------Specify the Doping (this includes static traps)----------------------------------------------
   # the GaN doping is most important when fitting to an experimental IV curve; it is usually p-type
   #GaN Doping units are /cm3, negative sign is p-type, 
   sel z= 2.0e14*Mater(GaN) name=GaN_Doping 
   
   #AlGaN Doping
   sel z= 1.0e0*Mater(AlGaN) name=AlGaN_Doping
   
   #AlN doping
   sel z= 1.0e12*Mater(AlN) name=AlN_Doping

   #SiC doping
   sel z= 1.0e12*Mater(SiC) name=SiC_Doping
   
   #GaN2 doping, GaN cap layer (no physics included)
   sel z= 1.0e0*Mater(GaN2) name=GaN2_Doping

   #for new gate lengths
   set SDedge [expr (1.9+$gateLinc)]

   #Source and Drain contact doping-from contact to 2DEG to make contacts ohmic (Gaussian profile is used in the expression)
   sel z=(1e19*(y>$SDedge)+(y<=$SDedge)*1.0e19*exp(-(y-$SDedge)*(y-$SDedge)/($SDedge*0.02*0.02)))*(exp(-(x*x)/($SDedge*0.03*0.03))) name=Drain_Doping
   sel z=(1e19*(y<-$SDedge)+(y>=-$SDedge)*1.0e19*exp(-(y+$SDedge)*(y+$SDedge)/($SDedge*0.02*0.02)))*(exp(-(x*x)/($SDedge*0.03*0.03))) name=Source_Doping

   #Total doping
   sel z=GaN_Doping+AlGaN_Doping+Drain_Doping+Source_Doping+SiC_Doping+AlN_Doping+GaN2_Doping name=Doping 
   
   #----------------------------------------------------------------------------------------------------------------------------------------  
   #------------------------------------------Specify the AlN ratio in the AlGaN------------------------------------------------------------
   # this is used in the modelfile to calculate AlGaN material parameters
   sel z=0.25 name=AlN_Ratio

   #----------------------------------------------------------------------------------------------------------------------------------------
   #------------------------------------------Species Concentrations and Mobilities---------------------------------------------------------
   #concentrations
   set Av 6.02e23    						;# Avogadro's number, ions/mol
   set cCl  [expr (100.0e-3*$Av)*1e-3]  	;# M*ions/mol*l/cm3 = ions/cm3
   set cNa  [expr (100.0e-3*$Av)*1e-3]  	;# M*ions/mol*m3/cm3 = ions/cm3
   set cH   [expr (1.0e-1*$Av)*1e-3]  		;# ions/cm3
   set cA   [expr (1.0e-1*$Av)*1e-3]  		;# ions/cm3
   set Ns 	 1.0e15						 	;# surface interface site density /cm2

   #mobilities and diffusivity constants 
   set D_Na 1.96e-5                    	;# cm2/s from Lopreore
   set D_Cl 2.03e-5                    	;# cm2/s from Lopreore
   set D_H  8.24e-5                    	;# cm2/s from Lopreore
   set D_A  2.00e-5                    	;# cm2/s from Lopreore
   set Namob [expr {$D_Na/$Vt}]
   set Clmob [expr {$D_Cl/$Vt}]
   set Hmob [expr {$D_H/$Vt}]
   set Amob [expr {$D_A/$Vt}]

   #assign mobility variable for each ion species
   pdbSetDouble Electrolyte Cl mob $Clmob
   pdbSetDouble Electrolyte Na mob $Namob
   pdbSetDouble Electrolyte H mob $Hmob 
   pdbSetDouble Electrolyte A mob $Amob 

   #-----------------------------------------------------------------------------------------------------------------------------------------	
   #-------------------------------------------Poisson in Electrolyte------------------------------------------------------------------------
   proc Poisson_Electro {Mat} {   
   global k q eps0 Vt 
   
   pdbSetDouble $Mat DevPsi DampValue 0.025
   pdbSetDouble $Mat DevPsi Abs.Error 0.001
   pdbSetDouble $Mat DevPsi Rel.Error 0.01
   
   set eqn " ($eps0 * [pdbDelayDouble $Mat DevPsi RelEps] * grad(DevPsi) / $q) + (- Cl + H - A + Na)"
   pdbSetString $Mat DevPsi Equation $eqn
   } ;#end Poisson proc
   #-----------------------------------------------------------------------------------------------------------------------------------------
   #------------------------------------------Continuity in Electrolyte----------------------------------------------------------------------
   proc Continuity2 {Mat} {
   global Vt Namob Clmob Hmob Amob

   #set eqn "ddt($species) - (([pdbDelayDouble $Mat $species mob]) * sgrad($species, -DevPsi / $Vt))"
   set eqnNa "ddt(Na) - ($Namob) * sgrad(Na, -DevPsi/$Vt)"
   set eqnCl "ddt(Cl) - ($Clmob) * sgrad(Cl, DevPsi/$Vt)"
   set eqnH "ddt(H) - ($Hmob) * sgrad(H, -DevPsi/$Vt)"
   set eqnA "ddt(A) - ($Amob) * sgrad(A, DevPsi/$Vt)"

   pdbSetString $Mat Na Equation $eqnNa
   pdbSetDouble $Mat Na Abs.Error 1.0e-5
   pdbSetDouble $Mat Na Rel.Error 1.0e-3
   pdbSetString $Mat Cl Equation $eqnCl
   pdbSetDouble $Mat Cl Abs.Error 1.0e-5
   pdbSetDouble $Mat Cl Rel.Error 1.0e-3
   pdbSetString $Mat H Equation $eqnH
   pdbSetDouble $Mat H Abs.Error 1.0e-5
   pdbSetDouble $Mat H Rel.Error 1.0e-3
   pdbSetString $Mat A Equation $eqnA
   pdbSetDouble $Mat A Abs.Error 1.0e-5
   pdbSetDouble $Mat A Rel.Error 1.0e-3

   }	;#end Continuity proc
   #-----------------------------------------------------------------------------------------------------------------------------------------
   #------------------------------------------Call procedures describing Physics-------------------------------------------------------------
   #Poisson equation and trap ionization
   Poisson GaN2
   Poisson AlGaN
   Poisson_Ins Nitride
   Poisson_Ins GaO
   Poisson SiC
   Poisson AlN
   Poisson GaN
   Poisson_Electro Electrolyte

   #Continuity equations for electron, hole and transient trap simulation
   ElecContinuity GaN
   ElecContinuity AlGaN
   ElecContinuity SiC
   ElecContinuity AlN

   HoleContinuity GaN
   HoleContinuity AlGaN
   HoleContinuity SiC
   HoleContinuity AlN
   Continuity2 Electrolyte			;#call proc for electrolyte continuity 
   #-----------------------------------------------------------------------------------------------------------------------------------------

   #--------------------------------------------Change in Ion Concentrations at Contact Proc-------------------------------------------------
   proc Dir.contact2 {Contact} {
   global cCl cNa cH cA

   pdbSetBoolean $Contact Cl Flux 0
   pdbSetBoolean $Contact Cl Fixed 1
   pdbSetBoolean $Contact Na Flux 0
   pdbSetBoolean $Contact Na Fixed 1
   pdbSetBoolean $Contact H Flux 0
   pdbSetBoolean $Contact H Fixed 1
   pdbSetBoolean $Contact A Flux 0
   pdbSetBoolean $Contact A Fixed 1

   pdbSetString $Contact Cl Equation "Cl-$cCl"
   pdbSetString $Contact Na Equation "Na-$cNa"
   pdbSetString $Contact H Equation "H-$cH"
   pdbSetString $Contact A Equation "A-$cA"
   }	;#end electrolyte ion conc proc

   #-----------------------------------------------------------------------------------------------------------------------------------------
   #---------------------------------------------Initialize all contacts to 0V---------------------------------------------------------------
   #contact name=G Electrolyte voltage supply=3.0
   contact name=B SiC voltage supply=0.0
   contact name=S AlGaN voltage supply=0.0
   contact name=D AlGaN voltage supply=0.0
   
   pdbSetBoolean B Temp Fixed 0
   pdbSetString B Temp Equation "(2.7*(Temp-300.0))/0.00125"

   #-----------------------------------------------------------------------------------------------------------------------------------------
   #----------------------------------------------Run the Initial Guess Procedure------------------------------------------------------------
   #Initial guess procedure, call for all materials 
   proc Init {Mat} {
   newton $Mat eqn=Doping+Donor-Acceptor+Hole-Elec var=DevPsi damp=0.025
   }
   Init GaN
   Init AlGaN
   Init AlN
   Init SiC

   sel z=$cH*Mater(Electrolyte) name=H1
   sel z=0 name=H2
   sel z=H1+H2 name=H
   sel z=$cNa*Mater(Electrolyte) name=Na
   sel z=$cCl*Mater(Electrolyte) name=Cl
   sel z=$cA*Mater(Electrolyte) name=A
   sel z=6.1 name=DevPsi

   Dir.contact2 G			;#call for ion concentration change at contact (electrolyte) proc
   #-----------------------------------------------------------------------------------------------------------------------------------------
   #--------------------------------------------Run and Plot---------------------------------------------------------------------------------

   set WinA [CreateGraphWindow]   	;# absorbed charge v. pH (H+, A- and total)
   set WinB [CreateGraphWindow]	;# semiconductor potential v. pH
   set WinC [CreateGraphWindow]	;# Id-Vds curve (point measurement)
   set WinD [CreateGraphWindow]	;# net adsorbed charge v. position (y-direction) across pH range 
   device													;# multiple device commands to help with convergence
   device
   device 
   device

   pdbSetString AlGaN_GaO DevPsi Equation " -3.26e12"   	;# polarization charge for AlGaN-DL
   pdbSetString AlN_GaN DevPsi Equation "2.2e13"			;# polarization charge for AlN-GaN
   pdbSetString AlGaN_Nitride DevPsi Equation " -3.26e13"  ;# polarization charge for AlGaN-Nitride
   pdbSetString AlGaN_GaN DevPsi Equation  "1.06e13" 		;# polarization charge for AlGaN-GaN

   device init
   device
   device 
   device

   #--------------------------------------------Energy Band Diagram v. Position (X-Direction, DL/AlGaN)--------------------------------------------------------------------------			
   sel z=Econd
   plot.1d y.v=0.0 plot_name=Band label=Ec !cle
   sel z=Eval
   plot.1d y.v=0.0 plot_name=Band label=Ev !cle
   sel z=Qfp
   plot.1d y.v=0.0 plot_name=Band label=Qfp !cle
   sel z=Qfn
   plot.1d y.v=0.0 plot_name=Band label=Qfn !cle
   sel z=DevPsi
   plot.1d y.v=0.0 plot_name=Band label=DevPsi !cle
   sel z=Doping
   plot.1d x.v=0.0 plot_name=DopingX label=Doping !cle
   sel z=Doping
   plot.1d y.v=0.0 plot_name=DopingY label=Doping !cle
   #specify solution variable for negative oxide surface site SO
   solution add name=SO solve !negative damp pde         

   #-------------------------------------------pH increment and adsorbed charge derivation----------------------------------------------------
   if {1} {
   #Kfa is a the forward reaction rate coefficient for the SO equation. 
   set Kfa 2.58867e-11 											
   #Kra is a the reverse reaction rate coefficient for the SO equation. 
   set Kra 1.0e5										
   #Kfb is a the forward reaction rate coefficient for the SOH2 equation. 
   set Kfb 2.58867e-11									
   #Krb is a the reverse reaction rate coefficient for the SOH2 equation. 
   set Krb 1.0e2											
   #Ns is the number of available surface sites
   set Ns 1.0e15

   # H interface equations for Electrolyte-DL 
   pdbSetString Electrolyte_GaO SO Equation "ddt(SO) - $Kra*($Ns-SO-SOH2) + $Kfa*(H(Electrolyte))*(SO)"
   pdbSetString Electrolyte_GaO SOH2 Equation "ddt(SOH2) - $Krb*(SOH2) + $Kfb*(H(Electrolyte))*($Ns-SOH2-SO)"
   pdbSetString Electrolyte_GaO H Electrolyte Equation "+ $Krb*(SOH2) - $Kfb*(H(Electrolyte))*($Ns-SOH2-SO) + $Kra*($Ns-SO-SOH2) - $Kfa*(H(Electrolyte))*(SO)"
   pdbSetString Electrolyte_GaO DevPsi Equation "SOH2 - SO"

   #ph increment
   for {set inc 0} {$inc<10.5} {set inc [expr $inc+1.0]} {
   set mult [expr pow(10,$inc)]    ;# need to put in expr to evaluate right away to avoid syntax error in pdb equation
   pdbSetString G H Equation "H-($cH/$mult)"
   pdbSetString G A Equation "A-($cA/$mult)"

   #drain bias ramp
   for {set Vds 0.05} {$Vds<1.1} {set Vds [expr $Vds+0.1]} {
   contact name=D supply=$Vds
   device init
   device
   device
   set cur [expr (1.0e6*(-[contact name=D sol=Qfn flux] + [contact name=D sol=Qfp flux]))]  ;# plots the current versus time (mA/mm) 
   AddtoLine $WinC IdVd.$inc $Vds $cur 
   }	

   # ion distribution (X-direction, electrolyte/DL/AlGaN), p[ion] = -log10[ion conc]
   sel z=log10(Na+1.0)
   plot.1d y.v=0.0 !cle plot_name=electrolyte label=Na
   sel z=log10(Cl+1.0)
   plot.1d y.v=0.0 !cle plot_name=electrolyte label=Cl
   sel z=log10(H+1.0)
   plot.1d y.v=0.0 !cle plot_name=electrolyte label=H
   sel z=log10(A+1.0)
   plot.1d y.v=0.0 !cle plot_name=electrolyte label=A

   #adsorbed charge derivation for electrolyte-DL 
   sel z=SOH2
   set nit_plus [interface Electrolyte /GaO y.v=0.0 val]
   puts $nit_plus
   sel z=SO
   set nit_minus [interface Electrolyte /GaO y.v=0.0 val]
   puts $nit_minus

   set int_chrg [expr ($nit_plus-$nit_minus)]
   set pH 1+$inc
   AddtoLine $WinA Int_chrg $pH $int_chrg 
   AddtoLine $WinA SOH2 $pH $nit_plus
   AddtoLine $WinA SO $pH $nit_minus  
   #-------------------------------------------Electrostatic Potential and pH (Y-Direction)----------------------------------------------------
   sel z=DevPsi
   set potential [interface Electrolyte /GaO y.v=0.0 val]
   AddtoLine $WinB Potential $pH $potential
   set potential [interface Electrolyte /GaO y.v=-0.25 val]
   AddtoLine $WinB Potential2 $pH $potential
   set potential [interface Electrolyte /GaO y.v=0.25 val]
   AddtoLine $WinB Potential3 $pH $potential

   #-------------------------------------------Net Adsorbed Charge v Position (Y-Direction)----------------------------------------------------	
   #these points are defined for 1 um gate length. For any other gate length, the y.v values must be changed. 
   if {1} {
   #positive oxide surface sites
   sel z=SOH2 
   set t1 [interface Electrolyte /GaO y.v=-0.5 val]
   set t2 [interface Electrolyte /GaO y.v=-0.375 val]
   set t3 [interface Electrolyte /GaO y.v=-0.25 val]
   set t4 [interface Electrolyte /GaO y.v=-0.125 val]
   set t5 [interface Electrolyte /GaO y.v=0.0 val]
   set t6 [interface Electrolyte /GaO y.v=0.125 val]
   set t7 [interface Electrolyte /GaO y.v=0.25 val]
   set t8 [interface Electrolyte /GaO y.v=0.375 val]
   set t9 [interface Electrolyte /GaO y.v=0.4 val]
   set t10 [interface Electrolyte /GaO y.v=0.5 val]
   set t11 [interface Electrolyte /GaO y.v=0.48 val]
   set t12 [interface Electrolyte /GaO y.v=0.5 val]

   sel z=SO 													;#negative oxide surface site
   set m1 [interface Electrolyte /GaO y.v=-0.5 val]
   set m2 [interface Electrolyte /GaO y.v=-0.375 val]
   set m3 [interface Electrolyte /GaO y.v=-0.25 val]
   set m4 [interface Electrolyte /GaO y.v=-0.125 val]
   set m5 [interface Electrolyte /GaO y.v=0.0 val]
   set m6 [interface Electrolyte /GaO y.v=0.125 val]
   set m7 [interface Electrolyte /GaO y.v=0.25 val]
   set m8 [interface Electrolyte /GaO y.v=0.375 val]
   set m9 [interface Electrolyte /GaO y.v=0.40 val]
   set m10 [interface Electrolyte /GaO y.v=0.45 val]
   set m11 [interface Electrolyte /GaO y.v=0.48 val]
   set m12 [interface Electrolyte /GaO y.v=0.5 val]

   set T1 [expr $t1-$m1]
   set T2 [expr $t2-$m2]
   set T3 [expr $t3-$m3]
   set T4 [expr $t4-$m4]
   set T5 [expr $t5-$m5]
   set T6 [expr $t6-$m6]
   set T7 [expr $t7-$m7]
   set T8 [expr $t8-$m8]
   set T9 [expr $t9-$m9]
   set T10 [expr $t10-$m10]
   set T11 [expr $t11-$m11]
   set T12 [expr $t12-$m12]

   AddtoLine $WinD Intchrg.$inc -0.5 $T1
   AddtoLine $WinD Intchrg.$inc -0.375 $T2
   AddtoLine $WinD Intchrg.$inc -0.25 $T3
   AddtoLine $WinD Intchrg.$inc -0.125 $T4
   AddtoLine $WinD Intchrg.$inc 0.0 $T5
   AddtoLine $WinD Intchrg.$inc 0.125 $T6
   AddtoLine $WinD Intchrg.$inc 0.25 $T7
   AddtoLine $WinD Intchrg.$inc 0.375 $T8
   AddtoLine $WinD Intchrg.$inc 0.4 $T9
   AddtoLine $WinD Intchrg.$inc 0.45 $T10
   AddtoLine $WinD Intchrg.$inc 0.48 $T11
   AddtoLine $WinD Intchrg.$inc 0.5 $T12
   }
   #-------------------------------------------SO/DevPsi v. Position Across pH Range----------------------------------------------
   # SO v position (y-direction) 
   sel z=SO
   plot.1d x.v=0.0005 plot_name=Intchrg label=SO.$inc !cle
   set nitT_minus [interface Electrolyte /GaO x.v=0.0005 val]		
   puts $nitT_minus
   sel z=Econd
   plot.1d y.v=0.0 plot_name=Band2 label=Ec.$inc !cle
   sel z=Eval
   plot.1d y.v=0.0 plot_name=Band2 label=Ev.$inc !cle
   sel z=Qfp
   plot.1d y.v=0.0 plot_name=Band2 label=Qfp.$inc !cle
   sel z=Qfn
   plot.1d y.v=0.0 plot_name=Band2 label=Qn.$inc !cle

   #DevPsi v position (x-direction, DL/AlGaN) 
   sel z=DevPsi
   plot.1d y.v=0.0 plot_name=Band2 label=DevPsi.$inc !cle

   #-------------------------------------------Current v. time and Carrier (E) Conc v. Position (X-Direction, DL/AlGaN)----------------------------------------------------------------------------
   plot.1d y.v=0.0 plot_name=elec_conc label=elec.$inc !cle		
   }	;#end for loop, increment pH
   }	;#end if{1}