SLiMv2.1: SLiMv2.1/extract_rho.f90 Source File

00001   program flatten_all
00002   character*80 fn
00003   fn="coeff_wavesim.fits"
00004   call flatten_cube(fn)
00005   end program
00006 
00007 !-------------------------------------------------------------------------------------------------
00008 
00009   subroutine flatten_cube(fn)
00010 ! takes a slice file and reduces it to a 3-D cube
00011   implicit none
00012 !  include 'param.f90'
00013   integer nx,ny,nz_ghostglobal,j
00014   parameter(nx=200,ny=50,nz_ghostglobal=556) 
00015   real dx,dy,dz
00016   parameter(dx=145.041E8/(nx), dy=2.*dx, dz=14.e8/(nz_ghostglobal-1))
00017 
00018 
00019   character*80 fn,comment,err
00020   real rdinfull(nx,ny,nz_ghostglobal)
00021   integer nsamp
00022   parameter(nsamp=227)
00023   real rd(nx,ny,nz_ghostglobal),bd(nx,ny,nz_ghostglobal),flux(nx,nz_ghostglobal),r,r1,fout(nx,ny,nz_ghostglobal)
00024   integer readwrite,group,fpixel,nelements,status,unit,junk,naxes(4),naxis,xc,rint
00025   integer dx1,dy1,dz1,bitpix,yc
00026   
00027   logical er,simple,extend
00028   integer i,k,blocksize,ol
00029   status=0
00030 
00031 
00032   yc=24
00033   blocksize=1
00034   group=1
00035   fpixel=1
00036   readwrite=0
00037   nelements=nx*ny*nz_ghostglobal
00038   extend=.true.
00039   call ftgiou(unit,status)
00040   call ftopen(unit,fn,readwrite,blocksize,status)
00041 !load density
00042      i=1
00043      write(*,*) i,status
00044         call FTMAHD(unit,i,junk,status)
00045         call ftgpve(unit,group,fpixel,nelements,0.,rdinfull,er,status)
00046      write(*,*) i,status   
00047      rd(:,:,:)=rdinfull
00048   
00049 !load bz
00050      i=17
00051      write(*,*) i,status
00052         call FTMAHD(unit,i,junk,status)
00053         call ftgpve(unit,group,fpixel,nelements,0.,rdinfull,er,status)
00054      write(*,*) i,status
00055      bd(:,:,:)=rdinfull
00056 
00057   call ftclos(unit, status)
00058   call ftfiou(unit, status)
00059 
00060   If(status .ne. 0) then
00061         call ftgerr(status,err)
00062         write(*,*) "read error: ",err
00063         stop
00064   endif
00065 
00066 !calculate flux
00067   xc=40*200./145.2
00068   flux(1,:)=0
00069   do i=2,nx-xc
00070      do j=1,nz_ghostglobal
00071         flux(i,j)=flux(i-1,j)+(i)*bd(xc+i-1,24,j)
00072      enddo
00073   enddo
00074   do i=nx-xc+1,nx
00075         flux(i,:)=flux(i-1,:)
00076   enddo
00077 ! load fout
00078    do i=1,nx
00079      do j=1,ny
00080      r=sqrt(float(i-xc)**2+(2.*float(j-yc))**2)
00081      rint=int(r)
00082      r1=r-rint
00083      rint=max(min(rint,nx-1),1)
00084      fout(i,j,:)=flux(rint,:)*(1-r1)+flux(rint+1,:)*r1
00085      enddo
00086    enddo
00087 
00088 
00089  write(*,*) "out" 
00090 !now OUTPUT
00091   fn="flat_"//fn
00092 !write_rho
00093   fn="rho_sim.fits"    
00094   call ftgiou(unit,status)
00095   call ftinit(unit,fn,blocksize,status)
00096   simple=.true.
00097   bitpix=-32
00098   naxis=3
00099   naxes(1)=nx   
00100   naxes(2)=ny
00101   naxes(3)=nz_ghostglobal
00102   dx1=dx/1e2        
00103   dy1=dy/1e2
00104   dz1=dz/1e2         
00105  write(*,*) status,nx,ny,nz_ghostglobal 
00106  call ftphpr(unit,simple,bitpix,naxis,naxes,0,1,extend,status)
00107   call ftpkyj(unit,'dx',dx1,'m',status)
00108   call ftpkyj(unit,'dy',dy1,'m',status)
00109   call ftpkyj(unit,'dz',dz1,'m',status)
00110 !     write the array to the FITS file
00111   group=1
00112   fpixel=1
00113 !  nelements=naxes(1)*naxes(2)*naxes(3) 
00114   call ftppre(unit,group,fpixel,nelements,rd,status) 
00115   call ftclos(unit, status)
00116   call ftfiou(unit, status)
00117 
00118  write(*,*) "1 done"
00119   fn="bz_sim.fits"
00120   call ftgiou(unit,status)
00121   write(*,*) 1,status
00122   call ftinit(unit,fn,blocksize,status)
00123   simple=.true.
00124   bitpix=-32
00125   naxis=3
00126   naxes(1)=nx
00127   naxes(2)=ny
00128   naxes(3)=nz_ghostglobal
00129  
00130   dx1=dx/1e2
00131   dy1=dy/1e2
00132   dz1=dz/1e2
00133   call ftphpr(unit,simple,bitpix,naxis,naxes,0,1,extend,status)
00134   call ftpkyj(unit,'dx',dx1,'m',status)
00135   call ftpkyj(unit,'dy',dy1,'m',status)
00136   call ftpkyj(unit,'dz',dz1,'m',status)
00137 !     write the array to the FITS file
00138   group=1
00139   fpixel=1
00140   call ftppre(unit,group,fpixel,nelements,bd,status)
00141   call ftclos(unit, status)
00142   call ftfiou(unit, status)
00143 
00144 !flux out
00145   fn="flux_sim.fits"
00146   call ftgiou(unit,status)
00147   write(*,*) 1,status
00148   call ftinit(unit,fn,blocksize,status)
00149   simple=.true.
00150   bitpix=-32
00151   naxis=3
00152   naxes(1)=nx
00153   naxes(2)=ny
00154   naxes(3)=nz_ghostglobal
00155 
00156   dx1=dx/1e2
00157   dy1=dy/1e2
00158   dz1=dz/1e2
00159   call ftphpr(unit,simple,bitpix,naxis,naxes,0,1,extend,status)
00160   call ftpkyj(unit,'dx',dx1,'m',status)
00161   call ftpkyj(unit,'dy',dy1,'m',status)
00162   call ftpkyj(unit,'dz',dz1,'m',status)
00163 !     write the array to the FITS file
00164   group=1
00165   fpixel=1
00166   call ftppre(unit,group,fpixel,nelements,fout,status)
00167   call ftclos(unit, status)
00168   call ftfiou(unit, status)
00169 
00170 
00171 
00172   if (status .ne. 0) then
00173      call ftgerr(status,err)
00174      write(*,*) "write err: ", err
00175      stop        
00176   endif     
00177 end subroutine