      SUBROUTINE starte(NHW,NHE,NVN,NVS)
C
C     THIS SUBROUTINE INITIALIZES ALL SOIL CHEMISTRY VARIABLES
C
      include "parameters.h"
      include "filec.h"
      include "files.h"
      include "blkc.h"
      include "blk2a.h"
      include "blk2b.h"
      include "blk2c.h"
      include "blk8a.h"
      include "blk8b.h"
      include "blk11a.h"
      include "blk11b.h"
      include "blk13a.h"
      include "blk13b.h"
      include "blk13c.h"
      include "blk19a.h"
      include "blk19b.h"
      include "blk19c.h"
      include "blk19d.h"
      CHARACTER*16 DATA(30)
      CHARACTER*3 CHOICE(102,20)
      CHARACTER*8 CDATE
      PARAMETER (SCO2X=7.391E-01,SCH4X=3.156E-02,SOXYX=2.925E-02
     2,SN2GX=1.510E-02,SN2OX=5.241E-01,SNH3X=2.852E+02,SH2GX=3.156E-02
     3,ACO2X=0.14,ACH4X=0.14,AOXYX=0.31,AN2GX=0.23,AN2OX=0.23
     4,ANH3X=0.07,AH2GX=0.14)
      PARAMETER (DPH2O=6.5E-09,SPALO=6.5E-22,SPFEO=6.5E-27
     2,SPCAC=3.8E-03,SPCAS=1.4E+01,SPALP=1.0E-15,SPFEP=1.0E-20
     3,SPCAM=7.0E+07,SPCAD=1.0E-01,SPCAH=2.3E-31,SXOH2=4.5E-05
     4,SXOH1=1.1E-06,SXH2P=2.0E+07,SXH1P=2.0E+07
     5,DPCO2=4.2E-04,DPHCO=5.6E-08,DPN4=5.7E-07
     6,DPAL1=4.6E-07,DPAL2=7.3E-07,DPAL3=1.8E-05
     7,DPAL4=1.2E-05,DPALS=0.16,DPFE1=4.6E-07,DPFE2=7.3E-07
     8,DPFE3=1.8E-05,DPFE4=1.2E-05,DPFES=7.1E-02,DPCAO=12.5
     9,DPCAC=4.2E-02,DPCAH=13.5,DPCAS=1.2,DPMGO=0.7,DPMGC=0.3
     1,DPMGH=67.0,DPMGS=2.1,DPNAC=0.45,DPNAS=3.3E+02,DPKAS=5.0E+01
     2,DPH1P=4.5E-10,DPH2P=6.3E-05,DPH3P=7.1,DPF1P=4.5E-02
     3,DPF2P=3.7E-03,DPC0P=3.5E-04,DPC1P=1.82,DPC2P=40.0
     4,DPM1P=1.23,DPCOH=1.0E-02,DPALO=6.3E+04,DPFEO=6.3E+04)
      PARAMETER (DPCO3=DPCO2*DPHCO,SHALO=SPALO/DPH2O**3
     2,SYAL1=SPALO/DPAL1,SHAL1=SYAL1/DPH2O**2,SYAL2=SYAL1/DPAL2
     3,SHAL2=SYAL2/DPH2O,SPAL3=SYAL2/DPAL3,SYAL4=SPAL3/DPAL4
     4,SHAL4=SYAL4*DPH2O,SHFEO=SPFEO/DPH2O**3,SYFE1=SPFEO/DPFE1
     5,SHFE1=SYFE1/DPH2O**2,SYFE2=SYFE1/DPFE2,SHFE2=SYFE2/DPH2O
     6,SPFE3=SYFE2/DPFE3,SYFE4=SPFE3/DPFE4,SHFE4=SYFE4*DPH2O
     7,SHCAC1=SPCAC/DPHCO,SYCAC1=SHCAC1*DPH2O,SHCAC2=SHCAC1/DPCO2
     8,SYCAC2=SHCAC2*DPH2O**2,SHA0P1=SPALP/DPH1P,SYA0P1=SHA0P1*DPH2O
     9,SPA1P1=SYA0P1/DPAL1,SYA2P1=SPA1P1/DPAL2,SHA2P1=SYA2P1*DPH2O
     1,SYA3P1=SYA2P1/DPAL3,SHA3P1=SYA3P1*DPH2O**2,SYA4P1=SYA3P1/DPAL4
     2,SHA4P1=SYA4P1*DPH2O**3,SHA0P2=SHA0P1/DPH2P
     3,SYA0P2=SHA0P2*DPH2O**2,SYA1P2=SYA0P2/DPAL1,SHA1P2=SYA1P2/DPH2O
     4,SPA2P2=SYA1P2/DPAL2,SYA3P2=SPA2P2/DPAL3,SHA3P2=SYA3P2*DPH2O
     5,SYA4P2=SYA3P2/DPAL4,SHA4P2=SYA4P2*DPH2O**2)
      PARAMETER (SHF0P1=SPFEP/DPH1P,SYF0P1=SHF0P1*DPH2O
     2,SPF1P1=SYF0P1/DPFE1,SYF2P1=SPF1P1/DPFE2,SHF2P1=SYF2P1*DPH2O
     3,SYF3P1=SYF2P1/DPFE3,SHF3P1=SYF3P1*DPH2O**2,SYF4P1=SYF3P1/DPFE4
     4,SHF4P1=SYF4P1*DPH2O**3,SHF0P2=SHF0P1/DPH2P,SYF0P2=SHF0P2*DPH2O**2
     5,SYF1P2=SYF0P2/DPFE1,SHF1P2=SYF1P2/DPH2O,SPF2P2=SYF1P2/DPFE2
     6,SYF3P2=SPF2P2/DPFE3,SHF3P2=SYF3P2*DPH2O,SYF4P2=SYF3P2/DPFE4
     7,SHF4P2=SYF4P2*DPH2O**2,SHCAD2=SPCAD/DPH2P,SYCAD2=SHCAD2*DPH2O
     8,SHCAH1=SPCAH/(DPH2O*DPH1P**3),SYCAH1=SHCAH1*DPH2O**4
     9,SHCAH2=SHCAH1/DPH2P**3,SYCAH2=SHCAH2*DPH2O**7)
      PARAMETER (TPD=2.5E-02,TAD=5.0E-02,TSL=0.5,A0=1.0,COOH=2.5E-02
     2,CALMX=10.0,CFEMX=10.0,CCAMX=10.0,MRXN=1000)
C
C     INITIALIZE CATION AND ANION CONCENTRATIONS
C     IN PRECIPITATION (K=1), IRRIGATION (K=2) AND SOIL (K=3)
C     FROM WEATHER, IRRIGATION AND SOIL FILES IN 'READS'
C
      DO 9995 NX=NHW,NHE
      DO 9990 NY=NVN,NVS
      CCO2M=CCO2EI(NY,NX)/12.0
      CCH4M=CCH4E(NY,NX)/12.0
      COXYM=COXYE(NY,NX)/32.0
      CZ2GM=CZ2GE(NY,NX)/14.0
      CZ2OM=CZ2OE(NY,NX)/14.0
      DO 1200 I=1,366
      DO 1200 L=NU(NY,NX),NL(NY,NX)
      DO 2000 K=1,3
C
C     INITIALIZE RAINFALL
C
      IF(K.EQ.1.AND.I.EQ.1.AND.L.EQ.1)THEN
      PH1=PHR(NY,NX)
      CHY1=10.0**(-(PHR(NY,NX)-3.0))
      COH1=DPH2O/CHY1
      CN4Z=CN4R(NY,NX)
      CNOZ=CNOR(NY,NX)
      CPOZ=CPOR(NY,NX)
      CALZ=CALR(NY,NX)
      CFEZ=CFER(NY,NX)
      CCAZ=CCAR(NY,NX)
      CMGZ=CMGR(NY,NX)
      CNAZ=CNAR(NY,NX)
      CKAZ=CKAR(NY,NX)
      CSOZ=CSOR(NY,NX)
      CCLZ=CCLR(NY,NX)
C
C     INITIALIZE IRRIGATION WATER
C
      ELSEIF(K.EQ.2.AND.L.EQ.1)THEN
      PH1=PHQ(I,NY,NX)
      CHY1=10.0**(-(PHQ(I,NY,NX)-3.0))
      COH1=DPH2O/CHY1
      CN4Z=CN4Q(I,NY,NX)
      CNOZ=CNOQ(I,NY,NX)
      CPOZ=CPOQ(I,NY,NX)
      CALZ=CALQ(I,NY,NX)
      CFEZ=CFEQ(I,NY,NX)
      CCAZ=CCAQ(I,NY,NX)
      CMGZ=CMGQ(I,NY,NX)
      CNAZ=CNAQ(I,NY,NX)
      CKAZ=CKAQ(I,NY,NX)
      CSOZ=CSOQ(I,NY,NX)
      CCLZ=CCLQ(I,NY,NX)
C
C     INITIALIZE SOIL WATER
C
      ELSEIF(K.EQ.3.AND.I.EQ.1.AND.DATA(20).EQ.'NO'.AND.IGO.EQ.0)THEN
      IF(BKVL(L,NY,NX).GT.ZEROS(NY,NX))THEN
      BKVLX=BKVL(L,NY,NX)
      ELSE
      BKVLX=VOLA(L,NY,NX)
      ENDIF
      CN4X=CNH4(L,NY,NX)
      CNOX=CNO3(L,NY,NX) 
      CPOX=CPO4(L,NY,NX)
      CALX=CAL(L,NY,NX)
      CFEX=CFE(L,NY,NX)
      CCAX=CCA(L,NY,NX)
      CMGX=CMG(L,NY,NX)
      CNAX=CNA(L,NY,NX)
      CKAX=CKA(L,NY,NX)
      CSOX=CSO4(L,NY,NX)
      CCLX=CCL(L,NY,NX)
      CALPOX=CALPO(L,NY,NX)
      CFEPOX=CFEPO(L,NY,NX)
      CCAPDX=CCAPD(L,NY,NX)
      CCAPHX=CCAPH(L,NY,NX)
      CALOHX=CALOH(L,NY,NX)
      CFEOHX=CFEOH(L,NY,NX)
      CCACOX=CCACO(L,NY,NX)
      CCASOX=CCASO(L,NY,NX)
      XCEC(L,NY,NX)=AMAX1(CN4X,CEC(L,NY,NX))*BKVLX
      XAEC(L,NY,NX)=AMAX1(CPOX,AEC(L,NY,NX))*BKVLX
      CHY1=10.0**(-(PH(L,NY,NX)-3.0))
      COH1=DPH2O/CHY1
      CN4Z=CN4X
      CNOZ=CNOX
      CPOZ=CPOX
      CALZ=CALX
      CFEZ=CFEX
      CCAZ=CCAX
      CMGZ=CMGX
      CNAZ=CNAX
      CKAZ=CKAX
      CSOZ=CSOX
      CCLZ=CCLX
      ELSE
      GO TO 2000
      ENDIF
C
C     INITIALIZE SOLUTE EQUILIBRIA
C
      CC3=AMAX1(0.0,CALZ)+AMAX1(0.0,CFEZ)
      CA3=0.0
      CC2=AMAX1(0.0,CCAZ)+CMGZ
      CA2=CSOZ
      CC1=CN4Z+CNAZ+CKAZ
      CA1=CNOZ+CCLZ+CPOZ
      CX2=0.0
      CX1=0.0
      CN=0.0
C
C     INITIALIZE ION STRENGTH AND ACTIVITIES
C
      CION2=AMAX1(0.0,CC3+CA3+CC2+CA2+CC1+CA1+CN)
      CSTR1=0.5E-03*(9.0*(CC3+CA3)+4.0*(CC2+CA2)+CC1+CA1)
      CSTRZ=0.5E-03*(9.0*(CC3+CA3)+4.0*(CC2+CX2)+CC1+CX1)
      CSTR2=SQRT(CSTR1)
      FSTR2=CSTR2/(1.0+CSTR2)
      FH2O=5.56E+04/(5.56E+04+CION2)
      IF(ISALTG.NE.0)THEN
      A1=AMIN1(1.0,10.0**(-0.509*1.0*FSTR2+0.20*CSTR2))
      A2=AMIN1(1.0,10.0**(-0.509*4.0*FSTR2+0.20*CSTR2))
      A3=AMIN1(1.0,10.0**(-0.509*9.0*FSTR2+0.20*CSTR2))
      ELSE
      A1=1.0
      A2=1.0
      A3=1.0
      ENDIF
C
C     INITIALIZE GASES
C
      CCO2X=CCO2M*SCO2X/(EXP(ACO2X*CSTRZ))
     2*EXP(0.843-0.0281*ATCA(NY,NX))*FH2O
      CCO2Y=LOG(CCO2X)
      CCO2Z=ABS(CCO2Y)
      CCO21=CCO2X
      FCO3=DPCO3*A0/(CHY1**2*A2)
      FHCO=DPCO2*A0/(CHY1*A1)
      Z=ACO2X*(2.0E-03*FCO3+0.5E-03*FHCO)
      DO 3020 MM=1,25
      R=(LOG(CCO21)+Z*CCO21-CCO2Y)/CCO2Z
      IF(R.LT.1.0E-03)GO TO 3030
      CCO21=CCO21/SQRT(1.0+R)
3020  CONTINUE
3030  CONTINUE
      CCH41=CCH4M*SCH4X/(EXP(ACH4X*CSTR1))
     2*EXP(0.597-0.0199*ATCA(NY,NX))*FH2O
      COXY1=COXYM*SOXYX/(EXP(AOXYX*CSTR1))
     2*EXP(0.516-0.0172*ATCA(NY,NX))*FH2O
      CZ2G1=CZ2GM*SN2GX/(EXP(AN2GX*CSTR1))
     2*EXP(0.456-0.0152*ATCA(NY,NX))*FH2O
      CZ2O1=CZ2OM*SN2OX/(EXP(AN2OX*CSTR1))
     2*EXP(0.897-0.0299*ATCA(NY,NX))*FH2O
      CCO31=CCO21*DPCO3*A0/(CHY1**2*A2)
      CHCO31=CCO21*DPCO2*A0/(CHY1*A1)
      CNO1=CNOZ
