text:
 THE WHAT
The Charity Classic is a flag football tournament for philanthropyminded UC Davis students.

 THE WHEN AND WHERE

  The Second Annual Charity Classic
  Sunday, April 30th 2006
  Russell Athletic Fields (Next to Toomey Field)
  Mens and coed leagues > only UC Davis students may register. (Coed teams must have at least 3 female players.)
  The cost is $40 per team of 710 individuals. (Registration fees due the day of the tournament.)

 HOW TO REGISTER
  Obtain a http://www.alphagammaomega.net/cc/CC%2006%20%20Team%20Roster%20Form.pdf Team Roster form, and turn it in to our information table on the MU Quad from 121:30, April 10th  21st

  Mail your team roster to CHARITY CLASSIC, 310 Parkway Circle, Davis, CA 95616.

  Register through email by sending your team roster to MailTo(charityclassic@gmail.com)

  IMPORTANT: Your team roster must include your team members names, ID numbers (to be sure all team members are students), and a contact email address.

 THE WHY

The purpose of our program is to encourage the student population of UC Davis to look beyond themselves by learning what it means to sacrifice of their own resources to benefit the global community at large. Our goal is for the student population of UCD to gain a better understanding of service, compassion, and selflessness by giving them the opportunity to give a little of themselves to meet an urgent need facing countless people around the world.

In the world in which we live today, nearly 1 billion people live without access to clean water  a terrible reality which contributes to death, disease, and poverty all over the globe. In response to this need, The 1000 Wells organization is a group that drills wells for impoverished communities in Africa, bringing hope and abundant life to people lacking this essential resource. It is this organization that we have decided to partner with this year through our philanthropy event, The 2nd Annual Charity Classic. This event consists of a flag football tournament, in which we hope to increase the awareness within the student body of great needs of people all over the world while having a great time in healthy competition and great fun. It is being put on by Alpha Gamma Omega, a ChristCentered Fraternity, and Alpha Delta Chi, a Christian sorority. In regards to the tournament, each team will pay a registration fee of $40, all of which will go directly to the 1000 Wells Project Davis 1000 Wells Project. The tournament will take place all day on April 30th 2006, and our desire is to raise enough money to build one well. How amazing it would be, how fulfilling and enriching to be able to see, touch, and experience what can happen when a body of people decides to look beyond their own needs and begins to care for the needs of others.

 DOWNLOADS
  http://www.alphagammaomega.net/cc/CC%2006%20%20Team%20Roster%20Form.pdf Team Roster Form
  http://www.alphagammaomega.net/cc/CC%2006%20%20Full%20Page%20Flyer%20.pdf Full Page Flyer
  http://www.alphagammaomega.net/cc/CC%2006%20%20Powerpoint%20Slide.pdf Powerpoint Slide
  http://www.alphagammaomega.net/cc/CC%2006%20%20Quarter%20Sheet%20Flyer.pdf Quarter Sheet Flyer

 ON THE WEB
  http://www.1000wells.org 1000 Wells
  http://iota.ago.org/ Alpha Gamma Omega
  http://www.alphadeltachi.net/ Alpha Delta Chi

 IN THE NEWS
  http://www.sfgate.com/cgibin/article.cgi?f/c/a/2006/04/04/MNGQGI34MS1.DTL&feedrss.news Water Scarce in Kenya (4.4.06 SF Times)
  http://www.sfgate.com/cgibin/article.cgi?f/c/a/2006/04/04/MNG6II33N51.DTL&feedrss.news Cost of Water in Africa Soaring (4.4.06 SF Times)

20060406 01:34:43 nbsp A water shortage in Africa? What about in California? Water shortages in Africa are just as much a result of poor civil planning in Africa as ours are here. In 1 Timothy 5:8, King James Version, it says But if any provide not for his own, and specially for those of his own house, he hath denied the faith, and is worse than an infidel. That stadium picture doesnt look like anywhere in Davis either. Users/NickSchmalenberger

20060408 00:12:28 nbsp Does water come out of your tap when you turn it on?  Until the time that you have to walk 90 miles for water, only to let your camel drink first,  I believe that we should first consider those who are truly threatened by water shortages and famines.  What happens when we have a water shortage in California, the price goes up what, a dollar or two?  In africa people die.  Not by the tens or hundreds, but by the thousands.  For one as blessed as yourself to even compare the problems is very disheartening.  And its a flag football tournament, thats a football field in the picture... Users/LeviRyan
 That is so patronizing and condescending. Do you think all Africans have camels? What is the farthest you have ever walked? Maybe, just maybe, most of them drive cars, ride bicycles, or take buses and trains like people in California. They dont need your pity. There are places in California where, if people choose to live there, they would have to walk much farther than 90 miles for water. People living in similar places in Africa also make that choice, or more likely just dont live there, like in California. Most cities anywhere in the world are built around existing sources of water like rivers, lakes, or oceans, and if they are contaminated or dry that is again a result of poor civil planning and overpopulation. Again, California has these issues too. Have you ever heard of Owens Lake or Mono Lake? Owens Lake doesnt exist anymore because Los Angeles drank it all and now they are starting in on Mono Lake and the entire Colorado River if the rest of the watershed will let them. I dont blame them either, I just blame the planners that oversold the land then had to bring in water from outside. As for the picture, it is not only a football field, but a large stadium, and so completely irrelevant to this event besides not being in Davis which already would make it irrelevant. Users/NickSchmalenberger

20060408 01:22:23 nbsp $40 or even 1000 wells isnt going to do very much good in Africa. What happens if a city grows too large and the whole water table dries up? A much better plan helping prevent that problem would be to invest that $40 in a business and grow the business until you are rich. Then donate money to universities for the establishment and improvement of civil planning departments, and not just in Africa, these problems exist in every part of the world. Users/NickSchmalenberger
   20060408 11:05:04 nbsp Are you kidding? Users/RyanPrendiville
     No. Are you? Users/NickSchmalenberger

20060420 13:24:06 nbsp Check this out: http://www.ug.edu.gh/deptdetails.php?recordID27 , the University of Ghana has a Department of Geography Resource and Development. Here: http://www.udsm.ac.tz/department_s/index.html , the University of Dar es Salaam has a Department of Urban and Rural Planning, a Department of Geography, a Department of Land Managment and Valuation, a Department of Environmental Engineering, and other departments that would be relevant to improvement of civil planning in Africa. Now tell me Africa doesnt know what civil planning is. Users/NickSchmalenberger

20060420 16:38:20 nbsp The 1000 Well Project does not take Civil Planning power away from the people of Africa.  It provides the capital necessary to implement the infrastructure supported by local governments.  Knowing where, how, and why to build a well is very different than being financially able to build that well.  The organizations funded by Blood: Water Mission act in partnership with the community.  BWM does what the local people think is best, forcing some other idea on them would be wrong.  Currently, local and national governments want wells.  Thats what BWM is going to give them. Users/JillWeinstein
  More water resource infrastructure is what Los Angeles built and look what happened. Civil planning doesnt just mean knowing how to extract resources, but also how to conserve them and what kind of population an environment can sustainably support. More wells dont make more water. Users/NickSchmalenberger
 The difference between NO water and SOME water is huge. Many of these communities are filled with refugees that have been forced onto marginalized lands with inadequate resources.  They have no choice when it comes to where they are living, and in matters of life and death I see no reason why we should allow innocent people to die.  The Los Angeles water system created by Mulholland was motivated by greed and a complete diregard for the Owens Valley, you have that correct.  However, Los Angeles grew as a result of the water projects, not the other way around.  In Africa you have population growth(partly b/c of individuals escaping genocide) in dessicated areas with water miles and miles away.  Also, in LA the water projects were primarily used to promote agriculture, a luxury.  BWM targets communities that lack DRINKING water, a BASIC NEED.  Im sorry, but if you think that completely sustainable water supplies are possible in communities that dont even have water to begin with, youre dreaming.  BWM does all that it can to promote ecological awareness and sustainability, but areas that lack basic infrastructure arent capable of maintaining sophisticated and costly systems that arent even considered costeffective in the first world. Users/JillWeinstein
  Im just saying building wells will only create worse problems later. The fact that there is a shortage now means that there are already too many people there. Building wells will encourage people to live there permanently, have families, and then there will be even more demand for water. As you say, these people are living on marginalized lands with inadequate resources and this is because the whole world is overpopulated. If nothing is done about that, people will eventually realize how expensive it is to have children, have fewer of them, and the problem will solve itself. But meanwhile, it will be expensive to have children, in the form of famines, wars, genocides, and pandemics. The only things that can be done to reduce the population faster without hurting people (as with genocide) will still take a long time to implement, so I think that some suffering is inevitable. It is really the same as any other economic shortage, except much worse. Users/NickSchmalenberger
  In response to Jill, with LA, it is a cycle, of LA needs water, LA gets more water, LA grows even larger, and then LA needs more water. While in the past the LA water projects may have been primarily for agriculture, but since that initial period, LAs desire for water has been driven by the aforementioned cycle, not by farmers. As a response to Nick, I understand your economic argument. I do agree with your assertion that if we increase water production, water demand will also increase in the form of increased population. However, with increased water production there will be increased agricultural production, which will be increased economic activity for the region. Additionally, at this time the affected populations and their governments are not capable of thinking in the long term that you suggest, because they are still attempting to cope on a day to day basis. In this manner, the increased economic activity and giving the populace less to worry about in the short term will allow them to try to address the long term problem. Let me reiterate that I agree with your ideas and logic to some extent economically, but I feel that this as a stopgap measure will allow the government to focus on the larger issues you brought up. This is why the world needs economists. Users/MaxMikalonis
   There will not be increased water production, only increased water use. Relative to the watershed, wells use water, even though they produce it relative to the well users. Even if some sort of facility for chemical synthesis or purification of otherwise undrinkable water were built, such as a desalination plant, it would require energy. wiki:WikiPedia:TANSTAAFL There is no such thing as a free lunch. As I said above, more wells dont make more water. Users/NickSchmalenberger


text:
C SUBROUTINE TO GET PHYSICAL MEMORY SIZE
        SUBROUTINE GET_MEM_SIZE(MEMSIZE)
        INTEGER MEMSIZE,WSEXT,WSLIM,WSA ! WSA IS ADJUSTMENT FOR BUFFERS ETC

	CALL XVPARM('WSA',WSA,I,K,0) ! GET WORKING SET ADJUSTMENT FOR BUFFERS

	WSEXT=0 ! INITIALIZE WSEXT
5	CALL SYS$ADJWSL(3000,WSLIM) ! TRY TO INCREASE WORKING SET
	IF(WSLIM.EQ.WSEXT) GOTO 10 ! IF NO CHANGE, GOTO 10
	WSEXT=WSLIM
	GOTO 5
10	MEMSIZE = WSEXT*512-WSA    !ALLOW WSA BYTES FOR BUFFERS ETC.
        RETURN
        END

C SUBROUTINE TO EXPLAIN WHY LGEOM CAN'T HANDLE IMAGE.
        SUBROUTINE TOOBIG
	CALL MABEND('** PICTURE IS TOO LARGE FOR THIS WORKING SET')
        RETURN
        END

C SUBROUTINE TO CHECK FOR SMALL PAGE FILE QUOTA UNDER VMS
        SUBROUTINE CHECK_VMS_PAGE_FILE_QUOTA
	INTEGER*4 RPFQ
	CALL GRPFQ(RPFQ) ! GET REMAINING PAGE FILE QUOTA
	IF(RPFQ.EQ.0)
     .     CALL XVMESSAGE('PAGE FILE QUOTA MAY BE TOO SMALL',' ')
        RETURN
        END

C	GWSEXT - GET WORKING SET EXTENT (MAX PAGES OF MEMORY)
C	GRPFQ  - GET REMAINING PAGE FILE QUOTA

	SUBROUTINE GWSEXT(WSEXT)
	INTEGER*4 WSEXT,RPFQ
	INCLUDE '($JPIDEF)'
	CALL LIB$GETJPI(JPI$_WSEXTENT,,,WSEXT)
	RETURN

	ENTRY GRPFQ(RPFQ)
	CALL LIB$GETJPI(JPI$_PAGFILCNT,,,RPFQ)
	RETURN

	END

ext: f

max_stars_repo_name: NASA-AMMOS/VICAR

max_issues_repo_name: NASA-AMMOS/VICAR

max_forks_repo_name: NASA-AMMOS/VICAR


text:
C Copyright(C) 1999-2020 National Technology & Engineering Solutions
C of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C See packages/seacas/LICENSE for details

C========================================================================
*     DECK, SETON0
      SUBROUTINE SETON0(ICONA,NELTN,SOLEA,SOLENA,IDBLK,XA,YA,ZA,ISTP,
     &     ITT,iblk)
C
C     *********************************************************************
C
C     Subroutine SETON0 extracts nodal values of shell element variables by
C     looping over each element and summing the value of the variable
C     in that element to each node in the connectivity list for that
C     element. Then the nodal summation of element variables is divided
C     by the number of elements that contributed to that node (resulting
C     in a nodal average of the element value.) This is done for the old
C     mesh elements and nodes to facilitate interpolation.
C
C     Each element block must be processed independently in order to
C     avoid averaging element variables across material boundaries.
C     Note: the last set of DO loops acts over all nodes; to make sense
C     one element block must be completely processed before another
C     element block is sent into this subroutine.
C
C     Calls subroutine VOL
C
c     Called by MAPVAR
C
C     *********************************************************************
C
C     ICONA       mesh-A connectivity (1:nelnda,1:numeba)
C     NELTN       number of elements tied to each node (1:nodesa)
C     SOLEA       element variables (1:numeba,1:nvarel)
C     SOLENA      element variables at nodes (1:nodesa,1:nvarel)
C     IDBLK       current element block I.D.
C     XA,YA,ZA    coordinates
C     XX,YY,ZZ    vector of coordinates of nodes for an element
C     ISTP        current time step
C     ITT         truth table
C     iblk        element block being processed (not ID)
C
C     *********************************************************************
C
      include 'aexds1.blk'
      include 'aexds2.blk'
      include 'amesh.blk'
      include 'ebbyeb.blk'
      include 'ex2tp.blk'
      include 'tapes.blk'
