Документ взят из кэша поисковой машины. Адрес оригинального документа : http://hea.iki.rssi.ru/conf/hea2007/presentations/after_dinner_speech/4-Disk/16C.FFF Дата изменения: Sun Mar 18 17:28:28 2001 Дата индексирования: Tue Oct 2 02:04:52 2012 Кодировка:

c**********************************************************
c Main program for sphere
c**********************************************************
c Set of dimensions :
c Number of concentrations : (0 - 4)
c Number of materics : (0 - 16)
c**********************************************************
common /blc00/ni,nj,nk,nc,nm,nt,di,dj,dk,dt,tt,ro
* /blc01/npm(30),wpm(10),eem(12),cnd(20)
* /blmvi/vi(&mv) /blmvj/vj(&mv) /blmvk/vk(&mv)
* /blmui/ui(&mv) /blmuj/uj(&mv) /blmuk/uk(&mv)
* /blmwi/wi(&mv) /blmwj/wj(&mv) /blmwk/wk(&mv)
* /blmqi/qi(&mv) /blmqj/qj(&mv) /blmqk/qk(&mv)
* /blmth/th(&mt) /blmph/ph(&mp)
* /blmsv/sv(&ms) /blmdg/dg(&md)
* /blmda/da(&md) /blmdb/db(&md) /blmdc/dc(&md)
* /blmdd/dd(&md) /blmde/de(&md) /blmdf/df(&md)
* /blism/ism(&ni) /bljsm/jsm(&nj) /blksm/ksm(&nk)
* /blsim/sim(&li) /blsjm/sjm(&lj) /blskm/skm(&lk)
common /blname/nfaa,nfbb,namea,nameb,namec,named
character*12 namea,nameb,namec,named,namef
c****************'data file '****************************
data namef/'16c .dat '/
data nvers/5/ ! - Version # 5 from Aug-01-1999
c**********************************************************
nt=nvers; call fildsph(&ni,&nj,&nk,&nc,&nm,namef)
10 call cdirsph ! - Looking for result files
if(nt.eq.0.and.tt.ne.0.) goto 10 ! - Read new file
20 call stepsph ! - Go to the next time step
call heatsph ! - Pz, Th, Ci ( t~=t+dt/2 )
call mtrasph ! - Continents ( t~=t+dt/2 )
30 call globsph ! - Exiliary vector (Ui,Uj,Uk)
if(nt.eq.0) goto 10 ! - Read new file
if(npm(20).lt.0) goto 40 ! - Global iteration
call ellpsph ! - Pressure, correction (Ui,Uj,Uk)
call ellisph ! - Additional vector component (Wi)
call elljsph ! - Additional vector component (Wj)
call ellksph ! - Additional vector component (Wk)
call optpsph ! - Calculation optimum parameters
if(npm(20).gt.0) goto 30 ! - Global iteration
40 call heatsph ! - Pz, Th, Ci ( t=t~+dt/2 )
call mtrasph ! - Continents ( t=t~+dt/2 )
call soutsph(1) ! - Write & output
if(nt.eq.0) goto 10 ! - Read new file
if(npm(30).gt.1) goto 20 ! - Next step
stop; end

function funcrat(ff)
c************************************************************
c Calculation of the right hand function
c ff(15) - dimension of parameters :
c (ff(1),ff(2),ff(3)) = (r,q,f) - point coordinats;
c ff(4) - Ph (pressure); ff(5) - Pz (phase);
c ff(6) - Th (themperature); ff(12) - t (time);
c ff(7) : ff(10) - Ca : Cd (concentrations A,B,C and D);
c ff(11) - Sij*Sij (Tensor); ff(13) - Vr (radial veloc.)
c ff(14) - Reserv; ff(15) - Reserv
c------------------------------------------------------------
c Rat, Rac - Railegh numbers for Temperature, Concentration
c************************************************************
common /blc00/ni,nj,nk,nc,nm,nt,di,dj,dk,dt,tt,ro
common /blc01/npm(30),wpm(10),eem(12),cnd(20)
dimension ff(15); Rat=cnd(1); Rac=cnd(2)
c-------------------------------------------------------------
r=ff(1); q=ff(2); f=ff(3) ! Coordinats of point
Ph=ff(4); Pz=ff(5) ! Pressure and Phase
Th=ff(6); Ss=ff(11) ! Temperature and Tensor inv.
Ca=ff(7); Cb=ff(8) ! Concentrations A and B
Cc=ff(9); Cd=ff(10) ! Concentrations C and D
t=ff(12); Vr=ff(13) ! Time and Radial Velocity
To=functho(r) ! Initial Temperature distribution
c-------------------------------------------------------------
c funcrat=Rat*Th-Rac*Pz+1.e+5*Ca+2.e+6*Cb
funcrat=Rat*Th-Rac*Pz
return; end

