text:
c> \file cauchy.f

c> \brief Compute the Generalized Cauchy Point along the projected gradient direction.
c>
c> For given x, l, u, g (with sbgnrm > 0), and a limited memory
c> BFGS matrix B defined in terms of matrices WY, WS, WT, and
c> scalars head, col, and theta, this subroutine computes the
c> generalized Cauchy point (GCP), defined as the first local
c> minimizer of the quadratic
c>
c>            Q(x + s) = g's + 1/2 s'Bs
c>
c> along the projected gradient direction P(x-tg,l,u).
c> The routine returns the GCP in xcp.
c>
c> @param n On entry n is the dimension of the problem.

c>          On exit n is unchanged.
c>
c> @param x On entry x is the starting point for the GCP computation.

c>          On exit x is unchanged.
c>
c> @param l On entry l is the lower bound of x.

c>          On exit l is unchanged.
c>
c> @param u On entry u is the upper bound of x.

c>          On exit u is unchanged.
c>
c> @param nbd On entry nbd represents the type of bounds imposed on the
c>               variables, and must be specified as follows:
c>               nbd(i)=

0 if x(i) is unbounded,

c>
1 if x(i) has only a lower bound,

c>
2 if x(i) has both lower and upper bounds, and

c>
3 if x(i) has only an upper bound.

c>            On exit nbd is unchanged.
c>
c> @param g On entry g is the gradient of f(x). g must be a nonzero vector.

c>          On exit g is unchanged.
c>
c> @param iorder iorder will be used to store the breakpoints in the piecewise
c>               linear path and free variables encountered.

c>               On exit,

iorder(1),...,iorder(nleft)
c>                               are indices of breakpoints
c>                               which have not been encountered;

c>
iorder(nleft+1),...,iorder(nbreak)
c>                               are indices of
c>                               encountered breakpoints; and

c>
iorder(nfree),...,iorder(n)
c>                               are indices of variables which
c>                               have no bound constraits along the search direction.

c>
c> @param iwhere On entry iwhere indicates only the permanently fixed (iwhere=3)
c>                  or free (iwhere= -1) components of x.

c>               On exit iwhere records the status of the current x variables.
c>                  iwhere(i)=

-3  if x(i) is free and has bounds, but is not moved

c>
 0   if x(i) is free and has bounds, and is moved

c>
 1   if x(i) is fixed at l(i), and l(i) .ne. u(i)

c>
 2   if x(i) is fixed at u(i), and u(i) .ne. l(i)

c>
 3   if x(i) is always fixed, i.e.,  u(i)=x(i)=l(i)

c>
-1  if x(i) is always free, i.e., it has no bounds.

c>
c> @param t working array; will be used to store the break points.
c>
c> @param d the Cauchy direction P(x-tg)-x
c>
c> @param xcp returns the GCP on exit
c>
c> @param m On entry m is the maximum number of variable metric corrections
c>             used to define the limited memory matrix.

c>          On exit m is unchanged.
c>
c> @param ws On entry this stores S, a set of s-vectors, that defines the
c>              limited memory BFGS matrix.

c>           On exit this array is unchanged.
c>
c> @param wy On entry this stores Y, a set of y-vectors, that defines the
c>              limited memory BFGS matrix.

c>           On exit this array is unchanged.
c>
c> @param sy On entry this stores S'Y, that defines the
c>              limited memory BFGS matrix.

c>           On exit this array is unchanged.
c>
c> @param wt On entry this stores the
c>              Cholesky factorization of (theta*S'S+LD^(-1)L'), that defines the
c>              limited memory BFGS matrix.

c>           On exit this array is unchanged.
c>
c> @param theta On entry theta is the scaling factor specifying B_0 = theta I.

c>              On exit theta is unchanged.
c>
c> @param col On entry col is the actual number of variable metric
c>               corrections stored so far.

c>            On exit col is unchanged.
c>
c> @param head On entry head is the location of the first s-vector (or y-vector)
c>                in S (or Y).

c>             On exit col is unchanged.
c>
c> @param p will be used to store the vector p = W^(T)d.
c>
c> @param c will be used to store the vector c = W^(T)(xcp-x).
c>
c> @param wbp will be used to store the row of W corresponding
c>              to a breakpoint.
c>
c> @param v working array
c>
c> @param nseg On exit nseg records the number of quadratic segments explored
c>                in searching for the GCP.
c>
c> @param iprint variable that must be set by the user.

c>       It controls the frequency and type of output generated:
c>

iprint<0    no output is generated;

c>
iprint=0    print only one line at the last iteration;