C
C     INITIALIZE ION PAIR EQUILIBRIA
C
      IF(K.NE.3)THEN
      CN41=CN4Z/(1.0+DPN4*A1/(CHY1*A0))
      CN31=CN41*DPN4*A1/(CHY1*A0)
      ELSE
      CN41=ZERO
      CN31=ZERO
      ENDIF
      IF(CALZ.LT.0.0)THEN
      CAL1=AMIN1(CALMX,SPALO/(COH1**3*A3))
      ELSE
      CAL1=AMIN1(CALZ,SPALO/(COH1**3*A3))
      ENDIF
      IF(CFEZ.LT.0.0)THEN
      CFE1=AMIN1(CFEMX,SPFEO/(COH1**3*A3))
      ELSE
      CFE1=AMIN1(CFEZ,SPFEO/(COH1**3*A3))
      ENDIF
      IF(CCAZ.LT.0.0)THEN
      CCA1=AMIN1(CCAMX,SPCAC/(CCO31*A2**2))
      ELSE
      CCA1=AMIN1(CCAZ,SPCAC/(CCO31*A2**2))
      ENDIF
      CMG1=CMGZ
      CNA1=CNAZ
      CKA1=CKAZ
      CSO41=CSOZ
      CCL1=CCLZ
      CALO1=CAL1*COH1*A3/(DPAL1*A2)
      CALO2=CAL1*COH1**2*A3/(DPAL1*DPAL2*A1)
      CALO3=CAL1*COH1**3*A3/(DPAL1*DPAL2*DPAL3*A0)
      CALO4=CAL1*COH1**4*A3/(DPAL1*DPAL2*DPAL3*DPAL4*A1)
      CALS1=0.0
      CFEO1=CFE1*COH1*A3/(DPFE1*A2)
      CFEO2=CFE1*COH1**2*A3/(DPFE1*DPFE2*A1)
      CFEO3=CFE1*COH1**3*A3/(DPFE1*DPFE2*DPFE3*A0)
      CFEO4=CFE1*COH1**4*A3/(DPFE1*DPFE2*DPFE3*DPFE4*A1)
      CFES1=0.0
      CCAO1=CCA1*COH1*A2/(DPCAO*A1)
      CCAC1=CCA1*CCO31*A2**2/(DPCAC*A0)
      CCAH1=CCA1*CHCO31*A2/DPCAH
      CCAS1=0.0
      CMGO1=CMG1*COH1*A2/(DPMGO*A1)
      CMGC1=CMG1*CCO31*A2**2/(DPMGC*A0)
      CMGH1=CMG1*CHCO31*A2/DPMGH
      CMGS1=0.0
      CNAC1=CNA1*CCO31*A2/DPNAC
      CNAS1=0.0
      CKAS1=0.0
      CF1P1=0.0
      CF2P1=0.0
      CC0P1=0.0
      CC1P1=0.0
      CC2P1=0.0
      CM1P1=0.0
C
C     INITIALIZE PHOSPHORUS EQUILIBRIA AMONG SOLUBLE, ADSORBED
C     AND PRECIPITATED FORMS
C
      IF(K.NE.3)THEN
      CH3P1=CPOZ/(1.0+DPH3P*A0/(CHY1*A1)+DPH3P*DPH2P*A0
     2/(CHY1**2*A2)+DPH3P*DPH2P*DPH1P*A0/(CHY1**3*A3))
      CH2P1=CH3P1*DPH3P*A0/(CHY1*A1)
      CH1P1=CH3P1*DPH3P*DPH2P*A0/(CHY1**2*A2)
      CH0P1=CH3P1*DPH3P*DPH2P*DPH1P*A0/(CHY1**3*A3)
      ELSE
      XHP=CPOZ 
      XOH=XAEC(L,NY,NX)/BKVLX
      FHP3=1.0/(1.0+DPH3P*A0/(CHY1*A1)+DPH3P*DPH2P*A0
     2/(CHY1**2*A2)+DPH3P*DPH2P*DPH1P*A0/(CHY1**3*A3))
      FHP2=FHP3*DPH3P*A0/(CHY1*A1)
      FHP1=FHP3*DPH3P*DPH2P*A0/(CHY1**2*A2)
      FHP0=FHP3*DPH3P*DPH2P*DPH1P*A0/(CHY1**3*A3)
      SPOH2=SXOH2/A1
      SPOH1=SXOH1/A1
      SPH2P=SXH2P*DPH2O/A1
      SPH1P=SXH1P*DPH2O*A1/A2
      FXH2=1.0/(1.0+SPOH2/CHY1+SPOH2*SPOH1/CHY1**2)
      FXH1=FXH2*SPOH2/CHY1
      FXH0=FXH1*SPOH1/CHY1
      FXP2=1.0/(1.0+SXH2P*DPH2P/(SXH1P*CHY1))
      FXP1=FXP2*SXH2P*DPH2P/(SXH1P*CHY1)
      FHPA=FHP2*A1
      XPT=(XOH+XHP+SXH2P*FXP2*COH1/(FXH1*FHPA)-SQRT(XOH**2*FXH1**2
     2*FHPA**2-2.0*XOH*FXH1**2*XHP*FHPA**2+FXH1**2*XHP**2*FHPA**2
     3+2.0*XOH*FXH1*FHPA*SXH2P*FXP2*COH1+2.0*FXH1*XHP*FHPA*SXH2P
     4*FXP2*COH1+SXH2P**2*FXP2**2*COH1**2)/(FXH1*FHPA))/2.0
      XOH21=(XOH-XPT)*FXH2
      XOH11=(XOH-XPT)*FXH1
      XOH01=(XOH-XPT)*FXH0
      XH1P1=XPT*FXP1
      XH2P1=XPT*FXP2
      CH3P1=(XHP-XPT)*FHP3
      CH2P1=(XHP-XPT)*FHP2
      CH1P1=(XHP-XPT)*FHP1
      CH0P1=(XHP-XPT)*FHP0
C
C     INITIALIZE CATION EQILIBRIA BETWEEN SOLUBLE
C     AND EXCHANGEABLE FORMS
C
      XCECQ=AMAX1(CN4X,CEC(L,NY,NX))
      XN4Q=CN4X 
      XHYQ=0.0
      XALQ=0.0
      XFEQ=0.0
      XCAQ=0.0
      XMGQ=0.0
      XNAQ=0.0
      XKAQ=0.0
      XHC1=0.0
      XALO21=0.0
      XFEO21=0.0
      XCOOH=AMAX1(0.0,COOH*ORGC(L,NY,NX))
      ENDIF
      CC3=CAL1+CFE1
      CA3=CH0P1
      CC2=CCA1+CMG1+CALO1+CFEO1+CF2P1
      CA2=CSO41+CCO31+CH1P1
      CC1=CN41+CHY1+CNA1+CKA1+CALO2+CFEO2+CALS1+CFES1+CCAO1
     2+CCAH1+CMGO1+CMGH1+CF1P1+CC2P1
      CA1=CNO1+COH1+CHCO31+CCL1+CALO4+CFEO4+CNAC1+CNAS1+CKAS1
     2+CH2P1+CC0P1
      CN=CCO21+CCH41+COXY1+CZ2G1+CZ2O1+CN31+CALO3+CFEO3+CCAC1+CCAS1
     2+CMGC1+CMGS1+CH3P1+CC1P1+CM1P1
      CX2=CA2-CCO31
      CX1=CA1-CHCO31
C
C     INITIALIZE EQUILIBRIA BETWEEN SOLUBLE AND PRECIPITATED FORMS
C
      IF(K.EQ.3)THEN
      PALOH1=CALOHX
      PFEOH1=CFEOHX
      PCACO1=CCACOX
      PCASO1=CCASOX
      PALPO1=CALPOX*VLPO4(L,NY,NX)
      PFEPO1=CFEPOX*VLPO4(L,NY,NX)
      PCAPD1=CCAPDX*VLPO4(L,NY,NX)
      PCAPH1=CCAPHX*VLPO4(L,NY,NX)
      CCEC=AMAX1(ZERO,XCEC(L,NY,NX)/BKVLX)
      CALX=CAL1**0.333
      CFEX=CFE1**0.333
      CCAX=CCA1**0.500
      CMGX=CMG1**0.500
      XCAX=CCEC/(1.0+GKC4(L,NY,NX)*CN41/CCAX
     3+GKCH(L,NY,NX)*CHY1/CCAX+GKCA(L,NY,NX)*CALX/CCAX
     4+GKCA(L,NY,NX)*CFEX/CCAX+GKCA(L,NY,NX)*CMG1/CCAX
     5+GKCN(L,NY,NX)*CNA1/CCAX+GKCK(L,NY,NX)*CKA1/CCAX)
      XN4Q=CN4X
      XHYQ=XCAX*GKCH(L,NY,NX) 
      XALQ=XCAX*GKCA(L,NY,NX) 
      XFEQ=XCAX*GKCA(L,NY,NX) 
      XCAQ=XCAX*CCAX
      XMGQ=XCAX*GKCM(L,NY,NX) 
      XNAQ=XCAX*GKCN(L,NY,NX) 
      XKAQ=XCAX*GKCK(L,NY,NX) 
      XTLQ=XN4Q+XHYQ+XALQ+XFEQ+XCAQ+XMGQ+XNAQ+XKAQ
      IF(XTLQ.GT.ZERO)THEN
      FX=CCEC/XTLQ
      ELSE
      FX=0.0
      ENDIF
      XN41=CN4X
      XHY1=FX*XHYQ
      XAL1=FX*XALQ/3.0
      XFE1=FX*XFEQ/3.0
      XCA1=FX*XCAQ/2.0
      XMG1=FX*XMGQ/2.0
      XNA1=FX*XNAQ
      XKA1=FX*XKAQ
C     WRITE(*,2222)'XN4S',K,L,CN4X,XN4Q,XN41,CNH4(L,NY,NX)
C    2,BKVL(L,NY,NX),XCECQ,XCEC(L,NY,NX),BKVLX,FX
C     WRITE(*,2222)'XOH2S',K,L,XOH21,XOH11,XOH01,XH2P1,XH1P1,XOH,XPT
C    2,CH3P1,CH2P1,CH1P1,CH0P1,CHY1,XHP,FHP3,FHP2,FHP1,FHP0
C    3,FXP2,FXP1,PALPO1,PFEPO1,PCAPD1,PCAPH1 
2222  FORMAT(A8,2I4,40E12.4)
      ENDIF
C
C     CONVERGE TOWARDS ALL SOLUBILITY EQUILIBRIA
C     IF SALT OPTION IS SELECTED
C
      IF(ISALTG.NE.0)THEN
      DO 1000 M=1,MRXN
      CCO21=AMAX1(ZERO,CCO21)
      CCO31=CCO21*DPCO3*A0/(CHY1**2*A2)
      CHCO31=CCO21*DPCO2*A0/(CHY1*A1)
      CN41=AMAX1(ZERO,CN41)
      CN31=AMAX1(ZERO,CN31)
      CAL1=AMAX1(ZERO,CAL1)
      CFE1=AMAX1(ZERO,CFE1)
      CCA1=AMAX1(ZERO,CCA1)
      CCA1=AMIN1(CCA1,SPCAC/(CCO31*A2**2))
      CMG1=AMAX1(ZERO,CMG1)
      CNA1=AMAX1(ZERO,CNA1)
      CKA1=AMAX1(ZERO,CKA1)
      CSO41=AMAX1(ZERO,CSO41)
      CALO1=AMAX1(ZERO,CALO1)
      CALO2=AMAX1(ZERO,CALO2)
      CALO3=AMAX1(ZERO,CALO3)
      CALO4=AMAX1(ZERO,CALO4)
      CALS1=AMAX1(ZERO,CALS1)
      CFEO1=AMAX1(ZERO,CFEO1)
      CFEO2=AMAX1(ZERO,CFEO2)
      CFEO3=AMAX1(ZERO,CFEO3)
      CFEO4=AMAX1(ZERO,CFEO4)
      CFES1=AMAX1(ZERO,CFES1)
      CCAO1=AMAX1(ZERO,CCAO1)
      CCAC1=AMAX1(ZERO,CCAC1)
      CCAH1=AMAX1(ZERO,CCAH1)
      CCAS1=AMAX1(ZERO,CCAS1)
      CMGO1=AMAX1(ZERO,CMGO1)
      CMGC1=AMAX1(ZERO,CMGC1)
      CMGH1=AMAX1(ZERO,CMGH1)
      CMGS1=AMAX1(ZERO,CMGS1)
      CNAC1=AMAX1(ZERO,CNAC1)
      CNAS1=AMAX1(ZERO,CNAS1)
      CKAS1=AMAX1(ZERO,CKAS1)
      CH0P1=AMAX1(ZERO,CH0P1)
      CH1P1=AMAX1(ZERO,CH1P1)
      CH2P1=AMAX1(ZERO,CH2P1)
      CH3P1=AMAX1(ZERO,CH3P1)
      CF1P1=AMAX1(ZERO,CF1P1)
      CF2P1=AMAX1(ZERO,CF2P1)
      CC0P1=AMAX1(ZERO,CC0P1)
      CC1P1=AMAX1(ZERO,CC1P1)
      CC2P1=AMAX1(ZERO,CC2P1)
      CM1P1=AMAX1(ZERO,CM1P1)
C
C     ION ACTIVITY COEFFICIENTS
C
      CC3=CAL1+CFE1
      CA3=CH0P1
      CC2=CCA1+CMG1+CALO1+CFEO1+CF2P1
      CA2=CSO41+CCO31+CH1P1
      CC1=CN41+CHY1+CNA1+CKA1+CALO2+CFEO2+CALS1+CFES1+CCAO1
     2+CCAH1+CMGO1+CMGH1+CF1P1+CC2P1
      CA1=CNO1+COH1+CHCO31+CCL1+CALO4+CFEO4+CNAC1+CNAS1+CKAS1
     2+CH2P1+CC0P1
      CN=CCO21+CCH41+COXY1+CZ2G1+CZ2O1+CN31+CALO3+CFEO3+CCAC1+CCAS1
     2+CMGC1+CMGS1+CH3P1+CC1P1+CM1P1
      CX2=CA2-CCO31
      CX1=CA1-CHCO31
      CION2=AMAX1(0.0,CC3+CA3+CC2+CA2+CC1+CA1+CN)
      CSTR1=0.5E-03*(9.0*(CC3+CA3)+4.0*(CC2+CA2)+CC1+CA1)
      CSTRZ=0.5E-03*(9.0*(CC3+CA3)+4.0*(CC2+CX2)+CC1+CX1)
      CSTR2=SQRT(CSTR1)
      FSTR2=CSTR2/(1.0+CSTR2)
      FH2O=5.56E+04/(5.56E+04+CION2)
      A1=AMIN1(1.0,10.0**(-0.509*1.0*FSTR2+0.20*CSTR2))
      A2=AMIN1(1.0,10.0**(-0.509*4.0*FSTR2+0.20*CSTR2))
      A3=AMIN1(1.0,10.0**(-0.509*9.0*FSTR2+0.20*CSTR2))