function funqqq(ff,ic)
c************************************************************
c Internal sourcees ic=0 - for Temperature
c ic=1,2,3,4 - for Concentrations
c************************************************************
c ff(15) - dimension of parameters :
c-------------------------------------------------------------
dimension ff(15)
r=ff(1); q=ff(2); f=ff(3) ! Coordinats of point
Ph=ff(4); Pz=ff(5) ! Pressure and Phase
Th=ff(6); Ss=ff(11) ! Temperature and Tensor inv.
Ca=ff(7); Cb=ff(8) ! Concentrations A and B
Cc=ff(9); Cd=ff(10) ! Concentrations C and D
t=ff(12); Vr=ff(13) ! Time and Radial Velocity
c-------------------------------------------------------------
funqqq=0. ! QT=F(r,q,f,Th,Ci,Vr)
c if(ic.eq.0) funqqq=5.
c if(ic.eq.0) funqqq=5.0*ss
c if(ic.eq.0) funqqq=Th*Vr
c------- Internal sourcees for Concentrations ----------------
if(ic.eq.1) funqqq=(1.-exp(-t)) ! For Conc.A
if(ic.eq.2) funqqq=0. ! For Conc.B
if(ic.eq.3) funqqq=0. ! For Conc.C
if(ic.eq.4) funqqq=0. ! For Conc.D
return; end

function funcsv(ff)
c************************************************************
c Viscosity coefficient
c************************************************************
c ff(15) - dimension of parameters :
c------------------------------------------------------------
dimension ff(15)
r=ff(1); q=ff(2); f=ff(3) ! Coordinats of point
Ph=ff(4); Pz=ff(5) ! Pressure and Phase
Th=ff(6); Ss=ff(11) ! Temperature and Tensor inv.
Ca=ff(7); Cb=ff(8) ! Concentrations A and B
Cc=ff(9); Cd=ff(10) ! Concentrations C and D
t=ff(12); Vr=ff(13) ! Time and Radial Velocity
c-------------------------------------------------------------
c funcsv=exp(2.0+1.5*(1.-r))-exp(2.0*(Th+0.1))
c XXX=(1.0-r)/0.454
c aa=1.0+20.*(1.0-r)
c bb=1.0*(Th+0.15+XXX)
c cc=aa/bb
c funcsv=0.01*exp(cc)
funcsv=1.0
ic=ifhmtr(r,q,f)
if(ic.eq.1) funcsv=10.
if(ic.eq.2) funcsv=10.
if(ic.eq.3) funcsv=10.
if(ic.eq.4) funcsv=10.
if(ic.eq.5) funcsv=10.
if(ic.eq.6) funcsv=10.
if(ic.eq.7) funcsv=10.
if(ic.eq.8) funcsv=10.
if(ic.eq.9) funcsv=10.
if(ic.eq.10) funcsv=10.
if(ic.eq.11) funcsv=10.
if(ic.eq.12) funcsv=10.
if(ic.eq.13) funcsv=10.
if(ic.eq.14) funcsv=10.
if(ic.eq.15) funcsv=10.
if(ic.eq.16) funcsv=10.
return; end

function funcsq(ff,ic)
c************************************************************
c Termoconductivity coefficient (if ic=0)
c Diffusion coeff. for concentr. (if ic=1,2,3,4)
c************************************************************
c ff(15) - dimension of parameters :
c-------------------------------------------------------------
dimension ff(15)
r=ff(1); q=ff(2); f=ff(3) ! Coordinats of point
Ph=ff(4); Pz=ff(5) ! Pressure and Phase
Th=ff(6); Ss=ff(11) ! Temperature and Tensor inv.
Ca=ff(7); Cb=ff(8) ! Concentrations A and B
Cc=ff(9); Cd=ff(10) ! Concentrations C and D
t=ff(12); Vr=ff(13) ! Time and Radial Velocity
c-------------------------------------------------------------
if(ic.eq.0) funcsq=(1.+0.*Th)
if(ic.ge.1) funcsq=1.
return; end

