text:
      subroutine decomp(nowcrp)
c
c     + + + PURPOSE + + +
c
c           XXX-NOTE:  The computation of POP looks like a real problem that
c           needs to be corrected!  The approach seems unsound,
c           and the location of this code in this subroutine seems
c           inappropriate.  (See note at end of DECOMP.)
c           -- CRM -- 12/01/92.
c
c     Estimates changes in residue and root masses, computes remaining
c     flat, standing, and buried residue, and root masses.  Computes
c     residue cover (COVCAL), and makes adjustments for tillage oper-
c     ations.  Also estimates mass and cover for flat and standing
c     residues.  Predicts buried residue mass and non-living root mass.
c     Considers up to three successive types of fallen and buried
c     residues: only one type of standing residue.
c
c
c  ADDENDUM:  Default parameter for wheat    ==> ORATEA = 0.0110
c                                   corn     ==>        = 0.0100
c                                   soybeans ==>        = 0.0200
c             For now, ORATER should continue to be set equal to
c             ORATEA, unless data determines otherwise.
c             -- D.E. Stott 9/2/92
c
c  ********************************************************************
c  *  NOTE: In the User's Documentation equation 8.5.3c is computed   *
c  *        on a weight basis.  OPTWAT (used in that equation) is     *
c  *        also given on a weight basis.  However, the computations  *
c  *        in equation 8.5.3c are done on a volume basis.  This      *
c  *        needs to be corrected. -- CRM -- 11/13/92.                *
c  ********************************************************************
c
c           XXX-NOTE:  Why are CRPPRM and PEREN separate common blocks?  They
c           seem to contain the same sorts of variables -- some of which are
c           appropriate to BOTH perenials AND annuals -- variables like
c           JDPLT & JDHARV are in CRPPRM, whereas JDCUT, JDBURN, etc. are
c           in PEREN. -- CRM -- 2/19/92
c
c     Called from WATBAL (If the land use is cropland)
c     Author(s): Alberts, Stott, Ferris, Meyer
c     Reference in User Guide: Chapter 8
c
c     Changes:
c          1) Local variable NRES deleted, since it was always set to
c             '4'.
c          2) Parameter TILTYP not used.  Include file PTILTY.INC
c             dereferenced.
c          3) The variables in include files CCRPVR5.INC, CCOVER.INC,
c             CCRPOUT.INC, and CSOLVAR.INC were not used.  They were
c             dereferenced.
c          4) The variables in the SAVE statement below never have an
c             opportunity to change value after the DATA initialization.
c             The statement was deleted.
c                save optwat,tfc,opttmp,varind,ferind,pszind
c          5) Deleted second calculation of ENV  in "do 999" loop,
c             since it was identical to the first.
c          6) Deleted use of NOWRES as a local variable with a con-
c             stant value of '1'.
c          7) Created following local variables to avoid repeated recal-
c             culation of:
c               TMPVR1:
c                     oratea(iresd(nowres,iplane))*varind*ferind*pszind
c               TMPVR2:
c                     orater(iresd(nowres,iplane))*varind*ferind*pszind
c               TMPVR3:
c                     1.0 - fbrnog(4,iplane)
c                     (It was being computed 3 times for each trip
c                      through the "do 15" loop.)
c               TMPVR4:
c                     rmagy-rmagt(iplane)
c               TMPVR5:
c                     1.0 - frmove(4,iplane)
c          8) TMPVR* calculations taken outside their respective
c             loops.
c          9) VARIND was deleted at request of D.E. Stott.
c             CRM -- 3/27/92
c         10) Calculation of RMAGT was determined to be incorrect.
c             Calculations of SENVIN, FENVIN, & BENVIN were changed
c             from cumulative to daily to correct for this, as per
c             instructions from D.E. Stott.  -- CRM -- 3/27/92
c         11) Changed the "residue management" IF's to ELSE-IF's.
c         12) Since RAIN is expressed in *meters* (not millimeters),
c             Equation 8.3.5a was incorrect by a factor of 1000.
c             This was corrected by D.E. Stott -- 9/2/92.
c             Changed from:
c                 if(rain.lt.4.0) then
c                   swatfc = rain/4.0
c             to:
c                 if(rain.lt.0.004) then
c                   swatfc = rain/0.004
c               CRM -- 10/30/92
c         13) Since SOILW is the volume of soil water (meters) in the
c             layer, and Equation 8.3.5c uses percent saturation by
c             volume; ie, soil water (meters) / layer thickness
c             (meters), Equation 8.3.5c was incorrect.  This correct-
c             ed by D.E. Stott -- 9/2/92.
c             Changed by D. Stott from:
c                 if(soilw(1,iplane).lt.optwat) then
c                   fwatfc = soilw(1,iplane)/optwat
c             to:
c                 if((soilw(1,iplane)/dg(1,iplane)).lt.optwat) then
c                   fwatfc = (soilw(1,iplane)/dg(1,iplane))/optwat
c               CRM -- 10/30/92
c -- XXX -- Or WAS it fully corrected?  The note with the equation
c           indicates it needs percent by WEIGHT.  -- CRM -- 11/06/92.
c         14) Equation to calculate root decomposition was *supposed* to
c             use ORATE*R*.  It was being calculated using ORATE*A*.
c             this was corrected by D.E. Stott -- 9/2/92.
c         15) OPTWAT, TFC, & OPTTMP were changed from 0.35, 6.1, & 30.0
c             to 0.30, 6.4, & 33.0 respectively. -- CRM -- 11/09/92.
c         16) The following line was added to the "Flat residue removed"
c             section:
c                 rilrm(nowres,iplane)=rilrm(nowres,iplane)*tmpvr5
c             and the following lines were removed:
c                 if((tilseq(nowcrp,iplane).gt.0) .and.
c            1       (iridge(tilseq(nowcrp,iplane)).ne.1))
c            2         rilrm(nowres,iplane)=rigrm(nowres,iplane)
c         17) Since residue removal only deals the FLAT residue (not
c             standing), the following line was deleted from the residue
c             removal management section:
c                 rmagt(iplane)=rmagt(iplane)*tmpvr5
c         18) While vainly attempting to convert cover calculations
c             to weight calculations, I noticed we had residue falling
c             UP the furrow!  Added the following statement to correct
c             the problem. -- CRM -- 11/30/92
c                 if(deltrm .lt. 0.0) deltrm = 0.0
c         19) For ridge-till systems, residue movement is supposed to
c             be simulated for 60 days after harvest (not just 59).
c             the following line:
c                 elseif(dah.gt.0 .and. dah.lt.60)then
c             changed to:
c                 elseif(dah.gt.0 .and. dah.le.60)then
c         20) Updated the calculation of residue flat on the ground
c             today (RMOGT) to reflect actual interrill/rill fractions
c             rather than assuming a 50:50 split as per Livingston.
c         21) References to JDBURN(4,iplane), JDCUT(4,iplane), and
c             JDMOVE(4,iplane) were changed to JDBURN(NOWCRP,iplane),
c             JDCUT(NOWCRP,iplane), and JDMOVE(NOWCRP,iplane), since
c             these are the values actually acquired in RESUP.
c         22) The line:
c                 if(tilseq(nowcrp,iplane).gt.0) then
c             was added back to the M3 & N1 loops to check for
c             positive subscripts on IRIDGE. -- CRM -- 12/10/92.
c         23) The variable WGHT1 was added at the beginning of the
c             residue management/shifting section to handle S. Living-
c             ston's new calculation of RMOGT. -- CRM -- 12/10/92.
c         24) The lines:
c                rmogt(nowres,iplane)=wght1*rilrm(nowres,iplane) +
c               1                     (1.0-wght1)*rigrm(nowres,iplane)
c             were changed to:
c                rmogt(1,iplane)=wght1*rilrm(1,iplane) +
c               1                     (1.0-wght1)*rigrm(1,iplane)
c             This was to correct an oversight.  NOWRES was not defined.
c             CRM -- 12/10/92.
c         25) The lines:
c                rmogt(1,iplane)=wght1*rilrm(1,iplane) +
c               1                     (1.0-wght1)*rigrm(1,iplane)
c             were changed to:                      ^
c                rmogt(1,iplane)=wght1*rilrm(1,iplane) +
c               1                     (1.0-wght1)*rilrm(1,iplane)
c             CRM -- 12/10/92.
c         26) The lines:
c                rmogt(1,iplane)=wght1*rilrm(1,iplane) +
c               1                     (1.0-wght1)*rigrm(1,iplane)
c             were incorrect (D. Flanagan) They were changed to:
c                rmogt(1,iplane)=wght1*riGrm(1,iplane) +
c               1                     (1.0-wght1)*riLrm(1,iplane)
c             CRM -- 12/31/92.
c         27) Moved call to DECOMP from within SR WATBAL to SR CONTIN,
c             because we need to perform residue management and tillage
c             and update residue cover BEFORE we compute the
c             infiltration parameters for today.  Some questions remain
c             as to how this move will impact the actual decomposition
c             itself (will be using yesterday's soil moisture in soil
c             layer 1, with today's rainfall and temperature info).
c             Best solution would be to split DECOMP into two parts -
c             the first which has only decomposition and is called by
c             SR WATBAL at the end of the day, and the second which
c             would have the residue management code and would be
c             called by SR CONTIN before the current day's infiltration
c             and erodibility parameter adjustments have been made.
c             DCF -- 2/3/93
c         28) Changed 'optwat' from a constant to a percentage of the
c             soil porosity in the top soil layer, 'avpor' This allows
c             the optimum water content to be set equal to 60% of the
c             total porosity filled with water. This is an extablished
c             optimum for soil microbial activity.  This also allows
c             'optwat' to be dynamic, allowing for changes in bulk
c             density. 'avpor' is called from ccons.inc.
c                         optwat = avpor(iplane)*0.60
c             DES -- 4/3/93
c         29) Streamlined the calculation of the tmpfac.  See comments
c             and equations for tmpfac below.  DES -- 4/3/93
c         30) Deleted the M3 IF-ELSE and altered the code so that all
c             systems (ridged or non-ridged) use the same computations
c             for tillage-induced effects on surface residue.  This will
c             allow easy insertion of separate MFO values for interrill
c             and rill areas, as requested by SCS.   dcf  8/25/93
c         31) variable initialization in data statements moved to
c             code body  jca2  8/31/93
c         32) Altered the equation for calculating the watfac for
c             wet conditions to allow for the impact of oxygen
c             depletion.  For water contents (%WFP) > 60%
c
c              fwatfc = optwat/(soilw(1,iplane)/dg(1,iplane))
c
c             DES -- 1/30/94
c         33) Added code to use different MFO (tillage intensity)
c             values for interrill and rill areas.  dcf  2/3/94
c         34) Changed from using variable RAIN to variable PRCP
c             and now also use TAVE in calculation of SWATFC,
c             because new winter routines (contin.for) no longer
c             set RAIN to zero when TAVE is less than zero.
c             Perhaps we should use TMAX instead of TAVE??? dcf 6/3/94
c         37) Changed the calculation of the water function to
c             use soil water information for top 2 soil layers since
c             comparison was being made with AVPOR which is computed
c             based on these top 2 layers.  dcf  5/4/95
c         36) Changed limits for temperature function.
c             FROM  if (tave.gt.-1.5.and.tave.lt.44.9) then
c               TO  if (tave.gt.-6.1.and.tave.lt.49.2) then
c             based on newest DECOMP documentation and check
c             with Stott.   dcf  5/5/95
c
c
c     Version: This module recoded from WEPP version 91.50.  Later
c              updated to version 92.25.
c
c     Date recoded: 12/26/91 - 02/19/92.
c              and: 05/12/92 - 05/12/92.
c              and: 10/30/92 - 12/01/92.
c     Recoded by: Charles R. Meyer.
c     Changes dated in April 1993 made to WEPP version 93.0 (DES)
c     Changes dated in January 1994 made to WEPP version 93.13 (DES)
c
c     + + + KEYWORDS + + +
c
c     + + + PARAMETERS + + +
      include 'pmxcrp.inc'
      include 'pmxcut.inc'
      include 'pmxelm.inc'
      include 'pmxhil.inc'
      include 'pmxnsl.inc'
      include 'pmxpln.inc'
      include 'pmxres.inc'
      include 'pmxtil.inc'
      include 'pmxtls.inc'
      include 'pntype.inc'