C
C     PRECIPITATION-DISSOLUTION EQUILIBRIA
C
      AHY1=CHY1*A1
      AOH1=COH1*A1
      AAL1=CAL1*A3
      AALO1=CALO1*A2
      AALO2=CALO2*A1
      AALO3=CALO3
      AALO4=CALO4*A1
      AFE1=CFE1*A3
      AFEO1=CFEO1*A2
      AFEO2=CFEO2*A1
      AFEO3=CFEO3
      AFEO4=CFEO4*A1
      ACA1=CCA1*A2
      ACO31=CCO31*A2
      AHCO31=CHCO31*A1
      ACO21=CCO21*A0
      ASO41=CSO41*A2
      AH0P1=CH0P1*A3
      AH1P1=CH1P1*A2
      AH2P1=CH2P1*A1
      AH3P1=CH3P1*A0
      AF1P1=CF1P1*A2
      AF2P1=CF2P1*A1
      AC0P1=CC0P1*A1
      AC1P1=CC1P1*A0
      AC2P1=CC2P1*A1
      AM1P1=CM1P1*A0
      AN41=CN41*A1
      AN31=CN31*A0
      AMG1=CMG1*A2
      ANA1=CNA1*A1
      AKA1=CKA1*A1
      AALX=AAL1**0.333
      AFEX=AFE1**0.333
      ACAX=ACA1**0.500
      AMGX=AMG1**0.500
      AALS1=CALS1*A1
      AFES1=CFES1*A1
      ACAO1=CCAO1*A1
      ACAC1=CCAC1*A0
      ACAS1=CCAS1*A0
      ACAH1=CCAH1*A1
      AMGO1=CMGO1*A1
      AMGC1=CMGC1*A0
      AMGH1=CMGH1*A1
      AMGS1=CMGS1*A0
      ANAC1=CNAC1*A1
      ANAS1=CNAS1*A1
      AKAS1=CKAS1*A1
C
C     ALUMINUM HYDROXIDE (GIBBSITE)
C
      IF(K.EQ.3)THEN
      PX=AMAX1(AAL1,AALO1,AALO2,AALO3,AALO4)
      IF(PX.EQ.AAL1)THEN
      R1=AHY1
      P1=AAL1
      P2=AOH1
      NR1=3
      NP2=0
      SP=SHALO 
      ELSEIF(PX.EQ.AALO1)THEN
      R1=AHY1
      P1=AALO1
      P2=AOH1
      NR1=2
      NP2=0
      SP=SHAL1 
      ELSEIF(PX.EQ.AALO2)THEN
      R1=AHY1
      P1=AALO2
      P2=AOH1
      NR1=1
      NP2=0
      SP=SHAL2
      ELSEIF(PX.EQ.AALO3)THEN
      R1=AHY1
      P1=AALO3
      P2=AOH1
      NR1=0
      NP2=0
      SP=SPAL3
      ELSEIF(PX.EQ.AALO4)THEN
      R1=AOH1
      P1=AALO4
      P2=AHY1
      NR1=0
      NP2=1
      SP=SHAL4 
      ENDIF
      RHAL1=0.0
      RHALO1=0.0
      RHALO2=0.0
      RHALO3=0.0
      RHALO4=0.0
      R1=AMAX1(ZERO,R1)
      P1=AMAX1(ZERO,P1)
      P2=AMAX1(ZERO,P2)
      SPX=SP*R1**NR1/P2**NP2
      RPALOX=AMAX1(-PALOH1,TPD*(P1-SPX))
      IF(PX.EQ.AAL1)THEN
      RHAL1=RPALOX
      ELSEIF(PX.EQ.AALO1)THEN
      RHALO1=RPALOX
      ELSEIF(PX.EQ.AALO2)THEN
      RHALO2=RPALOX
      ELSEIF(PX.EQ.AALO3)THEN
      RHALO3=RPALOX
      ELSEIF(PX.EQ.AALO4)THEN
      RHALO4=RPALOX
      ENDIF
C     IF((M/25)*25.EQ.M)THEN 
C     WRITE(*,1112)'ALOHI',I,NX,NY,L,K,M,PALOH1,AAL1
C    2,AALO1,AALO2,AALO3,AALO4,AOH1,CAL1,CALO1,CALO2,CALO3,CALO4
C    2,COH1,R1,P1,P2,SP,SPX,RPALOX,RHAL1,RHALO1,RHALO2,RHALO3,RHALO4
C    3,AAL1*AOH1**3,SPALO,A1
C     ENDIF
C
C     IRON HYDROXIDE
C
      PX=AMAX1(AFE1,AFEO1,AFEO2,AFEO3,AFEO4)
      IF(PX.EQ.AFE1)THEN
      R1=AHY1
      P1=AFE1
      P2=AOH1
      NR1=3
      NP2=0
      SP=SHFEO 
      ELSEIF(PX.EQ.AFEO1)THEN
      R1=AHY1
      P1=AFEO1
      P2=AOH1
      NR1=2
      NP2=0
      SP=SHFE1 
      ELSEIF(PX.EQ.AFEO2)THEN
      R1=AHY1
      P1=AFEO2
      P2=AOH1
      NR1=1
      NP2=0
      SP=SHFE2
      ELSEIF(PX.EQ.AFEO3)THEN
      R1=AHY1
      P1=AFEO3
      P2=AOH1
      NR1=0
      NP2=0
      SP=SPFE3
      ELSEIF(PX.EQ.AFEO4)THEN
      R1=AOH1
      P1=AFEO4
      P2=AHY1
      NR1=0
      NP2=1
      SP=SHFE4 
      ENDIF
      RHFE1=0.0
      RHFEO1=0.0
      RHFEO2=0.0
      RHFEO3=0.0
      RHFEO4=0.0
      R1=AMAX1(ZERO,R1)
      P1=AMAX1(ZERO,P1)
      P2=AMAX1(ZERO,P2)
      SPX=SP*R1**NR1/P2**NP2
      RPFEOX=AMAX1(-PFEOH1,TPD*(P1-SPX))
      IF(PX.EQ.AFE1)THEN
      RHFE1=RPFEOX
      ELSEIF(PX.EQ.AFEO1)THEN
      RHFEO1=RPFEOX
      ELSEIF(PX.EQ.AFEO2)THEN
      RHFEO2=RPFEOX
      ELSEIF(PX.EQ.AFEO3)THEN
      RHFEO3=RPFEOX
      ELSEIF(PX.EQ.AFEO4)THEN
      RHFEO4=RPFEOX
      ENDIF
C     IF((M/25)*25.EQ.M)THEN 
C     WRITE(*,1112)'FEOHI',I,NX,NY,L,K,M,PFEOH1,AFE1
C    2,AFEO1,AFEO2,AFEO3,AFEO4
C    2,AOH1,R1,P1,P2,SP,SPX,RPFEOX,RHFE1,RHFEO1,RHFEO2,RHFEO3,RHFEO4
C    3,AFE1*AOH1**3,SPFEO 
C     ENDIF
C
C     CALCITE
C
      PX=AMAX1(ACO31,AHCO31,ACO21)
      R1=AHY1
      P1=ACA1
      IF(PX.EQ.ACO31)THEN
      P2=ACO31
      NR1=0
      SP=SPCAC 
      ELSEIF(PX.EQ.AHCO31)THEN
      P2=AHCO31
      NR1=1
      SP=SHCAC1 
      ELSEIF(PX.EQ.ACO21)THEN
      P2=ACO21
      NR1=2
      SP=SHCAC2 
      ENDIF
      RHCAC3=0.0
      RHCACH=0.0
      RHCACO=0.0
      R1=AMAX1(ZERO,R1)
      P1=AMAX1(ZERO,P1)
      P2=AMAX1(ZERO,P2)
      SPX=SP*R1**NR1 
      S0=P1+P2
      S1=AMAX1(0.0,S0**2-4.0*(P1*P2-SPX))
      RPCACX=AMAX1(-PCACO1,TPD*(S0-SQRT(S1)))
      IF(PX.EQ.ACO31)THEN
      RHCAC3=RPCACX
      ELSEIF(PX.EQ.AHCO31)THEN
      RHCACH=RPCACX
      ELSEIF(PX.EQ.ACO21)THEN
      RHCACO=RPCACX
      ENDIF
C
C     GYPSUM
C
      P1=ACA1
      P2=ASO41
      P1=AMAX1(ZERO,P1)
      P2=AMAX1(ZERO,P2)
      SPX=SPCAS 
      S0=P1+P2
      S1=AMAX1(0.0,S0**2-4.0*(P1*P2-SPX))
      RPCASO=AMAX1(-PCASO1,TPD*(S0-SQRT(S1)))
C     IF((M/10)*10.EQ.M)THEN
C     WRITE(*,1112)'CALC',I,J,NX,NY,,L,K,M,PCASO1
C    2,ACO31,AHCO31,ACO21,CHY1
C    2,COH1,R1,P1,P2,P3,SP,Z,TX,RPCACX,RHCAC3,RHCACH,RHCACO
C    3,CCA1*A2*CCO3*A2,SPCAC
C     ENDIF
C
C     PHOSPHORUS PRECIPITATION-DISSOLUTION IN NON-BAND SOIL ZONE
C
C
C     ALUMINUM PHOSPHATE (VARISCITE)
C
      PX=AMAX1(AAL1,AALO1,AALO2,AALO3,AALO4)
      PY=AMAX1(AH1P1,AH2P1)
      R1=AHY1
      P3=AHY1
      IF(PY.EQ.AH1P1)THEN
      P2=AH1P1
      IF(PX.EQ.AAL1)THEN
      P1=AAL1
      NR1=1
      NP3=0
      SP=SHA0P1 
      ELSEIF(PX.EQ.AALO1)THEN
      P1=AALO1
      NR1=0
      NP3=0
      SP=SPA1P1 
      ELSEIF(PX.EQ.AALO2)THEN
      P1=AALO2
      NR1=0
      NP3=1
      SP=SHA2P1 
      ELSEIF(PX.EQ.AALO3)THEN
      P1=AALO3
      NR1=0
      NP3=2
      SP=SHA3P1 
      ELSEIF(PX.EQ.AALO4)THEN
      P1=AALO4
      NR1=0
      NP3=3
      SP=SHA4P1 
      ENDIF
      ELSE
      P2=AH2P1
      IF(PX.EQ.AAL1)THEN
      P1=AAL1
      NR1=2
      NP3=0
      SP=SHA0P2 
      ELSEIF(PX.EQ.AALO1)THEN
      P1=AALO1
      NR1=1
      NP3=0
      SP=SHA1P2
      ELSEIF(PX.EQ.AALO2)THEN
      P1=AALO2
      NR1=0
      NP3=0
      SP=SPA2P2 
      ELSEIF(PX.EQ.AALO3)THEN
      P1=AALO3
      NR1=0
      NP3=1
      SP=SHA3P2 
      ELSEIF(PX.EQ.AALO4)THEN
      P1=AALO4
      NR1=0
      NP3=2
      SP=SHA4P2 
      ENDIF
      ENDIF
      RHA0P1=0.0
      RHA1P1=0.0
      RHA2P1=0.0
      RHA3P1=0.0
      RHA4P1=0.0
      RHA0P2=0.0
      RHA1P2=0.0
      RHA2P2=0.0
      RHA3P2=0.0
      RHA4P2=0.0
      R1=AMAX1(ZERO,R1)
      P1=AMAX1(ZERO,P1)
      P2=AMAX1(ZERO,P2)
      P3=AMAX1(ZERO,P3)
      SPX=SP*R1**NR1/P3**NP3
      S0=P1+P2
      S1=AMAX1(0.0,S0**2-4.0*(P1*P2-SPX))
      RPALPX=AMAX1(-PALPO1,TPD*(S0-SQRT(S1)))
      IF(PY.EQ.AH1P1)THEN
      IF(PX.EQ.AAL1)THEN
      RHA0P1=RPALPX
      ELSEIF(PX.EQ.AALO1)THEN
      RHA1P1=RPALPX
      ELSEIF(PX.EQ.AALO2)THEN
      RHA2P1=RPALPX
      ELSEIF(PX.EQ.AALO3)THEN
      RHA3P1=RPALPX
      ELSEIF(PX.EQ.AALO4)THEN
      RHA4P1=RPALPX
      ENDIF
      ELSE
      IF(PX.EQ.AAL1)THEN
      RHA0P2=RPALPX
      ELSEIF(PX.EQ.AALO1)THEN
      RHA1P2=RPALPX
      ELSEIF(PX.EQ.AALO2)THEN
      RHA2P2=RPALPX
      ELSEIF(PX.EQ.AALO3)THEN
      RHA3P2=RPALPX
      ELSEIF(PX.EQ.AALO4)THEN
      RHA4P2=RPALPX
      ENDIF
      ENDIF