function funcph(r,Tx,fx)
c************************************************************
c Calculation of phase transition
c r - radius of the point, Tx - temperature.
c Function returns cp; if ep.lt.0 - no transition.
c------------------------------------------------------------
c cp(r)=cp(0)*(1.+ca*Tx*(F-F**2))
c F=0.5*(1.+thn((r-rp+at*(Tx-Tp))/ep))
c G=1./(F-F**2)
c------------------------------------------------------------
common /blc01/npm(30),wpm(10),eem(12),cnd(20)
data cp/1.0/ ! base heat capacity (1.00)
data ca/0.00/ ! changihg heat capacity (0.01)
data at/-0.03/ ! temperature coefficient (-0.003)
data rp/0.8965/ ! phase transition lavel (0.8965)
data Tp/0.468/ ! phase transition temperature (0.4??)
c if Tp<0 - calculate Tp=To(rp)
data ep/0.005/ ! phase transition thickness (0.0016)
c if ep<0 - no transition
c************************************************************
funcph=cp; eem(11)=ep; eem(12)=rp; if(ep.le.0.) return
if(Tp.lt.0.) Tp=functho(rp); x=(r-rp+at*(Tx-Tp))/ep
fx=x; x=abs(x); if(x.le.10.) then; s=exp(x)
G=(s+1./s)**2; funcph=cp*(1.+ca*Tx/G); endif
return; end

function functh(r,q,f)
c************************************************************
c Initial and boundary themperature distribution
c (r,q,f) - coordinats of the point
c************************************************************
common /blc00/ni,nj,nk,nc,nm,nt,di,dj,dk,dt,tt,ro
x=(r-ro)/(1.-ro); y=q+f; To=functho(r); Tq=0.
if(x.gt.0.1.and.x.lt.0.9) Tq=0.5*To*sin(q)
functh=To+Tq*cos(f)
return; end

function functho(r)
c************************************************************
c Static themperature distribution; r - radius
c Ta - down Temperature; Tb - up Temperature
c************************************************************
common /blc00/ni,nj,nk,nc,nm,nt,di,dj,dk,dt,tt,ro
Ta=1.; Tb=0.; x=(r-ro)/(r-ro*r)
functho=Ta+(Tb-Ta)*x
return; end

function funcch(r,q,f,ic)
c************************************************************
c Initial and boundary concentration distribution
c (r,q,f) - coordinats of the point
c ic - number of concentration (ic=1,2,3,4)
c if(ic=0) calls functh(r,q,f) - temperature distr.
c************************************************************
common /blc00/ni,nj,nk,nc,nm,nt,di,dj,dk,dt,tt,ro
if(ic.eq.0) then; funcch=functh(r,f,q); return; endif
x=(r-ro)/(1.-ro); sq=sin(q); sf=sin(f); funcch=0.
if(ic.eq.1) then; f=1.+x**2 ! Concentr A
if(x.le.0.001) f=f+tt ! Down Boundary
if(x.gt.0.999) f=f+tt**2 ! Upper Boundary
endif
if(ic.eq.2) f=1.-x ! Concentr B
funcch=f
return; end

function funcvs(rx,ic)
c************************************************************
c Sedimantation velocity, rx - radius
c ic=0 - for Temperature; ic=1,2,3,4 - for Concentrations
c************************************************************
r=rx; i=ic
funcvs=0.
return; end

function functs(Th,rx)
c************************************************************
c Tg=functs - Graphics Temperature (adiabatic)
c Th - real (superadiabatic) Temperature ; rx - radius
c************************************************************
c common /blc00/ni,nj,nk,nc,nm,nt,di,dj,dk,dt,tt,ro
c dimension txx(11); data nnn/11/
c data txx/0.,0.1,0.2,0.4,0.4,0.4,
c * 0.4,0.4,0.6,0.7,0.8/
c x=(1.-rx)/(1.-ro); ss=abs(x*(nnn-1)-0.0001)
c i=1+int(ss); s=(ss+1.)-i
c Tg=((1.-s)*txx(i)+s*txx(i+1))
c functs=0.*Th+Tg
c return; end

common /blc00/ni,nj,nk,nc,nm,nt,di,dj,dk,dt,tt,ro
c pi=3.14159;
x=(1.-rx)/(1.-ro); Tx=Th
c Tg=-0.5*sin(0.5*pi*x)
Tg=0.43*x
functs=Tg
return; end

function funcptn(ff)
c************************************************************
c Potential function for curent point
c ff(15) - dimension of parameters :
c (ff(1),ff(2),ff(3)) = (r,q,f) - point coordinats;
c ff(4) - Ph (pressure); ff(5) - Pz (phase);
c ff(6) - Th (themperature); ff(12) - t (time);
c ff(7) : ff(10) - Ca : Cd (concentrations A,B,C and D);
c------------------------------------------------------------
c************************************************************
common /blc00/ni,nj,nk,nc,nm,nt,di,dj,dk,dt,tt,ro
dimension ff(15); Th=ff(6); rx=ff(1)
funcptn=(1.-Th)*(1.+0.*rx)
return; end