c
c     + + + ARGUMENT DECLARATIONS + + +
      integer nowcrp
c
c     + + + ARGUMENT DEFINITIONS + + +
c     nowcrp - index of the current crop
c
c     + + + COMMON BLOCKS + + +
      include 'ccons.inc'
c       read: avpor(mxplan)
c
      include 'ccrpvr1.inc'
c       read: mfo(mxtill,mxtlsq,mxcrop),rmfo(mxtill,mxtlsq,mxcrop)
c     modify: rmogt(mxres,mxplan),  rmagt(mxplan), smrm(mxres,mxplan),
c             rtm(mxres,mxplan).
c
      include 'ccrpvr2.inc'
c       read: cn(ntype), aca(ntype), cf(ntype), ar(ntype)
c
      include 'cclim.inc'
c       read: tave
c
      include 'ccontcv.inc'
c       read: tilseq(mxcrop,mxplan)
c
      include 'ccrpprm.inc'
c       read: jdplt(mxcrop,mxplan), jdharv(mxcrop,mxplan),
c             iresd(mxres,mxplan), iroot(mxres,mxplan)
c
      include 'cdecvar.inc'
c       read: oratea(ntype)
c     modify: senvin(mxplan), fenvin(mxres,mxplan), benvin(mxres,mxplan)
      include 'cerrid.inc'
c       read: crpnam(iplane)
c
      include 'chydrol.inc'
c
      include 'cnew.inc'
c       read: manver
c       read: prcp
      include 'csenes.inc'
c       modify csenes(iplane) if residue addition
c
      include 'cperen.inc'
c     modify: pop(mxplan)
c       read: popmat(mxplan), jdburn(mxcrop,mxplan),
c             jdcut(mxcrop,mxplan), jdmove(mxcrop,mxplan),
c             fbrnog(mxcrop,mxplan), fbrnag(mxcrop,mxplan),
c             frmove(mxcrop,mxplan), srmhav(mxplan), fact(ntype)
c
      include 'cridge.inc'
c       read: iridge(mxtlsq)
c     modify: rilrm(mxplan),rigrm(mxres,mxplan)
c
      include 'cstruc.inc'
c       read: iplane
      include 'ctillge.inc'
c       read: resman
c
      include 'cupdate.inc'
c       read: sdate, mdate(mxtill,mxtlsq)
c     modify: indxy(mxplan)
c
      include 'cwater.inc'
c       read: soilw(mxnsl,mxplan)
c
      include 'cends4.inc'
c       read: rspace(mxplan), width(mxplan)
c
c     + + + LOCAL VARIABLES + + +
      real flrcov, rmogy(mxres), optwat, atemp,
     1    t1, t2, tmpfac, swatfc, fwatfc, ferind, pszind,
     1    tmpvr1, tmpvr2, tmpvr3, tmpvr4, tmpvr5, wght1,
     1    env, rmagy, envinx,sumwat
      integer i1, i2, numres, i, nowres
c
c     + + + LOCAL DEFINITIONS + + +
c     flrcov - surface cover provided by flat residues
c     rcinit - initial residue cover on ridges at harvest
c     rigcvy - residue cover on ridges before shifting today
c     rigcvt - residue cover remaining on ridges after shifting today
c     rmogy  - "residue mass on ground yesterday"; ie, yesterday's value
c              of RMOGT
c     optwat - optimum soil water content; calculated as 60%
c              of the total soil porosity.
c     atemp  - constant used in temperature function; now set to 6.1
c     opttmp - optimum temperature used in function; now set to 33 (C)
c     t1     - (avg. daily temp. [Celsius] + atemp) ** 2
c     t2     - (opttmp [Celsius] + atemp) ** 2 = 1528.81
c     tmpfac - temperature factor for residues (same factor is used for
c                standing, flat, buried residues & roots); 0<tmpfac<1
c     swatfc - water content factor for standing residues; 0<swatfc<1
c     fwatfc - water content factor for flat residues; 0<fwatfc<1
c     ferind - index for soil fertility; now set to 1.0
c     wght1  - fraction of space not in rills (ie, fraction of
c              interrill space).
c     pszind - index for residue particle size influence; now = 1.0
c     i1     - index of operation within tillage sequence
c     i2     - index of tillage sequence
c     numres - number of residue types being tracked on current OFE
c
c --------- Temporary Variables to enhance computational efficiency:
c     envinx - temporary variable used for computing FENVIN & BENVIN
c     tmpvr1 - temporary variable = oratea(iresd(1,iplane))*ferind*pszind
c     tmpvr2 - temporary variable = orater(iresd(1,iplane))*ferind*pszind
c     tmpvr3 - temporary variable = 1.0 - fbrnog(4,iplane)
c     tmpvr4 - temporary variable = rmagy-rmagt(iplane)
c     tmpvr5 - temporary variable = 1.0 - frmove(4,iplane)
c
c     + + + SAVES + + +
c     ... not used at this time 12-15-94 02:04pm sjl
c      save rcinit
c
c     + + + SUBROUTINES CALLED + + +
c     covcal
c
c     + + + DATA INITIALIZATIONS + + +
c
c
c     + + + END SPECIFICATIONS + + +
c
c
c               *********************************
c                *  Begin DECOMPOSITION MODEL  *
c               *********************************
c
c     Initialize local variables
c
      atemp = 6.1
      t2 = 1528.81
      ferind = 1.0
      pszind = 1.0
c
c     Determine number of residue groups currently being tracked
c
      if (iresd(3,iplane).ne.0) then
        numres = 3
      else if (iresd(2,iplane).ne.0) then
        numres = 2
      else
        numres = 1
      end if
c
c     Daily water coefficient for standing residues
c     *************************************************************
c     **    ( WEPP Equation 8.3.5a)                              **
c     **                                                         **
c     **                   Daily Rainfall (mm)                   **
c     **      WCF(s(t)) = ---------------------                  **
c     **                           4                             **
c     **                                                         **
c     **  Assumption: 4 mm of rainfall (per day) will saturate   **
c     **              the standing residues.                     **
c     *************************************************************
c
c     ----- Use daily rainfall to determine saturation of standing residues.
c
      if (tave.gt.0.0 .and. prcp.lt.0.004) then
        swatfc = prcp / 0.004
      elseif (tave.le.0.0) then
        swatfc = 0.0
      else
        swatfc = 1.0
      end if
c
c     Daily water coefficient for 1st soil layer.  (Used
c     for fallen residues, AND buried residues & roots.)
c     Note: This is a corruption of Equations 8.3.5b and
c     8.3.5c. -- CRM -- 3/26/92
c     *************************************************************
c     **    (WEPP Equation 8.3.5b)                               **
c     **                                                         **
c     **                       Residue Water Content             **
c     **      WCF(f(t)) = -------------------------------        **
c     **                   Optimum Residue Water Content         **
c     **                                                         **
c     **                                                         **
c     **    (WEPP Equation 8.3.5c)                               **
c     **                                                         **
c     **  For each                                               **
c     **  soil layer:                                            **
c     **                       Soil Layer Water Content          **
c     **      WCF(g(t)) = ----------------------------------     **
c     **                   Optimum Soil Layer Water Content      **
c     **                                                         **
c     **                                                         **
c     **  NOTE: Ideally, "soil layer water content" and "opti-   **
c     **        mum soil layer water content" are computed as    **
c     **        percent by *weight*.  In the equations below,    **
c     **        they are computed by *volume*.  If the para-     **
c     **        meters become available to compute by weight,    **
c     **        these equations should be converted.             **
c     **        -- CRM -- 11/8/92.                               **
c     **                                                         **
c     **  Currently soil water content is used for both flat     **
c     **  and buried residues because the WEPP WATBAL routines   **
c     **  do not yet track moisture within the surface residue   **
c     **  layer.  When residue moisture is added to WATBAL,      **
c     **  the WCF(f(t)) equation needs to be altered.            **
c     **  == D.E. Stott 6/30/92                                  **
c     **                                                         **
c     **  % Water filled pore space is actually a preferred way  **
c     **  to relate microbial activity to water content of the   **
c     **  soil as it takes bulk density into consideration.      **
c     **  This is especially true for soils that are wet enough  **
c     **  for oxygen to become limiting.  Water potential is     **
c     **  better for dry conditions.  == D.E. Stott 12/31/93     **
c     **                                                         **
c     **  For water contents (%WFP) greater than the optimal     **
c     **  water content (60%), WCF(f(t)) and WCF(g(t)) equal     **
c     **                                                         **
c     **                       Optimum Water Content             **
c     **      WCF(g(t)) = ----------------------------------     **
c     **                        Actual Water Content             **
c     **                                                         **
c     **    Per Linn and Doran (1984)  == D.E. Stott 1/30/93     **
c     **                                                         **
c     *************************************************************
c
c     ----- Calculate the value for 60% water filled pore space.
      optwat = avpor(iplane) * 0.60
c
c     ----- Use soil water content (surface layer) as percent moisture
c     of flat residue.
c

      sumwat=(soilw(1,iplane) + soilw(2,iplane))/
     1       (dg(1,iplane) + dg(2,iplane))

      if(sumwat.le.optwat)then
        fwatfc = sumwat/optwat
      else
        fwatfc = optwat/sumwat
      end if
c
c     Daily temperature index
c     *************************************************************
c     **    (WEPP Equation 8.3.6)                                **
c     **                                                         **
c     **             2*(T(avg)+A)^2 * (T(m)+A)^2 - (T(avg)+A)^4  **
c     **      TCF = -------------------------------------------- **
c     **                              (T(m)+A)^4                 **
c     **                                                         **
c     **  If T(m) = 33 and A = atemp = 6.1, then a temperature   **
c     **    constant, ktemp = 33+6.1 = 39.1.  Therefore:         **
c     **                                                         **
c     **        2*(T(avg)+Atemp)^2 * ktemp^2 - (T(avg)+Atemp)^4  **
c     **   TCF = --------------------------------------------    **
c     **                           ktemp^4                       **
c     **                                                         **
c     **  if (T(avg)+Atemp)^2 = t1 and ktemp^2 = t2, then        **
c     **                                                         **
c     **         (2*t1*t2) -  t1^2              t1[(t2*2)-t1]    **
c     **   TCF = ------------------ ==> TFC = -----------------  **
c     **               t2^2                         t2^2         **
c     **                                                         **
c     *************************************************************
c
c
      if (tave.gt.-6.1.and.tave.lt.49.2) then
        t1 = (tave+atemp) ** 2
c       The following is now set as a constant in the DATA INITIALIZATION
c       statements above.
c       t2 = (opttmp + atemp) ** 2
        tmpfac = t1 * (2*t2-t1) / t2 ** 2
      else
        tmpfac = 0.0
      end if
c
c     **********************************************************************
c     *   NOTE: SENVIN, FENVIN, & BENVIN are set by either *water* or      *
c     *         *temperature* -- whichever is limiting to decomposition.   *
c     *  senvin - daily environmental index for standing residues.         *
c     *  fenvin - daily environmental index for flat residues              *
c     *  benvin - daily environmental index for buried residues and roots  *
c     **********************************************************************
c     Calculate environmental index: separate indices for
c     standing, flat, and buried residues; later calculations
c     for roots use the same index as the buried residues.
c
c     *************************************************************
c     ** (Revised form of Equation 8.3.4, altered to work        **
c     **  with *daily* values of "M" in equation 8.3.7)          **
c     **  -- for standing residues --                            **
c     **                                                         **
c     **      ENVIND(t)   =   Min(TCF, WCF)                      **
c     *************************************************************
      senvin(iplane) = min(tmpfac,swatfc)
c
c     *************************************************************
c     ** (Revised form of Equation 8.3.4, altered to work        **
c     **  with *daily* values of "M" in equation 8.3.7)          **
c     **  -- for flat & buried residues --                       **
c     **  -- for up to 3 residue types. --                       **
c     **                                                         **
c     **      ENVIND(t)   =   Min(TCF, WCF)                      **
c     *************************************************************
      envinx = min(tmpfac,fwatfc)
      do 10 i = 1, numres
        fenvin(i,iplane) = envinx
        benvin(i,iplane) = envinx
   10 continue
