SLiMv2.1: SLiMv2.1/extract_flux.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_global),r,fout(nx,ny,nz_ghostglobal)
00024   integer readwrite,group,fpixel,nelements,status,unit,junk,naxes(4),naxis,xc
00025   integer dx1,dy1,dz1,bitpix
00026   
00027   logical er,simple,extend
00028   integer i,k,blocksize,ol
00029   status=0
00030 
00031   blocksize=1
00032   group=1
00033   fpixel=1
00034   readwrite=0
00035   nelements=nx*ny*nz_ghostglobal
00036   extend=.true.
00037   call ftgiou(unit,status)
00038   call ftopen(unit,fn,readwrite,blocksize,status)
00039 !load density
00040      i=1
00041      write(*,*) i,status
00042         call FTMAHD(unit,i,junk,status)
00043         call ftgpve(unit,group,fpixel,nelements,0.,rdinfull,er,status)
00044      write(*,*) i,status   
00045      rd(:,:,:)=rdinfull
00046   
00047 !load bz
00048      i=17
00049      write(*,*) i,status
00050         call FTMAHD(unit,i,junk,status)
00051         call ftgpve(unit,group,fpixel,nelements,0.,rdinfull,er,status)
00052      write(*,*) i,status
00053      bd(:,:,:)=rdinfull
00054 
00055   call ftclos(unit, status)
00056   call ftfiou(unit, status)
00057 
00058   If(status .ne. 0) then
00059         call ftgerr(status,err)
00060         write(*,*) "read error: ",err
00061         stop
00062   endif
00063 
00064 !calculate flux
00065   xc=40*200./145.2
00066   flux(1,:)=0
00067   do i=2,nx-xc
00068      do j=1,nz_ghostglobal
00069         flux(i,j)=flux(i-1,j)+(i-1)*bd(i-xc,24,j)
00070      enddo
00071   enddo
00072   do i=nx-xc+1:nx
00073         flux(i,:)=flux(i-1,:)
00074   enddo
00075 ! load fout
00076    do i=1,nx
00077      do j=1,ny
00078      r=sqrt(i**2+j**2)
00079      fout(i,j,:)=flux(max(min(int(r),nx),1),:)
00080      enddo
00081    enddo
00082 
00083 
00084  write(*,*) "out" 
00085 !now OUTPUT
00086   fn="flat_"//fn
00087 !write_rho
00088   fn="rho_sim.fits"    
00089   call ftgiou(unit,status)
00090   call ftinit(unit,fn,blocksize,status)
00091   simple=.true.
00092   bitpix=-32
00093   naxis=3
00094   naxes(1)=nx   
00095   naxes(2)=ny
00096   naxes(3)=nz_ghostglobal
00097   dx1=dx/1e2        
00098   dy1=dy/1e2
00099   dz1=dz/1e2         
00100  write(*,*) status,nx,ny,nz_ghostglobal 
00101  call ftphpr(unit,simple,bitpix,naxis,naxes,0,1,extend,status)
00102   call ftpkyj(unit,'dx',dx1,'m',status)
00103   call ftpkyj(unit,'dy',dy1,'m',status)
00104   call ftpkyj(unit,'dz',dz1,'m',status)
00105 !     write the array to the FITS file
00106   group=1
00107   fpixel=1
00108 !  nelements=naxes(1)*naxes(2)*naxes(3) 
00109   call ftppre(unit,group,fpixel,nelements,rd,status) 
00110   call ftclos(unit, status)
00111   call ftfiou(unit, status)
00112 
00113  write(*,*) "1 done"
00114   fn="bz_sim.fits"
00115   call ftgiou(unit,status)
00116   write(*,*) 1,status
00117   call ftinit(unit,fn,blocksize,status)
00118   simple=.true.
00119   bitpix=-32
00120   naxis=3
00121   naxes(1)=nx
00122   naxes(2)=ny
00123   naxes(3)=nz_ghostglobal
00124  
00125   dx1=dx/1e2
00126   dy1=dy/1e2
00127   dz1=dz/1e2
00128   call ftphpr(unit,simple,bitpix,naxis,naxes,0,1,extend,status)
00129   call ftpkyj(unit,'dx',dx1,'m',status)
00130   call ftpkyj(unit,'dy',dy1,'m',status)
00131   call ftpkyj(unit,'dz',dz1,'m',status)
00132 !     write the array to the FITS file
00133   group=1
00134   fpixel=1
00135   call ftppre(unit,group,fpixel,nelements,bd,status)
00136   call ftclos(unit, status)
00137   call ftfiou(unit, status)
00138 
00139 !flux out
00140   fn="flux_sim.fits"
00141   call ftgiou(unit,status)
00142   write(*,*) 1,status
00143   call ftinit(unit,fn,blocksize,status)
00144   simple=.true.
00145   bitpix=-32
00146   naxis=3
00147   naxes(1)=nx
00148   naxes(2)=ny
00149   naxes(3)=nz_ghostglobal
00150 
00151   dx1=dx/1e2
00152   dy1=dy/1e2
00153   dz1=dz/1e2
00154   call ftphpr(unit,simple,bitpix,naxis,naxes,0,1,extend,status)
00155   call ftpkyj(unit,'dx',dx1,'m',status)
00156   call ftpkyj(unit,'dy',dy1,'m',status)
00157   call ftpkyj(unit,'dz',dz1,'m',status)
00158 !     write the array to the FITS file
00159   group=1
00160   fpixel=1
00161   call ftppre(unit,group,fpixel,nelements,fout,status)
00162   call ftclos(unit, status)
00163   call ftfiou(unit, status)
00164 
00165 
00166 
00167   if (status .ne. 0) then
00168      call ftgerr(status,err)
00169      write(*,*) "write err: ", err
00170      stop        
00171   endif     
00172 end subroutine