C
      DIMENSION ICONA(NELNDA,*), NELTN(*)
      DIMENSION SOLEA(NUMEBA,*), SOLENA(NODESA,NVAREL), ITT(NVAREL,*)
      DIMENSION XA(*), YA(*), ZA(*), XX(27), YY(27), ZZ(27)
C
C     *********************************************************************
C
      NNODES = 4
      DO I = 1, NODESA
         NELTN(I) = 0
         DO J = 1, NVAREL
            SOLENA(I,J) = 0.
         end do
      end do
C
      DO NEL = 1, NUMEBA
         DO I = 1, NNODES
C
C     number of elements associated with each node - used for
C     computing an average later on
C
            NELTN(ICONA(I,NEL)) = NELTN(ICONA(I,NEL)) + 1
         end do
      end do
C
      DO IVAR = 1, NVAREL
         IF (ITT(IVAR,iblk) .EQ. 0)GO TO 40
         CALL EXGEV(NTP2EX,ISTP,IVAR,IDBLK,NUMEBA,SOLEA(1,IVAR),IERR)
C
         IF (NAMVAR(nvargp+IVAR)(1:6) .EQ. 'ELMASS') THEN
C
C     replace element mass with nodal density for interpolation
C
            DO IEL = 1, NUMEBA
               DO I = 1, NNODES
                  XX(I) = XA(ICONA(I,IEL))
                  YY(I) = YA(ICONA(I,IEL))
                  ZZ(I) = ZA(ICONA(I,IEL))
               end do
               CALL VOL(ITYPE,XX,YY,ZZ,VOLUME)
               SOLEA(IEL,IVAR) = SOLEA(IEL,IVAR) / VOLUME
            end do
         END IF
C
C     accumulate element variables to nodes
C
         DO NEL = 1, NUMEBA
            DO I = 1, NNODES
               SOLENA(ICONA(I,NEL),IVAR) =
     &              SOLENA(ICONA(I,NEL),IVAR) + SOLEA(NEL,IVAR)
            end do
         end do
C
C     divide by number of elements contributing to each node (average)
C
         DO I = 1, NODESA
            IF(NELTN(I) .NE. 0)THEN
               SOLENA(I,IVAR) = SOLENA(I,IVAR) / dble(NELTN(I))
            END IF
         end do
 40      CONTINUE
      end do
      RETURN
      END

max_stars_repo_name: tokusanya/seacas

max_issues_repo_name: tokusanya/seacas

max_forks_repo_name: tokusanya/seacas


text:
*> \brief  DGTSVX computes the solution to system of linear equations A * X = B for GT matrices
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DGTSVX + dependencies
*>
*> [TGZ]
*>
*> [ZIP]
*>
*> [TXT]
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE DGTSVX( FACT, TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF,
*                          DU2, IPIV, B, LDB, X, LDX, RCOND, FERR, BERR,
*                          WORK, IWORK, INFO )
*
*       .. Scalar Arguments ..
*       CHARACTER          FACT, TRANS
*       INTEGER            INFO, LDB, LDX, N, NRHS
*       DOUBLE PRECISION   RCOND
*       ..
*       .. Array Arguments ..
*       INTEGER            IPIV( * ), IWORK( * )
*       DOUBLE PRECISION   B( LDB, * ), BERR( * ), D( * ), DF( * ),
*      $                   DL( * ), DLF( * ), DU( * ), DU2( * ), DUF( * ),
*      $                   FERR( * ), WORK( * ), X( LDX, * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DGTSVX uses the LU factorization to compute the solution to a real
*> system of linear equations A * X = B or A**T * X = B,
*> where A is a tridiagonal matrix of order N and X and B are N-by-NRHS
*> matrices.
*>
*> Error bounds on the solution and a condition estimate are also
*> provided.
*> \endverbatim
*
*> \par Description:
*  =================
*>
*> \verbatim
*>
*> The following steps are performed:
*>
*> 1. If FACT = 'N', the LU decomposition is used to factor the matrix A
*>    as A = L * U, where L is a product of permutation and unit lower
*>    bidiagonal matrices and U is upper triangular with nonzeros in
*>    only the main diagonal and first two superdiagonals.
*>
*> 2. If some U(i,i)=0, so that U is exactly singular, then the routine
*>    returns with INFO = i. Otherwise, the factored form of A is used
*>    to estimate the condition number of the matrix A.  If the
*>    reciprocal of the condition number is less than machine precision,
*>    INFO = N+1 is returned as a warning, but the routine still goes on
*>    to solve for X and compute error bounds as described below.
*>
*> 3. The system of equations is solved for X using the factored form
*>    of A.
*>
*> 4. Iterative refinement is applied to improve the computed solution
*>    matrix and calculate error bounds and backward error estimates
*>    for it.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] FACT
*> \verbatim
*>          FACT is CHARACTER*1
*>          Specifies whether or not the factored form of A has been
*>          supplied on entry.
*>          = 'F':  DLF, DF, DUF, DU2, and IPIV contain the factored
*>                  form of A; DL, D, DU, DLF, DF, DUF, DU2 and IPIV
*>                  will not be modified.
*>          = 'N':  The matrix will be copied to DLF, DF, and DUF
*>                  and factored.
*> \endverbatim
*>
*> \param[in] TRANS
*> \verbatim
*>          TRANS is CHARACTER*1
*>          Specifies the form of the system of equations:
*>          = 'N':  A * X = B     (No transpose)
*>          = 'T':  A**T * X = B  (Transpose)
*>          = 'C':  A**H * X = B  (Conjugate transpose = Transpose)
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrix A.  N >= 0.
*> \endverbatim
*>
*> \param[in] NRHS
*> \verbatim
*>          NRHS is INTEGER
*>          The number of right hand sides, i.e., the number of columns
*>          of the matrix B.  NRHS >= 0.
*> \endverbatim
*>
*> \param[in] DL
*> \verbatim
*>          DL is DOUBLE PRECISION array, dimension (N-1)
*>          The (n-1) subdiagonal elements of A.
*> \endverbatim
*>
*> \param[in] D
*> \verbatim
*>          D is DOUBLE PRECISION array, dimension (N)
*>          The n diagonal elements of A.
*> \endverbatim
*>
*> \param[in] DU
*> \verbatim
*>          DU is DOUBLE PRECISION array, dimension (N-1)
*>          The (n-1) superdiagonal elements of A.
*> \endverbatim
*>
*> \param[in,out] DLF
*> \verbatim
*>          DLF is DOUBLE PRECISION array, dimension (N-1)
*>          If FACT = 'F', then DLF is an input argument and on entry
*>          contains the (n-1) multipliers that define the matrix L from
*>          the LU factorization of A as computed by DGTTRF.
*>
*>          If FACT = 'N', then DLF is an output argument and on exit
*>          contains the (n-1) multipliers that define the matrix L from
*>          the LU factorization of A.
*> \endverbatim
*>
*> \param[in,out] DF
*> \verbatim
*>          DF is DOUBLE PRECISION array, dimension (N)
*>          If FACT = 'F', then DF is an input argument and on entry
*>          contains the n diagonal elements of the upper triangular
*>          matrix U from the LU factorization of A.
*>
*>          If FACT = 'N', then DF is an output argument and on exit
*>          contains the n diagonal elements of the upper triangular
*>          matrix U from the LU factorization of A.
*> \endverbatim
*>
*> \param[in,out] DUF
*> \verbatim
*>          DUF is DOUBLE PRECISION array, dimension (N-1)
*>          If FACT = 'F', then DUF is an input argument and on entry
*>          contains the (n-1) elements of the first superdiagonal of U.
*>
*>          If FACT = 'N', then DUF is an output argument and on exit
*>          contains the (n-1) elements of the first superdiagonal of U.
*> \endverbatim
*>
*> \param[in,out] DU2
*> \verbatim
*>          DU2 is DOUBLE PRECISION array, dimension (N-2)
*>          If FACT = 'F', then DU2 is an input argument and on entry
*>          contains the (n-2) elements of the second superdiagonal of
*>          U.
*>
*>          If FACT = 'N', then DU2 is an output argument and on exit
*>          contains the (n-2) elements of the second superdiagonal of
*>          U.
*> \endverbatim
*>
*> \param[in,out] IPIV
*> \verbatim
*>          IPIV is INTEGER array, dimension (N)
*>          If FACT = 'F', then IPIV is an input argument and on entry
*>          contains the pivot indices from the LU factorization of A as
*>          computed by DGTTRF.
*>
*>          If FACT = 'N', then IPIV is an output argument and on exit
*>          contains the pivot indices from the LU factorization of A;
*>          row i of the matrix was interchanged with row IPIV(i).
*>          IPIV(i) will always be either i or i+1; IPIV(i) = i indicates
*>          a row interchange was not required.
*> \endverbatim
*>
*> \param[in] B
*> \verbatim
*>          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
*>          The N-by-NRHS right hand side matrix B.
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*>          LDB is INTEGER
*>          The leading dimension of the array B.  LDB >= max(1,N).
*> \endverbatim
*>
*> \param[out] X
*> \verbatim
*>          X is DOUBLE PRECISION array, dimension (LDX,NRHS)
*>          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
*> \endverbatim
*>
*> \param[in] LDX
*> \verbatim
*>          LDX is INTEGER
*>          The leading dimension of the array X.  LDX >= max(1,N).
*> \endverbatim
*>
*> \param[out] RCOND
*> \verbatim
*>          RCOND is DOUBLE PRECISION
*>          The estimate of the reciprocal condition number of the matrix
*>          A.  If RCOND is less than the machine precision (in
*>          particular, if RCOND = 0), the matrix is singular to working
*>          precision.  This condition is indicated by a return code of
*>          INFO > 0.
*> \endverbatim
*>
*> \param[out] FERR
*> \verbatim
*>          FERR is DOUBLE PRECISION array, dimension (NRHS)
*>          The estimated forward error bound for each solution vector
*>          X(j) (the j-th column of the solution matrix X).
*>          If XTRUE is the true solution corresponding to X(j), FERR(j)
*>          is an estimated upper bound for the magnitude of the largest
*>          element in (X(j) - XTRUE) divided by the magnitude of the
*>          largest element in X(j).  The estimate is as reliable as
*>          the estimate for RCOND, and is almost always a slight
*>          overestimate of the true error.
*> \endverbatim
*>
*> \param[out] BERR
*> \verbatim
*>          BERR is DOUBLE PRECISION array, dimension (NRHS)
*>          The componentwise relative backward error of each solution
*>          vector X(j) (i.e., the smallest relative change in
*>          any element of A or B that makes X(j) an exact solution).
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is DOUBLE PRECISION array, dimension (3*N)
*> \endverbatim
*>
*> \param[out] IWORK
*> \verbatim
*>          IWORK is INTEGER array, dimension (N)
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          = 0:  successful exit
*>          < 0:  if INFO = -i, the i-th argument had an illegal value
*>          > 0:  if INFO = i, and i is
*>                <= N:  U(i,i) is exactly zero.  The factorization
*>                       has not been completed unless i = N, but the
*>                       factor U is exactly singular, so the solution
*>                       and error bounds could not be computed.
*>                       RCOND = 0 is returned.
*>                = N+1: U is nonsingular, but RCOND is less than machine
*>                       precision, meaning that the matrix is singular
*>                       to working precision.  Nevertheless, the
*>                       solution and error bounds are computed because
*>                       there are a number of situations where the
*>                       computed solution can be more accurate than the
*>                       value of RCOND would suggest.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup doubleGTsolve
*
*  =====================================================================
      SUBROUTINE DGTSVX( FACT, TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF,
     $                   DU2, IPIV, B, LDB, X, LDX, RCOND, FERR, BERR,
     $                   WORK, IWORK, INFO )
*
*  -- LAPACK driver routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      CHARACTER          FACT, TRANS
      INTEGER            INFO, LDB, LDX, N, NRHS
      DOUBLE PRECISION   RCOND
*     ..
*     .. Array Arguments ..
      INTEGER            IPIV( * ), IWORK( * )
      DOUBLE PRECISION   B( LDB, * ), BERR( * ), D( * ), DF( * ),
     $                   DL( * ), DLF( * ), DU( * ), DU2( * ), DUF( * ),
     $                   FERR( * ), WORK( * ), X( LDX, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO
      PARAMETER          ( ZERO = 0.0D+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            NOFACT, NOTRAN
      CHARACTER          NORM
      DOUBLE PRECISION   ANORM
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      DOUBLE PRECISION   DLAMCH, DLANGT
      EXTERNAL           LSAME, DLAMCH, DLANGT
*     ..
*     .. External Subroutines ..
      EXTERNAL           DCOPY, DGTCON, DGTRFS, DGTTRF, DGTTRS, DLACPY,
     $                   XERBLA
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MAX
*     ..
*     .. Executable Statements ..
*
      INFO = 0
      NOFACT = LSAME( FACT, 'N' )
      NOTRAN = LSAME( TRANS, 'N' )
      IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
         INFO = -1
      ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
     $         LSAME( TRANS, 'C' ) ) THEN
         INFO = -2
      ELSE IF( N.LT.0 ) THEN
         INFO = -3
      ELSE IF( NRHS.LT.0 ) THEN
         INFO = -4
      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
         INFO = -14
      ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
         INFO = -16
      END IF
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'DGTSVX', -INFO )
         RETURN
      END IF
*
      IF( NOFACT ) THEN
*
*        Compute the LU factorization of A.
*
         CALL DCOPY( N, D, 1, DF, 1 )
         IF( N.GT.1 ) THEN
            CALL DCOPY( N-1, DL, 1, DLF, 1 )
            CALL DCOPY( N-1, DU, 1, DUF, 1 )
         END IF
         CALL DGTTRF( N, DLF, DF, DUF, DU2, IPIV, INFO )
*
*        Return if INFO is non-zero.
*
         IF( INFO.GT.0 )THEN
            RCOND = ZERO
            RETURN
         END IF
      END IF
*
*     Compute the norm of the matrix A.
*
      IF( NOTRAN ) THEN
         NORM = '1'
      ELSE
         NORM = 'I'
      END IF
      ANORM = DLANGT( NORM, N, DL, D, DU )
*
*     Compute the reciprocal of the condition number of A.
*
      CALL DGTCON( NORM, N, DLF, DF, DUF, DU2, IPIV, ANORM, RCOND, WORK,
     $             IWORK, INFO )
