cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c This is a program for numerical simulations of bubble growth c induced by diffusion of volatiles from the melt of finite volume. c Non-isothermal case. c c The physical model includes: c 1) limited volume of liquid around the bubble and spherical symmetry; c 2) diffusion of volatiles through the bubble walls; diffusion c coefficient is temperature and composition dependent c 3) varying saturation of the volitiles at the bubble interface c according to pressure inside it and Henry's law; c 4) mass balance between dissolved volatiles and gas inside the bubble; c 5) gas state law inside the bubble; c 6) hydrodynamics of liquid around growing bubble with continuity and c momentum equations for viscous Newtonian fluid; visocity coefficient c is temperature and composition dependent c 7) surface tension; c 8) heat of volatile evaporation; c 9) work of gas expansion; c 10) heat of viscous friction (dissipation energy); c 10) heat ballance at the bubble interface (energy conservation); c 11) thermal conductivity through the bubble wall. c c To solve the system of differential equations the following methods c are used: c 1) Lagrangian system of coordinates; c 2) non-dimensionalization; c 3) normalization. c All of this enabled to fix the computational domain for the c concentration profile. In this system they are not moving as at the c bubble interface as well as at the outer surface of the melt drop c (elementary cell with a single bubble inside). Then we used running c function of coordinate parameter to make integration gridding very c fine at the boundary without increasing number of intervals c (beta function). For computing,the concentration is also substituted c by potential function, which provided to eleminate the convection term c in diffusion equation. c c c Dr. Alexander Proussevitch alex.proussevitch@unh.edu c Prof. Dork L. Sahagian dork.sahagian@unh.edu c c Complex Systems Research Center, office (603)862-4796 c Institute for the Study of fax (603)862-0188 c Earth, Oceans, and Space, c University of New Hampshire, c Morse Hall, Room 445, c Durham, NH 03824-3525 c c Last updated: March 20, 1997 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c m1 number of nodal intervals to calculate the concentration c c0 initial gas concentration in the melt c cm heat capacity of the melt (P=const) c cp heat capacity of the gas (P=const) c denmelt melt density c depth0 initial depth of magma rise c dif dimensionless diffusion coefficient c dif0 initial diffusion coefficient c dpower diffusion power coefficient c dt dimensionless time step c dTvc vitrification or crystallization temperature interval c enact activation energy of viscous flow c eta dimensionless viscosity of the melt c eta0 initial viscosity of the melt c hen Henry's constant c heat dimensionless heat of volatile evaporation c heat0 initial heat of volatile evaporation c heatvc latent heat of vitrification or crystallization c pa ambient pressure c ph initial ambient pressure c r dimensionless bubble radius c r0 initial bubble radius in meter c rf final equilibrium bubble radius c rg universal gas constant c s0 initial sphere radius in meter c sigma surface tension of the gas-melt interface c td dimensionless time (R0*R0/D) c t total dimensional time c t1 total dimensionless time c tg dimensionless temperature of the gas c temp0 initial temperature of the system c tempvc liquidus temperature of vitrification or crystallization c vcoef coefficient of volatile effect on viscosity c vh rate of magma rise c wm molecular weight of the gas c xi temperature conductivity coefficient c ypar(i) dimensionless parameters of the process cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c implicit double precision(a-h,o-z) parameter(m1=500) parameter(kkf=10) character*1 yes character*14 fname,fpress,fprof,ftemp character*70 info dimension info(5) dimension y1(m1+1),y2(m1+1),yu(m1+1) dimension b(m1+1),y2t(m1+1),y2c(m1+1) dimension co(m1+1),cn(m1+1),con(m1+1),cont(m1+1) dimension dydyn1(m1+1),dydyn2(m1+1),dydyn3(m1+1),dydyn4(m1+1) dimension ypar(10) dimension aa(m1+1),bb(m1+1),cc(m1+1),rr(m1+1) dimension dif(m1+1),dtv(m1+1) dimension to(m1+1),tn(m1+1),temp(m1+1),tmp(m1+1) common enact,rg,vcoef,tempvc,dTvc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c Open output files for radius, pressures, concentration c profiles, and temperature terms and read information lines c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc open(unit=21,file='nbubble.dat',status='old') do 8 i=1,100 read(21,'(a1)') yes if(yes.ne.'#')goto 9 8 continue 9 backspace 21 read (21,'(a,t12)') fname fpress=fname fprof=fname ftemp=fname iend=index(fname,' ') fname(iend:iend+1)='.d' fpress(iend:iend+4)='-ps.d' fprof(iend:iend+4)='-pf.d' ftemp(iend:iend+4)='-tm.d' open(unit=77,file=fname,status='unknown') open(unit=78,file=fpress,status='unknown') open(unit=79,file=fprof,status='unknown') open(unit=80,file=ftemp,status='unknown') ccccccccccccccccccccccccccccc jinfo=0 do 10 i=1,20 read(21,'(a14)') fpress if(fpress(1:8).eq.'***start')then do 11 j=1,5 read (21,'(a70)') info(j) if(info(j)(1:6).eq.'***end')goto 12 jinfo=j 11 continue endif 10 continue cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Input of constants and writing formats cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 12 c0=dataread('c0') cm=dataread('cm') cp=dataread('cp') delta=0.0d+00 dcoef=dataread('dcoef') denact=dataread('denact') denmelt=dataread('denmelt') depth0=dataread('depth0') dTvc=dataread('dTvc') enact=dataread('enact') grav=9.8d+00 heat0=dataread('heat0') heatvc=dataread('heatvc')/dTvc hen=dataread('hen') pa=dataread('pa')*1.0d+06 r0=dataread('r0') rg=8.3143d+00 s0=dataread('s0') sigma=dataread('sigma') temp0=dataread('temp0')+2.7315d+02 tempvc=dataread('tempvc')+2.7315d+02 time=dataread('time') vcoef=dataread('vcoef') vh=dataread('vh') if(vh.eq.0.0d+00.and.pa.eq.0.0d+00)pa=depth0*denmelt*grav+1.0d+05 depth=(pa-1.0d+05)/denmelt/grav if(vh.ne.0.0d+00)pa=1.0d+05 p0=denmelt*grav*depth0+pa if(vh.ne.0.0d+00.and.c0.eq.0.0d+00) * c0=dsqrt(hen*(p0+2.0d+00*sigma/r0)) wm=dataread('wm') xi=dataread('xi') if(heat0.eq.0.0d+00.and.wm.eq.18.0d-03)then heat0=heatf(pa,temp0) iheat=1 else iheat=0 endif imonitor=ifix(sngl(dataread('monitor'))) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 50 format(10x,a1,i2,a3,e13.6) 51 format(1x,f13.6,4x,a) 52 format(1x,e13.6,4x,a) 53 format(1x,a,t40,f13.6,3x,a) 54 format(1x,a,t40,e13.6,3x,a) 55 format(1x,f14.7,f11.3,f12.4,2x,f9.3,2f11.5,f7.3) 551 format(1x,f14.3,f11.3,f12.4,2x,f9.3,2f11.5,f7.3) 56 format(1x,a,i1,a) 57 format(f14.5,f10.2,f12.4,f10.3,2x,3f8.4) 58 format(1x,f10.4,f13.6,f14.4,e17.6) 59 format(1x,3f12.5) 60 format(1x,f14.7,f11.3,2f12.6,f11.5,e14.3) 601 format(1x,f14.3,f11.3,2f12.6,f11.5,e14.3) 61 format(/,a7,e9.5,/,a37) 611 format("Time step is too big for the bubble. - No root found") 62 format(1x,f14.7,f10.2,5f10.3) 621 format(1x,f14.3,f10.2,5f10.3) 63 format("Time step is too big for the bubble. - Too many steps") 64 format("Bubble shrinks to",f7.3," %") pi=3.1415927d+00 crt=1.0d+00/3.0d+00 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Input of parameters cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ypar(1)=dataread('ypar(1)') ypar(2)=dataread('ypar(2)') ypar(3)=dataread('ypar(3)') ypar(4)=dataread('ypar(4)') ypar(5)=dataread('ypar(5)') ypar(6)=dataread('ypar(6)') ypar(7)=dataread('ypar(7)') ypar(8)=dataread('ypar(8)') ypar(9)=dataread('ypar(9)') ypar(10)=dataread('ypar(10)') close(unit=21) cccccccccccccccccccccccccccccccccccc test=1.0d+00 do 20 i=1,10 test=test*ypar(i) 20 continue cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc eta0=vis(temp0,c0) dif0=c0*dexp(-dcoef-denact/rg/temp0) ph1=p0/pa if(vh.eq.0.0d+00)ph1=1.0d+00 vp=denmelt*grav*vh*r0*r0/dif0/pa cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Calculating dimensionless parameters cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if(test.eq.0.0d+00)then ypar(1)=r0**3/(s0**3-r0**3) ypar(2)=denmelt*rg*temp0/wm/pa ypar(3)=2.0d+00*sigma/r0/pa ypar(4)=4.0d+00*eta0*dif0/pa/r0/r0 ypar(5)=dsqrt(hen*pa) ypar(6)=xi/dif0 ypar(7)=1.2d+01*eta0*dif0/r0/r0/denmelt/cm/temp0 ypar(8)=cm/cp ypar(9)=heat0/cm/temp0/wm ypar(10)=heatvc/cm cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Calculating system conditions cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc else pa=denmelt*rg*temp0/wm/ypar(2) r0=dsqrt(4.0d+00*eta0*dif0/pa/ypar(4)) s0=r0*(1.0d+00+1.0d+00/ypar(1))**crt sigma=ypar(3)*r0*pa/2.0d+00 hen=ypar(5)**2/pa xi=ypar(6)*dif0 cm=1.2d+01*eta0*dif0/r0/r0/denmelt/ypar(7)/temp0 cp=cm/ypar(8) heat0=ypar(9)*cm*temp0*wm heatvc=ypar(10)*cm iheat=0 endif cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Verification of the conditions for bubble growth cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if(vh.ne.0.0d+00)delta=1.0d-06 ratio=(c0*c0/ypar(5)/ypar(5)-ph1)/ypar(3)+delta if(ratio.le.1.0d+00)then write(77,*)'Bubble can not grow at this conditions (ratio < 1)' write(77,*)'Please, change parameters Y1-Y5' if(imonitor.eq.1)then write(6,*)'Bubble can not grow at this conditions (ratio < 1)' write(6,*)'Please, change parameters Y1-Y5' endif goto 2001 endif cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Equilibrium final bubble radius (isothermal) c c Using bisection, find the root of an implicit function (func) c for r(t=t->infinity). This root is known to lie between x1 and x2 c c Adopted subroutine from "Numerical Recipes" cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc x1=1.0d+00 x2=((1.0d+00-ypar(9)*(c0-ypar(5)))*(ypar(2)*(c0- * ypar(5))/ypar(1)+ypar(3)+1.0d+00))**crt fmid=f1(x2,ypar,c0) f=f1(x1,ypar,c0) if(f*fmid.ge.0.0d+00) pause 'Root of final bubble size must be br * acketed for bisection.' if(f.lt.0.0d+00)then rfinal=x1 dx=x2-x1 else rfinal=x2 dx=x1-x2 endif do 40 j=1,40 dx=dx/2.0d+00 xmid=rfinal+dx fmid=f1(xmid,ypar,c0) if(fmid.le.0.0d+00)rfinal=xmid if(dabs(dx).lt.1.0d-03.or.fmid.eq.0.0d+00)goto 45 40 continue write(77,*)'Too many bisections for final bubble size' if(imonitor.eq.1)then write(6,*)'Too many bisections for final bubble size' endif goto 2001 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Writing input dimensional and dimentionless parameters c to the output file and monitor cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 45 if(imonitor.eq.1)write(6,'(/,3a,//)') * 'Name of output file - "',fname(1:iend-1),'"' if(jinfo.eq.0)goto 69 do 68 i=1,jinfo if(imonitor.eq.1)write(6,*)info(i) write(77,*)info(i) write(78,*)info(i) write(79,*)info(i) write(80,*)info(i) 68 continue if(imonitor.eq.1)write(6,'(//)') write(77,'(//)') write(78,'(/)') write(79,'(/)') write(80,'(/)') 69 do 70 i=1,10 write(77,50)'Y',i,' = ',ypar(i) if(imonitor.eq.1)write(6,50)'Y',i,' = ',ypar(i) 70 continue write(77,'(/)') if(vh.ne.0.0d+00)then write(77,53)'Initial depth of the melt (*)',depth0,'m' write(77,53)'Initial ambient pressure (*)',p0*1.0d-06,'MPa' write(77,53)'Magma rise rate (*)',vh,'m/s' endif write(77,53)'Initial gas conc. in melt (*)',c0,'none' write(77,53)'Initial temperature (*)',temp0-2.7315d+02,'grad.C' write(77,*) write(77,53)'Molecular weight of the gas (*)',wm,'kg/mole' write(77,53)'Density (*)',denmelt,'kg/m**3' if(enact.eq.0.0d+00)then write(77,'(x,a,/,13x,a)')'Viscosity -','(Hess, Dingwell, 1996)' else write(77,53)'Activation energy of viscosity',enact,'J/mole' write(77,53)'Coefficient of volatile effect',vcoef,'none' endif write(77,53)'Activation energy of diffusion',denact,'J/mole' write(77,53)'Diffusion free coefficient',dcoef,'none' write(77,53)'Liquidus temperature of VC (*)',tempvc-2.7315d+02, * 'grad.C' write(77,53)'VC temperature interval (*)',dTvc,'grad.K' write(77,*) write(77,*)'(*) parameters are not calculated from ypar(i)' write(77,*)'VC - vitrification or crystallization' write(77,'(/)') write(77,53)'Ambient pressure (final)',pa*1.0d-06,'MPa' write(77,54)'Initial bubble radius',r0,'m' write(77,53)'Initial sphere radius',s0,'m' write(77,*) write(77,53)'Surface tension',sigma,'J/m**2' write(77,54)'Initial viscosity',eta0,'Pa s' write(77,54)'Initial diffusion coefficient',dif0,'m**2/s' write(77,*) write(77,54)'Henry"s constant',hen,'1/Pa' write(77,54)'Temperature conductivity',xi,'m**2/s' write(77,53)'Heat capacity for melt (cm)',cm,'J/(kg grad.K)' write(77,53)'Heat capacity for gas (cp)',cp,'J/(kg grad.K)' write(77,53)'Latent heat of VC',heatvc*dTvc,'J/kg' write(77,*) write(77,53)'Initial heat of evaporation (**)',heat0,'J/mole' write(77,'(/)') if(iheat.eq.1)then write(77,*)'(**) heat of evaporation is P-T dependent' else write(77,*)'(**) heat of evaporation is not P-T dependent' endif write(77,'(//)') write(77,*)'The ratio of r0 to nuclear radius is' write(77,'(10x,f13.6,//)')ratio write(77,*)'Equilibrium final bubble radius at t->infinity', * ' and heat=const is' write(77,'(10x,f13.6,////)')rfinal if(imonitor.eq.1)then write(6,'(/)') write(6,*)'The ratio of r0 to nuclear radius is' write(6,'(10x,f13.6,/)')ratio write(6,*)'Equilibrium final bubble radius at t->infinity', * ' and heat=const is' write(6,'(10x,f13.6,///)')rfinal endif write(77,*)'Time(s), depth(m), radius, gas temperature(C), volati *le oversaturation,' write(77,*)'gas fraction, and fraction of decompression in the gro *wth rate:' cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc write(78,*)'Time(s), depth(m), surface tension and dynamic pressu *res,' write(78,'(a,/)')'total volatile content in the melt, and average * melt viscosity:' write(79,*)'Concentration, temperature(C), and viscosity profile *for a real normalized bubble wall:' write(80,*)'Time(s), depth(m), average melt temperature, then' write(80,'(a,/)')'dissipation, vitrification, vaporization and pd *V heating terms (in C).' cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if(imonitor.eq.1)then write(6,'(a)')' START OF MAIN CALCULATIONS' write(6,'(a,//)')' Please, wait!!!!!' endif cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c Main calculations c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Initializing y-increment, and other variables cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c dyu=1.0d+00/dfloat(m1) m=m1+1 do 110 i=1,m yu(i)=dyu*dfloat(i-1) 110 continue cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc kp=1 kk=kkf yyy=9.0d+00*ypar(1)**(2.0d+00/3.0d+00) t1=0.0d+00 td=r0*r0/dif0 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c generating the initial conditions c c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c 9000 kgrid=kgrid+1 grid=4.0d+01/4.0d+00**kgrid if(kgrid.le.10)then if(dabs(ratioc-1.0d+00).lt.dabs(check-1.0d+00))then check=ratioc grid1=4.0d+01/4.0d+00**(kgrid-1) endif else grid=grid1 if(grid.eq.0.0d+00)pause 'grid is not good' endif c write(6,'(3f20.8)')grid,ratioc,check ccccccccccccccccccccccccccccc beta=ypar(1)*grid beta11=beta*0.9d+00 ratioc=1.01d+00 ratiot=1.01d+00 ratioc11=0.9d+00 ratiot11=0.9d+00 coef=beta coef11=beta11 ccccccccccccccccccccccccccccc do 120 i=1,m con(i)=c0 co(i)=0.0d+00 temp(i)=temp0 to(i)=0.0d+00 120 continue cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc r1=1.0d+00 r=r1 r3=r*r*r ph=ph1 pg=ph+ypar(3) pg1=pg tg=1.0d+00 tg1=1.0d+00 dt=1.0d-08*sqrt(eta0)/ypar(1)**crt vr=1.0d-04/dt cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Start of main cycle to solve the system of equations cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c Updating the time step ccccccccccccccccccccccccccccccccccccc 10000 par=vr*dt/r**0.8d+00 if(par.lt.0.8d+00*time)dt=dt*1.1d+00 if(par.gt.1.0d+00*time)dt=dt/1.1d+00 if(kgrid.eq.1.and.ph.lt.1.0d+00+vp*dt)then if(imonitor.eq.1)write(6,*)'Make "time" smaller' write(77,*)'Make "time" smaller' goto 2001 endif if(r.gt.0.98d+00*rfinal)goto 126 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Updating parameters to generate the non-uniform mesh c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ratioc12=ratioc11 ratioc11=ratioc beta12=beta11 beta11=beta if(beta11.eq.beta12)beta11=beta11*1.00001d+00 bcoef=(ratioc11-ratioc12)/(beta11-beta12) beta=(1.0d+00-ratioc11+bcoef*beta11)/bcoef delta=1.0d-01/dsqrt(r) c delta=1.0d-01/r**crt if(beta/beta11.gt.(1.0d+00+delta))beta=beta11*(1.0d+00+delta) if(beta/beta11.lt.(1.0d+00-delta))beta=beta11*(1.0d+00-delta) beta1=beta*r3+1.0d+00 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ratiot12=ratiot11 ratiot11=ratiot coef12=coef11 coef11=coef if(coef11.eq.coef12)coef11=coef11*1.00001d+00 bcoef=(ratiot11-ratiot12)/(coef11-coef12) coef=(1.0d+00-ratiot11+bcoef*coef11)/bcoef if(coef/coef11.gt.1.1d+00)coef=coef11*1.1d+00 if(coef/coef11.lt.0.9d+00)coef=coef11*0.9d+00 beta2=coef*(r**5+1.0d+02)+1.0d+00 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c Generating or updating the non-uniform mesh for concenttation c c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc do 130 i=1,m y1(i)=(beta1+1.0d+00)-(beta1-1.0d+00)*(((beta1+1.0d+00)/( * beta1-1.0d+00))**(1.0d+00-yu(i))) y1(i)=y1(i)/(((beta1+1.0d+00)/(beta1-1.0d+00))**(1.0d+00- * yu(i))+1.0d+00) 130 continue do 140 i=1,m dydyn1(i)=(2.0d+00*beta1)/((beta1*beta1-(1.0d+00-y1(i))* * (1.0d+00-y1(i)))*dlog((beta1+1.0d+00)/(beta1-1.0d+00))) dydyn2(i)=dydyn1(i)*((2.0d+00*(y1(i)-1.0d+00))/ * (beta1*beta1-(1.0d+00-y1(i))*(1.0d+00-y1(i)))) dydyn1(i)=dydyn1(i)*dydyn1(i) 140 continue cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Generating or updating the non-uniform mesh for temperature c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc do 131 i=1,m y2(i)=(beta2+1.0d+00)-(beta2-1.0d+00)*(((beta2+1.0d+00)/( * beta2-1.0d+00))**(1.0d+00-yu(i))) y2(i)=y2(i)/(((beta2+1.0d+00)/(beta2-1.0d+00))**(1.0d+00- * yu(i))+1.0d+00) 131 continue do 141 i=1,m dydyn3(i)=(2.0d+00*beta2)/((beta2*beta2-(1.0d+00-y2(i))* * (1.0d+00-y2(i)))*dlog((beta2+1.0d+00)/(beta2-1.0d+00))) dydyn4(i)=dydyn3(i)*((2.0d+00*(y2(i)-1.0d+00))/ * (beta2*beta2-(1.0d+00-y2(i))*(1.0d+00-y2(i)))) dydyn3(i)=dydyn3(i)*dydyn3(i) 141 continue cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c Interpolation of temperature for concentration profile gridding c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 126 call spline(y2,temp,m,1.0d+30,1.0d+30,y2t) call spline(y1,con,m,1.0d+30,1.0d+30,y2c) do 142 i=1,m call splint(y2,temp,y2t,m,y1(i),tmp(i)) call splint(y1,con,y2c,m,y2(i),cont(i)) 142 continue cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c Calculating the viscosity integral c c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc vint=0.0d+00 do 145 i=1,m-1 vint=vint+vis(temp(i),cont(i))/eta0*dyu/dsqrt(dydyn3(i))/(y2(i)+ * ypar(1)*r3)**2 dif(i)=con(i)*dexp(-dcoef-denact/rg/tmp(i))/dif0 145 continue dif(m)=con(m)*dexp(-dcoef-denact/rg/tmp(m))/dif0 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c Calculating the heat of evaporation c c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if(iheat.eq.1)then heat=heatf(pg*pa,tg*temp0)/heat0 else heat=1.0d+00 endif cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c Calculating the correction to temperature profile c c due to vitrification or crystallization c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 10001 do 146 i=1,m dtv(i)=dt if(temp(i).lt.tempvc.and.temp(i).gt.tempvc-dTvc) * dtv(i)=dt/(1.0d+00+ypar(10)) 146 continue cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Calculating new ambient pressure c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if(ph.gt.1.0d+00+vp*dt)then ph=ph-vp*dt else ph=1.0d+00 endif cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c Start of bisection procedure c c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ick=0 pcoef=1.0d+00 194 do 195 j=1,100 if(j.le.3.or.fmid*fmid0.lt.0.0d+00)then fmid00=fmid0 r00=rd0 endif fmid0=fmid rd0=r if(j.eq.1)r=r*(1.0d+00+1.0d-10) if(j.eq.2)r=r+dabs(vr)*dt*pcoef if(j.eq.3.and.fmid0*fmid00.gt.0.0d+00)then if(ick.eq.1.and.r.gt.r1)then ccccc*ccccccccccccccccccccc*********Replacing the bubble if(dt.lt.1.1d-08/dabs(r1)/ypar(1)**crt)then if(imonitor.eq.1)write(6,*)'Bubble radius is not found...' goto 2001 endif ccccc*ccccccccccccccccccccc*********Reducing time step if(imonitor.eq.1)write(6,611) r=r1 dt=dt*5.0d-01 goto 10001 endif ccccc*ccccccccccccccccccccc*********Shrinking if(r.gt.r1*1.3d+00)then if(imonitor.eq.1)write(6,'("Trying to shrink the bubble")') r=r1*0.9d+00 ick=1 pcoef=1.0d+00 goto 194 endif ccccc*ccccccccccccccccccccc r=r/(1.0d+00+1.0d-10) pcoef=pcoef*2.0d+00 goto 194 endif if(j.gt.2)then if(dabs(fmid0/fmid00).lt.1.0d-02)then fder=5.0d-02 else if(dabs(fmid0/fmid00).lt.1.0d+02)then fder=fmid0/(fmid0-fmid00) else fder=9.5d-01 endif endif r=rd0-(rd0-r00)*fder endif r3=r*r*r ccccc*cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c Calculating the pressure c c c ccccc*cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc a1=ph+ypar(3)/r a2=ypar(1)*ypar(4)*r*r*vint pg=a1+a2*(r-r1)/dt ccccc*cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c Calculating the concentration c c c ccccc*cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc der=ypar(5)*dsqrt(pg)-c0 do 150 i=1,m b(i)=yyy*dif(i)*(y1(i)+ypar(1)*r3)**(4.0d+00/3.0d+00) 150 continue do 160 i=2,m-1 aa(i)=dydyn2(i)*dyu/2.0d+00-dydyn1(i) bb(i)=2.0d+00*dydyn1(i)+dyu*dyu/dt/b(i) cc(i)=-(dydyn2(i)*dyu/2.0d+00+dydyn1(i)) rr(i)=co(i)*dyu*dyu/dt/b(i) 160 continue bb(1)=-1.0d+00 cc(1)=1.0d+00 rr(1)=dyu*der/dsqrt(dydyn1(1)) bb(m)=1.0d+00 cc(m)=0.0d+00 rr(m)=0.0d+00 c call tridag(aa,bb,cc,rr,cn,m) c dcdy=(cn(1)-co(1))/b(1)/dt dcdy1=dydyn2(1)*((-3.0d+00*cn(1)+4.0d+00*cn(2)-cn(3))/ * (2.0d+00*dyu))+dydyn1(1)*((5.0d+00*cn(1)-1.1d+01*cn(2)+ * 7.0d+00*cn(3)-cn(4))/(4.0d+00*dyu*dyu)) ccccc*cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c Calculating the temperature c c c ccccc*cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c a4=9.0d+00*ypar(1)*ypar(2)*dif(1)*r3*r*dcdy tg=pg*r3/(pg1*r1*r1*r1/tg1+a4*dt) c ccccc*cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c Calculating the temperature profile c c c ccccc*cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc der=tg-1.0d+00 do 170 i=1,m b(i)=yyy*ypar(6)*(y2(i)+ypar(1)*r3)**(4.0d+00/3.0d+00) 170 continue ymul=ypar(1)*ypar(1)*ypar(7)*r3*r*((r-r1)/dt)**2*vint do 180 i=2,m-1 aa(i)=dydyn4(i)*dyu/2.0d+00-dydyn3(i) bb(i)=2.0d+00*dydyn3(i)+dyu*dyu/dtv(i)/b(i) cc(i)=-(dydyn4(i)*dyu/2.0d+00+dydyn3(i)) rr(i)=dyu*dyu/b(i)*(to(i)/dtv(i)-ymul) 180 continue bb(1)=-1.0d+00 cc(1)=1.0d+00 rr(1)=dyu*der/dsqrt(dydyn3(1)) bb(m)=1.0d+00 cc(m)=0.0d+00 rr(m)=0.0d+00 c call tridag(aa,bb,cc,rr,tn,m) c dtdy=((tn(1)-to(1))/dtv(1)+ymul)/b(1) dtdy1=dydyn4(1)*((-3.0d+00*tn(1)+4.0d+00*tn(2)-tn(3))/ * (2.0d+00*dyu))+dydyn3(1)*((5.0d+00*tn(1)-1.1d+01*tn(2)+ * 7.0d+00*tn(3)-tn(4))/(4.0d+00*dyu*dyu)) ccccc*cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c End part of the bisection method c ccccc*cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc a3=9.0d+00*ypar(1)*ypar(2)*ypar(8)*tg*r/pg*(ypar(6)*dtdy- * ypar(9)*heat*dif(1)*dcdy+ypar(7)/2.7d+01/ypar(1)/ypar(4)/r* * (pg-pg1)/dt) fmid=tg-tg1-a3*dt if(dabs(fmid)+dabs(r-rd0).lt.1.0d-05.or.fmid.eq.0.0d+00)goto 196 195 continue if(kgrid.le.10)goto 196 if(dt.lt.1.1d-08/dabs(r1)/ypar(1)**crt)then if(imonitor.eq.1)write(6,*)'Bubble radius is not found...' goto 2001 endif if(imonitor.eq.1)write(6,63) r=r1 dt=dt*5.0d-01 goto 10001 196 if(ick.eq.1.and.r.lt.r1.and.imonitor.eq.1)write(6,64) * r/r1*1.0d+02 ratioc=dcdy/dcdy1 ratiot=dtdy/dtdy1 if(kgrid.le.10)goto 9000 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Verification growth conditions c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc vr=(r-r1)/dt if(vh.ne.0.0d+00.and.ic.eq.10.and.ratio.gt.1.0d+00)cn(1)= * (r*r*r*pg/tg-ph-ypar(3))*ypar(1)/ypar(2) oversat=c0-cn(1)-dsqrt(ypar(3)/r+ph)*ypar(5) if(r.lt.r1)pause 'Radius decrease' if(ph.eq.1.0d+00.and.oversat.lt.1.0d-08.or.vr.lt.1.0d-08)then write(77,*)'Very low oversaturation and growth rate' if(imonitor.eq.1)write(6,*) * 'Very low oversaturation and growth rate' goto 2000 endif cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Calculating temperature terms cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cmg=cm*(1.0d+00+cn(1)) tmelt=(1.0d+00-tn(1))*temp0 tevap=tevap+heat*heat0/wm*(cn(1)-co(1))/cmg tpdv=tpdv+ypar(1)*r**3*pa*(pg1-pg)/cmg/denmelt tdis=tdis+(pg-a1)*pa*ypar(1)*(r**3-r1**3)/cmg/denmelt if(tmelt.gt.tempvc)tvit=0.0d+00 if(tmelt.lt.tempvc)tvit=heatvc*(tempvc-tmelt)/cmg if(tmelt.lt.tempvc-dTvc)tvit=heatvc*dTvc/cmg cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Calculating decompression and diffusion terms in growth rate cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if(vh.ne.0.0d+00)then drdif=r1*((1.0d+00+(cn(1)-co(1))*ypar(2)/ypar(1)/ph/r3)**crt- * 1.0d+00)/(r-r1) if(ndec.ne.1)then drdec=r1*(((ph+vp*dt)/ph)**crt-1.0d+00)/(r-r1) if(ph-vp*dt.lt.1.0d+00)then deck=pg*r3 deca=drdec/ypar(1)/ypar(4)/vint/r3/r/r/vr ndec=1 endif else drdec=ddim((deck-r3)*deca,0.0d+00) endif decfr=drdec/(drdec+drdif) decint=decint+decfr*(r3-r1**3) endif cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c Updating the variables c c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc t1=t1+dt tg1=tg pg1=pg r1=r do 200 i=1,m co(i)=cn(i) to(i)=tn(i) 200 continue if(vh.ne.0.0d+00)depth=(ph-1.0d+00)*pa/denmelt/grav cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Calculating actual concentration cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc con(1)=dsqrt(dydyn1(1))/(2.0d+00*dyu)*(-3.0d+00*cn(1)+ * 4.0d+00*cn(2)-cn(3))+c0 do 165 i=2,m-1 con(i)=dsqrt(dydyn1(i))/(2.0d+00*dyu)*(cn(i+1)-cn(i-1))+c0 165 continue con(m)=dsqrt(dydyn1(m))/dyu*(cn(m)-cn(m-1))+c0 ratioc=dcdy/dcdy1 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Calculating actual temperature cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc temp(1)=(dsqrt(dydyn3(1))/(2.0d+00*dyu)*(-3.0d+00*tn(1)+ * 4.0d+00*tn(2)-tn(3))+1.0d+00)*temp0 do 190 i=2,m-1 temp(i)=(dsqrt(dydyn3(i))/(2.0d+00*dyu)*(tn(i+1)-tn(i-1)) * +1.0d+00)*temp0 190 continue temp(m)=(dsqrt(dydyn3(m))/dyu*(tn(m)-tn(m-1))+1.0d+00)*temp0 ratiot=dtdy/dtdy1 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c Calculating the bubble fraction in the system c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc v=r3/(r3+1.0d+00/ypar(1)) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Writing pressure on the screen cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if(r-1.0d+00.gt.(rfinal-1.0d+00)*dfloat(ic)*1.0d-01)then ic=ic+1 if(imonitor.eq.1)then write(6,*) write(6,*)'Time= ',t1*td write(6,*) write(6,*)'Surface tension and dynamic pressures, gas', * ' temperature:' write(6,*)ypar(3)/r,a2*vr,tg*temp0-2.7315d+02 endif cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Writing concentration and temperature profile c in the bubble wall normalized to 1.0 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if(ic.gt.1)then write(79,*) write(79,*)'Time, depth and radius (',ic-1,'0% of Rf) are ' write(79,59)t1*td,depth,r write(79,*) if(imonitor.eq.1)then write(6,*) write(6,*)'Concentration, temperature, and viscosity profile for' write(6,*)'a real normalized bubble wall (R is',ic-1,'0% of Rf):' endif do 210 i=0,20 x=dfloat(i)/2.0d+01 z=((((1.0d+00/ypar(1)+r3)**crt-r)*x+r)**3-r3)*ypar(1) j=int(z*dfloat(m1))+1 visc=vis(tmp(j),con(j)) write(79,58)y1(j),con(j),tmp(j)-2.7315d+02,visc if(imonitor.eq.1)write(6,58)y1(j),con(j),tmp(j)-2.7315d+02,visc 210 continue endif if(imonitor.eq.1)then write(6,*) print *,'Time(s), depth(m), radius, temperature(C) and three prec *ision indicators (=1.0)' endif endif cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c Writting the output c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if(((r-1.0d+00).lt.1.0d-02*(rfinal-1.0d+00)).or. * ((r-1.0d+00).gt.9.7d-01*(rfinal-1.0d+00)))then if(kk.eq.kkf)then kk=0 goto 1999 else kk=kk+1 endif endif if(r.gt.rfinal*dfloat(kp)/1.0d+02)then kp=kp+1 goto 1999 else goto 10000 endif 1999 if(t1*td.lt.1.0d+06)then write(77,55) t1*td,depth,r,tg*temp0-2.7315d+02,oversat,v,decfr else write(77,551)t1*td,depth,r,tg*temp0-2.7315d+02,oversat,v,decfr endif avis=0.0d+00 do 220 i=1,m1 avis=avis+vis(tmp(i),con(i))*(y1(i+1)-y1(i)) 220 continue if(t1*td.lt.1.0d+06)then write(78,60)t1*td,depth,ypar(3)/r,a2*vr,c0-cn(1),avis else write(78,601)t1*td,depth,ypar(3)/r,a2*vr,c0-cn(1),avis endif if(t1*td.lt.1.0d+06)then write(80,62)t1*td,depth,tmelt-2.7315d+02,tdis,tvit,-tevap,-tpdv else write(80,621)t1*td,depth,tmelt-2.7315d+02,tdis,tvit,-tevap,-tpdv endif ratio=r/((co(1)*ypar(2)/ypar(1)+ph+ypar(3))*tg/pg)**crt if(dtdy1.eq.0.0d+00)then dtdy1=1.0d+00 dtdy=1.0d+00 endif if(imonitor.eq.1)write(6,57) t1*td,depth,r,tg*temp0-2.7315d+02, * ratioc,ratiot,ratio if(r.lt.rfinal*1.2d+00)goto 10000 if(imonitor.eq.1)print *,'raduis is gt.1.2 rfinal' 2000 if(t1*td.lt.1.0d+06)then write(77,55) t1*td,depth,r,tg*temp0-2.7315d+02,oversat,v,decfr else write(77,551)t1*td,depth,r,tg*temp0-2.7315d+02,oversat,v,decfr endif write(77,'(a,f7.2)')'Integrated decompression bubble growth (%) = *',decint/(r3-1.0d+00)*1.0d+02 if(imonitor.eq.1)write(6,'(/,3a,/)') * 'Name of output file - "',fname(1:iend-1),'"' cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 2001 close(unit=77) close(unit=78) close(unit=79) end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c Subroutines c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Function for reading data cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c double precision function dataread(poisk) implicit double precision (a-h,o-z) character*(*) poisk character*30 line rewind 21 do 32 i=1,1000 read(21,'(a30)',end=33) line ind=index(line,'=') if(line(1:ind-1).eq.poisk) goto 34 32 continue 33 write(6,*)'Parameter "',poisk,'" is not found in "erupt.dat"' call abort 34 read(line(ind+1:30),*) dataread return end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Function for viscosity cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c double precision function vis(tem,conc) implicit double precision(a-h,o-z) common enact,rg,vcoef,tempvc,dTvc if(tem.gt.tempvc)then if(enact.eq.0.0d+00)then vln=dlog(conc*1.0d+02) vis=1.0d+01**(-3.545d+00+8.33d-01*vln+(9.601d+03- * 2.368d+03*vln)/(tem-195.7d+00-32.25d+00*vln)) else vis=10**(-4.5d+00)*dexp(enact*(1.0d+00-vcoef*conc)/rg/tem) endif else vis1=1.0d+01**(-4.5d+00)*dexp(enact*(1.0d+00-vcoef*conc)/rg/ * tempvc) enact2=rg*tempvc*(tempvc-dTvc)/dTvc*dlog(1.0d+12/vis1) hvis02=dlog(vis1)-enact2/rg/tempvc vis=dexp(hvis02+enact2/rg/tem) vis=dmin1(vis,1.0d+12) endif return end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Function for heat of vaporization cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c double precision function heatf(pres,tem) implicit double precision(a-h,o-z) dimension ak(6,3) data ak/6.736d+04,2.47d+03,-2.349d+03,-1.131d+03,2.591d+02, * -1.089d+01,-7.185d+01,-3.534d+00,3.65d+00,1.639d+00,-4.865d-01, * 2.939d-02,2.365d-02,1.744d-03,-1.594d-03,-6.038d-04,1.984d-04, * -1.27d-05/ press=pres*1.0d-06 heatf=0.0d+00 do 41 j=1,6 atem=0.0d+00 do 42 i=1,3 atem=atem+ak(j,i)*tem**(i-1) 42 continue heatf=heatf+atem*(dlog(press))**(j-1) 41 continue return end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Function "f1(x1,ypar,c0)" is an implicit function for equilibrium c bubble radius when time intends to infinity cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc double precision function f1(x,ypar,c0) implicit double precision(a-h,o-z) dimension ypar(9) dc=c0-ypar(5)*dsqrt(1.0d+00+ypar(3)/x) f1=x*x*x*(1.0d+00+ypar(3)/x)/(1.0d+00-ypar(9)*dc)-ypar(3)-ypar(2) * /ypar(1)*dc-1.0d+00 return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Subroutine TRIDAG solve the tridiagonal matrix which has nonzero c elements only on the diagonal plus and minus one column. c This subroutine is taken from the book "Numerical Recipes" cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine tridag(a,b,c,r,u,n) implicit double precision(a-h,o-z) parameter (nmax=1000) dimension gam(nmax),a(n),b(n),c(n),r(n),u(n) if(b(1).eq.0.)pause bet=b(1) u(1)=r(1)/bet do 300 j=2,n gam(j)=c(j-1)/bet bet=b(j)-a(j)*gam(j) if(bet.eq.0.)pause u(j)=(r(j)-a(j)*u(j-1))/bet 300 continue do 310 j=n-1,1,-1 u(j)=u(j)-gam(j+1)*u(j+1) 310 continue return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Subroutine SPLINE is supplimentary for interpolation cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine spline(x,y,n,yp1,ypn,y2) implicit double precision(a-h,o-z) parameter (nmax=1000) dimension x(n),y(n),y2(n),u(nmax) if (yp1.gt.0.99d+30) then y2(1)=0.0d+00 u(1)=0.0d+00 else y2(1)=-0.5d+00 u(1)=(3.0d+00/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1) endif do 411 i=2,n-1 sig=(x(i)-x(i-1))/(x(i+1)-x(i-1)) p=sig*y2(i-1)+2.0d+00 y2(i)=(sig-1.0d+00)/p u(i)=(6.0d+00*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1)) * /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p 411 continue if (ypn.gt.0.99d+30) then qn=0. un=0. else qn=0.5d+00 un=(3.0d+00/(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1))) endif y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.) do 412 k=n-1,1,-1 y2(k)=y2(k)*y2(k+1)+u(k) 412 continue return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Subroutine SPLINT is cubic spline interpolation cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine splint(xa,ya,y2a,n,x,y) implicit double precision(a-h,o-z) dimension xa(n),ya(n),y2a(n) klo=1 khi=n 401 if (khi-klo.gt.1) then k=(khi+klo)/2 if(xa(k).gt.x)then khi=k else klo=k endif goto 401 endif h=xa(khi)-xa(klo) if (h.eq.0.) pause 'bad xa input.' a=(xa(khi)-x)/h b=(x-xa(klo))/h y=a*ya(klo)+b*ya(khi)+ * ((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**2)/6.0d+00 return end