C     IF((M/25)*25.EQ.M)THEN 
C     WRITE(*,1112)'ALPO4I',I,NX,NY,L,K,M,PALPO1,AAL1
C    2,AALO1,AALO2,AALO3,AALO4
C    2,AH0P1,AH1P1,AH2P1,AHY1,AOH1,RPALPX,RHA0P1,RHA1P1,RHA2P1,RHA3P1 
C    3,RHA4P1,RHA0P2,RHA1P2,RHA2P2,RHA3P2,RHA4P2,SP,SPX,AAL1*AH0P1 
C    4,SPALP 
1112  FORMAT(A8,6I5,80E12.4)
C     ENDIF
C
C     IRON PHOSPHATE (STRENGITE)
C
      PX=AMAX1(AFE1,AFEO1,AFEO2,AFEO3,AFEO4)
      PY=AMAX1(AH1P1,AH2P1)
      R1=AHY1
      P3=AHY1
      IF(PY.EQ.AH1P1)THEN
      P2=AH1P1
      IF(PX.EQ.AFE1)THEN
      P1=AFE1
      NR1=1
      NP3=0
      SP=SHF0P1 
      ELSEIF(PX.EQ.AFEO1)THEN
      P1=AFEO1
      NR1=0
      NP3=0
      SP=SPF1P1 
      ELSEIF(PX.EQ.AFEO2)THEN
      P1=AFEO2
      NR1=0
      NP3=1
      SP=SHF2P1 
      ELSEIF(PX.EQ.AFEO3)THEN
      P1=AFEO3
      NR1=0
      NP3=2
      SP=SHF3P1 
      ELSEIF(PX.EQ.AFEO4)THEN
      P1=AFEO4
      NR1=0
      NP3=3
      SP=SHF4P1 
      ENDIF
      ELSE
      P2=AH2P1
      IF(PX.EQ.AFE1)THEN
      P1=AFE1
      NR1=2
      NP3=0
      SP=SHF0P2 
      ELSEIF(PX.EQ.AFEO1)THEN
      P1=AFEO1
      NR1=1
      NP3=0
      SP=SHF1P2 
      ELSEIF(PX.EQ.AFEO2)THEN
      P1=AFEO2
      NR1=0
      NP3=0
      SP=SPF2P2 
      ELSEIF(PX.EQ.AFEO3)THEN
      P1=AFEO3
      NR1=0
      NP3=1
      SP=SHF3P2 
      ELSEIF(PX.EQ.AFEO4)THEN
      P1=AFEO4
      NR1=0
      NP3=2
      SP=SHF4P2 
      ENDIF
      ENDIF
      RHF0P1=0.0
      RHF1P1=0.0
      RHF2P1=0.0
      RHF3P1=0.0
      RHF4P1=0.0
      RHF0P2=0.0
      RHF1P2=0.0
      RHF2P2=0.0
      RHF3P2=0.0
      RHF4P2=0.0
      R1=AMAX1(ZERO,R1)
      P1=AMAX1(ZERO,P1)
      P2=AMAX1(ZERO,P2)
      P3=AMAX1(ZERO,P3)
      SPX=SP*R1**NR1/P3**NP3
      S0=P1+P2
      S1=AMAX1(0.0,S0**2-4.0*(P1*P2-SPX))
      RPFEPX=AMAX1(-PFEPO1,TPD*(S0-SQRT(S1)))
      IF(PY.EQ.AH1P1)THEN
      IF(PX.EQ.AFE1)THEN
      RHF0P1=RPFEPX
      ELSEIF(PX.EQ.AFEO1)THEN
      RHF1P1=RPFEPX
      ELSEIF(PX.EQ.AFEO2)THEN
      RHF2P1=RPFEPX
      ELSEIF(PX.EQ.AFEO3)THEN
      RHF3P1=RPFEPX
      ELSEIF(PX.EQ.AFEO4)THEN
      RHF4P1=RPFEPX
      ENDIF
      ELSE
      IF(PX.EQ.AFE1)THEN
      RHF0P2=RPFEPX
      ELSEIF(PX.EQ.AFEO1)THEN
      RHF1P2=RPFEPX
      ELSEIF(PX.EQ.AFEO2)THEN
      RHF2P2=RPFEPX
      ELSEIF(PX.EQ.AFEO3)THEN
      RHF3P2=RPFEPX
      ELSEIF(PX.EQ.AFEO4)THEN
      RHF4P2=RPFEPX
      ENDIF
      ENDIF
C     IF((M/25)*25.EQ.M)THEN 
C     WRITE(*,1112)'FEPO4I',I,NX,NY,L,K,M,PFEPO1,AFE1
C    2,AFEO1,AFEO2,AFEO3,AFEO4
C    2,AH0P1,AH1P1,AH2P1,AHY1,AOH1,RPFEPX,RHF0P1,RHF1P1,RHF2P1,RHF3P1 
C    3,RHF4P1,RHF0P2,RHF1P2,RHF2P2,RHF3P2,RHF4P2,SP,SPX,AFE1*AH0P1 
C    4,SPFEP
C     ENDIF 
C
C     DICALCIUM PHOSPHATE
C
      PX=AMAX1(AH1P1,AH2P1)
      R1=AHY1
      P1=ACA1
      IF(PX.EQ.AH1P1)THEN
      P2=AH1P1
      NR1=0
      SP=SPCAD 
      ELSEIF(PX.EQ.AH2P1)THEN
      P2=AH2P1
      NR1=1
      SP=SHCAD2 
      ENDIF
      RPCAD1=0.0
      RHCAD2=0.0
      R1=AMAX1(ZERO,R1)
      P1=AMAX1(ZERO,P1)
      P2=AMAX1(ZERO,P2)
      SPX=SP*R1**NR1 
      S0=P1+P2
      S1=AMAX1(0.0,S0**2-4.0*(P1*P2-SPX))
      RPCADX=AMAX1(-PCAPD1,TPD*(S0-SQRT(S1)))
      IF(PX.EQ.AH1P1)THEN
      RPCAD1=RPCADX
      ELSEIF(PX.EQ.AH2P1)THEN
      RHCAD2=RPCADX
      ENDIF
C     IF((M/10)*10.EQ.M)THEN
C     WRITE(*,1112)'CAPO4',I,J,NX,NY,L,M,PCAPM1,PCAPD1,CCA1
C    2,CH1P1,CH2P1,CHY1,COH1,RPCADX,RPCAD1,RHCAD2,R1,P1,P2,P3
C    3,SP,Z,FX,Y,X,TX,A2,CCA1*A2*CH1P1*A2,SPCAD
C     ENDIF
C
C     HYDROXYAPATITE
C
      PX=AMAX1(AH1P1,AH2P1)
      R1=AHY1
      P1=ACA1
      IF(PX.EQ.AH1P1)THEN
      P2=AH1P1
      NR1=4
      SP=SHCAH1 
      ELSEIF(PX.EQ.AH2P1)THEN
      P2=AH2P1
      NR1=7
      SP=SHCAH2 
      ENDIF
      RHCAH1=0.0
      RHCAH2=0.0
      R1=AMAX1(ZERO,R1)
      P1=AMAX1(ZERO,P1)
      P2=AMAX1(ZERO,P2)
      SPX=(SP*R1**NR1/P1**5)**0.333
      RPCAHX=AMAX1(-PCAPH1,TPD*(P2-SPX))
      IF(PX.EQ.AH1P1)THEN
      RHCAH1=RPCAHX
      ELSEIF(PX.EQ.AH2P1)THEN
      RHCAH2=RPCAHX
      ENDIF
      PALOH1=PALOH1+RPALOX
      PFEOH1=PFEOH1+RPFEOX
      PCACO1=PCACO1+RPCACX
      PCASO1=PCASO1+RPCASO
      PALPO1=PALPO1+RPALPX
      PFEPO1=PFEPO1+RPFEPX
      PCAPD1=PCAPD1+RPCADX
      PCAPH1=PCAPH1+RPCAHX
C
C     ANION EXCHANGE EQILIBRIA
C
      IF(VOLW(L,NY,NX).GT.ZEROS(NY,NX))THEN
      VOLWBK=AMIN1(1.0,BKVLX/VOLW(L,NY,NX))
      ELSE
      VOLWBK=1.0
      ENDIF
      IF(XAEC(L,NY,NX).GT.ZEROS(NY,NX))THEN
      RXOH2=TAD*(XOH11*AHY1-SXOH2*XOH21)/(XOH11+SPOH2)*VOLWBK
      RXOH1=TAD*(XOH01*AHY1-SXOH1*XOH11)/(XOH01+SPOH1)*VOLWBK
      SPH2P=SXH2P*DPH2O 
      RXH2P=TAD*(XOH21*AH2P1-SPH2P*XH2P1)/(XOH21+SPH2P)*VOLWBK
      RYH2P=TAD*(XOH11*AH2P1-SXH2P*XH2P1*AOH1)
     2/(XOH11+SXH2P*AOH1)*VOLWBK 
C
C     HPO4 EXCHANGE
C
      SPH1P=SXH1P*DPH2O/DPH2P
      RXH1P=TAD*(XOH11*AH2P1-SPH1P*XH1P1)/(XOH11+SPH1P)*VOLWBK 
      XOH01=XOH01-RXOH1
      XOH11=XOH11+RXOH1-RXOH2-RYH2P-RXH1P 
      XOH21=XOH21+RXOH2-RXH2P
      XH1P1=XH1P1+RXH1P 
      XH2P1=XH2P1+RXH2P+RYH2P 
      ELSE
      RXOH2=0.0
      RXOH1=0.0
      RXH2P=0.0
      RYH2P=0.0
      RXH1P=0.0
      ENDIF
C
C     CATION EXCHANGE
C
      IF(XCEC(L,NY,NX).GT.ZEROS(NY,NX))THEN
      AALX=AAL1**0.333
      AFEX=AFE1**0.333
      ACAX=ACA1**0.500
      AMGX=AMG1**0.500
      XCAX=CCEC/(1.0+GKC4(L,NY,NX)*AN41/ACAX
     3+GKCH(L,NY,NX)*AHY1/ACAX+GKCA(L,NY,NX)*AALX/ACAX
     4+GKCA(L,NY,NX)*AFEX/ACAX+GKCM(L,NY,NX)*AMGX/ACAX
     5+GKCN(L,NY,NX)*ANA1/ACAX+GKCK(L,NY,NX)*AKA1/ACAX)
      XN4Q=XCAX*AN41*GKC4(L,NY,NX)
      XHYQ=XCAX*AHY1*GKCH(L,NY,NX) 
      XALQ=XCAX*AALX*GKCA(L,NY,NX) 
      XFEQ=XCAX*AFEX*GKCA(L,NY,NX) 
      XCAQ=XCAX*ACAX
      XMGQ=XCAX*AMGX*GKCM(L,NY,NX) 
      XNAQ=XCAX*ANA1*GKCN(L,NY,NX) 
      XKAQ=XCAX*AKA1*GKCK(L,NY,NX) 
      XTLQ=XN4Q+XHYQ+XALQ+XFEQ+XCAQ+XMGQ+XNAQ+XKAQ
      IF(XTLQ.GT.ZERO)THEN
      FX=CCEC/XTLQ
      ELSE
      FX=0.0
      ENDIF
      XN4Q=FX*XN4Q
      XHYQ=FX*XHYQ
      XALQ=FX*XALQ/3.0
      XFEQ=FX*XFEQ/3.0
      XCAQ=FX*XCAQ/2.0
      XMGQ=FX*XMGQ/2.0
      XNAQ=FX*XNAQ
      XKAQ=FX*XKAQ
      RXN4=TAD*AMIN1((XN4Q-XN41)*AN41/XN4Q,CN41)
      RXHY=TAD*AMIN1((XHYQ-XHY1)*AHY1/XHYQ,CHY1)
      RXAL=TAD*AMIN1((XALQ-XAL1)*AALX/XALQ,CAL1)
      RXFE=TAD*AMIN1((XFEQ-XFE1)*AFEX/XFEQ,CFE1)
      RXCA=TAD*AMIN1((XCAQ-XCA1)*ACAX/XCAQ,CCA1)
      RXMG=TAD*AMIN1((XMGQ-XMG1)*AMGX/XMGQ,CMG1)
      RXNA=TAD*AMIN1((XNAQ-XNA1)*ANA1/XNAQ,CNA1)
      RXKA=TAD*AMIN1((XKAQ-XKA1)*AKA1/XKAQ,CKA1)
      XN41=XN41+RXN4
      XHY1=XHY1+RXHY 
      XAL1=XAL1+RXAL 
      XFE1=XFE1+RXFE 
      XCA1=XCA1+RXCA 
      XMG1=XMG1+RXMG 
      XNA1=XNA1+RXNA 
      XKA1=XKA1+RXKA
      ELSE
      RXN4=0.0
      RXHY=0.0
      RXAL=0.0
      RXFE=0.0
      RXCA=0.0
      RXMG=0.0
      RXNA=0.0
      RXKA=0.0
      ENDIF
C     IF((M/25)*25.EQ.M)THEN 
C     WRITE(*,1112)'RXALI',I,NX,NY,L,K,M
C    2,RXAL,TAD,XALQ,XAL1,AALX,CCEC,XCAX,FX
C    3,GKCA(L,NY,NX),XCEC(L,NY,NX),BKVLX,VOLA(L,NY,NX)
C     ENDIF
C
C     ORGANIC MATTER
C
      DP=DPCOH*DPALO
      XHC1=AHY1*(XCOOH-XALO21-XFEO21)/(AHY1+DPCOH)
      XALO21=AALO2*(XCOOH-XHC1)/(AALO2+DPALO)
      XFEO21=AFEO2*(XCOOH-XHC1)/(AFEO2+DPFEO)
      XCOO=AMAX1(0.0,XCOOH-XHC1-XALO21-XFEO21)
      ELSE
      RHAL1=0.0
      RHALO1=0.0
      RHALO2=0.0
      RHALO3=0.0
      RHALO4=0.0
      RHFE1=0.0
      RHFEO1=0.0
      RHFEO2=0.0
      RHFEO3=0.0
      RHFEO4=0.0
      RHCAC3=0.0
      RHCACH=0.0
      RHCACO=0.0
      RPCACX=0.0
      RPCASO=0.0
      RPCADX=0.0
      RPCAHX=0.0
      RHA0P1=0.0
      RHA1P1=0.0
      RHA2P1=0.0
      RHA3P1=0.0
      RHA4P1=0.0
      RHA0P2=0.0
      RHA1P2=0.0
      RHA2P2=0.0
      RHA3P2=0.0
      RHA4P2=0.0
      RHF0P1=0.0
      RHF1P1=0.0
      RHF2P1=0.0
      RHF3P1=0.0
      RHF4P1=0.0
      RHF0P2=0.0
      RHF1P2=0.0
      RHF2P2=0.0
      RHF3P2=0.0
      RHF4P2=0.0
      RPCAD1=0.0
      RHCAD2=0.0
      RHCAH1=0.0
      RHCAH2=0.0
      RXOH2=0.0
      RXOH1=0.0
      RXH2P=0.0
      RYH2P=0.0
      RXH1P=0.0
      RXN4=0.0
      RXHY=0.0
      RXAL=0.0
      RXFE=0.0
      RXCA=0.0
      RXMG=0.0 
      RXNA=0.0 
      RXKA=0.0
      ENDIF
