      SUBROUTINE solute(I,J,NHW,NHE,NVN,NVS)
C
C     THIS SUBROUTINE CALCULATES ALL SOLUTE TRANSFORMATIONS
C     FROM THERMODYNAMIC EQUILIBRIA
C
      include "parameters.h"
      include "blkc.h"
      include "blk2a.h"
      include "blk2b.h"
      include "blk2c.h"
      include "blk8a.h"
      include "blk8b.h"
      include "blk10.h"
      include "blk11a.h"
      include "blk11b.h"
      include "blk13a.h"
      include "blk13b.h"
      include "blk13c.h"
      include "blk15a.h"
      include "blk15b.h"
      include "blk18a.h"
      include "blk18b.h"
      include "blk19a.h"
      include "blk19b.h"
      include "blk19c.h"
      include "blk19d.h"
      include "blk21a.h"
      include "blk21b.h"
C
C     EQUILIBRIUM CONSTANTS
C
      DIMENSION RNHUI(0:2)
C
C     DISSOLUTION AND DISSOCIATION CONSTANTS
C
C     DPH20=water,SPALO=AlOH3s,SPFEO=FeOH3s,SPCAC=CaCO3s,SPCAS=CaSO4s
C     SPALP=AlPO4s,SPFEP=FePO4s,SPCAM=Ca(H2PO4)2s,SPCAD=CAHPO4s
C     SPCAH=hydroxyapatite,SXOX2=R-OH2,SXOX1=R-OH,SXH2P=R-H2PO4
C     SXH1P=R-HPO4,DPCO2=CO2,DPHCO=HCO3,DPN4=NH4,DPAL1=ALOH,DPAL2=ALOH2
C     DPAL3=AlOH3,DPAL4=ALOH4,DPALS=AlSO4,DPFE1=FeOH,DPFE2=FeOH2,DPFE3=FeOH3
C     DPFE4=FeOH4,DPFES=FeSO4,DPCAO=CaOH,DPCAC=CaCO3,DPCAH=CaHCO3
C     DPCAS=CaSO4,DPMGO=MgOH,DPMGC=MgCO3,DPMGH=MgHCO3,DPMGS=MgSO4
C     DPNAC=NaCO3,DPNAS=NaSO4,DPKAS=KSO4,DPH1P=HPO4,DPH2P=H2PO4,DPH3P=H3PO4
C     DPF1P=FeHPO4,DPF2P=FeH2PO4,DPC0P=CaPO4,DPC1P=CaHPO4,DPC2P=CaH2PO4
C     DPM1P=MgHPO4,DPCOH=R-COO,DPALO=R-AlOH2,DPFEO=R-FeOH2
C
      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
     4,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)
C
C     MRXN=number of cycles for solving reaction equilibria
C     TPD,TADA,TADC,TSL=rate constants for dissoln,adsorpn,solute equilibria
C
      PARAMETER (MRXN=1,TPD=2.5E-03,TPDX=TPD/MRXN,TADA=5.0E-02
     2,TADAX=TADA/MRXN,TADC=5.0E-02,TADCX=TADC/MRXN
     3,TADC0=TADC*1.0E-02,TSL=5.0E-01,TSLX=TSL/MRXN)
C
C     DUKM,DUKI=Km,Ki for urea hydrolysis (SOIL SCI 136:56)
C     COOH=CEC of SOC
C     CCAMX=maximum Ca concentration
C     SPNH4,SPNH3,SPNHU,SPNO3,SPPO4=specific rate constants for 
C     NH4,NH3,urea,NO3,H2PO4 fertilizer dissolution
C     RNHUI=rate constants for decline in urea hydrolysis inhibition
C 
      PARAMETER (DUKM=1.0,DUKI=2.5,A0=1.0,COOH=2.5E-04
     2,CCAMX=10.0)
      PARAMETER (SPNH4=1.0E-00,SPNH3=1.0E-00,SPNHU=1.0E-01
     2,SPNO3=1.0E-00,SPPO4=5.0E-02)
      DATA RNHUI/10.0E-02,1.0E-02,0.5E-02/
      NPI=INT(NPH/2)
      DO 9995 NX=NHW,NHE
      DO 9990 NY=NVN,NVS
      DO 9985 L=NU(NY,NX),NL(NY,NX)
      IF(VOLX(L,NY,NX).GT.ZEROS2(NY,NX)
     2.AND.VOLWM(NPH,L,NY,NX).GT.ZEROS2(NY,NX))THEN
C
C     WATER VOLUME IN NON-BAND AND BAND SOIL ZONES
C
C     VOLWM=soil water volume
C     VLNH4,VLNHB=fractions of soil volume in NH4 non-band,band
C     VLNO3,VLNOB=fractions of soil volume in N03 non-band,band
C     VLPO4,VLPOB=fractions of soil volume in H2PO4 non-band,band
C     BKVL=soil mass
C
      VOLWNH=VOLWM(NPH,L,NY,NX)*VLNH4(L,NY,NX)
      VOLWNB=VOLWM(NPH,L,NY,NX)*VLNHB(L,NY,NX)
      VOLWNO=VOLWM(NPH,L,NY,NX)*VLNO3(L,NY,NX)
      VOLWNZ=VOLWM(NPH,L,NY,NX)*VLNOB(L,NY,NX)
      VOLWPO=VOLWM(NPH,L,NY,NX)*VLPO4(L,NY,NX)
      VOLWPB=VOLWM(NPH,L,NY,NX)*VLPOB(L,NY,NX)
      IF(BKVL(L,NY,NX).GT.ZEROS(NY,NX))THEN
      BKVLX=BKVL(L,NY,NX)
      BKVLNH=BKVL(L,NY,NX)*VLNH4(L,NY,NX)
      BKVLNB=BKVL(L,NY,NX)*VLNHB(L,NY,NX)
      BKVLNO=BKVL(L,NY,NX)*VLNO3(L,NY,NX)
      BKVLNZ=BKVL(L,NY,NX)*VLNOB(L,NY,NX)
      BKVLPO=BKVL(L,NY,NX)*VLPO4(L,NY,NX)
      BKVLPB=BKVL(L,NY,NX)*VLPOB(L,NY,NX)
      ELSE
      BKVLX=VOLA(L,NY,NX)
      BKVLNH=VOLWNH
      BKVLNB=VOLWNB
      BKVLNO=VOLWNO
      BKVLNZ=VOLWNZ
      BKVLPO=VOLWPO
      BKVLPB=VOLWPB
      ENDIF
C
C     UREA HYDROLYSIS IN BAND AND NON-BAND SOIL ZONES
C
C     VOLQ=biologically active soil water volume from nitro.f
C     COMA=concentration of active biomass
C     TOQCK=total microbial activity from nitro.f
C     DUKD=Km for urea hydrolysis
C
      IF(VOLQ(L,NY,NX).GT.ZEROS2(NY,NX))THEN
      COMA=AMIN1(0.1E+06,TOQCK(L,NY,NX)/VOLQ(L,NY,NX))
      ELSE
      COMA=0.1E+06
      ENDIF
      DUKD=DUKM*(1.0+COMA/DUKI)
C
C     UREA HYDROLYSIS INHIBITION
C
C     ZNHU0,ZNHUI=initial,current inhibition activity
C     RNHUI=rate constants for decline in urea hydrolysis inhibition
C
      IF(ZNHU0(L,NY,NX).GT.ZEROS(NY,NX)
     2.AND.ZNHUI(L,NY,NX).GT.ZEROS(NY,NX))THEN
      ZNHUI(L,NY,NX)=ZNHUI(L,NY,NX)
     2-RNHUI(IUTYP(NY,NX))*ZNHUI(L,NY,NX) 
     3*AMAX1(RNHUI(IUTYP(NY,NX)),1.0-ZNHUI(L,NY,NX)/ZNHU0(L,NY,NX))
      ELSE
      ZNHUI(L,NY,NX)=0.0
      ENDIF
C
C     UREA CONCENTRATION AND HYDROLYSIS IN NON-BAND
C
C     ZNHUFA=urea fertilizer in non-band
C     CNHUA=concentration of urea fertilizer in non-band
C     DFNSA=effect of microbial concentration on urea hydrolysis in non-band
C     RSNUA=rate of urea hydrolysis in non-band
C     SPNHU=specific rate constant for urea hydrolysis
C     TFNQ=temperature effect on microbial activity from nitro.f
C     ZNHUI=current inhibition activity
C
      IF(ZNHUFA(L,NY,NX).GT.ZEROS(NY,NX)
     2.AND.BKVL(L,NY,NX).GT.ZEROS(NY,NX))THEN
      CNHUA=ZNHUFA(L,NY,NX)/BKVL(L,NY,NX)
      ELSEIF(VOLW(L,NY,NX).GT.ZEROS2(NY,NX))THEN
      CNHUA=ZNHUFA(L,NY,NX)/VOLW(L,NY,NX)
      ELSE
      CNHUA=0.0
      ENDIF
      DFNSA=CNHUA/(CNHUA+DUKD)
      RSNUA=AMIN1(ZNHUFA(L,NY,NX)
     2,SPNHU*TOQCK(L,NY,NX)*DFNSA*TFNQ(L,NY,NX))*(1.0-ZNHUI(L,NY,NX))
C
C     UREA CONCENTRATION AND HYDROLYSIS IN BAND
C
C     ZNHUFB=urea fertilizer in band
C     CNHUB=concentration of urea fertilizer in band
C     DFNSB=effect of microbial concentration on urea hydrolysis in band
C     RSNUB=rate of urea hydrolysis in non-band
C     SPNHU=specific rate constant for urea hydrolysis
C     TFNQ=temperature effect on microbial activity from nitro.f
C
      IF(ZNHUFB(L,NY,NX).GT.ZEROS(NY,NX)
     2.AND.BKVL(L,NY,NX).GT.ZEROS(NY,NX))THEN
      CNHUB=ZNHUFB(L,NY,NX)/BKVL(L,NY,NX)
      ELSEIF(VOLW(L,NY,NX).GT.ZEROS2(NY,NX))THEN
      CNHUB=ZNHUFB(L,NY,NX)/VOLW(L,NY,NX)
      ELSE
      CNHUB=0.0
      ENDIF 
      DFNSB=CNHUB/(CNHUB+DUKD)
      RSNUB=AMIN1(ZNHUFB(L,NY,NX)
     2,SPNHU*TOQCK(L,NY,NX)*DFNSB*TFNQ(L,NY,NX))*(1.0-ZNHUI(L,NY,NX))
C     IF(ZNHUFA(L,NY,NX).GT.ZEROS(NY,NX))THEN
C     WRITE(*,8888)'UREA',I,J,L,IUTYP(NY,NX)
C    2,ZNHUFA(L,NY,NX),ZNHUFB(L,NY,NX),RSNUA,RSNUB
C    2,DFNSA,DFNSB,TFNQ(L,NY,NX),CNHUA,DUKD,DUKM,DUKI,TOQCK(L,NY,NX)
C    3,BKVL(L,NY,NX),SPNHU,ZNHU0(L,NY,NX),ZNHUI(L,NY,NX)
C    4,RNHUI(IUTYP(NY,NX)),VLNH4(L,NY,NX),VLNHB(L,NY,NX)
C    5,THETW(L,NY,NX)
8888  FORMAT(A8,4I4,40E12.4)
C     ENDIF
C
C     NH4, NH3, UREA, NO3 DISSOLUTION IN BAND AND NON-BAND
C     SOIL ZONES FROM FIRST-ORDER FUNCTIONS OF REMAINING
C     FERTILIZER (NOTE: SUPERPHOSPHATE AND ROCK PHOSPHATE
C     ARE REPRESENTED AS MONOCALCIUM PHOSPHATE AND HYDROXYAPATITE
C     MODELLED IN PHOSPHORUS REACTIONS BELOW)
C
C     RSN4AA,RSN4BA=rate of broadcast NH4 fertilizer dissoln in non-band,band
C     RSN3AA,RSN3BA=rate of broadcast NH3 fertilizer dissoln in non-band,band
C     RSNUAA,RSNUBA=rate of broadcast urea fertr dissoln in non-band,band
C     RSNOAA,RSNOBA=rate of broadcast NO3 fertilizer dissoln in non-band,band
C     RSN4BB= rate of banded NH4 fertilizer dissolution in band
C     RSN3BB= rate of banded NH3 fertilizer dissolution in band
C     RSNUBB= rate of banded urea fertilizer dissolution in band
C     RSNOBB= rate of banded NO3 fertilizer dissolution in band
C     ZNH4FA,ZNH3FA,ZNHUFA,ZNO3FA=broadcast NH4,NH3,urea,NO3 fertilizer
C     ZNH4FB,ZNH3FB,ZNHUFB,ZNO3FB=banded NH4,NH3,urea,NO3 fertilizer
C     VLNH4,VLNHB=fractions of soil volume in NH4 non-band,band
C     VLNO3,VLNOB=fractions of soil volume in N03 non-band,band
C     VLPO4,VLPOB=fractions of soil volume in H2PO4 non-band,band
C     THETW=soil water concentration
C
      RSN4AA=SPNH4*ZNH4FA(L,NY,NX)*VLNH4(L,NY,NX)
     2*THETW(L,NY,NX)
      RSN3AA=SPNH3*ZNH3FA(L,NY,NX)*VLNH4(L,NY,NX)
      RSNUAA=RSNUA*VLNH4(L,NY,NX)*THETW(L,NY,NX) 
      RSNOAA=SPNO3*ZNO3FA(L,NY,NX)*VLNO3(L,NY,NX)
     2*THETW(L,NY,NX)
      RSN4BA=SPNH4*ZNH4FA(L,NY,NX)*VLNHB(L,NY,NX)
     2*THETW(L,NY,NX)
      RSN3BA=SPNH3*ZNH3FA(L,NY,NX)*VLNHB(L,NY,NX)
      RSNUBA=RSNUA*VLNHB(L,NY,NX)*THETW(L,NY,NX) 
      RSNOBA=SPNO3*ZNO3FA(L,NY,NX)*VLNOB(L,NY,NX)
     2*THETW(L,NY,NX)
      RSN4BB=SPNH4*ZNH4FB(L,NY,NX)*THETW(L,NY,NX)
      RSN3BB=SPNH3*ZNH3FB(L,NY,NX) 
      RSNUBB=RSNUB*VLNHB(L,NY,NX)*THETW(L,NY,NX) 
      RSNOBB=SPNO3*ZNO3FB(L,NY,NX)*THETW(L,NY,NX)
C
C     SOLUBLE AND EXCHANGEABLE NH4 CONCENTRATIONS
C     IN NON-BAND AND BAND SOIL ZONES
C
C     VOLWNH,VOLWNB=water volume in NH4 non-band,band
C     RN4X,RN3X=NH4,NH3 input from uptake, mineraln, dissoln in non-band 
C     RNBX,R3BX=NH4,NH3 input from uptake, mineraln, dissoln in band 
C     TUPNH4,TUPNH3=soil-root exchange of NH4,NH3 in non-band from uptake.f
C     TUPNHB,TUPN3B=soil-root exchange of NH4,NH3 in band from uptake.f
C     XNH4S,XNH4B=net change in NH4 in band,non-band from nitro.f
C     CN41,CN31,CN4B,CN3B=total NH4,NH3 concentration in non-band,band
C     XN41,XN4B=adsorbed NH4 concentration in non-band,band
C
      IF(VOLWNH.GT.ZEROS2(NY,NX))THEN
      VOLWNX=14.0*VOLWNH 
      RN4X=(-TUPNH4(L,NY,NX)+XNH4S(L,NY,NX)+14.0*RSN4AA)/VOLWNX
      RN3X=(-TUPN3S(L,NY,NX)+14.0*RSNUAA)/VOLWNX
      CN41=AMAX1(0.0,ZNH4S(L,NY,NX)/VOLWNX+RN4X)
      CN31=AMAX1(0.0,ZNH3S(L,NY,NX)/VOLWNX+RN3X)
      XN41=AMAX1(0.0,XN4(L,NY,NX)/BKVLNH)
      ELSE
      RN4X=0.0
      RN3X=0.0
      CN41=0.0
      CN31=0.0
      XN41=0.0
      ENDIF
      IF(VOLWNB.GT.ZEROS2(NY,NX))THEN
      VOLWNX=14.0*VOLWNB 
      RNBX=(-TUPNHB(L,NY,NX)+XNH4B(L,NY,NX)+14.0*(RSN4BA+RSN4BB))
     2/VOLWNX
      R3BX=(-TUPN3B(L,NY,NX)+14.0*(RSNUBA+RSNUBB))
     2/VOLWNX 
      CN4B=AMAX1(0.0,ZNH4B(L,NY,NX)/VOLWNX+RNBX)
      CN3B=AMAX1(0.0,ZNH3B(L,NY,NX)/VOLWNX+R3BX)
      XN4B=AMAX1(0.0,XNB(L,NY,NX)/BKVLNB)
      ELSE
      RNBX=0.0
      R3BX=0.0
      CN4B=0.0
      CN3B=0.0
      XN4B=0.0
      ENDIF
C     IF(IYRC.EQ.2012.AND.I.EQ.151.AND.NX.EQ.1)THEN
C     WRITE(*,4141)'RN4X',I,J,NX,NY,L,RN4X,RN3X,RNBX,R3BX
C    2,CN41,CN31,CN4B,CN3B,ZNH4S(L,NY,NX),ZNH3S(L,NY,NX)
C    3,RSN4AA,TUPN3S(L,NY,NX),RSNUAA,TUPNHB(L,NY,NX)
C    4,XNH4B(L,NY,NX),RSN4BA,RSN4BB,TUPN3B(L,NY,NX)
C    5,RSNUBA,RSNUBB,ZNH4S(L,NY,NX),ZNH3S(L,NY,NX)
C    6,VOLWNX,BKVLNH 
4141  FORMAT(A8,5I4,30E12.4)
C     ENDIF
C
C     SOLUBLE, EXCHANGEABLE AND PRECIPITATED PO4 CONCENTRATIONS IN
C     NON-BAND AND BAND SOIL ZONES
C
C     VOLWPO,VOLWPB=water volume in H2PO4 non-band,band
C     RH1PX,RH2PX=HPO4,H2PO4 inputs from mineraln, uptake in non-band
C     RH1BX,RH2BX=HPO4,H2PO4 inputs from mineraln, uptake in band
C     XH1PS,XH1BS=net change in HPO4 in band,non-band from nitro.f 
C     TUPH1P,TUPH2P=soil-root exch of HPO4,H2PO4 in non-band from uptake.f
C     TUPH1B,TUPH2B=soil-root exch of HPO4,H2PO4 in band from uptake.f
C     CH1P1,CH2P1=HPO4,H2PO4 concentrations in non-band
C     CH1PB,CH2PB=HPO4,H2PO4 concentrations in band
C     XOH11,XOH11,XOH21=concn of adsorption sites R-,R-OH,R-OH2 in non-band
C     XOH1B,XH11B,XH21B=concn of adsorption sites R-,R-OH,R-OH2 in band
C     XH1P1,XH2P1=concentration of adsorbed HPO4,H2PO4 in non-band
C     XH11B,X2P1B=concentration of adsorbed HPO4,H2PO4 in band
C     PALPO1,PFEPO1=concn of precip AlPO4,FEPO4 in non-band
C     PALPOB,PFEPOB=concn of precip AlPO4,FEPO4 in band
C     PCAPM1,PCAPD1,PCAPH1=concn of precip CaH2PO4,CaHPO4,apatite in non-band
C     PCAPMB,PCAPDB,PCAPHB=concn of precip CaH2PO4,CaHPO4,apatite in band
C
      IF(VOLWPO.GT.ZEROS2(NY,NX))THEN
      VOLWPX=31.0*VOLWPO
      RH1PX=(XH1PS(L,NY,NX)-TUPH1P(L,NY,NX))/VOLWPX
      RH2PX=(XH2PS(L,NY,NX)-TUPH2P(L,NY,NX))/VOLWPX
      CH1P1=AMAX1(0.0,H1PO4(L,NY,NX)/VOLWPX+RH1PX)
      CH2P1=AMAX1(0.0,H2PO4(L,NY,NX)/VOLWPX+RH2PX)
      XOH01=AMAX1(0.0,XOH0(L,NY,NX))/BKVLPO
      XOH11=AMAX1(0.0,XOH1(L,NY,NX))/BKVLPO
      XOH21=AMAX1(0.0,XOH2(L,NY,NX))/BKVLPO
      XH1P1=AMAX1(0.0,XH1P(L,NY,NX))/BKVLPO
      XH2P1=AMAX1(0.0,XH2P(L,NY,NX))/BKVLPO
      PALPO1=AMAX1(0.0,PALPO(L,NY,NX))/BKVLPO
      PFEPO1=AMAX1(0.0,PFEPO(L,NY,NX))/BKVLPO
      PCAPM1=AMAX1(0.0,PCAPM(L,NY,NX))/BKVLPO
      PCAPD1=AMAX1(0.0,PCAPD(L,NY,NX))/BKVLPO
      PCAPH1=AMAX1(0.0,PCAPH(L,NY,NX))/BKVLPO
C     WRITE(*,8642)'CH2P1',I,J,L,CH2P1,H2PO4(L,NY,NX)
C    2,VOLWPX,RH2PX,XH2PS(L,NY,NX),TUPH2P(L,NY,NX)
8642  FORMAT(A8,3I4,20E12.4)
      ELSE
      RH1PX=0.0
      RH2PX=0.0
      CH1P1=0.0
      CH2P1=0.0
      XOH01=0.0
      XOH11=0.0
      XOH21=0.0
      XH1P1=0.0
      XH2P1=0.0
      PALPO1=0.0
      PFEPO1=0.0
      PCAPM1=0.0
      PCAPD1=0.0
      PCAPH1=0.0
      ENDIF
      IF(VOLWPB.GT.ZEROS2(NY,NX))THEN
      VOLWPX=31.0*VOLWPB
      RH1BX=(XH1BS(L,NY,NX)-TUPH1B(L,NY,NX))/VOLWPX
      RH2BX=(XH2BS(L,NY,NX)-TUPH2B(L,NY,NX))/VOLWPX
      CH1PB=AMAX1(0.0,H1POB(L,NY,NX)/VOLWPX+RH1BX)
      CH2PB=AMAX1(0.0,H2POB(L,NY,NX)/VOLWPX+RH2BX)
      XH01B=AMAX1(0.0,XOH0B(L,NY,NX))/BKVLPB
      XH11B=AMAX1(0.0,XOH1B(L,NY,NX))/BKVLPB
      XH21B=AMAX1(0.0,XOH2B(L,NY,NX))/BKVLPB
      X1P1B=AMAX1(0.0,XH1PB(L,NY,NX))/BKVLPB
      X2P1B=AMAX1(0.0,XH2PB(L,NY,NX))/BKVLPB
      PALPOB=AMAX1(0.0,PALPB(L,NY,NX))/BKVLPB
      PFEPOB=AMAX1(0.0,PFEPB(L,NY,NX))/BKVLPB
      PCAPMB=AMAX1(0.0,PCPMB(L,NY,NX))/BKVLPB
      PCAPDB=AMAX1(0.0,PCPDB(L,NY,NX))/BKVLPB
      PCAPHB=AMAX1(0.0,PCPHB(L,NY,NX))/BKVLPB
      ELSE
      RH1BX=0.0
      RH2BX=0.0
      CH1PB=0.0
      CH2PB=0.0
      XH01B=0.0
      XH11B=0.0
      XH21B=0.0
      X1P1B=0.0
      X2P1B=0.0
      PALPOB=0.0
      PFEPOB=0.0
      PCAPMB=0.0
      PCAPDB=0.0
      PCAPHB=0.0
      ENDIF
C
C     IF SALT OPTION SELECTED IN SITE FILE
C     THEN SOLVE FULL SET OF EQUILIBRIA REACTIONS
C
C     ISALTG=salt flag from site file
C
      IF(ISALTG.NE.0)THEN