c
c
c     ***************************************************************
c     *  rmagt - residue [bio]mass above ground (standing) today    *
c     *  rmogt - residue [bio]mass on the ground (flat on the soil  *
c     *           surface) today                                    *
c     *  smrm  - submerged (buried) residue (from above ground)     *
c     *           [bio]mass today                                   *
c     *  rtm   - non-living root [bio]mass                          *
c     *  rigrm - flat residue mass on ridges                        *
c     *  rilrm - flat residue mass in rills/furrows                 *
c     ***************************************************************
c
c     Compute Changes in Residue and Root Masses.
c
c     ... Compute decomposition of standing material (residue mass
c     standing above ground)  (expanded version of Equation 8.3.7.
c     TMPVR1 takes the place of "ORateA * FI * PS" in Equation 8.3.7.)
c     *************************************************************
c     **  (Corrected form of Equation 8.3.7)                     **
c     **   -- for standing residues --                           **
c     **                                                         **
c     **      M(t) = M(t-1) * exp(-ENVIND(t) * ORateA * FI * PS) **
c     **                                                         **
c     *************************************************************
c
      tmpvr1 = oratea(iresd(1,iplane)) * ferind * pszind
      rmagt(iplane) = rmagt(iplane) * exp(-senvin(iplane)*tmpvr1)
      if (rmagt(iplane).lt.0.00001) rmagt(iplane) = 0.00001
c
c     ----- fraction of interrill area
      wght1 = (rspace(iplane)-width(iplane)) / rspace(iplane)
c
      do 20 nowres = 1, numres
c
c       ...... Compute decomposition of flat residue
c       *************************************************************
c       **  (Corrected form of Equation 8.3.7)                     **
c       **   -- for flat residues --                               **
c       **                                                         **
c       **      M(t) = M(t-1) * exp(-ENVIND(t) * ORateA * FI * PS) **
c       **                                                         **
c       *************************************************************
c
        env = exp(-fenvin(nowres,iplane)*tmpvr1)
c
        rilrm(nowres,iplane) = rilrm(nowres,iplane) * env
        rigrm(nowres,iplane) = rigrm(nowres,iplane) * env
        if (rilrm(nowres,iplane).lt.0.00001) rilrm(nowres,iplane) =
     1      0.00001
        if (rigrm(nowres,iplane).lt.0.00001) rigrm(nowres,iplane) =
     1      0.00001
        rmogt(nowres,iplane) = wght1 * rigrm(nowres,iplane) + (1.0-wght1
     1      ) * rilrm(nowres,iplane)
c
c
c       ...... Compute decomposition of buried mass
c       *************************************************************
c       **  (Corrected form of Equation 8.3.7)                     **
c       **   -- for buried residues --                             **
c       **                                                         **
c       **      M(t) = M(t-1) * exp(-ENVIND(t) * ORateR * FI * PS) **
c       **                                                         **
c       *************************************************************
c
        smrm(nowres,iplane) = smrm(nowres,iplane) * env
        if (smrm(nowres,iplane).lt.0.00001) smrm(nowres,iplane) =
     1      0.00001
c
c
c       ...... Compute decomposition of root mass
c       TMPVR2 takes the place of "ORateR * FI * PS" in Equation 8.3.7.)
c       *************************************************************
c       **  (Corrected form of Equation 8.3.7)                     **
c       **   -- for root residues --                               **
c       **                                                         **
c       **      M(t) = M(t-1) * exp(-ENVIND(t) * ORateA * FI * PS) **
c       **                                                         **
c       *************************************************************
c
        tmpvr2 = orater(iresd(1,iplane)) * ferind * pszind
        env = exp(-benvin(nowres,iplane)*tmpvr2)
        rtm(nowres,iplane) = rtm(nowres,iplane) * env
        if (rtm(nowres,iplane).lt.0.00001) rtm(nowres,iplane) = 0.00001
c
   20 continue
c
c
c     ******************************************************************
c     *  CHANGE IN STANDING & FLAT RESIDUE MASSES DUE TO MANAGEMENT  *
c     ******************************************************************
c
      rmagy = rmagt(iplane)
      rmagt(iplane) = rmagt(iplane) * fact(iresd(1,iplane))
c
c     **************************************************************
c     *  NOTE: In the next 2 statements, it may appear that *all*  *
c     *        the "delta-residue" is being used *twice*; ie, in   *
c     *        *both* RILRM and RIGRM.  However, this is not the   *
c     *        case since the residue masses are expressed *per    *
c     *        unit area*.   CRM -- 10/13/92.                      *
c     **************************************************************
c
c     ---- update ridge & rill residue masses
      rilrm(1,iplane) = rilrm(1,iplane) + (rmagy-rmagt(iplane))
      rigrm(1,iplane) = rigrm(1,iplane) + (rmagy-rmagt(iplane))
c     ---- update residue mass on the ground
      rmogt(1,iplane) = wght1 * rigrm(1,iplane) + (1.0-wght1) *
     1    rilrm(1,iplane)
c
c     *** L1 IF ***
c     Burning
c
      if (sdate.eq.jdburn(nowcrp,iplane).and.jdburn(nowcrp,iplane).ge.
     1    jdharv(nowcrp,iplane)) then
c       ------ fraction of residue mass remaining
        tmpvr3 = 1.0 - fbrnog(nowcrp,iplane)
        do 30 nowres = 1, numres
c         -------- update ridge & rill residue masses
          rigrm(nowres,iplane) = rigrm(nowres,iplane) * tmpvr3
          rilrm(nowres,iplane) = rilrm(nowres,iplane) * tmpvr3
c         -------- update residue mass on the ground
c         (WEPP Equation 8.4.12)
          rmogt(nowres,iplane) = rmogt(nowres,iplane) * tmpvr3
   30   continue
c       ------ update residue mass above the ground
c       (WEPP Equation 8.4.11)
        rmagt(iplane) = rmagt(iplane) * (1.-fbrnag(nowcrp,iplane))
c     ***********************************************************************
c     * NOTE: In the User Documentation AND TILAGE.F (where they are read   *
c     *       in) AND in common block PEREN (where they are stored) FBRNAG  *
c     *       & FBRNOG are BACKWARDS!!  In the equations immediately above, *
c     *       they follow the convention of 'A' designating "above the      *
c     *       ground; ie, standing" and 'O' designating "on the ground; ie, *
c     *       flat". -- CRM -- 11/16/92.                                    *
c     ***********************************************************************
c
c     *** L1 ELSE-IF ***
c     Standing residue cut
c
      else if (sdate.eq.jdcut(nowcrp,iplane)) then
c       ------ "remember" RMAGT
        rmagy = rmagt(iplane)
c       ------ re-compute standing residue mass
c       (WEPP Equation 8.4.13)
        rmagt(iplane) = rmagt(iplane) * (1.-frcut(nowcrp,iplane))
c       ------ compute decrease in standing residue mass
        tmpvr4 = rmagy - rmagt(iplane)
        do 40 nowres = 1, numres
c         -------- add mass of cut residue to ridge & rill residue masses
          rilrm(nowres,iplane) = rilrm(nowres,iplane) + tmpvr4
          rigrm(nowres,iplane) = rigrm(nowres,iplane) + tmpvr4
c         -------- add mass of cut residue to residue mass on the ground
c         (WEPP Equation 8.4.14)
          rmogt(nowres,iplane) = rmogt(nowres,iplane) + tmpvr4
   40   continue