C
C     ION SPECIATION
C
      S0=AHY1+AN31+DPN4
      S1=S0**2-4.0*(AHY1*AN31-DPN4*AN41)
      RNH4=TSL*(S0-SQRT(S1))
      S0=AAL1+AOH1+DPAL1
      S1=S0**2-4.0*(AAL1*AOH1-DPAL1*AALO1)
      RALO1=TSL*(S0-SQRT(S1))
      S0=AALO1+AOH1+DPAL2
      S1=S0**2-4.0*(AALO1*AOH1-DPAL2*AALO2)
      RALO2=TSL*(S0-SQRT(S1))
      S0=CALO2+COH1+DPAL3
      S1=S0**2-4.0*(AALO2*AOH1-DPAL3*AALO3)
      RALO3=TSL*(S0-SQRT(S1))
      S0=AALO3+AOH1+DPAL4
      S1=S0**2-4.0*(AALO3*AOH1-DPAL4*AALO4)
      RALO4=TSL*(S0-SQRT(S1))
      S0=AAL1+ASO41+DPALS
      S1=S0**2-4.0*(AAL1*ASO41-DPALS*AALS1)
      RALS=TSL*(S0-SQRT(S1))
      S0=AFE1+AOH1+DPFE1
      S1=S0**2-4.0*(AFE1*AOH1-DPFE1*AFEO1)
      RFEO1=TSL*(S0-SQRT(S1))
      S0=AFEO1+AOH1+DPFE2
      S1=S0**2-4.0*(AFEO1*AOH1-DPFE2*AFEO2)
      RFEO2=TSL*(S0-SQRT(S1))
      S0=AFEO2+AOH1+DPFE3
      S1=S0**2-4.0*(AFEO2*AOH1-DPFE3*AFEO3)
      RFEO3=TSL*(S0-SQRT(S1))
      S0=AFEO3+AOH1+DPFE4
      S1=S0**2-4.0*(AFEO3*AOH1-DPFE4*AFEO4)
      RFEO4=TSL*(S0-SQRT(S1))
      S0=AFE1+ASO41+DPFES
      S1=S0**2-4.0*(AFE1*ASO41-DPFES*AFES1)
      RFES=TSL*(S0-SQRT(S1))
      S0=ACA1+AOH1+DPCAO
      S1=S0**2-4.0*(ACA1*AOH1-DPCAO*ACAO1)
      RCAO=TSL*(S0-SQRT(S1))
      S0=ACA1+ACO31+DPCAC
      S1=S0**2-4.0*(ACA1*ACO31-DPCAC*ACAC1)
      RCAC=TSL*(S0-SQRT(S1))
      S0=ACA1+AHCO31+DPCAH
      S1=S0**2-4.0*(ACA1*AHCO31-DPCAH*ACAH1)
      RCAH=TSL*(S0-SQRT(S1))
      S0=ACA1+ASO41+DPCAS
      S1=S0**2-4.0*(ACA1*ASO41-DPCAS*ACAS1)
      RCAS=TSL*(S0-SQRT(S1))
      S0=AMG1+AOH1+DPMGO
      S1=S0**2-4.0*(AMG1*AOH1-DPMGO*AMGO1)
      RMGO=TSL*(S0-SQRT(S1))
      S0=AMG1+ACO31+DPMGC
      S1=S0**2-4.0*(AMG1*ACO31-DPMGC*AMGC1)
      RMGC=TSL*(S0-SQRT(S1))
      S0=AMG1+AHCO31+DPMGH
      S1=S0**2-4.0*(AMG1*AHCO31-DPMGH*AMGH1)
      RMGH=TSL*(S0-SQRT(S1))
      S0=AMG1+ASO41+DPMGS
      S1=S0**2-4.0*(AMG1*ASO41-DPMGS*AMGS1)
      RMGS=TSL*(S0-SQRT(S1))
      S0=ANA1+ACO31+DPNAC
      S1=S0**2-4.0*(ANA1*ACO31-DPNAC*ANAC1)
      RNAC=TSL*(S0-SQRT(S1))
      S0=ANA1+ASO41+DPNAS
      S1=S0**2-4.0*(ANA1*ASO41-DPNAS*ANAS1)
      RNAS=TSL*(S0-SQRT(S1))
      S0=AKA1+ASO41+DPKAS
      S1=S0**2-4.0*(AKA1*ASO41-DPKAS*AKAS1)
      RKAS=TSL*(S0-SQRT(S1))
      S0=AH0P1+AHY1+DPH1P
      S1=S0**2-4.0*(AH0P1*AHY1-DPH1P*AH1P1)
      RH1P=TSL*(S0-SQRT(S1))
      S0=AH1P1+AHY1+DPH2P
      S1=S0**2-4.0*(AH1P1*AHY1-DPH2P*AH2P1)
      RH2P=TSL*(S0-SQRT(S1))
      S0=AH2P1+AHY1+DPH3P
      S1=S0**2-4.0*(AH2P1*AHY1-DPH3P*AH3P1)
      RH3P=TSL*(S0-SQRT(S1))
      S0=AFE1+AH1P1+DPF1P
      S1=S0**2-4.0*(AFE1*AH1P1-DPF1P*AF1P1)
      RF1P=TSL*(S0-SQRT(S1))
      S0=AFE1+AH2P1+DPF2P
      S1=S0**2-4.0*(AFE1*AH2P1-DPF2P*AF2P1)
      RF2P=TSL*(S0-SQRT(S1))
      S0=ACA1+AH0P1+DPC0P
      S1=S0**2-4.0*(ACA1*AH0P1-DPC0P*AC0P1)
      RC0P=TSL*(S0-SQRT(S1))
      S0=ACA1+AH1P1+DPC1P
      S1=S0**2-4.0*(ACA1*AH1P1-DPC1P*AC1P1)
      RC1P=TSL*(S0-SQRT(S1))
      S0=ACA1+AH2P1+DPC2P
      S1=S0**2-4.0*(ACA1*AH2P1-DPC2P*AC2P1)
      RC2P=TSL*(S0-SQRT(S1))
      S0=AMG1+AH1P1+DPM1P
      S1=S0**2-4.0*(AMG1*AH1P1-DPM1P*AM1P1)
      RM1P=TSL*(S0-SQRT(S1))
C
C     TOTAL ION FLUX FOR EACH ION SPECIES
C
      RN4S=RNH4-RXN4
      RN3S=-RNH4
      RAL=-RHAL1-RHA0P1-RHA0P2-RALO1-RALS-RXAL
      RFE=-RHFE1-RHF0P1-RHF0P2-RFEO1-RFES-RXFE
     2-RF1P-RF2P
      RCA=-RPCACX-RPCASO-RPCADX-5.0*RPCAHX-RXCA
     2-RCAO-RCAC-RCAH-RCAS-RC0P-RC1P-RC2P
      RMG=-RMGO-RMGC-RMGH-RMGS-RM1P-RXMG
      RNA=-RNAC-RNAS-RXNA
      RKA=-RKAS-RXKA
      RSO4=-RPCASO-RALS-RFES-RCAS-RMGS-RNAS-RKAS
      RAL1=-RHALO1-RHA1P1-RHA1P2+RALO1-RALO2
      RAL2=-RHALO2-RHA2P1-RHA2P2+RALO2-RALO3
      RAL3=-RHALO3-RHA3P1-RHA3P2+RALO3-RALO4
      RAL4=-RHALO4-RHA4P1-RHA4P2+RALO4
      RFE1=-RHFEO1-RHF1P1-RHF1P2+RFEO1-RFEO2
      RFE2=-RHFEO2-RHF2P1-RHF2P2+RFEO2-RFEO3
      RFE3=-RHFEO3-RHF3P1-RHF3P2+RFEO3-RFEO4
      RFE4=-RHFEO4-RHF4P1-RHF4P2+RFEO4
      RHP0=-RH1P-RC0P
      RHP1=-RHA0P1-RHA1P1-RHA2P1-RHA3P1
     2-RHA4P1-RHF0P1-RHF1P1-RHF2P1 
     3-RHF3P1-RHF4P1-RPCAD1-3.0*RHCAH1 
     4-RXH1P+RH1P-RH2P-RF1P-RC1P-RM1P
      RHP2=-RHA0P2-RHA1P2-RHA2P2-RHA3P2
     2-RHA4P2-RHF0P2-RHF1P2-RHF2P2 
     3-RHF3P2-RHF4P2-RHCAD2-3.0*RHCAH2
     4-RXH2P-RYH2P+RH2P-RH3P-RF2P-RC2P
      RHP3=RH3P
C
C     ION CONCENTRATIONS
C
      CCO2X=CCO2M*SCO2X/(EXP(ACO2X*CSTRZ))
     2*EXP(0.843-0.0281*ATCA(NY,NX))*FH2O
      CCO2Y=LOG(CCO2X)
      CCO2Z=ABS(CCO2Y)
      CCO21=CCO2X
      FCO3=DPCO3*A0/(AHY1**2*A2)
      FHCO=DPCO2*A0/(AHY1*A1)
      Z=ACO2X*(2.0E-03*FCO3+0.5E-03*FHCO)
      DO 2120 MM=1,25
      R=(LOG(CCO21)+Z*CCO21-CCO2Y)/CCO2Z
      IF(R.LT.1.0E-03)GO TO 2130
      CCO21=CCO21/SQRT(1.0+R)
2120  CONTINUE
2130  CONTINUE
      CCH41=CCH4M*SCH4X/(EXP(ACH4X*CSTR1))
     2*EXP(0.597-0.0199*ATCA(NY,NX))*FH2O
      COXY1=COXYM*SOXYX/(EXP(AOXYX*CSTR1))
     2*EXP(0.516-0.0172*ATCA(NY,NX))*FH2O
      CZ2G1=CZ2GM*SN2GX/(EXP(AN2GX*CSTR1))
     2*EXP(0.456-0.0152*ATCA(NY,NX))*FH2O
      CZ2O1=CZ2OM*SN2OX/(EXP(AN2OX*CSTR1))
     2*EXP(0.897-0.0299*ATCA(NY,NX))*FH2O
      CN41=CN41+RN4S 
      CN31=CN31+RN3S
      CAL1=CAL1+RAL
      CFE1=CFE1+RFE
      CCA1=CCA1+RCA
      CMG1=CMG1+RMG
      CNA1=CNA1+RNA
      CKA1=CKA1+RKA
      CSO41=CSO41+RSO4
      CCO31=CCO21*DPCO3*A0/(AHY1**2*A2)
      CHCO31=CCO21*DPCO2*A0/(AHY1*A1)
      CALO1=CALO1+RAL1
      CALO2=CALO2+RAL2
      CALO3=CALO3+RAL3
      CALO4=CALO4+RAL4
      CALS1=CALS1+RALS
      CFEO1=CFEO1+RFE1
      CFEO2=CFEO2+RFE2
      CFEO3=CFEO3+RFE3
      CFEO4=CFEO4+RFE4
      CFES1=CFES1+RFES
      CCAO1=CCAO1+RCAO
      CCAC1=CCAC1+RCAC
      CCAH1=CCAH1+RCAH
      CCAS1=CCAS1+RCAS
      CMGO1=CMGO1+RMGO
      CMGC1=CMGC1+RMGC
      CMGH1=CMGH1+RMGH
      CMGS1=CMGS1+RMGS
      CNAC1=CNAC1+RNAC
      CNAS1=CNAS1+RNAS
      CKAS1=CKAS1+RKAS
      CH0P1=CH0P1+RHP0
      CH1P1=CH1P1+RHP1
      CH2P1=CH2P1+RHP2
      CH3P1=CH3P1+RHP3
      CF1P1=CF1P1+RF1P
      CF2P1=CF2P1+RF2P
      CC0P1=CC0P1+RC0P
      CC1P1=CC1P1+RC1P
      CC2P1=CC2P1+RC2P
      CM1P1=CM1P1+RM1P
      IF(K.EQ.3.AND.(M/1)*1.EQ.M)THEN
C     WRITE(*,1112)'A1I',I,NX,NY,L,K,M,A1,A2,A3,FSTR2,CSTR1
C    2,CSTR2,CC3,CA3,CC2,CA2,CC1,CA1,VOLW(L,NY,NX) 
C     WRITE(*,1112)'ALPO4I',I,NX,NY,L,K,M,PALPO1,AAL1
C    2,AALO1,AALO2,AALO3,AALO4
C    2,AH0P1,AH1P1,AH2P1,AHY1,AOH1,RPALPX,RHA0P1,RHA1P1,RHA2P1,RHA3P1 
C    3,RHA4P1,RHA0P2,RHA1P2,RHA2P2,RHA3P2,RHA4P2,SP,SPX,AAL1*AH0P1 
C    4,SPALP,CH0P1,CH1P1,CH2P1,CH3P1,RHP0,RHP1,RHP2,RHP3
C    5,RAL,RHAL1,RHA0P1,RHA0P2,RALO1,RALS,RXAL 
C     WRITE(*,1112)'FEPO4I',I,NX,NY,L,K,M,PFEPO1,AFE1
C    2,AFEO1,AFEO2,AFEO3,AFEO4
C    2,AH0P1,AH1P1,AH2P1,AHY1,AOH1,RPFEPX,RHF0P1,RHF1P1,RHF2P1,RHF3P1 
C    3,RHF4P1,RHF0P2,RHF1P2,RHF2P2,RHF3P2,RHF4P2,SP,SPX,AFE1*AH0P1 
C    4,SPFEP
C     WRITE(*,1112)'APATITEI',I,NX,NY,L,K,M,PCAPH1,ACA1,XCA1 
C    2,AH0P1,AH1P1,AH2P1,AHY1,AOH1,RPCAHX,RHCAH1,RHCAH2 
C    3,SP,SPX,ACA1**5*AH0P1**3*AOH1,SPCAH,SHCAH1,SHCAH2
C    3,CH0P1,CH1P1,CH2P1,XOH01,XOH11,XOH21,XH1P1,XH2P1
C    4,RHA0P1,RHA1P1,RHA2P1,RHA3P1
C    2,RHA4P1,RHF0P1,RHF1P1,RHF2P1 
C    3,RHF3P1,RHF4P1,RPCAD1,3.0*RHCAH1 
C    4,RXH1P,RH1P,RH2P,RF1P,RC1P,RM1P
C    5,RHA0P2,RHA1P2,RHA2P2,RHA3P2
C    2,RHA4P2,RHF0P2,RHF1P2,RHF2P2 
C    3,RHF3P2,RHF4P2,RHCAD2,3.0*RHCAH2
C    4,RXH2P,RYH2P,RH2P,RH3P,RF2P,RC2P,RH3P 
      ENDIF
1000  CONTINUE
C
C     CONVERGE TOWARDS ALL SOLUBILITY EQUILIBRIA
C     IF SALT OPTION IS NOT SELECTED
C
      ELSE
      DO 1100 M=1,MRXN
      CCO21=AMAX1(ZERO,CCO21)
      CCO31=CCO21*DPCO3*A0/(CHY1**2*A2)
      CHCO31=CCO21*DPCO2*A0/(CHY1*A1)
      CN41=AMAX1(ZERO,CN41)
      CN31=AMAX1(ZERO,CN31)
      CAL1=AMAX1(ZERO,CAL1)
      CFE1=AMAX1(ZERO,CFE1)
      CCA1=AMAX1(ZERO,CCA1)
      CCA1=AMIN1(CCA1,SPCAC/(CCO31*A2**2))
      CMG1=AMAX1(ZERO,CMG1)
      CNA1=AMAX1(ZERO,CNA1)
      CKA1=AMAX1(ZERO,CKA1)
      CH1P1=AMAX1(ZERO,CH1P1)
      CH2P1=AMAX1(ZERO,CH2P1)