*
*     Compute the solution vectors X.
*
      CALL DLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
      CALL DGTTRS( TRANS, N, NRHS, DLF, DF, DUF, DU2, IPIV, X, LDX,
     $             INFO )
*
*     Use iterative refinement to improve the computed solutions and
*     compute error bounds and backward error estimates for them.
*
      CALL DGTRFS( TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF, DU2, IPIV,
     $             B, LDB, X, LDX, FERR, BERR, WORK, IWORK, INFO )
*
*     Set INFO = N+1 if the matrix is singular to working precision.
*
      IF( RCOND.LT.DLAMCH( 'Epsilon' ) )
     $   INFO = N + 1
*
      RETURN
*
*     End of DGTSVX
*
      END

max_stars_repo_name: sergey-v-kuznetsov/lapack

max_issues_repo_name: sergey-v-kuznetsov/lapack

max_forks_repo_name: sergey-v-kuznetsov/lapack


text:
!
!               Tennessee Eastman Process Control Test Problem
!
!                    James J. Downs and Ernest F. Vogel
!
!                  Process and Control Systems Engineering
!                        Tennessee Eastman Company
!                              P.O. Box 511
!                          Kingsport,TN  37662
!
!  Reference:
!    "A Plant-Wide Industrial Process Control Problem"
!    Presented at the AIChE 1990 Annual Meeting
!    Industrial Challenge Problems in Process Control,Paper #24a
!    Chicago,Illinois,November 14,1990
!
!  Subroutines:
!
!    TEFUNC - Function evaluator to be called by integrator
!    TEINIT - Initialization
!    TESUBi - Utility subroutines, i=1,2,..,8
!
!
!  The process simulation has 50 states (NN=50).  If the user wishes to
!  integrate additional states, NN must be increased accordingly in the
!  calling program.  The additional states should be appended to the end
!  of the YY vector, e.g. YY(51),...  The additional derivatives should
!  be appended to the end of the YP vector, e.g. YP(51),...  To initialize
!  the new states and to calculate derivatives for them, we suggest
!  creating new function evaluator and initialization routines as follows.
!
!          C-----------------------------------------------
!          C
!                SUBROUTINE FUNC(NN,TIME,YY,YP)
!          C
!                INTEGER NN
!                DOUBLE PRECISION TIME, YY(NN), YP(NN)
!          C
!          C  Call the function evaluator for the process
!          C
!                CALL TEFUNC(NN,TIME,YY,YP)
!          C
!          C  Calculate derivatives for additional states
!          C
!                YP(51) = ....
!                YP(52) = ....
!                   .
!                   .
!                   .
!                YP(NN) = ....
!          C
!                RETURN
!                END
!          C
!          C-----------------------------------------------
!          C
!                SUBROUTINE INIT(NN,TIME,YY,YP)
!          C
!                INTEGER NN
!                DOUBLE PRECISION TIME, YY(NN), YP(NN)
!          C
!          C  Call the initialization for the process
!          C
!                CALL TEINIT(NN,TIME,YY,YP)
!          C
!          C  Initialize additional states
!          C
!                YY(51) = ....
!                YY(52) = ....
!                   .
!                   .
!                   .
!                YY(NN) = ....
!          C
!                RETURN
!                END
!          C
!          C-----------------------------------------------
!
!  Differences between the code and its description in the paper:
!
!  1.  Subroutine TEINIT has TIME in the argument list.  TEINIT sets TIME
!      to zero.
!
!  2.  There are 8 utility subroutines (TESUBi) rather than 5.
!
!  3.  Process disturbances 14 through 20 do NOT need to be used in
!      conjunction with another disturbance as stated in the paper.  All
!      disturbances can be used alone or in any combination.
!
!
!  Manipulated Variables
!
!    XMV(1)     A Feed Flow (stream 1)
!    XMV(2)     D Feed Flow (stream 2)
!    XMV(3)     E Feed Flow (stream 3)
!    XMV(4)     A and C Feed Flow (stream 4)
!    XMV(5)     Compressor Recycle Valve
!    XMV(6)     Purge Valve (stream 9)
!    XMV(7)     Separator Pot Liquid Flow (stream 10)
!    XMV(8)     Stripper Liquid Product Flow (stream 11)
!    XMV(9)     Stripper Steam Valve
!    XMV(10)    Reactor Cooling Water Flow
!    XMV(11)    Condenser Cooling Water Flow
!    XMV(12)    Agitator Speed
!
!  Continuous Process Measurements
!
!    XMEAS(1)   A Feed  (stream 1)                    kscmh
!    XMEAS(2)   D Feed  (stream 2)                    kg/hr
!    XMEAS(3)   E Feed  (stream 3)                    kg/hr
!    XMEAS(4)   A and C Feed  (stream 4)              kscmh
!    XMEAS(5)   Recycle Flow  (stream 8)              kscmh
!    XMEAS(6)   Reactor Feed Rate  (stream 6)         kscmh
!    XMEAS(7)   Reactor Pressure                      kPa gauge
!    XMEAS(8)   Reactor Level                         %
!    XMEAS(9)   Reactor Temperature                   Deg C
!    XMEAS(10)  Purge Rate (stream 9)                 kscmh
!    XMEAS(11)  Product Sep Temp                      Deg C
!    XMEAS(12)  Product Sep Level                     %
!    XMEAS(13)  Prod Sep Pressure                     kPa gauge
!    XMEAS(14)  Prod Sep Underflow (stream 10)        m3/hr
!    XMEAS(15)  Stripper Level                        %
!    XMEAS(16)  Stripper Pressure                     kPa gauge
!    XMEAS(17)  Stripper Underflow (stream 11)        m3/hr
!    XMEAS(18)  Stripper Temperature                  Deg C
!    XMEAS(19)  Stripper Steam Flow                   kg/hr
!    XMEAS(20)  Compressor Work                       kW
!    XMEAS(21)  Reactor Cooling Water Outlet Temp     Deg C
!    XMEAS(22)  Separator Cooling Water Outlet Temp   Deg C
!
!  Sampled Process Measurements
!
!    Reactor Feed Analysis (Stream 6)
!        Sampling Frequency = 0.1 hr
!        Dead Time = 0.1 hr
!        Mole %
!    XMEAS(23)   Component A
!    XMEAS(24)   Component B
!    XMEAS(25)   Component C
!    XMEAS(26)   Component D
!    XMEAS(27)   Component E
!    XMEAS(28)   Component F
!
!    Purge Gas Analysis (Stream 9)
!        Sampling Frequency = 0.1 hr
!        Dead Time = 0.1 hr
!        Mole %
!    XMEAS(29)   Component A
!    XMEAS(30)   Component B
!    XMEAS(31)   Component C
!    XMEAS(32)   Component D
!    XMEAS(33)   Component E
!    XMEAS(34)   Component F
!    XMEAS(35)   Component G
!    XMEAS(36)   Component H
!
!    Product Analysis (Stream 11)
!        Sampling Frequency = 0.25 hr
!        Dead Time = 0.25 hr
!        Mole %
!    XMEAS(37)   Component D
!    XMEAS(38)   Component E
!    XMEAS(39)   Component F
!    XMEAS(40)   Component G
!    XMEAS(41)   Component H
!
!  Process Disturbances
!
!    IDV(1)   A/C Feed Ratio, B Composition Constant (Stream 4)          Step
!    IDV(2)   B Composition, A/C Ratio Constant (Stream 4)               Step
!    IDV(3)   D Feed Temperature (Stream 2)                              Step
!    IDV(4)   Reactor Cooling Water Inlet Temperature                    Step
!    IDV(5)   Condenser Cooling Water Inlet Temperature                  Step
!    IDV(6)   A Feed Loss (Stream 1)                                     Step
!    IDV(7)   C Header Pressure Loss - Reduced Availability (Stream 4)   Step
!    IDV(8)   A, B, C Feed Composition (Stream 4)            Random Variation
!    IDV(9)   D Feed Temperature (Stream 2)                  Random Variation
!    IDV(10)  C Feed Temperature (Stream 4)                  Random Variation
!    IDV(11)  Reactor Cooling Water Inlet Temperature        Random Variation
!    IDV(12)  Condenser Cooling Water Inlet Temperature      Random Variation
!    IDV(13)  Reaction Kinetics                                    Slow Drift
!    IDV(14)  Reactor Cooling Water Valve                            Sticking
!    IDV(15)  Condenser Cooling Water Valve                          Sticking
!    IDV(16)  Unknown
!    IDV(17)  Unknown
!    IDV(18)  Unknown
!    IDV(19)  Unknown
!    IDV(20)  Unknown
!
!
!=============================================================================
!

!-----------------------------------------------------------------------
!> Function Evaluator.
!!
!!  ## Input
!!      @param[in]  `nn`    Number of differential equations
!!      @param[in]  `time`  Current time(hrs)
!!      @param[in]  `yy`    Current state values
!!      @param[out] `yp`    Current derivative values
!!
!!  ### common:
!!      + /pv/
!!      + /dvec/
!!      + /teproc/
!!      + /wlk/
!!      + /const/
!!
Subroutine tefunc(nn, time, yy, yp)
    Double Precision yy(nn), yp(nn)
!
!  MEASUREMENT AND VALVE COMMON BLOCK
!
    Double Precision xmeas, xmv
    Integer isd
    Common /pv/ xmeas(41), xmv(12), isd
!
!   DISTURBANCE VECTOR COMMON BLOCK
!
    Integer idv
    Common /dvec/ idv(20)
    Double Precision uclr, ucvr, utlr, utvr, xlr, xvr, &
        etr, esr, tcr, tkr, dlr, vlr, vvr, vtr, ptr, ppr, &
        crxr, rr, rh, fwr, twr, qur, hwr, uar, ucls, ucvs, &
        utls, utvs, xls, xvs, ets, ess, tcs, tks, dls, vls, &
        vvs, vts, pts, pps, fws, tws, qus, hws, uclc, utlc, &
        xlc, etc, esc, tcc, dlc, vlc, vtc, quc, ucvv, utvv, &
        xvv, etv, esv, tcv, tkv, vtv, ptv, vcv, vrng, vtau, &
        ftm, fcm, xst, xmws, hst, tst, sfr, cpflmx, cpprmx, &
        cpdh, tcwr, tcws, htr, agsp, xdel, xns, tgas, tprod, vst
    Integer ivst
    Common /teproc/ uclr(8), ucvr(8), utlr, utvr, xlr(8), xvr(8), &
        etr, esr, tcr, tkr, dlr, vlr, vvr, vtr, ptr, &
        ppr(8), crxr(8), rr(4), rh, fwr, twr, qur, hwr, uar, &
        ucls(8), ucvs(8), utls, utvs, xls(8), xvs(8), ets, ess, &
        tcs, tks, dls, vls, vvs, vts, pts, pps(8), fws, tws, &
        qus, hws, uclc(8), utlc, xlc(8), etc, esc, tcc, dlc, &
        vlc, vtc, quc, ucvv(8), utvv, xvv(8), etv, esv, &
        tcv, tkv, vtv, ptv, vcv(12), vrng(12), vtau(12), ftm(13), &
        fcm(8, 13), xst(8, 13), xmws(13), hst(13), tst(13), sfr(8), &
        cpflmx, cpprmx, cpdh, tcwr, tcws, htr(3), agsp, &
        xdel(41), xns(41), tgas, tprod, vst(12), ivst(12)
    Integer idvwlk
    Double Precision adist, bdist, cdist, ddist, tlast, tnext, &
        hspan, hzero, sspan, szero, spspan
    Common /wlk/ adist(12), bdist(12), cdist(12), ddist(12), &
        tlast(12), tnext(12), hspan(12), hzero(12), &
        sspan(12), szero(12), spspan(12), idvwlk(12)
    Double Precision avp, bvp, cvp, ah, bh, ch, ag, bg, cg, &
        av, ad, bd, cd, xmw
    Common /const/ avp(8), bvp(8), cvp(8), ah(8), bh(8), ch(8), &
        ag(8), bg(8), cg(8), av(8), ad(8), bd(8), cd(8), xmw(8)
    Integer nn, i
    Double Precision rg, vpr, fin(8), time, flms, dlp, pr, &
        flcoef, uas, uac, vovrl, uarlev, vpos(12), &
        xmns, xcmp(41), tmpfac, r1f, r2f, hwlk, swlk, spwlk, &
        tesub7, tesub8

!!  Function Body --------------------------------
!
! Limiting of Disturbance Activations
    Do i = 1, 20
        If (idv(i)>0) Then
            idv(i) = 1
        Else
            idv(i) = 0
        End If
    End Do

! Assignment of Disturbance Activations
    idvwlk(1) = idv(8)
    idvwlk(2) = idv(8)
    idvwlk(3) = idv(9)
    idvwlk(4) = idv(10)
    idvwlk(5) = idv(11)
    idvwlk(6) = idv(12)
    idvwlk(7) = idv(13)
    idvwlk(8) = idv(13)
    idvwlk(9) = idv(16)
    idvwlk(10) = idv(17)
    idvwlk(11) = idv(18)
    idvwlk(12) = idv(20)

! Recalculation of Disturbance Process Parameters -
!   Update (1 - 9/13 - 20)
!
    Do i = 1, 9
        If (time>=tnext(i)) Then
            hwlk = tnext(i) - tlast(i)
            swlk  = adist(i) + hwlk * ( &
                bdist(i)+hwlk*(cdist(i)+hwlk*ddist(i)) &
            )
            spwlk = bdist(i) + hwlk * ( &
                2.D0*cdist(i)+3.D0*hwlk*ddist(i) &
            )
            tlast(i) = tnext(i)
            Call tesub5(swlk, spwlk, &
                    adist(i), bdist(i), cdist(i), ddist(i), &
                    tlast(i), tnext(i), &
                    hspan(i), hzero(i), sspan(i), szero(i), &
                    spspan(i), idvwlk(i))
        End If
    End Do