c
c     *** L1 ELSE-IF ***
c     Flat residue removed option modified to be a tillage operation
c     handled below (for management files version 98.3 and up
c
      else if (manver.lt.98.3.and.sdate.eq.jdmove(nowcrp,iplane)) then
c       ------ fraction of flat residue mass remaining
        tmpvr5 = 1.0 - frmove(nowcrp,iplane)
        do 45 nowres = 1, numres
c         -------- update ridge & rill residue masses
          rilrm(nowres,iplane) = rilrm(nowres,iplane) * tmpvr5
          rigrm(nowres,iplane) = rigrm(nowres,iplane) * tmpvr5
c         -------- update residue mass on the ground
c         (WEPP Equation 8.4.15)
          rmogt(nowres,iplane) = rmogt(nowres,iplane) * tmpvr5
   45   continue
c
c     *** L1 END-IF ***
      end if
c
c
c     *********************************************
c     *  CHANGE IN RESIDUE MASS DUE TO TILLAGE  *
c     *********************************************
c
c     ---- index of tillage implement within sequence (scenario)
      i1 = indxy(iplane)
c     ---- index of tillage sequence (scenario)
      i2 = tilseq(nowcrp,iplane)
c
c     *** M1 IF ***
c     Check whether there is tillage on this crop, this OFE.
      if (i2.gt.0.and.i1.gt.0) then
c
c       *** M2 IF ***
c       Check whether this is a management date.
        if (sdate.eq.mdate(i1,i2)) then
          if(manver.lt.98.3.or.resman(i1,i2).le.4.)then
c   ********************************************************************
c   * OLD DECOMP CODE PRIOR TO MULTIPLE RESIDUE ADDITIONS AND REMOVALS *
c   * OR NO RESIDUE REMOVALS OR ADDITIONS TO A NEWER MANAGEMENT FILE   *
c   ********************************************************************
            do 60 nowres = 1, numres
              rmogy(nowres) = rmogt(nowres,iplane)
   60       continue
            rmagy = rmagt(iplane)
c           -------- Calculate standing residue mass today (RMAGT)
c           based on yesterday's value.
c           (WEPP Equation 8.4.6)
            rmagt(iplane) = rmagt(iplane) *
     1          exp(-8.535*mfo(i1,i2,iresd(1,iplane))**2)
c
c           Calculate change in residue mass due to tillage for
c           ridged or non-ridged systems
c
c           (WEPP Equation 8.4.7 applied to ridge and rill areas)
            rigrm(1,iplane) = rigrm(1,iplane) + (rmagy-rmagt(iplane))
            rilrm(1,iplane) = rilrm(1,iplane) + (rmagy-rmagt(iplane))
c
            rmogy(1) = wght1 * rigrm(1,iplane) + (1-wght1) *
     1          rilrm(1,iplane)
c
            do 70 nowres = 1, numres
c             ------------ cover from flat residues (on ridges)
c
              flrcov = (1.0-exp(-cf(iresd(nowres,iplane))*
     1            rigrm(nowres,iplane))) * (1.0-
     1            mfo(i1,i2,iresd(nowres,iplane)))
              if (flrcov.lt.0.0) flrcov = 0.0
              if (flrcov.gt.0.999) flrcov = 0.999
c             ------------ residue mass on ridges
c
              rigrm(nowres,iplane) = log(1-flrcov) / (-
     1            cf(iresd(nowres,iplane)))
c
c             ------------ cover from flat residues (in rills)
              flrcov = (1.0-exp(-cf(iresd(nowres,iplane))*
     1            rilrm(nowres,iplane))) * (1.0-
     1            rmfo(i1,i2,iresd(nowres,iplane)))
              if (flrcov.lt.0.0) flrcov = 0.0
              if (flrcov.gt.0.999) flrcov = 0.999
c
c             ------------ Compute residue mass in rills (RILRM) from flat
c             residue cover (FLRCOV) and flat residue cover coef. (CF).
c             (Revised WEPP Equation 8.3.13)
              rilrm(nowres,iplane) = log(1-flrcov) / (-
     1            cf(iresd(nowres,iplane)))
              rmogt(nowres,iplane) = wght1 * rigrm(nowres,iplane) +
     1            (1.0 - wght1) * rilrm(nowres,iplane)
   70       continue
c
c           Computation of new submerged residue masses due to tillage
c
            do 80 nowres = 1, numres
              smrm(nowres,iplane) = smrm(nowres,iplane) +
     1            rmogy(nowres) - rmogt(nowres,iplane)
   80       continue
          else
c           *********************************************************************
c           * NEW RESIDUE CODE ALLOWING MULTIPLE RESIDUE ADDITIONS AND REMOVALS *
c           *********************************************************************
c
c           resman(mxtill,mxtilsq): residue management option for this operation
c                             10 = residue addition with no disturbance
c                             11 = residue removal with no disturbance
c                             12 = residue addition with disturbance
c                             13 = residue removal with disturbance
c
c           if this is a residue addition date for this OFE
c
            if(resman(i1,i2).eq.10.or.resman(i1,i2).eq.12)then
c
c           update residue tracking
c
c

              do 90 nowres = 1, numres
                rmogy(nowres) = rmogt(nowres,iplane)
   90         continue
              iresd(1,iplane)=iresad(i1,i2)
              call resup(nowcrp,resman(i1,i2))
              if(resman(i1,i2).eq.10)
     1            write(6,1000)resad(i1,i2),crpnam(iresd(1,iplane)),
     1                iplane,mdate(i1,i2)
              if(resman(i1,i2).eq.12)
     1            write(6,1010)resad(i1,i2),crpnam(iresd(1,iplane)),
     1            iplane,mdate(i1,i2)
              rmogt(1,iplane)=resad(i1,i2)
              rilrm(1,iplane)=rmogt(1,iplane)
              rigrm(1,iplane)=rmogt(1,iplane)
c
c           if this is a residue removal date for this OFE
            else if(resman(i1,i2).eq.11.or.resman(i1,i2).eq.13)then
c           amount of residue removed if any
              if(resman(i1,i2).eq.11)
     1            write(6,1050)frmove(i1,i2)*100,iplane,mdate(i1,i2)
              if(resman(i1,i2).eq.13)
     1            write(6,1100)frmove(i1,i2)*100,iplane,mdate(i1,i2)

c             ------ fraction of flat residue mass remaining
              tmpvr5 = 1.0 - frmove(i1,i2)
              do 100 nowres = 1, numres
c           -------- update ridge & rill residue masses
                rmogt(nowres,iplane) = rmogt(nowres,iplane) * tmpvr5
                rilrm(nowres,iplane) = rilrm(nowres,iplane) * tmpvr5
                rigrm(nowres,iplane) = rigrm(nowres,iplane) * tmpvr5
c               -------- update residue mass on the ground
c               (WEPP Equation 8.4.15)
  100         continue
            end if
c
          rmagy = rmagt(iplane)
c         -------- Calculate standing residue mass today (RMAGT)
c         based on yesterday's value.
c         (WEPP Equation 8.4.6)
c
c         if residue removal with disturbance
          if(resman(i1,i2).eq.13)
     1        rmagt(iplane) = rmagt(iplane) *
     1        exp(-8.535*mfo(i1,i2,iresd(1,iplane))**2)
c
          if(resman(i1,i2).eq.12)
     1        rmagt(iplane) = rmagt(iplane) *
     1        exp(-8.535*mfo(i1,i2,nowcrp)**2)
c
c         Calculate change in residue mass due to tillage for
c         ridged or non-ridged systems
c
c         (WEPP Equation 8.4.7 applied to ridge and rill areas)
          rigrm(1,iplane) = rigrm(1,iplane) + (rmagy-rmagt(iplane))
          rilrm(1,iplane) = rilrm(1,iplane) + (rmagy-rmagt(iplane))
          rmogy(1) = wght1 * rigrm(1,iplane) + (1-wght1) *
     1        rilrm(1,iplane)

c
c
c         Computation of new submerged residue masses due to tillage
c
          do 110 nowres = 1, numres
            if(resman(i1,i2).eq.10.or.resman(i1,i2).eq.11)then
              smrm(nowres,iplane) = smrm(nowres,iplane)
            else
              smrm(nowres,iplane) = smrm(nowres,iplane)
     1            + rmogy(nowres) - rmogt(nowres,iplane)
              if(smrm(nowres,iplane).lt.0.00001)
     1              smrm(nowres,iplane)=0.00001
            end if
  110     continue
cc       *** M2 END-IF ***
         end if
c        *** M1 END-IF ***
       end if
c
       end if
c
c XXX  Problem here with residue moving from interrill to rill area
c      when rills are narrow - get unrealistic amounts of residue
c      piling up in the rills.   Commented out this code  (N1 IF)
c      until problem is resolved, based on Nearing recommendation.
c      dcf  8/25/94
c
c     ******************************************************************
c     *  RESIDUE MOVEMENT FROM RIDGE TO FURROW IN RIDGE-TILLAGE SYSTEMS
c     ******************************************************************
c
c     *** N1 IF's (Note check for positive subscript.) ***
c     if (i2.gt.0) then
c       if (iridge(i2).eq.1) then
cc        ------ "days after harvest"
c         dah = sdate - jdharv(nowcrp,iplane)
cc        *** N2 IF ***
cc        If this is a harvest date....
c
cc        NOTE - change code to if(dah.eq.1) because change in order of
cc        subroutine calls in SR CONTIN has messed this up ----  we call
cc        DECOMP at the beginning of the simulation day before WATBAL
cc        and PTGRA and RESUP have been called - since RESUP has not yet
cc        been called on an sdate=jdharv  we have not yet incremented
cc        iresd by 1, so on dah=0  numres may only equal 1 when it
cc        should be 2.
c         if (dah.eq.1) then
c           do 90 nowres = 1, numres
cc            ---------- compute ridge cover at harvest, from residue
cc            mass on ridges
cc            (Note: "1 - exp(-cf)" converts from *mass* to *cover*.)
c             rcinit(nowres,iplane) = 1 -
c    1            exp(-cf(iresd(nowres,iplane))*rigrm(nowres,iplane))
c  90       continue
cc        *** N2 ELSE-IF ***
cc        If it's less than 60 days after harvest, calculate residue
cc        redistribution from ridges to furrows by wind, rain & gravity.
c         else if (dah.gt.1.and.dah.le.60) then
c           do 100 nowres = 1, numres
cc            ---------- residue cover on ridges before redistribution
c             rigcvy = 1 - exp(-cf(iresd(nowres,iplane))*
c    1            rigrm(nowres,iplane))
cc            ---------- residue cover remaining on the ridges today
cc            (prorated over the 60 days following harvest, down to
cc            30 percent).
cc            (WEPP Equation 8.3.12)
c             rigcvt = rcinit(nowres,iplane) - (rcinit(nowres,iplane)-
c    1            0.30) / 60.0 * dah
cc            ---------- mass of residue moved from the ridges into the
cc            furrows today
cc            (WEPP Equation 8.3.13)
c             deltrm = log(1.-(rigcvy-rigcvt)) / (-
c    1            cf(iresd(nowres,iplane)))
c             if (deltrm.lt.0.0) deltrm = 0.0
cc            ---------- residue mass in the rills
cc            (WEPP Equation 8.3.14)
c             rilrm(nowres,iplane) = rilrm(nowres,iplane) + deltrm
cc            ---------- residue mass on the ridges
cc            (WEPP Equation 8.3.16)
c             rigrm(nowres,iplane) = rigrm(nowres,iplane) - deltrm
c             rmogt(nowres,iplane) = wght1 * rigrm(nowres,iplane) + (1.0
c    1            -wght1) * rilrm(nowres,iplane)
c 100       continue
cc        *** N2 END-IF ***
c         end if
cc      *** N1 END-IF's ***
c       end if
c     end if
c XXX End of commented out code for problem with residue movement
c     from interrill to rills for ridge tillage systems. dcf 8/25/94
c
c
c     --- Compute the fraction of flat surface residue cover, standing
c     residue cover, and average surface cover
      call covcal(iresd,iplane)
c
c     Compute (POP) population of standing stalks(?) from the ratio of
c     standing residue mass today (RMAGT) to standing residue at
c     harvest (SRMHAV), and population at maturity (POPMAT).
c
c       added trap to prevent divide by zero if user has input a
c       value of zero for cutting height (making partcf = 0.0 )
c       dcf   3/8/94
c
      if(srmhav(iplane) .gt. 0.0)then
        pop(iplane) = (rmagt(iplane)/srmhav(iplane)) * popmat(iplane)
      else
        pop(iplane) = 0.0
      endif
c
      if (pop(iplane).gt.popmat(iplane)) pop(iplane) = popmat(iplane)
c
c     -- XXX -- This computation of POP looks like a real problem to me!
c     The approach seems unsound -- POP will max out before
c     maturity is reached.  The location of this code seems
c     inappropriate -- the code has nothing to do with
c     computing residue or decomposition.  I suggest it be
c     moved to the rangeland (or perennial) plant growth
c     routines.  -- CRM -- 12/01/92.
c
      return
 1000 format(1x,f4.1,'(kg/m^2) ',A10,' RESIDUE ADDED TO PLANE '
     1    ,i2,', NO DISTURBANCE, ON DAY',i4)
 1010 format(1x,f4.1,'(kg/m^2) ',A10,' RESIDUE ADDED TO PLANE '
     1    ,i2,', WITH DISTURBANCE, ON DAY',i4)
 1050 format(1x,f4.1,'% RESIDUE REMOVED FROM PLANE '
     1    ,i2,', NO DISTURBANCE, ON DAY',i4)
 1100 format(1x,f4.1,'% RESIDUE REMOVED FROM PLANE '
     1    ,i2,', WITH DISTURBANCE, ON DAY',i4)
      end

ext: for

max_stars_repo_name: jarad/dep

max_issues_repo_name: jarad/dep

max_forks_repo_name: akrherz/idep


text:
program main
  use iso_c_binding
  use blopex_fortran_driver
  implicit none
  integer(4) :: N, NZ
  integer(4) :: n_eigs, maxit, loglevel, prec_type, n_bc
  integer(4), allocatable :: index(:), item(:), i_bc(:)
  real(8) :: tol, t1, t2, t3, t4
  real(8), allocatable :: A(:), B(:), eigen_val(:), eigen_vec(:,:)
  character :: finA*100, finB*100, finBC*100
  logical :: is_B = .false.
  logical :: is_BC = .false.

  write(*,"(a)")"* blopex fortran driver"

  call cpu_time(t1)

  call get_input_arg(finA, finB, finBC, n_eigs, maxit, loglevel, tol, is_B, is_BC)

  if(is_BC)then
    call input_BC(finBC, n_bc, i_bc)
    call input_from_matrix_market_csr_bc(finA, N, NZ, index, item, A, n_bc, i_bc)
    if(is_B) call input_from_matrix_B_bc(finB, N, B, n_bc, i_bc)
  else
    call input_from_matrix_market_csr(finA, N, NZ, index, item, A)
    if(is_B) call input_from_matrix_B(finB, N, B)
  endif

  allocate(eigen_val(n_eigs), source = 0.0d0)
  allocate(eigen_vec(N,n_eigs), source = 0.0d0)

  call cpu_time(t2)

  if(is_B)then
    call blopex_lobpcg_solve(N, NZ, index, item, A, n_eigs, maxit, tol, loglevel, &
      & eigen_val, eigen_vec, prec_type, B)
  else
    call blopex_lobpcg_solve(N, NZ, index, item, A, n_eigs, maxit, tol, loglevel, &
      & eigen_val, eigen_vec, prec_type)
  endif

  call cpu_time(t3)

  if(is_BC)then
    call convet_eigval(N, n_bc, i_bc, n_eigs, eigen_vec)
    call output(N, n_eigs, eigen_val, eigen_vec)
  else
    call output(N, n_eigs, eigen_val, eigen_vec)
  endif

  call cpu_time(t4)

  write(*,"(a,1pe12.5)") "* input  time:", t2-t1
  write(*,"(a,1pe12.5)") "* LOBPCG time:", t3-t2
  write(*,"(a,1pe12.5)") "* output time:", t4-t3
  write(*,"(a,1pe12.5)") "* total  time:", t4-t1

contains

  subroutine convet_eigval(N, n_bc, i_bc, n_eigs, eigen_vec)
    implicit none
    integer(4) :: N, n_bc, i_bc(:), N_org, n_eigs, i, j, in, ip
    real(8), allocatable :: eigen_vec(:,:)
    integer(4), allocatable :: perm(:)
    real(8), allocatable :: temp(:,:)
    logical, allocatable :: is_use(:)

    allocate(temp(N,n_eigs), source = 0.0d0)
    temp = eigen_vec
    deallocate(eigen_vec)

    N_org = N + n_bc
    allocate(eigen_vec(N_org,n_eigs), source = 0.0d0)

    allocate(is_use(N_org), source = .true.)
    do i = 1, n_bc
      ip = i_bc(i)
      is_use(ip) = .false.
    enddo

    allocate(perm(N_org), source = 0)
    in = 1
    do i = 1, N_org
      if(.not. is_use(i)) cycle
      perm(i) = in
      in = in + 1
    enddo

    do j = 1, n_eigs
      do i = 1, N_org
        ip = perm(i)
        if(ip == 0) cycle
        eigen_vec(i,j) = temp(ip,j)
      enddo
    enddo

    N = N_org
  end subroutine convet_eigval

  subroutine output(N, n_eigs, eigen_val, eigen_vec)
    implicit none
    integer(4) :: N, n_eigs, i, j
    real(8) :: eigen_val(:), eigen_vec(:,:)
    character :: cnum*4

    call system('if [ ! -d output ]; then (echo "* create output"; mkdir -p output); fi')

    write(*,"(a)")"* eigen value"
    write(*,"(1pe12.5)")eigen_val

    open(20, file = "output/eigen_value.txt", status = "replace")
    write(20,"(i4)")n_eigs
    do i = 1, n_eigs
      write(20, "(i4,1pe12.4)")i, eigen_val(i)
    enddo
    close(20)

    do i = 1, n_eigs
      write(cnum,"(i0)")i
      open(20, file = "output/eigen_vector."//trim(cnum)//".txt", status = "replace")
      write(20,"(i8)")N
      do j = 1, N
        write(20, "(1pe12.4)")eigen_vec(j,i)
      enddo
      close(20)
    enddo
  end subroutine output

  subroutine get_input_arg(finA, finB, finBC, n_eigs, maxit, loglevel, tol, is_B, is_BC)
    implicit none
    integer(4) :: count
    integer(4) :: n_eigs, maxit, loglevel, prec
    real(8) :: tol
    character :: finA*100, finB*100, finBC*100
    logical :: is_B, is_BC

    count = iargc()
    if(count == 1)then
      call getarg(1, finA)
    elseif(count == 2)then
      call getarg(1, finA)
      call getarg(2, finBC)
      is_BC = .true.
    elseif(count == 3)then
      call getarg(1, finA)
      call getarg(2, finBC)
      call getarg(3, finB)
      is_BC = .true.
      is_B = .true.
    else
      stop "please enter input file name of matrix A and matrix B (optional)"
    endif

    open(10, file = "setting.dat", status = "old")
      read(10,*)n_eigs
      read(10,*)maxit
      read(10,*)tol
      read(10,*)loglevel
      read(10,*)prec_type
    close(10)

    write(*,"(a,i12)")     "* n_eigs  :", n_eigs
    write(*,"(a,i12)")     "* maxiter :", maxit
    write(*,"(a,1pe12.5)") "* tol     :", tol
    write(*,"(a,i12)")     "* loglevel:", loglevel
    write(*,"(a,l12)")     "* read BC :", is_BC
    write(*,"(a,l12)")     "* read B  :", is_B
    write(*,"(a,i12)")     "* apply M :", prec_type
  end subroutine get_input_arg

  subroutine input_from_matrix_market_mm(fin, N, A)
    implicit none
    integer(4) :: N, NZ, i, j, in, ierr
    real(8) :: val
    real(8), allocatable :: A(:,:)
    character :: fin*100, ctemp*10
    logical :: is_first

    is_first = .true.
    open(20, file = trim(fin), status = "old")
    in = 0
    do
      read(20, "(a)", iostat = ierr) ctemp
      if(0 /= ierr) exit

      if(ctemp(1:1) == "%") cycle

      if(is_first)then
        backspace(20)
        read(20,*) N, i, NZ
        allocate(A(N,N), source = 0.0d0)
        is_first = .false.
      else
        in = in + 1
        backspace(20)
        read(20,*) i, j, val
        A(i,j) = val
        A(j,i) = val
      endif
      if(in == NZ) exit
    enddo
    close(20)
  end subroutine input_from_matrix_market_mm

  subroutine input_from_matrix_market_csr(fin, N, NZ, index, item, A)
    implicit none
    integer(4) :: N, NZ, i, j, in
    integer(4), allocatable :: index(:), item(:)
    real(8), allocatable :: A(:), tmp(:,:)
    character :: fin*100

    call input_from_matrix_market_mm(fin, N, tmp)

    NZ = 0
    do i = 1, N
      do j = 1, N
        if(dabs(tmp(j,i)) > 1.0d-10) NZ = NZ + 1
      enddo
    enddo

    allocate(index(0:N), source = 0)
    allocate(item(NZ), source = 0)
    allocate(A(NZ), source = 0.0d0)

    in = 0
    do i = 1, N
      do j = 1, N
        if(dabs(tmp(j,i)) > 1.0d-10)then
          in = in + 1
          index(i) = index(i) + 1
          item(in) = j
          A(in) = tmp(j,i)
        endif
      enddo
    enddo

    do i = 1, N
      index(i) = index(i) + index(i-1)
    enddo
  end subroutine input_from_matrix_market_csr

  subroutine input_from_matrix_market_csr_bc(fin, N, NZ, index, item, A, n_bc, i_bc)
    implicit none
    integer(4) :: N, NZ, i, j, in, n_bc, ip, jp
    integer(4) :: N_new
    integer(4) :: i_bc(:)
    integer(4), allocatable :: index(:), item(:), perm(:)
    real(8), allocatable :: A(:), tmp(:,:)
    logical, allocatable :: is_use(:)
    character :: fin*100

    call input_from_matrix_market_mm(fin, N, tmp)

    allocate(is_use(N), source = .true.)
    do i = 1, n_bc
      in = i_bc(i)
      is_use(in) = .false.
    enddo

    N_new = N - n_bc
    allocate(perm(N), source = 0)
    in = 1
    do i = 1, N
      if(.not. is_use(i)) cycle
      perm(i) = in
      in = in + 1
    enddo

    NZ = 0
    do i = 1, N
      if(.not. is_use(i)) cycle
      do j = 1, N
        if(.not. is_use(j)) cycle
        if(dabs(tmp(j,i)) > 1.0d-10) NZ = NZ + 1
      enddo
    enddo

    allocate(index(0:N_new), source = 0)
    allocate(item(NZ), source = 0)
    allocate(A(NZ), source = 0.0d0)

    in = 0
    do i = 1, N
      if(.not. is_use(i)) cycle
      do j = 1, N
        if(.not. is_use(j)) cycle
        if(dabs(tmp(j,i)) > 1.0d-10)then
          in = in + 1
          ip = perm(i)
          index(ip) = index(ip) + 1
          jp = perm(j)
          item(in) = jp
          A(in) = tmp(j,i)
        endif
      enddo
    enddo

    do i = 1, N_new
      index(i) = index(i) + index(i-1)
    enddo

    N = N_new
  end subroutine input_from_matrix_market_csr_bc

  subroutine input_BC(fin, n_bc, i_bc)
    implicit none
    integer(4) :: n_bc, i, n, ndof, id, idof
    integer(4), allocatable :: i_bc(:)
    real(8) :: val
    character :: fin*100

    open(20, file = trim(fin), status = "old")
      read(20,*)n_bc, ndof
      allocate(i_bc(n_bc), source = 0)
      do i = 1, n_bc
        read(20,*) id, idof, val
        i_bc(i) = ndof*(id-1) + idof
      enddo
    close(20)
  end subroutine input_BC

  subroutine input_from_matrix_B_bc(fin, N, B, n_bc, i_bc)
    implicit none
    integer(4) :: N, L, i, n_bc, i_bc(:), in, ip
    real(8) :: val
    integer(4), allocatable :: perm(:)
    real(8), allocatable :: B(:)
    logical, allocatable :: is_use(:)
    character :: fin*100

    allocate(is_use(N+n_bc), source = .true.)
    do i = 1, n_bc
      in = i_bc(i)
      is_use(in) = .false.
    enddo

    allocate(perm(N+n_bc), source = 0)
    in = 1
    do i = 1, N+n_bc
      if(.not. is_use(i)) cycle
      perm(i) = in
      in = in + 1
    enddo

    allocate(B(N), source = 0.0d0)

    open(20, file = trim(fin), status = "old")
      read(20,*)L
      if(L /= N+n_bc) stop "input_from_matrix_B_bc: "
      do i = 1, L
        read(20,*) val
        ip = perm(i)
        if(ip == 0) cycle
        B(ip) = val
      enddo
    close(20)
  end subroutine input_from_matrix_B_bc

  subroutine input_from_matrix_B(fin, N, B)
    implicit none
    integer(4) :: N, i
    real(8), allocatable :: B(:)
    character :: fin*100

    open(20, file = trim(fin), status = "old")
      read(20,*)i
      if(i /= N) stop "input_from_matrix_B: "
      allocate(B(N), source = 0.0d0)
      do i = 1, N
        read(20,*) B(i)
      enddo
    close(20)
  end subroutine input_from_matrix_B
end program main

max_stars_repo_name: nqomorita/blopex

max_issues_repo_name: nqomorita/blopex

max_forks_repo_name: nqomorita/blopex


text:
*dk,try2to3i
      subroutine try2to3i(it1,ifpos1,it2,ibdytet,iflip,
     *                    npoints,ntets)
C
C ######################################################################
C
C     PURPOSE -
C
C        This routine checks whether it is feasible to flip a
C        connection on a material interface.  To consider such a flip,
C        the two tetrahedra, it1 and it2, must each have a boundary
C        face along the same connection, and that connection must be
C        only a two-material connection.  If the flip is feasible, that
C        boundary connection will be changed to an interior connection,
C        and a new boundary connection will be drawn between the two
C        points opposite the face common to it1 and it2.  The flip is
C        actually performed by the routine flp2to3i.
C
C     INPUT ARGUMENTS -
C
C        it1      - the first tet
C        ifpos1   - the position of the common face in it1
C
C     OUTPUT ARGUMENTS -
C
C        it2      - the second tet which shares the common face
C        ibdytet  - tells which of the three potential new tets will be
C                   the boundary tet
C        iflip    - 0 => flip is not indicated
C                   1 => flip is indicated
C
C     CHANGE HISTORY -
C
C        $Log: try2to3i.f,v $
C        Revision 2.00  2007/11/09 20:04:05  spchu
C        Import to CVS
C
CPVCS
CPVCS       Rev 1.7   05 Jan 2001 12:58:12   dcg
CPVCS    use dble in place of float so as not to lose digits
CPVCS
CPVCS       Rev 1.6   Mon Jul 14 15:49:32 1997   dcg
CPVCS    fix argument in call to getchain - previous argument
CPVCS    was imt1 which caused the mesh object attribute imt1
CPVCS    to be overwritten - changed to local array imt0
CPVCS
CPVCS       Rev 1.5   Mon Apr 14 17:05:08 1997   pvcs
CPVCS    No change.
CPVCS
CPVCS       Rev 1.4   12/02/94 16:19:20   dcg
CPVCS     changes required to compile on IBM RS6000
CPVCS
CPVCS       Rev 1.3   12/02/94 15:07:10   het
CPVCS    Added an option for the "cmo" access functions
CPVCS
CPVCS
CPVCS       Rev 1.2   12/01/94 18:49:44   het
CPVCS    Added a data type to the "cmo" calles
CPVCS       and added the "cmo.h" include file.
CPVCS
CPVCS       Rev 1.1   11/17/94 21:55:32   het
CPVCS    Added include files for chydro, neibor, cmerge, comdict. Added calles and
CPVCS    pointer statements for current_mesh_object database access.
CPVCS
CPVCS       Rev 1.0   11/10/94 12:19:50   pvcs
CPVCS    Original version.
C
C ######################################################################
C
      implicit none
C
      include "cmo.h"
      include "chydro.h"
      include "neibor.h"

C arguments (it1,ifpos1,it2,ibdytet,iflip,npoints,ntets)
      integer it1,ifpos1,it2,ibdytet,iflip,npoints,ntets

C variables

      integer nmulti
      parameter (nmulti = 200)

      integer kpts(3),ichain1(nmulti),imt0(nmulti)

      integer ierror,length,icmotype,lenitetclr,ier,lenitet,
     *        lenjtet,m,kpos1,kpt,kpos2,k,ktest,ict1,ipar

C     MACROS function statements
C
      integer iposfnd,i0,i1,i2,i3,i4
      iposfnd(i0,i1,i2,i3,i4)=min0(iabs(i1-i0)*4+1,iabs(i2-i0)*4+2,
     *                             iabs(i3-i0)*4+3,iabs(i4-i0)*4+4)

C
C ######################################################################
C BEGIN begin
C
C
C     ******************************************************************
C     FETCH MESH OBJECT NAME AND POINTER INFORMATION.
C
      if(icmoget.eq.1) then
C
      call cmo_get_name(cmo,ierror)
C
      call cmo_get_info('mbndry',cmo,mbndry,length,icmotype,ierror)
      call cmo_get_info('itetclr',cmo,ipitetclr,lenitetclr,icmotype,ier)
      call cmo_get_info('itet',cmo,ipitet,lenitet,icmotype,ierror)
      call cmo_get_info('jtet',cmo,ipjtet,lenjtet,icmotype,ierror)
C
      endif
C
C     ******************************************************************
C
C
      iflip=0
      it2=0.25*dble(jtet(ifpos1,it1))+0.9
C
C     ******************************************************************
C
C     DETERMINE IF TWO BOUNDARY FACES SHARE THE SAME CONNECTION.
C
      do 500 m=1,4
         if (jtet(m,it1) .le. mbndry) goto 500
         kpos1=m
         kpt=itet(kpos1,it1)
         kpos2=iposfnd(kpt,itet(1,it2),itet(2,it2),itet(3,it2),
     *                 itet(4,it2))
         if (jtet(kpos2,it2) .le. mbndry) goto 500
C
C        ***************************************************************
C
C        VERIFY THAT THE CONNECTION IS A TWO-MATERIAL CONNECTION.
C
         kpts(1)=itet(iflist(3*ifpos1-2),it1)
         kpts(3)=itet(iflist(3*ifpos1-1),it1)
         kpts(2)=itet(iflist(3*ifpos1  ),it1)
         do 300 k=1,3
            ktest=kpts(k)
            if (ktest .eq. kpt) then
C
C              .........................................................
C              DETERMINE WHICH OF THE THREE POTENTIAL NEW TETS WILL
C              BE THE BOUNDARY TET.
C
               if (k .eq. 1) then
                  ibdytet=2
               elseif (k .eq. 2) then
                  ibdytet=3
               else
                 ibdytet=1
               endif
               goto 300
            endif
            call getchain(ktest,ichain1,imt0,nmulti,ict1,ipar)
            if (ict1 .gt. 1)  goto 500
 300     continue
         iflip=1
         goto 9999
 500  continue
C
      goto 9999
 9999 continue
      end

ext: f

max_forks_repo_name: daniellivingston/LaGriT


text:
      SUBROUTINE SSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
*
*  -- LAPACK routine (version 1.0b) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     October 31, 1992
*
*     .. Scalar Arguments ..
      CHARACTER          COMPZ
      INTEGER            INFO, LDZ, N
*     ..
*     .. Array Arguments ..
      REAL               D( * ), E( * ), WORK( * ), Z( LDZ, * )
*     ..
*
*  Purpose
*  =======
*
*  SSTEQR computes all eigenvalues and, optionally, eigenvectors of a
*  symmetric tridiagonal matrix using the implicit QL or QR method.
*  The eigenvectors of a full or band symmetric matrix can also be found
*  if SSYTRD or SSPTRD or SSBTRD has been used to reduce this matrix to
*  tridiagonal form.
*
*  Arguments
*  =========
*
*  COMPZ   (input) CHARACTER*1
*          Specifies whether eigenvectors are to be computed
*          as follows
*
*             COMPZ = 'N' or 'n'     Compute eigenvalues only.
*
*             COMPZ = 'V' or 'v'     Compute eigenvectors of original
*                                    symmetric matrix also.
*                                    Array Z contains the orthogonal
*                                    matrix used to reduce the original
*                                    matrix to tridiagonal form.
*
*             COMPZ = 'I' or 'i'     Compute eigenvectors of
*                                    tridiagonal matrix also.
*
*  N       (input) INTEGER
*          The number of rows and columns in the matrix.  N >= 0.
*
*  D       (input/output) REAL array, dimension (N)
*          On entry, the diagonal elements of the tridiagonal matrix.
*          On exit, if INFO = 0, the eigenvalues in ascending order.
*
*  E       (input/output) REAL array, dimension (N-1)
*          On entry, the subdiagonal elements of the tridiagonal matrix.
*          On exit, E has been destroyed.
*
*  Z       (input/output) REAL array, dimension (LDZ, N)
*          If  COMPZ = 'V' or 'v', then:
*          On entry, Z contains the orthogonal matrix used in the
*          reduction to tridiagonal form.
*          If  COMPZ = 'V' or 'v' or 'I' or 'i', then:
*          On exit, if INFO = 0, Z contains the orthonormal eigenvectors
*          of the symmetric tridiagonal (or full) matrix.
*
*          If  COMPZ = 'N' or 'n', then Z is not referenced.
*
*  LDZ     (input) INTEGER
*          The leading dimension of the array Z.  If eigenvectors are
*          desired, then  LDZ >= max( 1, N ).  In any case,  LDZ >= 1.
*
*  WORK    (workspace) REAL array, dimension (max(1,2*N-2))
*          Workspace used in computing eigenvectors.
*          If  COMPZ = 'N' or 'n', then WORK is not referenced.
*
*  INFO    (output) INTEGER
*          = 0:  successful exit.
*          < 0:  if INFO = -i, the i-th argument had an illegal value.
*          > 0:  the algorithm has failed to find all the eigenvalues in
*                a total of 30*N iterations; if INFO = i, then i
*                elements of E have not converged to zero; on exit, D
*                and E contain the elements of a symmetric tridiagonal
*                matrix which is orthogonally similar to the original
*                matrix.
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO, ONE, TWO
      PARAMETER          ( ZERO = 0.0, ONE = 1.0, TWO = 2.0 )
      INTEGER            MAXIT
      PARAMETER          ( MAXIT = 30 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, ICOMPZ, II, J, JTOT, K, L, L1, LEND, LENDM1,
     $                   LENDP1, LM1, M, MM, MM1, NM1, NMAXIT
      REAL               B, C, EPS, F, G, P, R, RT1, RT2, S, TST
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      REAL               SLAMCH, SLAPY2
      EXTERNAL           LSAME, SLAMCH, SLAPY2
*     ..
*     .. External Subroutines ..
      EXTERNAL           SLAE2, SLAEV2, SLARTG, SLASR, SLAZRO, SSWAP,
     $                   XERBLA
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, MAX, SIGN
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      INFO = 0
*
      IF( LSAME( COMPZ, 'N' ) ) THEN
         ICOMPZ = 0
      ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
         ICOMPZ = 1
      ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
         ICOMPZ = 2
      ELSE
         ICOMPZ = -1
      END IF
      IF( ICOMPZ.LT.0 ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      ELSE IF( ( LDZ.LT.1 ) .OR. ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1,
     $         N ) ) ) THEN
         INFO = -6
      END IF
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'SSTEQR', -INFO )
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( N.EQ.0 )
     $   RETURN
*
      IF( N.EQ.1 ) THEN
         IF( ICOMPZ.GT.0 )
     $      Z( 1, 1 ) = ONE
         RETURN
      END IF
*
*     Determine the unit roundoff for this environment.
*
      EPS = SLAMCH( 'E' )
*
*     Compute the eigenvalues and eigenvectors of the tridiagonal
*     matrix.
*
      IF( ICOMPZ.EQ.2 )
     $   CALL SLAZRO( N, N, ZERO, ONE, Z, LDZ )
*
      NMAXIT = N*MAXIT
      JTOT = 0
*
*     Determine where the matrix splits and choose QL or QR iteration
*     for each block, according to whether top or bottom diagonal
*     element is smaller.
*
      L1 = 1
      NM1 = N - 1
*
   10 CONTINUE
      IF( L1.GT.N )
     $   GO TO 160
      IF( L1.GT.1 )
     $   E( L1-1 ) = ZERO
      IF( L1.LE.NM1 ) THEN
         DO 20 M = L1, NM1
            TST = ABS( E( M ) )
            IF( TST.LE.EPS*( ABS( D( M ) )+ABS( D( M+1 ) ) ) )
     $         GO TO 30
   20    CONTINUE
      END IF
      M = N
*
   30 CONTINUE
      L = L1
      LEND = M
      IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN
         L = LEND
         LEND = L1
      END IF
      L1 = M + 1
*
      IF( LEND.GE.L ) THEN
*
*        QL Iteration
*
*        Look for small subdiagonal element.
*
   40    CONTINUE
         IF( L.NE.LEND ) THEN
            LENDM1 = LEND - 1
            DO 50 M = L, LENDM1
               TST = ABS( E( M ) )
               IF( TST.LE.EPS*( ABS( D( M ) )+ABS( D( M+1 ) ) ) )
     $            GO TO 60
   50       CONTINUE
         END IF
*
         M = LEND
*
   60    CONTINUE
         IF( M.LT.LEND )
     $      E( M ) = ZERO
         P = D( L )
         IF( M.EQ.L )
     $      GO TO 80
*
*        If remaining matrix is 2-by-2, use SLAE2 or SLAEV2
*        to compute its eigensystem.
*
         IF( M.EQ.L+1 ) THEN
            IF( ICOMPZ.GT.0 ) THEN
               CALL SLAEV2( D( L ), E( L ), D( L+1 ), RT1, RT2, C, S )
               WORK( L ) = C
               WORK( N-1+L ) = S
               CALL SLASR( 'R', 'V', 'B', N, 2, WORK( L ),
     $                     WORK( N-1+L ), Z( 1, L ), LDZ )
            ELSE
               CALL SLAE2( D( L ), E( L ), D( L+1 ), RT1, RT2 )
            END IF
            D( L ) = RT1
            D( L+1 ) = RT2
            E( L ) = ZERO
            L = L + 2
            IF( L.LE.LEND )
     $         GO TO 40
            GO TO 10
         END IF
*
         IF( JTOT.EQ.NMAXIT )
     $      GO TO 140
         JTOT = JTOT + 1
*
*        Form shift.
*
         G = ( D( L+1 )-P ) / ( TWO*E( L ) )
         R = SLAPY2( G, ONE )
         G = D( M ) - P + ( E( L ) / ( G+SIGN( R, G ) ) )
*
         S = ONE
         C = ONE
         P = ZERO
*
*        Inner loop
*
         MM1 = M - 1
         DO 70 I = MM1, L, -1
            F = S*E( I )
            B = C*E( I )
            CALL SLARTG( G, F, C, S, R )
            IF( I.NE.M-1 )
     $         E( I+1 ) = R
            G = D( I+1 ) - P
            R = ( D( I )-G )*S + TWO*C*B
            P = S*R
            D( I+1 ) = G + P
            G = C*R - B
*
*           If eigenvectors are desired, then save rotations.
*
            IF( ICOMPZ.GT.0 ) THEN
               WORK( I ) = C
               WORK( N-1+I ) = -S
            END IF
*
   70    CONTINUE
*
*        If eigenvectors are desired, then apply saved rotations.
*
         IF( ICOMPZ.GT.0 ) THEN
            MM = M - L + 1
            CALL SLASR( 'R', 'V', 'B', N, MM, WORK( L ), WORK( N-1+L ),
     $                  Z( 1, L ), LDZ )
         END IF
*
         D( L ) = D( L ) - P
         E( L ) = G
         GO TO 40
*
*        Eigenvalue found.
*
   80    CONTINUE
         D( L ) = P
*
         L = L + 1
         IF( L.LE.LEND )
     $      GO TO 40
         GO TO 10
*
      ELSE
*
*        QR Iteration
*
*        Look for small superdiagonal element.
*
   90    CONTINUE
         IF( L.NE.LEND ) THEN
            LENDP1 = LEND + 1
            DO 100 M = L, LENDP1, -1
               TST = ABS( E( M-1 ) )
               IF( TST.LE.EPS*( ABS( D( M ) )+ABS( D( M-1 ) ) ) )
     $            GO TO 110
  100       CONTINUE
         END IF
*
         M = LEND
*
  110    CONTINUE
         IF( M.GT.LEND )
     $      E( M-1 ) = ZERO
         P = D( L )
         IF( M.EQ.L )
     $      GO TO 130
*
*        If remaining matrix is 2-by-2, use SLAE2 or SLAEV2
*        to compute its eigensystem.
*
         IF( M.EQ.L-1 ) THEN
            IF( ICOMPZ.GT.0 ) THEN
               CALL SLAEV2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2, C, S )
               WORK( M ) = C
               WORK( N-1+M ) = S
               CALL SLASR( 'R', 'V', 'F', N, 2, WORK( M ),
     $                     WORK( N-1+M ), Z( 1, L-1 ), LDZ )
            ELSE
               CALL SLAE2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2 )
            END IF
            D( L-1 ) = RT1
            D( L ) = RT2
            E( L-1 ) = ZERO
            L = L - 2
            IF( L.GE.LEND )
     $         GO TO 90
            GO TO 10
         END IF