C
C     PRECIPITATION-DISSOLUTION FLUXES
C
      IF(K.EQ.3)THEN
      CH2PA=SYA0P2/(CAL1*COH1**2)
      RPALPX=AMAX1(-PALPO1,TPD*(CH2P1-CH2PA))
      CH2PF=SYF0P2/(CFE1*COH1**2)
      RPFEPX=AMAX1(-PFEPO1,TPD*(CH2P1-CH2PF))
      CH2PD=SYCAD2/(CCA1*COH1)
      RPCADX=AMAX1(-PCAPD1,TPD*(CH2P1-CH2PD))
      CH2PH=(SYCAH2/(CCA1**5*COH1**7))**0.333
      RPCAHX=AMAX1(-PCAPH1,TPD*(CH2P1-CH2PH))
      PALPO1=PALPO1+RPALPX
      PFEPO1=PFEPO1+RPFEPX
      PCAPD1=PCAPD1+RPCADX
      PCAPH1=PCAPH1+RPCAHX
C
C     ANION EXCHANGE FLUXES
C
      IF(VOLW(L,NY,NX).GT.ZEROS(NY,NX))THEN
      VOLWBK=AMIN1(1.0,BKVLX/VOLW(L,NY,NX))
      ELSE
      VOLWBK=1.0
      ENDIF
      IF(XAEC(L,NY,NX).GT.ZEROS(NY,NX))THEN
      SPH2P=SXH2P*DPH2O
      RXH2P=TAD*(XOH21*CH2P1-SPH2P*XH2P1)/(XOH21+SPH2P)*VOLWBK 
      RYH2P=TAD*(XOH11*CH2P1-SXH2P*COH1*XH2P1)
     2/(XOH11+SXH2P*COH1)*VOLWBK 
      SPH1P=SXH1P*DPH2O/DPH2P
      RXH1P=TAD*(XOH11*CH1P1-SPH1P*XH1P1)/(XOH11+SPH1P)*VOLWBK 
      XOH11=XOH11-RYH2P-RXH1P
      XOH21=XOH21-RXH2P
      XH1P1=XH1P1+RXH1P
      XH2P1=XH2P1+RXH2P+RYH2P
      ELSE
      RXH2P=0.0
      RYH2P=0.0
      RXH1P=0.0
      ENDIF
C
C     CATION EXCHANGE
C
      IF(XCEC(L,NY,NX).GT.ZEROS(NY,NX))THEN
      CALX=CAL1**0.333
      CFEX=CFE1**0.333
      CCAX=CCA1**0.500
      CMGX=CMG1**0.500
      XCAX=CCEC/(1.0+GKC4(L,NY,NX)*CN41/CCAX
     3+GKCH(L,NY,NX)*CHY1/CCAX+GKCA(L,NY,NX)*CALX/CCAX
     4+GKCA(L,NY,NX)*CFEX/CCAX+GKCM(L,NY,NX)*CMGX/CCAX
     5+GKCN(L,NY,NX)*CNA1/CCAX+GKCK(L,NY,NX)*CKA1/CCAX)
      XN4Q=XCAX*CN41*GKC4(L,NY,NX)
      XHYQ=XCAX*CHY1*GKCH(L,NY,NX)
      XALQ=XCAX*CALX*GKCA(L,NY,NX)
      XFEQ=XCAX*CFEX*GKCA(L,NY,NX)
      XCAQ=XCAX*CCAX
      XMGQ=XCAX*CMGX*GKCM(L,NY,NX)
      XNAQ=XCAX*CNA1*GKCN(L,NY,NX)
      XKAQ=XCAX*CKA1*GKCK(L,NY,NX)
      XTLQ=XN4Q+XHYQ+XALQ+XFEQ+XCAQ+XMGQ+XNAQ+XKAQ
      IF(XTLQ.GT.ZERO)THEN
      FX=CCEC/XTLQ
      ELSE
      FX=0.0
      ENDIF
      XN4Q=FX*XN4Q
      RXN4=TSL*AMIN1((XN4Q-XN41)*CN41/XN4Q,CN41)
      XN41=XN41+RXN4
      ELSE
      RXN4=0.0
      ENDIF 
      ELSE
      RPALPX=0.0
      RPFEPX=0.0
      RPCADX=0.0
      RPCAHX=0.0
      RXH2P=0.0
      RYH2P=0.0
      RXH1P=0.0
      RXN4=0.0
      ENDIF
C
C     NH4 – NH3
C
      S0=CHY1+CN31+DPN4
      S1=AMAX1(0.0,S0**2-4.0*(CHY1*CN31-DPN4*CN41))
      RNH4=TSL*(S0-SQRT(S1))
C
C     H2PO4 – HPO4
C
      S0=CH1P1+CHY1+DPH2P
      S1=AMAX1(0.0,S0**2-4.0*(CH1P1*CHY1-DPH2P*CH2P1))
      RH2P=TSL*(S0-SQRT(S1))
C
C     ION FLUXES
C
      RN4S=RNH4-RXN4
      RN3S=-RNH4
      RHP1=-RXH1P-RH2P 
      RHP2=-RXH2P-RYH2P+RH2P-RPALPX-RPFEPX-RPCADX-3.0*RPCAHX
      CN41=CN41+RN4S
      CN31=CN31+RN3S
      CH1P1=CH1P1+RHP1
      CH2P1=CH2P1+RHP2
      IF(K.EQ.3)THEN
C     WRITE(*,2222)'RNH4E',K,L,CN41,RN4S,RNH4,RXN4
C    2,TSL,XN4Q,XN41,CN41,FX,XTLQ,XHYQ,XALQ,XFEQ
C    2,XCAQ,XMGQ,XNAQ,XKAQ,CCEC,BKVLX,XCAX,CAL1,CFE1
C    2,CCA1,CMG1 
C     WRITE(*,2222)'RHP1E',K,L,CH2P1,CH1P1,RHP2,RHP1
C    2,XH2P1,XH1P1,RXH1P,RXH2P,RYH2P,RH2P,CHY1,COH1
C    3,XOH21,XOH11 
      ENDIF
1100  CONTINUE
      ENDIF
C
C     SOLUTE CONCENTRATIONS IN PRECIPITATION
C
      IF(K.EQ.1.AND.I.EQ.1.AND.L.EQ.1)THEN
      CCOR(NY,NX)=CCO21
      CCHR(NY,NX)=CCH41
      COXR(NY,NX)=COXY1
      CNNR(NY,NX)=CZ2G1
      CN2R(NY,NX)=CZ2O1
      CN4R(NY,NX)=CN41
      CN3R(NY,NX)=CN31
      CALR(NY,NX)=CAL1
      CFER(NY,NX)=CFE1
      CHYR(NY,NX)=CHY1
      CCAR(NY,NX)=CCA1
      CMGR(NY,NX)=CMG1
      CNAR(NY,NX)=CNA1
      CKAR(NY,NX)=CKA1
      COHR(NY,NX)=COH1
      CSOR(NY,NX)=CSO41
      CCLR(NY,NX)=CCL1
      CC3R(NY,NX)=CCO31
      CHCR(NY,NX)=CHCO31
      CAL1R(NY,NX)=CALO1
      CAL2R(NY,NX)=CALO2
      CAL3R(NY,NX)=CALO3
      CAL4R(NY,NX)=CALO4
      CALSR(NY,NX)=CALS1
      CFE1R(NY,NX)=CFEO1
      CFE2R(NY,NX)=CFEO2
      CFE3R(NY,NX)=CFEO3
      CFE4R(NY,NX)=CFEO4
      CFESR(NY,NX)=CFES1
      CCAOR(NY,NX)=CCAO1
      CCACR(NY,NX)=CCAC1
      CCAHR(NY,NX)=CCAH1
      CCASR(NY,NX)=CCAS1
      CMGOR(NY,NX)=CMGO1
      CMGCR(NY,NX)=CMGC1
      CMGHR(NY,NX)=CMGH1
      CMGSR(NY,NX)=CMGS1
      CNACR(NY,NX)=CNAC1
      CNASR(NY,NX)=CNAS1
      CKASR(NY,NX)=CKAS1
      CH0PR(NY,NX)=CH0P1
      CH1PR(NY,NX)=CH1P1
      CPOR(NY,NX)=CH2P1
      CH3PR(NY,NX)=CH3P1
      CF1PR(NY,NX)=CF1P1
      CF2PR(NY,NX)=CF2P1
      CC0PR(NY,NX)=CC0P1
      CC1PR(NY,NX)=CC1P1
      CC2PR(NY,NX)=CC2P1
      CM1PR(NY,NX)=CM1P1
      CSTRR(NY,NX)=CSTR1
C
C     SOLUTE CONCENTRATIONS IN IRRIGATION
C
      ELSEIF(K.EQ.2.AND.L.EQ.1)THEN
      CCOQ(NY,NX)=CCO21
      CCHQ(NY,NX)=CCH41
      COXQ(NY,NX)=COXY1
      CNNQ(NY,NX)=CZ2G1
      CN2Q(NY,NX)=CZ2O1
      CN4Q(I,NY,NX)=CN41
      CN3Q(I,NY,NX)=CN31
      CALQ(I,NY,NX)=CAL1
      CFEQ(I,NY,NX)=CFE1
      CHYQ(I,NY,NX)=CHY1
      CCAQ(I,NY,NX)=CCA1
      CMGQ(I,NY,NX)=CMG1
      CNAQ(I,NY,NX)=CNA1
      CKAQ(I,NY,NX)=CKA1
      COHQ(I,NY,NX)=COH1
      CSOQ(I,NY,NX)=CSO41
      CCLQ(I,NY,NX)=CCL1
      CC3Q(I,NY,NX)=CCO31
      CHCQ(I,NY,NX)=CHCO31
      CAL1Q(I,NY,NX)=CALO1
      CAL2Q(I,NY,NX)=CALO2
      CAL3Q(I,NY,NX)=CALO3
      CAL4Q(I,NY,NX)=CALO4
      CALSQ(I,NY,NX)=CALS1
      CFE1Q(I,NY,NX)=CFEO1
      CFE2Q(I,NY,NX)=CFEO2
      CFE3Q(I,NY,NX)=CFEO3
      CFE4Q(I,NY,NX)=CFEO4
      CFESQ(I,NY,NX)=CFES1
      CCAOQ(I,NY,NX)=CCAO1
      CCACQ(I,NY,NX)=CCAC1
      CCAHQ(I,NY,NX)=CCAH1
      CCASQ(I,NY,NX)=CCAS1
      CMGOQ(I,NY,NX)=CMGO1
      CMGCQ(I,NY,NX)=CMGC1
      CMGHQ(I,NY,NX)=CMGH1
      CMGSQ(I,NY,NX)=CMGS1
      CNACQ(I,NY,NX)=CNAC1
      CNASQ(I,NY,NX)=CNAS1
      CKASQ(I,NY,NX)=CKAS1
      CH0PQ(I,NY,NX)=CH0P1
      CH1PQ(I,NY,NX)=CH1P1
      CPOQ(I,NY,NX)=CH2P1
      CH3PQ(I,NY,NX)=CH3P1
      CF1PQ(I,NY,NX)=CF1P1
      CF2PQ(I,NY,NX)=CF2P1
      CC0PQ(I,NY,NX)=CC0P1
      CC1PQ(I,NY,NX)=CC1P1
      CC2PQ(I,NY,NX)=CC2P1
      CM1PQ(I,NY,NX)=CM1P1
      CSTRQ(I,NY,NX)=CSTR1
C
C     SOLUTE CONCENTRATIONS IN SOIL
C
      ELSEIF(K.EQ.3.AND.I.EQ.1.AND.DATA(20).EQ.'NO'.AND.IGO.EQ.0)THEN
      CCOU=CCO21
      CCHU=CCH41
      COXU=0.0
      CNNU=CZ2G1
      CN2U=CZ2O1
      CN4U(L,NY,NX)=CN41
      CN3U(L,NY,NX)=CN31
      CNOU(L,NY,NX)=CNOX
      CALU(L,NY,NX)=CAL1
      CFEU(L,NY,NX)=CFE1
      CHYU(L,NY,NX)=CHY1
      CCAU(L,NY,NX)=CCA1
      CMGU(L,NY,NX)=CMG1
      CNAU(L,NY,NX)=CNA1
      CKAU(L,NY,NX)=CKA1
      COHU(L,NY,NX)=COH1
      CSOU(L,NY,NX)=CSO41
      CCLU(L,NY,NX)=CCL1
      CC3U(L,NY,NX)=CCO31
      CHCU(L,NY,NX)=CHCO31
      CAL1U(L,NY,NX)=CALO1
      CAL2U(L,NY,NX)=CALO2
      CAL3U(L,NY,NX)=CALO3
      CAL4U(L,NY,NX)=CALO4
      CALSU(L,NY,NX)=CALS1
      CFE1U(L,NY,NX)=CFEO1
      CFE2U(L,NY,NX)=CFEO2
      CFE3U(L,NY,NX)=CFEO3
      CFE4U(L,NY,NX)=CFEO4
      CFESU(L,NY,NX)=CFES1
      CCAOU(L,NY,NX)=CCAO1
      CCACU(L,NY,NX)=CCAC1
      CCAHU(L,NY,NX)=CCAH1
      CCASU(L,NY,NX)=CCAS1
      CMGOU(L,NY,NX)=CMGO1
      CMGCU(L,NY,NX)=CMGC1
      CMGHU(L,NY,NX)=CMGH1
      CMGSU(L,NY,NX)=CMGS1
      CNACU(L,NY,NX)=CNAC1
      CNASU(L,NY,NX)=CNAS1
      CKASU(L,NY,NX)=CKAS1
      CH0PU(L,NY,NX)=CH0P1
      CH1PU(L,NY,NX)=CH1P1
      CH2PU(L,NY,NX)=CH2P1
      CH3PU(L,NY,NX)=CH3P1
      CF1PU(L,NY,NX)=CF1P1
      CF2PU(L,NY,NX)=CF2P1
      CC0PU(L,NY,NX)=CC0P1
      CC1PU(L,NY,NX)=CC1P1
      CC2PU(L,NY,NX)=CC2P1
      CM1PU(L,NY,NX)=CM1P1