! Recalculation of Disturbance Process Parameters -
!   Determination of Processes to be Updated (10 - 12)
!
    Do i = 10, 12
        If (time>=tnext(i)) Then
            hwlk = tnext(i) - tlast(i)
            swlk  = adist(i) + hwlk * ( &
                bdist(i)+hwlk*(cdist(i)+hwlk*ddist(i)) &
            )
            spwlk = bdist(i) + hwlk * ( &
                2.D0*cdist(i)+3.D0*hwlk*ddist(i) &
            )
            tlast(i) = tnext(i)
            If (swlk>0.1D0) Then
                adist(i) = swlk
                bdist(i) = spwlk
                cdist(i) = -(3.D0*swlk+0.2D0*spwlk)/0.01D0
                ddist(i) = (2.D0*swlk+0.1D0*spwlk)/0.001D0
                tnext(i) = tlast(i) + 0.1D0
            Else
                isd = -1
                hwlk = hspan(i)*tesub7(isd) + hzero(i)
                adist(i) = 0.D0
                bdist(i) = 0.D0
                cdist(i) = dble(idvwlk(i))/hwlk**2
                ddist(i) = 0.D0
                tnext(i) = tlast(i) + hwlk
            End If
        End If
    End Do

! Initilization of Disturbance Processes Parameters
    If (time==0.D0) Then
        Do i = 1, 12
            adist(i) = szero(i)
            bdist(i) = 0.D0
            cdist(i) = 0.D0
            ddist(i) = 0.D0
            tlast(i) = 0.0D0
            tnext(i) = 0.1D0
        End Do
    End If

! Determination of Disturbed Values
    xst(1, 4) = tesub8(1, time) - idv(1)*0.03D0 - idv(2)*2.43719D-3
    xst(2, 4) = tesub8(2, time) + idv(2)*0.005D0
    xst(3, 4) = 1.D0 - xst(1, 4) - xst(2, 4)
    tst(1) = tesub8(3, time) + idv(3)*5.D0
    tst(4) = tesub8(4, time)
    tcwr = tesub8(5, time) + idv(4)*5.D0
    tcws = tesub8(6, time) + idv(5)*5.D0
    r1f = tesub8(7, time)
    r2f = tesub8(8, time)

! Retrieving of Current States
    Do i = 1, 3
        ucvr(i) = yy(i)
        ucvs(i) = yy(i+9)
        uclr(i) = 0.0
        ucls(i) = 0.0
    End Do
    Do i = 4, 8
        uclr(i) = yy(i)
        ucls(i) = yy(i+9)
    End Do
    Do i = 1, 8
        uclc(i) = yy(i+18)
        ucvv(i) = yy(i+27)
    End Do
    etr = yy(9)
    ets = yy(18)
    etc = yy(27)
    etv = yy(36)
    twr = yy(37)
    tws = yy(38)
    Do i = 1, 12
        vpos(i) = yy(i+38)
    End Do

! Calculation of Collective Holdup
    utlr = 0.0
    utls = 0.0
    utlc = 0.0
    utvv = 0.0
    Do i = 1, 8
        utlr = utlr + uclr(i)
        utls = utls + ucls(i)
        utlc = utlc + uclc(i)
        utvv = utvv + ucvv(i)
    End Do
    Do i = 1, 8
        xlr(i) = uclr(i)/utlr
        xls(i) = ucls(i)/utls
        xlc(i) = uclc(i)/utlc
        xvv(i) = ucvv(i)/utvv
    End Do

    esr = etr/utlr
    ess = ets/utls
    esc = etc/utlc
    esv = etv/utvv

! Calculation of Temperatures
    Call tesub2(xlr, tcr, esr, 0)
    tkr = tcr + 273.15
    Call tesub2(xls, tcs, ess, 0)
    tks = tcs + 273.15
    Call tesub2(xlc, tcc, esc, 0)
    Call tesub2(xvv, tcv, esv, 2)
    tkv = tcv + 273.15

! Calculation of Densities
    Call tesub4(xlr, tcr, dlr)
    Call tesub4(xls, tcs, dls)
    Call tesub4(xlc, tcc, dlc)

! Calculation of Volume of Liquid and Vapor Phase
    vlr = utlr/dlr
    vls = utls/dls
    vlc = utlc/dlc
    vvr = vtr - vlr
    vvs = vts - vls

    rg = 998.9

! Calculation of Pressure
    ptr = 0.0
    pts = 0.0

    Do i = 1, 3
        ppr(i) = ucvr(i)*rg*tkr/vvr
        ptr = ptr + ppr(i)
        pps(i) = ucvs(i)*rg*tks/vvs
        pts = pts + pps(i)
    End Do
    Do i = 4, 8
        vpr = dexp(avp(i)+bvp(i)/(tcr+cvp(i)))
        ppr(i) = vpr*xlr(i)
        ptr = ptr + ppr(i)
        vpr = dexp(avp(i)+bvp(i)/(tcs+cvp(i)))
        pps(i) = vpr*xls(i)
        pts = pts + pps(i)
    End Do

! Calculation of Component Concentration in Vapor Phase
!   (Reactor and Separator)
    ptv = utvv*rg*tkv/vtv
    Do i = 1, 8
        xvr(i) = ppr(i)/ptr
        xvs(i) = pps(i)/pts
    End Do

! Calculation of Collective Holdup of Components in Vapor Phase
!   (Reactor and Separator)
    utvr = ptr*vvr/rg/tkr
    utvs = pts*vvs/rg/tks

! Calculation of Single Holdup of Components in Vapor Phase
!   (Reactor and Separator)
    Do i = 4, 8
        ucvr(i) = utvr*xvr(i)
        ucvs(i) = utvs*xvs(i)
    End Do

! Reaction kinetics
    rr(1) = dexp(31.5859536-40000.0/1.987/tkr)*r1f
    rr(2) = dexp(3.00094014-20000.0/1.987/tkr)*r2f
    rr(3) = dexp(53.4060443-60000.0/1.987/tkr)
    rr(4) = rr(3)*0.767488334D0
    If (ppr(1)>0.0 .And. ppr(3)>0.0) Then
        r1f = ppr(1)**1.1544
        r2f = ppr(3)**0.3735
        rr(1) = rr(1)*r1f*r2f*ppr(4)
        rr(2) = rr(2)*r1f*r2f*ppr(5)
    Else
        rr(1) = 0.0
        rr(2) = 0.0
    End If
    rr(3) = rr(3)*ppr(1)*ppr(5)
    rr(4) = rr(4)*ppr(1)*ppr(4)
    Do i = 1, 4
        rr(i) = rr(i)*vvr
    End Do

! Consumption and Creation of Components in Reactor
    crxr(1) = -rr(1) - rr(2) - rr(3)
    crxr(3) = -rr(1) - rr(2)
    crxr(4) = -rr(1) - 1.5D0*rr(4)
    crxr(5) = -rr(2) - rr(3)
    crxr(6) = rr(3) + rr(4)
    crxr(7) = rr(1)
    crxr(8) = rr(2)
    rh = rr(1)*htr(1) + rr(2)*htr(2)

    xmws(1) = 0.0
    xmws(2) = 0.0
    xmws(6) = 0.0
    xmws(8) = 0.0
    xmws(9) = 0.0
    xmws(10) = 0.0

    Do i = 1, 8
        xst(i, 6) = xvv(i)
        xst(i, 8) = xvr(i)
        xst(i, 9) = xvs(i)
        xst(i, 10) = xvs(i)
        xst(i, 11) = xls(i)
        xst(i, 13) = xlc(i)

        xmws(1) = xmws(1) + xst(i, 1)*xmw(i)
        xmws(2) = xmws(2) + xst(i, 2)*xmw(i)
        xmws(6) = xmws(6) + xst(i, 6)*xmw(i)
        xmws(8) = xmws(8) + xst(i, 8)*xmw(i)
        xmws(9) = xmws(9) + xst(i, 9)*xmw(i)
        xmws(10) = xmws(10) + xst(i, 10)*xmw(i)
    End Do

    tst(6) = tcv
    tst(8) = tcr
    tst(9) = tcs
    tst(10) = tcs
    tst(11) = tcs
    tst(13) = tcc

    Call tesub1(xst(1,1), tst(1), hst(1), 1)
    Call tesub1(xst(1,2), tst(2), hst(2), 1)
    Call tesub1(xst(1,3), tst(3), hst(3), 1)
    Call tesub1(xst(1,4), tst(4), hst(4), 1)
    Call tesub1(xst(1,6), tst(6), hst(6), 1)
    Call tesub1(xst(1,8), tst(8), hst(8), 1)
    Call tesub1(xst(1,9), tst(9), hst(9), 1)
    hst(10) = hst(9)
    Call tesub1(xst(1,11), tst(11), hst(11), 0)
    Call tesub1(xst(1,13), tst(13), hst(13), 0)

    ftm(1) = vpos(1)*vrng(1)/100.0
    ftm(2) = vpos(2)*vrng(2)/100.0
    ftm(3) = vpos(3)*(1.D0-idv(6))*vrng(3)/100.0
    ftm(4) = vpos(4)*(1.D0-idv(7)*0.2D0)*vrng(4)/100.0 + 1.D-10
    ftm(11) = vpos(7)*vrng(7)/100.0
    ftm(13) = vpos(8)*vrng(8)/100.0

    uac = vpos(9)*vrng(9)*(1.D0+tesub8(9,time))/100.0
    fwr = vpos(10)*vrng(10)/100.0
    fws = vpos(11)*vrng(11)/100.0
    agsp = (vpos(12)+150.0)/100.0

    dlp = ptv - ptr
    If (dlp<0.0) dlp = 0.0
    flms = 1937.6D0*dsqrt(dlp)
    ftm(6) = flms/xmws(6)
    dlp = ptr - pts
    If (dlp<0.0) dlp = 0.0
    flms = 4574.21D0*dsqrt(dlp)*(1.D0-0.25D0*tesub8(12,time))
    ftm(8) = flms/xmws(8)
    dlp = pts - 760.0
    If (dlp<0.0) dlp = 0.0
    flms = vpos(6)*0.151169D0*dsqrt(dlp)
    ftm(10) = flms/xmws(10)
    pr = ptv/pts
    If (pr<1.0) pr = 1.0
    If (pr>cpprmx) pr = cpprmx
    flcoef = cpflmx/1.197D0

    flms = cpflmx + flcoef*(1.0-pr**3)
    cpdh = flms * (tcs+273.15D0) * 1.8D-6 * 1.9872D0 * (ptv-pts) &
            / (xmws(9)*pts)
    dlp = ptv - pts
    If (dlp<0.0) dlp = 0.0
    flms = flms - vpos(5)*53.349D0*dsqrt(dlp)
    If (flms<1.D-3) flms = 1.D-3
    ftm(9) = flms/xmws(9)
    hst(9) = hst(9) + cpdh/ftm(9)

    Do i = 1, 8
        fcm(i, 1) = xst(i, 1)*ftm(1)
        fcm(i, 2) = xst(i, 2)*ftm(2)
        fcm(i, 3) = xst(i, 3)*ftm(3)
        fcm(i, 4) = xst(i, 4)*ftm(4)
        fcm(i, 6) = xst(i, 6)*ftm(6)
        fcm(i, 8) = xst(i, 8)*ftm(8)
        fcm(i, 9) = xst(i, 9)*ftm(9)
        fcm(i, 10) = xst(i, 10)*ftm(10)
        fcm(i, 11) = xst(i, 11)*ftm(11)
        fcm(i, 13) = xst(i, 13)*ftm(13)
    End Do

    If (ftm(11)>0.1) Then
        If (tcc>170.) Then
            tmpfac = tcc - 120.262
        Else If (tcc<5.292) Then
            tmpfac = 0.1
        Else
            tmpfac = 363.744/(177.-tcc) - 2.22579488
        End If
        vovrl = ftm(4)/ftm(11)*tmpfac
        sfr(4) = 8.5010*vovrl/(1.0+8.5010*vovrl)
        sfr(5) = 11.402*vovrl/(1.0+11.402*vovrl)
        sfr(6) = 11.795*vovrl/(1.0+11.795*vovrl)
        sfr(7) = 0.0480*vovrl/(1.0+0.0480*vovrl)
        sfr(8) = 0.0242*vovrl/(1.0+0.0242*vovrl)
    Else
        sfr(4) = 0.9999
        sfr(5) = 0.999
        sfr(6) = 0.999
        sfr(7) = 0.99
        sfr(8) = 0.98
    End If ! (ftm(11)>0.1)

    Do i = 1, 8
        fin(i) = 0.0
        fin(i) = fin(i) + fcm(i, 4)
        fin(i) = fin(i) + fcm(i, 11)
    End Do
    ftm(5) = 0.0
    ftm(12) = 0.0
    Do i = 1, 8
        fcm(i, 5) = sfr(i)*fin(i)
        fcm(i, 12) = fin(i) - fcm(i, 5)
        ftm(5) = ftm(5) + fcm(i, 5)
        ftm(12) = ftm(12) + fcm(i, 12)
    End Do
    Do i = 1, 8
        xst(i, 5) = fcm(i, 5)/ftm(5)
        xst(i, 12) = fcm(i, 12)/ftm(12)
    End Do
    tst(5) = tcc
    tst(12) = tcc
    Call tesub1(xst(1,5), tst(5), hst(5), 1)
    Call tesub1(xst(1,12), tst(12), hst(12), 0)
    ftm(7) = ftm(6)
    hst(7) = hst(6)
    tst(7) = tst(6)
    Do i = 1, 8
        xst(i, 7) = xst(i, 6)
        fcm(i, 7) = fcm(i, 6)
    End Do

! Calculation of Heat Transfer in Reactor
    If (vlr/7.8>50.0) Then
        uarlev = 1.0
    Else If (vlr/7.8<10.0) Then
        uarlev = 0.0
    Else
        uarlev = 0.025*vlr/7.8 - 0.25
    End If

    uar = uarlev*(-0.5*agsp**2+2.75*agsp-2.5)*855490.D-6
    qur = uar*(twr-tcr)*(1.D0-0.35D0*tesub8(10,time))

! Calculation of Heat Transfer in Condenser (Separator)
    uas = 0.404655*(1.0-1.0/(1.0+(ftm(8)/3528.73)**4))
    qus = uas*(tws-tst(8))*(1.D0-0.25D0*tesub8(11,time))

! Calculation of Heat Transfer in Stripper
    quc = 0.D0
    If (tcc<100.) quc = uac*(100.0-tcc)

