! This code generated by rnpl, a numerical programming language ! copyright (c) 1994,1995 by Robert L. Marsa and Matthew W. Choptuik c----------------------------------------------------------------------- c c M E M O R Y M A N A G E M E N T c c R O U T I N E S c c----------------------------------------------------------------------- c These routines, in conjunction with the file 'mm.inc' reproduced c here: c c integer maxblk c parameter ( maxblk = 1 000 ) c c integer mmmloc(0:maxblk), c * mmsize(0:maxblk), c * mmnext(0:maxblk) c logical mfree(0:maxblk) c c common / commm / mmmloc, mmsize, mmnext, mfree c c constitute a FORTRAN implementation of the ``First-fit'' c memory management scheme described in c c D. Knuth c The Art of Computer Programming c Vol. 1, Fundamental Algorithms c Section 2.5: Dynamic Storage Allocation c c Invoking routines are responsible for declaring actual storage c array (arena) from which blocks are allocated/deallocated. c c Example arena declaration: c c integer memsiz c parameter ( memsiz = 1 000 000 ) c real*8 q(memsiz) c c----------------------------------------------------------------------- c c Sample usage (given above declarations) c c Initialize memory management structure ... c c integer mmaloc, mmdeal c integer ptr1, ptr2 c c call mmini(memsiz) c c Allocate two blocks of storage ... c c ptr1 = mmaloc(1000) c if( ptr1 .lt. 0 ) then c write(*,*) 'Unable to allocate block of size 1000' c stop c end if c ptr2 = mmaloc(1000) c if( ptr2 .lt. 0 ) then c write(*,*) 'Unable to allocate block of size 1000' c stop c end if c c ... do something with the blocks ... c c call subroutine_which_wants_2_vectors(q(ptr1),q(ptr2),1000) c c ... and free them. Note that block size must be supplied to c deallocating routine. c c ptr1 = mmdeal(ptr1,1000) c ptr2 = mmdeal(ptr2,1000) c c c----------------------------------------------------------------------- c======================================================================= c c subroutine mmini(memsiz) c c I: memsiz integer Arena size c c O: None c c c Initializes memory management structure. c c======================================================================= subroutine mmini(memsiz) implicit none C----------------------------------------------------------------------- C C M E M O R Y M A N A G E M E N T C C S T R U C T U R E C C----------------------------------------------------------------------- C integer maxblk parameter ( maxblk = 1 000 ) C integer mmmloc(0:maxblk), * mmsize(0:maxblk), * mmnext(0:maxblk) logical mfree(0:maxblk) C common / commm / mmmloc, mmsize, mmnext, mfree C C---------------------------------------------------------------------- integer memsiz call ivls(mmmloc(0),0,maxblk+1) call ivls(mmsize(0),0,maxblk+1) call ivls(mmnext(0),0,maxblk+1) mmnext(0) = 1 mmmloc(1) = 1 mmsize(1) = memsiz return end c======================================================================= c c integer function mmfree() c c I: None c c O: None c c Returns first available location in list of memory block nodes, c -1 if none available. c c======================================================================= integer function mmfree() implicit none C----------------------------------------------------------------------- C C M E M O R Y M A N A G E M E N T C C S T R U C T U R E C C----------------------------------------------------------------------- C integer maxblk parameter ( maxblk = 1 000 ) C integer mmmloc(0:maxblk), * mmsize(0:maxblk), * mmnext(0:maxblk) logical mfree(0:maxblk) C common / commm / mmmloc, mmsize, mmnext, mfree C C---------------------------------------------------------------------- integer p, q do mmfree = 0 , maxblk mfree(mmfree) = .true. end do q = 0 200 continue p = mmnext(q) if( p .eq. 0 ) go to 300 mfree(p) = .false. q = p go to 200 300 continue do mmfree = 1 , maxblk if( mfree(mmfree) ) return end do mmfree = -1 return end c======================================================================= c c integer function mmaloc(reqsiz) c c I: reqsiz integer Requested block size. c c O: None c c c Allocates memory block of size REQSIZ, returning pointer to first c location, or -1 if unable to perform allocation. c c======================================================================= integer function mmaloc(reqsiz) implicit none C----------------------------------------------------------------------- C C M E M O R Y M A N A G E M E N T C C S T R U C T U R E C C----------------------------------------------------------------------- C integer maxblk parameter ( maxblk = 1 000 ) C integer mmmloc(0:maxblk), * mmsize(0:maxblk), * mmnext(0:maxblk) logical mfree(0:maxblk) C common / commm / mmmloc, mmsize, mmnext, mfree C C---------------------------------------------------------------------- integer k, p, q, reqsiz q = 0 100 continue p = mmnext(q) if( p .eq. 0 ) then mmaloc = -1 return end if if( mmsize(p) .ge. reqsiz ) then k = mmsize(p) - reqsiz if( k .eq. 0 ) then mmnext(q) = mmnext(p) else mmsize(p) = k end if mmaloc = mmmloc(p) + k return end if q = mmnext(q) go to 100 end c======================================================================= c c integer function mmdeal(p0loc,p0size) c c I: p0loc integer Pointer to block to be freed. c p0size integer Size of block to be freed. c c O: None c c c Returns block of storage of length P0SIZE, which begins at c location P0LOC to free storage list. c c======================================================================= integer function mmdeal(p0loc,p0size) implicit none C----------------------------------------------------------------------- C C M E M O R Y M A N A G E M E N T C C S T R U C T U R E C C----------------------------------------------------------------------- C integer maxblk parameter ( maxblk = 1 000 ) C integer mmmloc(0:maxblk), * mmsize(0:maxblk), * mmnext(0:maxblk) logical mfree(0:maxblk) C common / commm / mmmloc, mmsize, mmnext, mfree C C---------------------------------------------------------------------- integer mmfree integer p, p0, p0loc, p0size, q p0 = mmfree() if( p0 .lt. 0 ) then mmdeal = p0 return else mmmloc(p0) = p0loc mmsize(p0) = p0size end if q = 0 100 continue p = mmnext(q) if( mmmloc(p) .gt. mmmloc(p0) .or. p .eq. 0 ) then if( mmmloc(p0) + mmsize(p0) .eq. mmmloc(p) ) then mmnext(p0) = mmnext(p) mmsize(p0) = mmsize(p0) + mmsize(p) else mmnext(p0) = p end if if( mmmloc(q) + mmsize(q) .eq. mmmloc(p0) ) then mmnext(q) = mmnext(p0) mmsize(q) = mmsize(q) + mmsize(p0) else mmnext(q) = p0 end if mmdeal = 0 return end if q = p go to 100 end c======================================================================= c c subroutine mmdump(unit) c c I: unit integer Logical unit for dump. c c O: None c c c Dumps map of free list on UNIT. c c======================================================================= subroutine mmdump(unit) implicit none C----------------------------------------------------------------------- C C M E M O R Y M A N A G E M E N T C C S T R U C T U R E C C----------------------------------------------------------------------- C integer maxblk parameter ( maxblk = 1 000 ) C integer mmmloc(0:maxblk), * mmsize(0:maxblk), * mmnext(0:maxblk) logical mfree(0:maxblk) C common / commm / mmmloc, mmsize, mmnext, mfree C C---------------------------------------------------------------------- integer n, p, unit write(unit,100) 100 FORMAT(//T10,' <<< Free Block Structure. >>>'// * T11,'Block', T29,'Locations'/) p = 0 n = 0 200 continue p = mmnext(p) if( p .eq. 0 ) then return else n = n + 1 write(unit,300) n, mmmloc(p), mmmloc(p) + mmsize(p) - 1 300 FORMAT(T12,I3,T25,I7,' ... ',I8) end if go to 200 end !----------------------------------------------------------------------- ! Front ends to latinfo for extracting: ! ! getlatid(ilat): Lattice id of ilat'th lattice ! getlatcs(ilat): Coordinate system of ilat'th lattice (0..) ! getrank (ilat): Rank of ilat'th lat ! getpshape(ilat): Pointer into PTRS for shape of lat ! ptrs(getpshape(ilat) - 1 + j) := shape(j) ! getshape(ilat,j): j'th dim (element of shape) of ilat'th lat ! ptrs(getpshape(ilat) - 1 + j) := shape(j) ! getpbounds(ilat): Pointer into PTRS for bounds of lat ! ptrs(getpbounds(ilat) - 1 + j) := bounds(j) ! getbounds(ilat,j): j'th component of bounds ! getpcoord(ilat,j) Pointer into Q for j'th coords. of ilat'th lat ! getsize(ilat): Size of lattice (product reduction of shape) !----------------------------------------------------------------------- integer function getnlat() implicit none integer latinfo getnlat = latinfo(0,0,0) return end integer function getlatcs(ilat) implicit none integer latinfo integer ilat getlatcs = latinfo(ilat,0,9) return end integer function getlatid(ilat) implicit none integer latinfo integer ilat getlatid = latinfo(ilat,0,1) return end integer function getrank(ilat) implicit none integer latinfo integer ilat getrank = latinfo(ilat,0,2) return end integer function getpshape(ilat) implicit none integer latinfo integer ilat getpshape = latinfo(ilat,0,3) return end integer function getshape(ilat,j) implicit none integer latinfo integer ilat, j getshape = latinfo(ilat,j,4) return end integer function getpbounds(ilat) implicit none integer latinfo integer ilat getpbounds = latinfo(ilat,0,7) return end integer function getbounds(ilat,j) implicit none integer latinfo integer ilat, j getbounds = latinfo(ilat,j,8) return end integer function getpcoord(ilat,j) implicit none integer latinfo integer ilat, j getpcoord = latinfo(ilat,j,5) return end integer function getsize(ilat) implicit none integer latinfo integer ilat getsize = latinfo(ilat,0,6) return end !----------------------------------------------------------------------- ! Main decoding routine ! ! Requires 'ivprod' for product reduction of integer vector. ! ! option = 0: Number of lattices ! option = 1: Lattice id ! option = 2: Rank of ilat'th lattice ! option = 3: Pointer into ptrs for shape of ilat'th lattice ! option = 4: j'th dim of ilat'th lattice ! option = 5: Pointer into q for coordinates of j'th dim of ! ilat'th lattice ! option = 6: Size of lattice (product reduction of shape) ! option = 7: Pointer into ptrs for bounds of ilat'th lattice ! option = 8: j'th bound of ilat'th lattice ! option = 9: Coordinate system !----------------------------------------------------------------------- integer function latinfo(ilat,j,option) implicit none include 'sys_param.inc' !----------- from gfuni0.inc ---------------------------------- integer getnlat, getlatid, getrank, getpshape, & getshape, getpcoord, getsize, getpbounds, & getbounds, getlatcs external getnlat, getlatid, getrank, getpshape, & getshape, getpcoord, getsize, getpbounds, & getbounds, getlatcs real*8 q(memsiz) common / com_gfuni0 / q integer maxrank parameter ( maxrank = 5 ) integer MAXSTEP parameter ( MAXSTEP = 50 ) integer sizelatcb parameter ( sizelatcb = 3 + 4 * maxrank ) integer maxgfcn parameter ( maxgfcn = 100 ) integer ngfcn common / com_gfuni0 / ngfcn integer maxptrs parameter ( maxptrs = maxgfcn * (2 + sizelatcb) ) integer ptrs(maxptrs) common / com_gfuni0 / ptrs integer nptrs common / com_gfuni0 / nptrs integer ivprod integer ilat, j, option, cs integer nlat, base, latid, rank integer shape, bounds, cptr nlat = (nptrs - 2*ngfcn) / sizelatcb if( option .eq. 0 ) then latinfo = nlat return end if if( 1 .le. ilat .and. ilat .le. nlat ) then base = 2*ngfcn + (ilat - 1)*sizelatcb + 1 latid = ptrs(base) cs = ptrs(base+1) rank = ptrs(base+2) shape = base + 3 bounds = shape + maxrank cptr = bounds + 2*maxrank if( option .eq. 1 ) then latinfo = latid else if( option .eq. 2 ) then latinfo = rank else if( option .eq. 3 ) then latinfo = shape else if( option .eq. 4 ) then if( 1 .le. j .and. j .le. rank ) then latinfo = ptrs(shape + j - 1) else write(0,1000) j, rank, ilat 1000 format(' latinfo: j = ',i2,' out of bounds for ', & 'rank ',i2,' lattice ',i2) latinfo = -1 end if else if( option .eq. 5 ) then if( 1 .le. j .and. j .le. rank ) then latinfo = ptrs(cptr + j-1) else write(0,1000) j, rank, ilat latinfo = -1 end if else if( option .eq. 6 ) then latinfo = ivprod(ptrs(shape),rank) else if( option .eq. 7 ) then latinfo = bounds else if( option .eq. 8 ) then if( 1 .le. j .and. j .le. 2*rank ) then latinfo = ptrs(bounds + j - 1) else write(0,1000) j, rank, ilat latinfo = -1 end if else if( option .eq. 9 ) then latinfo = cs end if else write(0,1100) ilat, nlat 1100 format(' latinfo: ilat = ',i2,' out of bounds 1 - ',i2) latinfo = -1 end if return end subroutine res_Phi(Phi_res,Phi_np1,Phi_n,Phi_nm1,r_n,g1_shp,g1_bds & ,rho,z,dt,drho,dz) implicit none include 'globals.inc' integer g1_shp(2) integer g1_bds(4) real*8 Phi_res(g1_bds(1):g1_bds(2),g1_bds(3):g1_bds(4)) real*8 Phi_np1(g1_bds(1):g1_bds(2),g1_bds(3):g1_bds(4)) real*8 Phi_n(g1_bds(1):g1_bds(2),g1_bds(3):g1_bds(4)) real*8 Phi_nm1(g1_bds(1):g1_bds(2),g1_bds(3):g1_bds(4)) real*8 r_n(g1_bds(1):g1_bds(2),g1_bds(3):g1_bds(4)) real*8 rho(*) real*8 z(*) real*8 dt real*8 drho real*8 dz integer i integer j integer Nrho integer Nz Nrho = Nrho0 * 2**level + 1 Nz = Nz0 * 2**level + 1 i=1 do j=1, Nz, 1 Phi_res(i,j)=+Phi_np1(i,j)-4*Phi_np1(i+1,j)/3+Phi_np1(i+2,j)/3 end do i=Nrho do j=2, Nz-1, 1 Phi_res(i,j)=r_n(i,j)*(3*Phi_np1(i,j)-4*Phi_n(i,j)+Phi_nm1(i,j & ))/(2*dt)+rho(i)*(3*Phi_np1(i,j)-4*Phi_np1(i-1,j)+Phi_np1( & i-2,j))/(2*drho)+z(j)*(Phi_np1(i,j+1)-Phi_np1(i,j-1))/(2*dz) end do do i=2, Nrho-1, 1 do j=2, Nz-1, 1 Phi_res(i,j)=(Phi_np1(i,j)-2*Phi_n(i,j)+Phi_nm1(i,j))/dt**2-(( & Phi_n(i,j+1)-2*Phi_n(i,j)+Phi_n(i,j-1))/dz**2+(Phi_n(i+1,j)- & 2*Phi_n(i,j)+Phi_n(i-1,j))/drho**2+(Phi_n(i+1,j)-Phi_n(i-1,j & ))/(2*drho)/rho(i)) end do end do do i=2, Nrho-1, 1 j=1 Phi_res(i,j)=r_n(i,j)*(3*Phi_np1(i,j)-4*Phi_n(i,j)+Phi_nm1(i,j & ))/(2*dt)+rho(i)*(Phi_np1(i+1,j)-Phi_np1(i-1,j))/(2*drho)+ & z(j)*((-3)*Phi_np1(i,j)+4*Phi_np1(i,j+1)-Phi_np1(i,j+2))/(2* & dz) end do do i=2, Nrho-1, 1 j=Nz Phi_res(i,j)=r_n(i,j)*(3*Phi_np1(i,j)-4*Phi_n(i,j)+Phi_nm1(i,j & ))/(2*dt)+rho(i)*(Phi_np1(i+1,j)-Phi_np1(i-1,j))/(2*drho)+ & z(j)*(3*Phi_np1(i,j)-4*Phi_np1(i,j-1)+Phi_np1(i,j-2))/(2*dz) end do return end !---------------------------------------------------------------------- ! swap_top !---------------------------------------------------------------------- subroutine swap_top() implicit none include 'sys_param.inc' include 'gfuni0.inc' include 'globals.inc' integer tmp tmp=ptrs(gf_st(1)+1) ptrs(gf_st(1)+1)=ptrs(gf_st(1)+2) ptrs(gf_st(1)+2)=ptrs(gf_st(1)+3) ptrs(gf_st(1)+3)=tmp return end !---------------------------------------------------------------------- ! initial_guess !---------------------------------------------------------------------- subroutine initial_guess() implicit none include 'sys_param.inc' include 'gfuni0.inc' include 'other_glbs.inc' include 'globals.inc' integer i,j,k call rdvcpy(q(ptrs(gf_st(1)+1)), q(ptrs(gf_st(1)+2)),getsize(1)) return end !---------------------------------------------------------------------- ! calc_resid !---------------------------------------------------------------------- real*8 function calc_resid() implicit none include 'sys_param.inc' include 'gfuni0.inc' include 'globals.inc' real*8 u,l2norm real*8 u0,v0 u0=l2norm(getsize(1),q(ptrs(gf_st(3)+1))) v0=l2norm(getsize(1),q(ptrs(gf_st(1)+1))) if (v0 .ne. 0.0) then u0 = u0/v0 end if u=u0 calc_resid=u return end !---------------------------------------------------------------------- ! take_step !---------------------------------------------------------------------- subroutine take_step(t) implicit none include 'sys_param.inc' include 'gfuni0.inc' include 'other_glbs.inc' include 'globals.inc' integer mmaloc, mmdeal real*8 t call update0(q(ptrs(gf_st(1)+1)),q(ptrs(gf_st(1)+2)), & q(ptrs(gf_st(1)+3)),q(ptrs(gf_st(2)+1)),ptrs(getpshape(1)), & ptrs(getpbounds(1)),q(rho),q(z),dt,drho,dz) call res_Phi(q(ptrs(gf_st(3)+1)),q(ptrs(gf_st(1)+1)), & q(ptrs(gf_st(1)+2)),q(ptrs(gf_st(1)+3)),q(ptrs(gf_st(2)+1)), & ptrs(getpshape(1)),ptrs(getpbounds(1)),q(rho),q(z),dt,drho,dz) return end !---------------------------------------------------------------------- ! iter_gfuncs !---------------------------------------------------------------------- integer function iter_gfuncs(tt) implicit none real*8 tt include 'sys_param.inc' include 'gfuni0.inc' include 'other_glbs.inc' include 'globals.inc' real*8 t,u,d,calc_resid integer i,j,k,step u=0.0 t=tt d=dt dt=-dt/2 call rdvcpy(q(ptrs(gf_st(1)+2)),q(ptrs(gf_st(1)+3)),getsize(1)) call rdvcpy(q(ptrs(gf_st(1)+1)),q(ptrs(gf_st(1)+2)),getsize(1)) step=0 100 continue step=step+1 call rdvaddc(q(ptrs(gf_st(1)+2)),(0.5D0),q(ptrs(gf_st(1)+1)), & (0.5D0),q(ptrs(gf_st(1)+3)),getsize(1)) call take_step(t) u=calc_resid() if(u .gt. epsiter .AND. step .lt. MAXSTEP) then goto 100 end if call swap_top() if(step .eq. MAXSTEP) then write (*,*) 'Initialization did not converge within ', & MAXSTEP,' iterations.' end if dt=d iter_gfuncs = step return end !---------------------------------------------------------------------- ! init_params_attribs !---------------------------------------------------------------------- subroutine init_params_attribs() implicit none include 'sys_param.inc' include 'gfuni0.inc' include 'other_glbs.inc' include 'globals.inc' epsdis=(0D0) epsiter=(1D-05) fout=0 idstr='gts' iter=100 lambda=(0.5D0) level=0 Nrho0=3 Nz0=3 output(1)=1 output(2)=-1 output(3)=-1 output(4)=1 output(5)=0 rho0=(1D0) ser=0 start_t=(0D0) s_step=0 tag=' ' z0=(1D0) out_gf(1)=0 out_gf(2)=0 out_gf(3)=0 return end !---------------------------------------------------------------------- ! read_parameters !---------------------------------------------------------------------- subroutine read_parameters(p_file) implicit none character*(*) p_file include 'sys_param.inc' include 'gfuni0.inc' include 'other_glbs.inc' include 'globals.inc' call get_real_param(p_file,'amp',0,amp,1) call get_real_param(p_file,'delta',0,delta,1) call get_real_param(p_file,'epsdis',1,epsdis,1) call get_real_param(p_file,'epsiter',1,epsiter,1) call get_int_param(p_file,'fout',1,fout,1) call get_str_param(p_file,'idstr',1,idstr,1) call get_str_param(p_file,'in_file',0,in_file,1) call get_int_param(p_file,'iter',1,iter,1) call get_real_param(p_file,'lambda',1,lambda,1) call get_int_param(p_file,'level',1,level,1) call get_int_param(p_file,'Nrho0',1,Nrho0,1) call get_int_param(p_file,'Nz0',1,Nz0,1) call get_ivec_param(p_file,'output',1,output,5) call get_str_param(p_file,'out_file',0,out_file,1) call get_real_param(p_file,'r0',0,r0,1) call get_real_param(p_file,'rho0',1,rho0,1) call get_real_param(p_file,'rhomax',0,rhomax,1) call get_real_param(p_file,'rhomin',0,rhomin,1) call get_int_param(p_file,'ser',1,ser,1) call get_real_param(p_file,'start_t',1,start_t,1) call get_int_param(p_file,'s_step',1,s_step,1) call get_str_param(p_file,'tag',1,tag,1) call get_real_param(p_file,'z0',1,z0,1) call get_real_param(p_file,'zmax',0,zmax,1) call get_real_param(p_file,'zmin',0,zmin,1) return end !---------------------------------------------------------------------- ! read_attributes !---------------------------------------------------------------------- subroutine read_attributes(p_file) implicit none character*(*) p_file integer encode include 'sys_param.inc' include 'gfuni0.inc' include 'globals.inc' call get_int_param(p_file,'out_gf',1,out_gf,3) nout_gf=encode(lout_gf,out_gf,3) return end !---------------------------------------------------------------------- ! init_coord_difs !---------------------------------------------------------------------- subroutine init_coord_difs() implicit none include 'sys_param.inc' include 'gfuni0.inc' include 'other_glbs.inc' include 'globals.inc' 10 format (a,a,i5,a) 11 format (a,i5,a) if(level .gt. 0 .and. mod(Nrho0,2) .eq. 1) then write(*,11)'WARNING: Nrho0=',Nrho0,' is odd and level != 0' write(*,11)'WARNING: Adjusting Nrho0 to ',Nrho0-1, & '. If this is part of a convergence' write(*,10)'WARNING: test, rerun level 0 calculation with' & ,' Nrho0=',Nrho0-1,'.' Nrho0 = Nrho0 - 1 end if Nrho = Nrho0 * 2**level + 1 if(level .gt. 0 .and. mod(Nz0,2) .eq. 1) then write(*,11)'WARNING: Nz0=',Nz0,' is odd and level != 0' write(*,11)'WARNING: Adjusting Nz0 to ',Nz0-1, & '. If this is part of a convergence' write(*,10)'WARNING: test, rerun level 0 calculation with' & ,' Nz0=',Nz0-1,'.' Nz0 = Nz0 - 1 end if Nz = Nz0 * 2**level + 1 drho=(rhomax - rhomin)/(Nrho - 1) dz=(zmax - zmin)/(Nz - 1) dt=lambda*sqrt((drho*drho+dz*dz)/2) return end !---------------------------------------------------------------------- ! init_lats !---------------------------------------------------------------------- subroutine init_lats() implicit none include 'sys_param.inc' include 'gfuni0.inc' include 'other_glbs.inc' include 'globals.inc' integer i,mmaloc cname(1)='t' cname(2)='rho' cname(3)='z' ncnames = 6 ngfcn = 5 nptrs = sizelatcb * 1 + 2 * ngfcn ptrs(11)=1 ptrs(12)=0 ptrs(13)=2 ptrs(14)=Nrho-1+1 ptrs(15)=Nz-1+1 ptrs(16)=1 ptrs(17)=1 ptrs(18)=1 ptrs(19)=1 ptrs(20)=Nrho ptrs(21)=1 ptrs(22)=Nz ptrs(23)=1 ptrs(24)=1 ptrs(25)=1 ptrs(26)=1 ptrs(27)=1 ptrs(28)=1 rho=mmaloc(getsize(1)) call rdvramp(q(rho),getshape(1,1),rhomin,drho) ptrs(29)=rho z=mmaloc(getsize(1)) call rdvramp(q(z),getshape(1,2),zmin,dz) ptrs(30)=z return end !---------------------------------------------------------------------- ! init_gfuncs !---------------------------------------------------------------------- subroutine init_gfuncs() implicit none integer mmaloc character*(128) temp character*256 catsqz include 'sys_param.inc' include 'gfuni0.inc' include 'globals.inc' 901 format (a,i1) 902 format (a,i1,a) ptrs(ngfcn + 1) = 2*ngfcn + 0*sizelatcb + 1 ptrs(1) = mmaloc(getsize(1)) ptrs(ngfcn + 2) = 2*ngfcn + 0*sizelatcb + 1 ptrs(2) = mmaloc(getsize(1)) ptrs(ngfcn + 3) = 2*ngfcn + 0*sizelatcb + 1 ptrs(3) = mmaloc(getsize(1)) write(temp,901) 'Phi_',level attrfname(1) = catsqz(tag,temp) gf_ln(1)=3 gf_st(1)=0 ptrs(ngfcn + 4) = 2*ngfcn + 0*sizelatcb + 1 ptrs(4) = mmaloc(getsize(1)) write(temp,901) 'r_',level attrfname(2) = catsqz(tag,temp) gf_ln(2)=1 gf_st(2)=3 ptrs(ngfcn + 5) = 2*ngfcn + 0*sizelatcb + 1 ptrs(5) = mmaloc(getsize(1)) write(temp,901) 'Phi_res_',level attrfname(3) = catsqz(tag,temp) gf_ln(3)=1 gf_st(3)=4 return end !---------------------------------------------------------------------- ! gen_gfuncs !---------------------------------------------------------------------- subroutine gen_gfuncs(t) implicit none include 'sys_param.inc' include 'gfuni0.inc' include 'other_glbs.inc' include 'globals.inc' real*8 t integer i,j,k do i=1, Nrho, 1 do j=1, Nz, 1 q(ptrs(gf_st(2)+gf_ln(2))+(i-1)+(j-1)*getshape(1,1))=sqrt( & q(rho+i-1)**2+q(z+j-1)**2) end do end do do i=1, Nrho, 1 do j=1, Nz, 1 q(ptrs(gf_st(1)+gf_ln(1))+(i-1)+(j-1)*getshape(1,1))=amp*exp(- & ((sqrt(q(rho+i-1)**2+q(z+j-1)**2)-r0)/delta)**2) end do end do return end !---------------------------------------------------------------------- ! write_attribs !---------------------------------------------------------------------- subroutine write_attribs(f_name) implicit none character*64 f_name include 'sys_param.inc' include 'gfuni0.inc' include 'globals.inc' include 'other_glbs.inc' integer res, gft_write_id_int_p integer gft_write_id_float_p, gft_write_id_str_p res=gft_write_id_float_p(f_name,'amp',amp,1) res=gft_write_id_float_p(f_name,'delta',delta,1) res=gft_write_id_float_p(f_name,'epsdis',epsdis,1) res=gft_write_id_float_p(f_name,'epsiter',epsiter,1) res=gft_write_id_int_p(f_name,'fout',fout,1) res=gft_write_id_str_p(f_name,'idstr',idstr,1) res=gft_write_id_str_p(f_name,'in_file',in_file,1) res=gft_write_id_int_p(f_name,'iter',iter,1) res=gft_write_id_float_p(f_name,'lambda',lambda,1) res=gft_write_id_int_p(f_name,'level',level,1) res=gft_write_id_int_p(f_name,'Nrho0',Nrho0,1) res=gft_write_id_int_p(f_name,'Nz0',Nz0,1) res=gft_write_id_int_p(f_name,'output',output,5) res=gft_write_id_str_p(f_name,'out_file',out_file,1) res=gft_write_id_float_p(f_name,'r0',r0,1) res=gft_write_id_float_p(f_name,'rho0',rho0,1) res=gft_write_id_float_p(f_name,'rhomax',rhomax,1) res=gft_write_id_float_p(f_name,'rhomin',rhomin,1) res=gft_write_id_int_p(f_name,'ser',ser,1) res=gft_write_id_float_p(f_name,'start_t',start_t,1) res=gft_write_id_int_p(f_name,'s_step',s_step,1) res=gft_write_id_str_p(f_name,'tag',tag,1) res=gft_write_id_float_p(f_name,'z0',z0,1) res=gft_write_id_float_p(f_name,'zmax',zmax,1) res=gft_write_id_float_p(f_name,'zmin',zmin,1) return end !---------------------------------------------------------------------- ! dump_state !---------------------------------------------------------------------- subroutine dump_state(t, st, f_name) implicit none include 'sys_param.inc' include 'gfuni0.inc' include 'globals.inc' real*8 t integer st, res character*64 f_name integer gft_close, gft_write_idata res=gft_write_idata(f_name,'Phi[1]',ptrs(getpshape(1)), & getrank(1),q(ptrs(gf_st(1)+2))) res=gft_write_idata(f_name,'Phi[0]',ptrs(getpshape(1)), & getrank(1),q(ptrs(gf_st(1)+3))) res=gft_write_idata(f_name,'r[0]',ptrs(getpshape(1)), & getrank(1),q(ptrs(gf_st(2)+1))) start_t=t s_step=st call write_attribs(f_name) res=gft_close(f_name) return end !---------------------------------------------------------------------- ! cleanup !---------------------------------------------------------------------- subroutine cleanup implicit none include 'sys_param.inc' include 'gfuni0.inc' include 'globals.inc' call gft_close_all() return end !----------------------------------------------------------------------- ! ! This is the f77 init main ! !----------------------------------------------------------------------- program cyl2d_init implicit none include 'sys_param.inc' include 'globals.inc' include 'other_glbs.inc' real*8 q(memsiz) integer iargc, indlnb, qloc character*256 param_file logical exist if( iargc() .lt. 1 ) go to 900 call getarg(1,param_file) inquire(file=param_file(1:indlnb(param_file)),exist=exist) if( .not. exist ) go to 910 call init_params_attribs() call read_parameters(param_file(1:indlnb(param_file))) call init_coord_difs() call maninit_gfuncs(nrho,nz, & q(qloc(nrho*nz)),q(qloc(nrho*nz)),q(qloc(nrho*nz)), & q(qloc(nrho)),q(qloc(nz)) & ) stop 900 continue write(0,*) 'Usage: cyl2d_init ' stop 910 continue write(0,*) 'cyl2d_init: parameter file '''// & param_file(1:indlnb(param_file))// & ''' does not exist' stop end subroutine maninit_gfuncs(nrhol,nzl, & Phi_n,Phi_nm1,r_n, rhol,zl & ) implicit none integer indlnb integer nrhol, nzl real*8 Phi_n(nrhol,nzl), Phi_nm1(nrhol,nzl), & r_n(nrhol,nzl), * rhol(nrhol), zl(nzl) integer i, j include 'sys_param.inc' include 'globals.inc' include 'other_glbs.inc' call dvmesh(rhol,nrhol,rhomin,rhomax) call dvmesh(zl ,nzl ,zmin, zmax) if( idstr .eq. 'ts' ) then write(0,*) '*** Time symmetric initial data' do j = 1 , nzl do i = 1 , nrhol r_n(i,j) = sqrt(rhol(i)*rhol(i)+zl(j)*zl(j)) Phi_n(i,j) = amp*exp(-((r_n(i,j)-r0)/delta)**2) end do end do do j = 1 , nzl do i = 1 , nrhol Phi_nm1(i,j) = Phi_n(i,j) end do end do else if( idstr .eq. 'in' ) then write(0,*) '*** Ingoing initial data' do j = 1 , nzl do i = 1 , nrhol r_n(i,j) = sqrt(rhol(i)*rhol(i)+zl(j)*zl(j)) Phi_n(i,j) = amp*exp(-((r_n(i,j)-r0)/delta)**2) end do end do do j = 1 , nzl Phi_nm1(i,j) = Phi_n(i,j) do i = 2 , nrhol Phi_nm1(i,j) = Phi_n(i,j) * & (1.0 - dt * (-2 * (r_n(i,j) - r0) / delta**2 + & 1 /r_n(i,j) )) end do end do else if( idstr .eq. 'out' ) then write(0,*) '*** Outgoing initial data' do j = 1 , nzl do i = 1 , nrhol r_n(i,j) = sqrt(rhol(i)*rhol(i)+zl(j)*zl(j)) Phi_n(i,j) = amp*exp(-((r_n(i,j)-r0)/delta)**2) end do end do do j = 1 , nzl Phi_nm1(i,j) = Phi_n(i,j) do i = 2 , nrhol Phi_nm1(i,j) = Phi_n(i,j) * & (1.0 + dt * (-2 * (r_n(i,j) - r0) / delta**2 + & 1 /r_n(i,j) )) end do end do else if( idstr .eq. 'gts' ) then write(0,*) '*** Generalized time-symmetric data' do j = 1 , nzl do i = 1 , nrhol r_n(i,j) = sqrt((rhol(i)/rho0)**2 + & (zl(j)/z0)**2) Phi_n(i,j) = amp*exp(-((r_n(i,j)-r0)/delta)**2) Phi_nm1(i,j) = Phi_n(i,j) end do end do end if !---------------------------------------------------------------------- ! Code fragment for creating initial data file !---------------------------------------------------------------------- call rnpl_id_set_file(in_file(1:indlnb(in_file))) call rnpl_id_set_rank(2) call rnpl_id_set_shape(nrhol,nzl,nzl) call rnpl_id_write('Phi[1]',Phi_n) call rnpl_id_write('Phi[0]',Phi_nm1) call rnpl_id_write('r[0]',r_n) call rnpl_id_end() return end integer function qloc(size) implicit none logical first integer size, loc data first / .true. / save if( first ) then loc = 1 first = .false. end if qloc = loc loc = loc + size return end subroutine dvmesh(v,n,vmin,vmax) implicit none integer n real*8 v(n) real*8 vmin, vmax real*8 dv integer j v(1) = vmin if( n .ge. 2 ) then dv = (vmax - vmin) / (n - 1) do j = 2 , n v(j) = v(j-1) + dv end do end if return end