C
C     INITIAL STATE VARIABLES FOR GAS IN SOIL
C
      CO2G(L,NY,NX)=CCO2EI(NY,NX)*VOLP(L,NY,NX)
      CH4G(L,NY,NX)=CCH4E(NY,NX)*VOLP(L,NY,NX)
      OXYG(L,NY,NX)=COXYE(NY,NX)*VOLP(L,NY,NX)
      Z2GG(L,NY,NX)=CZ2GE(NY,NX)*VOLP(L,NY,NX)
      Z2OG(L,NY,NX)=CZ2OE(NY,NX)*VOLP(L,NY,NX)
      ZNH3G(L,NY,NX)=CNH3E(NY,NX)*VOLP(L,NY,NX)
      H2GG(L,NY,NX)=CH2GE(NY,NX)*VOLP(L,NY,NX)
      IF(CDPTH(L-1,NY,NX).LT.DTBLZ(NY,NX))THEN
      OXYS(L,NY,NX)=COXYE(NY,NX)*SOXYX/(EXP(AOXYX*CSTR1))
     2*EXP(0.516-0.0172*ATCA(NY,NX))*FH2O*VOLW(L,NY,NX)
      ELSE
      OXYS(L,NY,NX)=0.0
      ENDIF
      CO2S(L,NY,NX)=CCO2EI(NY,NX)*SCO2X/(EXP(ACO2X*CSTR1))
     2*EXP(0.843-0.0281*ATCA(NY,NX))*FH2O*VOLW(L,NY,NX)
      CH4S(L,NY,NX)=CCH4E(NY,NX)*SCH4X/(EXP(ACH4X*CSTR1))
     2*EXP(0.597-0.0199*ATCA(NY,NX))*FH2O*VOLW(L,NY,NX)
      Z2GS(L,NY,NX)=CZ2GE(NY,NX)*SN2GX/(EXP(AN2GX*CSTR1))
     2*EXP(0.456-0.0152*ATCA(NY,NX))*FH2O*VOLW(L,NY,NX)
      Z2OS(L,NY,NX)=CZ2OE(NY,NX)*SN2OX/(EXP(AN2OX*CSTR1))
     2*EXP(0.897-0.0299*ATCA(NY,NX))*FH2O*VOLW(L,NY,NX)
      H2GS(L,NY,NX)=CH2GE(NY,NX)*SH2GX/(EXP(AH2GX*CSTR1))
     2*EXP(0.597-0.0199*ATCA(NY,NX))*FH2O*VOLW(L,NY,NX)
C
C     INITIAL STATE VARIABLES FOR MINERAL N AND P IN SOIL
C
      ZNH4S(L,NY,NX)=CN4U(L,NY,NX)*VOLW(L,NY,NX)*VLNH4(L,NY,NX)*14.0
      ZNH4B(L,NY,NX)=CN4U(L,NY,NX)*VOLW(L,NY,NX)*VLNHB(L,NY,NX)*14.0
      ZNH3S(L,NY,NX)=CN3U(L,NY,NX)*VOLW(L,NY,NX)*VLNH4(L,NY,NX)*14.0
      ZNH3B(L,NY,NX)=CN3U(L,NY,NX)*VOLW(L,NY,NX)*VLNHB(L,NY,NX)*14.0
      ZNO3S(L,NY,NX)=CNOU(L,NY,NX)*VOLW(L,NY,NX)*VLNO3(L,NY,NX)*14.0
      ZNO3B(L,NY,NX)=CNOU(L,NY,NX)*VOLW(L,NY,NX)*VLNOB(L,NY,NX)*14.0
      H2PO4(L,NY,NX)=CH2PU(L,NY,NX)*VOLW(L,NY,NX)*VLPO4(L,NY,NX)*31.0
      H2POB(L,NY,NX)=CH2PU(L,NY,NX)*VOLW(L,NY,NX)*VLPOB(L,NY,NX)*31.0
      H1PO4(L,NY,NX)=CH1PU(L,NY,NX)*VOLW(L,NY,NX)*VLPO4(L,NY,NX)*31.0
      H1POB(L,NY,NX)=CH1PU(L,NY,NX)*VOLW(L,NY,NX)*VLPOB(L,NY,NX)*31.0
      ZNO2S(L,NY,NX)=0.0
      ZNO2B(L,NY,NX)=0.0
C     WRITE(*,444)'ZNH4S',NX,NY,L,ZNH4S(L,NY,NX),CN4U(L,NY,NX)
C    2,VOLW(L,NY,NX),VLNH4(L,NY,NX),H2PO4(L,NY,NX),CH2PU(L,NY,NX)
C    3,VLPO4(L,NY,NX)
444   FORMAT(A8,3I4,12E12.4)
C
C     INITIAL STATE VARIABLES FOR CATIONS, ANIONS AND ION PAIRS IN SOIL
C
      ZAL(L,NY,NX)=CALU(L,NY,NX)*VOLW(L,NY,NX)
      ZFE(L,NY,NX)=CFEU(L,NY,NX)*VOLW(L,NY,NX)
      ZHY(L,NY,NX)=CHYU(L,NY,NX)*VOLW(L,NY,NX)
      ZCA(L,NY,NX)=CCAU(L,NY,NX)*VOLW(L,NY,NX)
      ZMG(L,NY,NX)=CMGU(L,NY,NX)*VOLW(L,NY,NX)
      ZNA(L,NY,NX)=CNAU(L,NY,NX)*VOLW(L,NY,NX)
      ZKA(L,NY,NX)=CKAU(L,NY,NX)*VOLW(L,NY,NX)
      ZOH(L,NY,NX)=COHU(L,NY,NX)*VOLW(L,NY,NX)
      ZSO4(L,NY,NX)=CSOU(L,NY,NX)*VOLW(L,NY,NX)
      ZCL(L,NY,NX)=CCLU(L,NY,NX)*VOLW(L,NY,NX)
      ZCO3(L,NY,NX)=CC3U(L,NY,NX)*VOLW(L,NY,NX)
      ZHCO3(L,NY,NX)=CHCU(L,NY,NX)*VOLW(L,NY,NX)
      ZALOH1(L,NY,NX)=CAL1U(L,NY,NX)*VOLW(L,NY,NX)
      ZALOH2(L,NY,NX)=CAL2U(L,NY,NX)*VOLW(L,NY,NX)
      ZALOH3(L,NY,NX)=CAL3U(L,NY,NX)*VOLW(L,NY,NX)
      ZALOH4(L,NY,NX)=CAL4U(L,NY,NX)*VOLW(L,NY,NX)
      ZALS(L,NY,NX)=CALSU(L,NY,NX)*VOLW(L,NY,NX)
      ZFEOH1(L,NY,NX)=CFE1U(L,NY,NX)*VOLW(L,NY,NX)
      ZFEOH2(L,NY,NX)=CFE2U(L,NY,NX)*VOLW(L,NY,NX)
      ZFEOH3(L,NY,NX)=CFE3U(L,NY,NX)*VOLW(L,NY,NX)
      ZFEOH4(L,NY,NX)=CFE4U(L,NY,NX)*VOLW(L,NY,NX)
      ZFES(L,NY,NX)=CFESU(L,NY,NX)*VOLW(L,NY,NX)
      ZCAO(L,NY,NX)=CCAOU(L,NY,NX)*VOLW(L,NY,NX)
      ZCAC(L,NY,NX)=CCACU(L,NY,NX)*VOLW(L,NY,NX)
      ZCAH(L,NY,NX)=CCAHU(L,NY,NX)*VOLW(L,NY,NX)
      ZCAS(L,NY,NX)=CCASU(L,NY,NX)*VOLW(L,NY,NX)
      ZMGO(L,NY,NX)=CMGOU(L,NY,NX)*VOLW(L,NY,NX)
      ZMGC(L,NY,NX)=CMGCU(L,NY,NX)*VOLW(L,NY,NX)
      ZMGH(L,NY,NX)=CMGHU(L,NY,NX)*VOLW(L,NY,NX)
      ZMGS(L,NY,NX)=CMGSU(L,NY,NX)*VOLW(L,NY,NX)
      ZNAC(L,NY,NX)=CNACU(L,NY,NX)*VOLW(L,NY,NX)
      ZNAS(L,NY,NX)=CNASU(L,NY,NX)*VOLW(L,NY,NX)
      ZKAS(L,NY,NX)=CKASU(L,NY,NX)*VOLW(L,NY,NX)
      H0PO4(L,NY,NX)=CH0PU(L,NY,NX)*VOLW(L,NY,NX)*VLPO4(L,NY,NX)
      H3PO4(L,NY,NX)=CH3PU(L,NY,NX)*VOLW(L,NY,NX)*VLPO4(L,NY,NX)
      ZFE1P(L,NY,NX)=CF1PU(L,NY,NX)*VOLW(L,NY,NX)*VLPO4(L,NY,NX)
      ZFE2P(L,NY,NX)=CF2PU(L,NY,NX)*VOLW(L,NY,NX)*VLPO4(L,NY,NX)
      ZCA0P(L,NY,NX)=CC0PU(L,NY,NX)*VOLW(L,NY,NX)*VLPO4(L,NY,NX)
      ZCA1P(L,NY,NX)=CC1PU(L,NY,NX)*VOLW(L,NY,NX)*VLPO4(L,NY,NX)
      ZCA2P(L,NY,NX)=CC2PU(L,NY,NX)*VOLW(L,NY,NX)*VLPO4(L,NY,NX)
      ZMG1P(L,NY,NX)=CM1PU(L,NY,NX)*VOLW(L,NY,NX)*VLPO4(L,NY,NX)
      H0POB(L,NY,NX)=CH0PU(L,NY,NX)*VOLW(L,NY,NX)*VLPOB(L,NY,NX)
      H3POB(L,NY,NX)=CH3PU(L,NY,NX)*VOLW(L,NY,NX)*VLPOB(L,NY,NX)
      ZFE1PB(L,NY,NX)=CF1PU(L,NY,NX)*VOLW(L,NY,NX)*VLPOB(L,NY,NX)
      ZFE2PB(L,NY,NX)=CF2PU(L,NY,NX)*VOLW(L,NY,NX)*VLPOB(L,NY,NX)
      ZCA0PB(L,NY,NX)=CC0PU(L,NY,NX)*VOLW(L,NY,NX)*VLPOB(L,NY,NX)
      ZCA1PB(L,NY,NX)=CC1PU(L,NY,NX)*VOLW(L,NY,NX)*VLPOB(L,NY,NX)
      ZCA2PB(L,NY,NX)=CC2PU(L,NY,NX)*VOLW(L,NY,NX)*VLPOB(L,NY,NX)
      ZMG1PB(L,NY,NX)=CM1PU(L,NY,NX)*VOLW(L,NY,NX)*VLPOB(L,NY,NX)
C
C     INITIAL STATE VARIABLES FOR ALL MATERIAL IN SOIL MACROPORES
C
      CO2SH(L,NY,NX)=0.0
      CH4SH(L,NY,NX)=0.0
      OXYSH(L,NY,NX)=0.0
      Z2GSH(L,NY,NX)=0.0
      Z2OSH(L,NY,NX)=0.0
      H2GSH(L,NY,NX)=0.0
      ZNH4SH(L,NY,NX)=0.0
      ZNH4BH(L,NY,NX)=0.0
      ZNH3SH(L,NY,NX)=0.0
      ZNH3BH(L,NY,NX)=0.0
      ZNO3SH(L,NY,NX)=0.0
      ZNO3BH(L,NY,NX)=0.0
      ZNO2SH(L,NY,NX)=0.0
      ZNO2BH(L,NY,NX)=0.0
      H2PO4H(L,NY,NX)=0.0
      H2POBH(L,NY,NX)=0.0
      H1PO4H(L,NY,NX)=0.0
      H1POBH(L,NY,NX)=0.0
      ZALH(L,NY,NX)=0.0
      ZFEH(L,NY,NX)=0.0
      ZHYH(L,NY,NX)=0.0
      ZCCH(L,NY,NX)=0.0
      ZMAH(L,NY,NX)=0.0
      ZNAH(L,NY,NX)=0.0
      ZKAH(L,NY,NX)=0.0
      ZOHH(L,NY,NX)=0.0
      ZSO4H(L,NY,NX)=0.0
      ZCLH(L,NY,NX)=0.0
      ZCO3H(L,NY,NX)=0.0
      ZHCO3H(L,NY,NX)=0.0
      ZALO1H(L,NY,NX)=0.0
      ZALO2H(L,NY,NX)=0.0
      ZALO3H(L,NY,NX)=0.0
      ZALO4H(L,NY,NX)=0.0
      ZALSH(L,NY,NX)=0.0
      ZFEO1H(L,NY,NX)=0.0
      ZFEO2H(L,NY,NX)=0.0
      ZFEO3H(L,NY,NX)=0.0
      ZFEO4H(L,NY,NX)=0.0
      ZFESH(L,NY,NX)=0.0
      ZCAOH(L,NY,NX)=0.0
      ZCACH(L,NY,NX)=0.0
      ZCAHH(L,NY,NX)=0.0
      ZCASH(L,NY,NX)=0.0
      ZMGOH(L,NY,NX)=0.0
      ZMGCH(L,NY,NX)=0.0
      ZMGHH(L,NY,NX)=0.0
      ZMGSH(L,NY,NX)=0.0
      ZNACH(L,NY,NX)=0.0
      ZNASH(L,NY,NX)=0.0
      ZKASH(L,NY,NX)=0.0
      H0PO4H(L,NY,NX)=0.0
      H1PO4H(L,NY,NX)=0.0
      H3PO4H(L,NY,NX)=0.0
      ZFE1PH(L,NY,NX)=0.0
      ZFE2PH(L,NY,NX)=0.0
      ZCA0PH(L,NY,NX)=0.0
      ZCA1PH(L,NY,NX)=0.0
      ZCA2PH(L,NY,NX)=0.0
      ZMG1PH(L,NY,NX)=0.0
      H0POBH(L,NY,NX)=0.0
      H1POBH(L,NY,NX)=0.0
      H3POBH(L,NY,NX)=0.0
      ZFE1BH(L,NY,NX)=0.0
      ZFE2BH(L,NY,NX)=0.0
      ZCA0BH(L,NY,NX)=0.0
      ZCA1BH(L,NY,NX)=0.0
      ZCA2BH(L,NY,NX)=0.0
      ZMG1BH(L,NY,NX)=0.0