! Setting of Measured Values
    xmeas(1) = ftm(3)*0.359/35.3145
    xmeas(2) = ftm(1)*xmws(1)*0.454
    xmeas(3) = ftm(2)*xmws(2)*0.454
    xmeas(4) = ftm(4)*0.359/35.3145
    xmeas(5) = ftm(9)*0.359/35.3145
    xmeas(6) = ftm(6)*0.359/35.3145
    xmeas(7) = (ptr-760.0)/760.0*101.325
    xmeas(8) = (vlr-84.6)/666.7*100.0
    xmeas(9) = tcr
    xmeas(10) = ftm(10)*0.359/35.3145
    xmeas(11) = tcs
    xmeas(12) = (vls-27.5)/290.0*100.0
    xmeas(13) = (pts-760.0)/760.0*101.325
    xmeas(14) = ftm(11)/dls/35.3145
    xmeas(15) = (vlc-78.25)/vtc*100.0
    xmeas(16) = (ptv-760.0)/760.0*101.325
    xmeas(17) = ftm(13)/dlc/35.3145
    xmeas(18) = tcc
    xmeas(19) = quc*1.04D3*0.454
    xmeas(20) = cpdh*0.0003927D6
    xmeas(20) = cpdh*0.29307D3
    xmeas(21) = twr
    xmeas(22) = tws

! Checking of Shut-Down-Constraints
    isd = 0
    If (xmeas(7)>3000.0) isd = 1
    If (vlr/35.3145>24.0) isd = 1
    If (vlr/35.3145<2.0) isd = 1
    If (xmeas(9)>175.0) isd = 1
    If (vls/35.3145>12.0) isd = 1
    If (vls/35.3145<1.0) isd = 1
    If (vlc/35.3145>8.0) isd = 1
    If (vlc/35.3145<1.0) isd = 1

! Adding of Measurement Noise
    If (time>0.0 .And. isd==0) Then
        Do i = 1, 22
            Call tesub6(xns(i), xmns)
            xmeas(i) = xmeas(i) + xmns
        End Do
    End If

! Analyzer Outputs
    xcmp(23) = xst(1, 7)*100.0
    xcmp(24) = xst(2, 7)*100.0
    xcmp(25) = xst(3, 7)*100.0
    xcmp(26) = xst(4, 7)*100.0
    xcmp(27) = xst(5, 7)*100.0
    xcmp(28) = xst(6, 7)*100.0
    xcmp(29) = xst(1, 10)*100.0
    xcmp(30) = xst(2, 10)*100.0
    xcmp(31) = xst(3, 10)*100.0
    xcmp(32) = xst(4, 10)*100.0
    xcmp(33) = xst(5, 10)*100.0
    xcmp(34) = xst(6, 10)*100.0
    xcmp(35) = xst(7, 10)*100.0
    xcmp(36) = xst(8, 10)*100.0
    xcmp(37) = xst(4, 13)*100.0
    xcmp(38) = xst(5, 13)*100.0
    xcmp(39) = xst(6, 13)*100.0
    xcmp(40) = xst(7, 13)*100.0
    xcmp(41) = xst(8, 13)*100.0

    If (time==0.D0) Then
        Do i = 23, 41
            xdel(i) = xcmp(i)
            xmeas(i) = xcmp(i)
        End Do
        tgas = 0.1
        tprod = 0.25
    End If

    If (time>=tgas) Then
        Do i = 23, 36
            xmeas(i) = xdel(i)
            Call tesub6(xns(i), xmns)
            xmeas(i) = xmeas(i) + xmns
            xdel(i) = xcmp(i)
        End Do
        tgas = tgas + 0.1
    End If

    If (time>=tprod) Then
        Do i = 37, 41
            xmeas(i) = xdel(i)
            Call tesub6(xns(i), xmns)
            xmeas(i) = xmeas(i) + xmns
            xdel(i) = xcmp(i)
        End Do
        tprod = tprod + 0.25
    End If

! Calculation of State Derivative
    Do i = 1, 8
        yp(i)    = fcm(i, 7)  - fcm(i, 8) + crxr(i)
        yp(i+9)  = fcm(i, 8)  - fcm(i, 9) - fcm(i, 10) - fcm(i, 11)
        yp(i+18) = fcm(i, 12) - fcm(i, 13)
        yp(i+27) = fcm(i, 1)  + fcm(i, 2) + fcm(i, 3) &
                    + fcm(i, 5) + fcm(i, 9) - fcm(i, 6)
    End Do
    yp(9) = hst(7)*ftm(7) - hst(8)*ftm(8) + rh + qur

!   Here is the "correct" version of the separator energy balance:
!
!	    YP(18) = HST(8)*FTM(8) - (HST(9)*FTM(9)-cpdh)
!               - HST(10)*FTM(10) - HST(11)*FTM(11) + QUS
!
!   Here is the original version
!
    yp(18) = hst(8)*ftm(8) - hst(9)*ftm(9) &
            - hst(10)*ftm(10) - hst(11)*ftm(11) + qus

    yp(27) = hst(4)*ftm(4) + hst(11)*ftm(11) &
            - hst(5)*ftm(5) - hst(13)*ftm(13) + quc

    yp(36) = hst(1)*ftm(1) + hst(2)*ftm(2) + hst(3)*ftm(3) &
            + hst(5)*ftm(5) + hst(9)*ftm(9) - hst(6)*ftm(6)

    yp(37) = (fwr*500.53*(tcwr-twr)-qur*1.D6/1.8)/hwr
    yp(38) = (fws*500.53*(tcws-tws)-qus*1.D6/1.8)/ hws

    ivst(10) = idv(14)
    ivst(11) = idv(15)
    ivst(5) = idv(19)
    ivst(7) = idv(19)
    ivst(8) = idv(19)
    ivst(9) = idv(19)

    Do i = 1, 12
        If ( time == 0.D0 &
            .Or. dabs(vcv(i)-xmv(i)) > vst(i) * ivst(i) &
        ) vcv(i) = xmv(i)
        If (vcv(i) < 0.0)   vcv(i) = 0.0
        If (vcv(i) > 100.0) vcv(i) = 100.0

        yp(i+38) = (vcv(i)-vpos(i)) / vtau(i)
    End Do

    If (time>0.0 .And. isd/=0) Then
        Do i = 1, nn
            yp(i) = 0.0
        End Do
    End If

    Return
End Subroutine tefunc

!-----------------------------------------------------------------------
!> Initialization TEP.
!!
!!  ## Input
!!      @param[in]  `nn`    Number of differential equations
!!      @param[out] `time`  Current time(hrs)
!!      @param[out] `yy`    Current state values
!!      @param[out] `yp`    Current derivative values
!!
!!  ### common:
!!      + /pv/
!!      + /dvec/
!!      + /randsd/
!!      + /teproc/
!!      + /wlk/
!!      + /const/
!!
Subroutine teinit(nn, time, yy, yp)
    Double Precision xmeas, xmv
    Common /pv/ xmeas(41), xmv(12)
    Integer idv
    Common /dvec/ idv(20)
    Double Precision g
    Common /randsd/ g
    Double Precision uclr, ucvr, utlr, utvr, xlr, xvr, etr, esr, &
        tcr, tkr, dlr, vlr, vvr, vtr, ptr, ppr, crxr, rr, rh, &
        fwr, twr, qur, hwr, uar, ucls, ucvs, utls, utvs, &
        xls, xvs, ets, ess, tcs, tks, dls, vls, vvs, vts, &
        pts, pps, fws, tws, qus, hws, uclc, utlc, xlc, &
        etc, esc, tcc, dlc, vlc, vtc, quc, ucvv, utvv, xvv, &
        etv, esv, tcv, tkv, vtv, ptv, vcv, vrng, vtau, ftm, fcm, &
        xst, xmws, hst, tst, sfr, cpflmx, cpprmx, cpdh, tcwr, tcws, &
        htr, agsp, xdel, xns, tgas, tprod, vst
    Integer ivst
    Common /teproc/ uclr(8), ucvr(8), utlr, utvr, xlr(8), xvr(8), &
        etr, esr, tcr, tkr, dlr, vlr, vvr, vtr, ptr, ppr(8), &
        crxr(8), rr(4), rh, fwr, twr, qur, hwr, uar, &
        ucls(8), ucvs(8), utls, utvs, xls(8), xvs(8), &
        ets, ess, tcs, tks, dls, vls, vvs, vts, pts, pps(8), &
        fws, tws, qus, hws, uclc(8), utlc, xlc(8), etc, esc, &
        tcc, dlc, vlc, vtc, quc, ucvv(8), utvv, xvv(8), etv, &
        esv, tcv, tkv, vtv, ptv, vcv(12), vrng(12), vtau(12), &
        ftm(13), fcm(8, 13), xst(8, 13), xmws(13), hst(13), &
        tst(13), sfr(8), cpflmx, cpprmx, cpdh, tcwr, tcws, &
        htr(3), agsp, xdel(41), xns(41), tgas, tprod, &
        vst(12), ivst(12)
    Integer idvwlk
    Double Precision adist, bdist, cdist, ddist, tlast, tnext, &
        hspan, hzero, sspan, szero, spspan
    Common /wlk/ adist(12), bdist(12), cdist(12), ddist(12), &
        tlast(12), tnext(12), hspan(12), hzero(12), &
        sspan(12), szero(12), spspan(12), idvwlk(12)
    Double Precision avp, bvp, cvp, ah, bh, ch, &
        ag, bg, cg, av, ad, bd, cd, xmw
    Common /const/ avp(8), bvp(8), cvp(8), ah(8), bh(8), ch(8), &
        ag(8), bg(8), cg(8), av(8), ad(8), bd(8), cd(8), xmw(8)

    Integer i, nn
    Double Precision yy(nn), yp(nn), time

!! Common init seq
!   1. Common /const/
!   2. Common /pv/
!   3. Common /teproc/
!   4. Common /randsd/
!   5. Common /dvec/
!   6. Common /wlk/
!

! Init Common /const/ ---------------------------
!   xmw(8)
!   avp(8), bvp(8), cvp(8),
!   ad(8), bd(8), cd(8),
!   ah(8), bh(8), ch(8),
!   av(8),
!   ag(8), bg(8), cg(8),
!
    xmw(1) = 2.0
    xmw(2) = 25.4
    xmw(3) = 28.0
    xmw(4) = 32.0
    xmw(5) = 46.0
    xmw(6) = 48.0
    xmw(7) = 62.0
    xmw(8) = 76.0

    avp(1) = 0.0
    avp(2) = 0.0
    avp(3) = 0.0
    avp(4) = 15.92
    avp(5) = 16.35
    avp(6) = 16.35
    avp(7) = 16.43
    avp(8) = 17.21
    bvp(1) = 0.0
    bvp(2) = 0.0
    bvp(3) = 0.0
    bvp(4) = -1444.0
    bvp(5) = -2114.0
    bvp(6) = -2114.0
    bvp(7) = -2748.0
    bvp(8) = -3318.0
    cvp(1) = 0.0
    cvp(2) = 0.0
    cvp(3) = 0.0
    cvp(4) = 259.0
    cvp(5) = 265.5
    cvp(6) = 265.5
    cvp(7) = 232.9
    cvp(8) = 249.6

    ad(1) = 1.0
    ad(2) = 1.0
    ad(3) = 1.0
    ad(4) = 23.3
    ad(5) = 33.9
    ad(6) = 32.8
    ad(7) = 49.9
    ad(8) = 50.5
    bd(1) = 0.0
    bd(2) = 0.0
    bd(3) = 0.0
    bd(4) = -0.0700
    bd(5) = -0.0957
    bd(6) = -0.0995
    bd(7) = -0.0191
    bd(8) = -0.0541
    cd(1) = 0.0
    cd(2) = 0.0
    cd(3) = 0.0
    cd(4) = -0.0002
    cd(5) = -0.000152
    cd(6) = -0.000233
    cd(7) = -0.000425
    cd(8) = -0.000150

    ah(1) = 1.0D-6
    ah(2) = 1.0D-6
    ah(3) = 1.0D-6
    ah(4) = 0.960D-6
    ah(5) = 0.573D-6
    ah(6) = 0.652D-6
    ah(7) = 0.515D-6
    ah(8) = 0.471D-6
    bh(1) = 0.0
    bh(2) = 0.0
    bh(3) = 0.0
    bh(4) = 8.70D-9
    bh(5) = 2.41D-9
    bh(6) = 2.18D-9
    bh(7) = 5.65D-10
    bh(8) = 8.70D-10
    ch(1) = 0.0
    ch(2) = 0.0
    ch(3) = 0.0
    ch(4) = 4.81D-11
    ch(5) = 1.82D-11
    ch(6) = 1.94D-11
    ch(7) = 3.82D-12
    ch(8) = 2.62D-12

    av(1) = 1.0D-6
    av(2) = 1.0D-6
    av(3) = 1.0D-6
    av(4) = 86.7D-6
    av(5) = 160.D-6
    av(6) = 160.D-6
    av(7) = 225.D-6
    av(8) = 209.D-6

    ag(1) = 3.411D-6
    ag(2) = 0.3799D-6
    ag(3) = 0.2491D-6
    ag(4) = 0.3567D-6
    ag(5) = 0.3463D-6
    ag(6) = 0.3930D-6
    ag(7) = 0.170D-6
    ag(8) = 0.150D-6
    bg(1) = 7.18D-10
    bg(2) = 1.08D-9
    bg(3) = 1.36D-11
    bg(4) = 8.51D-10
    bg(5) = 8.96D-10
    bg(6) = 1.02D-9
    bg(7) = 0.D0
    bg(8) = 0.D0
    cg(1) = 6.0D-13
    cg(2) = -3.98D-13
    cg(3) = -3.93D-14
    cg(4) = -3.12D-13
    cg(5) = -3.27D-13
    cg(6) = -3.12D-13
    cg(7) = 0.D0
    cg(8) = 0.D0