c>
0<iprint<99 print also f and |proj g| every iprint iterations;

c>
iprint=99   print details of every iteration except n-vectors;

c>
iprint=100  print also the changes of active set and final x;

c>
iprint>100  print details of every iteration including x and g;

c>       When iprint > 0, the file iterate.dat will be created to
c>                        summarize the iteration.
c>
c> @param sbgnrm On entry sbgnrm is the norm of the projected gradient at x.

c>               On exit sbgnrm is unchanged.
c>
c> @param info On entry info is 0.
c>             On exit info =

0       for normal return,

c>
= nonzero for abnormal return when the the system
c>                                     used in routine bmv is singular.

c> @param epsmch machine precision epsilon
      subroutine cauchy(n, x, l, u, nbd, g, iorder, iwhere, t, d, xcp,
     +                  m, wy, ws, sy, wt, theta, col, head, p, c, wbp,
     +                  v, nseg, iprint, sbgnrm, info, epsmch)
      implicit none
      integer          n, m, head, col, nseg, iprint, info,
     +                 nbd(n), iorder(n), iwhere(n)
      double precision theta, epsmch,
     +                 x(n), l(n), u(n), g(n), t(n), d(n), xcp(n),
     +                 wy(n, col), ws(n, col), sy(m, m),
     +                 wt(m, m), p(2*m), c(2*m), wbp(2*m), v(2*m)