*
         IF( JTOT.EQ.NMAXIT )
     $      GO TO 140
         JTOT = JTOT + 1
*
*        Form shift.
*
         G = ( D( L-1 )-P ) / ( TWO*E( L-1 ) )
         R = SLAPY2( G, ONE )
         G = D( M ) - P + ( E( L-1 ) / ( G+SIGN( R, G ) ) )
*
         S = ONE
         C = ONE
         P = ZERO
*
*        Inner loop
*
         LM1 = L - 1
         DO 120 I = M, LM1
            F = S*E( I )
            B = C*E( I )
            CALL SLARTG( G, F, C, S, R )
            IF( I.NE.M )
     $         E( I-1 ) = R
            G = D( I ) - P
            R = ( D( I+1 )-G )*S + TWO*C*B
            P = S*R
            D( I ) = G + P
            G = C*R - B
*
*           If eigenvectors are desired, then save rotations.
*
            IF( ICOMPZ.GT.0 ) THEN
               WORK( I ) = C
               WORK( N-1+I ) = S
            END IF
*
  120    CONTINUE
*
*        If eigenvectors are desired, then apply saved rotations.
*
         IF( ICOMPZ.GT.0 ) THEN
            MM = L - M + 1
            CALL SLASR( 'R', 'V', 'F', N, MM, WORK( M ), WORK( N-1+M ),
     $                  Z( 1, M ), LDZ )
         END IF