C
C     INITIAL STATE VARIABLES FOR EXCHANGEABLE CATIONS AND ANIONS
C
      XN4(L,NY,NX)=XN41*BKVL(L,NY,NX)*VLNH4(L,NY,NX)
      XNB(L,NY,NX)=XN41*BKVL(L,NY,NX)*VLNHB(L,NY,NX)
      XHY(L,NY,NX)=XHY1*BKVL(L,NY,NX)
      XAL(L,NY,NX)=XAL1*BKVL(L,NY,NX) 
      XFE(L,NY,NX)=XFE1*BKVL(L,NY,NX) 
      XCA(L,NY,NX)=XCA1*BKVL(L,NY,NX) 
      XMG(L,NY,NX)=XMG1*BKVL(L,NY,NX) 
      XNA(L,NY,NX)=XNA1*BKVL(L,NY,NX)
      XKA(L,NY,NX)=XKA1*BKVL(L,NY,NX)
      XHC(L,NY,NX)=XHC1*BKVL(L,NY,NX)
      XALO2(L,NY,NX)=XALO21*BKVL(L,NY,NX)
      XFEO2(L,NY,NX)=XFEO21*BKVL(L,NY,NX)
      XOH0(L,NY,NX)=XOH01*BKVL(L,NY,NX)*VLPO4(L,NY,NX)
      XOH1(L,NY,NX)=XOH11*BKVL(L,NY,NX)*VLPO4(L,NY,NX)
      XOH2(L,NY,NX)=XOH21*BKVL(L,NY,NX)*VLPO4(L,NY,NX)
      XH1P(L,NY,NX)=XH1P1*BKVL(L,NY,NX)*VLPO4(L,NY,NX)
      XH2P(L,NY,NX)=XH2P1*BKVL(L,NY,NX)*VLPO4(L,NY,NX)
      XOH0B(L,NY,NX)=XOH01*BKVL(L,NY,NX)*VLPOB(L,NY,NX)
      XOH1B(L,NY,NX)=XOH11*BKVL(L,NY,NX)*VLPOB(L,NY,NX)
      XOH2B(L,NY,NX)=XOH21*BKVL(L,NY,NX)*VLPOB(L,NY,NX)
      XH1PB(L,NY,NX)=XH1P1*BKVL(L,NY,NX)*VLPOB(L,NY,NX)
      XH2PB(L,NY,NX)=XH2P1*BKVL(L,NY,NX)*VLPOB(L,NY,NX)
C
C     INITIAL STATE VARIABLES FOR PRECIPITATES
C
      PALOH(L,NY,NX)=PALOH1*BKVL(L,NY,NX)
      PFEOH(L,NY,NX)=PFEOH1*BKVL(L,NY,NX)
      PCACO(L,NY,NX)=PCACO1*BKVL(L,NY,NX)
      PCASO(L,NY,NX)=PCASO1*BKVL(L,NY,NX)
      PALPO(L,NY,NX)=PALPO1*BKVL(L,NY,NX)*VLPO4(L,NY,NX)
      PFEPO(L,NY,NX)=PFEPO1*BKVL(L,NY,NX)*VLPO4(L,NY,NX)
      PCAPD(L,NY,NX)=PCAPD1*BKVL(L,NY,NX)*VLPO4(L,NY,NX)
      PCAPH(L,NY,NX)=PCAPH1*BKVL(L,NY,NX)*VLPO4(L,NY,NX)
      PCAPM(L,NY,NX)=0.0
      PALPB(L,NY,NX)=PALPO1*BKVL(L,NY,NX)*VLPOB(L,NY,NX)
      PFEPB(L,NY,NX)=PFEPO1*BKVL(L,NY,NX)*VLPOB(L,NY,NX)
      PCPDB(L,NY,NX)=PCAPD1*BKVL(L,NY,NX)*VLPOB(L,NY,NX)
      PCPHB(L,NY,NX)=PCAPH1*BKVL(L,NY,NX)*VLPOB(L,NY,NX)
      PCPMB(L,NY,NX)=0.0
      ECND(L,NY,NX)=0.0
      CSTR(L,NY,NX)=0.0
      CION(L,NY,NX)=0.0
      ZNH4S(L,NY,NX)=ZNH4S(L,NY,NX)+0.5*OSN(1,2,L,NY,NX)
      OSN(1,2,L,NY,NX)=OSN(1,2,L,NY,NX)-0.5*OSN(1,2,L,NY,NX)
C     IF(K.EQ.3)THEN
C     WRITE(*,2222)'XN4F',K,L,XN4(L,NY,NX),XN4Q,BKVL(L,NY,NX)
C    2,VLNH4(L,NY,NX)
C     WRITE(*,2222)'XOH2F',K,L,XOH21,XOH11,XOH01,XH2P1,XH1P1,XOH,XPT
C    2,CH3P1,CH2P1,CH1P1,CH0P1,AHY1,XHP,FHP3,FHP2,FHP1,FHP0
C    3,FXP2,FXP1,PALPO1,PFEPO1,PCAPD1,PCAPH1 
C     ENDIF
      ENDIF
2000  CONTINUE
1200  CONTINUE
C
C     INITIAL STATE VARIABLES FOR MINERALS IN SURFACE RESIDUE
C
      IF(DATA(20).EQ.'NO'.AND.IGO.EQ.0)THEN
      ZNH4S(0,NY,NX)=0.0
      ZNH3S(0,NY,NX)=0.0
      ZNO3S(0,NY,NX)=0.0
      ZNO2S(0,NY,NX)=0.0
      H2PO4(0,NY,NX)=0.0
      H1PO4(0,NY,NX)=0.0
      ZNH4B(0,NY,NX)=0.0
      ZNH3B(0,NY,NX)=0.0
      ZNO3B(0,NY,NX)=0.0
      ZNO2B(0,NY,NX)=0.0
      H2POB(0,NY,NX)=0.0
      H1POB(0,NY,NX)=0.0
      XN4(0,NY,NX)=0.0
      XNB(0,NY,NX)=0.0
      XOH0(0,NY,NX)=0.0
      XOH1(0,NY,NX)=0.0
      XOH2(0,NY,NX)=0.0
      XH1P(0,NY,NX)=0.0
      XH2P(0,NY,NX)=0.0
      XNB(0,NY,NX)=0.0
      XOH0B(0,NY,NX)=0.0
      XOH1B(0,NY,NX)=0.0
      XOH2B(0,NY,NX)=0.0
      XH1PB(0,NY,NX)=0.0
      XH2PB(0,NY,NX)=0.0
      PALPO(0,NY,NX)=0.0
      PFEPO(0,NY,NX)=0.0
      PCAPD(0,NY,NX)=0.0
      PCAPH(0,NY,NX)=0.0
      PCAPM(0,NY,NX)=0.0
      PALPB(0,NY,NX)=0.0
      PFEPB(0,NY,NX)=0.0
      PCPDB(0,NY,NX)=0.0
      PCPHB(0,NY,NX)=0.0
      PCPMB(0,NY,NX)=0.0
C
C     INITIAL STATE VARIABLES FOR MINERAL N AND P IN SNOWPACK
C
      DO 9985 L=1,JS
      IF(VHCPW(L,NY,NX).GT.VHCPWX(NY,NX))THEN
      VOLWW=VOLWSL(L,NY,NX)+VOLSSL(L,NY,NX)+VOLISL(L,NY,NX)*DENSI
      CO2W(L,NY,NX)=VOLWW*CCOR(NY,NX)
      CH4W(L,NY,NX)=VOLWW*CCHR(NY,NX)
      OXYW(L,NY,NX)=VOLWW*COXR(NY,NX)
      ZNGW(L,NY,NX)=VOLWW*CNNR(NY,NX)
      ZN2W(L,NY,NX)=VOLWW*CN2R(NY,NX)
      ZN4W(L,NY,NX)=VOLWW*CN4R(NY,NX)*14.0
      ZN3W(L,NY,NX)=VOLWW*CN3R(NY,NX)*14.0
      ZNOW(L,NY,NX)=VOLWW*CNOR(NY,NX)*14.0
      Z1PW(L,NY,NX)=VOLWW*CH1PR(NY,NX)*31.0
      ZHPW(L,NY,NX)=VOLWW*CPOR(NY,NX)*31.0
C
C     INITIAL STATE VARIABLES FOR CATIONS AND ANIONS IN SNOWPACK
C
      IF(ISALTG.NE.0)THEN
      ZALW(L,NY,NX)=VOLWW*CALR(NY,NX)
      ZFEW(L,NY,NX)=VOLWW*CFER(NY,NX)
      ZHYW(L,NY,NX)=VOLWW*CHYR(NY,NX)
      ZCAW(L,NY,NX)=VOLWW*CCAR(NY,NX)
      ZMGW(L,NY,NX)=VOLWW*CMGR(NY,NX)
      ZNAW(L,NY,NX)=VOLWW*CNAR(NY,NX)
      ZKAW(L,NY,NX)=VOLWW*CKAR(NY,NX)
      ZOHW(L,NY,NX)=VOLWW*COHR(NY,NX)
      ZSO4W(L,NY,NX)=VOLWW*CSOR(NY,NX)
      ZCLW(L,NY,NX)=VOLWW*CCLR(NY,NX)
      ZCO3W(L,NY,NX)=VOLWW*CC3R(NY,NX)
      ZHCO3W(L,NY,NX)=VOLWW*CHCR(NY,NX)
      ZALH1W(L,NY,NX)=VOLWW*CAL1R(NY,NX)
      ZALH2W(L,NY,NX)=VOLWW*CAL2R(NY,NX)
      ZALH3W(L,NY,NX)=VOLWW*CAL3R(NY,NX)
      ZALH4W(L,NY,NX)=VOLWW*CAL4R(NY,NX)
      ZALSW(L,NY,NX)=VOLWW*CALSR(NY,NX)
      ZFEH1W(L,NY,NX)=VOLWW*CFE1R(NY,NX)
      ZFEH2W(L,NY,NX)=VOLWW*CFE2R(NY,NX)
      ZFEH3W(L,NY,NX)=VOLWW*CFE3R(NY,NX)
      ZFEH4W(L,NY,NX)=VOLWW*CFE4R(NY,NX)
      ZFESW(L,NY,NX)=VOLWW*CFESR(NY,NX)
      ZCAOW(L,NY,NX)=VOLWW*CCAOR(NY,NX)
      ZCACW(L,NY,NX)=VOLWW*CCACR(NY,NX)
      ZCAHW(L,NY,NX)=VOLWW*CCAHR(NY,NX)
      ZCASW(L,NY,NX)=VOLWW*CCASR(NY,NX)
      ZMGOW(L,NY,NX)=VOLWW*CMGOR(NY,NX)
      ZMGCW(L,NY,NX)=VOLWW*CMGCR(NY,NX)
      ZMGHW(L,NY,NX)=VOLWW*CMGHR(NY,NX)
      ZMGSW(L,NY,NX)=VOLWW*CMGSR(NY,NX)
      ZNACW(L,NY,NX)=VOLWW*CNACR(NY,NX)
      ZNASW(L,NY,NX)=VOLWW*CNASR(NY,NX)
      ZKASW(L,NY,NX)=VOLWW*CKASR(NY,NX)
      H0PO4W(L,NY,NX)=VOLWW*CH0PR(NY,NX)
      H3PO4W(L,NY,NX)=VOLWW*CH3PR(NY,NX)
      ZFE1PW(L,NY,NX)=VOLWW*CF1PR(NY,NX)
      ZFE2PW(L,NY,NX)=VOLWW*CF2PR(NY,NX)
      ZCA0PW(L,NY,NX)=VOLWW*CC0PR(NY,NX)
      ZCA1PW(L,NY,NX)=VOLWW*CC1PR(NY,NX)
      ZCA2PW(L,NY,NX)=VOLWW*CC2PR(NY,NX)
      ZMG1PW(L,NY,NX)=VOLWW*CM1PR(NY,NX)
      ENDIF
      ELSE
      CO2W(L,NY,NX)=0.0
      CH4W(L,NY,NX)=0.0
      OXYW(L,NY,NX)=0.0
      ZNGW(L,NY,NX)=0.0
      ZN2W(L,NY,NX)=0.0
      ZN4W(L,NY,NX)=0.0
      ZN3W(L,NY,NX)=0.0
      ZNOW(L,NY,NX)=0.0
      Z1PW(L,NY,NX)=0.0
      ZHPW(L,NY,NX)=0.0
C
C     INITIAL STATE VARIABLES FOR CATIONS AND ANIONS IN SNOWPACK
C
      IF(ISALTG.NE.0)THEN
      ZALW(L,NY,NX)=0.0
      ZFEW(L,NY,NX)=0.0
      ZHYW(L,NY,NX)=0.0
      ZCAW(L,NY,NX)=0.0
      ZMGW(L,NY,NX)=0.0
      ZNAW(L,NY,NX)=0.0
      ZKAW(L,NY,NX)=0.0
      ZOHW(L,NY,NX)=0.0
      ZSO4W(L,NY,NX)=0.0 
      ZCLW(L,NY,NX)=0.0 
      ZCO3W(L,NY,NX)=0.0 
      ZHCO3W(L,NY,NX)=0.0
      ZALH1W(L,NY,NX)=0.0
      ZALH2W(L,NY,NX)=0.0
      ZALH3W(L,NY,NX)=0.0
      ZALH4W(L,NY,NX)=0.0
      ZALSW(L,NY,NX)=0.0 
      ZFEH1W(L,NY,NX)=0.0 
      ZFEH2W(L,NY,NX)=0.0
      ZFEH3W(L,NY,NX)=0.0
      ZFEH4W(L,NY,NX)=0.0
      ZFESW(L,NY,NX)=0.0
      ZCAOW(L,NY,NX)=0.0
      ZCACW(L,NY,NX)=0.0
      ZCAHW(L,NY,NX)=0.0
      ZCASW(L,NY,NX)=0.0
      ZMGOW(L,NY,NX)=0.0
      ZMGCW(L,NY,NX)=0.0
      ZMGHW(L,NY,NX)=0.0
      ZMGSW(L,NY,NX)=0.0
      ZNACW(L,NY,NX)=0.0
      ZNASW(L,NY,NX)=0.0
      ZKASW(L,NY,NX)=0.0
      H0PO4W(L,NY,NX)=0.0
      H3PO4W(L,NY,NX)=0.0
      ZFE1PW(L,NY,NX)=0.0
      ZFE2PW(L,NY,NX)=0.0
      ZCA0PW(L,NY,NX)=0.0
      ZCA1PW(L,NY,NX)=0.0
      ZCA2PW(L,NY,NX)=0.0
      ZMG1PW(L,NY,NX)=0.0
      ENDIF
      ENDIF
9985  CONTINUE
      ENDIF
9990  CONTINUE
9995  CONTINUE
      RETURN
      END