c
c     References:
c
c       [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
c       memory algorithm for bound constrained optimization'',
c       SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.
c
c       [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN
c       Subroutines for Large Scale Bound Constrained Optimization''
c       Tech. Report, NAM-11, EECS Department, Northwestern University,
c       1994.
c
c       (Postscript files of these papers are available via anonymous
c        ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.)
c
c                           *  *  *
c
c     NEOS, November 1994. (Latest revision June 1996.)
c     Optimization Technology Center.
c     Argonne National Laboratory and Northwestern University.
c     Written by
c                        Ciyou Zhu
c     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
c
c
c     ************

      logical          xlower,xupper,bnded
      integer          i,j,col2,nfree,nbreak,pointr,
     +                 ibp,nleft,ibkmin,iter
      double precision f1,f2,dt,dtm,tsum,dibp,zibp,dibp2,bkmin,
     +                 tu,tl,wmc,wmp,wmw,ddot,tj,tj0,neggi,sbgnrm,
     +                 f2_org
      double precision one,zero
      parameter        (one=1.0d0,zero=0.0d0)

c     Check the status of the variables, reset iwhere(i) if necessary;
c       compute the Cauchy direction d and the breakpoints t; initialize
c       the derivative f1 and the vector p = W'd (for theta = 1).

      if (sbgnrm .le. zero) then
         if (iprint .ge. 0) write (6,*) 'Subgnorm = 0.  GCP = X.'
         call dcopy(n,x,1,xcp,1)
         return
      endif
      bnded = .true.
      nfree = n + 1
      nbreak = 0
      ibkmin = 0
      bkmin = zero
      col2 = 2*col
      f1 = zero
      if (iprint .ge. 99) write (6,3010)

c     We set p to zero and build it up as we determine d.

      do 20 i = 1, col2
         p(i) = zero
  20  continue

c     In the following loop we determine for each variable its bound
c        status and its breakpoint, and update p accordingly.
c        Smallest breakpoint is identified.

      do 50 i = 1, n
         neggi = -g(i)
         if (iwhere(i) .ne. 3 .and. iwhere(i) .ne. -1) then
c             if x(i) is not a constant and has bounds,
c             compute the difference between x(i) and its bounds.
            if (nbd(i) .le. 2) tl = x(i) - l(i)
            if (nbd(i) .ge. 2) tu = u(i) - x(i)

c           If a variable is close enough to a bound
c             we treat it as at bound.
            xlower = nbd(i) .le. 2 .and. tl .le. zero
            xupper = nbd(i) .ge. 2 .and. tu .le. zero

c              reset iwhere(i).
            iwhere(i) = 0
            if (xlower) then
               if (neggi .le. zero) iwhere(i) = 1
            else if (xupper) then
               if (neggi .ge. zero) iwhere(i) = 2
            else
               if (abs(neggi) .le. zero) iwhere(i) = -3
            endif
         endif
         pointr = head
         if (iwhere(i) .ne. 0 .and. iwhere(i) .ne. -1) then
            d(i) = zero
         else
            d(i) = neggi
            f1 = f1 - neggi*neggi
c             calculate p := p - W'e_i* (g_i).
            do 40 j = 1, col
               p(j) = p(j) +  wy(i,pointr)* neggi
               p(col + j) = p(col + j) + ws(i,pointr)*neggi
               pointr = mod(pointr,m) + 1
  40        continue
            if (nbd(i) .le. 2 .and. nbd(i) .ne. 0
     +                        .and. neggi .lt. zero) then
c                                 x(i) + d(i) is bounded; compute t(i).
               nbreak = nbreak + 1
               iorder(nbreak) = i
               t(nbreak) = tl/(-neggi)
               if (nbreak .eq. 1 .or. t(nbreak) .lt. bkmin) then
                  bkmin = t(nbreak)
                  ibkmin = nbreak
               endif
            else if (nbd(i) .ge. 2 .and. neggi .gt. zero) then
c                                 x(i) + d(i) is bounded; compute t(i).
               nbreak = nbreak + 1
               iorder(nbreak) = i
               t(nbreak) = tu/neggi
               if (nbreak .eq. 1 .or. t(nbreak) .lt. bkmin) then
                  bkmin = t(nbreak)
                  ibkmin = nbreak
               endif
            else
c                x(i) + d(i) is not bounded.
               nfree = nfree - 1
               iorder(nfree) = i
               if (abs(neggi) .gt. zero) bnded = .false.
            endif
         endif
  50  continue

c     The indices of the nonzero components of d are now stored
c       in iorder(1),...,iorder(nbreak) and iorder(nfree),...,iorder(n).
c       The smallest of the nbreak breakpoints is in t(ibkmin)=bkmin.

      if (theta .ne. one) then
c                   complete the initialization of p for theta not= one.
         call dscal(col,theta,p(col+1),1)
      endif

c     Initialize GCP xcp = x.

      call dcopy(n,x,1,xcp,1)

      if (nbreak .eq. 0 .and. nfree .eq. n + 1) then
c                  is a zero vector, return with the initial xcp as GCP.
         if (iprint .gt. 100) write (6,1010) (xcp(i), i = 1, n)
         return
      endif

c     Initialize c = W'(xcp - x) = 0.

      do 60 j = 1, col2
         c(j) = zero
  60  continue

c     Initialize derivative f2.

      f2 =  -theta*f1
      f2_org  =  f2
      if (col .gt. 0) then
         call bmv(m,sy,wt,col,p,v,info)
         if (info .ne. 0) return
         f2 = f2 - ddot(col2,v,1,p,1)
      endif
      dtm = -f1/f2
      tsum = zero
      nseg = 1
      if (iprint .ge. 99)
     +   write (6,*) 'There are ',nbreak,'  breakpoints '

c     If there are no breakpoints, locate the GCP and return.

      if (nbreak .eq. 0) goto 888

      nleft = nbreak
      iter = 1

      tj = zero

c------------------- the beginning of the loop -------------------------

 777  continue

c     Find the next smallest breakpoint;
c       compute dt = t(nleft) - t(nleft + 1).

      tj0 = tj
      if (iter .eq. 1) then
c         Since we already have the smallest breakpoint we need not do
c         heapsort yet. Often only one breakpoint is used and the
c         cost of heapsort is avoided.
         tj = bkmin
         ibp = iorder(ibkmin)
      else
         if (iter .eq. 2) then
c             Replace the already used smallest breakpoint with the
c             breakpoint numbered nbreak > nlast, before heapsort call.
            if (ibkmin .ne. nbreak) then
               t(ibkmin) = t(nbreak)
               iorder(ibkmin) = iorder(nbreak)
            endif
c        Update heap structure of breakpoints
c           (if iter=2, initialize heap).
         endif
         call hpsolb(nleft,t,iorder,iter-2)
         tj = t(nleft)
         ibp = iorder(nleft)
      endif

      dt = tj - tj0

      if (dt .ne. zero .and. iprint .ge. 100) then
         write (6,4011) nseg,f1,f2
         write (6,5010) dt
         write (6,6010) dtm
      endif

c     If a minimizer is within this interval, locate the GCP and return.

      if (dtm .lt. dt) goto 888

c     Otherwise fix one variable and
c       reset the corresponding component of d to zero.

      tsum = tsum + dt
      nleft = nleft - 1
      iter = iter + 1
      dibp = d(ibp)
      d(ibp) = zero
      if (dibp .gt. zero) then
         zibp = u(ibp) - x(ibp)
         xcp(ibp) = u(ibp)
         iwhere(ibp) = 2
      else
         zibp = l(ibp) - x(ibp)
         xcp(ibp) = l(ibp)
         iwhere(ibp) = 1
      endif
      if (iprint .ge. 100) write (6,*) 'Variable  ',ibp,'  is fixed.'
      if (nleft .eq. 0 .and. nbreak .eq. n) then
c                                             all n variables are fixed,
c                                                return with xcp as GCP.
         dtm = dt
         goto 999
      endif

c     Update the derivative information.

      nseg = nseg + 1
      dibp2 = dibp**2

c     Update f1 and f2.

c        temporarily set f1 and f2 for col=0.
      f1 = f1 + dt*f2 + dibp2 - theta*dibp*zibp
      f2 = f2 - theta*dibp2

      if (col .gt. 0) then
c                          update c = c + dt*p.
         call daxpy(col2,dt,p,1,c,1)

c           choose wbp,
c           the row of W corresponding to the breakpoint encountered.
         pointr = head
         do 70 j = 1,col
            wbp(j) = wy(ibp,pointr)
            wbp(col + j) = theta*ws(ibp,pointr)
            pointr = mod(pointr,m) + 1
  70     continue

c           compute (wbp)Mc, (wbp)Mp, and (wbp)M(wbp)'.
         call bmv(m,sy,wt,col,wbp,v,info)
         if (info .ne. 0) return
         wmc = ddot(col2,c,1,v,1)
         wmp = ddot(col2,p,1,v,1)
         wmw = ddot(col2,wbp,1,v,1)

c           update p = p - dibp*wbp.
         call daxpy(col2,-dibp,wbp,1,p,1)

c           complete updating f1 and f2 while col > 0.
         f1 = f1 + dibp*wmc
         f2 = f2 + 2.0d0*dibp*wmp - dibp2*wmw
      endif

      f2 = max(epsmch*f2_org,f2)
      if (nleft .gt. 0) then
         dtm = -f1/f2
         goto 777
c                 to repeat the loop for unsearched intervals.
      else if(bnded) then
         f1 = zero
         f2 = zero
         dtm = zero
      else
         dtm = -f1/f2
      endif

c------------------- the end of the loop -------------------------------

 888  continue
      if (iprint .ge. 99) then
         write (6,*)
         write (6,*) 'GCP found in this segment'
         write (6,4010) nseg,f1,f2
         write (6,6010) dtm
      endif
      if (dtm .le. zero) dtm = zero
      tsum = tsum + dtm

c     Move free variables (i.e., the ones w/o breakpoints) and
c       the variables whose breakpoints haven't been reached.

      call daxpy(n,tsum,d,1,xcp,1)

 999  continue

c     Update c = c + dtm*p = W'(x^c - x)
c       which will be used in computing r = Z'(B(x^c - x) + g).

      if (col .gt. 0) call daxpy(col2,dtm,p,1,c,1)
      if (iprint .gt. 100) write (6,1010) (xcp(i),i = 1,n)
      if (iprint .ge. 99) write (6,2010)

 1010 format ('Cauchy X =  ',/,(4x,1p,6(1x,d11.4)))
 2010 format (/,'---------------- exit CAUCHY----------------------',/)
 3010 format (/,'---------------- CAUCHY entered-------------------')
 4010 format ('Piece    ',i3,' --f1, f2 at start point ',1p,2(1x,d11.4))
 4011 format (/,'Piece    ',i3,' --f1, f2 at start point ',
     +        1p,2(1x,d11.4))
 5010 format ('Distance to the next break point =  ',1p,d11.4)
 6010 format ('Distance to the stationary point =  ',1p,d11.4)

      return

      end

max_stars_repo_name: jonathanschilling/L-BFGS-B

max_issues_repo_name: jonathanschilling/L-BFGS-B

max_forks_repo_name: jonathanschilling/L-BFGS-B


text:
\ Mforth compiler

\ Equates for terminal variables puts copies in header space. See tiff.c/tiff.h.
\ Some are omitted to save header space.

handler         equ handler
base            equ base                \ 6.1.0750
hp              equ hp
cp              equ cp
dp              equ dp
state           equ state               \ 6.1.2250
current         equ current
source-id       equ source-id           \ 6.2.2218
personality     equ personality
tibs            equ tibs
tibb            equ tibb
>in             equ >in                 \ 6.1.0560
c_wids          equ c_wids
c_casesens      equ c_casesens
\ head            equ head
context         equ context
forth-wordlist  equ forth-wordlist      \ 16.6.1.1595
hld             equ hld
blk             equ blk
tib             equ tib
\ maxtibsize    equ maxtibsize
|pad|           equ |pad|
pad             equ pad                 \ 6.2.2000
RAMsize         equ RAMsize
ROMsize         equ ROMsize
SPIflashBlocks  equ SPIflashBlocks

\ Define a compiler
\ char params:  c_colondef c_litpend c_slot c_called

defer FlushLit
defer NewGroup

: ClearIR   \ --                        \ Initialize the IR. slot=26, litpend=0
   0 iracc !  26 c_slot w!  0 c_called c!
;
: AppendIR  \ opcode imm --             \ insert opcode at slot, adding imm
   swap  c_slot c@ lshift  iracc @      \ shift it and get IR
   +  +                    iracc !      \ merge into instruction
;

: Implicit  \ opcode --                 \ add implicit opcode to group
   FlushLit
   c_slot c@ 0=  over 60 and  and if    \ doesn't fit in the last slot
      NewGroup
   then
   0 AppendIR
   c_slot c@  dup if
      6 -  dup 0< invert and            \ sequence: 26, 20, 14, 8, 2, 0
      c_slot c! exit
   then
   drop NewGroup                        \ slot 0 -> slot 26
;

: Explicit  \ imm opcode --
   swap FlushLit
   1  c_slot c@  dup >r  lshift         \ opcode imm maximum | slot
   over <  r> 0= or                     \ imm doesn't fit of slot=0
   if  NewGroup  then
   AppendIR
   c_slot c@  24 lshift
   cp @ +  calladdr !                   \ remember call address in case of ;
   0 c_slot c!  NewGroup
;

hex
: HardLit  \ n --                       \ compile a hard literal
   dup >r -if: negate |                 \ |n| | n
   1FFFFFF invert 2* over and  if       \ compile unsigned long
      drop r>  dup 18 rshift  op_lit Explicit   \ upper part
      FFFFFF and  op_litx Explicit      \ lower part
      exit
   then
   r> 0< if                             \ compile negative
      1- op_lit Explicit
      op_com Implicit  exit
   then
   op_lit Explicit                      \ compile positive
;

:noname \ FlushLit  \ --                \ compile a literal if it's pending
   c_litpend c@ if
      0 c_litpend c!
      nextlit @ HardLit
   then
   0 c_called c!                        \ everything clears the call tail
; is FlushLit

:noname \ NewGroup  \ --                \ finish the group and start a new one
   FlushLit                             \ if a literal is pending, compile it
   c_slot c@
   15 over - |-if 2drop exit |          \ already at first slot
   drop if
      op_no: Implicit                   \ skip unused slots
   then
   iracc @ ,c
   ClearIR
; is NewGroup

: literal,  \ n --                      \ compile a literal
   FlushLit
   1FFFFFF invert 2*  over and
   if HardLit exit then                 \ oversized literal, compile it now
   nextlit !  1 c_litpend c!            \ literal is pending
;
: compile,  \ xt --                     \ 6.2.0945  compile a call
   2/ 2/ op_call Explicit
   head @ 4 + c@  80 and  c_called c!   \ 0 = call-only
;
: ,exit  \ --
   c_called c@ if
      calladdr link>
      8  [ calladdr 3 + ] literal c@  lshift
      invert swap ROM!                  \ convert CALL to JUMP
      0 c_called c!
      0 calladdr ! exit
   then op_exit Implicit
;
decimal

ext: f

max_stars_repo_name: bradleyeckert/Mforth

max_issues_repo_name: bradleyeckert/Mforth

max_forks_repo_name: bradleyeckert/Mforth


!============================================================================
      SUBROUTINE grid_construction_TPM()
!============================================================================
! Region two-phase (TP)
!
! The specific internal energy along from critical point to triple
! point. The region is bounded by the saturation line LH and the triple
! pressure line, the critical internal, point E energy.
!
!
! The bilinear mapping is used for the interpolation from physical (v,e)
! doamin to the transformed domaine (X,Y).(Bicubic can be
! considered later by using also the derivatives)
!========================================================================
!
     USE def_constants, ONLY: pr,ord_spline,e_tri_R,e_cr,x_mesh_max,x_mesh_min,&
&                             e_tri_R,e_tri_L,v_tri_R,v_tri_L,P_tri,T_tri,c_cr,MMM_TPM,NNN_TPM,NNN_sat_TPM,&
&                             P_cr,T_cr,rho_cr
     USE def_variables, ONLY: x_mesh_TPM,y_mesh_TPM,y_mesh_sat_TPM,v_right_TPM,v_left_TPM,&
&        spline_right_TPM,spline_left_TPM,saturP,ppp_TPM,vvv_TPM,xxx_TPM,ccc_TPM,TTT_TPM
!&        vL_psat_spline, vV_psat_spline,uL_psat_spline, uV_psat_spline, Tsat_psat_spline
     USE non_linear_solvers,ONLY: New_Rap3D, BrentRoots
     USE properties, ONLY: pressure,sound_speed,satderiv,satprop
     USE grid_functions,ONLY: polyfit
!
      IMPLICIT NONE
!
      REAL(pr) :: u_max, u_min , delta, deltap, vvv_TPM_log
      REAL(pr) :: ev,el,vv,vl,qual,sound,pp,ratio,press,temp,Tsat
      INTEGER  :: i,k,j_tpl, Niter, flag,j
      REAL(pr) :: duL_dp, duV_dp, dvL_dp, dvV_dp,du_dp_x,dv_dp_x
      REAL(pr) :: pguess,res, T_guess,v_six,in_2,out_2,v_guess,vliq,vliq_guess,pguess1,&
&                 Tguess,vvguess,vlguess, res1,dummy,upper, lower

!--------------------------------------------------------------------------
!
! Construct vector y_mesh_TPM(physical domain), x_mesh_TPM(tansformed
! domain)
!
!-------------------------------------------------------------------------
        u_max = e_tri_R
        u_min = e_cr
!
        delta = (u_max - u_min) / (NNN_TPM - 1)
        y_mesh_TPM = u_min + (/(i*delta,i=0,NNN_TPM-1)/)
!
        delta = (x_mesh_max - x_mesh_min) / (MMM_TPM -1)
        x_mesh_TPM = x_mesh_min + (/(i*delta, i=0,MMM_TPM-1)/)
!
        delta             = (u_max - u_min)/(NNN_sat_TPM-1)
        y_mesh_sat_TPM   = u_min + (/(i*delta, i=0,NNN_sat_TPM-1)/)
!

! Right boundary = iso triple pressure line
!print*, 'stat right'
        ev = e_tri_R
        el = e_tri_L
        vv = v_tri_R
        vl = v_tri_L
        ppp_TPM(:,MMM_TPM) = P_tri
        TTT_TPM(:,MMM_TPM) = T_tri
!
        DO k=1,NNN_TPM
!
           xxx_TPM(k,MMM_TPM) = (y_mesh_TPM(k) - el) / (ev - el)
           qual = (y_mesh_TPM(k) - el) / (ev - el)
           vvv_TPM(k,MMM_TPM) = qual*vv + (1_pr-qual)*vl
! Speed of sound
!
           CALL satderiv(3, ppp_TPM(k,MMM_TPM), duL_dp, duV_dp, dvL_dp, dvV_dp)
!
           ratio = ((ev - el)/(vv - vl))   ! (J/kg)/(m3/kg)
           du_dp_x = qual * duV_dp + (1_pr - qual) * duL_dp
           dv_dp_x = qual * dvV_dp + (1_pr - qual) * dvL_dp
!
           sound = SQRT((P_tri + ratio)/(du_dp_x - ratio * dv_dp_x)) * vvv_TPM(k,MMM_TPM)
!
           ccc_TPM(k,MMM_TPM) = sound
!        print*,k,'v', vvv_TPM(k,MMM_TPM),'x',qual,'sound',ccc_TPM(k,MMM_TPM),y_mesh_TPM(k)
        END DO
!stop
!        vvv_TPM(1,MMM_TPM) = v_D
!        xxx_TPM(1,MMM_TPM) = (e_cr - el) / (ev - el)
!
! Left boundary == a part of right boundary of LH
!print*, 'stat left'
        vvv_TPM(1,1) = 1_pr/rho_cr
        ppp_TPM(1,1) = P_cr
        TTT_TPM(1,1) = T_cr
        ccc_TPM(1,1) = c_cr
        xxx_TPM(:,1) = 1_pr
!
!        v_six = 2.2530239912767237e-3_pr
!        DO i = 2, 18
!        ! Calculation of the volumes in the linear part for 100*100 mesh
!          vvv_TPM(i,1) = (y_mesh_TPM(i) - e_cr) / ( y_mesh_TPM(19) - e_cr) * &
!     &    (v_six - 1d0 / rho_cr) + 1d0 / rho_cr
!
!          CALL New_Rap1D(1, TTT_TPM(i,1), out_2, res, Niter, &
!     &         flag, y_mesh_TPM(i), T_cr-0.01,&
!     &         vvv_TPM(i,1), out_2)
!
!          CALL pressure(TTT_TPM(i,1),vvv_TPM(i,1),press)
!          ppp_TPM(i,1) = press
!!
!          CALL sound_speed(TTT_TPM(i,1),vvv_TPM(i,1), sound)
!          ccc_TPM(i,1) = sound
!
!          IF (res > 1e-12) THEN
!             print*, "Left in TPM", res, "iter", Niter,i
!             STOP
!          ENDIF
!
!  print*, i, vvv_TPM(i,1), ppp_TPM(i,1),TTT_TPM(i,1),ccc_TPM(i,1)
!        ENDDO
             T_guess = T_cr-0.001_pr
             v_guess = 1_pr / rho_cr + 1e-4_pr
             vliq_guess = 1_pr / rho_cr - 1e-4_pr
!             T_guess = TTT_TPM(18,1)
!             v_guess = vvv_TPM(18,1)
!             vliq_guess = 1_pr / rho_cr - 2e-4_pr
        DO i=2, NNN_TPM
!
         CALL New_Rap3D(4,vliq,vvv_TPM(i,1),TTT_TPM(i,1),&
     &        res, Niter, flag, y_mesh_TPM(i),in_2,&
     &        vliq_guess,v_guess,T_guess)
!
         v_guess = vvv_TPM(i,1)
         vliq_guess = vliq
         T_guess = TTT_TPM(i,1)
!
         CALL pressure(TTT_TPM(i,1),vvv_TPM(i,1),press)
         ppp_TPM(i,1) = press
!
         CALL sound_speed(TTT_TPM(i,1),vvv_TPM(i,1), sound)
         ccc_TPM(i,1) = sound
!
         IF (res > 1e-7) THEN
            print*, "Left in TPM", res, "iter", Niter
            STOP
         ENDIF
!      print*, i, vvv_TPM(i,1), ppp_TPM(i,1),TTT_TPM(i,1),ccc_TPM(i,1)
        END DO
!STOP
!
!
! Middle points
!
!
!print*, 'stat middle'
!         DO j=2,MMM_TPM-1
          DO i=1, NNN_TPM
!             pguess = ppp_TPM(i,1)
!          DO i=1, NNN_TPM
             DO j=2,MMM_TPM-1
!
!           vvv_TPM(i,j) = ((x_mesh_TPM(j) - x_mesh_min) /&
!     &    (x_mesh_max - x_mesh_min)) *(vvv_TPM(i,MMM_TPM)-vvv_TPM(i,1))&
!     &     +vvv_TPM(i,1)
           vvv_TPM_log = ((x_mesh_TPM(j) - x_mesh_min) /&
     &     (x_mesh_max - x_mesh_min)) *(log10(vvv_TPM(i,MMM_TPM))-log10(vvv_TPM(i,1)))+&
     &     log10(vvv_TPM(i,1))
           vvv_TPM(i,j) = 10_pr**vvv_TPM_log
!
!      print*, vvv_TPL(i,MMM_TPL), vvv_TPL(i,1)
!      print*, y_mesh_TPL(i),vvv_TPL(i,j),pguess
!           CALL New_Rap1D(2, press, qual, res1, Niter, flag,&
!     &           y_mesh_TPM(i), pguess, vvv_TPM(i,j), temp)
!
           lower = P_tri
           upper = P_cr - 0.01e6_pr
!           upper = pguess+2.0e6_pr
!print*, i,j,'lower',lower,'upper',upper
           CALL BrentRoots(1, press, temp, res, Niter, y_mesh_TPM(i), lower, upper, vvv_TPM(i,j), dummy)
!
!print*, 'v',vvv_TPM(i,j),'e',y_mesh_TPM(i), 'p,T',press,temp
!print*, ' '
           ppp_TPM(i,j) = press
           TTT_TPM(i,j) = temp
!           pguess       = press !- 1e4_pr
!
           IF (res  > 1e-10_pr) THEN
             print*, "middle points TPM", res, "iter", Niter,i,j
             STOP
           ENDIF
!
! Speed of sound
!
           CALL satprop(3, press, dummy, vv, vl, ev, el)
           CALL satderiv(3, press, duL_dp, duV_dp, dvL_dp, dvV_dp)
!
           qual = (vvv_TPM(i,j)-vl) / (vv-vl)
!
           ratio = ((ev - el)/(vv - vl))   ! (J/kg)/(m3/kg)
           du_dp_x = qual * duV_dp + (1_pr - qual) * duL_dp
           dv_dp_x = qual * dvV_dp + (1_pr - qual) * dvL_dp
!
           sound = SQRT((press + ratio)/(du_dp_x - ratio * dv_dp_x)) * vvv_TPM(i,j)
!
           ccc_TPM(i,j) = sound
           xxx_TPM(i,j) = qual
!
!      print*,j,i, ppp_TPM(i,j), TTT_TPM(i,j), xxx_TPM(i,j)
          END DO
         END DO
!STOP
!  print*, 'TPM','1,1', y_mesh_TPM(1),  vvv_TPM(1,1),&
!&                          TTT_TPM   (1,1),ppp_TPM(1,1),ccc_TPM(1,1)
!  print*, 'TPM','2,1   ', y_mesh_TPM(2),  vvv_TPM(2,1),&
!&                          TTT_TPM   (2,1),ppp_TPM(2,1),ccc_TPM(2,1)
!  print*, 'TPM','1,2   ', y_mesh_TPM(1),  vvv_TPM(1,2),&
!&                          TTT_TPM   (1,2),ppp_TPM(1,2),ccc_TPM(1,2)
!
! Spline line coefficient right boundary
!
!print*, 'stat spline right'
        ev = e_tri_R
        el = e_tri_L
        vv = v_tri_R
        vl = v_tri_L
       DO i=1,NNN_sat_TPM
           qual = (y_mesh_sat_TPM(i) - el) / (ev - el)
           v_right_TPM(i) = qual*vv + (1_pr-qual)*vl
       ENDDO
       DO i = 1, NNN_TPM-1
          j = ord_spline * (i-1) + 1
          CALL polyfit(spline_right_TPM(:,i), &
     &               y_mesh_sat_TPM(j:j+ord_spline),&
     &               v_right_TPM(j:j+ord_spline), ord_spline)
       ENDDO
!
! Spline line coefficient left boundary
!
!print*, 'stat spline left'
        v_left_TPM(1) = 1_pr/rho_cr

!        T_guess = T_cr - 0.0001_pr
!        v_guess =  1_pr/rho_cr + 1e-5_pr
!        vliq_guess = 1_pr/rho_cr - 1e-5_pr
!
! USE test_grid to compute T_guess, v_guess, vliq_guess
        T_guess = 304.1281
        v_guess =  2.1443612972628459E-003
        vliq_guess = 2.1322415401290890E-003

        DO i=2, NNN_sat_TPM
!              T_guess_sat = T_sat_sat(i-1)
!              v_guess_sat = v_Lsat_LH(i-1)
!              vliq_guess_sat = vliq_sat_LH(i-1)
!           IF( i == 2) THEN
!print*,i, 'in', y_mesh_sat_TPM(i),in_2,vliq_guess,v_guess, T_guess

          CALL New_Rap3D(4,vliq,v_left_TPM(i),Tsat,&
      &        res, Niter, flag, y_mesh_sat_TPM(i),in_2,&
      &        vliq_guess,v_guess, T_guess)
!
!print*, 'out', vliq,v_left_TPM(i),Tsat,res,Niter,flag

          CALL New_Rap3D(4,vliq,v_left_TPM(i),Tsat,&
      &        res,Niter,flag, y_mesh_sat_TPM(i),in_2,&
      &        vliq,v_left_TPM(i),Tsat)
!
!           ELSE
!
!               CALL New_Rap3D(4,vliq_sat_LH(i),v_Lsat_LH(i),T_sat_sat(i),&
!              & res_sat_LH(i), Niter, exitflag, y_mesh_sat_LH(i),in_2,&
!              & vliq_guess_sat,v_guess_sat, T_guess_sat)
!           ENDIF

              T_guess    = Tsat
              v_guess    = v_left_TPM(i)
              vliq_guess = vliq
!
              IF (res > 1e-8) THEN
               print*, "Left spline TPM", res, "i", i,&
               &        'T=',Tsat
               STOP
              ENDIF
!      print*, i, v_left_TPM(i)
        END DO
!
!
      DO i = 1, NNN_TPM-1
        j = 3 * (i-1) + 1
        CALL polyfit(spline_right_TPM(:,i), y_mesh_sat_TPM(j:j+3),v_right_TPM(j:j+3),ord_spline)
       CALL polyfit(spline_left_TPM(:,i), y_mesh_sat_TPM(j:j+3),v_left_TPM(j:j+3),ord_spline)
!        CALL polyfit(spline_left_TPM(:,i), y_mesh_sat_TPM(j:j+3),LOG10(v_left_TPM(j:j+3)),ord_spline)
      ENDDO
!
!     print*,'------------------------------------------------------'
     print*,           'TPM GRID FINISH'
!     print*,'------------------------------------------------------'
!
!
      END SUBROUTINE grid_construction_TPM