After looking at it… it is a small code that calls somethings I do not understand:

```
!gfortran Perplex_LitMod.for SUB* -o Perplex_LitMod -mpc80 -m64 -DD -O2
implicit none
include 'perplex_parameters.h'
integer i,idead,ii,j,k,nodes_inv,whocalled,loop1,useless,iki,
* icounting
logical meemum, nodata, bulk
character amount*6 !, yes*1
integer itri(4),jtri(4),ijpt,ierflag
double precision ctot, wt(3)
double precision props,psys,psys1
common/ cxt22 /props(i8,k5),psys(i8),psys1(i8)
integer kkp, np, ncpd, ntot
double precision cp3, ctot3
common/ cxt15 /cp3(k5,k5),ctot3(k5),kkp(k5),np,ncpd,ntot
integer iwt
common/ cst209 /iwt
double precision atwt
common/ cst45 /atwt(k0)
double precision v,tr,pr,r,ps
common/ cst5 /v(l2),tr,pr,r,ps
integer ipot,jv,iv
common / cst24 /ipot,jv(l2),iv(l2)
character*8 vname,xname
common/ csta2 /xname(k5),vname(l2)
cccccccccccccccccccccccccccccccccccccccccc
character pname*14
common/ cxt21a /pname(k5)
character*2 trim_name(k5)
double precision pcomp
common/ cst324 /pcomp(k5,k5)
ccccccccccccccccccccccccccccccccccccccccccc
character*5 cname
common/ csta4 /cname(k5)
integer jbulk
double precision cblk
common/ cst300 /cblk(k5),jbulk
double precision a,b
common/ cst313 /a(k5,k1),b(k5)
integer icomp,istct,iphct,icp
common/ cst6 /icomp,istct,iphct,icp
integer io3,io4,io9
common / cst41 /io3,io4,io9
c meemum flag used for output routines
c to distunguish meemum from werami
save meemum
data meemum/.true./
! REAL(8),ALLOCATABLE:: vs_sys(:),den_sys(:),vp_sys(:),
! * compres_p(:),alpha_p(:)
REAL(8) vs_sys,den_sys,vp_sys,
* compres_p,alpha_p
REAL(8) t_inv,p_inv, t_eq
!! ...intialize thermodynamic solver (Perple_X solver)
CALL iniprp
! Do
! WRITE(*,*)'Input T (C), P (bar)?, Al2O3,FeO'
READ(*,*) t_inv,p_inv,cblk(2),cblk(3)
t_inv=t_inv+273.15D0
! WRITE(*,*) 'Input bulk composition (wt%) SiO2,Al2O3,FeO,MgO,CaO '
! READ(*,*) cblk(1:5)
! WRITE(*,*) 'Input bulk composition (wt%) Al2O3,FeO'
! READ(*,*) cblk(2),cblk(3)
cblk(4)=49.369-4.106*cblk(2)+.343*cblk(2)**2D0
cblk(5)= -.164+.906*cblk(2)
cblk(1)=100D0-SUM(cblk(2:5))
bulk = .true.
icounting=0
! cblk=0.d0
ierflag=0
amount = 'weight'
idead=0
!T in K, P in bar
if(t_inv<773.1) then
t_eq=773.1
else
t_eq=t_inv
endif
v(jv(1))=t_eq
v(jv(2))=p_inv
IF (bulk) THEN
if (iwt.eq.1) then
! Convert mass to molar
do i = 1, jbulk
cblk(i) = cblk(i)/atwt(i)
enddo
endif
! normalize the composition vector, this
!is necessary for reasons of stupidity (lpopt0).
ctot = 0.d0
do i = 1, icp
ctot = ctot + cblk(i)
enddo
do i = 1, icp
b(i) = cblk(i)/ctot
enddo
! lpopt does the minimization and outputs
! the results to the print file.
call lpopt0 (idead)
IF (idead>0) THEN
write (*,*) 'Energy Minimization failed at this point'
do iki=1,jbulk
write(*,*) 'node cmp T(K) P(bar) '
write(*,*) ii, cblk(iki)*atwt(iki),v(jv(1)),v(jv(2))
enddo
write (*,*) 'I will ignore this point and continue'
icounting=1+icounting ;idead=0;ierflag=1 ;nodes_inv=useless
ELSE
! compute derivative properties
! v(jv(1))=t_inv
call getloc (itri,jtri,ijpt,wt,nodata,meemum)
! den_sys(nodes_inv)=psys(10) ! kg m-3
! vs_sys(nodes_inv)=psys(8) ! km s-1
! vp_sys(nodes_inv)=psys(7) ! km s-1
! compres_p(nodes_inv)=psys(14) / 1.0D5 ! (1/Pa)
! alpha_p(nodes_inv)=psys(13) ! (1/K)
den_sys=psys(10) ! kg m-3
vs_sys=psys(8) ! km s-1
vp_sys=psys(7) ! km s-1
compres_p=psys(14) / 1.0D5 ! (1/Pa)
alpha_p=psys(13)
ENDIF
ENDIF
! WRITE(*,*)'Density (kg/m3)', den_sys
! WRITE(*,*)'Vp (km/s)', vp_sys
! WRITE(*,*)'Vs (km/s)', vs_sys
WRITE(*,*) vp_sys, vs_sys, den_sys
! WRITE(*,*)'ntot',ntot
! DO i=1,ntot
! WRITE(*,*) 'pname(i)',pname(i)
! WRITE(*,*) ' props(17,i)*props(16,i)/psys(17)*1d2',
! * props(17,i)*props(16,i)/psys(17)*1d2
! write(*,*)'props para la fase',props(:,i)
! write(*,*) 'pcomp', pcomp(:,i)
! ENDDO
! enddo
END
```