c
c   THIS PROGRAM SOLVES A SYSTEM OF DIFFERENTIAL EQUATIONS WITH BOUNDARY
c   CONDITONS AT TWO POINTS USING THE FINITE DIFFERENCE METHOD
c   (IMSL SUBROUTINE DBVPFD)
c
        program ode_O2
c                                                                               
	integer neq, mxgrid, ninit, ldyfin, ldyini, maxKz
c
        parameter (neq = 2)		! Number (Order) of diff eqns
        parameter (mxgrid = 10001)	! Max number grid pts
        parameter (ninit = mxgrid)	! Number of initial grid points
	parameter (maxKz = 10)		! Max number of eddy diffusivity constraints
        parameter (ldyfin = neq, ldyini = neq)
c
        real*8 x_bc1, x_bc2		! Boundaries conditions
        real*8 C_bc1, C_bc2 
	real*8 x_Res			! Resolution transition depth
c
	real*8 xinit(mxgrid),Cinit(neq,mxgrid)		! Initial grid (input)
	real*8 xfinal(mxgrid),Cfinal(neq,mxgrid)	! Final grid (output)
c
	real*8 T, S, DEPTH, P		! Temperature, salinity, wc depth, pressure
	real*8 p_w, p_sw		! Density of pure water, seawater 
	real*8 MU, V			! Dynamic viscosity, kinematic viscosity
	real*8 C_d			! Drag coefficient
	real*8 U_100			! Current velocity 100 cm above bottom
	real*8 U_star			! Friction velocity
	real*8 Sc			! Schmidt number
	real*8 LBL, VBL			! Logarithmic, viscous, and diffusive boundary layers
	real*8 eDBL			! Effective diffusive boundary layer
	real*8 DBL			! Diffusive boundary layer (where Kz = D)
	real*8 TLP, TL, Q1, Q2		! Transition layer parameter
	real*8 sigma_DBL, D_patch
c
	real*8 D			! O2 molecular diffusivity
	real*8 E			! Eddy diffusivity
	real*8 O_thresh			! Threshold conc. for O2 limitation
	real*8 O_TGP			! Threshold gradient parameter for O2
	real*8 R			! Zero-order O2 consumption rate constant
	real*8 REF			! Rate enhancement factor (linear function of depth)
	real*8 IGP			! Interface gradient parameter
c
	real*8 por_0, porbeta		! Porosity
	real*8 porinf, PGP		!   parameters
c
	real*8 O2(mxgrid)		! O2 concentration
	real*8 break_O2(mxgrid)		! Spline
	real*8 cscoef_O2(4,mxgrid)	!   parameters
	real*8 O2_DBL			! O2 concentration at top of DBL
	real*8 O2_lin			! O2 concentration in linear region of DBL
c
	real*8 dO2_dx(mxgrid)		! O2 gradient
	real*8 break_dO2_dx(mxgrid)	! Spline
	real*8 cscoef_dO2_dx(4,mxgrid)	!   parameters
	real*8 dO2_dx_bc1		! O2 gradient at x_bc1
	real*8 dO2_dx_DBL		! O2 gradient at top of DBL
	real*8 dO2_dx_eDBL		! O2 gradient at effective DBL
	real*8 dO2_dx_if		! O2 gradient at sed-water interface
	real*8 dO2_dx_lin		! O2 gradient in linear region of DBL
c
	real*8 POR(mxgrid)		! Porosity
	real*8 break_POR(mxgrid)	! Spline
	real*8 cscoef_POR(4,mxgrid)	!   parameters
	real*8 POR_bc1			! Porosity at x_bc1
	real*8 POR_DBL			! Porosity at top of DBL
	real*8 POR_eDBL			! Porosity at effective DBL
	real*8 POR_if			! Porosity at sed-water interface
