subroutine therm_init(ngtabrem,tempref,cpref,href,cgref,wt,r0,nsi, 1 temp1tabl,odtabl,temptabl,cgsptabl, 2 hsptabl,cpsptabl,cpsptb_tl,hsptb_tl, 3 cgsptb_tl) IMPLICIT NONE include 'parameter.par' include 'therm.com' integer ngtabrem,nsi double precision tempref,cpref,href,cgref,r0,tempv real*8 wt(nsi) character*80 fmt integer ncol,j,k double precision temp1tabl double precision odtabl real*8 temptabl(ntmx) real*8 cgsptabl(nsmx,ntmx) real*8 hsptabl(nsmx,ntmx) real*8 cpsptabl(nsmx,ntmx) real*8 cpsptb_tl(nsmx,ntmx) real*8 hsptb_tl(nsmx,ntmx) real*8 cgsptb_tl(nsmx,ntmx) c------------------------------------------------------------------------------- c ... fill up ns for the therm common blocks, and to be used below n_species = nsi c------------------------------------------------------------------------------- c ... read species specific heats cp_k(T) table open(unit=50,file='cpsp.tab',status='old') read(50,*) read(50,'(2x,e22.14)')dtab read(50,'(2x,2i9)')ngtab,ncol nstab = ncol-1 ngtabrem = ngtab if(ngtab .gt. ntmx)stop 'ngtab > ntmx' if(nstab .gt. nsmx)stop 'nstab > nsmx' dtab = dtab/tempref write(fmt,'(a,i3,a)') '(',nstab+1,'e22.14)' do j=1,ngtab read(50,fmt)temptab(j),(cpsptab(k,j),k=1,nstab) temptab(j) = temptab(j)/tempref temptabl(j) = temptab(j) do k=1,nstab cpsptab(k,j) = cpsptab(k,j)/cpref cpsptabl(k,j) = cpsptab(k,j) end do end do write(*,40) nstab, ngtab, dtab 40 format('cpsp_k table has ',i3,' species, ',i5, + ' points, interval',1p,e10.3) if(nstab.ne.n_species)stop 1 'therm_init: cpsp_k: nstab_vs_ns mismatch' close(50) c------------------------------------------------------------------------------- odtab = 1./dtab temp1tab = temptab(1) odtabl = odtab temp1tabl = temp1tab c write(*,*)' 1/dtab = odtab = ',odtab c write(*,*)' temptab(1) = temp1tab = ',temp1tab c------------------------------------------------------------------------------- c ... read species enthalpies h_k(T) table open(unit=50,file='hk.tab',status='old') read(50,*) read(50,'(2x,e22.14)')dtab read(50,'(2x,2i9)')ngtab,ncol nstab = ncol-1 if(ngtab .gt. ntmx)stop 'ngtab > ntmx' if(nstab .gt. nsmx)stop 'nstab > nsmx' dtab = dtab/tempref write(fmt,'(a,i3,a)') '(',nstab+1,'e22.14)' do j=1,ngtab read(50,fmt)tempv,(hsptab(k,j),k=1,nstab) do k=1,nstab hsptab(k,j) = hsptab(k,j)/href hsptabl(k,j) = hsptab(k,j) end do end do write(*,50) nstab, ngtab, dtab 50 format('hsp_k table has ',i3,' species, ',i5, + ' points, interval',1p,e10.3) if(nstab.ne.n_species)stop 'therm_init: hsp_k: nstab_vs_ns mismatch' close(50) c------------------------------------------------------------------------------- c ... read species gibbs free energies g_k(T) = h_k(T)-T*s_k(T) table c ... then fill up the table with g_k*wt(k) or G_k, (per mole rather c ... than per kg), for computational convenience c ... so: cgsptab(j,k) is G(j,k)=g(j,k)*wt(k) non-dimensionalized by c ... cgref=gref*Wref=r0*Tref open(unit=50,file='gk.tab',status='old') read(50,*) read(50,'(2x,e22.14)')dtab read(50,'(2x,2i9)')ngtab,ncol nstab = ncol-1 if(ngtab .gt. ntmx)stop 'ngtab > ntmx' if(nstab .gt. nsmx)stop 'nstab > nsmx' dtab = dtab/tempref write(fmt,'(a,i3,a)') '(',nstab+1,'e22.14)' do j=1,ngtab read(50,fmt)tempv,(cgsptab(k,j),k=1,nstab) do k=1,nstab cgsptab(k,j) = cgsptab(k,j)*wt(k)/cgref cgsptabl(k,j) = cgsptab(k,j) end do end do write(*,60) nstab, ngtab, dtab 60 format('cgsp_k table has ',i3,' species, ',i5, + ' points, interval',1p,e10.3) if(nstab.ne.n_species)stop 'therm_init: cgsp_k: nstab_vs_ns mismatch' close(50) c------------------------------------------------------------------------------- c ... read d(species_specific_heats)/dT: cp_k_T(T) table open(unit=50,file='cpsp_T.tab',status='old') read(50,*) read(50,'(2x,e22.14)')dtab read(50,'(2x,2i9)')ngtab,ncol nstab = ncol-1 if(ngtab .gt. ntmx)stop 'ngtab > ntmx' if(nstab .gt. nsmx)stop 'nstab > nsmx' dtab = dtab/tempref write(fmt,'(a,i3,a)') '(',nstab+1,'e22.14)' do j=1,ngtab read(50,fmt)tempv,(cpsptb_t(k,j),k=1,nstab) do k=1,nstab cpsptb_t(k,j) = cpsptb_t(k,j)*(tempref/cpref) cpsptb_tl(k,j) = cpsptb_t(k,j) end do end do write(*,70) nstab, ngtab, dtab 70 format('cpsp_k_T table has ',i3,' species, ',i5, + ' points, interval',1p,e10.3) if(nstab.ne.n_species)stop 1 'therm_init: cpsp_k_T:nstab_vs_ns mismatch' close(50) c------------------------------------------------------------------------------- c ... read species enthalpies h_k_T(T) table open(unit=50,file='hk_T.tab',status='old') read(50,*) read(50,'(2x,e22.14)')dtab read(50,'(2x,2i9)')ngtab,ncol nstab = ncol-1 if(ngtab .gt. ntmx)stop 'ngtab > ntmx' if(nstab .gt. nsmx)stop 'nstab > nsmx' dtab = dtab/tempref write(fmt,'(a,i3,a)') '(',nstab+1,'e22.14)' do j=1,ngtab read(50,fmt)tempv,(hsptb_t(k,j),k=1,nstab) do k=1,nstab hsptb_t(k,j) = hsptb_t(k,j)*(tempref/href) hsptb_tl(k,j) = hsptb_t(k,j) end do end do write(*,80) nstab, ngtab, dtab 80 format('hsp_k_T table has ',i3,' species, ',i5, + ' points, interval',1p,e10.3) if(nstab.ne.n_species)stop 1 'therm_init: hsp_k_T: nstab_vs_ns mismatch' close(50) c------------------------------------------------------------------------------- c ... read species gibbs free energies g_k_T(T) =d(h_k(T)-T*s_k(T))/dT table c ... then fill up the table with g_k_T*wt(k) or G_k_T, (per mole rather c ... than per kg), for computational convenience c ... so: cgsptb_T(j,k) is G_T(j,k)=d(g(j,k)*wt(k))/dT non-dimensionalized by c ... cgref/Tref=gref*Wref/Tref=r0*Tref/Tref = r0 open(unit=50,file='gk_T.tab',status='old') read(50,*) read(50,'(2x,e22.14)')dtab read(50,'(2x,2i9)')ngtab,ncol nstab = ncol-1 if(ngtab .gt. ntmx)stop 'ngtab > ntmx' if(nstab .gt. nsmx)stop 'nstab > nsmx' dtab = dtab/tempref write(fmt,'(a,i3,a)') '(',nstab+1,'e22.14)' do j=1,ngtab read(50,fmt)tempv,(cgsptb_t(k,j),k=1,nstab) do k=1,nstab cgsptb_t(k,j) = cgsptb_t(k,j)*wt(k)/r0 cgsptb_tl(k,j) = cgsptb_t(k,j) end do end do write(*,90) nstab, ngtab, dtab 90 format('cgsp_k_T table has ',i3,' species, ',i5, + ' points, interval',1p,e10.3) if(nstab.ne.n_species)stop 1 'therm_init: cgsp_k_T:nstab_vs_ns mismatch' close(50) c------------------------------------------------------------------------------- return end c******************************************************************************