!
! END Init Common /const/ -----------------------

    yy(1) = 10.40491389
    yy(2) = 4.363996017
    yy(3) = 7.570059737
    yy(4) = 0.4230042431
    yy(5) = 24.15513437
    yy(6) = 2.942597645
    yy(7) = 154.3770655
    yy(8) = 159.1865960
    yy(9) = 2.808522723
    yy(10) = 63.75581199
    yy(11) = 26.74026066
    yy(12) = 46.38532432
    yy(13) = 0.2464521543
    yy(14) = 15.20484404
    yy(15) = 1.852266172
    yy(16) = 52.44639459
    yy(17) = 41.20394008
    yy(18) = 0.5699317760
    yy(19) = 0.4306056376
    yy(20) = 7.9906200783D-03
    yy(21) = 0.9056036089
    yy(22) = 1.6054258216D-02
    yy(23) = 0.7509759687
    yy(24) = 8.8582855955D-02
    yy(25) = 48.27726193
    yy(26) = 39.38459028
    yy(27) = 0.3755297257
    yy(28) = 107.7562698
    yy(29) = 29.77250546
    yy(30) = 88.32481135
    yy(31) = 23.03929507
    yy(32) = 62.85848794
    yy(33) = 5.546318688
    yy(34) = 11.92244772
    yy(35) = 5.555448243
    yy(36) = 0.9218489762
    yy(37) = 94.59927549
    yy(38) = 77.29698353
    yy(39) = 63.05263039
    yy(40) = 53.97970677
    yy(41) = 24.64355755
    yy(42) = 61.30192144
    yy(43) = 22.21000000
    yy(44) = 40.06374673
    yy(45) = 38.10034370
    yy(46) = 46.53415582
    yy(47) = 47.44573456
    yy(48) = 41.10581288
    yy(49) = 18.11349055
    yy(50) = 50.00000000

! Init ----------------------------------------------------
!   1. Common /pv/ xmv(12)
!   2. Common /teproc/
!       vcv(12), vst(12), ivst(12)
!
    Do i = 1, 12
        ! Init Common /pv/ xmv(12)
        xmv(i) = yy(i+38)
        ! Init Common /teproc/ vcv(12), vst(12), ivst(12)
        vcv(i) = xmv(i)
        vst(i) = 2.0D0
        ivst(i) = 0
    End Do

! Init Common /teproc/
!   vrng(12), vtr, vts, vtc, vtv,
!   htr(3), hwr, hws, sfr(8),
!   xst(8, 13), tst(13),
!   cpflmx, cpprmx,
!   vtau(12), xns(41),
!
    vrng(1) = 400.00
    vrng(2) = 400.00
    vrng(3) = 100.00
    vrng(4) = 1500.00
    vrng(7) = 1500.00
    vrng(8) = 1000.00
    vrng(9) = 0.03
    vrng(10) = 1000.
    vrng(11) = 1200.0

    vtr = 1300.0
    vts = 3500.0
    vtc = 156.5
    vtv = 5000.0

    htr(1) = 0.06899381054D0
    htr(2) = 0.05D0
    hwr = 7060.
    hws = 11138.

    sfr(1) = 0.99500
    sfr(2) = 0.99100
    sfr(3) = 0.99000
    sfr(4) = 0.91600
    sfr(5) = 0.93600
    sfr(6) = 0.93800
    sfr(7) = 5.80000D-02
    sfr(8) = 3.01000D-02

    xst(1, 1) = 0.0
    xst(2, 1) = 0.0001
    xst(3, 1) = 0.0
    xst(4, 1) = 0.9999
    xst(5, 1) = 0.0
    xst(6, 1) = 0.0
    xst(7, 1) = 0.0
    xst(8, 1) = 0.0
    tst(1) = 45.

    xst(1, 2) = 0.0
    xst(2, 2) = 0.0
    xst(3, 2) = 0.0
    xst(4, 2) = 0.0
    xst(5, 2) = 0.9999
    xst(6, 2) = 0.0001
    xst(7, 2) = 0.0
    xst(8, 2) = 0.0
    tst(2) = 45.

    xst(1, 3) = 0.9999
    xst(2, 3) = 0.0001
    xst(3, 3) = 0.0
    xst(4, 3) = 0.0
    xst(5, 3) = 0.0
    xst(6, 3) = 0.0
    xst(7, 3) = 0.0
    xst(8, 3) = 0.0
    tst(3) = 45.

    xst(1, 4) = 0.4850
    xst(2, 4) = 0.0050
    xst(3, 4) = 0.5100
    xst(4, 4) = 0.0
    xst(5, 4) = 0.0
    xst(6, 4) = 0.0
    xst(7, 4) = 0.0
    xst(8, 4) = 0.0
    tst(4) = 45.

    cpflmx = 280275.
    cpprmx = 1.3

    vtau(1) = 8.
    vtau(2) = 8.
    vtau(3) = 6.
    vtau(4) = 9.
    vtau(5) = 7.
    vtau(6) = 5.
    vtau(7) = 5.
    vtau(8) = 5.
    vtau(9) = 120.
    vtau(10) = 5.
    vtau(11) = 5.
    vtau(12) = 5.
    Do i = 1, 12
        vtau(i) = vtau(i)/3600.
    End Do

!
! Init Common /randsd/ g
    g = 1431655765.D0
! END Common /randsd/ g
!

    xns(1) = 0.0012D0
    xns(2) = 18.000D0
    xns(3) = 22.000D0
    xns(4) = 0.0500D0
    xns(5) = 0.2000D0
    xns(6) = 0.2100D0
    xns(7) = 0.3000D0
    xns(8) = 0.5000D0
    xns(9) = 0.0100D0
    xns(10) = 0.0017D0
    xns(11) = 0.0100D0
    xns(12) = 1.0000D0
    xns(13) = 0.3000D0
    xns(14) = 0.1250D0
    xns(15) = 1.0000D0
    xns(16) = 0.3000D0
    xns(17) = 0.1150D0
    xns(18) = 0.0100D0
    xns(19) = 1.1500D0
    xns(20) = 0.2000D0
    xns(21) = 0.0100D0
    xns(22) = 0.0100D0
    xns(23) = 0.250D0
    xns(24) = 0.100D0
    xns(25) = 0.250D0
    xns(26) = 0.100D0
    xns(27) = 0.250D0
    xns(28) = 0.025D0
    xns(29) = 0.250D0
    xns(30) = 0.100D0
    xns(31) = 0.250D0
    xns(32) = 0.100D0
    xns(33) = 0.250D0
    xns(34) = 0.025D0
    xns(35) = 0.050D0
    xns(36) = 0.050D0
    xns(37) = 0.010D0
    xns(38) = 0.010D0
    xns(39) = 0.010D0
    xns(40) = 0.500D0
    xns(41) = 0.500D0
! END Init Common /teproc/ ----------------------

! Init Common /dvec/ idv(20)
    Do i = 1, 20
        idv(i) = 0
    End Do
! END Init Common /dvec/ idv(20)

! Init Common /wlk/ --------------------------------------------------
!   hspan(12), hzero(12), sspan(12), szero(12), spspan(12),
!   tlast(12), tnext(12), adist(12), bdist(12), cdist(12), ddist(12),
!
! NOT Init: idvwlk(12),
!
    hspan(1) = 0.2D0
    hzero(1) = 0.5D0
    sspan(1) = 0.03D0
    szero(1) = 0.485D0
    spspan(1) = 0.D0

    hspan(2) = 0.7D0
    hzero(2) = 1.0D0
    sspan(2) = .003D0
    szero(2) = .005D0
    spspan(2) = 0.D0

    hspan(3) = 0.25D0
    hzero(3) = 0.5D0
    sspan(3) = 10.D0
    szero(3) = 45.D0
    spspan(3) = 0.D0

    hspan(4) = 0.7D0
    hzero(4) = 1.0D0
    sspan(4) = 10.D0
    szero(4) = 45.D0
    spspan(4) = 0.D0

    hspan(5) = 0.15D0
    hzero(5) = 0.25D0
    sspan(5) = 10.D0
    szero(5) = 35.D0
    spspan(5) = 0.D0

    hspan(6) = 0.15D0
    hzero(6) = 0.25D0
    sspan(6) = 10.D0
    szero(6) = 40.D0
    spspan(6) = 0.D0

    hspan(7) = 1.D0
    hzero(7) = 2.D0
    sspan(7) = 0.25D0
    szero(7) = 1.0D0
    spspan(7) = 0.D0

    hspan(8) = 1.D0
    hzero(8) = 2.D0
    sspan(8) = 0.25D0
    szero(8) = 1.0D0
    spspan(8) = 0.D0

    hspan(9) = 0.4D0
    hzero(9) = 0.5D0
    sspan(9) = 0.25D0
    szero(9) = 0.0D0
    spspan(9) = 0.D0

    hspan(10) = 1.5D0
    hzero(10) = 2.0D0
    sspan(10) = 0.0D0
    szero(10) = 0.0D0
    spspan(10) = 0.D0

    hspan(11) = 2.0D0
    hzero(11) = 3.0D0
    sspan(11) = 0.0D0
    szero(11) = 0.0D0
    spspan(11) = 0.D0

    hspan(12) = 1.5D0
    hzero(12) = 2.0D0
    sspan(12) = 0.0D0
    szero(12) = 0.0D0
    spspan(12) = 0.D0

    Do i = 1, 12
        tlast(i) = 0.D0
        tnext(i) = 0.1D0
        adist(i) = szero(i)
        bdist(i) = 0.D0
        cdist(i) = 0.D0
        ddist(i) = 0.D0
    End Do
! END Init Common /wlk/ -----------------------

    time = 0.0
    Call tefunc(nn, time, yy, yp)

    Return
End Subroutine teinit

!=============================================================================
!     TESUBi - Utility subroutines, i=1,2,..,8
!
!-----------------------------------------------------------------------
!> Update `h`.
!!
!!  ## Input
!!      @param[in]  `z`
!!      @param[in]  `t`
!!      @param[out] `h`
!!      @param[in]  `ity`
!!  ### common:
!!      + /const/
!!
Subroutine tesub1(z, t, h, ity)
    Double Precision avp, bvp, cvp, &
        ah, bh, ch, &
        ag, bg, cg,  av, &
        ad, bd, cd,  xmw
    Common /const/ avp(8), bvp(8), cvp(8), &
        ah(8), bh(8), ch(8), &
        ag(8), bg(8), cg(8),  av(8), &
        ad(8), bd(8), cd(8),  xmw(8)
    Double Precision z(8), r, t, h, hi
    Integer ity, i

    If (ity==0) Then
        h = 0.0D0
        Do i = 1, 8
            hi = t*(ah(i)+bh(i)*t/2.D0+ch(i)*t**2/3.D0)
            hi = 1.8D0*hi
            h = h + z(i)*xmw(i)*hi
        End Do
    Else
        h = 0.0D0
        Do i = 1, 8
            hi = t*(ag(i)+bg(i)*t/2.D0+cg(i)*t**2/3.D0)
            hi = 1.8D0*hi
            hi = hi + av(i)
            h = h + z(i)*xmw(i)*hi
        End Do
    End If

    If (ity==2) Then
        r = 3.57696D0/1.D6
        h = h - r*(t+273.15)
    End If

    Return
End Subroutine tesub1

!-----------------------------------------------------------------------
!> Update `t`.
!!
!!  ## Input
!!      @param[in]  `z`
!!      @param[inout] `t`
!!      @param[in]  `h`
!!      @param[in]  `ity`
!!  ### common:
!!      + /const/
!!
Subroutine tesub2(z, t, h, ity)
    Double Precision avp, bvp, cvp, ah, bh, ch, &
        ag, bg, cg, av, ad, bd, cd, xmw
    Common /const/ avp(8), bvp(8), cvp(8), ah(8), bh(8), ch(8), &
        ag(8), bg(8), cg(8), av(8), ad(8), bd(8), cd(8), xmw(8)
    Integer ity, j
    Double Precision z(8), t, h, tin, htest, err, dh, dt

    tin = t
    Do j = 1, 100
        Call tesub1(z, t, htest, ity)
        err = htest - h
        Call tesub3(z, t, dh, ity)
        dt = -err/dh
        t = t + dt
        If (dabs(dt)<1.D-12) Goto 300
    End Do
    t = tin

300 Return
End Subroutine tesub2

!-----------------------------------------------------------------------
!> Update `dh`.
!!
!!  ## Input
!!      @param[in] `z`
!!      @param[in] `t`
!!      @param[inout] `dh`
!!      @param[in] `ity`
!!  ### common:
!!      + /const/
!!
Subroutine tesub3(z, t, dh, ity)
    Double Precision avp, bvp, cvp, &
        ah, bh, ch, &
        ag, bg, cg,  av, &
        ad, bd, cd,  xmw
    Common /const/ avp(8), bvp(8), cvp(8), &
        ah(8), bh(8), ch(8), &
        ag(8), bg(8), cg(8),  av(8), &
        ad(8), bd(8), cd(8),  xmw(8)
    Integer ity, i
    Double Precision z(8), r, t, dh, dhi

    If (ity==0) Then
        dh = 0.0D0
        Do i = 1, 8
            dhi = ah(i) + bh(i)*t + ch(i)*t**2
            dhi = 1.8D0*dhi
            dh = dh + z(i)*xmw(i)*dhi
        End Do
    Else
        dh = 0.0D0
        Do i = 1, 8
            dhi = ag(i) + bg(i)*t + cg(i)*t**2
            dhi = 1.8D0*dhi
            dh = dh + z(i)*xmw(i)*dhi
        End Do
    End If

    If (ity==2) Then
        r = 3.57696D0/1.D6
        dh = dh - r
    End If

    Return
End Subroutine tesub3

!-----------------------------------------------------------------------
!> Update `r`.
!!
!!  ## Input
!!      @param[in] `x`
!!      @param[in] `t`
!!      @param[inout] `r`
!!  ### common:
!!      + /const/
!!
Subroutine tesub4(x, t, r)
    Double Precision avp, bvp, cvp, ah, bh, ch, ag, bg, cg, av, &
        ad, bd, cd,  xmw
    Common /const/ avp(8), bvp(8), cvp(8), ah(8), bh(8), ch(8), &
        ag(8), bg(8), cg(8),  av(8), &
        ad(8), bd(8), cd(8),  xmw(8)
    Double Precision v, r, x(8), t
    Integer i

    v = 0.0
    Do i = 1, 8
        v = v + x(i)*xmw(i)/(ad(i)+(bd(i)+cd(i)*t)*t)
    End Do
    r = 1.0/v

    Return
End Subroutine tesub4