c
	real*8 Ds(mxgrid)		! Sediment diffusivity
	real*8 break_Ds(mxgrid)		! Spline
	real*8 cscoef_Ds(4,mxgrid)	!   parameters
	real*8 Ds_bc1			! Sediment diffusivity at x_bc1
	real*8 Ds_DBL			! Sediment diffusivity at top of DBL
	real*8 Ds_eDBL			! Sediment diffusivity at effective DBL
	real*8 Ds_if			! Sediment diffusivity at sed-water interface
c
	real*8 x_Kz(maxKz), Kz(maxKz)	! Eddy diffusivity constraints
	real*8 break_Kz(maxKz)		! Spline
	real*8 cscoef_Kz(4,maxKz)	!   parameters
	real*8 EDDY(mxgrid)		! Eddy diffusivity
	real*8 Kz_bc1			! Eddy diffusivity at x_bc1
	real*8 Kz_DBL			! Eddy diffusivity at top of DBL
	real*8 Kz_eDBL			! Eddy diffusivity at effective DBL
	real*8 Kz_if			! Eddy diffusivity at sed-water interface
c
	real*8 FLUX_bc1			! O2 flux at x_bc1
	real*8 FLUX_DBL			! O2 flux at top of DBL
	real*8 FLUX_eDBL		! O2 flux at effective DBL
	real*8 FLUX_if			! O2 flux at sed-water interface
c
	real*8 OCR(mxgrid)		! O2 consumption rate
	real*8 OCR_ws(mxgrid)		! O2 consumption rate
	real*8 break_OCR_ws(mxgrid)	! Spline
	real*8 cscoef_OCR_ws(4,mxgrid)	!   parameters
	real*8 OCR_ws_itg		! Integrated OCR
c
	real*8 tol			! Relative error control parameter
        real*8 errest(neq)		! Error estimate
        real*8 pistep			! Normally set to zero
c
	real*8 xplot
c
        integer nfinal			! Number of final grid points 
	integer nleft			! Number of initial conditions
	integer ncupbc			! Number of coupled boundary conditions
	integer nKz			! Number of eddy diffusivity constraints

	integer nskip, i_toggle
c
	logical linear, print
c                                                                               
	common/intsed/nKz
	common/realsed/D,V,O_thresh,O_TGP,R,REF,IGP,DBL,TLP,sigma_DBL,D_patch,break_Kz,cscoef_Kz
	common/realpor/por_0,porbeta,porinf,PGP
	common/realbound/C_bc1,C_bc2
c
 	external dbvpfd
	real*8 fcneqn, fcnjac, fcnbc
        external fcneqn, fcnjac, fcnbc
	real*8 derf
	external derf
	real*8 dcsval, dcsitg
	external dcsint, dcsval, dcsitg, dcsakm     
	real*8 dabs, dsqrt
	intrinsic dabs, dsqrt
	integer int
	intrinsic int