C
C     salt code: *HY*=H+,*OH*=OH-,*AL*=Al3+,*FE*=Fe3+,*CA*=Ca2+,*MG*=Mg2+
C          :*NA*=Na+,*KA*=K+,*SO4*=SO42-,*CL*=Cl-,*CO3*=CO32-,*HCO3*=HCO3-
C          :*CO2*=CO2,*ALO1*=AlOH2-,*ALOH2=AlOH2-,*ALOH3*=AlOH3
C          :*ALOH4*=AlOH4+,*ALS*=AlSO4+,*FEO1*=FeOH2-,*FEOH2=F3OH2-
C          :*FEOH3*=FeOH3,*FEOH4*=FeOH4+,*FES*=FeSO4+,*CAO*=CaOH
C          :*CAC*=CaCO3,*CAH*=CaHCO3-,*CAS*=CaSO4,*MGO*=MgOH,*MGC*=MgCO3
C          :*MHG*=MgHCO3-,*MGS*=MgSO4,*NAC*=NaCO3-,*NAS*=NaSO4-,*KAS*=KSO4-
C     phosphorus code: *H0P*=PO43-,*H3P*=H3PO4,*F1P*=FeHPO42-,*F2P*=F1H2PO4-
C          :*C0P*=CaPO4-,*C1P*=CaHPO4,*C2P*=CaH2PO4+,*M1P*=MgHPO4,*COO*=COOH-
C          :*1=non-band,*B=band
C     C*,X*,Z*=soluble,exchangeable concentration, mass
C
C     SOLUBLE NO3 CONCENTRATIONS
C     IN NON-BAND AND BAND SOIL ZONES
C
C     VOLWNO,VOLWNZ=soil water volume in NO3 non-band,band
C     ZNO3S,ZNO3B=NO3 mass in non-band,band
C     CNO1,CNOB=NO3 concentrations in non-band,band
C
      IF(VOLWNO.GT.ZEROS2(NY,NX))THEN
      CNO1=AMAX1(0.0,ZNO3S(L,NY,NX)/(14.0*VOLWNO))
      ELSE
      CNO1=0.0
      ENDIF
      IF(VOLWNZ.GT.ZEROS2(NY,NX))THEN
      CNOB=AMAX1(0.0,ZNO3B(L,NY,NX)/(14.0*VOLWNZ))
      ELSE
      CNOB=0.0
      ENDIF
C
C     H CONCENTRATION
C
C     XZHYS=total H+ production from nitro.f
C
      CHY1=AMAX1(0.0,ZHY(L,NY,NX)+XZHYS(L,NY,NX))
     2/VOLWM(NPH,L,NY,NX)