!-----------------------------------------------------------------------
!> Update `[a-d]dist` and `tnext`.
!!
!!  ## Input
!!      @param[in] `s, sp`
!!      @param[in] `hspan, hzero`
!!      @param[in] `sspan, szero`
!!      @param[in] `spspan`
!!      @param[in] `idvflag`
!!      @param[in] `tlast`
!!      @param[out] `adist, bdist, cdist, ddist`
!!      @param[out] `tnext`
Subroutine tesub5(s, sp,  adist, bdist, cdist, ddist,   tlast, tnext, &
        hspan, hzero,   sspan, szero,   spspan, idvflag)
! ----------------------------------------------------------------------
    Double Precision s, sp, h, s1, s1p,   adist, bdist, cdist, ddist, &
        tlast, tnext,   hspan, hzero,   sspan, szero,   spspan, tesub7
    Integer i, idvflag

    i = -1
    h = hspan*tesub7(i) + hzero
    s1 = sspan*tesub7(i)*idvflag + szero
    s1p = spspan*tesub7(i)*idvflag

    adist = s
    bdist = sp
    cdist = (3.D0*(s1-s)-h*(s1p+2.D0*sp))/h**2
    ddist = (2.D0*(s-s1)+h*(s1p+sp))/h**3
    tnext = tlast + h

    Return
End Subroutine tesub5

!-----------------------------------------------------------------------
!> Update `x`.
!!
!!  ## Input
!!      @param[in] `std`
!!      @param[out] `x`
!!
Subroutine tesub6(std, x)
    Integer i
    Double Precision std, x, tesub7

    x = 0.D0
    Do i = 1, 12
        x = x + tesub7(i)
    End Do
    x = (x-6.D0)*std

    Return
End Subroutine tesub6

!-----------------------------------------------------------------------
!> Return random number.
!!
!!  ## Input
!!      @param[in] `i` flag
!!  ### commom
!!      + `/randsd/ g`
!!
Double Precision Function tesub7(i)
    Integer i
    Double Precision g, dmod
    Common /randsd/g

    g = dmod(g*9228907.D0, 4294967296.D0)
    If (i>=0) tesub7 = g/4294967296.D0
    If (i<0) tesub7 = 2.D0*g/4294967296.D0 - 1.D0

    Return
End Function tesub7

!-----------------------------------------------------------------------
!> Return dist?
!!
!!  ## Input
!!      @param[in] `i`
!!      @param[in] `t`
!!  ### commom
!!      + `/wlk/`
!!
Double Precision Function tesub8(i, t)
    Integer i
    Double Precision h, t
    Integer idvwlk
    Double Precision adist, bdist, cdist, ddist,  tlast, &
        tnext, hspan, hzero, sspan, szero, spspan
    Common /wlk/ adist(12), bdist(12), cdist(12), ddist(12), &
        tlast(12), tnext(12), hspan(12), hzero(12), sspan(12), &
        szero(12), spspan(12), idvwlk(12)

    h = t - tlast(i)
    tesub8 = adist(i) + h*(bdist(i)+h*(cdist(i)+h*ddist(i)))

    Return
End Function tesub8

max_stars_repo_name: inkydragon/Tennessee-Eastman-Challenge

max_issues_repo_name: inkydragon/Tennessee-Eastman-Challenge

max_forks_repo_name: inkydragon/Tennessee-Eastman-Challenge


text:
subroutine chkcom(lundia    ,error     ,neffil    ,ftype     ,gdp       )
!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation version 3.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: chkcom.f90 5986 2016-04-07 13:44:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/engines_gpl/flow2d3d/packages/kernel/src/inichk/chkcom.f90 $
!!--description-----------------------------------------------------------------
!
!    Function: Check the comm-file version by reading the element
!              FILE-VERSION or
!              check the existence of the group DWQTIM
! Method used:
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
    use precision
    use globaldata
    use string_module
    !
    implicit none
    !
    type(globdat),target :: gdp
    !
    ! The following list of pointer parameters is used to point inside the gdp structure
    !
!
! Local parameters
!
    integer, parameter :: start = 1, stopp = 2, incr = 3
!
! Global variables
!
    integer     , intent(in)  :: lundia !  Description and declaration in inout.igs
    logical     , intent(out) :: error  !  Flag=TRUE if an error is encountered
    character(*), intent(in)  :: neffil !  File name for FLOW NEFIS output files: tri"h/m/d"-"casl""labl" or for Comm. file com-"casl""labl"
    character(3), intent(in)  :: ftype  !  String containing to which output file version group or to diagnostic file should be written
!
! Local variables
!
    integer                        :: celidt                ! Array with inidices identifying the cell from where the reading and the writing starts
    integer                        :: datlen                ! Length of 'data file name'
    integer                        :: deflen                ! Length of 'definition file name'
    integer                        :: fds                   ! File description of NEFIS file set
    integer                        :: ierror               ! Local errorflag for NEFIS files
    integer                        :: ind
    integer, dimension(3)          :: uindex                !  Lowest allowed reference number
    integer, external              :: clsnef
    integer, external              :: crenef
    integer, external              :: getels
    integer, external              :: inqdat
    integer, external              :: neferr
    character(1)                   :: access                ! NEFIS file access type, here read only
    character(1)                   :: coding                ! NEFIS file representation, here NEUTRAL
    character(11)                  :: refnr                 ! Lowest allowed reference number
    character(16)                  :: elnms4                ! Element name defined for the NEFIS-files
    character(16)                  :: grdef5                ! Data-group name defined for the NEFIS-files as defined in definition part
    character(16)                  :: grnam4                ! Data-group name defined for the NEFIS-files as defined in data part
    character(16)                  :: grnam5                ! Data-group name defined for the NEFIS-files as defined in data part
    character(16), dimension(1)    :: hlptxt                ! Help character string to check comm file version number
    character(256)                 :: datnam                ! Data file name
    character(256)                 :: defnam                ! Definition file name
    character(256)                 :: errmsg                ! Character var. containing the errormessage to be written to file. The message depend on the error.
    character(256)                 :: filnam                ! Help var. for FLOW file name
    character(4)                   :: errnr                 ! Character var. containing the errormessage number corresponding to errormessage in ERRFIL
!
!
!! executable statements -------------------------------------------------------
!
    !
    !
    !
    !-----Initialize local variables
    !
    coding = 'N'
    error = .false.
    errnr = ' '
    fds   = -1
    refnr = '03.20.00.00'
    celidt = 1
    hlptxt(1) = refnr
    filnam = neffil
    !
    grnam4 = ftype // '-version'
    elnms4 = 'FILE-VERSION'
    !
    grnam5 = 'DWQTIM'
    !
    uindex(start) = celidt
    uindex(stopp) = celidt
    uindex(incr) = 1
    !
    !
    !-----aggregate file names
    !
    ind = index(filnam, ' ')
    !
    datnam = filnam
    datnam(ind:ind + 3) = '.dat'
    call remove_leading_spaces(datnam    ,datlen    )
    !
    defnam = filnam
    defnam(ind:ind + 3) = '.def'
    call remove_leading_spaces(defnam    ,deflen    )
    !
    !-----Open NEFIS files with read access
    !
    access = 'r'
    !
    ierror = crenef(fds, datnam(1:datlen), defnam(1:deflen), coding, access)
    if (ierror/=0) goto 9999
    !
    !     If NEFIS exist then it should be the right version
    !     group 'com-version', element 'FILE-VERSION'
    !
    ierror = getels(fds, grnam4, elnms4, uindex, 1, 16, hlptxt)
    if (ierror ==0) then
       if (hlptxt(1)(1:7)>=refnr(2:8)) then
       !----------file-version read is higher than lowest allowed number
       else
          !----------file-version read is lower than lowest allowed number
          error = .true.
          errnr = 'S027'
          errmsg = hlptxt(1)(1:11)
       endif
    else
       !
       !        group com-version does not exist
       !        check existence of flow-group DWQTIM
       !
       ierror = inqdat(fds, grnam5, grdef5)
       if (ierror ==0) then
          !           group DWQTIM exists
          error = .true.
          errnr = 'S027'
          errmsg = 'group not defined, but DWQTIM is defined'
       !           group DWQTIM does not exist
       else
            error = .true.
            errnr = 'U021'
            ierror = neferr(0, errmsg)
            call prterr(lundia    ,errnr     ,errmsg    )
            write (errmsg, '(5a)') 'Incomplete COM-files #',trim(datnam),'/',trim(defnam),'# exist'
            call prterr(lundia    ,errnr     ,errmsg    )
            errmsg = 'Try removing the COM-files first'
       endif
    endif
    !
    !
    !-----Close data and definition file
    !
    ierror = clsnef(fds)
    if (ierror/=0 .and. .not.error) then
       error = .true.
       errnr = 'U021'
       ierror = neferr(0, errmsg)
    endif
    !
    !-----Report error
    !
 9999 continue
    if (error) then
       call prterr(lundia    ,errnr     ,errmsg    )
    endif