c
	print 5
    5	format(//,10x,'***** Program "ode_O2" *****',/)
c
c BVPFD PARAMETERS 
c
	tol = 0.01
	linear = .false.
	print = .false.
        pistep = 0.0
c
c TYPES OF BOUNDARY CONDITIONS
c
        nleft = 1  
        ncupbc = 0
c
c INPUT BOUNDARY CONDITIONS
c
	x_bc1 = -1.0		! cm (0 = sediment-water interface)
	x_bc2 = 2.0
	x_Res = 0.5		! set to AVG(x_bc1,x_bc2) for evenly spaced grid
c
	print 7
    7	format(' Enter Bottom Water O2 concentration (uM): ',$)
	read *, C_bc1
	C_bc1 = C_bc1/1000	! mM
	i_toggle = 2		! 1 = C_bc1 measured 10 cm above bottom (Hydrolab)
				! 2 = C_bc1 measured 2 mm above bottom (microelectrode)
	if(i_toggle.eq.1) then
		C_bc1 = C_bc1*0.975	!([x_bc1 = -1.0] corrects for dC/dx from 1-10 cm) 
	else
		C_bc1 = C_bc1/0.985	!([x_bc1 = -1.0] corrects for dC/dx from 0.2-1 cm) 
	end if	
	C_bc2 = 0.0		! mM
c
c INPUT HYDRODYNAMIC PARAMETERS
c
	print 10
   10	format(' Enter Bottom Current (cm/s): ',$)
	read *, U_100
c
	T = 27.0		! C
	S = 11.0		! psu
	DEPTH = 0.0		! m
c
	P = 1. + DEPTH/10.	! bar
c
c INPUT MOLECULAR AND EDDY DIFFUSIVITIES
c
	D = 2.30e-5 * (60*60*24)	! cm2/d (for O2 at 0 psu, 25C) 
	E = 1.0     * (60*60*24)	! 1.0 cm2/d
c
c CALCULATE DYNAMIC (MOLECULAR) VISCOSITY
c
	MU = (1.7910 - (6.144e-2)*T + (1.4510e-3)*T**2 - (1.6826e-5)*T**3
     &		- (1.5290e-4)*P + (8.3885e-8)*P**2 + (2.4727e-3)*S
     &		+ T*(6.0574e-6*P - 2.670e-9*(P**2)) + S*(4.8429e-5*T
     &		- (4.7172e-6)*T**2 + (7.5986e-8)*T**3))*0.01			! g/cm s
c
c CORRECT MOLECULAR DIFFUSIVITY FOR T AND S
c
	D = D * (0.0089897/MU) * ((T+273.15)/298.15)
c
c EQUATION OF STATE FOR SEAWATER
c
	p_w = 999.842594 + 6.793952e-2*T - 9.095290e-3*T**2 + 1.001685e-4*T**3
     &		- 1.120083e-6*T**4 + 6.536332e-9*T**5				! kg/m3
	p_sw = (p_w + (8.24493e-1 - 4.0899e-3*T + 7.6438e-5*T**2
     &		- 8.2467e-7*T**3 + 5.3875e-9*T**4)*S + (-5.72466e-3
     &		+ 1.0227e-4*T -1.6546e-6*t**2)*S**1.5 + 4.8314e-4*S**2)/1000	! g/cm3
c
c CALCULATE KINEMATIC VISCOSITY, DRAG COEFFICIENT, FRICTION VELOCITY, & SCHMIDT NUMBER
c
	V = MU/p_sw * (60*60*24)					! cm2/d
	C_d = 1.e-3 * (2.33 - (0.0700*U_100) + (0.000646*U_100**2))	! unitless
	U_star = (dsqrt(C_d) * U_100) * (60*60*24)			! cm/d
	Sc = V/D							! unitless
c
c CALCULATE LOGARITHMIC, VISCOUS, AND DIFFUSIVE BOUNDARY LAYERS (cm)
c
	LBL = 100.					! Note that boundary layers are
	VBL = 10*V/U_star				! expressed as positive numbers
	eDBL = D/(0.0889*U_star*(Sc**(-0.704)))		! even though they represent 
							! distances above the sediment surface.
c
	print 15
   15	format(' Enter Q1 (TL = Q1*DBL): ',$)
   	read *, Q1
	Q2 = 0.5					! eDBL = DBL - Q2*TL (0.lt.Q2.lt.1)
	DBL = eDBL/(1-(Q1*Q2))				! normally Q2 = 0.5
	TL = Q1 * DBL
c
	TLP = TL/2	
	sigma_DBL = DBL/100.
	D_patch = 5*D
c
c INPUT KNOTS FOR EDDY DIFFUSION COEFFICIENT
c
	nKz = 5
c	
	x_Kz(1) = -2.00 * LBL
	x_Kz(2) = -1.75 * LBL
	x_Kz(3) = -1.50 * LBL
	x_Kz(4) = -VBL
	x_Kz(5) = -DBL
c
	Kz(1) = E
	Kz(2) = E
	Kz(3) = E
	Kz(4) = V
	Kz(5) = D_patch
c
	call dcsakm(nKz,x_Kz,Kz,break_Kz,cscoef_Kz)
c
c ENTER OXYGEN CONSUMPTION RATE AND FACTOR FOR LINEAR DEPTH DEPENDENCE
c
	print 20
   20	format(' Enter O2 Consumption Rate (mM/d): ',$)
	read *, R
  	print 25
   25	format(' Enter REF (OCR = R[1+REF*x]): ',$)
	read *, REF
	IGP = 0.001		! Interface gradient parameter (10 um)
c
c INPUT THRESHOLD AND GRADIENT PARAMETER VALUES
c
	O_thresh = 0.002	! Oxygen threshold (0.002 mM) [min. O_thresh=0.0002]
	O_TGP = O_thresh/2	! Oxygen threshold gradient parameter (mM)
c
c INPUT POROSITY PARAMETERS
c
c
c Porosity and solid matter density
c
	i_toggle = 1	! 1 = default values; 2 = interactive
	if(i_toggle.eq.1) then
		por_0 = 0.970		! por_0 = porinf for constant porosity
		porinf = 0.904
		porbeta = 0.742		! cm-1
	else
		print 26
		read *, por_0
		print 27
		read *, porinf
		print 28
		read *, porbeta
	end if 
	PGP = 0.001			! Porosity gradient parameter (10 um)
c
   26	format(' Enter porosity at sediment-water interface: ',$)
   27	format(' Enter porosity below zone of compaction: ',$)
   28	format(' Enter depth-attenuation coefficient for porosity (cm-1): ',$)
c
c SET UP INITIAL GRID
c
	jinit = int(ninit/2)
	kinit = ninit-jinit
	do i=1,ninit
		if(i.le.jinit) then
			xinit(i) = x_bc1 + (i-1)*(x_Res - x_bc1)/(jinit)
		else
			xinit(i) = x_Res + (i-kinit)*(x_bc2 - x_Res)/(kinit-1)
		end if
	end do
c
c  Initial guess of Cinit values
c
	do j=1,neq
		do i=1,ninit
			Cinit(j,i) = 0.0
		end do
	end do
c
c SOLVE ODE
c                                                                            
        call dbvpfd(fcneqn,fcnjac,fcnbc,fcneqn,fcnbc,neq,nleft,ncupbc,
     &         x_bc1,x_bc2,pistep,tol,ninit,xinit,cinit,ldyini,linear,
     &		print,mxgrid,nfinal,xfinal,cfinal,ldyfin,errest)
c
c CALCULATE O2 AND dO2/dx
c
	do i=1,nfinal
		O2(i) = Cfinal(1,i)
		dO2_dx(i) = Cfinal(2,i)
	end do
	call dcsint(nfinal,xfinal,O2,break_O2,cscoef_O2)
	O2_DBL = dcsval(-DBL,mxgrid-1,break_O2,cscoef_O2)
	O2_lin = dcsval(-(DBL-TL)/2,mxgrid-1,break_O2,cscoef_O2)
c
	call dcsint(nfinal,xfinal,dO2_dx,break_dO2_dx,cscoef_dO2_dx)
	dO2_dx_bc1 = dcsval(x_bc1,mxgrid-1,break_dO2_dx,cscoef_dO2_dx)
	dO2_dx_DBL = dcsval(-DBL,mxgrid-1,break_dO2_dx,cscoef_dO2_dx)
	dO2_dx_eDBL = dcsval(-eDBL,mxgrid-1,break_dO2_dx,cscoef_dO2_dx)
	dO2_dx_if = dcsval(0.,mxgrid-1,break_dO2_dx,cscoef_dO2_dx)
	dO2_dx_lin = dcsval(-(DBL-TL)/2,mxgrid-1,break_dO2_dx,cscoef_dO2_dx)
c
c CALCULATE POROSITY
c
	do i=1,nfinal
		POR(i) = 0.5*(1-derf(xfinal(i)/PGP)) 
     &			+ 0.5*(1+derf(xfinal(i)/PGP))
     &			*((por_0-porinf)*dexp(-porbeta*xfinal(i))+porinf)
     	end do
	call dcsint(nfinal,xfinal,POR,break_POR,cscoef_POR)
	POR_bc1 = dcsval(x_bc1,mxgrid-1,break_POR,cscoef_POR)
	POR_DBL = dcsval(-DBL,mxgrid-1,break_POR,cscoef_POR)
	POR_eDBL = dcsval(-eDBL,mxgrid-1,break_POR,cscoef_POR)
	POR_if = dcsval(0.,mxgrid-1,break_POR,cscoef_POR)
c
c CALCULATE DIFFUSION COEFFICIENTS (Ds and Kz)
c
	do i=1,nfinal
		Ds(i) = POR(i)**2 * D
		EDDY(i) = dcsval(xfinal(i),nKz-1,break_Kz,cscoef_Kz)
     &				* 0.5*(1-derf((xfinal(i)+DBL)/sigma_DBL))
     &				+ 0.5*(1+derf((xfinal(i)+DBL)/sigma_DBL))
     &				* D_Patch*dexp(-((xfinal(i)+DBL)/TLP)**2)
	end do
c
	call dcsint(nfinal,xfinal,Ds,break_Ds,cscoef_Ds)
	Ds_bc1 = dcsval(x_bc1,mxgrid-1,break_Ds,cscoef_Ds)
	Ds_DBL = dcsval(-DBL,mxgrid-1,break_Ds,cscoef_Ds)
	Ds_eDBL = dcsval(-eDBL,mxgrid-1,break_Ds,cscoef_Ds)
	Ds_if = dcsval(0.,mxgrid-1,break_Ds,cscoef_Ds)
c
	Kz_bc1 = dcsval(x_bc1,nKz-1,break_Kz,cscoef_Kz)
	Kz_DBL = dcsval(-DBL,nKz-1,break_Kz,cscoef_Kz)
	Kz_eDBL = dcsval(-eDBL,nKz-1,break_Kz,cscoef_Kz)
	Kz_if = dcsval(0.,nKz-1,break_Kz,cscoef_Kz)
c
c CALCULATE O2 FLUXES
c
	FLUX_bc1 = -POR_bc1 * dO2_dx_bc1
     &			* (Ds_bc1+(Kz_bc1*0.5*(1-derf((x_bc1+DBL)/sigma_DBL)))
     &			+ 0.5*(1+derf((x_bc1+DBL)/sigma_DBL))
     &			* D_Patch*dexp(-((x_bc1+DBL)/TLP)**2))
	FLUX_DBL = -POR_DBL * dO2_dx_DBL
     &			* (Ds_DBL + (Kz_DBL*0.5)
     &			+ 0.5*D_Patch)
	FLUX_eDBL = -POR_eDBL * dO2_dx_eDBL
     &			* (Ds_eDBL+(Kz_eDBL*0.5*(1-derf((-eDBL+DBL)/sigma_DBL)))
     &			+ 0.5*(1+derf((-eDBL+DBL)/sigma_DBL))
     &			* D_Patch*dexp(-((-eDBL+DBL)/TLP)**2))
	FLUX_if = -POR_if * dO2_dx_if
     &			* (Ds_if+(Kz_if* 0.5*(1-derf(DBL/sigma_DBL)))
     &			+ 0.5*(1+derf(DBL/sigma_DBL))
     &			* D_Patch*dexp(-(DBL/TLP)**2))
c
c INTEGRATE OXYGEN CONSUMPTION RATE
c
	do i=1,nfinal
        	OCR(i) = R*(1.+REF*xfinal(i))*0.5*(1+derf(xfinal(i)/IGP))
     &				*0.5*(1+derf((cfinal(1,i)-O_thresh)/O_TGP))
     		OCR_ws(i) = POR(i) * OCR(i)
	end do
	call dcsint(nfinal,xfinal,OCR_ws,break_OCR_ws,cscoef_OCR_ws)
	OCR_ws_itg = dcsitg(b_bc1,x_bc2,nfinal-1,break_OCR_ws,cscoef_OCR_ws)
c
c PRINT RESULTS ON SCREEN
c
	print *, ' '
	print 100, x_bc1*10, FLUX_bc1*(0.32)	! g O2/m2 d
 	print 110, -DBL*10, FLUX_DBL*(0.32)
	print 115, -eDBL*10, FLUX_eDBL*(0.32)
 	print 120, FLUX_if*(0.32)
	print *, ' '
  	print 130, OCR_ws_itg*(0.32)
	print *, ' '
	print 140, (O2_lin - dO2_dx_lin*(eDBL-((DBL-TL)/2)))*1000
	print *, ' '
 100	format(' O2 Flux at upper boundary (',f6.2,' mm) = ',f5.3,' g/m2 d')
 110	format(' O2 Flux at top of DBL (',f6.2,' mm) = ',f5.3,' g/m2 d')
 115	format(' O2 Flux at effective DBL (',f6.2,' mm) = ',f5.3,' g/m2 d')
 120	format(' O2 Flux sed-water interface = ',f5.3,' g/m2 d')
 130	format(' Integrated OCR = ',f5.3,' g/m2 d')
 140	format(' Predicted O2 in bottom water = ',f4.0,' uM')
c
c OUTPUT SOLUTION TO FILES
c
	if(mxgrid.gt.2000) then
		nskip = int(mxgrid/2000)
	else
		nskip = 1
	end if
c
	xplot = 2.0	! cm
c
        open (unit=50,file='O2.mod',status='unknown') 
        do i=1,nfinal,nskip
                if(xfinal(i).lt.xplot) write (50,600) xfinal(i)*10, Cfinal(1,i)*1000
        end do   
        close (unit=50,status='keep')
c       
        open (unit=50,file='OCR.mod',status='unknown') 
        do i=1,nfinal,nskip 
                if(xfinal(i).lt.xplot) write (50,600) xfinal(i)*10, OCR(i)
        end do   
        close (unit=50,status='keep')
c       
        open (unit=50,file='POR.mod',status='unknown') 
        do i=1,nfinal,nskip
                if(xfinal(i).lt.xplot) write (50,600) xfinal(i)*10, POR(i)
        end do   
        close (unit=50,status='keep')
c       
        open (unit=50,file='EDDY.mod',status='unknown') 
        do i=1,nfinal,nskip
                if(xfinal(i).lt.xplot) write (50,600) xfinal(i)*10, EDDY(i)
        end do   
        close (unit=50,status='keep')
c       
        open (unit=50,file='Ds.mod',status='unknown') 
        do i=1,nfinal,nskip
		if(xfinal(i).lt.xplot) write (50,600) xfinal(i)*10, Ds(i)
        end do   
        close (unit=50,status='keep')
c       
        open (unit=50,file='DBL.mod',status='unknown') 
        write (50,610) -DBL*10, 0.
	write (50,610) -DBL*10, 1000000.
        close (unit=50,status='keep')
c       
        open (unit=50,file='eDBL.mod',status='unknown') 
        write (50,610) -eDBL*10, 0.
	write (50,610) -eDBL*10, 200.
        close (unit=50,status='keep')
c       
        open (unit=50,file='VBL.mod',status='unknown') 
        write (50,610) -VBL*10, 0.
	write (50,610) -VBL*10, 1000000.
        close (unit=50,status='keep')
c       
        open (unit=50,file='TL.mod',status='unknown') 
        write (50,610) -(DBL-TL)*10, 0.
	write (50,610) -(DBL-TL)*10, 10000.
        close (unit=50,status='keep')
c       
        open (unit=50,file='Param.mod',status='unknown') 
        write (50,*) (O2_lin - (dO2_dx_lin*DBL/2)), D, eDBL, R, por_0
        close (unit=50,status='keep')
c
 600	format(' ',f9.4,'	',f12.6)
 610	format(' ',f8.3,'	',f16.3)
        stop    
        end    
c       
c SUBROUTINE FCNEQN 
c        
        subroutine fcneqn(neq,x,C,p,dCdx)
c
	integer maxKz
	parameter (maxKz = 10)
c
	real*8 SQPI
	parameter (SQPI = 1.772453851)
c
        real*8 x, C(neq), dCdx(neq)
        real*8 p
        real*8 D, V, O_thresh, O_TGP, R, REF, IGP
	real*8 DBL, TLP, sigma_DBL, D_patch, break_Kz(maxKz), cscoef_Kz(4,maxKz)
	real*8 por_0, porbeta, porinf, PGP
	real*8 GAM1, GAM2, GAM3, GAM5, GAM6, GAM7
	real*8 OCR, POR, Ds, Kz, DIFF 
	real*8 dGAM1_dx, dGAM2_dx, dGAM5_dx, dGAM6_dx
	real*8 dPOR_dx, dDs_dx, dKz_dx, dDIFF_dx
c
	real*8 derf, dcsder, dcsval
	external derf, dcsder, dcsval
	real*8 dexp
	intrinsic dexp
c
	integer neq, nKz
c
	common/intsed/nKz
	common/realsed/D,V,O_thresh,O_TGP,R,REF,IGP,DBL,TLP,sigma_DBL,D_patch,break_Kz,cscoef_Kz
	common/realpor/por_0,porbeta,porinf,PGP
c
	GAM1 = 0.5*(1+derf(x/PGP))			! =0 if x.lt.0 else =1
	GAM2 = 0.5*(1-derf(x/PGP))			! =1 if x.lt.0 else =0
	GAM3 = 0.5*(1+derf(x/IGP))			! =0 if x.lt.0 else =1
c
	GAM5 = 0.5*(1+derf((x+DBL)/sigma_DBL))		! =0 if x.lt.DBL else =1 
	GAM6 = 0.5*(1-derf((x+DBL)/sigma_DBL))		! =1 if x.lt.DBL else =0 
	GAM7 = 0.5*(1+derf((C(1)-O_thresh)/O_TGP))	! =0 if C(1).lt.O_thresh else =1
c
	OCR = R * (1.+REF*x) * GAM3 * GAM7
c
	POR = GAM2 + GAM1*((por_0-porinf)*dexp(-porbeta*x)+porinf)
	Ds = POR**2 * D
	Kz = dcsval(x,nKz-1,break_Kz,cscoef_Kz)
	DIFF = Ds + GAM6*Kz + GAM5*D_patch*dexp(-((x+DBL)/TLP)**2)
c
	dGAM1_dx = dexp(-(x/PGP)**2)/(SQPI*PGP)
	dGAM2_dx = -dexp(-(x/PGP)**2)/(SQPI*PGP)
	dGAM5_dx = dexp(-((x+DBL)/sigma_DBL)**2)/(SQPI*sigma_DBL)
	dGAM6_dx = -dexp(-((x+DBL)/sigma_DBL)**2)/(SQPI*sigma_DBL)
c
	dPOR_dx = dGAM2_dx - GAM1*(por_0-porinf)*porbeta*dexp(-porbeta*x)
     &			+ ((por_0-porinf)*dexp(-porbeta*x)+porinf)*dGAM1_dx
	dDs_dx = 2*POR*D*dPOR_dx
	dKz_dx = dcsder(1,x,nKz-1,break_Kz,cscoef_Kz)
	dDIFF_dx = dDs_dx + Kz*dGAM6_dx + GAM6*dKz_dx 
     &			+ D_patch*dexp(-((x+DBL)/TLP)**2)*(dGAM5_dx-(2*(x+DBL)*GAM5/TLP**2))
c
        dCdx(1) = C(2)
        dCdx(2) = (OCR - ((dDIFF_dx+(DIFF*dPOR_dx/POR))*C(2)))/DIFF
c              
        return  
        end  
c  
c SUBROUTINE FCNJAC - EVALUATES JACOBIAN MATRIX 
c              
        subroutine fcnjac(neq,x,C,p,dCpdC)
c
	integer maxKz
	parameter (maxKz = 10)
c
	real*8 SQPI
	parameter (SQPI = 1.772453851)
c
        real*8 x, C(neq), dCpdC(neq,neq)
        real*8 p
        real*8 D, V, O_thresh, O_TGP, R, REF, IGP
	real*8 DBL, TLP, sigma_DBL, D_patch, break_Kz(maxKz), cscoef_Kz(4,maxKz)
	real*8 por_0, porbeta, porinf, PGP
	real*8 GAM1, GAM2, GAM3, GAM5, GAM6
	real*8 POR, Ds, Kz, DIFF 
	real*8 dGAM1_dx, dGAM2_dx, dGAM5_dx, dGAM6_dx, dGAM7_dC1
	real*8 dPOR_dx, dDs_dx, dKz_dx, dDIFF_dx
c
	real*8 derf, dcsder, dcsval
	external derf, dcsder, dcsval
	real*8 dexp
	intrinsic dexp
c
	integer neq, nKz                                                            
c
	common/intsed/nKz
	common/realsed/D,V,O_thresh,O_TGP,R,REF,IGP,DBL,TLP,sigma_DBL,D_patch,break_Kz,cscoef_Kz
	common/realpor/por_0,porbeta,porinf,PGP
c
	GAM1 = 0.5*(1+derf(x/PGP))			! =0 if x.lt.0 else =1
	GAM2 = 0.5*(1-derf(x/PGP))			! =1 if x.lt.0 else =0
	GAM3 = 0.5*(1+derf(x/IGP))			! =0 if x.lt.0 else =1
c
	GAM5 = 0.5*(1+derf((x+DBL)/sigma_DBL))		! =0 if x.lt.DBL else =1 
	GAM6 = 0.5*(1-derf((x+DBL)/sigma_DBL))		! =1 if x.lt.DBL else =0 
c
	POR = GAM2 + GAM1*((por_0-porinf)*dexp(-porbeta*x)+porinf)
	Ds = POR**2 * D
	Kz = dcsval(x,nKz-1,break_Kz,cscoef_Kz)
	DIFF = Ds + GAM6*Kz + GAM5*D_patch*dexp(-((x+DBL)/TLP)**2)
c
	dGAM1_dx = dexp(-(x/PGP)**2)/(SQPI*PGP)
	dGAM2_dx = -dexp(-(x/PGP)**2)/(SQPI*PGP)
	dGAM5_dx = dexp(-((x+DBL)/sigma_DBL)**2)/(SQPI*sigma_DBL)
	dGAM6_dx = -dexp(-((x+DBL)/sigma_DBL)**2)/(SQPI*sigma_DBL)
	dGAM7_dC1 = dexp(-((C(1)-O_thresh)/O_TGP)**2)/(SQPI*O_TGP)
c
	dPOR_dx = dGAM2_dx - GAM1*(por_0-porinf)*porbeta*dexp(-porbeta*x)
     &			+ ((por_0-porinf)*dexp(-porbeta*x)+porinf)*dGAM1_dx
	dDs_dx = 2*POR*dPOR_dx
	dKz_dx = dcsder(1,x,nKz-1,break_Kz,cscoef_Kz)
	dDIFF_dx = dDs_dx + Kz*dGAM6_dx + GAM6*dKz_dx 
     &			+ D_patch*dexp(-((x+DBL)/TLP)**2)*(dGAM5_dx-(2*(x+DBL)*GAM5/TLP**2))
c
        dCpdC(1,1) = 0.0
        dCpdC(1,2) = 1.0  
        dCpdC(2,1) = R*(1.+REF*x)*GAM3*dGAM7_dC1/DIFF
        dCpdC(2,2) = -(dDIFF_dx+(DIFF*dPOR_dx/POR))/DIFF
c
        return 
        end     
c                                                                               
c SUBROUTINE FCNBC - EVALUATES BOUNDARY CONDITIONS 
c                                                                               
        subroutine fcnbc(neq,ya,yb,p,f)
c
        real*8 ya(neq), yb(neq), f(neq)
        real*8 p 
        real*8 C_bc1, C_bc2 
c
        integer neq 
c
	common/realbound/C_bc1,C_bc2
c      
        f(1) = ya(1) - C_bc1
        f(2) = yb(1) - C_bc2
c      
        return  
        end