*
         D( L ) = D( L ) - P
         E( LM1 ) = G
         GO TO 90
*
*        Eigenvalue found.
*
  130    CONTINUE
         D( L ) = P
*
         L = L - 1
         IF( L.GE.LEND )
     $      GO TO 90
         GO TO 10
*
      END IF
*
*     Set error -- no convergence to an eigenvalue after a total
*     of N*MAXIT iterations.
*
  140 CONTINUE
      DO 150 I = 1, N - 1
         IF( E( I ).NE.ZERO )
     $      INFO = INFO + 1
  150 CONTINUE
      RETURN
*
*     Order eigenvalues and eigenvectors.
*
  160 CONTINUE
      DO 180 II = 2, N
         I = II - 1
         K = I
         P = D( I )
         DO 170 J = II, N
            IF( D( J ).LT.P ) THEN
               K = J
               P = D( J )
            END IF
  170    CONTINUE
         IF( K.NE.I ) THEN
            D( K ) = D( I )
            D( I ) = P
            IF( ICOMPZ.GT.0 )
     $         CALL SSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 )
         END IF
  180 CONTINUE
*
      RETURN
*
*     End of SSTEQR
*
      END

max_stars_repo_name: khinsen/HYDROLIB

max_issues_repo_name: khinsen/HYDROLIB

