subroutine init C Self-Similar Driven Wave in 2 or 3 dimensions C Hardwired for spherical geometry C Boundary conditions are hardwired into sweepX.f C DATA is initially blocked in j C aug00 blondin C======================================================================= implicit none C GLOBALS include 'global.h' include 'zone.h' include 'mpif.h' C LOCALS integer i, j, k, mpierr, ibuffer, ndata, n, seed, icontact & ,nhi, nb, nlo real ridt, rodt, xmin, xmax, ymin, ymax, zmin, zmax, xvel & ,yvel, zvel, width, widthy, widthz, ratio & ,delx, ran3, dleft, fluct, sound, log_r_max, log_r_min & ,del_log_r, angle, delphi, bubr, Bub_radius, Bubc2 & ,radi1, radi2, radi3, radi4, radi5, radi6, radi7, radi8 & ,radi9, rado, cphi, sphi, cphip, cphim, sphip, sphim & ,cthetap, cthetam, sthetap, sthetam, stheta, ctheta & ,theta, phi, radius, density, pressure, velocity & ,frac, omf, contact real input_x(2101), input_r(2101) & ,input_v(2101), input_p(2101) C------------------------------------------------------------------------ pi = 2.0 * asin(1.0) C Define the problem... ndim = 2 C====================================================================== C Read input data open(26,file='x120d') ndata = 2101 npower = 12 spower = 0 gam = 1.1 ! 4. / 3. do n = 1, ndata read(26,*) input_x(n), input_r(n), input_p(n), input_v(n) enddo contact = input_x(1001) close(26) C====================================================================== xmin = 5.0 * (input_x(1) - 1.0) + 1.0 xmax = 1.01 ymin = 0.5 * pi - 0.5*(xmax-xmin) ymax = 0.5 * pi + 0.5*(xmax-xmin) zmin = 0.0 zmax = (ymax-ymin) courant= 0.4 gamm = gam - 1.0 yspan = ymax - ymin zspan = zmax - zmin C Set up grid coordinates and parabolic coefficients call grid(imax,xmin,xmax,zxa,zxc,zdx,zparax) call grid(jmax,ymin,ymax,zya,zyc,zdy,zparay) call grid(kmax,zmin,zmax,zza,zzc,zdz,zparaz) if(ndim.eq.2) then zzc(1) = 0.0 zspan = 0.0 endif C====================================================================== C Log parameters of problem in history file if(mype.eq.0) then write (8,*) 'SSDW in ', ndim, ' dimensions.' write (8,*) ' on a moving grid' write (8,*) write (8,*) 'Adiabatic index, gamma = ', gam write (8,*) 'n = ',npower write (8,*) 's = ',spower write (8,*) '------------------------------------------' endif C initialize grid: icontact= (contact-zxc(1))/zdx(1) seed = -431 - mype*145 time = input_x(1) / input_v(1) do i = 1, imax if(zxc(i).lt.input_x(1)) then ratio = (input_x(1)/zxc(i))**npower density = input_r(1)*ratio pressure = input_p(1) velocity = zxc(i) / time else nhi = 2 do while(input_x(nhi).lt.zxc(i)) nhi = nhi + 1 enddo nlo = nhi - 1 frac = (zxc(i)-input_x(nlo))/(input_x(nhi)-input_x(nlo)) omf = 1.0 - frac density = omf * input_r(nlo) + frac * input_r(nhi) pressure = omf * input_p(nlo) + frac * input_p(nhi) velocity = omf * input_v(nlo) + frac * input_v(nhi) endif do k = 1, kmax do j = 1, js zro(i,j,k) = density zpr(i,j,k) = pressure zux(i,j,k) = velocity zuy(i,j,k) = 0. zuz(i,j,k) = 0. zcl(i,j,k) = 1.0 if(zxc(i).gt.contact) zcl(i,j,k) = 0.0 enddo enddo enddo do k = 1, kmax do j = 1, js do i = icontact-20, icontact+20 sound = sqrt(gam*zpr(i,j,k)/zro(i,j,k)) fluct = sound * 2.0*(ran3(seed)-0.5) zux(i,j,k) = zux(i,j,k) + fluct fluct = sound * 2.0*(ran3(seed)-0.5) zuy(i,j,k) = fluct enddo enddo enddo C set parameter for powerlaw ejecta ejecta = zro(1,1,1)*time**3 * (zxc(1)/time)**npower pscale = zpr(imax-10,1,1)*time**2 C set outer boundary conditions dotflo = zro(imax,1,1)*zxc(imax)**spower potflo = zpr(imax,1,1)*zxc(imax)**spower uotflo = 0.0 votflo = 0.0 wotflo = 0.0 C******************************************************************* C set value of cmin & cmax for normalizing density in movie images cmax = dleft * 0.2 cmin = 0.0 C######################################################################## C Compute Courant-limited timestep ridt = 0. if(ndim.eq.2) then k = 1 do j = 1, js do i = 1, imax widthy = zdy(mype*js+j) * zxc(i) width = min(zdx(i),widthy) svel = sqrt(gam*zpr(i,j,k)/zro(i,j,k))/width xvel = abs(zux(i,j,k)) / zdx(i) yvel = abs(zuy(i,j,k)) / widthy ridt = max(xvel,yvel,svel,ridt) enddo enddo else do k = 1, kmax do j = 1, js do i = 1, imax widthy = zdy(mype*js+j) * zxc(i) widthz = zdz(k) * zxc(i) * sin(zyc(mype*js+j)) width = min(zdx(i),widthy,widthz) svel = sqrt(gam*zpr(i,j,k)/zro(i,j,k))/width xvel = abs(zux(i,j,k)) / zdx(i) yvel = abs(zuy(i,j,k)) / widthy zvel = abs(zuz(i,j,k)) / widthz ridt = max(xvel,yvel,zvel,svel,ridt) enddo enddo enddo endif call MPI_ALLREDUCE( ridt, rodt, 1, MPI_DOUBLE_PRECISION, MPI_MAX, & MPI_COMM_WORLD, mpierr ) dt = courant / rodt return end c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! function ran3(idum) c Returns a uniform random deviate between 0.0 and 1.0. Set IDUM to any c negative value to initialize or reinitialize the sequence. integer idum, mbig, mseed, mz real fac c real ran3 parameter (mbig=1000000000,mseed=161803398,mz=0,fac=1./mbig) integer i,iff,ii,inext,inextp,k, mj,mk,ma(55) save iff,inext,inextp,ma data iff /0/ if(idum.lt.0.or.iff.eq.0) then iff=1 mj=mseed-iabs(idum) mj=mod(mj,mbig) ma(55)=mj mk=1 do 11 i=1,54 ii=mod(21*i,55) ma(ii)=mk mk=mj-mk if(mk.lt.mz)mk=mk+mbig mj=ma(ii) 11 continue do 13 k=1,4 do 12 i=1,55 ma(i)=ma(i)-ma(1+mod(i+30,55)) if(ma(i).lt.mz)ma(i)=ma(i)+mbig 12 continue 13 continue inext=0 inextp=31 idum=1 endif inext=inext+1 if(inext.eq.56)inext=1 inextp=inextp+1 if(inextp.eq.56)inextp=1 mj=ma(inext)-ma(inextp) if(mj.lt.mz)mj=mj+mbig ma(inext)=mj ran3=mj*fac return end