C
C     SOLUTE ION AND ION PAIR CONCENTRATIONS
C
C     CCEC,XCEC=cation exchange concentration,capacity
C     BKVLX=soil mass
C     VOLWM=soil water volume
C
      IF(BKVLX.GT.ZEROS(NY,NX))THEN
      CCEC=AMAX1(ZERO,XCEC(L,NY,NX)/BKVLX)
      ELSE
      CCEC=ZERO
      ENDIF
      IF(VOLWM(NPH,L,NY,NX).GT.ZEROS2(NY,NX))THEN
      COH1=AMAX1(0.0,ZOH(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CAL1=AMAX1(0.0,ZAL(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CFE1=AMAX1(0.0,ZFE(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CCA1=AMAX1(0.0,ZCA(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CMG1=AMAX1(0.0,ZMG(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CNA1=AMAX1(0.0,ZNA(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CKA1=AMAX1(0.0,ZKA(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CSO41=AMAX1(0.0,ZSO4(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CCL1=AMAX1(0.0,ZCL(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CCO31=AMAX1(0.0,ZCO3(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CHCO31=AMAX1(0.0,ZHCO3(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CCO21=AMAX1(0.0,CO2S(L,NY,NX)/(12.0*VOLWM(NPH,L,NY,NX)))
      CALO1=AMAX1(0.0,ZALOH1(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CALO2=AMAX1(0.0,ZALOH2(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CALO3=AMAX1(0.0,ZALOH3(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CALO4=AMAX1(0.0,ZALOH4(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CALS1=AMAX1(0.0,ZALS(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CFEO1=AMAX1(0.0,ZFEOH1(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CFEO2=AMAX1(0.0,ZFEOH2(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CFEO3=AMAX1(0.0,ZFEOH3(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CFEO4=AMAX1(0.0,ZFEOH4(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CFES1=AMAX1(0.0,ZFES(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CCAO1=AMAX1(0.0,ZCAO(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CCAC1=AMAX1(0.0,ZCAC(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CCAH1=AMAX1(0.0,ZCAH(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CCAS1=AMAX1(0.0,ZCAS(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CMGO1=AMAX1(0.0,ZMGO(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CMGC1=AMAX1(0.0,ZMGC(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CMGH1=AMAX1(0.0,ZMGH(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CMGS1=AMAX1(0.0,ZMGS(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CNAC1=AMAX1(0.0,ZNAC(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CNAS1=AMAX1(0.0,ZNAS(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CKAS1=AMAX1(0.0,ZKAS(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      ELSE
      COH1=0.0
      CAL1=0.0
      CFE1=0.0
      CCA1=0.0
      CMG1=0.0
      CNA1=0.0
      CKA1=0.0
      CSO41=0.0
      CCL1=0.0
      CCO31=0.0
      CHCO31=0.0
      CCO21=0.0
      CALO1=0.0
      CALO2=0.0
      CALO3=0.0
      CALO4=0.0
      CALS1=0.0
      CFEO1=0.0
      CFEO2=0.0
      CFEO3=0.0
      CFEO4=0.0
      CFES1=0.0
      CCAO1=0.0
      CCAC1=0.0
      CCAH1=0.0
      CCAS1=0.0
      CMGO1=0.0
      CMGC1=0.0
      CMGH1=0.0
      CMGS1=0.0
      CNAC1=0.0
      CNAS1=0.0
      CKAS1=0.0
      ENDIF
C
C     PO4 CONCENTRATIONS IN NON-BAND AND BAND SOIL ZONES
C
C     VOLWPO,VOLWPB=water volume in PO4 non-band,band
C
      IF(VOLWPO.GT.ZEROS2(NY,NX))THEN
      VOLWPX=31.0*VOLWPO
      CH0P1=AMAX1(0.0,H0PO4(L,NY,NX)/VOLWPO)
      CH3P1=AMAX1(0.0,H3PO4(L,NY,NX)/VOLWPO)
      CF1P1=AMAX1(0.0,ZFE1P(L,NY,NX)/VOLWPO)
      CF2P1=AMAX1(0.0,ZFE2P(L,NY,NX)/VOLWPO)
      CC0P1=AMAX1(0.0,ZCA0P(L,NY,NX)/VOLWPO)
      CC1P1=AMAX1(0.0,ZCA1P(L,NY,NX)/VOLWPO)
      CC2P1=AMAX1(0.0,ZCA2P(L,NY,NX)/VOLWPO)
      CM1P1=AMAX1(0.0,ZMG1P(L,NY,NX)/VOLWPO)
      ELSE
      CH0P1=0.0
      CH3P1=0.0
      CF1P1=0.0
      CF2P1=0.0
      CC0P1=0.0
      CC1P1=0.0
      CC2P1=0.0
      CM1P1=0.0
      ENDIF
      IF(VOLWPB.GT.ZEROS2(NY,NX))THEN
      CH0PB=AMAX1(0.0,H0POB(L,NY,NX)/VOLWPB)
      CH3PB=AMAX1(0.0,H3POB(L,NY,NX)/VOLWPB)
      CF1PB=AMAX1(0.0,ZFE1PB(L,NY,NX)/VOLWPB)
      CF2PB=AMAX1(0.0,ZFE2PB(L,NY,NX)/VOLWPB)
      CC0PB=AMAX1(0.0,ZCA0PB(L,NY,NX)/VOLWPB)
      CC1PB=AMAX1(0.0,ZCA1PB(L,NY,NX)/VOLWPB)
      CC2PB=AMAX1(0.0,ZCA2PB(L,NY,NX)/VOLWPB)
      CM1PB=AMAX1(0.0,ZMG1PB(L,NY,NX)/VOLWPB)
      ELSE
      CH0PB=0.0
      CH3PB=0.0
      CF1PB=0.0
      CF2PB=0.0
      CC0PB=0.0
      CC1PB=0.0
      CC2PB=0.0
      CM1PB=0.0
      ENDIF
C
C     EXCHANGEABLE ION CONCENTRATIONS
C
      IF(BKVLX.GT.ZEROS(NY,NX))THEN
      XHY1=AMAX1(0.0,XHY(L,NY,NX)/BKVLX)
      XAL1=AMAX1(0.0,XAL(L,NY,NX)/BKVLX)
      XFE1=AMAX1(0.0,XFE(L,NY,NX)/BKVLX)
      XCA1=AMAX1(0.0,XCA(L,NY,NX)/BKVLX)
      XMG1=AMAX1(0.0,XMG(L,NY,NX)/BKVLX)
      XNA1=AMAX1(0.0,XNA(L,NY,NX)/BKVLX)
      XKA1=AMAX1(0.0,XKA(L,NY,NX)/BKVLX)
      XHC1=AMAX1(0.0,XHC(L,NY,NX)/BKVLX)
      XALO21=AMAX1(0.0,XALO2(L,NY,NX)/BKVLX)
      XFEO21=AMAX1(0.0,XFEO2(L,NY,NX)/BKVLX)
      XCOOH=AMAX1(0.0,COOH*ORGC(L,NY,NX)/BKVLX)
C
C     PRECIPITATE CONCENTRATIONS
C
      PALOH1=AMAX1(0.0,PALOH(L,NY,NX)/BKVLX)
      PFEOH1=AMAX1(0.0,PFEOH(L,NY,NX)/BKVLX)
      PCACO1=AMAX1(0.0,PCACO(L,NY,NX)/BKVLX)
      PCASO1=AMAX1(0.0,PCASO(L,NY,NX)/BKVLX)
      ELSE
      XHY1=0.0
      XAL1=0.0
      XFE1=0.0
      XCA1=0.0
      XMG1=0.0
      XNA1=0.0
      XKA1=0.0
      XHC1=0.0
      XALO21=0.0
      XFEO21=0.0
      XCOOH=0.0
      PALOH1=0.0
      PFEOH1=0.0
      PCACO1=0.0
      PCASO1=0.0
      ENDIF
C
C     CONVERGENCE TOWARDS SOLUTE EQILIBRIA
C
      DO 1000 M=1,MRXN
C
C     SOLUTE CONCENTRATIONS
C
      CN41=AMAX1(ZERO,CN41)
      CN4B=AMAX1(ZERO,CN4B)
      CN31=AMAX1(ZERO,CN31)
      CN3B=AMAX1(ZERO,CN3B)
      CHY1=AMAX1(ZERO,10.0**(-(PH(L,NY,NX)-3.0)))
      COH1=AMAX1(ZERO,DPH2O/CHY1)
      CCO31=AMAX1(ZERO,CCO31)
      CAL1=AMAX1(ZERO,CAL1)
      CFE1=AMAX1(ZERO,CFE1)
      CCA1=AMAX1(ZERO,AMIN1(CCAMX,CCA1))
      CMG1=AMAX1(ZERO,CMG1)
      CNA1=AMAX1(ZERO,CNA1)
      CKA1=AMAX1(ZERO,CKA1)
      CSO41=AMAX1(ZERO,CSO41)
      CHCO31=AMAX1(ZERO,CHCO31)
      CCO21=AMAX1(ZERO,CCO21)
      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)
      CH0PB=AMAX1(ZERO,CH0PB)
      CH1PB=AMAX1(ZERO,CH1PB)
      CH2PB=AMAX1(ZERO,CH2PB)
      CH3PB=AMAX1(ZERO,CH3PB)
      CF1PB=AMAX1(ZERO,CF1PB)
      CF2PB=AMAX1(ZERO,CF2PB)
      CC0PB=AMAX1(ZERO,CC0PB)
      CC1PB=AMAX1(ZERO,CC1PB)
      CC2PB=AMAX1(ZERO,CC2PB)
      CM1PB=AMAX1(ZERO,CM1PB)
      XCOO=AMAX1(0.0,XCOOH-XHC1-XALO21-XFEO21)
C
C     IONIC STRENGTH FROM SUMS OF ION CONCENTRATIONS
C
C     CC3,CA3,CC2,CA2,CC1,CA1=total tri-,di-,univalent cations C,anions A
C     CSTR1=ion strength
C
      CC3=CAL1+CFE1
      CA3=CH0P1*VLPO4(L,NY,NX)+CH0PB*VLPOB(L,NY,NX)
      CC2=CCA1+CMG1+CALO1+CFEO1+CF2P1*VLPO4(L,NY,NX)
     2+CF2PB*VLPOB(L,NY,NX)
      CA2=CSO41+CCO31+CH1P1*VLPO4(L,NY,NX)+CH1PB*VLPOB(L,NY,NX)
      CC1=CN41*VLNH4(L,NY,NX)+CN4B*VLNHB(L,NY,NX)+CHY1+CNA1+CKA1
     2+CALO2+CFEO2+CALS1+CFES1+CCAO1+CCAH1+CMGO1+CMGH1
     3+(CF1P1+CC2P1)*VLPO4(L,NY,NX)+(CF1PB+CC2PB)*VLPOB(L,NY,NX)
      CA1=CNO1*VLNO3(L,NY,NX)+CNOB*VLNOB(L,NY,NX)+COH1+CHCO31+CCL1
     2+CALO4+CFEO4+CNAC1+CNAS1+CKAS1+(CH2P1+CC0P1)*VLPO4(L,NY,NX)
     3+(CH2PB+CC0PB)*VLPOB(L,NY,NX)
      CSTR1=AMAX1(0.0,0.5E-03*(9.0*(CC3+CA3)+4.0*(CC2+CA2)
     2+CC1+CA1))
      CSTR2=SQRT(CSTR1)
      FSTR2=CSTR2/(1.0+CSTR2)
C
C     ACTIVITY COEFFICIENTS CALCULATED FROM ION STRENGTH
C
      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 REACTIONS IN NON-BAND, BAND
C     CALCULATED FROM ACTIVITIES OF REACTANTS AND PRODUCTS THROUGH 
C     CONVERGENCE SOLUTIONS FOR THEIR EQUILIBRIUM CONCENTRATIONS USING 
C     SOLUTE FORMS CURRENTLY AT HIGHEST CONCENTRATIONS
C
C     for all reactions:
C     A*=ion activity
C     PX,PY=solute forms with greatest activity
C     R*,P*=reactant,product
C     NR*,NP*=reactant,product stoichiometry
C     SP=solubility product of PX from parameters above
C     SPX=equilibrium product concentration
C     R*X=precipitation-dissolution rate
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
      AH0PB=CH0PB*A3
      AH1PB=CH1PB*A2
      AH2PB=CH2PB*A1
      AH3PB=CH3PB*A0
      AF1PB=CF1PB*A2
      AF2PB=CF2PB*A1
      AC0PB=CC0PB*A1
      AC1PB=CC1PB*A0
      AC2PB=CC2PB*A1
      AM1PB=CM1PB*A0
      AN41=CN41*A1
      AN4B=CN4B*A1
      AN31=CN31*A0
      AN3B=CN3B*A0
      AMG1=CMG1*A2
      ANA1=CNA1*A1
      AKA1=CKA1*A1
      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
      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
C     R1=AMAX1(ZERO,R1)
C     P1=AMAX1(ZERO,P1)
C     P2=AMAX1(ZERO,P2)
      SPX=SP*R1**NR1/P2**NP2
      RPALOX=AMAX1(-PALOH1,TPDX*(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(I.EQ.256)THEN
C     WRITE(*,1112)'ALOH',I,J,L,M,PALOH1,AAL1,AALO1,AALO2,AALO3,AALO4
C    2,AOH1,AHY1,PX,R1,P1,P2,SP,SPX,RPALOX,RHAL1,RHALO1,RHALO2
C    3,RHALO3,RHALO4,AAL1*AOH1**3,SPALO 
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,TPDX*(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(I.EQ.180.AND.J.EQ.12)THEN
C     WRITE(*,1112)'FEOH',I,J,L,M,PFEOH1,AFE1,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,TPDX*(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,TPDX*(S0-SQRT(S1)))
C     IF((M/10)*10.EQ.M)THEN
C     WRITE(*,1112)'CALC',I,J,L,M,PCASO1,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
      IF(VOLWPO.GT.ZEROS2(NY,NX))THEN
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,TPDX*(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(I.EQ.180.AND.J.EQ.12)THEN
C     WRITE(*,1112)'ALPO4',I,J,L,M,PALPO1,AAL1,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 
C     ENDIF
1112  FORMAT(A8,4I5,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,TPDX*(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(I.EQ.180.AND.J.EQ.12)THEN
C     WRITE(*,1112)'FEPO4',I,J,L,M,PFEPO1,AFE1,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,TPDX*(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,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,TPDX*(P2-SPX))
      IF(PX.EQ.AH1P1)THEN
      RHCAH1=RPCAHX
      ELSEIF(PX.EQ.AH2P1)THEN
      RHCAH2=RPCAHX
      ENDIF
C     IF((I/10)*10.EQ.I.AND.J.EQ.12)THEN
C     WRITE(*,1112)'A1',I,L,K,M,A1,A2,A3,FSTR2,CSTR1
C    2,CSTR2,CC3,CA3,CC2,CA2,CC1,CA1,VOLWM(NPH,L,NY,NX) 
C     WRITE(*,1112)'APATITE',I,J,L,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
C     ENDIF
C
C     MONOCALCIUM PHOSPHATE
C
      P1=ACA1
      P2=AH2P1
      P1=AMAX1(ZERO,P1)
      P2=AMAX1(ZERO,P2)
      SPX=SPCAM 
      S0=P1+P2
      S1=AMAX1(0.0,S0**2-4.0*(P1*P2-SPX))
      RPCAMX=AMAX1(-PCAPM1,TPDX*(S0-SQRT(S1)))
      ELSE
      RPALPX=0.0
      RPFEPX=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
      RPCAMX=0.0
      ENDIF
C
C     PHOSPHORUS PRECIPITATION-DISSOLUTION IN BAND SOIL ZONE
C
      IF(VOLWPB.GT.ZEROS2(NY,NX))THEN
C
C     ALUMINUM PHOSPHATE (VARISCITE)
C
      PX=AMAX1(AAL1,AALO1,AALO2,AALO3,AALO4)
      PY=AMAX1(AH1PB,AH2PB)
      R1=AHY1
      P3=AHY1
      IF(PY.EQ.AH1PB)THEN
      P2=AH1PB
      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=AH2PB
      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
      RHA0B1=0.0
      RHA1B1=0.0
      RHA2B1=0.0
      RHA3B1=0.0
      RHA4B1=0.0
      RHA0B2=0.0
      RHA1B2=0.0
      RHA2B2=0.0
      RHA3B2=0.0
      RHA4B2=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))
      RPALBX=AMAX1(-PALPOB,TPDX*(S0-SQRT(S1)))
      IF(PY.EQ.AH1PB)THEN
      IF(PX.EQ.AAL1)THEN
      RHA0B1=RPALBX
      ELSEIF(PX.EQ.AALO1)THEN
      RHA1B1=RPALBX
      ELSEIF(PX.EQ.AALO2)THEN
      RHA2B1=RPALBX
      ELSEIF(PX.EQ.AALO3)THEN
      RHA3B1=RPALBX
      ELSEIF(PX.EQ.AALO4)THEN
      RHA4B1=RPALBX
      ENDIF
      ELSE
      IF(PX.EQ.AAL1)THEN
      RHA0B2=RPALBX
      ELSEIF(PX.EQ.AALO1)THEN
      RHA1B2=RPALBX
      ELSEIF(PX.EQ.AALO2)THEN
      RHA2B2=RPALBX
      ELSEIF(PX.EQ.AALO3)THEN
      RHA3B2=RPALBX
      ELSEIF(PX.EQ.AALO4)THEN
      RHA4B2=RPALBX
      ENDIF
      ENDIF
C
C     IRON PHOSPHATE (STRENGITE)
C
      PX=AMAX1(AFE1,AFEO1,AFEO2,AFEO3,AFEO4)
      PY=AMAX1(AH1PB,AH2PB)
      R1=AHY1
      P3=AHY1
      IF(PY.EQ.AH1PB)THEN
      P2=AH1PB
      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=AH2PB
      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
      RHF0B1=0.0
      RHF1B1=0.0
      RHF2B1=0.0
      RHF3B1=0.0
      RHF4B1=0.0
      RHF0B2=0.0
      RHF1B2=0.0
      RHF2B2=0.0
      RHF3B2=0.0
      RHF4B2=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))
      RPFEBX=AMAX1(-PFEPOB,TPDX*(S0-SQRT(S1)))
      IF(PY.EQ.AH1PB)THEN
      IF(PX.EQ.AFE1)THEN
      RHF0B1=RPFEBX
      ELSEIF(PX.EQ.AFEO1)THEN
      RHF1B1=RPFEBX
      ELSEIF(PX.EQ.AFEO2)THEN
      RHF2B1=RPFEBX
      ELSEIF(PX.EQ.AFEO3)THEN
      RHF3B1=RPFEBX
      ELSEIF(PX.EQ.AFEO4)THEN
      RHF4B1=RPFEBX
      ENDIF
      ELSE
      IF(PX.EQ.AFE1)THEN
      RHF0B2=RPFEBX
      ELSEIF(PX.EQ.AFEO1)THEN
      RHF1B2=RPFEBX
      ELSEIF(PX.EQ.AFEO2)THEN
      RHF2B2=RPFEBX
      ELSEIF(PX.EQ.AFEO3)THEN
      RHF3B2=RPFEBX
      ELSEIF(PX.EQ.AFEO4)THEN
      RHF4B2=RPFEBX
      ENDIF
      ENDIF
C
C     DICALCIUM PHOSPHATE
C
      PX=AMAX1(AH1PB,AH2PB)
      R1=AHY1
      P1=ACA1
      IF(PX.EQ.AH1PB)THEN
      P2=AH1PB
      NR1=0
      SP=SPCAD 
      ELSEIF(PX.EQ.AH2PB)THEN
      P2=AH2PB
      NR1=1
      SP=SHCAD2 
      ENDIF
      RPCDB1=0.0
      RHCDB2=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))
      RPCDBX=AMAX1(-PCAPDB,TPDX*(S0-SQRT(S1)))
      IF(PX.EQ.AH1PB)THEN
      RPCDB1=RPCDBX
      ELSEIF(PX.EQ.AH2PB)THEN
      RHCDB2=RPCDBX
      ENDIF
C
C     HYDROXYAPATITE
C
      PX=AMAX1(AH1PB,AH2PB)
      R1=AHY1
      P1=ACA1
      IF(PX.EQ.AH1PB)THEN
      P2=AH1PB
      NR1=4
      SP=SHCAH1 
      ELSEIF(PX.EQ.AH2PB)THEN
      P2=AH2PB
      NR1=7
      SP=SHCAH2 
      ENDIF
      RHCHB1=0.0
      RHCHB2=0.0
      R1=AMAX1(ZERO,R1)
      P1=AMAX1(ZERO,P1)
      P2=AMAX1(ZERO,P2)
      SPX=(SP*R1**NR1/P1**5)**0.333
      RPCHBX=AMAX1(-PCAPHB,TPDX*(P2-SPX))
      IF(PX.EQ.AH1PB)THEN
      RHCHB1=RPCHBX
      ELSEIF(PX.EQ.AH2PB)THEN
      RHCHB2=RPCHBX
      ENDIF
C
C     MONOCALCIUM PHOSPHATE
C
      P1=ACA1
      P2=AH2PB
      P1=AMAX1(ZERO,P1)
      P2=AMAX1(ZERO,P2)
      SPX=SPCAM 
      S0=P1+P2
      S1=AMAX1(0.0,S0**2-4.0*(P1*P2-SPX))
      RPCMBX=AMAX1(-PCAPMB,TPDX*(S0-SQRT(S1)))
      ELSE
      RPALBX=0.0
      RPFEBX=0.0
      RPCDBX=0.0
      RPCHBX=0.0
      RPCMBX=0.0
      RHA0B1=0.0
      RHA1B1=0.0
      RHA2B1=0.0
      RHA3B1=0.0
      RHA4B1=0.0
      RHA0B2=0.0
      RHA1B2=0.0
      RHA2B2=0.0
      RHA3B2=0.0
      RHA4B2=0.0
      RHF0B1=0.0
      RHF1B1=0.0
      RHF2B1=0.0
      RHF3B1=0.0
      RHF4B1=0.0
      RHF0B2=0.0
      RHF1B2=0.0
      RHF2B2=0.0
      RHF3B2=0.0
      RHF4B2=0.0
      RPCDB1=0.0
      RHCDB2=0.0
      RHCHB1=0.0
      RHCHB2=0.0
      ENDIF
C
C     PHOSPHORUS ANION EXCHANGE IN NON-BAND SOIL ZONE
C     CALCULATED FROM EXCHANGE EQUILIBRIA AMONG H2PO4-,
C     HPO4--, H+, OH- AND PROTONATED AND NON-PROTONATED -OH
C     EXCHANGE SITES
C
C     BKVL=soil mass
C     VOLWM=soil water volume
C     XAEC=anion exchange capacity
C     VOLWPO=soil water volume in non-band
C     TADAX=adsorption rate constant
C     RXOH2,RXOH1=OH2,OH exchange with R-OH2,R-OH in non-band
C     SXOH2,SXOH1=equilibrium constant for OH2,OH exchange with R-OH2,R-OH
C
      IF(VOLWM(NPH,L,NY,NX).GT.ZEROS2(NY,NX))THEN
      VOLWBK=AMIN1(1.0,BKVL(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      ELSE
      VOLWBK=1.0
      ENDIF
      IF(VOLWPO.GT.ZEROS2(NY,NX)
     2.AND.XAEC(L,NY,NX).GT.ZEROS(NY,NX))THEN
      RXOH2=TADAX*(XOH11*AHY1-SXOH2*XOH21)/(XOH11+SXOH2)*VOLWBK
      RXOH1=TADAX*(XOH01*AHY1-SXOH1*XOH11)/(XOH01+SXOH1)*VOLWBK 
C
C     H2PO4 EXCHANGE IN NON-BAND SOIL ZONE FROM CONVERGENCE
C     SOLUTION FOR EQUILIBRIUM AMONG H2PO4-, H+, OH-, X-OH
C     AND X-H2PO4
C
C     SPH2P,SXH2P=equilibrium constant for H2PO4 exchange with R-OH2,R-OH
C     RXH2P,RYH2P=H2PO4 exchange with R-OH2,R-OH in non-band
C
      SPH2P=SXH2P*DPH2O 
      RXH2P=TADAX*(XOH21*AH2P1-SPH2P*XH2P1)/(XOH21+SPH2P)*VOLWBK 
      RYH2P=TADAX*(XOH11*AH2P1-SXH2P*XH2P1*AOH1)/(XOH11+SXH2P)*VOLWBK 
C
C     HPO4 EXCHANGE IN NON-BAND SOIL ZONE FROM CONVERGENCE
C     SOLUTION FOR EQUILIBRIUM AMONG HPO4--, H+, OH-, X-OH
C     AND X-HPO4
C
C     SPH1P=equilibrium constant for HPO4 exchange with R-OH
C     RXH1P=HPO4 exchange with R-OH in non-band
C
      SPH1P=SXH1P*DPH2O/DPH2P
      RXH1P=TADAX*(XOH11*AH1P1-SPH1P*XH1P1)/(XOH11+SPH1P)*VOLWBK 
      ELSE
      RXOH2=0.0
      RXOH1=0.0
      RXH2P=0.0
      RYH2P=0.0
      RXH1P=0.0
      ENDIF
C
C     PHOSPHORUS ANION EXCHANGE IN BAND SOIL ZONE
C     CALCULATED FROM EXCHANGE EQUILIBRIA AMONG H2PO4-,
C     HPO4--, H+, OH- AND PROTONATED AND NON-PROTONATED -OH
C     EXCHANGE SITES
C
      IF(VOLWPB.GT.ZEROS2(NY,NX)
     2.AND.XAEC(L,NY,NX).GT.ZEROS(NY,NX))THEN
C
C     RXO2B,RXO1B=OH2,OH exchange with R-OH2,R-OH in band
C     SXOH2,SXOH1=equilibrium constant for OH2,OH exchange with R-OH2,R-OH
C
      RXO2B=TADAX*(XH11B*AHY1-SXOH2*XH21B)/(XH11B+SXOH2)*VOLWBK 
      RXO1B=TADAX*(XH01B*AHY1-SXOH1*XH11B)/(XH01B+SXOH1)*VOLWBK 
C
C     H2PO4 EXCHANGE IN BAND SOIL ZONE FROM CONVERGENCE
C     SOLUTION FOR EQUILIBRIUM AMONG H2PO4-, H+, OH-, X-OH
C     AND X-H2PO4
C
C     SPH2P,SXH2P=equilibrium constant for H2PO4 exchange with R-OH2,R-OH
C     RXH2B,RYH2B=H2PO4 exchange with R-OH2,R-OH in band
C
      SPH2P=SXH2P*DPH2O 
      RXH2B=TADAX*(XH21B*AH2PB-SPH2P*X2P1B)/(XH21B+SPH2P)*VOLWBK 
      RYH2B=TADAX*(XH11B*AH2PB-SXH2P*X2P1B*AOH1)/(XH11B+SXH2P)*VOLWBK 
C
C     HPO4 EXCHANGE IN BAND SOIL ZONE FROM CONVERGENCE
C     SOLUTION FOR EQUILIBRIUM AMONG HPO4--, H+, OH-, X-OH
C     AND X-HPO4
C
C     SPH1P=equilibrium constant for HPO4 exchange with R-OH
C     RXH1B=HPO4 exchange with R-OH in band
C
      SPH1P=SXH1P*DPH2O/DPH2P
      RXH1B=TADAX*(XH11B*AH1PB-SPH1P*X1P1B)/(XH11B+SPH1P)*VOLWBK 
C     WRITE(*,2226)'RXH1B',I,J,L,M,RXH1B,XH11B,XH21B,CH1PB
C    2,SPH1P,X1P1B,RYH2B,CH2PB,SXH2P,X2P1B,COH1,CHY1,ROH
2226  FORMAT(A8,4I4,20E12.4)
      ELSE
      RXO2B=0.0
      RXO1B=0.0
      RXH2B=0.0
      RYH2B=0.0
      RXH1B=0.0
      ENDIF
C
C     CATION EXCHANGE FROM GAPON SELECTIVITY COEFFICIENTS
C     FOR CA-NH4, CA-H, CA-AL, CA-MG, CA-NA, CA-K
C
      IF(XCEC(L,NY,NX).GT.ZEROS(NY,NX))THEN
C
C     CATION CONCENTRATIONS
C
C     EQUILIBRIUM X-CA CONCENTRATION FROM CEC, GAPON COEFFICIENTS
C     AND CATION CONCENTRATIONS
C
C     CCEC,XCEC=cation exchange concentration,capacity
C     XCAX=equilibrium R-Ca concentration
C     GKC4,GKCH,GKCA,GKCM,GKCN,GKCK=Gapon selectivity coefficients for
C     Ca-NH4,Ca-H,Ca-Al,Ca-Mg,Ca-Na,Ca-K
C     X*Q=equilibrium exchangeable concentrations
C     XTLQ=total equilibrium exchangeable concentration 
C
      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*VLNH4(L,NY,NX)
     2+GKC4(L,NY,NX)*AN4B/ACAX*VLNHB(L,NY,NX)
     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)
      XNBQ=XCAX*AN4B*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*VLNH4(L,NY,NX)+XNBQ*VLNHB(L,NY,NX)
     2+XHYQ+XALQ+XFEQ+XCAQ+XMGQ+XNAQ+XKAQ
      IF(XTLQ.GT.ZERO)THEN
      FX=CCEC/XTLQ
      ELSE
      FX=0.0
      ENDIF
      XN4Q=FX*XN4Q
      XNBQ=FX*XNBQ
      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
C
C     NH4 EXCHANGE IN NON-BAND AND BAND SOIL ZONES
C
C     RXN4,RXNB=NH4 adsorption in non-band,band
C     TADCX=adsorption rate constant
C
      RXN4=TADCX*AMAX1(AMIN1((XN4Q-XN41)*AN41/XN4Q,CN41),-XN41) 
      RXNB=TADCX*AMAX1(AMIN1((XNBQ-XN4B)*AN4B/XNBQ,CN4B),-XN4B)
C
C     H,AL,FE,CA,MG,NA,K EXCHANGE
C
C     RX*=ion adsorption
C
      RXHY=TADCX*AMIN1((XHYQ-XHY1)*AHY1/XHYQ,CHY1) 
      RXAL=TADCX*AMIN1((XALQ-XAL1)*AALX/XALQ,CAL1) 
      RXFE=TADCX*AMIN1((XFEQ-XFE1)*AFEX/XFEQ,CFE1) 
      RXCA=TADCX*AMIN1((XCAQ-XCA1)*ACAX/XCAQ,CCA1) 
      RXMG=TADCX*AMIN1((XMGQ-XMG1)*AMGX/XMGQ,CMG1) 
      RXNA=TADCX*AMIN1((XNAQ-XNA1)*ANA1/XNAQ,CNA1)
      RXKA=TADCX*AMIN1((XKAQ-XKA1)*AKA1/XKAQ,CKA1) 
C     IF(I.EQ.256.AND.L.EQ.1)THEN
C     WRITE(*,1112)'RXAL',I,J,L,M,RXAL,TADCX,XALQ,XAL1,CAL1
C    2,AAL1,AALX,CCEC,XCAX,GKCA(L,NY,NX),FX
C     ENDIF
      ELSE
      RXN4=0.0
      RXNB=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     SOLUTE DISSOCIATION REACTIONS
C
C     for all reactions:
C     S0,S1=equilibrium constant,equilibrium solute concentration**2
C     TSLX=dissociation rate constant
C
C     DISSOCIATION OF CARBOXYL RADICALS
C     AND ADSORPTION OF AL AND FE OH2
C
C     RXHC=COOH-COO+H dissociation
C     RXALO2,RXFLO2=Al(OH2)2+,FeOH2+ adsorption
C
      S0=AHY1+XCOO+DPCOH 
      S1=AMAX1(0.0,S0**2-4.0*(AHY1*XCOO-DPCOH*XHC1))
      RXHC=TADCX*(S0-SQRT(S1))
      S0=AALO2+XCOO+DPALO 
      S1=AMAX1(0.0,S0**2-4.0*(AALO2*XCOO-DPALO*XALO21))
      RXALO2=TADAX*(S0-SQRT(S1))
      S0=AFEO2+XCOO+DPFEO 
      S1=AMAX1(0.0,S0**2-4.0*(AFEO2*XCOO-DPFEO*XFEO21))
      RXFEO2=TADAX*(S0-SQRT(S1))
C
C     RNH4,RNHB-NH4-NH3+H dissociation in non-band,band
C     DPN4=NH4 dissociation constant
C
      IF(VOLWNH.GT.ZEROS2(NY,NX))THEN
      RNH4=TSLX*(AHY1*AN31-DPN4*AN41)/(DPN4+AHY1)
      ELSE
      RNH4=0.0
      ENDIF
      IF(VOLWNB.GT.ZEROS2(NY,NX))THEN
      RNHB=TSLX*(AHY1*AN3B-DPN4*AN4B)/(DPN4+AHY1)
      ELSE
      RNHB=0.0
      ENDIF
C
C     RCO2Q=CO2-HCO3+H dissociation
C
      S0=AHY1+AHCO31+DPCO2
      S1=AMAX1(0.0,S0**2-4.0*(AHY1*AHCO31-DPCO2*ACO21))
      RCO2Q=TSLX*(S0-SQRT(S1))
C
C     RHCO3=HCO3-CO3+H dissociation
C
      S0=AHY1+ACO31+DPHCO
      S1=AMAX1(0.0,S0**2-4.0*(AHY1*ACO31-DPHCO*AHCO31))
      RHCO3=TSLX*(S0-SQRT(S1))
C
C     RALO1=ALOH-AL+OH dissociation
C
      RALO1=TSLX*(AAL1*AOH1-DPAL1*AALO1)/(AOH1+DPAL1)
C
C     RALO2=ALOH2-ALOH+OH dissociation 
C
      RALO2=TSLX*(AALO1*AOH1-DPAL2*AALO2)/(AOH1+DPAL2)
C
C     RALO3=ALOH3-ALOH2+OH dissociation 
C
      RALO3=TSLX*(AALO2*AOH1-DPAL3*AALO3)/(AOH1+DPAL3)
C
C     RALO4=ALOH4-ALOH3+OH dissociation 
C
      RALO4=TSLX*(AALO3*AOH1-DPAL4*AALO4)/(AOH1+DPAL4)
C
C     RALS=ALSO4-AL+SO4 dissociation 
C
      S0=AAL1+ASO41+DPALS
      S1=AMAX1(0.0,S0**2-4.0*(AAL1*ASO41-DPALS*AALS1))
      RALS=TSLX*(S0-SQRT(S1))
C
C     RFEO1=FEOH-FE+OH dissociation 
C
      RFEO1=TSLX*(AFE1*AOH1-DPFE1*AFEO1)/(AOH1+DPFE1)
C
C     RFEO2=FEOH2-FEOH+OH dissociation 
C
      RFEO2=TSLX*(AFEO1*AOH1-DPFE2*AFEO2)/(AOH1+DPFE2)
C
C     RFEO3=FEOH3-FEOH2+OH dissociation 
C
      RFEO3=TSLX*(AFEO2*AOH1-DPFE3*AFEO3)/(AOH1+DPFE3)
C
C     RFE04=ALOH4-ALOH3+OH dissociation 
C
      RFEO4=TSLX*(AFEO3*AOH1-DPFE4*AFEO4)/(AOH1+DPFE4)
C
C     RFES-FE+SO4 dissociation 
C
      S0=AFE1+ASO41+DPFES
      S1=AMAX1(0.0,S0**2-4.0*(AFE1*ASO41-DPFES*AFES1))
      RFES=TSLX*(S0-SQRT(S1))
C
C     RCAO=CAOH-CA+OH dissociation 
C
      RCAO=TSLX*(ACA1*AOH1-DPCAO*ACAO1)/(AOH1+DPCAO)
C
C     RCAC=CACO3-CA+CO3 dissociation 
C
      S0=ACA1+ACO31+DPCAC
      S1=AMAX1(0.0,S0**2-4.0*(ACA1*ACO31-DPCAC*ACAC1))
      RCAC=TSLX*(S0-SQRT(S1))
C
C     RCAH=CAHCO3-CA+HCO3 dissociation 
C
      S0=ACA1+AHCO31+DPCAH
      S1=AMAX1(0.0,S0**2-4.0*(ACA1*AHCO31-DPCAH*ACAH1))
      RCAH=TSLX*(S0-SQRT(S1))
C
C     RCAS=CASO4-CA+SO4 dissociation 
C
      S0=ACA1+ASO41+DPCAS
      S1=AMAX1(0.0,S0**2-4.0*(ACA1*ASO41-DPCAS*ACAS1))
      RCAS=TSLX*(S0-SQRT(S1))
C
C     RMGO=MGOH-MG+OH dissociation 
C
      RMGO=TSLX*(AMG1*AOH1-DPMGO*AMGO1)/(AOH1+DPMGO)
C
C     RMGC=MGCO3-MG+CO3 dissociation 
C
      S0=AMG1+ACO31+DPMGC
      S1=AMAX1(0.0,S0**2-4.0*(AMG1*ACO31-DPMGC*AMGC1))
      RMGC=TSLX*(S0-SQRT(S1))
C
C     RMGH=MGHCO3-MG+HCO3 dissociation 
C
      S0=AMG1+AHCO31+DPMGH
      S1=AMAX1(0.0,S0**2-4.0*(AMG1*AHCO31-DPMGH*AMGH1))
      RMGH=TSLX*(S0-SQRT(S1))
C
C     RMGS=MGSO4-MG+SO4 dissociation 
C
      S0=AMG1+ASO41+DPMGS
      S1=AMAX1(0.0,S0**2-4.0*(AMG1*ASO41-DPMGS*AMGS1))
      RMGS=TSLX*(S0-SQRT(S1))
C
C     RNAC=NACO3-NA+CO3 dissociation 
C
      S0=ANA1+ACO31+DPNAC
      S1=AMAX1(0.0,S0**2-4.0*(ANA1*ACO31-DPNAC*ANAC1))
      RNAC=TSLX*(S0-SQRT(S1))
C
C     RNAS=NASO4-NA+SO4 dissociation 
C
      S0=ANA1+ASO41+DPNAS
      S1=AMAX1(0.0,S0**2-4.0*(ANA1*ASO41-DPNAS*ANAS1))
      RNAS=TSLX*(S0-SQRT(S1))
C
C     RKAS=KSO4-K+SO4 dissociation 
C
      S0=AKA1+ASO41+DPKAS
      S1=AMAX1(0.0,S0**2-4.0*(AKA1*ASO41-DPKAS*AKAS1))
      RKAS=TSLX*(S0-SQRT(S1))
C
C     PHOSPHORUS IN NON-BAND SOIL ZONE
C
      IF(VOLWPO.GT.ZEROS2(NY,NX))THEN
C
C     RH1P=HPO4-H+PO4 dissociation in non-band 
C
      RH1P=TSLX*(AH0P1*AHY1-DPH1P*AH1P1)/(DPH1P+AHY1)
C
C     RH2P=H2PO4-H+HPO4 dissociation in non-band 
C
      RH2P=TSLX*(AH1P1*AHY1-DPH2P*AH2P1)/(DPH2P+AHY1)
C     IF(NY.EQ.5.AND.L.EQ.10)THEN
C     WRITE(*,22)'RH2P',I,J,NX,NY,L,M,RH2P,TSLX,S0,S1,DP,DPH2P,A2
C    2,CH1P1,CHY1,CH2P1,H2PO4(L,NY,NX),VOLWPX,RH2PX,XH2PS(L,NY,NX)
C    3,TUPH2P(L,NY,NX)
22    FORMAT(A8,6I4,60E12.4)
C     ENDIF
C
C     RH3P=H3PO4-H+H2PO4 dissociation in non-band 
C
      RH3P=TSLX*(AH2P1*AHY1-DPH3P*AH3P1)/(DPH3P+AHY1) 
C
C     RF1P=FEHPO4-FE+HPO4 dissociation in non-band 
C
      S0=AFE1+AH1P1+DPF1P
      S1=AMAX1(0.0,S0**2-4.0*(AFE1*AH1P1-DPF1P*AF1P1))
      RF1P=TSLX*(S0-SQRT(S1))
C
C     RF2P=FEH2PO4-FE+H2PO4 dissociation in non-band 
C
      S0=AFE1+AH2P1+DPF2P
      S1=AMAX1(0.0,S0**2-4.0*(AFE1*AH2P1-DPF2P*AF2P1))
      RF2P=TSLX*(S0-SQRT(S1))
C
C     RC0P=CAPO4-CA+PO4 dissociation in non-band 
C
      S0=ACA1+AH0P1+DPC0P
      S1=AMAX1(0.0,S0**2-4.0*(ACA1*AH0P1-DPC0P*AC0P1))
      RC0P=TSLX*(S0-SQRT(S1))
C
C     RC1P=CAHPO4-CA+HPO4 dissociation in non-band 
C
      S0=ACA1+AH1P1+DPC1P
      S1=AMAX1(0.0,S0**2-4.0*(ACA1*AH1P1-DPC1P*AC1P1))
      RC1P=TSLX*(S0-SQRT(S1))
C
C     RC2P=CAH2PO4-CA+H2PO4 dissociation in non-band 
C
      S0=ACA1+AH2P1+DPC2P
      S1=AMAX1(0.0,S0**2-4.0*(ACA1*AH2P1-DPC2P*AC2P1))
      RC2P=TSLX*(S0-SQRT(S1))
C
C     RM1P=MGHPO4-MG+HPO4 dissociation in non-band 
C
      S0=AMG1+AH1P1+DPM1P
      S1=AMAX1(0.0,S0**2-4.0*(AMG1*AH1P1-DPM1P*AM1P1))
      RM1P=TSLX*(S0-SQRT(S1))
      ELSE
      RH1P=0.0
      RH2P=0.0
      RH3P=0.0
      RF1P=0.0
      RF2P=0.0
      RC0P=0.0
      RC1P=0.0
      RC2P=0.0
      RM1P=0.0
      ENDIF
C
C     PHOSPHORUS IN BAND SOIL ZONE
C
      IF(VOLWPB.GT.ZEROS2(NY,NX))THEN
C
C     RH1B=HPO4-H+PO4 dissociation in band
C
      RH1B=TSLX*(AH0PB*AHY1-DPH1P*AH1PB)/(AHY1+DPH1P)
C
C     RH2B=H2PO4-H+HPO4 dissociation in band 
C
      RH2B=TSLX*(AH1PB*AHY1-DPH2P*AH2PB)/(AHY1+DPH2P)
C
C     RH3B=H3PO4-H+H2PO4 dissociation in band 
C
      RH3B=TSLX*(AH2PB*AHY1-DPH3P*AH3PB)/(AHY1+DPH3P)
C
C     RF1B=FEHPO4-FE+HPO4 dissociation in band 
C
      S0=AFE1+AH1PB+DPF1P
      S1=AMAX1(0.0,S0**2-4.0*(AFE1*AH1PB-DPF1P*AF1PB))
      RF1B=TSLX*(S0-SQRT(S1))
C
C     RF2B=FEH2PO4-FE+H2PO4 dissociation in band 
C
      S0=AFE1+AH2PB+DPF2P
      S1=AMAX1(0.0,S0**2-4.0*(AFE1*AH2PB-DPF2P*AF2PB))
      RF2B=TSLX*(S0-SQRT(S1))
C
C     RC0B=CAPO4-CA+PO4 dissociation in band 
C
      S0=ACA1+AH0PB+DPC0P
      S1=AMAX1(0.0,S0**2-4.0*(ACA1*AH0PB-DPC0P*AC0PB))
      RC0B=TSLX*(S0-SQRT(S1))
C
C     RC1B=CAHPO4-CA+HPO4 dissociation in band 
C
      S0=ACA1+AH1PB+DPC1P
      S1=AMAX1(0.0,S0**2-4.0*(ACA1*AH1PB-DPC1P*AC1PB))
      RC1B=TSLX*(S0-SQRT(S1))
C
C     RC2B=CAH2PO4-CA+H2PO4 dissociation in band 
C
      S0=ACA1+AH2PB+DPC2P
      S1=AMAX1(0.0,S0**2-4.0*(ACA1*AH2PB-DPC2P*AC2PB))
      RC2B=TSLX*(S0-SQRT(S1))
C
C     RM1B=MGHPO4-MG+HPO4 dissociation in band 
C
      S0=AMG1+AH1PB+DPM1P
      S1=AMAX1(0.0,S0**2-4.0*(AMG1*AH1PB-DPM1P*AM1PB))
      RM1B=TSLX*(S0-SQRT(S1))
      ELSE
      RH1B=0.0
      RH2B=0.0
      RH3B=0.0
      RF1B=0.0
      RF2B=0.0
      RC0B=0.0
      RC1B=0.0
      RC2B=0.0
      RM1B=0.0
      ENDIF
      
C
C     TOTAL ION FLUXES FOR CURRENT ITERATION
C     FROM ALL REACTIONS ABOVE
C
C     RN4S,RN4B=net NH4 flux in non-band,band
C     RN3S,RN3B=net NH3 flux in non-band,band
C     RAL,RFE,RHY,RCA,RMG,RNA,RKA,ROH=net Al,Fe,H,Ca,Mg,Na,K,OH flux
C     RSO4,RCO3,RHCO,RCO2=net SO4,CO3,HCO3,CO2 flux
C     RAL1,RAL2,RAL3,RAL4,RALS=net AlOH,AlOH2,AlOH3,AlOH4,AlSO4
C     RFE1,RFE2,RFE3,RFE4,RFES=net FeOH,FeOH2,FeOH3,FeOH4,FeSO4
C     RHP0,RHP1,RHP2,RHP3=net PO4,HPO4,H2PO4,H3PO4 flux in non-band
C     RXH0,RXH1,RXH2,RX1P,RX2P=net R-O,R-OH,R-OH2,R-HPO4,R-H2PO4 in non-band
C     RHB0,RHB1,RHB2,RHB3=net PO4,HPO4,H2PO4,H3PO4 flux in band
C     RBH0,RBH1,RBH2,RB1P,RB2P=net R-O,R-OH,R-OH2,R-HPO4,R-H2PO4 in band
C
      RN4S=RNH4-RXN4
      RN4B=RNHB-RXNB
      RN3S=-RNH4
      RN3B=-RNHB
      RAL=-RHAL1-RXAL-RALO1-RALS
     2-(RHA0P1+RHA0P2)*VLPO4(L,NY,NX)
     3-(RHA0B1+RHA0B2)*VLPOB(L,NY,NX)
      RFE=-RHFE1-RXFE-RFEO1-RFES
     2-(RHF0P1+RHF0P2+RF1P+RF2P)*VLPO4(L,NY,NX)
     2-(RHF0B1+RHF0B2+RF1B+RF2B)*VLPOB(L,NY,NX)
      RHY=-RNH4*VLNH4(L,NY,NX)-RNHB*VLNHB(L,NY,NX)
     2-RXHY-RXHC+2.0*(RHALO1+RHFEO1+RHCACO
     2+(RHA0P2+RHF0P2-RHA3P1-RHA4P2-RHF3P1-RHF4P2)*VLPO4(L,NY,NX)
     3+(RHA0B2+RHF0B2-RHA3B1-RHA4B2-RHF3B1-RHF4B2)*VLPOB(L,NY,NX))
     4+3.0*(RHAL1+RHFE1
     5-(RHA4P1+RHF4P1)*VLPO4(L,NY,NX)
     6-(RHF4B1+RHA4B1)*VLPOB(L,NY,NX))
     7+4.0*(RHCAH1*VLPO4(L,NY,NX)+RHCHB1*VLPOB(L,NY,NX))
     8+7.0*(RHCAH2*VLPO4(L,NY,NX)+RHCHB2*VLPOB(L,NY,NX))
     9+RHALO2+RHFEO2-RHALO4-RHFEO4+RHCACH-RCO2Q-RHCO3
     1+(RHA0P1-RHA2P1+RHA1P2-RHA3P2+RHF0P1-RHF2P1+RHF1P2-RHF3P2
     2+RHCAD2-RXOH2-RXOH1-RH1P-RH2P-RH3P)*VLPO4(L,NY,NX)
     3+(RHA0B1-RHA2B1+RHA1B2-RHA3B2+RHF0B1-RHF2B1+RHF1B2-RHF3B2
     4+RHCDB2-RXO2B-RXO1B-RH1B-RH2B-RH3B)*VLPOB(L,NY,NX)
      RCA=-RPCACX-RPCASO-RXCA-RCAO-RCAC-RCAH-RCAS
     2-(RPCADX+RPCAMX+RC0P+RC1P+RC2P)*VLPO4(L,NY,NX)
     3-(RPCDBX+RPCMBX+RC0B+RC1B+RC2B)*VLPOB(L,NY,NX)
     4-5.0*(RPCAHX*VLPO4(L,NY,NX)+RPCHBX*VLPOB(L,NY,NX))
      RMG=-RXMG-RMGO-RMGC-RMGH-RMGS
     2-RM1P*VLPO4(L,NY,NX)-RM1B*VLPOB(L,NY,NX)
      RNA=-RXNA-RNAC-RNAS
      RKA=-RXKA-RKAS
      ROH=-RCAO-RMGO-RALO1
     2-RALO2-RALO3-RALO4-RFEO1-RFEO2-RFEO3-RFEO4
     3-(-RYH2P-RXH1P)*VLPO4(L,NY,NX)
     4-(-RYH2B-RXH1B)*VLPOB(L,NY,NX)
      RSO4=-RPCASO-RALS-RFES-RCAS-RMGS-RNAS-RKAS
      RCO3=-RHCAC3-RHCO3-RCAC-RMGC-RNAC
      RHCO=-RHCACH-RCO2Q-RCAH-RMGH+RHCO3
      RCO2=-RHCACO+RCO2Q
      RAL1=-RHALO1+RALO1-RALO2
     2-(RHA1P1+RHA1P2)*VLPO4(L,NY,NX)
     3-(RHA1B1+RHA1B2)*VLPOB(L,NY,NX)
      RAL2=-RHALO2+RALO2-RALO3
     2-(RHA2P1+RHA2P2)*VLPO4(L,NY,NX)
     3-(RHA2B1+RHA2B2)*VLPOB(L,NY,NX) 
      RAL3=-RHALO3+RALO3-RALO4
     2-(RHA3P1+RHA3P2)*VLPO4(L,NY,NX)
     3-(RHA3B1+RHA3B2)*VLPOB(L,NY,NX)
      RAL4=-RHALO4+RALO4
     2-(RHA4P1+RHA4P2)*VLPO4(L,NY,NX)
     3-(RHA4B1+RHA4B2)*VLPOB(L,NY,NX)
      RFE1=-RHFEO1+RFEO1-RFEO2
     2-(RHF1P1+RHF1P2)*VLPO4(L,NY,NX)
     3-(RHF1B1+RHF1B2)*VLPOB(L,NY,NX)
      RFE2=-RHFEO2+RFEO2-RFEO3
     2-(RHF2P1+RHF2P2)*VLPO4(L,NY,NX)
     3-(RHF2B1+RHF2B2)*VLPOB(L,NY,NX)
      RFE3=-RHFEO3+RFEO3-RFEO4
     2-(RHF3P1+RHF3P2)*VLPO4(L,NY,NX)
     3-(RHF3B1+RHF3B2)*VLPOB(L,NY,NX)
      RFE4=-RHFEO4+RFEO4
     2-(RHF4P1+RHF4P2)*VLPO4(L,NY,NX)
     3-(RHF4B1+RHF4B2)*VLPOB(L,NY,NX)
      RHP0=-RH1P-RC0P
      RHP1=-RHA0P1-RHA1P1-RHA2P1-RHA3P1
     2-RHA4P1-RHF0P1-RHF1P1-RHF2P1 
     3-RHF3P1-RHF4P1-RPCAD1-3.0*RHCAH1-RXH1P 
     4+RH1P-RH2P-RF1P-RC1P-RM1P
      RHP2=-RHA0P2-RHA1P2-RHA2P2-RHA3P2
     2-RHA4P2-RHF0P2-RHF1P2-RHF2P2 
     3-RHF3P2-RHF4P2-RHCAD2-3.0*RHCAH2
     4-2.0*RPCAMX-RXH2P-RYH2P+RH2P-RH3P-RF2P-RC2P
      RHP3=RH3P
      RXH0=-RXOH1
      RXH1=RXOH1-RXOH2-RYH2P-RXH1P 
      RXH2=RXOH2-RXH2P
      RX1P=RXH1P 
      RX2P=RXH2P+RYH2P 
      RHB0=-RH1B-RC0B
      RHB1=-RHA0B1-RHA1B1-RHA2B1-RHA3B1
     2-RHA4B1-RHF0B1-RHF1B1-RHF2B1 
     3-RHF3B1-RHF4B1-RPCDB1-3.0*RHCHB1-RXH1B 
     4+RH1B-RH2B-RF1B-RC1B-RM1B
      RHB2=-RHA0B2-RHA1B2-RHA2B2-RHA3B2
     2-RHA4B2-RHF0B2-RHF1B2-RHF2B2 
     3-RHF3B2-RHF4B2-RHCDB2-3.0*RHCHB2
     4-2.0*RPCMBX-RXH2B-RYH2B+RH2B-RH3B-RF2B-RC2B
      RHB3=RH3B
      RBH0=-RXO1B
      RBH1=RXO1B-RXO2B-RYH2B-RXH1B 
      RBH2=RXO2B-RXH2B
      RB1P=RXH1B 
      RB2P=RXH2B+RYH2B 
C     IF(NY.EQ.5.AND.L.EQ.10)THEN
C     WRITE(20,23)'RHP1',I,J,NX,NY,L,M,RHP1,RHA0P1 
C    2,RHA1P1,RHA2P1,RHA3P1,RHA4P1 
C    3,RHF0P1,RHF1P1,RHF2P1,RHF3P1 
C    4,RHF4P1,RPCAD1,3.0*( RHCAH1),RXH1P,RH1P,RH2P,RF1P,RC1P,RM1P
C     WRITE(20,23)'RHP2',I,J,NX,NY,L,M,RHP2,RHA0P2,RHA1P2
C    2,RHA2P2,RHA3P2,RHA4P2,RHF0P2 
C    3,RHF1P2,RHF2P2,RHF3P2,RHF4P2,RHCAD2 
C    4,RHCAH2,RPCAMX,RXH2P,RYH2P,RH2P,RH3P,RF2P,RC2P
23    FORMAT(A8,6I4,60E12.4)
C     ENDIF
C
C     UPDATE ION CONCENTRATIONS FOR CURRENT ITERATION
C     FROM TOTAL ION FLUXES
C
      CN41=CN41+RN4S
      CN4B=CN4B+RN4B
      CN31=CN31+RN3S
      CN3B=CN3B+RN3B
      CAL1=CAL1+RAL
      CFE1=CFE1+RFE
      CHY1=CHY1+RHY
      CCA1=CCA1+RCA
      CMG1=CMG1+RMG
      CNA1=CNA1+RNA
      CKA1=CKA1+RKA
      COH1=COH1+ROH
      CSO41=CSO41+RSO4
      CCO31=CCO31+RCO3
      CHCO31=CHCO31+RHCO
      CCO21=CCO21+RCO2
      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
      CH0PB=CH0PB+RHB0
      CH1PB=CH1PB+RHB1
      CH2PB=CH2PB+RHB2
      CH3PB=CH3PB+RHB3
      CF1PB=CF1PB+RF1B
      CF2PB=CF2PB+RF2B
      CC0PB=CC0PB+RC0B
      CC1PB=CC1PB+RC1B
      CC2PB=CC2PB+RC2B
      CM1PB=CM1PB+RM1B
C
C     RHHY,RHOH=H2O-H+OH equilibration
C
      CHY2=10.0**(-PH(L,NY,NX))*1.0E+03
      COH2=DPH2O/CHY2
      RHHY=CHY2-CHY1
      RHOH=COH2-COH1
      CHY1=CHY1+RHHY
      COH1=COH1+RHOH
C     IF(I.EQ.256.AND.L.EQ.1)THEN 
C     WRITE(*,1111)'CO3',I,J,L,M,CCO31,CHCO31,CCO21,DPHCO,DPCO2
C    2,A1,A2,RCO3,RHCAC3,RHCO3,RCAC,RMGC,RNAC
C    3,RHCO,RHCACH,RCO2Q,RCAH,RMGH,RHCO3
C    4,RCO2,RHCACO,AHY1,AHCO31,ACO21,TRCO2(L,NY,NX) 
C     WRITE(*,1111)'CCA1',I,J,L,M,CCA1,ACA1,AHY1,AH1P1,AH2P1,ACO31 
C    2,XCA1,AHCO31,RCA,RPCACX,RPCASO,RXCA,RCAO,RCAC,RCAH,RCAS
C    2,(RPCADX+RPCAMX+RC0P+RC1P+RC2P)*VLPO4(L,NY,NX)
C    3,(RPCDBX+RPCMBX+RC0B+RC1B+RC2B)*VLPOB(L,NY,NX)
C    4,5.0*(RPCAHX*VLPO4(L,NY,NX)+RPCHBX*VLPOB(L,NY,NX))
C     WRITE(*,1111)'CAL1',I,J,L,M,CAL1,A3,AAL1
C    2,RAL,RHAL1,RXAL,RALO1,RALS
C    2,RHA0P1,RHA0P2,VLPO4(L,NY,NX)
C    3,RHA0B1,RHA0B2,VLPOB(L,NY,NX)
C     WRITE(*,1111)'CFEO2',I,J,L,M,CFEO2,CFEO2*A1
C    2,RFE2,RHFEO2,RHF2P1,RHF2P2,RHF2B1
C    2,RHF2B2,RFEO2,RFEO3
C     WRITE(20,1111)'CHY1',I,J,L,M,CHY1,RHY,RHHY,CN31,CN41,RNH4 
C    2,RXHY,RXHC
C    2,2.0*(RHALO1+RHFEO1+RHCACO
C    2+(RHA0P2+RHF0P2-RHA3P1-RHA4P2-RHF3P1-RHF4P2)*VLPO4(L,NY,NX)
C    3+(RHA0B2+RHF0B2-RHA3B1-RHA4B2-RHF3B1-RHF4B2)*VLPOB(L,NY,NX))
C    4,3.0*(RHAL1+RHFE1
C    5-(RHA4P1+RHF4P1)*VLPO4(L,NY,NX)
C    6-(RHF4B1+RHA4B1)*VLPOB(L,NY,NX))
C    7,4.0*(RHCAH1*VLPO4(L,NY,NX)+RHCHB1*VLPOB(L,NY,NX))
C    8,7.0*(RHCAH2*VLPO4(L,NY,NX)+RHCHB2*VLPOB(L,NY,NX))
C    9,RHALO2,RHFEO2,RHALO4,RHFEO4,RHCACH,RCO2Q,RHCO3
C    1,(RHA0P1-RHA2P1+RHA1P2-RHA3P2+RHF0P1-RHF2P1+RHF1P2-RHF3P2
C    2+RHCAD2-RXOH2-RXOH1-RH1P-RH2P-RH3P)*VLPO4(L,NY,NX)
C    3,(RHA0B1-RHA2B1+RHA1B2-RHA3B2+RHF0B1-RHF2B1+RHF1B2-RHF3B2
C    4+RHCDB2-RXO2B-RXO1B-RH1B-RH2B-RH3B)*VLPOB(L,NY,NX)
C     WRITE(20,1111)'COH1',I,J,L,M,COH1,ROH,RHOH 
C    2,RYH2P,RYH2B,RXH1P,RXH1B,RPALPX,RCAO,RMGO
C    2,RPCAHX,RALO1,RALO2,RALO3,RALO4,RFEO1,RFEO2,RFEO3,RFEO4
1111  FORMAT(A8,4I4,80E12.4)
C     ENDIF
C
C     UPDATE EXCHANGEABLE ION CONCENTRATIONS IN CURRENT
C     ITERATION FROM TOTAL ION FLUXES
C
      XN41=XN41+RXN4
      XN4B=XN4B+RXNB
      XHY1=XHY1+RXHY
      XAL1=XAL1+RXAL
      XFE1=XFE1+RXFE
      XCA1=XCA1+RXCA
      XMG1=XMG1+RXMG
      XNA1=XNA1+RXNA
      XKA1=XKA1+RXKA
      XHC1=XHC1+RXHC
      XALO21=XALO21+RXALO2
      XFEO21=XFEO21+RXFEO2
      XOH01=XOH01+RXH0
      XOH11=XOH11+RXH1
      XOH21=XOH21+RXH2
      XH1P1=XH1P1+RX1P
      XH2P1=XH2P1+RX2P
      XH01B=XH01B+RBH0
      XH11B=XH11B+RBH1
      XH21B=XH21B+RBH2
      X1P1B=X1P1B+RB1P
      X2P1B=X2P1B+RB2P
C
C     UPDATE PRECIPITATE CONCENTRATIONS IN CURRENT
C     ITERATION FROM TOTAL ION FLUXES
C
      PALOH1=PALOH1+RPALOX
      PFEOH1=PFEOH1+RPFEOX
      PCACO1=PCACO1+RPCACX
      PCASO1=PCASO1+RPCASO
      PALPO1=PALPO1+RPALPX
      PFEPO1=PFEPO1+RPFEPX
      PCAPD1=PCAPD1+RPCADX
      PCAPH1=PCAPH1+RPCAHX
      PCAPM1=PCAPM1+RPCAMX
      PALPOB=PALPOB+RPALBX
      PFEPOB=PFEPOB+RPFEBX
      PCAPDB=PCAPDB+RPCDBX
      PCAPHB=PCAPHB+RPCHBX
      PCAPMB=PCAPMB+RPCMBX
C
C     ACCUMULATE TOTAL ION FLUXES FOR ALL ITERATIONS
C
C     TRN4S,TRN4B=total NH4 flux in non-band,band
C     TRN3S,TRN3B=total NH3 flux in non-band,band
C     TRAL,TRFE,TRHY,TRCA,TRMG,TRNA,TRKA,TROH=totalAl,Fe,H,Ca,Mg,Na,K,OH flux
C     TRSO4,TRCO3,TRHCO,TRCO2=total SO4,CO3,HCO3,CO2 flux
C     TRAL1,TRAL2,TRAL3,TRAL4,TRALS=total AlOH,AlOH2,AlOH3,AlOH4,AlSO4
C     TRFE1,TRFE2,TRFE3,TRFE4,TRFES=total FeOH,FeOH2,FeOH3,FeOH4,FeSO4
C     TRCAO,TRCAC,TRCAH,TRCAS=total CaOH,CaCO3,CaHCO3,CaSO4 flux
C     TRMGO,TRMGC,TRMGH,TRMGS=total MgOH,MgCO3,MgHCO3,MgSO4 flux
C     TRNAC,TRNAS,TRKAS=total NaCO3,NaSO4,KSO4 flux
C     TRH0P,TRH1P,TRH2P,TRH3P=net PO4,HPO4,H2PO4,H3PO4 flux in non-band
C     TRF1P,TRF2P,TRC0P,TRC1P,TRC2P,TRM1P
C     =total FeHPO4,FeH2PO4,CaPO4,CaHPO4,CaH2PO4,MgHPO4 in non-band
C     TRH0B,TRH1B,TRH2B,TRH3B=net PO4,HPO4,H2PO4,H3PO4 flux in band
C     TRF1B,TRF2B,TRC0B,TRC1B,TRC2B,TRM1B
C     =total FeHPO4,FeH2PO4,CaPO4,CaHPO4,CaH2PO4,MgHPO4 in band
C     TRNX4,TRNXB=total NH4 adsorption in non-band,band
C     TRXHY,TRXAL,TRXFE,TRXCA,TRXMG,TRXNA,TRXKA=total H,Al,Fe,Ca,Mg,Na,K adsorpn
C     TRXHC,TRXAL2,TRXFE2=total HCO3,AlOH2,FeOH2 adsorption
C     TRXH0,TRXH1,TRXH2,TRX1P,TRX2P
C     =total R-O,R-OH,R-OH2,R-HPO4,R-H2PO4 adsorption in non-band
C     TRBH0,TRBH1,TRBH2,TRB1P,TRB2P
C     =total R-O,R-OH,R-OH2,R-HPO4,R-H2PO4 adsorption in band
C     TRALOH,TRFEOH,TRCACO,TRCASO=total AlOH3,FeOH3,CaCO3,CaSO4 precipitation
C     TRALPO,TRFEPO,TRCAPD,TRCAPH,TRCAPM
C     =total AlPO4,FePO4,CaHPO4,apatite,Ca(H2PO4)2 precipitation in non-band
C     TRALPB,TRFEPB,TRCPDB,TRCPHB,TRCPMB
C     =total AlPO4,FePO4,CaHPO4,apatite,Ca(H2PO4)2 precipitation in band
C
      TRN4S(L,NY,NX)=TRN4S(L,NY,NX)+RN4S
      TRN4B(L,NY,NX)=TRN4B(L,NY,NX)+RN4B
      TRN3S(L,NY,NX)=TRN3S(L,NY,NX)+RN3S
      TRN3B(L,NY,NX)=TRN3B(L,NY,NX)+RN3B
      TRAL(L,NY,NX)=TRAL(L,NY,NX)+RAL
      TRFE(L,NY,NX)=TRFE(L,NY,NX)+RFE
      TRHY(L,NY,NX)=TRHY(L,NY,NX)+RHY+RHHY
      TRCA(L,NY,NX)=TRCA(L,NY,NX)+RCA
      TRMG(L,NY,NX)=TRMG(L,NY,NX)+RMG
      TRNA(L,NY,NX)=TRNA(L,NY,NX)+RNA
      TRKA(L,NY,NX)=TRKA(L,NY,NX)+RKA
      TROH(L,NY,NX)=TROH(L,NY,NX)+ROH+RHOH
      TRSO4(L,NY,NX)=TRSO4(L,NY,NX)+RSO4
      TRCO3(L,NY,NX)=TRCO3(L,NY,NX)+RCO3
      TRHCO(L,NY,NX)=TRHCO(L,NY,NX)+RHCO
      TRCO2(L,NY,NX)=TRCO2(L,NY,NX)+RCO2
      TRAL1(L,NY,NX)=TRAL1(L,NY,NX)+RAL1
      TRAL2(L,NY,NX)=TRAL2(L,NY,NX)+RAL2 
      TRAL3(L,NY,NX)=TRAL3(L,NY,NX)+RAL3
      TRAL4(L,NY,NX)=TRAL4(L,NY,NX)+RAL4
      TRALS(L,NY,NX)=TRALS(L,NY,NX)+RALS
      TRFE1(L,NY,NX)=TRFE1(L,NY,NX)+RFE1
      TRFE2(L,NY,NX)=TRFE2(L,NY,NX)+RFE2
      TRFE3(L,NY,NX)=TRFE3(L,NY,NX)+RFE3
      TRFE4(L,NY,NX)=TRFE4(L,NY,NX)+RFE4
      TRFES(L,NY,NX)=TRFES(L,NY,NX)+RFES
      TRCAO(L,NY,NX)=TRCAO(L,NY,NX)+RCAO
      TRCAC(L,NY,NX)=TRCAC(L,NY,NX)+RCAC
      TRCAH(L,NY,NX)=TRCAH(L,NY,NX)+RCAH
      TRCAS(L,NY,NX)=TRCAS(L,NY,NX)+RCAS
      TRMGO(L,NY,NX)=TRMGO(L,NY,NX)+RMGO
      TRMGC(L,NY,NX)=TRMGC(L,NY,NX)+RMGC
      TRMGH(L,NY,NX)=TRMGH(L,NY,NX)+RMGH
      TRMGS(L,NY,NX)=TRMGS(L,NY,NX)+RMGS
      TRNAC(L,NY,NX)=TRNAC(L,NY,NX)+RNAC
      TRNAS(L,NY,NX)=TRNAS(L,NY,NX)+RNAS
      TRKAS(L,NY,NX)=TRKAS(L,NY,NX)+RKAS
      TRH0P(L,NY,NX)=TRH0P(L,NY,NX)+RHP0
      TRH1P(L,NY,NX)=TRH1P(L,NY,NX)+RHP1
      TRH2P(L,NY,NX)=TRH2P(L,NY,NX)+RHP2
      TRH3P(L,NY,NX)=TRH3P(L,NY,NX)+RHP3
      TRF1P(L,NY,NX)=TRF1P(L,NY,NX)+RF1P
      TRF2P(L,NY,NX)=TRF2P(L,NY,NX)+RF2P
      TRC0P(L,NY,NX)=TRC0P(L,NY,NX)+RC0P
      TRC1P(L,NY,NX)=TRC1P(L,NY,NX)+RC1P
      TRC2P(L,NY,NX)=TRC2P(L,NY,NX)+RC2P
      TRM1P(L,NY,NX)=TRM1P(L,NY,NX)+RM1P
      TRH0B(L,NY,NX)=TRH0B(L,NY,NX)+RHB0
      TRH1B(L,NY,NX)=TRH1B(L,NY,NX)+RHB1
      TRH2B(L,NY,NX)=TRH2B(L,NY,NX)+RHB2
      TRH3B(L,NY,NX)=TRH3B(L,NY,NX)+RHB3
      TRF1B(L,NY,NX)=TRF1B(L,NY,NX)+RF1B
      TRF2B(L,NY,NX)=TRF2B(L,NY,NX)+RF2B
      TRC0B(L,NY,NX)=TRC0B(L,NY,NX)+RC0B
      TRC1B(L,NY,NX)=TRC1B(L,NY,NX)+RC1B
      TRC2B(L,NY,NX)=TRC2B(L,NY,NX)+RC2B
      TRM1B(L,NY,NX)=TRM1B(L,NY,NX)+RM1B
      TRXN4(L,NY,NX)=TRXN4(L,NY,NX)+RXN4
      TRXNB(L,NY,NX)=TRXNB(L,NY,NX)+RXNB
      TRXHY(L,NY,NX)=TRXHY(L,NY,NX)+RXHY
      TRXAL(L,NY,NX)=TRXAL(L,NY,NX)+RXAL
      TRXFE(L,NY,NX)=TRXFE(L,NY,NX)+RXFE
      TRXCA(L,NY,NX)=TRXCA(L,NY,NX)+RXCA
      TRXMG(L,NY,NX)=TRXMG(L,NY,NX)+RXMG
      TRXNA(L,NY,NX)=TRXNA(L,NY,NX)+RXNA
      TRXKA(L,NY,NX)=TRXKA(L,NY,NX)+RXKA
      TRXHC(L,NY,NX)=TRXHC(L,NY,NX)+RXHC
      TRXAL2(L,NY,NX)=TRXAL2(L,NY,NX)+RXALO2
      TRXFE2(L,NY,NX)=TRXFE2(L,NY,NX)+RXFEO2
      TRXH0(L,NY,NX)=TRXH0(L,NY,NX)+RXH0
      TRXH1(L,NY,NX)=TRXH1(L,NY,NX)+RXH1
      TRXH2(L,NY,NX)=TRXH2(L,NY,NX)+RXH2
      TRX1P(L,NY,NX)=TRX1P(L,NY,NX)+RX1P
      TRX2P(L,NY,NX)=TRX2P(L,NY,NX)+RX2P
      TRBH0(L,NY,NX)=TRBH0(L,NY,NX)+RBH0
      TRBH1(L,NY,NX)=TRBH1(L,NY,NX)+RBH1
      TRBH2(L,NY,NX)=TRBH2(L,NY,NX)+RBH2
      TRB1P(L,NY,NX)=TRB1P(L,NY,NX)+RB1P
      TRB2P(L,NY,NX)=TRB2P(L,NY,NX)+RB2P
      TRALOH(L,NY,NX)=TRALOH(L,NY,NX)+RPALOX
      TRFEOH(L,NY,NX)=TRFEOH(L,NY,NX)+RPFEOX
      TRCACO(L,NY,NX)=TRCACO(L,NY,NX)+RPCACX
      TRCASO(L,NY,NX)=TRCASO(L,NY,NX)+RPCASO
      TRALPO(L,NY,NX)=TRALPO(L,NY,NX)+RPALPX
      TRFEPO(L,NY,NX)=TRFEPO(L,NY,NX)+RPFEPX
      TRCAPD(L,NY,NX)=TRCAPD(L,NY,NX)+RPCADX
      TRCAPH(L,NY,NX)=TRCAPH(L,NY,NX)+RPCAHX
      TRCAPM(L,NY,NX)=TRCAPM(L,NY,NX)+RPCAMX
      TRALPB(L,NY,NX)=TRALPB(L,NY,NX)+RPALBX
      TRFEPB(L,NY,NX)=TRFEPB(L,NY,NX)+RPFEBX
      TRCPDB(L,NY,NX)=TRCPDB(L,NY,NX)+RPCDBX
      TRCPHB(L,NY,NX)=TRCPHB(L,NY,NX)+RPCHBX
      TRCPMB(L,NY,NX)=TRCPMB(L,NY,NX)+RPCMBX
C
C     GO TO NEXT ITERATION
C
1000  CONTINUE
C
C     CONVERGENCE ITERATIONS COMPLETED
C
C     IF(J.EQ.24)THEN
C     WRITE(*,1119)'GAPON',I,J,L,M,CH0P1,CAL1,CFE1,CH0P1*A3*CAL1*A3
C    2,SPALP,CH0P1*A3*CFE1*A3,SPFEP
C    6,SPOH2,XOH11*CHY1*A1/XOH21,SPOH1,XOH01*CHY1*A1/XOH11
C    7,SPH2P,XOH21*CH2P1*A1/XH2P1,SXH2P,XOH11*CH2P1/(XH2P1*COH1)
C    8,SPH1P,XOH11*CH1P1*A2/(XH1P1*COH1*A1)
C    9,COH1*A1,CHY1*A1
1119  FORMAT(A8,4I4,24E11.3)
C     WRITE(*,1119)'CATION',I,J,L,M,CCEC,XN41+XHY1+3*XAL1+2*(XCA1+XMG1)
C    2+XNA1+XKA1,XN41,XHY1,XAL1,XCA1,XMG1,XNA1,XKA1,CN41,CHY1,CAL1,CCA1
C    2,CMG1,CNA1,CKA1,(CCA1*A2)**0.5*XN41/(CN41*A1*XCA1*2)
C    3,(CCA1*A2)**0.5*XHY1/(CHY1*A1*XCA1*2)
C    2,(CCA1*A2)**0.5*XAL1*3/((CAL1*A3)**0.333*XCA1*2)
C    3,(CCA1*A2)**0.5*XMG1*2/((CMG1*A2)**0.5*XCA1*2)
C    3,(CCA1*A2)**0.5*XNA1/(CNA1*A1*XCA1*2)
C    5,(CCA1*A2)**0.5*XKA1/(CKA1*A1*XCA1*2)
C    6,CHY1*A1*XCOO/XHC1,CALO2*A1*XCOO/XALO21
C     ENDIF
C
C     CONVERT TOTAL ION FLUXES FROM CHANGES IN CONCENTRATION
C     TO CHANGES IN MASS PER UNIT AREA FOR USE IN 'REDIST'
C
      TRN4S(L,NY,NX)=TRN4S(L,NY,NX)*VOLWNH
      TRN4B(L,NY,NX)=TRN4B(L,NY,NX)*VOLWNB
      TRN3S(L,NY,NX)=TRN3S(L,NY,NX)*VOLWNH
      TRN3B(L,NY,NX)=TRN3B(L,NY,NX)*VOLWNB
      TRAL(L,NY,NX)=TRAL(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRFE(L,NY,NX)=TRFE(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRHY(L,NY,NX)=TRHY(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRCA(L,NY,NX)=TRCA(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRMG(L,NY,NX)=TRMG(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRNA(L,NY,NX)=TRNA(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRKA(L,NY,NX)=TRKA(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TROH(L,NY,NX)=TROH(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRSO4(L,NY,NX)=TRSO4(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRCO3(L,NY,NX)=TRCO3(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRHCO(L,NY,NX)=TRHCO(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRCO2(L,NY,NX)=TRCO2(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRAL1(L,NY,NX)=TRAL1(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRAL2(L,NY,NX)=TRAL2(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRAL3(L,NY,NX)=TRAL3(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRAL4(L,NY,NX)=TRAL4(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRALS(L,NY,NX)=TRALS(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRFE1(L,NY,NX)=TRFE1(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRFE2(L,NY,NX)=TRFE2(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRFE3(L,NY,NX)=TRFE3(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRFE4(L,NY,NX)=TRFE4(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRFES(L,NY,NX)=TRFES(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRCAO(L,NY,NX)=TRCAO(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRCAC(L,NY,NX)=TRCAC(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRCAH(L,NY,NX)=TRCAH(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRCAS(L,NY,NX)=TRCAS(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRMGO(L,NY,NX)=TRMGO(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRMGC(L,NY,NX)=TRMGC(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRMGH(L,NY,NX)=TRMGH(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRMGS(L,NY,NX)=TRMGS(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRNAC(L,NY,NX)=TRNAC(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRNAS(L,NY,NX)=TRNAS(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRKAS(L,NY,NX)=TRKAS(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRH0P(L,NY,NX)=TRH0P(L,NY,NX)*VOLWPO
      TRH1P(L,NY,NX)=TRH1P(L,NY,NX)*VOLWPO
      TRH2P(L,NY,NX)=TRH2P(L,NY,NX)*VOLWPO
      TRH3P(L,NY,NX)=TRH3P(L,NY,NX)*VOLWPO
      TRF1P(L,NY,NX)=TRF1P(L,NY,NX)*VOLWPO
      TRF2P(L,NY,NX)=TRF2P(L,NY,NX)*VOLWPO
      TRC0P(L,NY,NX)=TRC0P(L,NY,NX)*VOLWPO
      TRC1P(L,NY,NX)=TRC1P(L,NY,NX)*VOLWPO
      TRC2P(L,NY,NX)=TRC2P(L,NY,NX)*VOLWPO
      TRM1P(L,NY,NX)=TRM1P(L,NY,NX)*VOLWPO
      TRH0B(L,NY,NX)=TRH0B(L,NY,NX)*VOLWPB
      TRH1B(L,NY,NX)=TRH1B(L,NY,NX)*VOLWPB
      TRH2B(L,NY,NX)=TRH2B(L,NY,NX)*VOLWPB
      TRH3B(L,NY,NX)=TRH3B(L,NY,NX)*VOLWPB
      TRF1B(L,NY,NX)=TRF1B(L,NY,NX)*VOLWPB
      TRF2B(L,NY,NX)=TRF2B(L,NY,NX)*VOLWPB
      TRC0B(L,NY,NX)=TRC0B(L,NY,NX)*VOLWPB
      TRC1B(L,NY,NX)=TRC1B(L,NY,NX)*VOLWPB
      TRC2B(L,NY,NX)=TRC2B(L,NY,NX)*VOLWPB
      TRM1B(L,NY,NX)=TRM1B(L,NY,NX)*VOLWPB
      TRXN4(L,NY,NX)=TRXN4(L,NY,NX)*VOLWNH
      TRXNB(L,NY,NX)=TRXNB(L,NY,NX)*VOLWNB
      TRXHY(L,NY,NX)=TRXHY(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRXAL(L,NY,NX)=TRXAL(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRXFE(L,NY,NX)=TRXFE(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRXCA(L,NY,NX)=TRXCA(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRXMG(L,NY,NX)=TRXMG(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRXNA(L,NY,NX)=TRXNA(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRXKA(L,NY,NX)=TRXKA(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRXHC(L,NY,NX)=TRXHC(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRXAL2(L,NY,NX)=TRXAL2(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRXFE2(L,NY,NX)=TRXFE2(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRXH0(L,NY,NX)=TRXH0(L,NY,NX)*VOLWPO
      TRXH1(L,NY,NX)=TRXH1(L,NY,NX)*VOLWPO
      TRXH2(L,NY,NX)=TRXH2(L,NY,NX)*VOLWPO
      TRX1P(L,NY,NX)=TRX1P(L,NY,NX)*VOLWPO
      TRX2P(L,NY,NX)=TRX2P(L,NY,NX)*VOLWPO
      TRBH0(L,NY,NX)=TRBH0(L,NY,NX)*VOLWPB
      TRBH1(L,NY,NX)=TRBH1(L,NY,NX)*VOLWPB
      TRBH2(L,NY,NX)=TRBH2(L,NY,NX)*VOLWPB
      TRB1P(L,NY,NX)=TRB1P(L,NY,NX)*VOLWPB
      TRB2P(L,NY,NX)=TRB2P(L,NY,NX)*VOLWPB
      TRALOH(L,NY,NX)=TRALOH(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRFEOH(L,NY,NX)=TRFEOH(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRCACO(L,NY,NX)=TRCACO(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRCASO(L,NY,NX)=TRCASO(L,NY,NX)*VOLWM(NPH,L,NY,NX)
      TRALPO(L,NY,NX)=TRALPO(L,NY,NX)*VOLWPO
      TRFEPO(L,NY,NX)=TRFEPO(L,NY,NX)*VOLWPO
      TRCAPD(L,NY,NX)=TRCAPD(L,NY,NX)*VOLWPO
      TRCAPH(L,NY,NX)=TRCAPH(L,NY,NX)*VOLWPO
      TRCAPM(L,NY,NX)=TRCAPM(L,NY,NX)*VOLWPO
      TRALPB(L,NY,NX)=TRALPB(L,NY,NX)*VOLWPB
      TRFEPB(L,NY,NX)=TRFEPB(L,NY,NX)*VOLWPB
      TRCPDB(L,NY,NX)=TRCPDB(L,NY,NX)*VOLWPB
      TRCPHB(L,NY,NX)=TRCPHB(L,NY,NX)*VOLWPB
      TRCPMB(L,NY,NX)=TRCPMB(L,NY,NX)*VOLWPB
C
C     BOUNDARY SALT FLUXES FOR C, H, OH, P, AL+FE, CA, OH
C     USED TO CHECK MATERIAL BALANCES IN REDIST.F
C
C     TBCO2=CO2 net change from all solute equilibria
C     TRH2O=H2O net change from all solute equilibria
C     TBION=total solute net change from all solute equilibria
C
      TBCO2(L,NY,NX)=TRCO3(L,NY,NX)
     2+TRCAC(L,NY,NX)+TRMGC(L,NY,NX)+TRNAC(L,NY,NX)+TRCACO(L,NY,NX)
     3+2.0*(TRHCO(L,NY,NX)+TRCAH(L,NY,NX)+TRMGH(L,NY,NX)) 
      TRH2O(L,NY,NX)=TRHY(L,NY,NX)+TROH(L,NY,NX)+TRXHY(L,NY,NX)
     2+TRXHC(L,NY,NX)
      TBION(L,NY,NX)=4.0*(TRH3P(L,NY,NX)+TRH3B(L,NY,NX))
     2+3.0*(TRF2P(L,NY,NX)+TRC2P(L,NY,NX)
     2+TRF2B(L,NY,NX)+TRC2B(L,NY,NX))
     3+2.0*(TRF1P(L,NY,NX)+TRC1P(L,NY,NX)+TRM1P(L,NY,NX)
     4+TRF1B(L,NY,NX)+TRC1B(L,NY,NX)+TRM1B(L,NY,NX)) 
     5+TRH0P(L,NY,NX)+TRC0P(L,NY,NX)+TRH0B(L,NY,NX)+TRC0B(L,NY,NX)
     6-(TRALPO(L,NY,NX)+TRFEPO(L,NY,NX)
     6+TRALPB(L,NY,NX)+TRFEPB(L,NY,NX))
     7-(TRCAPD(L,NY,NX)+TRCAPM(L,NY,NX)
     6+TRCPDB(L,NY,NX)+TRCPMB(L,NY,NX)
     8+5.0*(TRCAPH(L,NY,NX)+TRCPHB(L,NY,NX)))
     9+TRAL1(L,NY,NX)+TRFE1(L,NY,NX)
     1+2.0*(TRAL2(L,NY,NX)+TRFE2(L,NY,NX))
     1+3.0*(TRAL3(L,NY,NX)+TRFE3(L,NY,NX))
     2+4.0*(TRAL4(L,NY,NX)+TRFE4(L,NY,NX))
     3+TRCAO(L,NY,NX)+TRMGO(L,NY,NX)
     4+3.0*(TRALOH(L,NY,NX)+TRFEOH(L,NY,NX))
C     IF(L.EQ.11)THEN
C     WRITE(*,1111)'TRCO2',I,J,L,M,TRCO2(L,NY,NX),TRCO3(L,NY,NX)
C    2,TRHCO(L,NY,NX),TRCAC(L,NY,NX),TRMGC(L,NY,NX)
C    2,TRNAC(L,NY,NX),TRCAH(L,NY,NX)
C    2,TRMGH(L,NY,NX),TRCACO(L,NY,NX),VOLWM(NPH,L,NY,NX),RCO2
C    3,RHCO,RHCACH,RCO2Q,RCAH,RMGH,RHCO3,AHY1,AHCO31,ACO21,DPCO2
C     WRITE(*,1111)'TBION',I,J,L,M,TBION(L,NY,NX)
C     ENDIF
C
C     IF NO SALTS IS SELECTED IN SITE FILE THEN A SUBSET
C     OF THE EQUILIBRIA REACTIONS ARE SOLVED: MOSTLY THOSE
C     FOR PHOSPHORUS AND CO-REACTANTS
C
      ELSE
C
C     PRECIPITATION-DISSOLUTION CALCULATED FROM ACTIVITIES
C     OF REACTANTS AND PRODUCTS THROUGH SOLUTIONS
C     FOR THEIR EQUILIBRIUM CONSTANTS USING CURRENT
C     ION CONCENTRATION
C
C     CCEC,XCEC=cation exchange concentration,capacity
C     BLVLX=soil mass
C     C*,Z*=solute concentration, mass
C     DP*=dissociation constant from PARAMETER above
C     SP*=solubility product from PARAMETER above
C     C*<0.0=solve for C* from equilibrium with other solutes
C
      IF(BKVLX.GT.ZEROS(NY,NX))THEN
      CCEC=AMAX1(ZERO,XCEC(L,NY,NX)/BKVLX)
      ELSE
      CCEC=ZERO
      ENDIF
      CHY1=AMAX1(ZERO,10.0**(-(PH(L,NY,NX)-3.0)))
      COH1=AMAX1(ZERO,DPH2O/CHY1)
      IF(CAL(L,NY,NX).LT.0.0)THEN
      CAL1=AMAX1(ZERO,SPALO/COH1**3)
      ELSE
      CAL1=AMAX1(ZERO,AMIN1(CAL(L,NY,NX),SPALO/COH1**3))
      ENDIF
      IF(CFE(L,NY,NX).LT.0.0)THEN
      CFE1=AMAX1(ZERO,SPFEO/COH1**3)
      ELSE
      CFE1=AMAX1(ZERO,AMIN1(CFE(L,NY,NX),SPFEO/COH1**3))
      ENDIF
      CMG1=AMAX1(ZERO,ZMG(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CNA1=AMAX1(ZERO,ZNA(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      CKA1=AMAX1(ZERO,ZKA(L,NY,NX)/VOLWM(NPH,L,NY,NX))
C
C     CA CONCENTRATION FROM CURRENT CO2 CONCENTRATION
C
      CCO21=AMAX1(ZERO,CCO2S(L,NY,NX)/12.0)
      CCO31=AMAX1(ZERO,CCO21*DPCO3/CHY1**2)
      IF(CCA(L,NY,NX).LT.0.0)THEN
      CCA1=AMAX1(ZERO,AMIN1(CCAMX,SPCAC/CCO31))
      ELSE
      CCA1=AMAX1(ZERO,AMIN1(CCA(L,NY,NX),SPCAC/CCO31))
      ENDIF
C
C     PHOSPHORUS TRANSFORMATIONS IN NON-BAND SOIL ZONE
C
      IF(VOLWPO.GT.ZEROS2(NY,NX))THEN
C
C     ALUMINUM PHOSPHATE (VARISCITE)
C
C     CH2PA,CH2P1=equilibrium,current H2PO4 concentration in non-band
C     SYA0P2=solubility product derived from SPALO
C     RPALPX=H2PO4 dissolution from AlPO4 in non-band
C
      CH2PA=SYA0P2/(CAL1*COH1**2)
      RPALPX=AMAX1(-PALPO1,TPD*(CH2P1-CH2PA))
C     IF((I/30)*30.EQ.I.AND.J.EQ.12)THEN
C     WRITE(*,1117)'RPALPX',I,J,L,CH2P1,CH2PA,SYA0P2,CAL1,COH1,PALPO1
C    2,RPALPX,CAL1*CH2P1*COH1**2
C     ENDIF
C
C     IRON PHOSPHATE (STRENGITE)
C
C     CH2PF,CH2P1=equilibrium,current H2PO4 concentration in non-band
C     SYF0P2=solubility product derived from SPALO
C     RPFEPX=H2PO4 dissolution from FePO4 in non-band
C
      CH2PF=SYF0P2/(CFE1*COH1**2)
      RPFEPX=AMAX1(-PFEPO1,TPD*(CH2P1-CH2PF))
C     IF((I/30)*30.EQ.I.AND.J.EQ.12)THEN
C     WRITE(*,1117)'RPFEPX',I,J,L,CH2P1,CH2PF,SYF0P2,CFE1,COH1,PFEPO1
C    3,RPFEPX,CFE1*CH2P1*COH1**2
C     ENDIF
C
C     DICALCIUM PHOSPHATE
C
C     CH2PD,CH2P1=equilibrium,current H2PO4 concentration in non-band
C     SYCAD2=solubility product derived from SPALO
C     RPCADX=H2PO4 dissolution from CaHPO4 in non-band
C
      CH2PD=SYCAD2/(CCA1*COH1)
      RPCADX=AMAX1(-PCAPD1,TPD*(CH2P1-CH2PD))
C
C     HYDROXYAPATITE
C
C     CH2PH,CH2P1=equilibrium,current H2PO4 concentration in non-band
C     SYCAH2=solubility product derived from SPALO
C     RPCAHX=H2PO4 dissolution from apatite in non-band
C
      CH2PH=(SYCAH2/(CCA1**5*COH1**7))**0.333
      RPCAHX=AMAX1(-PCAPH1,TPD*(CH2P1-CH2PH))
C     IF((I/30)*30.EQ.i.AND.J.EQ.12)THEN
C     WRITE(*,1117)'RPCAHX',I,J,L,CH2P1,CH2PH,SYCAH2,CCA1,COH1,SPCAC
C    2,DPCO3,CCO31,CCO21,CHY1,PH(L,NY,NX),PCAPH1,RPCAHX
C    3,CCA1**5*CH2P1**3*COH1**7
C     ENDIF
C
C     MONOCALCIUM PHOSPHATE
C
C     CH2PM,CH2P1=equilibrium,current H2PO4 concentration in non-band
C     SPCAM=solubility product for Ca(H2PO4)2 
C     RPCAMX=H2PO4 dissolution from Ca(H2PO4)2 in non-band
C
      CH2PM=SQRT(SPCAM/CCA1)
      RPCAMX=AMAX1(-PCAPM1*SPPO4,TPD*(CH2P1-CH2PM))
C     IF(I.GT.315)THEN
C     WRITE(*,1117)'RPPO4',I,J,L,RPCADX,CH2P1,CH2PD,PCAPD1,RPCAHX 
C    2,CH2PA,CH2PH,SYA0P2,CAL1,COH1,SYCAH2,CCA1,CCO21,CCO31,PCAPH1 
C    3,VOLWPO,SPCAC/CCO31,CCA(L,NY,NX),H2PO4(L,NY,NX)
C    4,VOLWM(NPH,L,NY,NX),ZCA(L,NY,NX),CCO2S(L,NY,NX)
1117  FORMAT(A8,3I4,30E12.4)
C     ENDIF
C
C     PHOSPHORUS ANION EXCHANGE IN NON-BAND SOIL ZONE
C     CALCULATED FROM EXCHANGE EQUILIBRIA AMONG H2PO4-,
C     HPO4--, H+, OH- AND PROTONATED AND NON-PROTONATED -OH
C     EXCHANGE SITES
C
      IF(VOLWM(NPH,L,NY,NX).GT.ZEROS2(NY,NX))THEN
      VOLWBK=AMIN1(1.0,BKVL(L,NY,NX)/VOLWM(NPH,L,NY,NX))
      ELSE
      VOLWBK=1.0
      ENDIF
      IF(XAEC(L,NY,NX).GT.ZEROS(NY,NX))THEN
C
C     H2PO4 EXCHANGE IN NON-BAND SOIL ZONE FROM CONVERGENCE
C     SOLUTION FOR EQUILIBRIUM AMONG H2PO4-, H+, OH-, X-OH
C     AND X-H2PO4
C
C     SPH2P,SXH2P=equilibrium constant for H2PO4 exchange with R-OH2,R-OH
C     RXH2P,RYH2P=H2PO4 exchange with R-OH2,R-OH in non-band
C
      SPH2P=SXH2P*DPH2O
      RXH2P=TADA*(XOH21*CH2P1-SPH2P*XH2P1)/(XOH21+SPH2P)*VOLWBK 
      RYH2P=TADA*(XOH11*CH2P1-SXH2P*COH1*XH2P1)/(XOH11+SXH2P)*VOLWBK 
C
C     HPO4 EXCHANGE IN NON-BAND SOIL ZONE FROM CONVERGENCE
C     SOLUTION FOR EQUILIBRIUM AMONG HPO4--, H+, OH-, X-OH
C     AND X-HPO4
C
C     SPH1P=equilibrium constant for HPO4 exchange with R-OH
C     RXH1P=HPO4 exchange with R-OH in non-band
C
      SPH1P=SXH1P*DPH2O/DPH2P
      RXH1P=TADA*(XOH11*CH1P1-SPH1P*XH1P1)/(XOH11+SPH1P)*VOLWBK 
      ELSE
      RXH2P=0.0
      RYH2P=0.0
      RXH1P=0.0
      ENDIF
C     IF((I/120)*120.EQ.I.AND.J.EQ.24.AND.L.LE.6)THEN
C     WRITE(*,1116)'RXH2P',I,J,NX,NY,L,RXH2P
C    2,XOH21,CH2P1,XH2P1,XOH21*(CH2P1-RXH2P)/(XH2P1+RXH2P),SPH2P
C    3,H2PO4(L,NY,NX),RH2PX,VOLWPO
C     WRITE(*,1116)'RYH2P',I,J,NX,NY,L,RYH2P
C    2,XOH11,CH2P1,XH2P1,COH1,(XOH11*(CH2P1-RYH2P))
C    3/((XH2P1+RYH2P)*COH1),SXH2P
C     WRITE(*,1116)'RXH1P',I,J,NX,NY,L,RXH1P 
C    2,XOH11,CH1P1,XH1P1,XOH11*(CH1P1-RXH1P)/(XH1P1+RXH1P),SPH1P
C    3,SYH1P,DPH2O,DPH2P,XOH1(L,NY,NX),VLPO4(L,NY,NX),VLPOB(L,NY,NX)
C    4,TKS(L,NY,NX),XOH21,XOH01
1116  FORMAT(A8,5I4,40E12.4)
C     ENDIF
C
C     H2PO4-H+HPO4
C
C     DPH2P=dissociation constant
C     S1=equilibrium concentration in non-band
C     RH2P=H2PO4-H+HPO4 dissociation in non-band 
C     
      DP=DPH2P 
      S0=CH1P1+CHY1+DP
      S1=AMAX1(0.0,S0**2-4.0*(CH1P1*CHY1-DP*CH2P1))
      RH2P=TSL*(S0-SQRT(S1))
      ELSE
      RPALPX=0.0
      RPFEPX=0.0
      RPCADX=0.0
      RPCAHX=0.0
      RPCAMX=0.0
      RXH2P=0.0
      RYH2P=0.0
      RXH1P=0.0
      RH2P=0.0
      ENDIF
C     IF(J.EQ.1)THEN
C     WRITE(*,2222)'PO4',I,J,L,CH2P1,PALPO1,PFEPO1,PCAPD1,PCAPH1,PCAPM1
C    2,CH2PA,CH2PF,CH2PD,CH2PH,CH2PM,RPALPX,RPFEPX,RPCADX,RPCAHX,RPCAMX
C    3,XH2P1,RXH2P,RYH2P
C    3,CAL1,CFE1,CCA1,CHY1,COH1
2222  FORMAT(A8,3I4,40E12.4)
C     ENDIF
C
C     PHOSPHORUS PRECIPITATION-DISSOLUTION IN BAND SOIL ZONE
C
      IF(VOLWPB.GT.ZEROS2(NY,NX))THEN
C
C     ALUMINUM PHOSPHATE (VARISCITE)
C
C     CH2PA,CH2PB=equilibrium,current H2PO4 concentration in band
C     SYA0P2=solubility product derived from SPALO
C     RPALBX=H2PO4 dissolution from AlPO4 in band
C
      CH2PA=SYA0P2/(CAL1*COH1**2)
      RPALBX=AMAX1(-PALPOB,TPD*(CH2PB-CH2PA))
C
C     IRON PHOSPHATE (STRENGITE)
C
C     CH2PF,CH2PB=equilibrium,current H2PO4 concentration in band
C     SYF0P2=solubility product derived from SPALO
C     RPFEBX=H2PO4 dissolution from FePO4 in band
C
      CH2PF=SYF0P2/(CFE1*COH1**2)
      RPFEBX=AMAX1(-PFEPOB,TPD*(CH2PB-CH2PF))
C
C     DICALCIUM PHOSPHATE
C
C     CH2PD,CH2PB=equilibrium,current H2PO4 concentration in band
C     SYCAD2=solubility product derived from SPALO
C     RPCDBX=H2PO4 dissolution from CaHPO4 in band
C
      CH2PD=SYCAD2/(CCA1*COH1)
      RPCDBX=AMAX1(-PCAPDB,TPD*(CH2PB-CH2PD))
C
C     HYDROXYAPATITE
C
C     CH2PH,CH2PB=equilibrium,current H2PO4 concentration in band
C     SYCAH2=solubility product derived from SPALO
C     RPCHBX=H2PO4 dissolution from apatite in band
C
      CH2PH=(SYCAH2/(CCA1**5*COH1**7))**0.333
      RPCHBX=AMAX1(-PCAPHB,TPD*(CH2PB-CH2PH))
C
C     MONOCALCIUM PHOSPHATE
C
C     CH2PM,CH2PB=equilibrium,current H2PO4 concentration in band
C     SPCAM=solubility product for Ca(H2PO4)2 
C     RPCMBX=H2PO4 dissolution from Ca(H2PO4)2 in band
C
      CH2PM=SQRT(SPCAM/CCA1)
      RPCMBX=AMAX1(-PCAPMB*SPPO4,TPD*(CH2PB-CH2PM))
C     IF(I.GT.315)THEN
C     WRITE(*,1117)'RPPOB',I,J,L,RPCMBX,CH2PM,CH2PB,SPCAM,CCA1
C    2,PCAPMB,SPPO4,TPD,PCPMB(L,NY,NX),BKVLPB
C     ENDIF
C
C     PHOSPHORUS ANION EXCHANGE IN BAND SOIL ZONE
C     CALCULATED FROM EXCHANGE EQUILIBRIA AMONG H2PO4-,
C     HPO4--, H+, OH- AND PROTONATED AND NON-PROTONATED -OH
C     EXCHANGE SITES
C
      IF(XAEC(L,NY,NX).GT.ZEROS(NY,NX))THEN
C
C     H2PO4 EXCHANGE IN BAND SOIL ZONE FROM CONVERGENCE
C     SOLUTION FOR EQUILIBRIUM AMONG H2PO4-, H+, OH-, X-OH
C     AND X-H2PO4
C
C     SPH2P,SXH2P=equilibrium constant for H2PO4 exchange with R-OH2,R-OH
C     RXH2B,RYH2B=H2PO4 exchange with R-OH2,R-OH in band
C
      SPH2P=SXH2P*DPH2O
      RXH2B=TADA*(XH21B*CH2PB-SPH2P*X2P1B)/(XH21B+SPH2P)*VOLWBK 
      RYH2B=TADA*(XH11B*CH2PB-SXH2P*X2P1B*COH1)/(XH11B+SXH2P)*VOLWBK 
C
C     HPO4 EXCHANGE IN BAND SOIL ZONE FROM CONVERGENCE
C     SOLUTION FOR EQUILIBRIUM AMONG HPO4--, H+, OH-, X-OH
C     AND X-HPO4
C
C     SPH1P=equilibrium constant for HPO4 exchange with R-OH
C     RXH1B=HPO4 exchange with R-OH in band
C
      SPH1P=SXH1P*DPH2O/DPH2P
      RXH1B=TADA*(XH11B*CH1PB-SPH1P*X1P1B)/(XH11B+SPH1P)*VOLWBK 
      ELSE
      RXH2B=0.0
      RYH2B=0.0
      RXH1B=0.0
      ENDIF
C     WRITE(*,2224)'RXH1B',I,J,L,RXH1B,XH11B,CH1PB,SPH1P,X1P1B
2224  FORMAT(A8,3I4,40E12.4)
C
C     H2PO4-H+HPO4
C
C     DPH2P=dissociation constant
C     S1=eqilibriunm concentration in band
C     RH2B=H2PO4-H+HPO4 dissociation in band 
C
      DP=DPH2P 
      S0=CH1PB+CHY1+DP
      S1=AMAX1(0.0,S0**2-4.0*(CH1PB*CHY1-DP*CH2PB))
      RH2B=TSLX*(S0-SQRT(S1))
      ELSE
      RPALBX=0.0
      RPFEBX=0.0
      RPCDBX=0.0
      RPCHBX=0.0
      RPCMBX=0.0
      RXH2B=0.0
      RYH2B=0.0
      RXH1B=0.0
      RH2B=0.0
      ENDIF
C
C     CATION EXCHANGE FROM GAPON SELECTIVITY COEFFICIENTS
C     FOR CA-NH4, CA-H, CA-AL
C
      IF(XCEC(L,NY,NX).GT.ZEROS(NY,NX))THEN
C
C     CATION CONCENTRATIONS
C
C     EQUILIBRIUM X-CA CONCENTRATION FROM CEC, GAPON COEFFICIENTS
C     AND CATION CONCENTRATIONS
C
C     CCEC,XCEC=cation exchange concentration,capacity
C     XCAX=equilibrium R-Ca concentration
C     GKC4,GKCH,GKCA,GKCM,GKCN,GKCK=Gapon selectivity coefficients for
C     CA-NH4,CA-H,CA-AL,CA-MG,CA-NA,CA-K
C     X*Q=equilibrium exchangeable concentrations
C     XTLQ=total equilibrium exchangeable concentration 
C
      CN41=AMAX1(ZERO,CN41)
      CN4B=AMAX1(ZERO,CN4B)
      CALX=AMAX1(ZERO,CAL1)**0.333
      CFEX=AMAX1(ZERO,CFE1)**0.333
      CCAX=AMAX1(ZERO,CCA1)**0.500
      CMGX=AMAX1(ZERO,CMG1)**0.500
      CNA1=AMAX1(ZERO,CNA1)
      CKA1=AMAX1(ZERO,CKA1)
C
C     EQUILIBRIUM X-CA CONCENTRATION FROM CEC AND CATION
C     CONCENTRATIONS
C
      XCAX=CCEC/(1.0+GKC4(L,NY,NX)*CN41/CCAX*VLNH4(L,NY,NX)
     2+GKC4(L,NY,NX)*CN4B/CCAX*VLNHB(L,NY,NX)
     3+GKCH(L,NY,NX)*CHY1/CCAX+GKCA(L,NY,NX)*CALX/CCAX 
     3+GKCA(L,NY,NX)*CFEX/CCAX+GKCM(L,NY,NX)*CMGX/CCAX
     3+GKCN(L,NY,NX)*CNA1/CCAX+GKCK(L,NY,NX)*CKA1/CCAX)
      XN4Q=XCAX*CN41*GKC4(L,NY,NX)
      XNBQ=XCAX*CN4B*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*VLNH4(L,NY,NX)+XNBQ*VLNHB(L,NY,NX)
     2+XHYQ+XALQ+XFEQ+XCAQ+XMGQ+XNAQ+XKAQ
      IF(XTLQ.GT.ZERO)THEN
      FX=CCEC/XTLQ
      ELSE
      FX=0.0
      ENDIF
      XN4Q=FX*XN4Q
      XNBQ=FX*XNBQ
C
C     NH4 EXCHANGE IN NON-BAND AND BAND SOIL ZONES
C
C     RXN4,RXNB=NH4 adsorption in non-band,band
C     TADC=adsorption rate constant
C
      RXN4=TADC*AMAX1(AMIN1((XN4Q-XN41)*CN41/XN4Q,CN41),-XN41) 
      RXNB=TADC*AMAX1(AMIN1((XNBQ-XN4B)*CN4B/XNBQ,CN4B),-XN4B)
      ELSE
      RXN4=0.0
      RXNB=0.0
      ENDIF
C     IF(J.EQ.12.AND.L.EQ.0)THEN
C     WRITE(*,2222)'RXN4',I,J,L,RXN4,CN41,XN41,CCAX,CCA1,XCAQ 
C    2,CCEC,XCAX,FN4X,FCAQ,GKC4(L,NY,NX),PH(L,NY,NX),VOLWBK
C    3,XN4Q,XTLQ,FX
C    3,(CCA1)**0.5*XN41/(CN41*XCAQ),ZCA(L,NY,NX),BKVLX
C    4,CN4B,CHY1,CALX,CFEX,CMGX,CNA1,CKA1
C     ENDIF
C
C     NH4-NH3+H IN NON-BAND AND BAND SOIL ZONES
C
C     RNH4,RNHB=NH4-NH3+H dissociation in non-band,band
C     DPN4=NH4 dissociation constant
C
      IF(VOLWNH.GT.ZEROS2(NY,NX))THEN
      RNH4=(CHY1*CN31-DPN4*CN41)/(DPN4+CHY1)
      ELSE
      RNH4=0.0
      ENDIF
      IF(VOLWNB.GT.ZEROS2(NY,NX))THEN
      RNHB=(CHY1*CN3B-DPN4*CN4B)/(DPN4+CHY1)
      ELSE
      RNHB=0.0
      ENDIF
C     IF(IYRC.EQ.2012.AND.I.EQ.151.AND.NX.EQ.1)THEN
C     WRITE(*,2222)'RNH4',I,J,L,RNH4,CHY1,CN31,DPN4,CN41
C    2,RXN4,XN41,VOLWNH,RNHB,CN3B,CN4B,VOLWNB,RXNB,XN4B,FN4X
C    2,CN41*VOLWNH,XN41*VOLWNH,CN4B*VOLWNB,XN4B*VOLWNB 
C    3,(CCA1)**0.5*XN41/(CN41*XCAQ),(CCA1)**0.5*XN4B/(CN4B*XCAQ)
C    4,RN4X,RN3X,RNBX,R3BX,ZEROS2(NY,NX)
C     ENDIF
C
C     TOTAL ION FLUXES FOR ALL REACTIONS ABOVE
C
C     RN4S,RN4B=net NH4 flux in non-band,band
C     RN3S,RN3B=net NH3 flux in non-band,band
C     RAL,RFE,RHY,RCA,RMG,RNA,RKA,ROH=net Al,Fe,H,Ca,Mg,Na,K,OH flux
C     RHP1,RHP2=net HPO4,H2PO4 flux in non-band
C     RXH1,RXH2,RX1P,RX2P=net R-OH,R-OH2,R-HPO4,R-H2PO4 in non-band
C     RHB1,RHB2=net HPO4,H2PO4 flux in band
C     RBH1,RBH2,RB1P,RB2P=net R-OH,R-OH2,R-HPO4,R-H2PO4 in band
C
      RN4S=RNH4-RXN4
      RN4B=RNHB-RXNB
      RN3S=-RNH4
      RN3B=-RNHB
      RHP1=-RH2P-RXH1P 
      RHP2=RH2P-RXH2P-RYH2P
     2-RPALPX-RPFEPX-RPCADX-2.0*RPCAMX-3.0*RPCAHX
      RHB1=-RH2B-RXH1B 
      RHB2=RH2B-RXH2B-RYH2B
     2-RPALBX-RPFEBX-RPCDBX-2.0*RPCMBX-3.0*RPCHBX
      RXH1=-RYH2P-RXH1P
      RXH2=-RXH2P
      RX1P=RXH1P
      RX2P=RXH2P+RYH2P
      RBH1=-RYH2B-RXH1B
      RBH2=-RXH2B
      RB1P=RXH1B
      RB2P=RXH2B+RYH2B
C
C     CONVERT TOTAL ION FLUXES FROM CHANGES IN CONCENTRATION
C     TO CHANGES IN MASS PER UNIT AREA FOR USE IN 'REDIST'
C
C     TRN4S,TRN4B=total NH4 flux in non-band,band
C     TRN3S,TRN3B=total NH3 flux in non-band,band
C     TRH1P,TRH2P=net HPO4,H2PO4 flux in non-band
C     TRH1B,TRH2B=net HPO4,H2PO4 flux in band
C     TRNX4,TRNXB=total NH4 adsorption in non-band,band
C     TRXH1,TRXH2,TRX1P,TRX2P
C     =total R-OH,R-OH2,R-HPO4,R-H2PO4 adsorption in non-band
C     TRBH1,TRBH2,TRB1P,TRB2P
C     =total R-OH,R-OH2,R-HPO4,R-H2PO4 adsorption in band
C     TRALPO,TRFEPO,TRCAPD,TRCAPH,TRCAPM
C     =total AlPO4,FePO4,CaHPO4,apatite,Ca(H2PO4)2 precipitation in non-band
C     TRALPB,TRFEPB,TRCPDB,TRCPHB,TRCPMB
C     =total AlPO4,FePO4,CaHPO4,apatite,Ca(H2PO4)2 precipitation in band
C
      TRN4S(L,NY,NX)=TRN4S(L,NY,NX)+RN4S*VOLWNH
      TRN4B(L,NY,NX)=TRN4B(L,NY,NX)+RN4B*VOLWNB
      TRN3S(L,NY,NX)=TRN3S(L,NY,NX)+RN3S*VOLWNH
      TRN3B(L,NY,NX)=TRN3B(L,NY,NX)+RN3B*VOLWNB
      TRH1P(L,NY,NX)=TRH1P(L,NY,NX)+RHP1*VOLWPO
      TRH2P(L,NY,NX)=TRH2P(L,NY,NX)+RHP2*VOLWPO
      TRH1B(L,NY,NX)=TRH1B(L,NY,NX)+RHB1*VOLWPB
      TRH2B(L,NY,NX)=TRH2B(L,NY,NX)+RHB2*VOLWPB
      TRXN4(L,NY,NX)=TRXN4(L,NY,NX)+RXN4*VOLWNH
      TRXNB(L,NY,NX)=TRXNB(L,NY,NX)+RXNB*VOLWNB
      TRXH1(L,NY,NX)=TRXH1(L,NY,NX)+RXH1*VOLWPO
      TRXH2(L,NY,NX)=TRXH2(L,NY,NX)+RXH2*VOLWPO
      TRX1P(L,NY,NX)=TRX1P(L,NY,NX)+RX1P*VOLWPO
      TRX2P(L,NY,NX)=TRX2P(L,NY,NX)+RX2P*VOLWPO
      TRBH1(L,NY,NX)=TRBH1(L,NY,NX)+RBH1*VOLWPB
      TRBH2(L,NY,NX)=TRBH2(L,NY,NX)+RBH2*VOLWPB
      TRB1P(L,NY,NX)=TRB1P(L,NY,NX)+RB1P*VOLWPB
      TRB2P(L,NY,NX)=TRB2P(L,NY,NX)+RB2P*VOLWPB
      TRALPO(L,NY,NX)=TRALPO(L,NY,NX)+RPALPX*VOLWPO
      TRFEPO(L,NY,NX)=TRFEPO(L,NY,NX)+RPFEPX*VOLWPO
      TRCAPD(L,NY,NX)=TRCAPD(L,NY,NX)+RPCADX*VOLWPO
      TRCAPH(L,NY,NX)=TRCAPH(L,NY,NX)+RPCAHX*VOLWPO
      TRCAPM(L,NY,NX)=TRCAPM(L,NY,NX)+RPCAMX*VOLWPO
      TRALPB(L,NY,NX)=TRALPB(L,NY,NX)+RPALBX*VOLWPB
      TRFEPB(L,NY,NX)=TRFEPB(L,NY,NX)+RPFEBX*VOLWPB
      TRCPDB(L,NY,NX)=TRCPDB(L,NY,NX)+RPCDBX*VOLWPB
      TRCPHB(L,NY,NX)=TRCPHB(L,NY,NX)+RPCHBX*VOLWPB
      TRCPMB(L,NY,NX)=TRCPMB(L,NY,NX)+RPCMBX*VOLWPB
C     IF(IYRC.EQ.2012.AND.I.EQ.151.AND.NX.EQ.1)THEN
C     WRITE(*,24)'RN4S',I,J,L,RN4S,RN3S,RNH4,RXN4,VOLWNH 
C     WRITE(*,24)'RHP1',I,J,L,RHP1,RH2P,RXH1P
C    2,TRX1P(L,NY,NX),TRH2P(L,NY,NX)
C     WRITE(*,24)'RHP2',I,J,L,RHP2,RH2P,RXH2P,RYH2P
C    2,RPALPX,RPFEPX,RPCADX,2.0*RPCAMX,3.0*RPCAHX 
C    3,TRX2P(L,NY,NX)
24    FORMAT(A8,3I4,60E12.4)
C     ENDIF
      ENDIF
C
C     CHANGE IN WIDTHS AND DEPTHS OF FERTILIZER BANDS FROM
C     VERTICAL AND HORIZONTAL DIFFUSION DRIVEN BY CONCENTRATION
C     DIFFERENCES BETWEEN BAND AND NON-BAND SOIL ZONES
C
C     ROWI=band row width
C     FLWD=net vertical flow relative to area  
C
C     IF(ROWI(I,NY,NX).GT.0.0)THEN
      FLWD=0.5*(FLW(3,L,NY,NX)+FLW(3,L+1,NY,NX))/AREA(3,L,NY,NX)
C
C     NH4 FERTILIZER BAND
C
C     IFNHB=banded NH4 fertilizer flag
C     ROWN=NH4 fertilizer band row width
C     DPNH4=NH4 fertilizer band depth 
C
      IF(IFNHB(NY,NX).EQ.1.AND.ROWN(NY,NX).GT.0.0)THEN
      IF(L.EQ.NU(NY,NX).OR.CDPTH(L-1,NY,NX).LT.DPNH4(NY,NX))THEN
C
C     NH4 BAND WIDTH
C
C     DWNH4=change in NH4 fertilizer band width
C     WDNHB=layer NH4 fertilizer band width
C     ZNSGL=NH4 diffusivity
C     TORT=tortuosity
C
      DWNH4=0.5*SQRT(ZNSGL(L,NY,NX))*TORT(NPH,L,NY,NX)
      WDNHB(L,NY,NX)=AMIN1(ROWN(NY,NX)
     2,AMAX1(0.025,WDNHB(L,NY,NX))+DWNH4)
C
C     NH4 BAND DEPTH
C
C     DPFLW=change in NH4 fertilizer band depth
C     DPNH4,DPNHB=total,layer NH4 fertilizer band depth
C
      IF(CDPTH(L,NY,NX).GE.DPNH4(NY,NX))THEN
      DPFLW=FLWD+DWNH4
      DPNH4(NY,NX)=DPNH4(NY,NX)+DPFLW
      DPNHB(L,NY,NX)=DPNHB(L,NY,NX)+DPFLW
      IF(DPNHB(L,NY,NX).GT.DLYR(3,L,NY,NX))THEN
      DPNHB(L+1,NY,NX)=DPNHB(L+1,NY,NX)+(DPNHB(L,NY,NX)-DLYR(3,L,NY,NX))
      WDNHB(L+1,NY,NX)=WDNHB(L,NY,NX)
      DPNHB(L,NY,NX)=DLYR(3,L,NY,NX)
      ELSEIF(DPNHB(L,NY,NX).LT.0.0)THEN
      DPNHB(L-1,NY,NX)=DPNHB(L-1,NY,NX)+DPNHB(L,NY,NX)
      DPNHB(L,NY,NX)=0.0
      WDNHB(L,NY,NX)=0.0
      ENDIF
      ENDIF
C
C     FRACTION OF SOIL LAYER OCCUPIED BY NH4 BAND
C     FROM BAND WIDTH X DEPTH
C
C     VLNH4,VLNHB=fraction of soil volume in NH4 non-band,band
C     DLYR=soil layer thickness
C     FVLNH4=relative change in VLNH4
C
      XVLNH4=VLNH4(L,NY,NX)
      IF(DLYR(3,L,NY,NX).GT.ZERO)THEN
      VLNHB(L,NY,NX)=AMAX1(0.0,AMIN1(0.999,WDNHB(L,NY,NX) 
     2/ROWN(NY,NX)*DPNHB(L,NY,NX)/DLYR(3,L,NY,NX)))
      ELSE
      VLNHB(L,NY,NX)=0.0
      ENDIF
      VLNH4(L,NY,NX)=1.0-VLNHB(L,NY,NX)
      FVLNH4=AMIN1(0.0,(VLNH4(L,NY,NX)-XVLNH4)/XVLNH4)
C
C     TRANSFER NH4, NH3 FROM NON-BAND TO BAND
C     DURING BAND GROWTH
C
C     DNH4S,DNH3S,DXNH4=transfer of NH4,NH3,exchangeable NH4
C
      DNH4S=FVLNH4*ZNH4S(L,NY,NX)/14.0
      DNH3S=FVLNH4*ZNH3S(L,NY,NX)/14.0
      DXNH4=FVLNH4*XN4(L,NY,NX)
      TRN4S(L,NY,NX)=TRN4S(L,NY,NX)+DNH4S
      TRN4B(L,NY,NX)=TRN4B(L,NY,NX)-DNH4S
      TRN3S(L,NY,NX)=TRN3S(L,NY,NX)+DNH3S
      TRN3B(L,NY,NX)=TRN3B(L,NY,NX)-DNH3S
      TRXN4(L,NY,NX)=TRXN4(L,NY,NX)+DXNH4
      TRXNB(L,NY,NX)=TRXNB(L,NY,NX)-DXNH4
      ELSE
C
C     AMALGAMATE NH4 BAND WITH NON-BAND IF BAND NO LONGER EXISTS
C
      DPNHB(L,NY,NX)=0.0
      WDNHB(L,NY,NX)=0.0
      VLNH4(L,NY,NX)=1.0
      VLNHB(L,NY,NX)=0.0
      ZNH4S(L,NY,NX)=ZNH4S(L,NY,NX)+ZNH4B(L,NY,NX)
      ZNH3S(L,NY,NX)=ZNH3S(L,NY,NX)+ZNH3B(L,NY,NX)
      ZNH4B(L,NY,NX)=0.0
      ZNH3B(L,NY,NX)=0.0
      XN4(L,NY,NX)=XN4(L,NY,NX)+XNB(L,NY,NX)
      XNB(L,NY,NX)=0.0
      ENDIF
      ENDIF
C
C     NO3 FERTILIZER BAND
C
C     IFNOB=banded NO3 fertilizer flag
C     ROWO=NO3 fertilizer band row width
C     DPNO3=NO3 fertilizer band depth 
C
      IF(IFNOB(NY,NX).EQ.1.AND.ROWO(NY,NX).GT.0.0)THEN
      IF(L.EQ.NU(NY,NX).OR.CDPTH(L-1,NY,NX).LT.DPNO3(NY,NX))THEN
C
C     NO3 BAND WIDTH
C
C     DWNO3=change in NO3 fertilizer band width
C     WDNOB=layer NO3 fertilizer band width
C     ZOSGL=NO3 diffusivity
C     TORT=tortuosity
C
      DWNO3=0.5*SQRT(ZOSGL(L,NY,NX))*TORT(NPH,L,NY,NX)
      WDNOB(L,NY,NX)=AMIN1(ROWO(NY,NX),WDNOB(L,NY,NX)+DWNO3)
C
C     NO3 BAND DEPTH
C
C     DPFLW=change in NO3 fertilizer band depth
C     DPNO3,DPNOB=total,layer NO3 fertilizer band depth
C
      IF(CDPTH(L,NY,NX).GE.DPNO3(NY,NX))THEN
      DPFLW=FLWD+DWNO3
      DPNO3(NY,NX)=DPNO3(NY,NX)+DPFLW
      DPNOB(L,NY,NX)=DPNOB(L,NY,NX)+DPFLW
      IF(DPNOB(L,NY,NX).GT.DLYR(3,L,NY,NX))THEN
      DPNOB(L+1,NY,NX)=DPNOB(L+1,NY,NX)+(DPNOB(L,NY,NX)-DLYR(3,L,NY,NX))
      WDNOB(L+1,NY,NX)=WDNOB(L,NY,NX)
      DPNOB(L,NY,NX)=DLYR(3,L,NY,NX)
      ELSE IF(DPNOB(L,NY,NX).LT.0.0)THEN
      DPNOB(L-1,NY,NX)=DPNOB(L-1,NY,NX)+DPNOB(L,NY,NX)
      DPNOB(L,NY,NX)=0.0
      WDNOB(L,NY,NX)=0.0
      ENDIF
      ENDIF
C
C     FRACTION OF SOIL LAYER OCCUPIED BY NO3 BAND
C     FROM BAND WIDTH X DEPTH
C
C     VLNO3,VLNOB=fraction of soil volume in NO3 non-band,band
C     DLYR=soil layer thickness
C     FVLNO3=relative change in VLNO3
C
      XVLNO3=VLNO3(L,NY,NX)
      IF(DLYR(3,L,NY,NX).GT.ZERO)THEN
      VLNOB(L,NY,NX)=AMAX1(0.0,AMIN1(0.999,WDNOB(L,NY,NX) 
     2/ROWO(NY,NX)*DPNOB(L,NY,NX)/DLYR(3,L,NY,NX)))
      ELSE
      VLNOB(L,NY,NX)=0.0
      ENDIF
      VLNO3(L,NY,NX)=1.0-VLNOB(L,NY,NX)
      FVLNO3=AMIN1(0.0,(VLNO3(L,NY,NX)-XVLNO3)/XVLNO3)
C
C     TRANSFER NO3 FROM NON-BAND TO BAND
C     DURING BAND GROWTH
C
C     DNO3S,DNO2S=transfer of NO3,NO2
C
      DNO3S=FVLNO3*ZNO3S(L,NY,NX)/14.0
      DNO2S=FVLNO3*ZNO2S(L,NY,NX)/14.0
      TRNO3(L,NY,NX)=TRNO3(L,NY,NX)+DNO3S
      TRNO2(L,NY,NX)=TRNO2(L,NY,NX)+DNO2S
      TRNOB(L,NY,NX)=TRNOB(L,NY,NX)-DNO3S
      TRN2B(L,NY,NX)=TRN2B(L,NY,NX)-DNO2S
      ELSE
C
C     AMALGAMATE NO3 BAND WITH NON-BAND IF BAND NO LONGER EXISTS 
C
      DPNOB(L,NY,NX)=0.0
      WDNOB(L,NY,NX)=0.0
      VLNO3(L,NY,NX)=1.0
      VLNOB(L,NY,NX)=0.0
      ZNO3S(L,NY,NX)=ZNO3S(L,NY,NX)+ZNO3B(L,NY,NX)
      ZNO2S(L,NY,NX)=ZNO2S(L,NY,NX)+ZNO2B(L,NY,NX)
      ZNO3B(L,NY,NX)=0.0
      ZNO2B(L,NY,NX)=0.0
      ENDIF
      ENDIF
C
C     PO4 FERTILIZER BAND
C
C     IFPOB=banded H2PO4 fertilizer flag
C     ROWP=H2PO4 fertilizer band row width
C     DPPO4=H2PO4 fertilizer band depth 
C
      IF(IFPOB(NY,NX).EQ.1.AND.ROWP(NY,NX).GT.0.0)THEN
      IF(L.EQ.NU(NY,NX).OR.CDPTH(L-1,NY,NX).LT.DPPO4(NY,NX))THEN
C
C     PO4 BAND WIDTH
C
C     DWPO4=change in H2PO4 fertilizer band width
C     WDPO4=layer H2PO4 fertilizer band width
C     POSGL=H2PO4 diffusivity
C     TORT=tortuosity
C
      DWPO4=0.5*SQRT(POSGL(L,NY,NX))*TORT(NPH,L,NY,NX)
      WDPOB(L,NY,NX)=AMIN1(ROWP(NY,NX),WDPOB(L,NY,NX)+DWPO4)
C
C     PO4 BAND DEPTH
C
C     DPFLW=change in H2PO4 fertilizer band depth
C     DPPO4,DPPOB=total,layer H2PO4 fertilizer band depth
C
      IF(CDPTH(L,NY,NX).GE.DPPO4(NY,NX))THEN
      DPFLW=FLWD+DWPO4
      DPPO4(NY,NX)=DPPO4(NY,NX)+DPFLW
      DPPOB(L,NY,NX)=DPPOB(L,NY,NX)+DPFLW
      IF(DPPOB(L,NY,NX).GT.DLYR(3,L,NY,NX))THEN
      DPPOB(L+1,NY,NX)=DPPOB(L+1,NY,NX)+(DPPOB(L,NY,NX)-DLYR(3,L,NY,NX))
      WDPOB(L+1,NY,NX)=WDPOB(L,NY,NX)
      DPPOB(L,NY,NX)=DLYR(3,L,NY,NX)
      ELSEIF(DPPOB(L,NY,NX).LT.0.0)THEN
      DPPOB(L-1,NY,NX)=DPPOB(L-1,NY,NX)+DPPOB(L,NY,NX)
      DPPOB(L,NY,NX)=0.0
      WDPOB(L,NY,NX)=0.0
      ENDIF
      ENDIF
C
C     FRACTION OF SOIL LAYER OCCUPIED BY PO4 BAND
C     FROM BAND WIDTH X DEPTH
C
C     VLPO43,VLPOB=fraction of soil volume in H2PO4 non-band,band
C     DLYR=soil layer thickness
C     FVLPO4=relative change in VLPO4
C
      XVLPO4=VLPO4(L,NY,NX)
      IF(DLYR(3,L,NY,NX).GT.ZERO)THEN
      VLPOB(L,NY,NX)=AMAX1(0.0,AMIN1(0.999,WDPOB(L,NY,NX)
     2/ROWP(NY,NX)*DPPOB(L,NY,NX)/DLYR(3,L,NY,NX)))
      ELSE
      VLPOB(L,NY,NX)=0.0
      ENDIF
      VLPO4(L,NY,NX)=1.0-VLPOB(L,NY,NX)
      FVLPO4=AMIN1(0.0,(VLPO4(L,NY,NX)-XVLPO4)/XVLPO4)
C
C     TRANSFER HPO4,H2PO4 FROM NON-BAND TO BAND
C     DURING BAND GROWTH DEPENDING ON SALT
C     VS. NON-SALT OPTION
C
C     DZ*,DX*,DP*=transfer of solute,adsorbed,precipitated HPO4,H2PO4
C
      IF(ISALTG.NE.0)THEN
      DZH0P=FVLPO4*H0PO4(L,NY,NX)
      DZH1P=FVLPO4*H1PO4(L,NY,NX)/31.0
      DZH2P=FVLPO4*H2PO4(L,NY,NX)/31.0
      DZH3P=FVLPO4*H3PO4(L,NY,NX)
      DZF1P=FVLPO4*ZFE1P(L,NY,NX)
      DZF2P=FVLPO4*ZFE2P(L,NY,NX)
      DZC0P=FVLPO4*ZCA0P(L,NY,NX)
      DZC1P=FVLPO4*ZCA1P(L,NY,NX)
      DZC2P=FVLPO4*ZCA2P(L,NY,NX)
      DZM1P=FVLPO4*ZMG1P(L,NY,NX)
      DXOH0=FVLPO4*XOH0(L,NY,NX)
      DXOH1=FVLPO4*XOH1(L,NY,NX)
      DXOH2=FVLPO4*XOH2(L,NY,NX)
      DXH1P=FVLPO4*XH1P(L,NY,NX)
      DXH2P=FVLPO4*XH2P(L,NY,NX)
      DPALP=FVLPO4*PALPO(L,NY,NX)
      DPFEP=FVLPO4*PFEPO(L,NY,NX)
      DPCDP=FVLPO4*PCAPD(L,NY,NX)
      DPCHP=FVLPO4*PCAPH(L,NY,NX)
      DPCMP=FVLPO4*PCAPM(L,NY,NX)
      TRH0P(L,NY,NX)=TRH0P(L,NY,NX)+DZH0P
      TRH1P(L,NY,NX)=TRH1P(L,NY,NX)+DZH1P
      TRH2P(L,NY,NX)=TRH2P(L,NY,NX)+DZH2P
      TRH3P(L,NY,NX)=TRH3P(L,NY,NX)+DZH3P
      TRF1P(L,NY,NX)=TRF1P(L,NY,NX)+DZF1P
      TRF2P(L,NY,NX)=TRF2P(L,NY,NX)+DZF2P
      TRC0P(L,NY,NX)=TRC0P(L,NY,NX)+DZC0P
      TRC1P(L,NY,NX)=TRC1P(L,NY,NX)+DZC1P
      TRC2P(L,NY,NX)=TRC2P(L,NY,NX)+DZC2P
      TRM1P(L,NY,NX)=TRM1P(L,NY,NX)+DZM1P
      TRH0B(L,NY,NX)=TRH0B(L,NY,NX)-DZH0P
      TRH1B(L,NY,NX)=TRH1B(L,NY,NX)-DZH1P
      TRH2B(L,NY,NX)=TRH2B(L,NY,NX)-DZH2P
      TRH3B(L,NY,NX)=TRH3B(L,NY,NX)-DZH3P
      TRF1B(L,NY,NX)=TRF1B(L,NY,NX)-DZF1P
      TRF2B(L,NY,NX)=TRF2B(L,NY,NX)-DZF2P
      TRC0B(L,NY,NX)=TRC0B(L,NY,NX)-DZC0P
      TRC1B(L,NY,NX)=TRC1B(L,NY,NX)-DZC1P
      TRC2B(L,NY,NX)=TRC2B(L,NY,NX)-DZC2P
      TRM1B(L,NY,NX)=TRM1B(L,NY,NX)-DZM1P
      TRXH0(L,NY,NX)=TRXH0(L,NY,NX)+DXOH0
      TRXH1(L,NY,NX)=TRXH1(L,NY,NX)+DXOH1
      TRXH2(L,NY,NX)=TRXH2(L,NY,NX)+DXOH2
      TRX1P(L,NY,NX)=TRX1P(L,NY,NX)+DXH1P
      TRX2P(L,NY,NX)=TRX2P(L,NY,NX)+DXH2P
      TRBH0(L,NY,NX)=TRBH0(L,NY,NX)-DXOH0
      TRBH1(L,NY,NX)=TRBH1(L,NY,NX)-DXOH1
      TRBH2(L,NY,NX)=TRBH2(L,NY,NX)-DXOH2
      TRB1P(L,NY,NX)=TRB1P(L,NY,NX)-DXH1P
      TRB2P(L,NY,NX)=TRB2P(L,NY,NX)-DXH2P
      TRALPO(L,NY,NX)=TRALPO(L,NY,NX)+DPALP
      TRFEPO(L,NY,NX)=TRFEPO(L,NY,NX)+DPFEP
      TRCAPD(L,NY,NX)=TRCAPD(L,NY,NX)+DPCDP
      TRCAPH(L,NY,NX)=TRCAPH(L,NY,NX)+DPCHP
      TRCAPM(L,NY,NX)=TRCAPM(L,NY,NX)+DPCMP
      TRALPB(L,NY,NX)=TRALPB(L,NY,NX)-DPALP
      TRFEPB(L,NY,NX)=TRFEPB(L,NY,NX)-DPFEP
      TRCPDB(L,NY,NX)=TRCPDB(L,NY,NX)-DPCDP
      TRCPHB(L,NY,NX)=TRCPHB(L,NY,NX)-DPCHP
      TRCPMB(L,NY,NX)=TRCPMB(L,NY,NX)-DPCMP
      ELSE
      DZH1P=FVLPO4*H1PO4(L,NY,NX)/31.0
      DZH2P=FVLPO4*H2PO4(L,NY,NX)/31.0
      DXOH1=FVLPO4*XOH1(L,NY,NX)
      DXOH2=FVLPO4*XOH2(L,NY,NX)
      DXH2P=FVLPO4*XH2P(L,NY,NX)
      DPALP=FVLPO4*PALPO(L,NY,NX)
      DPFEP=FVLPO4*PFEPO(L,NY,NX)
      DPCDP=FVLPO4*PCAPD(L,NY,NX)
      DPCHP=FVLPO4*PCAPH(L,NY,NX)
      DPCMP=FVLPO4*PCAPM(L,NY,NX)
      TRH1P(L,NY,NX)=TRH1P(L,NY,NX)+DZH1P
      TRH2P(L,NY,NX)=TRH2P(L,NY,NX)+DZH2P
      TRXH1(L,NY,NX)=TRXH1(L,NY,NX)+DXOH1
      TRXH2(L,NY,NX)=TRXH2(L,NY,NX)+DXOH2
      TRX2P(L,NY,NX)=TRX2P(L,NY,NX)+DXH2P
      TRH1B(L,NY,NX)=TRH1B(L,NY,NX)-DZH1P
      TRH2B(L,NY,NX)=TRH2B(L,NY,NX)-DZH2P
      TRBH1(L,NY,NX)=TRBH1(L,NY,NX)-DXOH1
      TRBH2(L,NY,NX)=TRBH2(L,NY,NX)-DXOH2
      TRB2P(L,NY,NX)=TRB2P(L,NY,NX)-DXH2P
      TRALPO(L,NY,NX)=TRALPO(L,NY,NX)+DPALP
      TRFEPO(L,NY,NX)=TRFEPO(L,NY,NX)+DPFEP
      TRCAPD(L,NY,NX)=TRCAPD(L,NY,NX)+DPCDP
      TRCAPH(L,NY,NX)=TRCAPH(L,NY,NX)+DPCHP
      TRCAPM(L,NY,NX)=TRCAPM(L,NY,NX)+DPCMP
      TRALPB(L,NY,NX)=TRALPB(L,NY,NX)-DPALP
      TRFEPB(L,NY,NX)=TRFEPB(L,NY,NX)-DPFEP
      TRCPDB(L,NY,NX)=TRCPDB(L,NY,NX)-DPCDP
      TRCPHB(L,NY,NX)=TRCPHB(L,NY,NX)-DPCHP
      TRCPMB(L,NY,NX)=TRCPMB(L,NY,NX)-DPCMP
      ENDIF
      ELSE
C
C     AMALGAMATE PO4 BAND WITH NON-BAND IF BAND NO LONGER EXISTS 
C
      DPPOB(L,NY,NX)=0.0
      WDPOB(L,NY,NX)=0.0
      VLPOB(L,NY,NX)=0.0
      VLPO4(L,NY,NX)=1.0
      H0PO4(L,NY,NX)=H0PO4(L,NY,NX)+H0POB(L,NY,NX)
      H1PO4(L,NY,NX)=H1PO4(L,NY,NX)+H1POB(L,NY,NX)
      H2PO4(L,NY,NX)=H2PO4(L,NY,NX)+H2POB(L,NY,NX)
      H3PO4(L,NY,NX)=H3PO4(L,NY,NX)+H3POB(L,NY,NX)
      ZFE1P(L,NY,NX)=ZFE1P(L,NY,NX)+ZFE1PB(L,NY,NX)
      ZFE2P(L,NY,NX)=ZFE2P(L,NY,NX)+ZFE2PB(L,NY,NX)
      ZCA0P(L,NY,NX)=ZCA0P(L,NY,NX)+ZCA0PB(L,NY,NX)
      ZCA1P(L,NY,NX)=ZCA1P(L,NY,NX)+ZCA1PB(L,NY,NX)
      ZCA2P(L,NY,NX)=ZCA2P(L,NY,NX)+ZCA2PB(L,NY,NX)
      ZMG1P(L,NY,NX)=ZMG1P(L,NY,NX)+ZMG1PB(L,NY,NX)
      H0POB(L,NY,NX)=0.0
      H1POB(L,NY,NX)=0.0
      H2POB(L,NY,NX)=0.0
      H3POB(L,NY,NX)=0.0
      ZFE1PB(L,NY,NX)=0.0
      ZFE2PB(L,NY,NX)=0.0
      ZCA0PB(L,NY,NX)=0.0
      ZCA1PB(L,NY,NX)=0.0
      ZCA2PB(L,NY,NX)=0.0
      ZMG1PB(L,NY,NX)=0.0
      XOH0(L,NY,NX)=XOH0(L,NY,NX)+XOH0B(L,NY,NX)
      XOH1(L,NY,NX)=XOH1(L,NY,NX)+XOH1B(L,NY,NX)
      XOH2(L,NY,NX)=XOH2(L,NY,NX)+XOH2B(L,NY,NX)
      XH1P(L,NY,NX)=XH1P(L,NY,NX)+XH1PB(L,NY,NX)
      XH2P(L,NY,NX)=XH2P(L,NY,NX)+XH2PB(L,NY,NX)
      XOH0B(L,NY,NX)=0.0
      XOH1B(L,NY,NX)=0.0
      XOH2B(L,NY,NX)=0.0
      XH1PB(L,NY,NX)=0.0
      XH2PB(L,NY,NX)=0.0
      PALPO(L,NY,NX)=PALPO(L,NY,NX)+PALPB(L,NY,NX)
      PFEPO(L,NY,NX)=PFEPO(L,NY,NX)+PFEPB(L,NY,NX)
      PCAPD(L,NY,NX)=PCAPD(L,NY,NX)+PCPDB(L,NY,NX)
      PCAPH(L,NY,NX)=PCAPH(L,NY,NX)+PCPHB(L,NY,NX)
      PCAPM(L,NY,NX)=PCAPM(L,NY,NX)+PCPMB(L,NY,NX)
      PALPB(L,NY,NX)=0.0
      PFEPB(L,NY,NX)=0.0
      PCPDB(L,NY,NX)=0.0
      PCPHB(L,NY,NX)=0.0
      PCPMB(L,NY,NX)=0.0
      ENDIF
      ENDIF
C     ENDIF
C
C     SUBTRACT FERTILIZER DISSOLUTION FROM FERTILIZER POOLS
C
C     ZNH4FA,ZNH3FA,ZNHUFA,ZNO3FA=broadcast NH4,NH3,urea,NO3 fertilizer
C     ZNH4FB,ZNH3FB,ZNHUFB,ZNO3FB=banded NH4,NH3,urea,NO3 fertilizer
C     RSN4AA,RSN4BA=rate of broadcast NH4 fertilizer dissoln in non-band,band
C     RSN3AA,RSN3BA=rate of broadcast NH3 fertilizer dissoln in non-band,band
C     RSNUAA,RSNUBA=rate of broadcast urea fertr dissoln in non-band,band
C     RSNOAA,RSNOBA=rate of broadcast NO3 fertilizer dissoln in non-band,band
C     RSN4BB=rate of banded NH4 fertilizer dissolution in band
C     RSN3BB=rate of banded NH3 fertilizer dissolution in band
C     RSNUBB=rate of banded urea fertilizer dissolution in band
C     RSNOBB=rate of banded NO3 fertilizer dissolution in band
C
      ZNH4FA(L,NY,NX)=ZNH4FA(L,NY,NX)-RSN4AA-RSN4BA
      ZNH3FA(L,NY,NX)=ZNH3FA(L,NY,NX)-RSN3AA-RSN3BA
      ZNHUFA(L,NY,NX)=ZNHUFA(L,NY,NX)-RSNUAA-RSNUBA
      ZNO3FA(L,NY,NX)=ZNO3FA(L,NY,NX)-RSNOAA-RSNOBA
      ZNH4FB(L,NY,NX)=ZNH4FB(L,NY,NX)-RSN4BB
      ZNH3FB(L,NY,NX)=ZNH3FB(L,NY,NX)-RSN3BB
      ZNHUFB(L,NY,NX)=ZNHUFB(L,NY,NX)-RSNUBB
      ZNO3FB(L,NY,NX)=ZNO3FB(L,NY,NX)-RSNOBB
C
C     ADD FERTILIZER DISSOLUTION TO ION FLUXES AND CONVERT TO MASS
C
C     TRN3G=NH3 dissolution from NH3 
C     TRN4S,TRN4B=NH4 dissolution in non-band,band
C     TRN3S,TRN3B=NH3 dissolution from urea in non-band,band
C     TRNO3,TRNOB=NO3 dissolution in non-band,band
C
      TRN3G(L,NY,NX)=TRN3G(L,NY,NX)+RSN3AA+RSN3BA+RSN3BB 
      TRN4S(L,NY,NX)=TRN4S(L,NY,NX)+RSN4AA 
      TRN4B(L,NY,NX)=TRN4B(L,NY,NX)+RSN4BA+RSN4BB 
      TRN3S(L,NY,NX)=TRN3S(L,NY,NX)+RSNUAA
      TRN3B(L,NY,NX)=TRN3B(L,NY,NX)+RSNUBA+RSNUBB
      TRNO3(L,NY,NX)=TRNO3(L,NY,NX)+RSNOAA
      TRNOB(L,NY,NX)=TRNOB(L,NY,NX)+RSNOBA+RSNOBB
      TRCO2(L,NY,NX)=TRCO2(L,NY,NX)*12.0
      TRN3G(L,NY,NX)=TRN3G(L,NY,NX)*14.0
      TRN4S(L,NY,NX)=TRN4S(L,NY,NX)*14.0
      TRN4B(L,NY,NX)=TRN4B(L,NY,NX)*14.0
      TRN3S(L,NY,NX)=TRN3S(L,NY,NX)*14.0
      TRN3B(L,NY,NX)=TRN3B(L,NY,NX)*14.0
      TRNO3(L,NY,NX)=TRNO3(L,NY,NX)*14.0
      TRNOB(L,NY,NX)=TRNOB(L,NY,NX)*14.0
      TRNO2(L,NY,NX)=TRNO2(L,NY,NX)*14.0
      TRN2B(L,NY,NX)=TRN2B(L,NY,NX)*14.0
      TRH1P(L,NY,NX)=TRH1P(L,NY,NX)*31.0
      TRH2P(L,NY,NX)=TRH2P(L,NY,NX)*31.0
      TRH1B(L,NY,NX)=TRH1B(L,NY,NX)*31.0
      TRH2B(L,NY,NX)=TRH2B(L,NY,NX)*31.0
      ENDIF
C     IF(I.EQ.116)THEN
C     WRITE(*,9984)'TRN3S',I,J,L,TRN4S(L,NY,NX),TRN3S(L,NY,NX)
9984  FORMAT(A8,3I4,20F14.7)
C     ENDIF
9985  CONTINUE 
C
C     SURFACE RESIDUE
C
C     BKVL=litter mass
C
      IF(VOLWM(NPH,0,NY,NX).GT.ZEROS2(NY,NX))THEN
      BKVLX=BKVL(0,NY,NX)
C
C     UREA HYDROLYSIS IN SURFACE RESIDUE
C
C     VOLQ=biologically active litter water volume from nitro.f
C     COMA=concentration of active biomass
C     TOQCK=total microbial activity from nitro.f
C     DUKD=Km for urea hydrolysis
C
      IF(VOLQ(0,NY,NX).GT.ZEROS2(NY,NX))THEN
      COMA=AMIN1(0.1E+06,TOQCK(0,NY,NX)/VOLQ(0,NY,NX))
      ELSE
      COMA=0.1E+06
      ENDIF
      DUKD=DUKM*(1.0+COMA/DUKI)
C
C     UREA HYDROLYSIS INHIBITION
C
C     ZNHU0,ZNHUI=initial,current inhibition activity
C     RNHUI=rate constant for decline in urea hydrolysis inhibition
C
      IF(ZNHU0(0,NY,NX).GT.ZEROS(NY,NX)
     2.AND.ZNHUI(0,NY,NX).GT.ZEROS(NY,NX))THEN
      ZNHUI(0,NY,NX)=ZNHUI(0,NY,NX)-TFNQ(0,NY,NX)**0.25 
     2*RNHUI(IUTYP(NY,NX))*ZNHUI(0,NY,NX)
     3*AMAX1(RNHUI(IUTYP(NY,NX)),1.0-ZNHUI(0,NY,NX)/ZNHU0(0,NY,NX))
      ELSE
      ZNHUI(0,NY,NX)=0.0
      ENDIF
C
C     UREA CONCENTRATION AND HYDROLYSIS IN SURFACE RESIDUE
C
C     ZNHUFA=urea fertilizer 
C     CNHUA=concentration of urea fertilizer 
C     DFNSA=effect of microbial concentration on urea hydrolysis 
C     RSNUA=rate of urea hydrolysis 
C     SPNHU=specific rate constant for urea hydrolysis
C     TFNQ=temperature effect on microbial activity from nitro.f
C
      IF(ZNHUFA(0,NY,NX).GT.ZEROS(NY,NX)
     2.AND.BKVL(0,NY,NX).GT.ZEROS(NY,NX))THEN
      CNHUA=ZNHUFA(0,NY,NX)/BKVL(0,NY,NX)
      DFNSA=CNHUA/(CNHUA+DUKD)
      RSNUA=AMIN1(ZNHUFA(0,NY,NX)
     2,SPNHU*TOQCK(0,NY,NX)*DFNSA*TFNQ(0,NY,NX))*(1.0-ZNHUI(0,NY,NX))
      ELSE
      RSNUA=0.0
      ENDIF
C     IF(ZNHUFA(0,NY,NX).GT.ZEROS(NY,NX))THEN
C     WRITE(*,8778)'UREA0',I,J,IUTYP(NY,NX)
C    2,ZNHUFA(0,NY,NX),RSNUA
C    2,DFNSA,TFNQ(0,NY,NX),CNHUA,DUKD,DUKM,DUKI,TOQCK(0,NY,NX)
C    3,BKVL(0,NY,NX),TFNQ(0,NY,NX),SPNHU,ZNHU0(0,NY,NX),ZNHUI(0,NY,NX)
C    4,RNHUI(IUTYP(NY,NX))
8778  FORMAT(A8,3I4,40E12.4)
C     ENDIF
C
C     NH4, NH3, UREA, NO3 DISSOLUTION IN SURFACE RESIDUE
C     FROM FIRST-ORDER FUNCTIONS OF REMAINING
C     FERTILIZER (NOTE: SUPERPHOSPHATE AND ROCK PHOSPHATE
C     ARE REPRESENTED AS MONOCALCIUM PHOSPHATE AND HYDROXYAPATITE
C     MODELLED IN PHOSPHORUS REACTIONS BELOW)
C
C     RSN4AA=rate of broadcast NH4 fertilizer dissoln 
C     RSN3AA=rate of broadcast NH3 fertilizer dissoln 
C     RSNUAA=rate of broadcast urea fertr dissoln 
C     RSNOAA=rate of broadcast NO3 fertilizer dissoln 
C
      IF(VOLWRX(NY,NX).GT.ZEROS(NY,NX))THEN
      THETWR=AMIN1(1.0,VOLW(0,NY,NX)/VOLWRX(NY,NX))
      ELSE
      THETWR=1.0
      ENDIF
      RSN4AA=SPNH4*ZNH4FA(0,NY,NX)*THETWR
      RSN3AA=SPNH3*ZNH3FA(0,NY,NX)
      RSNUAA=RSNUA*THETWR
      RSNOAA=SPNO3*ZNO3FA(0,NY,NX)*THETWR
C
C     SOLUBLE AND EXCHANGEABLE NH4 CONCENTRATIONS
C
C     VOLWNH=water volume 
C     RN4X,RN3X=NH4,NH3 input from uptake, mineraln, dissoln 
C     XNH4S=net change in NH4 from nitro.f
C     CN41,CN31=total NH4,NH3 concentration 
C     XN41=adsorbed NH4 concentration 
C
      IF(VOLWM(NPH,0,NY,NX).GT.ZEROS2(NY,NX))THEN
      VOLWMX=14.0*VOLWM(NPH,0,NY,NX)
      RN4X=(XNH4S(0,NY,NX)+14.0*RSN4AA)/VOLWMX
      RN3X=14.0*RSNUAA/VOLWMX
      CN41=AMAX1(ZERO,ZNH4S(0,NY,NX)/VOLWMX+RN4X)
      CN31=AMAX1(ZERO,ZNH3S(0,NY,NX)/VOLWMX+RN3X)
      IF(BKVLX.GT.ZEROS(NY,NX))THEN
      XN41=AMAX1(ZERO,XN4(0,NY,NX)/BKVLX)
      ELSE
      XN41=0.0
      ENDIF
C
C     SOLUBLE, EXCHANGEABLE AND PRECIPITATED PO4 CONCENTRATIONS 
C
C     VOLWMP=water volume 
C     RH1PX,RH2PX=HPO4,H2PO4 inputs from mineraln, uptake 
C     XH1PS=net change in HPO4 from nitro.f 
C     CH1P1,CH2P1=HPO4,H2PO4 concentrations 
C
      VOLWMP=31.0*VOLWM(NPH,0,NY,NX)
      RH1PX=XH1PS(0,NY,NX)/VOLWMP
      RH2PX=XH2PS(0,NY,NX)/VOLWMP
      CH1P1=AMAX1(0.0,H1PO4(0,NY,NX)/VOLWMP+RH1PX)
      CH2P1=AMAX1(0.0,H2PO4(0,NY,NX)/VOLWMP+RH2PX)
      ELSE
      RN4X=0.0
      RN3X=0.0
      CN41=0.0
      CN31=0.0
      XN41=0.0
      RH1PX=0.0
      RH2PX=0.0
      CH1P1=0.0
      CH2P1=0.0
      ENDIF
C
C     PHOSPHORUS TRANSFORMATIONS IN SURFACE RESIDUE
C
C     PALPO1,PFEPO1=concn of precip AlPO4,FEPO4 
C     PCAPM1,PCAPD1,PCAPH1=concn of precip CaH2PO4,CaHPO4,apatite 
C
      IF(BKVLX.GT.ZEROS(NY,NX))THEN
      PALPO1=AMAX1(0.0,PALPO(0,NY,NX)/BKVLX)
      PFEPO1=AMAX1(0.0,PFEPO(0,NY,NX)/BKVLX)
      PCAPM1=AMAX1(0.0,PCAPM(0,NY,NX)/BKVLX)
      PCAPD1=AMAX1(0.0,PCAPD(0,NY,NX)/BKVLX)
      PCAPH1=AMAX1(0.0,PCAPH(0,NY,NX)/BKVLX)
      ELSE
      PALPO1=0.0
      PFEPO1=0.0
      PCAPM1=0.0
      PCAPD1=0.0
      PCAPH1=0.0
      ENDIF
C
C     CALCULATE H2PO4 COPRECIPITATES FRPM LITTER PH
C
      CHY1=AMAX1(ZERO,10.0**(-(PH(0,NY,NX)-3.0)))
      COH1=AMAX1(ZERO,DPH2O/CHY1)
      CAL1=AMAX1(ZERO,SPALO/COH1**3)
      CFE1=AMAX1(ZERO,SPFEO/COH1**3)
      CCO20=AMAX1(ZERO,CCO2S(0,NY,NX)/12.0)
      CCO31=AMAX1(ZERO,CCO20*DPCO3/CHY1**2)
      CCA1=AMAX1(ZERO,AMIN1(CCAMX,SPCAC/CCO31))
C
C     ALUMINUM PHOSPHATE (VARISCITE)
C
C     CH2PA,CH2P1=equilibrium,current H2PO4 concentration in litter 
C     SYA0P2=solubility product derived from SPALO
C     RPALPX=H2PO4 dissolution from AlPO4 in litter
C
      CH2PA=SYA0P2/(CAL1*COH1**2)
      RPALPX=AMIN1(AMAX1(0.0,4.0E-08*ORGC(0,NY,NX)-PALPO1)
     2,AMAX1(-PALPO1,TPD*(CH2P1-CH2PA)))
C
C     IRON PHOSPHATE (STRENGITE)
C
C     CH2PF,CH2P1=equilibrium,current H2PO4 concentration in litter 
C     SYF0P2=solubility product derived from SPALO
C     RPFEPX=H2PO4 dissolution from FePO4 in litter
C
      CH2PF=SYF0P2/(CFE1*COH1**2)
      RPFEPX=AMIN1(AMAX1(0.0,2.0E-06*ORGC(0,NY,NX)-PFEPO1)
     2,AMAX1(-PFEPO1,TPD*(CH2P1-CH2PF)))
C
C     DICALCIUM PHOSPHATE
C
C     CH2PD,CH2P1=equilibrium,current H2PO4 concentration in litter
C     SYCAD2=solubility product derived from SPALO
C     RPCADX=H2PO4 dissolution from CaHPO4 in litter
C
      CH2PD=SYCAD2/(CCA1*COH1)
      RPCADX=AMIN1(AMAX1(0-.0,5.0E-05*ORGC(0,NY,NX)-PCAPD1)
     2,AMAX1(-PCAPD1,TPD*(CH2P1-CH2PD)))
C
C     HYDROXYAPATITE
C
C     CH2PH,CH2P1=equilibrium,current H2PO4 concentration in litter
C     SYCAH2=solubility product derived from SPALO
C     RPCAHX=H2PO4 dissolution from apatite in litter
C
      CH2PH=(SYCAH2/(CCA1**5*COH1**7))**0.333
      RPCAHX=AMIN1(AMAX1(0.0,5.0E-05*ORGC(0,NY,NX)-PCAPH1)
     2,AMAX1(-PCAPH1,TPD*(CH2P1-CH2PH)))
C
C     MONOCALCIUM PHOSPHATE
C
C     CH2PM,CH2P1=equilibrium,current H2PO4 concentration in litter
C     SPCAM=solubility product for Ca(H2PO4)2 
C     RPCAMX=H2PO4 dissolution from Ca(H2PO4)2 in litter
C
      CH2PM=SQRT(SPCAM/CCA1)
      RPCAMX=AMIN1(AMAX1(0.0,5.0E-05*ORGC(0,NY,NX)-PCAPM1)
     2,AMAX1(-PCAPM1*SPPO4,TPD*(CH2P1-CH2PM)))
C     IF(I.GT.315)THEN
C     WRITE(*,2227)'RPPO4',I,J,L,RPCAHX,CH2P1,CH2PA,CH2PH
C    2,SYA0P2,CAL1,COH1,SYCAH2,CCA1,CCO21,CCO31,PCAPH1
C    3,VOLWM(NPH,0,NY,NX),SPCAC/CCO31,H2PO4(0,NY,NX)
C    4,CCO20,DPCO3,CHY1,CCO2S(0,NY,NX)
2227  FORMAT(A8,3I4,20E12.4)
C     ENDIF
C
C     PHOSPHORUS ANION EXCHANGE IN SURFACE REDISUE
C     CALCULATED FROM EXCHANGE EQUILIBRIA AMONG H2PO4-,
C     HPO4--, H+, OH- AND PROTONATED AND NON-PROTONATED -OH
C     EXCHANGE SITES (NOT CALCULATED)
C
C     H2PO4-H+HPO4
C
C     DPH2P=dissociation constant
C     S1=equilibrium concentration in litter
C     RH2P=H2PO4-H+HPO4 dissociation in litter 
C     
      DP=DPH2P 
      S0=CH1P1+CHY1+DP
      S1=AMAX1(0.0,S0**2-4.0*(CH1P1*CHY1-DP*CH2P1))
      RH2P=TSL*(S0-SQRT(S1))
C
C     EQUILIBRIUM X-CA CONCENTRATION FROM CEC, GAPON COEFFICIENTS 
C     AND CATION CONCENTRATIONS
C
C     CCEC,XCEC=cation exchange concentration,capacity
C     XCAX=equilibrium R-Ca concentration
C     GKC4,GKCH,GKCA,GKCM,GKCN,GKCK=Gapon selectivity coefficients for
C     CA-NH4,CA-H,CA-AL,CA-MG,CA-NA,CA-K
C     X*Q=equilibrium exchangeable concentrations
C     XTLQ=total equilibrium exchangeable concentration 
C
      IF(BKVLX.GT.ZEROS(NY,NX))THEN
      CCEC0=AMAX1(ZERO,COOH*ORGC(0,NY,NX)/BKVLX)
      ELSE
      CCEC0=ZERO
      ENDIF
      CALX=AMAX1(ZERO,CAL1)**0.333
      CFEX=AMAX1(ZERO,CFE1)**0.333
      CCAX=AMAX1(ZERO,CCA1)**0.500
C
C     EQUILIBRIUM X-CA CONCENTRATION FROM CEC AND CATION
C     CONCENTRATIONS
C
      XCAX=CCEC0/(1.0+GKC4(NU(NY,NX),NY,NX)*CN41/CCAX
     2+GKCH(NU(NY,NX),NY,NX)*CHY1/CCAX
     3+GKCA(NU(NY,NX),NY,NX)*CALX/CCAX
     3+GKCA(NU(NY,NX),NY,NX)*CFEX/CCAX)
      XN4Q=XCAX*CN41*GKC4(NU(NY,NX),NY,NX)
      XHYQ=XCAX*CHY1*GKCH(NU(NY,NX),NY,NX)
      XALQ=XCAX*CALX*GKCA(NU(NY,NX),NY,NX)
      XFEQ=XCAX*CFEX*GKCA(NU(NY,NX),NY,NX)
      XCAQ=XCAX*CCAX
      XTLQ=XN4Q+XHYQ+XALQ+XFEQ+XCAQ 
      IF(XTLQ.GT.ZERO)THEN
      FX=CCEC0/XTLQ
      ELSE
      FX=0.0
      ENDIF
      XN4Q=FX*XN4Q
C
C     NH4 AND NH3 EXCHANGE IN SURFACE RESIDUE
C
C     RXN4=NH4 adsorption in litter
C     TADC0=adsorption rate constant
C     RNH4=NH4-NH3+H dissociation in litter
C     DPN4=NH4 dissociation constant
C
      RXN4=TADC0*(XN4Q-XN41) 
      RNH4=(CHY1*CN31-DPN4*CN41)/(DPN4+CHY1)
C     IF(J.EQ.12)THEN
C     WRITE(*,2223)'RXN4',I,J,RXN4,CN41,XN41,TADC0,XN4Q
C    2,CCAX,CCA1,CCO20,CCO31 
C    2,XCAQ,CCEC0,FN4X,FCAQ,GKC4(NU(NY,NX),NY,NX)
C    3,PH(0,NY,NX),CHY1,RNH4
C    3,CN31,DPN4,ZNH4S(0,NY,NX),XN4(0,NY,NX),14.0*RSN4AA,RN4X,BKVLX
C    4,BKVL(0,NY,NX),VOLWM(NPH,0,NY,NX) 
2223  FORMAT(A8,2I4,30E12.4)
C     ENDIF
      ELSE
      RSN4AA=0.0
      RSN3AA=0.0
      RSNUAA=0.0
      RSNOAA=0.0
      RPALPX=0.0
      RPFEPX=0.0
      RPCADX=0.0
      RPCAHX=0.0
      RPCAMX=0.0
      RXN4=0.0
      RNH4=0.0
      RH2P=0.0
      RPALPX=0.0
      RPFEPX=0.0
      RPCADX=0.0
      RPCAMX=0.0
      RPCAHX=0.0 
      ENDIF
C
C     TOTAL ION FLUXES FOR ALL REACTIONS ABOVE
C
C     RN4S=net NH4 flux in litter
C     RN3S=net NH3 flux in litter
C     RHP1,RHP2=net HPO4,H2PO4 flux in litter
C
      RN4S=RNH4-RXN4
      RN3S=-RNH4
      RHP1=-RH2P
      RHP2=RH2P-RPALPX-RPFEPX-RPCADX-2.0*RPCAMX-3.0*RPCAHX
C
C     CONVERT TOTAL ION FLUXES FROM CHANGES IN CONCENTRATION
C     TO CHANGES IN MASS PER UNIT AREA FOR USE IN 'REDIST'
C
C     TRN4S=total NH4 flux 
C     TRN3S=total NH3 flux 
C     TRH1P,TRH2P=net HPO4,H2PO4 flux 
C     TRNX4=total NH4 adsorption 
C     TRALPO,TRFEPO,TRCAPD,TRCAPH,TRCAPM
C     =total AlPO4,FePO4,CaHPO4,apatite,Ca(H2PO4)2 precipitation 
C
      TRN4S(0,NY,NX)=TRN4S(0,NY,NX)+RN4S*VOLWM(NPH,0,NY,NX)
      TRN3S(0,NY,NX)=TRN3S(0,NY,NX)+RN3S*VOLWM(NPH,0,NY,NX)
      TRH1P(0,NY,NX)=TRH1P(0,NY,NX)+RHP1*VOLWM(NPH,0,NY,NX)
      TRH2P(0,NY,NX)=TRH2P(0,NY,NX)+RHP2*VOLWM(NPH,0,NY,NX)
      TRXN4(0,NY,NX)=TRXN4(0,NY,NX)+RXN4*VOLWM(NPH,0,NY,NX)
      TRALPO(0,NY,NX)=TRALPO(0,NY,NX)+RPALPX*VOLWM(NPH,0,NY,NX)
      TRFEPO(0,NY,NX)=TRFEPO(0,NY,NX)+RPFEPX*VOLWM(NPH,0,NY,NX)
      TRCAPD(0,NY,NX)=TRCAPD(0,NY,NX)+RPCADX*VOLWM(NPH,0,NY,NX)
      TRCAPH(0,NY,NX)=TRCAPH(0,NY,NX)+RPCAHX*VOLWM(NPH,0,NY,NX)
      TRCAPM(0,NY,NX)=TRCAPM(0,NY,NX)+RPCAMX*VOLWM(NPH,0,NY,NX)
      ZNH4FA(0,NY,NX)=ZNH4FA(0,NY,NX)-RSN4AA
      ZNH3FA(0,NY,NX)=ZNH3FA(0,NY,NX)-RSN3AA
      ZNHUFA(0,NY,NX)=ZNHUFA(0,NY,NX)-RSNUAA
      ZNO3FA(0,NY,NX)=ZNO3FA(0,NY,NX)-RSNOAA
      TRN4S(0,NY,NX)=TRN4S(0,NY,NX)+RSN4AA 
      TRN3S(0,NY,NX)=TRN3S(0,NY,NX)+RSN3AA+RSNUAA 
      TRNO3(0,NY,NX)=TRNO3(0,NY,NX)+RSNOAA
      TRN4S(0,NY,NX)=TRN4S(0,NY,NX)*14.0
      TRN3S(0,NY,NX)=TRN3S(0,NY,NX)*14.0
      TRNO3(0,NY,NX)=TRNO3(0,NY,NX)*14.0
      TRH1P(0,NY,NX)=TRH1P(0,NY,NX)*31.0
      TRH2P(0,NY,NX)=TRH2P(0,NY,NX)*31.0
C     WRITE(*,9989)'TRN4S',I,J,TRN4S(0,NY,NX)
C    2,RN4S,RNH4,RXN4,RSN4AA,VOLWM(NPH,0,NY,NX)
C    3,SPNH4,ZNH4FA(0,NY,NX),THETWR
9989  FORMAT(A8,2I4,12E12.4)
9990  CONTINUE
9995  CONTINUE
      RETURN
      END