max_forks_repo_name: khinsen/HYDROLIB


text:
! unix_errno.f90
module unix_errno
    use, intrinsic :: iso_c_binding
    implicit none

#if defined (__linux__)

    integer(kind=c_int), parameter, public :: EPERM           = 1
    integer(kind=c_int), parameter, public :: ENOENT          = 2
    integer(kind=c_int), parameter, public :: ESRCH           = 3
    integer(kind=c_int), parameter, public :: EINTR           = 4
    integer(kind=c_int), parameter, public :: EIO             = 5
    integer(kind=c_int), parameter, public :: ENXIO           = 6
    integer(kind=c_int), parameter, public :: E2BIG           = 7
    integer(kind=c_int), parameter, public :: ENOEXEC         = 8
    integer(kind=c_int), parameter, public :: EBADF           = 9
    integer(kind=c_int), parameter, public :: ECHILD          = 10
    integer(kind=c_int), parameter, public :: EAGAIN          = 11
    integer(kind=c_int), parameter, public :: ENOMEM          = 12
    integer(kind=c_int), parameter, public :: EACCES          = 13
    integer(kind=c_int), parameter, public :: EFAULT          = 14
    integer(kind=c_int), parameter, public :: ENOTBLK         = 15
    integer(kind=c_int), parameter, public :: EBUSY           = 16
    integer(kind=c_int), parameter, public :: EEXIST          = 17
    integer(kind=c_int), parameter, public :: EXDEV           = 18
    integer(kind=c_int), parameter, public :: ENODEV          = 19
    integer(kind=c_int), parameter, public :: ENOTDIR         = 20
    integer(kind=c_int), parameter, public :: EISDIR          = 21
    integer(kind=c_int), parameter, public :: EINVAL          = 22
    integer(kind=c_int), parameter, public :: ENFILE          = 23
    integer(kind=c_int), parameter, public :: EMFILE          = 24
    integer(kind=c_int), parameter, public :: ENOTTY          = 25
    integer(kind=c_int), parameter, public :: ETXTBSY         = 26
    integer(kind=c_int), parameter, public :: EFBIG           = 27
    integer(kind=c_int), parameter, public :: ENOSPC          = 28
    integer(kind=c_int), parameter, public :: ESPIPE          = 29
    integer(kind=c_int), parameter, public :: EROFS           = 30
    integer(kind=c_int), parameter, public :: EMLINK          = 31
    integer(kind=c_int), parameter, public :: EPIPE           = 32
    integer(kind=c_int), parameter, public :: EDOM            = 33
    integer(kind=c_int), parameter, public :: ERANGE          = 34
    integer(kind=c_int), parameter, public :: EDEADLK         = 35
    integer(kind=c_int), parameter, public :: ENAMETOOLONG    = 36
    integer(kind=c_int), parameter, public :: ENOLCK          = 37
    integer(kind=c_int), parameter, public :: ENOSYS          = 38
    integer(kind=c_int), parameter, public :: ENOTEMPTY       = 39
    integer(kind=c_int), parameter, public :: ELOOP           = 40
    integer(kind=c_int), parameter, public :: EWOULDBLOCK     = EAGAIN
    integer(kind=c_int), parameter, public :: ENOMSG          = 42
    integer(kind=c_int), parameter, public :: EIDRM           = 43
    integer(kind=c_int), parameter, public :: ECHRNG          = 44
    integer(kind=c_int), parameter, public :: EL2NSYNC        = 45
    integer(kind=c_int), parameter, public :: EL3HLT          = 46
    integer(kind=c_int), parameter, public :: EL3RST          = 47
    integer(kind=c_int), parameter, public :: ELNRNG          = 48
    integer(kind=c_int), parameter, public :: EUNATCH         = 49
    integer(kind=c_int), parameter, public :: ENOCSI          = 50
    integer(kind=c_int), parameter, public :: EL2HLT          = 51
    integer(kind=c_int), parameter, public :: EBADE           = 52
    integer(kind=c_int), parameter, public :: EBADR           = 53
    integer(kind=c_int), parameter, public :: EXFULL          = 54
    integer(kind=c_int), parameter, public :: ENOANO          = 55
    integer(kind=c_int), parameter, public :: EBADRQC         = 56
    integer(kind=c_int), parameter, public :: EBADSLT         = 57
    integer(kind=c_int), parameter, public :: EDEADLOCK       = EDEADLK
    integer(kind=c_int), parameter, public :: EBFONT          = 59
    integer(kind=c_int), parameter, public :: ENOSTR          = 60
    integer(kind=c_int), parameter, public :: ENODATA         = 61
    integer(kind=c_int), parameter, public :: ETIME           = 62
    integer(kind=c_int), parameter, public :: ENOSR           = 63
    integer(kind=c_int), parameter, public :: ENONET          = 64
    integer(kind=c_int), parameter, public :: ENOPKG          = 65
    integer(kind=c_int), parameter, public :: EREMOTE         = 66
    integer(kind=c_int), parameter, public :: ENOLINK         = 67
    integer(kind=c_int), parameter, public :: EADV            = 68
    integer(kind=c_int), parameter, public :: ESRMNT          = 69
    integer(kind=c_int), parameter, public :: ECOMM           = 70
    integer(kind=c_int), parameter, public :: EPROTO          = 71
    integer(kind=c_int), parameter, public :: EMULTIHOP       = 72
    integer(kind=c_int), parameter, public :: EDOTDOT         = 73
    integer(kind=c_int), parameter, public :: EBADMSG         = 74
    integer(kind=c_int), parameter, public :: EOVERFLOW       = 75
    integer(kind=c_int), parameter, public :: ENOTUNIQ        = 76
    integer(kind=c_int), parameter, public :: EBADFD          = 77
    integer(kind=c_int), parameter, public :: EREMCHG         = 78
    integer(kind=c_int), parameter, public :: ELIBACC         = 79
    integer(kind=c_int), parameter, public :: ELIBBAD         = 80
    integer(kind=c_int), parameter, public :: ELIBSCN         = 81
    integer(kind=c_int), parameter, public :: ELIBMAX         = 82
    integer(kind=c_int), parameter, public :: ELIBEXEC        = 83
    integer(kind=c_int), parameter, public :: EILSEQ          = 84
    integer(kind=c_int), parameter, public :: ERESTART        = 85
    integer(kind=c_int), parameter, public :: ESTRPIPE        = 86
    integer(kind=c_int), parameter, public :: EUSERS          = 87
    integer(kind=c_int), parameter, public :: ENOTSOCK        = 88
    integer(kind=c_int), parameter, public :: EDESTADDRREQ    = 89
    integer(kind=c_int), parameter, public :: EMSGSIZE        = 90
    integer(kind=c_int), parameter, public :: EPROTOTYPE      = 91
    integer(kind=c_int), parameter, public :: ENOPROTOOPT     = 92
    integer(kind=c_int), parameter, public :: EPROTONOSUPPORT = 93
    integer(kind=c_int), parameter, public :: ESOCKTNOSUPPORT = 94
    integer(kind=c_int), parameter, public :: EOPNOTSUPP      = 95
    integer(kind=c_int), parameter, public :: EPFNOSUPPORT    = 96
    integer(kind=c_int), parameter, public :: EAFNOSUPPORT    = 97
    integer(kind=c_int), parameter, public :: EADDRINUSE      = 98
    integer(kind=c_int), parameter, public :: EADDRNOTAVAIL   = 99
    integer(kind=c_int), parameter, public :: ENETDOWN        = 100
    integer(kind=c_int), parameter, public :: ENETUNREACH     = 101
    integer(kind=c_int), parameter, public :: ENETRESET       = 102
    integer(kind=c_int), parameter, public :: ECONNABORTED    = 103
    integer(kind=c_int), parameter, public :: ECONNRESET      = 104
    integer(kind=c_int), parameter, public :: ENOBUFS         = 105
    integer(kind=c_int), parameter, public :: EISCONN         = 106
    integer(kind=c_int), parameter, public :: ENOTCONN        = 107
    integer(kind=c_int), parameter, public :: ESHUTDOWN       = 108
    integer(kind=c_int), parameter, public :: ETOOMANYREFS    = 109
    integer(kind=c_int), parameter, public :: ETIMEDOUT       = 110
    integer(kind=c_int), parameter, public :: ECONNREFUSED    = 111
    integer(kind=c_int), parameter, public :: EHOSTDOWN       = 112
    integer(kind=c_int), parameter, public :: EHOSTUNREACH    = 113
    integer(kind=c_int), parameter, public :: EALREADY        = 114
    integer(kind=c_int), parameter, public :: EINPROGRESS     = 115
    integer(kind=c_int), parameter, public :: ESTALE          = 116
    integer(kind=c_int), parameter, public :: EUCLEAN         = 117
    integer(kind=c_int), parameter, public :: ENOTNAM         = 118
    integer(kind=c_int), parameter, public :: ENAVAIL         = 119
    integer(kind=c_int), parameter, public :: EISNAM          = 120
    integer(kind=c_int), parameter, public :: EREMOTEIO       = 121
    integer(kind=c_int), parameter, public :: EDQUOT          = 122
    integer(kind=c_int), parameter, public :: ENOMEDIUM       = 123
    integer(kind=c_int), parameter, public :: EMEDIUMTYPE     = 124
    integer(kind=c_int), parameter, public :: ECANCELED       = 125
    integer(kind=c_int), parameter, public :: ENOKEY          = 126
    integer(kind=c_int), parameter, public :: EKEYEXPIRED     = 127
    integer(kind=c_int), parameter, public :: EKEYREVOKED     = 128
    integer(kind=c_int), parameter, public :: EKEYREJECTED    = 129
    integer(kind=c_int), parameter, public :: EOWNERDEAD      = 130
    integer(kind=c_int), parameter, public :: ENOTRECOVERABLE = 131
    integer(kind=c_int), parameter, public :: ERFKILL         = 132
    integer(kind=c_int), parameter, public :: EHWPOISON       = 133

