anchor.diff:  C     reset count of nRed, nFull
                 ireset = 1
              endif
  -	elseif(jsolv.eq.1) then !model already selected for these conditions
  -	   imodNum = iModel(igrid)
  -	   kr = 0
  -	   do k=1,kk
  -	      yfull(k) = 0.0d0
  -	      if(nus(k,imodNum) .eq. 1) then
  -		 kr = kr + 1
  -		 yfull(k) = y(kr)
  -		 yinert(k,igrid) = 0.0d0 !frees up some memory
  -	      else
  -		 yfull(k) = yinert(k,igrid)
  -	      endif
  -	   enddo
  -	   call ackytcp(p,t,yfull,ickwrk,rckwrk,cfull)
  -
  -	endif ! "if jsolv=0/1" conditional
  +      elseif(jsolv.eq.1) then !model already selected for these conditions
  +         imodNum = iModel(igrid)
  +         kr = 0
  +         do k=1,kk
  +            yfull(k) = 0.0d0
  +            if(nus(k,imodNum) .eq. 1) then
  +                kr = kr + 1
  +                yfull(k) = y(kr)
  +                yinert(k,igrid) = 0.0d0 !frees up some memory
  +            else
  +                yfull(k) = yinert(k,igrid)
  +            endif
  +         enddo
  +         call ackytcp(p,t,yfull,ickwrk,rckwrk,cfull)
  +
  +      endif ! "if jsolv=0/1" conditional
   C     <initialize parameters>
   C     On first call, initialize, read inputs, and save data
  -	if(iCall(1).ne.2) call adapchinit
  +      if(iCall(1).ne.2) call adapchinit
        & (igridTot,icount,igrTotm1,icountm1,ii,kk,kfuel,resid_tol,
        &  jprint,nprint,nsteps,iicount,iModel,nFull,nRed,rmaxfull)
   C********NOTE: OLDER VERSIONS ASSUME ARRAYS (eg. ICAL()) ARE ALWAYS INITIALIZED
   C              TO 0 BY YOUR COMPILER. I.E. ICAL(I) = 0 FORALL I
   C     </initialize>
  -C       
  +C
   C     is a new model needed?
  -	if(jsolv .eq. -1) then
  +      if(jsolv .eq. -1) then
   C     this call is to load reduced models only; then return
  -	   jfindmod = 1
  -	   iwritesize = 0
  -	else	      
  -	   call needmodel(igrid,igridTot,igrTotm1,nsteps,icount,icountm1,
  +         jfindmod = 1
  +         iwritesize = 0
  +      else
  +         call needmodel(igrid,igridTot,igrTotm1,nsteps,icount,icountm1,
        & jfindmod,iwritesize,jprint,nprint,iicount)
  -	endif
  +      endif
   C     record model usage statistics to date
  -	if(iwritesize.eq.1) then
  -!Luwi--10/15/10--don't write stats!	   call sizewrite(nFull,nRed,iSizeTot,icountTot,ii,kk)
  +      if(iwritesize.eq.1) then
  +!Luwi--10/15/10--don't write stats! call sizewrite(nFull,nRed,iSizeTot,icountTot,ii,kk)
   C     <modelLibrary>
   C     skip the following if a new model library is being constructed
  -	   if(inewlib.gt.0) goto 007
  -  007	   continue
  +         if(inewlib.gt.0) goto 007
  +  007    continue
   C     </modelLibrary>
  -	endif
  -C	reset iSizeTot and icountTot
  +      endif
  +C     reset iSizeTot and icountTot
         if(ireset .eq. 1) then
           if(jfindmod.eq.1) then
  -	      iSizeTot = 0
  -	      icountTot = 0
  -	      nFull = 0
  -	      nRed = 0
  -	   endif
  -	   icountm1 = icount
  -	endif
  +            iSizeTot = 0
  +            icountTot = 0
  +            nFull = 0
  +            nRed = 0
  +         endif
  +         icountm1 = icount
  +      endif
   C     <Update models: only if necessary>
   C      Find smallest valid reduced model
  -	if(jfindmod.eq.1) then
  -	call findModel(jsolv,1,
  +      if(jfindmod.eq.1) then
  +      call findModel(jsolv,1,
        & iCall(1), kk,ii, kfuel, T, Cfull, imodNum,nFull,
        & nRed)
  -	if (jsolv .eq. -1) return !Initialization call only
  -	   iModel(igrid) = imodNum
  -	   iSizeTot = iSizeTot + modelSize(imodNum,2)!computes avg number of species
  -	   icountTot = icountTot + 1
  -	endif
  -C     Once models obtained for all grid points, don't search for new models 
  +      if (jsolv .eq. -1) return !Initialization call only
  +         iModel(igrid) = imodNum
  +         iSizeTot = iSizeTot + modelSize(imodNum,2)!computes avg number of species
  +         icountTot = icountTot + 1
  +      endif
  +C     Once models obtained for all grid points, don't search for new models
   C     again until retired models or new mesh.
   C     </Update models>
   
  -C	<Calculate rates using selected model>
  -	if(jsolv .eq. 0) then
  +C     <Calculate rates using selected model>
  +      if(jsolv .eq. 0) then
   C     this call was for determining reduced model only; AND
   C     it is possible reduced model already determined on prior call
   C     e.g. ODE solver reducing stepsize and repeating integration
   C     after prior failed attempt
  -	   imodNum = iModel(igrid)
  -	   do k=1,kk
  -	      ispecies(k) = nus(k,imodNum)
  -	   enddo
  -	elseif(jsolv.eq.1) then
  - 	   imodNum = iModel(igrid)
  -	   imodSize = modelSize(imodNum,1)
  -	   kmodSize = modelSize(imodNum,2)
  -     
  +         imodNum = iModel(igrid)
  +         do k=1,kk
  +            ispecies(k) = nus(k,imodNum)
  +         enddo
  +      elseif(jsolv.eq.1) then
  +         imodNum = iModel(igrid)
  +         imodSize = modelSize(imodNum,1)
  +         kmodSize = modelSize(imodNum,2)
  +
   C       chemkin call uses species from reduced-order ODE mapped into full concentration array
   C
  -	   call ackwc(imodSize,p,t,cfull,nur(1,imodNum),nus(1,imodNum),
  +         call ackwc(imodSize,p,t,cfull,nur(1,imodNum),nus(1,imodNum),
        & ickwrk,rckwrk,ispecies,wdot)
   C     return the reduced model species binary array to CFD code
  -	   do k=1,kk
  -	      ispecies(k) = nus(k,imodNum)
  -	   enddo
  -	elseif(jsolv .gt. 1) then
  -	   print*,'Error: Invalid jsolv value:',jsolv
  -	   print*,'jsolv must be -1, 0 or 1'
  -	   stop
  -	endif
  +         do k=1,kk
  +            ispecies(k) = nus(k,imodNum)
  +         enddo
  +      elseif(jsolv .gt. 1) then
  +         print*,'Error: Invalid jsolv value:',jsolv
  +         print*,'jsolv must be -1, 0 or 1'
  +         stop
  +      endif
   C     </Calc. rates>
   
  -	return
  -	end
  +      return
  +      end
   
   C*************************************************************************
  -	subroutine needmodel(igrid,igridTot,igrTotm1,nsteps,icount,
  +      subroutine needmodel(igrid,igridTot,igrTotm1,nsteps,icount,
        & icountm1,jfindmod,iwritesize,jprint,nprint,iicount)
  -	implicit double precision (a-h, o-z), integer (i-n)
  -	dimension iicount(*)
  +      implicit double precision (a-h, o-z), integer (i-n)
  +      dimension iicount(*)
   C     Keep track of age of current model distribution
   C     icount changes on a new iteration, remains constant for all calls
   C     on the same iteration (e.g calls for jacobian calculations)
   C     Force update if number of gridpoints has changed
  -	if(igridTot.ne.igrTotm1) then
  -	   do j = 1, igridTot
  -	      iicount(j) = -nsteps
  -	   enddo
  +      if(igridTot.ne.igrTotm1) then
  +         do j = 1, igridTot
  +            iicount(j) = -nsteps
  +         enddo
   C     <Update>
  -	   igrTotm1 = igridTot
  +         igrTotm1 = igridTot
   C     </Update>
  -	endif
  +      endif
   C     default: don't need to update models
  -	jfindmod = 0
  +      jfindmod = 0
   C     default: don't write model use stats to file
  -	iwritesize = 0
  +      iwritesize = 0
   C     iicount(igrid) = iteration number of most recent model update for igrid
  -	iicountI = iicount(igrid)
  -	if(iicountI.ne.icount) then
  +      iicountI = iicount(igrid)
  +      if(iicountI.ne.icount) then
   C     new iteration/integration step at this grid point;
  -	   if(icount.eq.0) then
  +         if(icount.eq.0) then
   C     user wishes to force model update
  -	      istepI = 0
  -	   else
  +            istepI = 0
  +         else
   C     determine whether update is needed
  -	      istepI = mod(iicountI,nsteps)
  -	   endif
  -	   if(istepI.eq.0) then
  +            istepI = mod(iicountI,nsteps)
  +         endif
  +         if(istepI.eq.0) then
   C     Models at retirement age, or first call, or new mesh
  -	      jfindmod = 1
  -	   endif
  +            jfindmod = 1
  +         endif
   C     update iicount(igrid)
  -	   iicount(igrid) = icount
  -	endif
  +         iicount(igrid) = icount
  +      endif
   C     <modelUsage statistics>
  -	if(icount.ne.icountm1) then
  +      if(icount.ne.icountm1) then
   C     new step
   C     increment count of number of steps between 'modelUsage' writes
  -	   jprint = jprint + 1
  -	   if(jprint.eq.nprint) then
  -	      jprint = 0
  -	      iwritesize = 1
  -	   endif
  -	endif
  +         jprint = jprint + 1
  +         if(jprint.eq.nprint) then
  +            jprint = 0
  +            iwritesize = 1
  +         endif
  +      endif
   C     </modelUsage>
   
  -	return
  -	end
  +      return
  +      end
   
   C************************************************************************
  -	subroutine adapchinit(igridTot,icount,igrTotm1,icountm1,ii,kk,
  +      subroutine adapchinit(igridTot,icount,igrTotm1,icountm1,ii,kk,
        & kfuel,resid_tol,jprint,nprint,nsteps,iicount,iModel,nFull,nRed,
        & rmaxfull)
  -	implicit double precision (a-h,o-z), integer (i-n)
  -	dimension iModel(*), iicount(*), iCall(1)
  -	character line*5
  +      implicit double precision (a-h,o-z), integer (i-n)
  +      dimension iModel(*), iicount(*), iCall(1)
  +      character line*5
   
   C       initialize and set defaults on number of steps between modelUsage write
  -	jprint = 0
  -	nprint = 1
  +      jprint = 0
  +      nprint = 1
   C     for calculating average model size; printed to avgModelSize
  -	icountTot = 0
  -	iSizeTot = 0
  +      icountTot = 0
  +      iSizeTot = 0
   C     modelUsage:
  -	nFull = 0
  -	nRed = 0
  -	rmaxfull = 1.
  -	igrTotm1 = igridTot
  -	icountm1 = icount
  -	do j=1, igridTot
  -	   iModel(j) = 0
  -	   iicount(j) = 0
  -	end do
  -	print*,'reading adapchem.in...'
  -	open(7, file='adapchem.in', status='old')
  -	do while(7.gt.0)
  -	   read(7,*) line
  -	   if(line(:3).eq.'END') goto 005
  -	   if(line.eq.'TOLSS') read(7,*) tolss
  -	   if(line.eq.'TOLMR') read(7,*) tolmr
  -	   if(line.eq.'NSTEP') read(7,*) nsteps
  -	   if(line.eq.'NPRIN') read(7,*) nprint
  +      nFull = 0
  +      nRed = 0
  +      rmaxfull = 1.
  +      igrTotm1 = igridTot
  +      icountm1 = icount
  +      do j=1, igridTot
  +         iModel(j) = 0
  +         iicount(j) = 0
  +      end do
  +      print*,'reading adapchem.in...'
  +      open(7, file='adapchem.in', status='old')
  +      do while(7.gt.0)
  +         read(7,*) line
  +         if(line(:3).eq.'END') goto 005
  +         if(line.eq.'TOLSS') read(7,*) tolss
  +         if(line.eq.'TOLMR') read(7,*) tolmr
  +         if(line.eq.'NSTEP') read(7,*) nsteps
  +         if(line.eq.'NPRIN') read(7,*) nprint
   C     <modelLib>
  -	   if(line.eq.'MAXFU') read(7,*) rmaxfull
  +         if(line.eq.'MAXFU') read(7,*) rmaxfull
   C     </modelLib>
  -	enddo
  - 005	continue
  -	close(7)
  -	resid_tol = tolss - tolmr
  -	print*,'resid_tol= ',resid_tol
  -	print*,'nsteps= ',nsteps
  -	print*,'rmaxfull= ',rmaxfull
  -	print*,'reading models.chem...'
  -	
  -	open(1, file='models.chem', status='old')
  -	read (1,*) ii, kk!Luwi-9/2/10-kfuel not needed/used!, kfuel
  +      enddo
  + 005  continue
  +      close(7)
  +      resid_tol = tolss - tolmr
  +      print*,'resid_tol= ',resid_tol
  +      print*,'nsteps= ',nsteps
  +      print*,'rmaxfull= ',rmaxfull
  +      print*,'reading models.chem...'
  +
  +      open(1, file='models.chem', status='old')
  +      read (1,*) ii, kk!Luwi-9/2/10-kfuel not needed/used!, kfuel
   C     closed in subroutine findModel
   !       <Luwi--9/3/2010-- avoiding multiple file copies in parallel code>
   CC     Write headers for model usage documentation
  -C	open(unit=2,file='avgModelSize',status='unknown')
  -C	write(2,'(A14)') '======================'
  -C	write(2,'(A4,X,A9)') 'nRxns','% Of Full'
  -C	write(2,'(A14)') '======================'
  +C     open(unit=2,file='avgModelSize',status='unknown')
  +C     write(2,'(A14)') '======================'
  +C     write(2,'(A4,X,A9)') 'nRxns','% Of Full'
  +C     write(2,'(A14)') '======================'
   CC
  -C	open(unit=3,file='modelUsage',status='unknown')
  -C	write(3,'(A22)')        '======================'
  -C	write(3,'(A10,2X,A10)') '     nFull','  nReduced'
  -C	write(3,'(A22)')        '======================'
  +C     open(unit=3,file='modelUsage',status='unknown')
  +C     write(3,'(A22)')        '======================'
  +C     write(3,'(A10,2X,A10)') '     nFull','  nReduced'
  +C     write(3,'(A22)')        '======================'
   !       <Luwi>
  -	return
  -	end
  +      return
  +      end
   
   C**********************************************************************
  -	subroutine sizewrite(nFull,nRed,iSizeTot,icountTot,ii,kk)
  -	implicit none
  -	integer iSizeTot, icountTot, ii, nFull, nRed,kk
  -	write(2,1) nint(dble(iSizeTot)/dble(icountTot)),
  +      subroutine sizewrite(nFull,nRed,iSizeTot,icountTot,ii,kk)
  +      implicit none
  +      integer iSizeTot, icountTot, ii, nFull, nRed,kk
  +      write(2,1) nint(dble(iSizeTot)/dble(icountTot)),
        &dble(iSizeTot)/dble(icountTot*kk)*100.
  -	write(3,2) nFull, nRed
  -	
  -    1	format(I4,2X,F5.1)
  -    2	format(I10,2X,I10)
  -
  -	return
  -	end
  +      write(3,2) nFull, nRed
  +
  +    1 format(I4,2X,F5.1)
  +    2 format(I10,2X,I10)
  +
  +      return
  +      end
   
   C**********************************************************************
  -	subroutine readModels(lin,kk,ii,rangeLo,rangeHi,nModels)
  +      subroutine readModels(lin,kk,ii,rangeLo,rangeHi,nModels)
   C     O.Oluwole 6/9/2004
   C     This subroutine reads in the reduced models in the model library
   C     for use during the adaptive chemistry calculation
  -	implicit double precision (a-h,o-z), integer (i-n)
  -	include 'ackparams.i'
  -	dimension rangeLo(kmax+1,nm+1), rangeHi(kmax+1,nm+1)
  -	dimension modelSize(nm+1,2)
  +      implicit double precision (a-h,o-z), integer (i-n)
  +      include 'ackparams.i'
  +      dimension rangeLo(kmax+1,nm+1), rangeHi(kmax+1,nm+1)
  +      dimension modelSize(nm+1,2)
         dimension NUR(imax,nm+1), NUS(kmax,nm+1)
  -	common /saveNUS/ NUS
  +      common /saveNUS/ NUS
         common /saveNUR/ NUR, modelSize
  -	character line*12
  +      character line*12
   C     Columns of NUR(,) contain all models including full model
  -C     Columns of rangeLo and rangeHi contain valid ranges of reduced 
  -C     models beginning in column 2       
  +C     Columns of rangeLo and rangeHi contain valid ranges of reduced
  +C     models beginning in column 2
   C     all reduced models are stored in order of increasing size
   C     <O.Oluwole--6/26/2008>
   C     Aerodyne Research, Inc.
  
positive.code: C     starting a new iteration (or integration) step.
  C     'resid_tol' is the steady-state convergence tolerance:
  C     steady state when max(|dT/dt|,|dC/dt|)<resid_tol
  C     This should give you an idea of where to insert the common block 
  C     in your CFD solver 
  C     resid_tol should ONLY be changed by adapchem routine
  C *********************************************************************
  C
  
negative.code: C
        DO 25 N = 1, NLAN
           I = ILAN(N)
         if(nur(i).eq.0) goto 25
           TFAC = PLT(1,N)/T**(1.0/3.0) + PLT(2,N)/T**(2.0/3.0)
           RKFT(I) = RKFT(I) * EXP(TFAC)
     25 CONTINUE
  C
        CALL ACKSMH (T, ICKWRK, RCKWRK, SMH)
        DO 50 I = 1, II
         if(nur(i).eq.0) goto 50
            SUMSMH = 0.0
            DO 40 N = 1, MAXSP
               IF (NUNK(N,I).NE.0)
  