00001 subroutine outputcube(fd,rd,rpress,rrho,rdb,time,mark,ol)
00002
00003 implicit none
00004 include 'param.f90'
00005 integer ol,i,mark
00006 complex fdin(nx/2+1,ny,nz_max,nvar)
00007 complex fd(nx/2+1,ny,nz_max,nvar)
00008 real rd(nx,ny,nz_max,nvar), rpress(nx,ny,nz_max), rrho(nx,ny,nz_max)
00009 real rdb(nx,ny,nz_max,3), rdfull(nx,ny,nz_ghostglobal),temp(nx,ny,nz_max)
00010 real time
00011 character*80 fn,stringol,zeros
00012
00013 write(stringol,*) mark
00014 zeros=""
00015 if (mark .lt. 10000) zeros="_0"
00016 if (mark .lt. 1000 ) zeros="_00"
00017 if (mark .lt. 100 ) zeros="_000"
00018 if (mark .lt. 10 ) zeros="_0000"
00019 fn="rd.output."//trim(adjustl(zeros))//trim(adjustl(stringol)) &
00020 //".fits"
00021 do i=1,6
00022 temp=rd(:,:,:,i)
00023 call collectall(rdfull,temp)
00024 if(my_rank .eq. 0) then
00025 call fw_cube(fn,rdfull,nx,ny,nz_ghostglobal,dx,dy,dz,time,ol,i)
00026 endif
00027 enddo
00028 end subroutine outputcube
00029
00030
00031
00032
00033
00034
00035
00036 subroutine outputslice(fd,rd,rpress,rrho,rdb,time,ol_in)
00037
00038 implicit none
00039 include 'param.f90'
00040 integer ol,i,ol_in
00041 complex fdin(nx/2+1,ny,nz_max,nvar)
00042 complex fd(nx/2+1,ny,nz_max,nvar)
00043 real rd(nx,ny,nz_max,nvar), rpress(nx,ny,nz_max), rrho(nx,ny,nz_max)
00044 real rdb(nx,ny,nz_max,3), rdfull(nx,ny,nz_ghostglobal),temp(nx,ny,nz_max)
00045 real time
00046 real hslice(nx,ny),vslice(nx,nz_ghostglobal)
00047 character*80 fn
00048 character*60 fns(11)
00049 fns(1)="xdisp.fits"
00050 fns(2)="ydisp.fits"
00051 fns(3)="zdisp.fits"
00052 fns(4)="xdisp_dot.fits"
00053 fns(5)="ydisp_dot.fits"
00054 fns(6)="zdisp_dot.fits"
00055 fns(7)="bx.fits"
00056 fns(8)="by.fits"
00057 fns(9)="bz.fits"
00058 fns(10)="rho.fits"
00059 fns(11)="p.fits"
00060 ol=ol_in+1
00061 do i=1,6
00062 temp=rd(:,:,:,i)
00063 call collectall(rdfull,temp)
00064 if(my_rank .eq. 0) then
00065 fn="hslice_"//fns(i)
00066 hslice=rdfull(:,:,nz_max-100)
00067 call fw_cube(fn,hslice,nx,ny,1,dx,dy,dz,time,ol,ol)
00068 fn="vslice_"//fns(i)
00069 vslice=rdfull(:,yc,:)
00070 call fw_cube(fn,vslice,nx,1,nz_ghostglobal,dx,dy,dz,time,ol,ol)
00071 endif
00072 enddo
00073
00074 call collectall(rdfull,rrho)
00075 if(my_rank .eq. 0) then
00076 fn="hslice_"//fns(10)
00077 hslice=rdfull(:,:,nz_max-100)
00078 call fw_cube(fn,hslice,nx,ny,1,dx,dy,dz,time,ol,ol)
00079 fn="vslice_"//fns(10)
00080 vslice=rdfull(:,yc,:)
00081 call fw_cube(fn,vslice,nx,1,nz_ghostglobal,dx,dy,dz,time,ol,ol)
00082 endif
00083
00084 call collectall(rdfull,rpress)
00085 if(my_rank .eq. 0) then
00086 fn="hslice_"//fns(11)
00087 hslice=rdfull(:,:,nz_max-100)
00088 call fw_cube(fn,hslice,nx,ny,1,dx,dy,dz,time,ol,ol)
00089 fn="vslice_"//fns(11)
00090 vslice=rdfull(:,yc,:)
00091 call fw_cube(fn,vslice,nx,1,nz_ghostglobal,dx,dy,dz,time,ol,ol)
00092 endif
00093
00094 do i=1,3
00095 temp=rdb(:,:,:,i)
00096 call collectall(rdfull,temp)
00097 if(my_rank .eq. 0) then
00098 fn="hslice_"//fns(i+6)
00099 hslice=rdfull(:,:,nz_max-100)
00100 call fw_cube(fn,hslice,nx,ny,1,dx,dy,dz,time,ol,ol)
00101 fn="vslice_"//fns(i+6)
00102 vslice=rdfull(:,yc,:)
00103 call fw_cube(fn,vslice,nx,1,nz_ghostglobal,dx,dy,dz,time,ol,ol)
00104 endif
00105 enddo
00106
00107 end subroutine outputslice
00108
00109
00110
00111
00112
00113
00114
00115 subroutine input(ol,time)
00116
00117 implicit none
00118 include 'param.f90'
00119 include 'terms3d.f90'
00120 character*80 fn,comment,err
00121 character*8 kywrd
00122 real rdinfull(nx,ny,nz_ghostglobal)
00123 real rdin(nx,ny,nz_max)
00124 integer readwrite,group,fpixel,nelements,status,unit,junk
00125 logical er
00126 integer i,k,blocksize,ol,time_int
00127 real time
00128 status=0
00129 comment="Step number"
00130
00131
00132 fn="IC.fits"
00133 blocksize=1
00134 group=1
00135 fpixel=1
00136 readwrite=0
00137 nelements=nx*ny*nz_ghostglobal
00138 if(my_rank .eq. 0) then
00139 call ftgiou(unit,status)
00140 call ftopen(unit,fn,readwrite,blocksize,status)
00141 endif
00142 do i=1,nvar
00143 if(my_rank .eq. 0) then
00144 call FTMAHD(unit,i,junk,status)
00145 call ftgpve(unit,group,fpixel,nelements,0.,rdinfull,er,status)
00146 kywrd="STEP"
00147 call ftgkyj(unit,kywrd,ol,comment,status)
00148 kywrd="TIME"
00149 call ftgkyj(unit,kywrd,time_int,comment,status)
00150 time=time_int/1.e4
00151 endif
00152 call distributeall(rdinfull,rdin)
00153 rd(:,:,:,i)=rdin
00154 enddo
00155 if(my_rank .eq. 0) then
00156 call ftclos(unit, status)
00157 call ftfiou(unit, status)
00158 endif
00159 if(status .ne. 0) then
00160 call ftgerr(status,err)
00161 write(*,*) err
00162 stop
00163 endif
00164
00165
00166
00167
00168 call sfftw_execute(plan_r2fd_1)
00169 call sfftw_execute(plan_r2fd_2)
00170 call sfftw_execute(plan_r2fd_3)
00171 call sfftw_execute(plan_r2fd_4)
00172 call sfftw_execute(plan_r2fd_5)
00173 call sfftw_execute(plan_r2fd_6)
00174 end subroutine input
00175
00176
00177 subroutine distributeall(rdinfull,rdin)
00178 implicit none
00179 include 'param.f90'
00180 real rdinfull(nx,ny,nz_ghostglobal)
00181 real rdin(nx,ny,nz_max)
00182 real temp(nx,ny,nz_max)
00183 integer*8 i
00184 if(my_rank .eq. 0) then
00185
00186 rdin(:,:,1:nz_max)=rdinfull(:,:,1:nz_max)
00187 do i=1,nprocs-1
00188 temp(:,:,1:ends(i+1)-start(i+1)+5)=rdinfull(:,:,start(i+1)-2:ends(i+1)+2)
00189 call mpi_send(temp,nx*ny*nz_max,MPI_REAL,&
00190 & i,9,MPI_COMM_WORLD,ierr)
00191 enddo
00192 endif
00193 if(my_rank .ne. 0) then
00194 call mpi_recv(rdin,nx*ny*nz_max,MPI_REAL,&
00195 & 0,9,MPI_COMM_WORLD,istatus,ierr)
00196 endif
00197 end subroutine distributeall
00198
00199 subroutine collectall(rdinfull,rdin)
00200 implicit none
00201 include 'param.f90'
00202
00203 real rdinfull(nx,ny,nz_ghostglobal)
00204 real rdin(nx,ny,nz_max)
00205 real temp(nx,ny,nz_max)
00206 integer*8 i
00207 if(my_rank .ne. 0) then
00208 call mpi_send(rdin,nx*ny*nz_max,MPI_REAL,&
00209 & 0,2,MPI_COMM_WORLD,ierr)
00210 endif
00211 if(my_rank .eq. 0) then
00212
00213 rdinfull(:,:,1:nz)=rdin(:,:,1:nz)
00214 do i=1,nprocs-1
00215 call mpi_recv(temp,nx*ny*nz_max,MPI_REAL,&
00216 & i,2,MPI_COMM_WORLD,istatus,ierr)
00217 rdinfull(:,:,start(i+1):ends(i+1)+2)=temp(:,:,3:ends(i+1)-start(i+1)+5)
00218 enddo
00219 endif
00220 end subroutine collectall
00221
00222
00223
00224
00225
00226
00227
00228 subroutine fw_cube(fn,data,nx,ny,nz,dx,dy,dz,time,step,k)
00229 implicit none
00230 character*80 fn
00231 character*30 errmsg
00232 real data(nx,ny,nz)
00233 real dx,dy,dz,time
00234 integer dx1,dy1,dz1,step,junk
00235 integer status,unit,blocksize,bitpix,naxis,naxes(3)
00236 integer i,j,group,fpixel,nelements,nx,ny,nz
00237 integer n,k,time_int,temp
00238 logical simple,extend
00239
00240 group=1
00241 fpixel=1
00242 simple=.true.
00243 status=0
00244 bitpix=-32
00245 naxis=3
00246 naxes(1)=nx
00247 naxes(2)=ny
00248 naxes(3)=nz
00249 extend=.true.
00250 status=0
00251 blocksize=1
00252 nelements=naxes(1)*naxes(2)*naxes(3)
00253 call ftgiou(unit,status)
00254 if(k .eq. 1) then
00255 call ftinit(unit,fn,blocksize,status)
00256 endif
00257 if(k .ne. 1) then
00258 call FTOPEN(unit,fn,1,blocksize,status)
00259 call FTMAHD(unit,k-1,junk,status)
00260 write(*,*) "attempting", k,fn
00261 call FTGHDN(unit,temp)
00262 write(*,*) "now at ", temp
00263 call ftthdu(unit,temp,status)
00264 write(*,*) "Deleting", k,":",temp
00265 do i=temp,k,-1
00266 write(*,*) i
00267 call FTMAHD(unit,i,junk,status)
00268 write(*,*) "attempting delete", i
00269 call FTDHDU(unit,junk,status)
00270 enddo
00271 call FTCRHD(unit,status)
00272 call FTGHDN(unit,temp)
00273 write(*,*) "now at ",temp
00274 endif
00275
00276 call ftphpr(unit,simple,bitpix,naxis,naxes,0,1,extend,status)
00277 time_int=1.e4*time
00278 call ftpkyj(unit,'time',time_int,'restart time (10^-4 s)',status)
00279 call ftpkyj(unit,'step',step,'Step number',status)
00280 dx1=dx/1e2
00281 dy1=dy/1e2
00282 dz1=dz/1e2
00283 call ftpkyj(unit,'dx',dx1,'m',status)
00284 call ftpkyj(unit,'dy',dy1,'m',status)
00285 call ftpkyj(unit,'dz',dz1,'m',status)
00286
00287 call ftppre(unit,group,fpixel,nelements,data,status)
00288
00289
00290 call ftclos(unit, status)
00291 call ftfiou(unit, status)
00292 if (status .ne. 0) then
00293 call ftgerr(status,errmsg)
00294 write(*,*) "err: ", errmsg
00295 stop
00296 endif
00297 end subroutine fw_cube
00298
00299
00300