end subroutine chkcom

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
      subroutine compute_th(ctfm,ctfmp,n_a1,n_alpha,n_b1,n_beta,np1_a1,n
     &p1_alpha,np1_b1,np1_beta,x,Nx,ht,hx,myzero,phys_bdy,res)
      implicit none
      integer i
      integer Nx
      real*8 ht
      real*8 hx
      real*8 myzero
      real*8 ctfm(Nx)
      real*8 ctfmp(Nx)
      real*8 n_a1(Nx)
      real*8 n_alpha(Nx)
      real*8 n_b1(Nx)
      real*8 n_beta(Nx)
      real*8 np1_a1(Nx)
      real*8 np1_alpha(Nx)
      real*8 np1_b1(Nx)
      real*8 np1_beta(Nx)
      real*8 x(Nx)
      integer phys_bdy(2)
      real*8 res(Nx)
      real*8 qb
      if (phys_bdy(1) .eq. 1) then
      do i=1, 1, 1
      qb = 0.1D1 + myzero * x(i)
      res(i)=qb
      end do
      endif
      do i=2, Nx-1, 1
      qb = ((0.5000000000000000D0 / np1_alpha(i) + 0.5000000000000000D0
     #/ n_alpha(i)) * (-0.1D1 * n_b1(i) + np1_b1(i)) / ht + 0.5000000000
     #000000D0 * (0.1D1 / np1_a1(i) - 0.1D1 * np1_beta(i) / np1_alpha(i)
     #) * (-0.5000000000000000D0 * (np1_b1(i - 1) - 0.1D1 * np1_b1(i + 1
     #)) / hx / ctfmp(i) + np1_b1(i) / ctfm(i)) + 0.5000000000000000D0 *
     # (0.1D1 / n_a1(i) - 0.1D1 * n_beta(i) / n_alpha(i)) * (-0.50000000
     #00000000D0 * (n_b1(i - 1) - 0.1D1 * n_b1(i + 1)) / hx / ctfmp(i) +
     # n_b1(i) / ctfm(i))) / (0.5000000000000000D0 * np1_b1(i) + 0.50000
     #00000000000D0 * n_b1(i))
      res(i)=qb
      end do
      if (phys_bdy(2) .eq. 1) then
      do i=Nx, Nx, 1
      qb = myzero * x(i)
      res(i)=qb
      end do
      endif
      END


text:
      Program Ex0204

      Implicit None

*   Initial, Final Point
      Double Precision  xMin
      Parameter         (xMin = 0.0D0)
      Double Precision  xMax
      Parameter         (xMax = 3.0D0)
      Double Precision  yMin
      Parameter         (yMin = 1.0D0)
*   Data Dimension
      Integer           DNum
      Parameter         (DNum = 9)
*   Data Storage
      Integer           GridN(DNum)
      Double Precision  Gap(DNum), y1RK1(DNum), y3RK1(DNum),yrRK1(DNum),
     $                  y1RK2(DNum), y3RK2(DNum), yrRK2(DNum),
     $                  y1RK3(DNum), y3RK3(DNum), yrRK3(DNum)
*   Dummy
      Integer           I
*   Function
      Double Precision  Func, Ext
      External          Func, Ext

      Open(11, File = '02-04.txt')
 1000 Format(F6.3,9F14.7)
 1010 Format(A6,9A14)
 1020 Format(A6,3A42)

*   Data Declaration
      Data GridN /6,15,30,60,150,300,600,1500,3000/
      Do 10 I = 1, DNum
        Gap(I) = (xMax - xMin) / DBLE(GridN(I))
  10  Continue

      Do 110 I = 1, DNum
        Call RK1(y1RK1(I), y3RK1(I), yrRK1(I), xMin,xMax,yMin, GridN(I),
     $          Gap(I))
        Call RK2(y1RK2(I), y3RK2(I), yrRK2(I), xMin,xMax,yMin, GridN(I),
     $          Gap(I))
        Call RK3(y1RK3(I), y3RK3(I), yrRK3(I), xMin,xMax,yMin, GridN(I),
     $          Gap(I))
  110 Continue

      Write(11,1020) '------',
     $               '------------------------------------------',
     $               '------------------------------------------',
     $               '------------------------------------------'
      Write(11,1020) ' |',
     $               'Runge-Kutta 2-order       |',
     $               'Runge-Kutta 3-order       |',
     $               'Runge-Kutta 4-order       |'
      Write(11,1020) ' |',
     $               'Eq.(2.22)            |',
     $               'Eq.(2.24)            |',
     $               'Eq.(2.25)            |'
      Write(11,1010) 'h |',
     $               'y(1)  ', 'y(3)  ', 'Retrun |',
     $               'y(1)  ', 'y(3)  ', 'Return |',
     $               'y(1)  ', 'y(3)  ', 'Retrun |'
      Write(11,1020) '------',
     $               '------------------------------------------',
     $               '------------------------------------------',
     $               '------------------------------------------'
      Do 20 I = 1, DNum
        Write(11,1000) Gap(I), y1RK1(I), y3RK1(I), yrRK1(I), y1RK2(I),
     $                 y3RK2(I), yrRK2(I), y1RK3(I), y3RK3(I), yrRK3(I)
  20  Continue
      Write(11,1020) '------',
     $               '------------------------------------------',
     $               '------------------------------------------',
     $               '------------------------------------------'

      End Program Ex0204

*-----------------------------------------------------------------------
*
      Subroutine RK1(y1, y3, yr, xMin, xMax, yMin, GridN, Gap)

      Implicit None

      Double Precision  y1, y3, yr, xMin, xMax, yMin, Gap
      Integer           GridN
*   Dummy
      Integer           I
*   Value for n+1, n+0
      Double Precision  xn1, yn1, xn0, yn0
*   Function
      Double Precision  Func, Ext
      External          Func, Ext
*   Internal Variable
      Double Precision  k1, k2

      xn0 = xMin
      yn0 = yMin
      Do 10 I = 1, GridN
        xn1 = DBLE(I) / DBLE(GridN) * xMax +
     $        (1.0D0 - DBLE(I) / DBLE(GridN)) * xMin
        k1 = Gap * Func(xn0, yn0)
        k2 = Gap * Func(xn0 + 0.5D0*Gap, yn0 + 0.5D0*k1)
        yn1 = yn0 + k2
        If(DAbs(xn1-1.0D0).LT.1.0D-10) y1 = yn1 - Ext(xn1)
        If(DAbs(xn1-3.0D0).LT.1.0D-10) y3 = yn1 - Ext(xn1)
        xn0 = xn1
        yn0 = yn1
  10  Continue

      Do 20 I = 1, GridN
        xn1 = DBLE(I) / DBLE(GridN) * xMin +
     $        (1.0D0 - DBLE(I) / DBLE(GridN)) * xMax
        k1 = -Gap * Func(xn0, yn0)
        k2 = -Gap * Func(xn0 - 0.5D0*Gap, yn0 + 0.5D0*k1)
        yn1 = yn0 + k2
        xn0 = xn1
        yn0 = yn1
  20  Continue
      yr = Ext(xn1) - yn1

      End Subroutine RK1

*-----------------------------------------------------------------------
*
      Subroutine RK2(y1, y3, yr, xMin, xMax, yMin, GridN, Gap)

      Implicit None

      Double Precision  y1, y3, yr, xMin, xMax, yMin, Gap
      Integer           GridN
*   Dummy
      Integer           I
*   Value for n+1, n+0
      Double Precision  xn1, yn1, xn0, yn0
*   Function
      Double Precision  Func, Ext
      External          Func, Ext
*   Internal Variable
      Double Precision  k1, k2, k3

      xn0 = xMin
      yn0 = yMin
      Do 10 I = 1, GridN
        xn1 = DBLE(I) / DBLE(GridN) * xMax +
     $        (1.0D0 - DBLE(I) / DBLE(GridN)) * xMin
        k1 = Gap * Func(xn0, yn0)
        k2 = Gap * Func(xn0 + 0.5D0*Gap, yn0 + 0.5D0*k1)
        k3 = Gap * Func(xn0 + Gap, yn0 - k1 + 2.0D0*k2)
        yn1 = yn0 + 1.0D0/6.0D0 * (k1 + 4.0D0*k2 + k3)
        If(DAbs(xn1-1.0D0).LT.1.0D-10) y1 = yn1 - Ext(xn1)
        If(DAbs(xn1-3.0D0).LT.1.0D-10) y3 = yn1 - Ext(xn1)
        xn0 = xn1
        yn0 = yn1
  10  Continue

      Do 20 I = 1, GridN
        xn1 = DBLE(I) / DBLE(GridN) * xMin +
     $        (1.0D0 - DBLE(I) / DBLE(GridN)) * xMax
        k1 = -Gap * Func(xn0, yn0)
        k2 = -Gap * Func(xn0 - 0.5D0*Gap, yn0 + 0.5D0*k1)
        k3 = -Gap * Func(xn0 - Gap, yn0 - k1 + 2.0D0*k2)
        yn1 = yn0 + 1.0D0/6.0D0 * (k1 + 4.0D0*k2 + k3)
        xn0 = xn1
        yn0 = yn1
  20  Continue
      yr = Ext(xn1) - yn1

      End Subroutine RK2

*-----------------------------------------------------------------------
*
      Subroutine RK3(y1, y3, yr, xMin, xMax, yMin, GridN, Gap)

      Implicit None

      Double Precision  y1, y3, yr, xMin, xMax, yMin, Gap
      Integer           GridN
*   Dummy
      Integer           I
*   Value for n+1, n+0
      Double Precision  xn1, yn1, xn0, yn0
*   Function
      Double Precision  Func, Ext
      External          Func, Ext
*   Internal Variable
      Double Precision  k1, k2, k3, k4

      xn0 = xMin
      yn0 = yMin
      Do 10 I = 1, GridN
        xn1 = DBLE(I) / DBLE(GridN) * xMax +
     $        (1.0D0 - DBLE(I) / DBLE(GridN)) * xMin
        k1 = Gap * Func(xn0, yn0)
        k2 = Gap * Func(xn0 + 0.5D0*Gap, yn0 + 0.5D0*k1)
        k3 = Gap * Func(xn0 + 0.5D0*Gap, yn0 + 0.5D0*k2)
        k4 = Gap * Func(xn0 + Gap, yn0 + k3)
        yn1 = yn0 + 1.0D0/6.0D0 * (k1 + 2.0D0*k2 + 2.0D0*k3 + k4)
        If(DAbs(xn1-1.0D0).LT.1.0D-10) y1 = yn1 - Ext(xn1)
        If(DAbs(xn1-3.0D0).LT.1.0D-10) y3 = yn1 - Ext(xn1)
        xn0 = xn1
        yn0 = yn1
  10  Continue

      Do 20 I = 1, GridN
        xn1 = DBLE(I) / DBLE(GridN) * xMin +
     $        (1.0D0 - DBLE(I) / DBLE(GridN)) * xMax
        k1 = -Gap * Func(xn0, yn0)
        k2 = -Gap * Func(xn0 - 0.5D0*Gap, yn0 + 0.5D0*k1)
        k3 = -Gap * Func(xn0 - 0.5D0*Gap, yn0 + 0.5D0*k2)
        k4 = -Gap * Func(xn0 - Gap, yn0 + k3)
        yn1 = yn0 + 1.0D0/6.0D0 * (k1 + 2.0D0*k2 + 2.0D0*k3 + k4)
        xn0 = xn1
        yn0 = yn1
  20  Continue
      yr = Ext(xn1) - yn1

      End Subroutine RK3

*-----------------------------------------------------------------------
*
      Function Func(x,y)

      Implicit None
      Double Precision  Func, x, y

      Func = - x * y

      End Function Func

*-----------------------------------------------------------------------
*
      Function FuncDx(x,y)

      Implicit None
      Double Precision  FuncDx, x, y

      FuncDx = - y

      End Function FuncDx

*-----------------------------------------------------------------------
*
      Function FuncDy(x,y)

      Implicit None
      Double Precision  FuncDy, x, y

      FuncDy = - x

      End Function FuncDy

*-----------------------------------------------------------------------
*
      Function FuncG(x)

      Implicit None
      Double Precision  FuncG, x

      FuncG = - x

      End Function FuncG

*-----------------------------------------------------------------------
*
      Function Ext(x)

      Implicit None
      Double Precision  Ext, x

      Ext = DExp(-x*x/2.0D0)

      End Function Ext

ext: for


text:
module dummy_assimilator

  implicit none

contains

  ! Copyright (c) 2004-2018 Lars Nerger
  !
  ! This file is part of PDAF.
  !
  ! PDAF is free software: you can redistribute it and/or modify
  ! it under the terms of the GNU Lesser General Public License
  ! as published by the Free Software Foundation, either version
  ! 3 of the License, or (at your option) any later version.
  !
  ! PDAF is distributed in the hope that it will be useful,
  ! but WITHOUT ANY WARRANTY; without even the implied warranty of
  ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  ! GNU Lesser General Public License for more details.
  !
  ! You should have received a copy of the GNU Lesser General Public
  ! License along with PDAF.  If not, see <http://www.gnu.org/licenses/>.
  !
  !$Id: PDAF-D_lenkf_analysis_rsm.F90 29 2018-03-09 20:06:46Z lnerger $
  !BOP

  ! !ROUTINE: PDAF_lenkf_analysis_rsm --- Perform LEnKF analysis step
  !
  ! !INTERFACE:
  SUBROUTINE dummy_assimilator_assimilate( &
    step, ind_p, dim_p, dim_obs_p, &
    dim_obs, dim_ens, rank_ana, state_p, ens_p, predictions, innovations, &
    U_add_obs_err, U_localize, forget, flag, info)

    ! !DESCRIPTION:
    ! Analysis step of ensemble Kalman filter with
    ! representer-type formulation.  In this version
    ! HP is explicitly computed.  This variant is
    ! optimal if the number of observations is
    ! smaller than or equal to half of the ensemble
    ! size.
    ! The final ensemble update uses a block
    ! formulation to reduce memory requirements.
    !
    ! Variant for domain decomposition.
    !
    ! !  This is a core routine of PDAF and
    !    should not be changed by the user   !
    !
    ! !REVISION HISTORY:
    ! 2003-10 - Lars Nerger - Initial code
    ! Later revisions - see svn log
    !
    ! !USES:
    ! Include definitions for real type of different precision
    ! (Defines BLAS/LAPACK routines and MPI_REALTYPE)

    USE iso_c_binding

    IMPLICIT NONE

    abstract INTERFACE
      SUBROUTINE add_obs_err(step, ind_p, dim_obs, HPH, info)
        ! Add observation error covariance matrix
        USE iso_c_binding
        INTEGER(c_int32_t), INTENT(in), value :: step, ind_p, dim_obs
        REAL(c_double), INTENT(inout) :: HPH(dim_obs, dim_obs)
        class(*), intent(in)::info
      END SUBROUTINE add_obs_err
      SUBROUTINE localize(step, ind_p, dim_p, dim_obs, HP_p, HPH, info)
        ! Apply localization to HP and HPH^T
        USE iso_c_binding
        INTEGER(c_int32_t), INTENT(in), value :: step, ind_p, dim_p, dim_obs
        REAL(c_double), INTENT(inout) :: HP_p(dim_obs, dim_p)
        REAL(c_double), INTENT(inout) :: HPH(dim_obs, dim_obs)
        class(*), intent(in)::info
      END SUBROUTINE localize

    END INTERFACE

    ! !ARGUMENTS:
    INTEGER(c_int32_t), INTENT(in), value :: step      ! Iteration number
    INTEGER(c_int32_t), INTENT(in), value :: ind_p      ! Integer index of PE
    INTEGER(c_int32_t), INTENT(in), value  :: dim_p     ! PE-local dimension of model state
    INTEGER(c_int32_t), INTENT(in), value :: dim_obs_p  ! PE-local dimension of observation vector
    INTEGER(c_int32_t), INTENT(in), value :: dim_obs    ! Global dimension of observation vector
    INTEGER(c_int32_t), INTENT(in), value :: dim_ens   ! Size of state ensemble
    INTEGER(c_int32_t), INTENT(in), value :: rank_ana  ! Rank to be considered for inversion of HPH
    REAL(c_double), INTENT(inout)  :: state_p(dim_p)        ! PE-local ensemble mean state
    REAL(c_double), INTENT(inout)  :: ens_p(dim_p, dim_ens) ! PE-local state ensemble
    REAL(c_double), INTENT(inout)  :: predictions(dim_obs, dim_ens) ! PE-local state ensemble
    REAL(c_double), INTENT(in)     :: innovations(dim_obs, dim_ens) ! Global array of innovations
    REAL(c_double), INTENT(in), value     :: forget    ! Forgetting factor
    INTEGER(c_int32_t), INTENT(out) :: flag    ! Status flag
    class(*), intent(inout)::info

    procedure(add_obs_err) :: U_add_obs_err
    procedure(localize) :: U_localize

  END SUBROUTINE dummy_assimilator_assimilate

end module dummy_assimilator

max_stars_repo_name: LightDA-assim/lightda-core

max_issues_repo_name: LightDA-assim/lightda-core

max_forks_repo_name: LightDA-assim/lightda-core


text:
!*******************************************************************
!                                                                  *
      SUBROUTINE MES(I)
!                                                                  *
!   Written by G. Gaigalas,                                        *
!   Vilnius,  Lithuania                             December 1993  *
!   Transform to fortran 90/95 by G. Gaigalas       December 2012  *
!   The last modification made by G. Gaigalas       October  2017  *
!                                                                  *
!*******************************************************************
!
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: I
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: J
      CHARACTER(LEN=10), DIMENSION(6) :: STRING5
!
      DATA STRING5/'   I T L S',' I T L S 2',' I T L S 3', &
                   'AW P 1 L S','WA P 1 L S','       W 1'/
!-----------------------------------------------
      IF(I > 50) THEN
        J=I-50
        WRITE(6,'(A)') ' error in func./sub. '
        WRITE(6,'(20X,A10)') STRING5(J)
        WRITE(6,'(A)') ' susimaise f sluoksnio termu kodavimas  '
      ELSE
        WRITE(6,'(A)') ' yra daugiau nei 2 ele. sluoks. f,g,h,i,k,l,m'
        WRITE(6,'(3X,I5)') I
        IF(I == 1) THEN
          WRITE(6,'(A)') ' error in Subroutine   W 1 G '
        ELSEIF(I == 2) THEN
          WRITE(6,'(A)') ' error in Subroutine   W 1 '
        ELSEIF(I == 11) THEN
          WRITE(6,'(A)') ' error in Function     I T L S  '
        ELSEIF(I == 12) THEN
          WRITE(6,'(A)') ' error in Function     I T L S 2  '
        ELSEIF(I == 13) THEN
          WRITE(6,'(A)') ' error in Function     I T L S 3  '
        ELSEIF(I == 30) THEN
          WRITE(6,'(A)') ' error in Subroutine   A 1 A 2 A 3 A 4 L S '
        ELSEIF(I == 31) THEN
          WRITE(6,'(A)') ' error in Subroutine   A 1 A 2 L S '
        ELSEIF(I == 32) THEN
          WRITE(6,'(A)') ' error in Subroutine   A 1 A 2 W 3 L S '
        ELSEIF(I == 33) THEN
          WRITE(6,'(A)') ' error in Subroutine   A 1 A W 2 L S '
        ELSEIF(I == 34) THEN
          WRITE(6,'(A)') ' error in Subroutine   W A 1 A 2 L S '
        ELSEIF(I == 35) THEN
          WRITE(6,'(A)') ' error in Subroutine   W 1 W 2 L S '
        ELSE
          WRITE(6,'(A)') ' error in unknown Subroutine  '
        ENDIF
      ENDIF
      WRITE(6,'(A)') ' Contact to   G. Gediminas please ! ! ! '
      STOP
      END SUBROUTINE MES

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum