c
c
      subroutine invproc(nlay,thick,xcomp) 

*$noereference
      include 'p1werm.inc'
*$reference

c     + + + PURPOSE + + +
c     
c	  Invert the component passed to xcomp 
c	  
c     + + + KEYWORDS + + +
c     inversion, tillage 
c
c     + + + ARGUMENT DECLARATIONS + + +

      real xcomp(mnsz), thick(mnsz)
      integer nlay
c
c     + + + ARGUMENT DEFINITIONS + + +

c     nlay		- number of soil layers used
c     thick		- thickness of each layer in a subregion
c	  xcomp		- component that needs inverting
c
c     + + + LOCAL VARIABLES + + +

      integer i
	  real dum(mnsz),dum1(mnsz),dum2(mnsz) 
	  real inverse(mnsz), depth(mnsz), oldepth(0:mnsz)
	  real x, new,  p(mnsz)  
c
c     + + + LOCAL VARIABLE DEFINITIONS + + +
c
c     dum		- dummy variable used in making a property array
c     dum1		- dummy variable used in making a property array
c     dum2		- dummy variable used in making a property array
c     depth		- depth matrix containing inverted layer depths from the surface
c     i			- loop variable for soil layers
c     inverse	- inverse layer thickness matrix
c	  new		- new inverted property value based on interpolation or
c	  			  extrapolation
c     oldepth	- depth matrix containing original layer depths from the surface
c     p			- property matrix after inversion
c     x			- variable containing the original depths from the surface
c				  used in the call to polint
	   
c     Do the inversion process. 
c     Initialize the dummy variables to zero.

	  dum(1) = 0.0
	  dum1(1) = 0.0
	  dum2(1) = 0.0
	  oldepth(0) = 0.0
	  
	  do 100 i=1,nlay
c      invert the layers (layer thickness) 
	    inverse(i) = thick((nlay+1)-i)
c      invert the property passed to xcomp  
		p(i) = xcomp((nlay+1)-i)
c     form a property array for the depth (thichness) and component based
c     on the inverted matrix 
	  dum(i+1) = thick(i)
	  dum1(i+1) = inverse(i) 
	  depth(i)=inverse(i)/2.0+dum1(i)/2.0+dum2(i) 
	  dum2(i+1)= depth(i)
	  oldepth(i)=thick(i)/2.0+dum(i)/2.0+oldepth(i-1)
100   continue
c
c     make a call to subroutine polint which takes the current
c     property matrix and depth matrix and either interpolates
c     or extrapolates the property matrix to correspond to the
c     original layer thickness before the inversion process 
c     was performed.  Make call for each layer (nlay). 
c
      do 200 i=1,nlay
        x=oldepth(i) 
	    call polint(depth, p, nlay, x, new)
c
c     set the component for a layer equal to the interpolated
c     or extrapolated value calculated in polint.for 
c
	    xcomp(i) = new
200   continue   
	  end