c -*- mode: fortran -*- setup.F c Homogenous turbulence setup c---------------------------------------------------------------- c Setup globals - Here, for an entire memory manager's set of c tiles. This should be called once for each tile set. c c nodeSize: xyz zones on this MM c prbSize: xyz whole domain size, all MMs c xl/yl/zl: coordinate mesh for THE WHOLE DOMAIN, here initialized c others are physics parameters from rdparams() c subroutine setupturbglobals include 'iq.h' include 'globals.h' c 2601 wavenumbers is enough for kmax=8 integer MX_WAVE parameter (MX_WAVE=2604) real snn(MX_AXIS), css(MX_AXIS), aro(MX_WAVE), ape(MX_WAVE), & aror(MX_WAVE), aroi(MX_WAVE), aper(MX_WAVE), apei(MX_WAVE), & axr(MX_WAVE), ayr(MX_WAVE), azr(MX_WAVE), axi(MX_WAVE), & ayi(MX_WAVE), azi(MX_WAVE) common /setTurbcmn/ isetup, nlen, snn,css, & aro,ape,aror,aroi,aper,apei, & axr,ayr,azr, axi,ayi,azi integer nxzons,nyzons,nzzons c write(0,*) "entering setupturbglobals gamma=", gamma c write(0,*) "xxl/yyl/zzzl=" c do i= 1-nbdyf, 10 c write(0,*) i, xxl(i), yyl(i), zzl(i) c enddo c write(0,*) "xxl/yyl/zzzl=" c do i= nx-10, nx+nbdyf c write(0,*) i, xxl(i), yyl(i), zzl(i) c enddo cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Set sin and cosine arrays c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c NOTE: setvec as it is (15 July 1993) works for any arrangment of c c sub-meshes (identical or not, cubical or not). c c However, the entire simulation mesh must be uniform and cubical. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc pi = 3.1415926535897932 pi2 = 2. * pi cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c SET CONSTANTS c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c grav = acceleration do to gravity in +r direction (grav = 0) c c beta = initial temperature gradient at top. c c flux = vertical heat flux everywhere in the initial state, c c it is also imposed at the bottom boundary at all times c c eitop = internal energy in center of zone just above top boundary c c deidr0 = vertical slope of internal energy imposed at bottom = 0 c c heatcv = heat capacity at constant volume = 1 c c rhoavg = mean mass density c c pavg = mean presure c c c0 = mean sound speed c c ei0 = mean internal energy c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c0 = 1. rhoavg = 1. pavg = (rhoavg * c0**2) / gamma c ei0 = pavg / ((gamma - 1.) * rhoavg) c eitop = ei0 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c This routine sets the presure, density, and velocity fields for c c homogenous turbulence. Density and presure are constant. c c density = 1, speed of sound = 1. c c----------------------------------------------------------------------c c Inputs are: c c read (19,7) acard, dvamp, bvamp, tcamp c c read (19,9) acard,iseed,kmax,velamp,ratsc,wdthk,roamb,pgamb c c----------------------------------------------------------------------c c dvamp = Amplitude of Dynamic Viscousity c c bvamp = Amplitude of Bulk Viscousity c c tcamp = Amplitude of Thermal Conductivity c c iseed = Random number generator seed c c velamp = Root (Ensemble Mean) (Spatial Mean) Square velocity c c ratsc = (solinoidal power) / (compresible power) c c fkmax = Maximum mode excited c c wdthk = Width of excited modes c c roamb = fractional RMS desity fluctuations c c pgamb = fractional RMS presure fluctuations c c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Parameter c c Peak excited mode: fk0 = 0. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc fk0 = 0. c sea added ^^^^^^ nlen = max(nx,ny,nz) dxk = pi2 / float(nlen) do 60 i = 1, nlen arg = dxk * float(i-1) snn(i) = sin(arg) css(i) = cos(arg) 60 continue call setrmy(iseed) c Now setup ALL the required wavenumber amplitudes. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Initial Fourier modes are from a Gaussian distribution c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Set scaling factors for density, presure, and velocity c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c call setrmy(iseed) kmax2 = kmax ** 2 neglim = kmax * nlen wdthi = 1. / wdthk vel2 = 0. do 550 kz = -kmax, kmax do 550 ky = -kmax, kmax do 550 kx = 0, kmax cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c the following is for debug of 2-period waves only c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c do 550 kz = -kmax, kmax,2 c do 550 ky = -kmax, kmax,2 c do 550 kx = 0, kmax,2 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if(kx**2+ky**2+kz**2 .gt. kmax2) go to 550 if((kx .eq. 0) .and. (ky .lt. 0)) go to 550 if((kx .eq. 0) .and. (ky .eq. 0) .and. (kz .le. 0)) go to 550 fk = sqrt(float(kx**2 + ky**2 + kz**2)) vel22 = (fk * exp(-((fk - fk0)*wdthi) ** 2)) ** 2 vel2 = vel2 + vel22 550 continue c roamp = rhoavg * roamb * sqrt(4. / vel2) pgamp = pavg * pgamb * sqrt(4. / vel2) c ------------- change for 3-D was made here ----------------- vels = velamp * sqrt(2. * ratsc / (vel2 * (1. + ratsc))) velc = velamp * sqrt(4. / (vel2 * (1. + ratsc))) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Run through all of the wave vectors "k" c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c debugging: c ncount = 0 c as2tot = 0. c ac2tot = 0. c ro2tot = 0. c pe2tot = 0. c ax2tot = 0. c ay2tot = 0. c az2tot = 0. disp = 1. / float(2 * nz) c c // wavenumber: iw = 0 do 900 kz = -kmax, kmax do 800 ky = -kmax, kmax do 700 kx = 0, kmax iw = iw + 1 if ( iw .gt. MX_WAVE ) then print *, MX_WAVE,' wave table size exceeded' stop endif c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c the following is for debug of 2-period waves only c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c do 900 kz = -kmax, kmax,2 c do 800 ky = -kmax, kmax,2 c do 700 kx = 0, kmax,2 c kx = mod(iseed , 2) c ky = mod(iseed/2, 2) c kz = mod(iseed/4, 2) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if(kx**2+ky**2+kz**2 .gt. kmax2) go to 700 if((kx .eq. 0) .and. (ky .lt. 0)) go to 700 if((kx .eq. 0) .and. (ky .eq. 0) .and. (kz .le. 0)) go to 700 phase0 = disp * float(kx + ky + kz) fk = sqrt(float(kx**2 + ky**2 + kz**2)) c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c set up density and presure amplitudes (real + imaginary) c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc samp = fk * exp(-((fk - fk0)*wdthi) ** 2) aro(iw) = roamp * samp * erfinv(ranfmy(dumb)) ape(iw) = pgamp * samp * erfinv(ranfmy(dumb)) c$ phase = pi2 * (erfinv(ranfmy(dumb)) - phase0) phase = pi2 * (ranfmy(dumb) - phase0) aror(iw) = aro(iw) * cos(phase) aroi(iw) = aro(iw) * sin(phase) c$ phase = pi2 * (erfinv(ranfmy(dumb)) - phase0) phase = pi2 * (ranfmy(dumb) - phase0) aper(iw) = ape(iw) * cos(phase) apei(iw) = ape(iw) * sin(phase) c c ro2tot = ro2tot + aro ** 2 c pe2tot = pe2tot + ape ** 2 c write (6,*) 'ro2tot,aro', ro2tot,aro c write (6,*) 'roamp,samp', roamp, samp cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c set up velocity amplitudes (real + imaginary) c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc fkinv = 1. / fk fkx = fkinv * float(kx) fky = fkinv * float(ky) fkz = fkinv * float(kz) if((kx .ne. 0) .or. (ky .ne. 0)) then fkxyinv = 1. / sqrt(fkx**2 + fky**2) fk1x = fky * fkxyinv fk1y = -fkx * fkxyinv fk1z = 0. else fk1x = 1. fk1y = 0. fk1z = 0. endif fk2x = fky * fk1z - fkz * fk1y fk2y = fkz * fk1x - fkx * fk1z fk2z = fkx * fk1y - fky * fk1x c vel0 = fk * exp(-((fk - fk0)*wdthi) ** 2) as1 = vels * vel0 * erfinv(ranfmy(dumb)) as2 = vels * vel0 * erfinv(ranfmy(dumb)) ac = velc * vel0 * erfinv(ranfmy(dumb)) c$ phase = pi2 * (erfinv(ranfmy(dumb)) - phase0) phase = pi2 * (ranfmy(dumb) - phase0) as1r = as1 * cos(phase) as1i = as1 * sin(phase) c$ phase = pi2 * (erfinv(ranfmy(dumb)) - phase0) phase = pi2 * (ranfmy(dumb) - phase0) as2r = as2 * cos(phase) as2i = as2 * sin(phase) c$ phase = pi2 * (erfinv(ranfmy(dumb)) - phase0) phase = pi2 * (ranfmy(dumb) - phase0) acr = ac * cos(phase) aci = ac * sin(phase) c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c the following is for debug of single wave only c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c as1r = 0. c as1i = 0. c as2r = 0. c as2i = 0. c acr = velamp c aci = 0. c aror = roamb c aroi = 0. c aper = pgamb c apei = 0. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c the preceeding is for debug only c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c axr(iw) = fkx * acr + fk1x * as1r + fk2x * as2r ayr(iw) = fky * acr + fk1y * as1r + fk2y * as2r azr(iw) = fkz * acr + fk1z * as1r + fk2z * as2r axi(iw) = fkx * aci + fk1x * as1i + fk2x * as2i ayi(iw) = fky * aci + fk1y * as1i + fk2y * as2i azi(iw) = fkz * aci + fk1z * as1i + fk2z * as2i c dbging: c as2tot = as2tot + as1 ** 2 + as2 ** 2 c ac2tot = ac2tot + ac ** 2 c ax2tot = ax2tot + axr ** 2 + axi ** 2 c ay2tot = ay2tot + ayr ** 2 + ayi ** 2 c az2tot = az2tot + azr ** 2 + azi ** 2 c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Add one wave-number (sin + cos) onto mesh c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c dxk = pi2 / float(nlen) c dxkx = dxk * kx c dxky = dxk * ky c dxkz = dxk * kz 700 continue 800 continue 900 continue isetup = 1 c write(0,*) "setupturb complete, iw=", iw return end c ---------------------------------------------------------------- subroutine readblock(jtx,jty,jtz) implicit none INCLUDE 'iq.h' INCLUDE 'thingftn.h' integer itx,ity,itz, ierr integer status(NWSTATUS_THING) integer jtx,jty,jtz real rho(npx,npy,npz),prs(npx,npy,npz),uux(npx,npy,npz) real uuy(npx,npy,npz),uuz(npx,npy,npz) common /rst/ rho,prs,uux,uuy,uuz integer ioff, nbytes character*(256) name real dmin,dmax integer ix,iy,iz c convert to block number offsets itx= jtx/npx ity= jty/npy itz= jtz/npz name= "BOF2_h13---0120aaa" c 123456789012345678 name(9:11)= "???" name(16:16) = char( ichar('a') + itx) name(17:17) = char( ichar('a') + ity) write(0,*) "reading ", name(1:18), itx,ity, itz nbytes = npz * 4 * npx * npy ioff = (itz) * nbytes name(9:11)= "RHO" ierr= read_thing( name, ioff, rho, nbytes, status ) if ( ierr .ne. 0 ) goto 666 ierr= wait_thing( status, FOREVER ) if ( ierr .ne. nbytes ) goto 667 name(9:11)= "PRS" ierr= read_thing( name, ioff, prs, nbytes, status ) if ( ierr .ne. 0 ) goto 666 ierr= wait_thing( status, FOREVER ) if ( ierr .ne. nbytes ) goto 667 name(9:11)= "Vx_" ierr= read_thing( name, ioff, uux, nbytes, status ) if ( ierr .ne. 0 ) goto 666 ierr= wait_thing( status, FOREVER ) if ( ierr .ne. nbytes ) goto 667 name(9:11)= "Vy_" ierr= read_thing( name, ioff, uuy, nbytes, status ) if ( ierr .ne. 0 ) goto 666 ierr= wait_thing( status, FOREVER ) if ( ierr .ne. nbytes ) goto 667 name(9:11)= "Vz_" ierr= read_thing( name, ioff, uuz, nbytes, status ) if ( ierr .ne. 0 ) goto 666 ierr= wait_thing( status, FOREVER ) if ( ierr .ne. nbytes ) goto 667 write(0,*) "first bytes rho,prs,vx,vy,vz=",rho(1,1,1), & prs(1,1,1), uux(1,1,1), uuy(1,1,1), uuz(1,1,1) write(0,*) "last bytes rho,prs,vx,vy,vz=",rho(npx,npy,npz), & prs(npx,npy,npz), uux(npx,npy,npz), uuy(npx,npy,npz), uuz(npx,npy,npz) c More thorough error checking: dmin= 1.0e+6 dmax= -dmin do iz= 1, npx do iy= 1, npy do ix = 1, npx if ( rho(ix,iy,iz) .lt. dmin ) dmin= rho(ix,iy,iz) if ( rho(ix,iy,iz) .gt. dmax ) dmax= rho(ix,iy,iz) enddo enddo enddo write(0,*) "rho min,max=", dmin, dmax dmin= 1.0e+6 dmax= -dmin do iz= 1, npx do iy= 1, npy do ix = 1, npx if ( prs(ix,iy,iz) .lt. dmin ) dmin= prs(ix,iy,iz) if ( prs(ix,iy,iz) .gt. dmax ) dmax= prs(ix,iy,iz) enddo enddo enddo write(0,*) "prs min,max=", dmin, dmax dmin= 1.0e+6 dmax= -dmin do iz= 1, npx do iy= 1, npy do ix = 1, npx if ( uux(ix,iy,iz) .lt. dmin ) dmin= uux(ix,iy,iz) if ( uux(ix,iy,iz) .gt. dmax ) dmax= uux(ix,iy,iz) enddo enddo enddo write(0,*) "uux min,max=", dmin, dmax dmin= 1.0e+6 dmax= -dmin do iz= 1, npx do iy= 1, npy do ix = 1, npx if ( uuy(ix,iy,iz) .lt. dmin ) dmin= uuy(ix,iy,iz) if ( uuy(ix,iy,iz) .gt. dmax ) dmax= uuy(ix,iy,iz) enddo enddo enddo write(0,*) "uuy min,max=", dmin, dmax dmin= 1.0e+6 dmax= -dmin do iz= 1, npx do iy= 1, npy do ix = 1, npx if ( uuz(ix,iy,iz) .lt. dmin ) dmin= uuz(ix,iy,iz) if ( uuz(ix,iy,iz) .gt. dmax ) dmax= uuz(ix,iy,iz) enddo enddo enddo write(0,*) "uuz min,max=", dmin, dmax return 666 continue write(0,*) "can't read ", name(1:18), " at offset ",ioff," aborting." stop 667 continue write(0,*) "couldn't read",nbytes," bytes of ", & name(1:18)," at offset ",ioff," aborting with ", ierr stop end c ---------------------------------------------------------------- c This routine loads the supplied files which must be local. c as an initial state. See below for params (same as setupturbpencil) c subroutine loadturbpencil( pen0, pen1, tileOffset, & sumkntc, sumprs, dmin, dmax ) implicit none INCLUDE 'iq.h' INCLUDE 'thingftn.h' INCLUDE 'context.h' INCLUDE 'globals.h' integer pen0(3),pen1(3), tileOffset(3) REAL dmin(nvar), dmax(nvar) double precision sumkntc,sumprs integer ix,iy,iz,i real rho(npx,npy,npz),prs(npx,npy,npz),uux(npx,npy,npz) real uuy(npx,npy,npz),uuz(npx,npy,npz) common /rst/ rho,prs,uux,uuy,uuz integer nz0 integer mxstart,mxstop, mystart,mystop, mzstart,mzstop real zzz integer loadedblock(3) data loadedblock/-1,-1,-1/ c------- if ( tileOffset(1).ne. loadedblock(1) .or. & ( tileOffset(2).ne. loadedblock(2)) .or. & ( tileOffset(3).ne. loadedblock(3)) ) then call readblock(tileOffset(1),tileOffset(2),tileOffset(3)) loadedblock(1)= tileOffset(1) loadedblock(2)= tileOffset(2) loadedblock(3)= tileOffset(3) endif mxstart = pen0(1) mystart = pen0(2) mzstart = pen0(3) mxstop = pen1(1) mystop = pen1(2) mzstop = pen1(3) nz0= tileOffset(3) do iz= mzstart, mzstop zzz= zzl(iz+nz0) do iy= mystart, mystop do ix= mxstart, mxstop ddd0(1,ix,iy,iz) = rho(ix,iy,iz) ddd0(2,ix,iy,iz) = prs(ix,iy,iz) ddd0(3,ix,iy,iz) = uux(ix,iy,iz) ddd0(4,ix,iy,iz) = uuy(ix,iy,iz) ddd0(5,ix,iy,iz) = uuz(ix,iy,iz) ddd0(6,ix,iy,iz) = zzz sumkntc = sumkntc + ddd0(1,ix,iy,iz) * & ( ddd0(3,ix,iy,iz)**2 + ddd0(4,ix,iy,iz)**2 + ddd0(5,ix,iy,iz)**2 ) sumprs = sumprs + ddd0(2,ix,iy,iz) do i = 1, nvar if ( ddd0(i,ix,iy,iz) .lt. dmin(i) ) dmin(i) = ddd0(i,ix,iy,iz) if ( ddd0(i,ix,iy,iz) .gt. dmax(i) ) dmax(i) = ddd0(i,ix,iy,iz) enddo enddo enddo enddo return end c ---------------------------------------------------------------- c This routine is called once for each PENCIL, and uses the many c amplitude numbers set up by the initialization routine. c c There are (2kmax+1)*(2kmax+1)*(kmax+1) ==> 2601 for kmax=8. c c (All this is not only for cache performance, but to keep the pseudorandom c number sequence correct and repeatable.) c c tileOffset: Global coords for this tile c pen0,pen1 Limits in tile coords for this pencil c tileOffset Offset of this block in the global mesh (0,0,0) for llc c sumkntc,prs returned energy sums (accumulated w/input values) c dmin,dmax variable extrema returned (composite w/input) c subroutine setupturbpencil( pen0,pen1, & tileOffset, sumkntc,sumprs,dmin,dmax ) include 'iq.h' include 'globals.h' include 'context.h' REAL dmin(nvar), dmax(nvar) double precision sumkntc,sumprs integer pen0(3), pen1(3) integer tileOffset(3) c 2601 wavenumbers is enough for kmax=8 parameter (MX_WAVE=2604) REAL snn(MX_AXIS), css(MX_AXIS), aro(MX_WAVE), ape(MX_WAVE), & aror(MX_WAVE), aroi(MX_WAVE), aper(MX_WAVE), apei(MX_WAVE), & axr(MX_WAVE), ayr(MX_WAVE), azr(MX_WAVE), axi(MX_WAVE), & ayi(MX_WAVE), azi(MX_WAVE) common /setTurbcmn/ isetup, nlen, snn,css, & aro,ape,aror,aroi,aper,apei, & axr,ayr,azr, axi,ayi,azi c --- locals cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Set sin and cosine arrays c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c NOTE: setvec as it is (15 July 1993) works for any arrangment of c c sub-meshes (identical or not, cubical or not). c c However, the entire simulation mesh must be uniform and cubical. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc pi = 3.1415926535897932 pi2 = 2. * pi cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c SET CONSTANTS c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c grav = acceleration do to gravity in +r direction (grav = 0) c c beta = initial temperature gradient at top. c c flux = vertical heat flux everywhere in the initial state, c c it is also imposed at the bottom boundary at all times c c eitop = internal energy in center of zone just above top boundary c c deidr0 = vertical slope of internal energy imposed at bottom = 0 c c heatcv = heat capacity at constant volume = 1 c c rhoavg = mean mass density c c pavg = mean presure c c c0 = mean sound speed c c ei0 = mean internal energy c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c0 = 1. rhoavg = 1. pavg = (rhoavg * c0**2) / gamma c ei0 = pavg / ((gamma - 1.) * rhoavg) c eitop = ei0 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c This routine sets the presure, density, and velocity fields for c c homogeneius turbulence. Density and presure are constant. c c density = 1, speed of sound = 1. c c----------------------------------------------------------------------c c Inputs are: c c read (19,7) acard, dvamp, bvamp, tcamp c c read (19,9) acard,iseed,kmax,velamp,ratsc,wdthk,roamb,pgamb c c----------------------------------------------------------------------c c dvamp = Amplitude of Dynamic Viscousity c c bvamp = Amplitude of Bulk Viscousity c c tcamp = Amplitude of Thermal Conductivity c c iseed = Random number generator seed c c velamp = Root (Ensemble Mean) (Spatial Mean) Square velocity c c ratsc = (solinoidal power) / (compresible power) c c fkmax = Maximum mode excited c c wdthk = Width of excited modes c c roamb = fractional RMS desity fluctuations c c pgamb = fractional RMS presure fluctuations c c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Parameter c c Peak excited mode: fk0 = 0. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc fk0 = 0. c sea added ^^^^^^ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c SET mean values of field arrays c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Compute global pencil zone numbers: c cc nx0 = tileOffset(1) + pen0(1) - 1 cc ny0 = tileOffset(2) + pen0(2) - 1 cc nz0 = tileOffset(3) + pen0(3) - 1 nx0 = tileOffset(1) ny0 = tileOffset(2) nz0 = tileOffset(3) nxzons = nx nyzons = ny nzzons = nz mxstart = pen0(1) mystart = pen0(2) mzstart = pen0(3) mxstop = pen1(1) mystop = pen1(2) mzstop = pen1(3) do 10 iz = mzstart, mzstop zzz= zzl(iz+nz0) + 0.5 if ( zzz .gt. 1.0 ) zzz= zzz - 1.0 do 10 iy = mystart, mystop do 10 ix = mxstart, mxstop ddd0(1, ix,iy,iz ) = rhoavg ddd0(2, ix,iy,iz ) = pavg ddd0(3, ix,iy,iz ) = 0 ddd0(4, ix,iy,iz ) = 0 ddd0(5, ix,iy,iz ) = 0 ddd0(6, ix,iy,iz ) = zzz 10 continue if(kmax .eq. 0) then xlen1 = xxl(nxzons+1) - xxl(1) ylen1 = yyl(nyzons+1) - yyl(1) zlen1 = zzl(nzzons+1) - zzl(1) facxx = pi2 / xlen1 facyy = pi2 / ylen1 faczz = pi2 / zlen1 do 51 k= mzstart, mzstop zzz = zzl(k+nz0) do 51 j= mystart, mystop yyy = yyl(j+ny0) do 51 i= mxstart, mxstop xxx = xxl(i+nx0) arg = facxx * xxx + facyy * yyy + faczz * zzz vvvv = velamp * sin(arg) ddd0(3,i,j,k) = vvvv ddd0(4,i,j,k) = vvvv ddd0(5,i,j,k) = vvvv 51 continue endif c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Run through all of the wave vectors "k" c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c // wavenumber: iw = 0 kmax2 = kmax**2 neglim = kmax * nlen do 900 kz = -kmax, kmax do 800 ky = -kmax, kmax do 700 kx = 0, kmax cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c the following is for debug of 2-period waves only c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c do 900 kz = -kmax, kmax,2 c do 800 ky = -kmax, kmax,2 c do 700 kx = 0, kmax,2 c kx = mod(iseed , 2) c ky = mod(iseed/2, 2) c kz = mod(iseed/4, 2) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc iw = iw + 1 if ( iw.gt.MX_WAVE ) then print *, MX_WAVE,' wave number exceeded' stop endif if(kx**2+ky**2+kz**2 .gt. kmax2) go to 650 if((kx .eq. 0) .and. (ky .lt. 0)) go to 650 if((kx .eq. 0) .and. (ky .eq. 0) .and. (kz .le. 0)) go to 650 do 600 k = mzstart, mzstop kmod = 1 + mod((k+nz0)*kz+neglim, nlen) snk = snn(kmod) csk = css(kmod) do 600 j = mystart, mystop jmod = 1 + mod((j+ny0)*ky+neglim, nlen) snj = snn(jmod) csj = css(jmod) snjk = snj*csk + csj*snk csjk = csj*csk - snj*snk do 600 i = mxstart, mxstop imod = 1 + mod((i+nx0)*kx+neglim, nlen) sni = snn(imod) csi = css(imod) snijk = sni*csjk + csi*snjk csijk = csi*csjk - sni*snjk ddd0(1,i,j,k) = ddd0(1,i,j,k) + aror(iw) * csijk + aroi(iw) * snijk ddd0(2,i,j,k) = ddd0(2,i,j,k) + aper(iw) * csijk + apei(iw) * snijk ddd0(3,i,j,k) = ddd0(3,i,j,k) + axr(iw) * csijk + axi(iw) * snijk ddd0(4,i,j,k) = ddd0(4,i,j,k) + ayr(iw) * csijk + ayi(iw) * snijk ddd0(5,i,j,k) = ddd0(5,i,j,k) + azr(iw) * csijk + azi(iw) * snijk 600 continue 650 continue 700 continue 800 continue c write (6,*) 'kz,ro2,ac2,as2 = ', kz, ro2tot, ac2tot, as2tot 900 continue c if ( ithread.eq.1 ) then c write (6,*) 'Vel. squared: Comp. Soln: ', 0.5*ac2tot, 0.5*as2tot c write (6,*) 'Vx2,Vy2,Vz2 = ', ax2tot, ay2tot, az2tot c rmsdro = sqrt(0.5 * ro2tot) c rmsdpe = sqrt(0.5 * pe2tot) c write (6,*) 'RMS(ro-ro0), RMS(p-pavg) = ', rmsdro, rmsdpe c endif c Sweep through one more time collecting statistics c do iz = mzstart, mzstop do iy = mystart, mystop do ix = mxstart, mxstop sumkntc = sumkntc + ddd0(1,ix,iy,iz) * & ( ddd0(3,ix,iy,iz)**2 + ddd0(4,ix,iy,iz)**2 + ddd0(5,ix,iy,iz)**2 ) sumprs = sumprs + ddd0(2,ix,iy,iz) do i = 1, nvar if ( ddd0(i,ix,iy,iz) .lt. dmin(i) ) dmin(i) = ddd0(i,ix,iy,iz) if ( ddd0(i,ix,iy,iz) .gt. dmax(i) ) dmax(i) = ddd0(i,ix,iy,iz) enddo enddo enddo enddo return end c---------------------------------------------------------------------72 c c THIS FUNCTION CALCULATES THE INVERSE ERF(P) c FUNCTION ERFINV(P) c ************************************************************************** c * THIS CODE IS D36PPM3D * c * THIS CODE WAS WRITTEN BY: DAVID PORTER ALL RIGHTS RESERVED * c ************************************************************************** EPS=.00001 SP=0.8862269255 z=0. pp=p C 5 Z=Z+(PP-ERF1(Z))*SP*EXP(Z*Z) IF((ABS(PP-ERF1(Z))*SP*EXP(Z*Z)) .GT. EPS) GO TO 5 erfinv=z C RETURN END c---------------------------------------------------------------------72 c function ranfmy(dumb) c ************************************************************************** c * THIS CODE IS D36PPM3D * c * THIS CODE WAS WRITTEN BY: DAVID PORTER ALL RIGHTS RESERVED * c ************************************************************************** integer*4 i0, i1, i2 common /ran_sav/ i0, i1, i2 i0 = mod(i0*i1, i2) ranfmy = float(i0) / float(i2) return end c---------------------------------------------------------------- c subroutine setrmy(iseed) c ************************************************************************** c * THIS CODE IS D36PPM3D * c * THIS CODE WAS WRITTEN BY: DAVID PORTER ALL RIGHTS RESERVED * c ************************************************************************** integer*4 i0, i1, i2 common /ran_sav/ i0, i1, i2 i0 = iseed i1 = 11447 i2 = 8803 return end c ---------------------------------------------------------------- c FUNCTION ERF1(X1) c ************************************************************************** c * THIS CODE IS D36PPM3D * c * THIS CODE WAS WRITTEN BY: DAVID PORTER ALL RIGHTS RESERVED * c ************************************************************************** A1=0.0705230784 A2=0.0422820123 A3=0.0092705272 A4=0.0001520143 A5=0.0002765672 A6=0.0000430638 S=1. IF(X1 .LT. 0.) S=-1. X=ABS(X1) Y=1.+X*(A1+X*(A2+X*(A3+X*(A4+X*(A5+X*A6))))) ERF1=S*(1.-1./(Y**16)) RETURN END c---------------------------------------------------------------- c c subroutine printtrace( ncycle, stime, dtime, c & courmx, sumkntc,sumprs, dvol,gamma ) c REAL courmx, dtime, stime c double precision sumkntc,sumprs, dvol,gamma c logical virgin c save virgin c data virgin/.true./ c c if ( virgin ) then c write(6,100) c 100 format(" Step",t8,"Time",t22,"dTime",t35,"Courant",t44,"Ek", c & t59,"Eh",t75,"Etotal") c virgin = .false. c open(unit=36,file='trace',form='formatted',access='append') c write(36,100) c endif c c ek = dvol * 0.5*sumkntc c eh = dvol * sumprs / (gamma-1.0) c et = ek + eh c c write(6,101) ncycle,stime,dtime, courmx, ek,eh,et c write(36,101) ncycle,stime,dtime, courmx, ek,eh,et cc call flush(36) c c 101 format( i6,1p,e13.6,1x,e13.6,1x,0p,f7.5,1x,1p,3(e15.8)) c c return c end