#elif defined (__FreeBSD__)

    integer(kind=c_int), parameter, public :: EPERM           = 1
    integer(kind=c_int), parameter, public :: ENOENT          = 2
    integer(kind=c_int), parameter, public :: ESRCH           = 3
    integer(kind=c_int), parameter, public :: EINTR           = 4
    integer(kind=c_int), parameter, public :: EIO             = 5
    integer(kind=c_int), parameter, public :: ENXIO           = 6
    integer(kind=c_int), parameter, public :: E2BIG           = 7
    integer(kind=c_int), parameter, public :: ENOEXEC         = 8
    integer(kind=c_int), parameter, public :: EBADF           = 9
    integer(kind=c_int), parameter, public :: ECHILD          = 10
    integer(kind=c_int), parameter, public :: EDEADLK         = 11
    integer(kind=c_int), parameter, public :: ENOMEM          = 12
    integer(kind=c_int), parameter, public :: EACCES          = 13
    integer(kind=c_int), parameter, public :: EFAULT          = 14
    integer(kind=c_int), parameter, public :: ENOTBLK         = 15
    integer(kind=c_int), parameter, public :: EBUSY           = 16
    integer(kind=c_int), parameter, public :: EEXIST          = 17
    integer(kind=c_int), parameter, public :: EXDEV           = 18
    integer(kind=c_int), parameter, public :: ENODEV          = 19
    integer(kind=c_int), parameter, public :: ENOTDIR         = 20
    integer(kind=c_int), parameter, public :: EISDIR          = 21
    integer(kind=c_int), parameter, public :: EINVAL          = 22
    integer(kind=c_int), parameter, public :: ENFILE          = 23
    integer(kind=c_int), parameter, public :: EMFILE          = 24
    integer(kind=c_int), parameter, public :: ENOTTY          = 25
    integer(kind=c_int), parameter, public :: ETXTBSY         = 26
    integer(kind=c_int), parameter, public :: EFBIG           = 27
    integer(kind=c_int), parameter, public :: ENOSPC          = 28
    integer(kind=c_int), parameter, public :: ESPIPE          = 29
    integer(kind=c_int), parameter, public :: EROFS           = 30
    integer(kind=c_int), parameter, public :: EMLINK          = 31
    integer(kind=c_int), parameter, public :: EPIPE           = 32
    integer(kind=c_int), parameter, public :: EDOM            = 33
    integer(kind=c_int), parameter, public :: ERANGE          = 34
    integer(kind=c_int), parameter, public :: EAGAIN          = 35
    integer(kind=c_int), parameter, public :: EWOULDBLOCK     = EAGAIN
    integer(kind=c_int), parameter, public :: EINPROGRESS     = 36
    integer(kind=c_int), parameter, public :: EALREADY        = 37
    integer(kind=c_int), parameter, public :: ENOTSOCK        = 38
    integer(kind=c_int), parameter, public :: EDESTADDRREQ    = 39
    integer(kind=c_int), parameter, public :: EMSGSIZE        = 40
    integer(kind=c_int), parameter, public :: EPROTOTYPE      = 41
    integer(kind=c_int), parameter, public :: ENOPROTOOPT     = 42
    integer(kind=c_int), parameter, public :: EPROTONOSUPPORT = 43
    integer(kind=c_int), parameter, public :: ESOCKTNOSUPPORT = 44
    integer(kind=c_int), parameter, public :: EOPNOTSUPP      = 45
    integer(kind=c_int), parameter, public :: ENOTSUP         = EOPNOTSUPP
    integer(kind=c_int), parameter, public :: EPFNOSUPPORT    = 46
    integer(kind=c_int), parameter, public :: EAFNOSUPPORT    = 47
    integer(kind=c_int), parameter, public :: EADDRINUSE      = 48
    integer(kind=c_int), parameter, public :: EADDRNOTAVAIL   = 49
    integer(kind=c_int), parameter, public :: ENETDOWN        = 50
    integer(kind=c_int), parameter, public :: ENETUNREACH     = 51
    integer(kind=c_int), parameter, public :: ENETRESET       = 52
    integer(kind=c_int), parameter, public :: ECONNABORTED    = 53
    integer(kind=c_int), parameter, public :: ECONNRESET      = 54
    integer(kind=c_int), parameter, public :: ENOBUFS         = 55
    integer(kind=c_int), parameter, public :: EISCONN         = 56
    integer(kind=c_int), parameter, public :: ENOTCONN        = 57
    integer(kind=c_int), parameter, public :: ESHUTDOWN       = 58
    integer(kind=c_int), parameter, public :: ETOOMANYREFS    = 59
    integer(kind=c_int), parameter, public :: ETIMEDOUT       = 60
    integer(kind=c_int), parameter, public :: ECONNREFUSED    = 61
    integer(kind=c_int), parameter, public :: ELOOP           = 62
    integer(kind=c_int), parameter, public :: ENAMETOOLONG    = 63
    integer(kind=c_int), parameter, public :: EHOSTDOWN       = 64
    integer(kind=c_int), parameter, public :: EHOSTUNREACH    = 65
    integer(kind=c_int), parameter, public :: ENOTEMPTY       = 66
    integer(kind=c_int), parameter, public :: EPROCLIM        = 67
    integer(kind=c_int), parameter, public :: EUSERS          = 68
    integer(kind=c_int), parameter, public :: EDQUOT          = 69

#endif

    public :: c_errno

    ! Interface to `error_number()` in `errno.c`
    interface
        ! int error_number()
        function c_errno() bind(c, name='error_number')
            import :: c_int
            integer(c_int) :: c_errno
        end function c_errno
    end interface
end module unix_errno

max_stars_repo_name: interkosmos/fortran-unix

max_issues_repo_name: interkosmos/fortran-unix

max_forks_repo_name: interkosmos/fortran-unix


text: 20070905 10:02:32 nbsp Howdy, and welcome to the wiki!  While were all pretty opinionated around here, please try to keep it focused on the subject and avoid personal attacks on individuals.  Remember, the people on DavisWiki are real people... and your neighbors. Users/JabberWokky

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Copyright 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011
c Washington University, Mallinckrodt Institute of Radiology.
c All rights reserved.
c This software may not be reproduced, copied, or distributed without written
c permission of washington university. For further information contact A. Z. Snyder.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c $Header: /data/petsun4/data1/src_solaris/librms/RCS/dmatinv.f,v 1.6 2011/01/19 04:54:43 avi Exp $
c $Log: dmatinv.f,v $
c Revision 1.6  2011/01/19  04:54:43  avi
c restore malloc functionality
c
c Revision 1.5  2010/02/12  07:14:43  avi
c hard-coded matrix dimenssion limit 512->2048
c
c Revision 1.4  2007/05/05  20:17:03  avi
c typo
c
c Revision 1.3  2007/05/05  20:14:08  avi
c Revision 1.2  2007/05/05  20:05:19  avi
c gcc v3 compliant (no pointers, drand()->rand)
c
c Revision 1.1  2005/09/04  06:24:51  avi
c Initial revision
c
c     program dmatinvtst
      subroutine dmatinvtst
      implicit real*8 (a-h,o-z)
      real*4 rand
      parameter (n=384)
      real*8 q(n,n),t(n,n),r(n,n)
      do 1 i=1,n
      do 1 j=1,n
    1 t(i,j)=0.5*(rand(0)-.1)
      do 2 i=1,n
      do 2 j=1,n
      q(i,j)=0.0
      do 2 k=1,n
      q(i,j)=q(i,j)+t(i,k)*t(j,k)
    2 r(i,j)=q(i,j)
      call dmatinv(q,n,det)
      write(*,"('dimension=',i4,' det=',e12.6)")n,det
  101 format(3(3f16.12,4x))
      do 12 i=1,n
      do 12 j=1,n
      t(i,j)=0.0
      do 12 k=1,n
   12 t(i,j)=t(i,j)+q(i,k)*r(k,j)
      write(*,"()")
      if (n.eq.3)then
        do i=1,n
          write(*,101)(r(i,k),k=1,n),(q(i,k),k=1,n),(t(i,k),k=1,n)
         enddo
      endif
      derr=0.0
      ierr=0
      jerr=0
      do 6 i=1,n
      do 6 j=1,n
      x=t(i,j)
      if(i.eq.j)x=x-1.0
      if(dabs(x).gt.derr)then
        ierr=i
        jerr=j
        derr=dabs(x)
      endif
    6 continue
      write(*,"('maximum error = ',e12.6)")derr
      end
      subroutine dmatinv(a,n,det)
c     algorithm follows Bevington, Phillip R.,
c     Data Reduction and Error Analysis for the Physical Sciences,
c     Mcgraw-Hill, New York, 1969.
      implicit real*8 (a-h,o-z)
      real*8 a(n,n),det
      pointer (pik,ik),(pjk,jk)
      integer*4 ik(n),jk(n)
c     integer*4 ik(2048),jk(2048)
c     if(n.gt.2048)then
c       write (*,"('dmatinv matrix dimension limit (',i4,') exeeded')")2048
c       call exit (-1)
c     endif
      pik=malloc(4*n)
      pjk=malloc(4*n)
      det=1.
      do 100 k=1,n
      vmax=0.
   21 do 30 i=k,n
      do 30 j=k,n
      if(dabs(vmax).gt.dabs(a(i,j)))goto 30
      vmax=a(i,j)
      ik(k)=i
      jk(k)=j
   30 continue
      i=ik(k)
c     if(i-k)21,51,43
      if(i.eq.k)goto 51
      if(i.lt.k)goto 21
   43 do 50 j=1,n
      t=a(k,j)
      a(k,j)=a(i,j)
   50 a(i,j)=-t
   51 j=jk(k)
c     if(j-k)21,61,53
      if(j.eq.k)goto 61
      if(j.lt.k)goto 21
   53 do 60 i=1,n
      t=a(i,k)
      a(i,k)=a(i,j)
   60 a(i,j)=-t
   61 do 70 i=1,n
      if(i.ne.k)a(i,k)=-a(i,k)/vmax
   70 continue
   71 do 80 i=1,n
      do 80 j=1,n
      if(i.ne.k.and.j.ne.k)a(i,j)=a(i,j)+a(i,k)*a(k,j)
   80 continue
   81 do 90 j=1,n
      if(j.ne.k)a(k,j)=a(k,j)/vmax
   90 continue
      a(k,k)=1./vmax
  100 det=det*vmax
      do 130 l=1,n
      k=n-l+1
      j=ik(k)
c     if(j-k)111,111,105
      if(j.le.k)goto 111
  105 do 110 i=1,n
      t=a(i,k)
      a(i,k)=-a(i,j)
  110 a(i,j)=t
  111 i=jk(k)
c     if(i-k)130,130,113
      if(i.le.k)goto 130
  113 do 120 j=1,n
      t=a(k,j)
      a(k,j)=-a(i,j)
  120 a(i,j)=t
  130 continue
      call free(pik)
      call free(pjk)
      return
      end

max_stars_repo_name: llevitis/PUP

max_issues_repo_name: llevitis/PUP

max_forks_repo_name: llevitis/PUP