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

00001 program makeIC
00002   ! write out result
00003   implicit none
00004   include 'param.f90'
00005   include 'terms3d.f90'
00006   integer i,j,k,j1,n,l
00007   real z,rsun,sigma,kz,omega,a(nz_ghostglobal),amplitude,r,omega_star
00008   real mask(nx/2+1)
00009   real data(nx,ny,nz_ghostglobal,nvar)
00010   real tdata(nx,ny,nz_ghostglobal)
00011 
00012   character*80 fn
00013   parameter (rsun=6.98e10)
00014   
00015   fn="ICstuff/dispx.fits" 
00016   call fr(fn,tdata,nx,ny,nz_ghostglobal,dx,dy,dz) 
00017   write(*,*) "x-d in"
00018   data(:,:,:,1)=tdata
00019   fn="ICstuff/dispy.fits" 
00020   call fr(fn,tdata,nx,ny,nz_ghostglobal,dx,dy,dz) 
00021   data(:,:,:,2)=tdata
00022   fn="ICstuff/dispz.fits" 
00023   call fr(fn,tdata,nx,ny,nz_ghostglobal,dx,dy,dz) 
00024   data(:,:,:,3)=tdata
00025   fn="ICstuff/vx.fits" 
00026   call fr(fn,tdata,nx,ny,nz_ghostglobal,dx,dy,dz) 
00027   data(:,:,:,4)=tdata
00028   fn="ICstuff/vy.fits" 
00029   call fr(fn,tdata,nx,ny,nz_ghostglobal,dx,dy,dz) 
00030   data(:,:,:,5)=tdata
00031   fn="ICstuff/vz.fits" 
00032   call fr(fn,tdata,nx,ny,nz_ghostglobal,dx,dy,dz) 
00033   data(:,:,:,6)=tdata
00034 
00035 
00036   fn="IC.fits" 
00037   call fw(fn,data,nx,ny,nz_ghostglobal,dx,dy,dz,0.,0) 
00038   write(*,*) "Success" 
00039 end program makeIC
00040 
00041 !============================================================================
00042       subroutine fr(fn,data,nx,ny,nz)
00043       implicit none     
00044       character*80 fn
00045       character*30 errmsg
00046       real data(nx,ny,nz), data1(nx,ny,nz)    
00047       real dx,dy,dz,time
00048       integer dx1,dy1,dz1,step
00049       integer status,unit,blocksize,bitpix,naxis,naxes(3)
00050       integer i,j,group,fpixel,nelements,nx,ny,nz
00051       integer n,k
00052       logical simple,extend
00053       integer readwrite
00054       logical er
00055 
00056       status=0
00057       readwrite=0
00058       blocksize=1
00059       group=1
00060       fpixel=1   
00061       nelements=nx*ny*nz
00062       
00063       call ftgiou(unit,status)
00064       call ftopen(unit,fn,readwrite,blocksize,status)
00065       call ftgpve(unit,group,fpixel,nelements,0.,data,er,status)
00066       call ftclos(unit, status)
00067       call ftfiou(unit, status)
00068       write(*,*) fn,maxval(data),status 
00069       if (status .ne. 0) then
00070             call ftgerr(status,errmsg)
00071             write(*,*) fn,nx,ny,nz,"err: ", errmsg  
00072             stop        
00073       endif     
00074       end subroutine fr
00075 !=============================================================================
00076       subroutine fw(fn,data,nx,ny,nz,dx,dy,dz,time,step)
00077       implicit none     
00078       character*80 fn
00079       character*30 errmsg
00080       real data(nx,ny,nz,6), data1(nx,ny,nz)    
00081       real dx,dy,dz,time
00082       integer dx1,dy1,dz1,step
00083       integer status,unit,blocksize,bitpix,naxis,naxes(3)
00084       integer i,j,group,fpixel,nelements,nx,ny,nz
00085       integer n,k
00086       logical simple,extend
00087 
00088       status=0
00089       blocksize=1
00090   
00091       call ftgiou(unit,status)
00092       
00093       blocksize=1
00094       call ftinit(unit,fn,blocksize,status)
00095       simple=.true.
00096       bitpix=-32
00097       naxis=3
00098       naxes(1)=nx
00099       naxes(2)=ny
00100       naxes(3)=nz       
00101       extend=.true.
00102 
00103 !     write the required header keywords
00104       call ftphpr(unit,simple,bitpix,naxis,naxes,0,1,extend,status)
00105       call ftpkyj(unit,'time',0,'restart time',status)
00106       call ftpkyj(unit,'step',step,'Step number',status)
00107       call ftpkyj(unit,'W-MODE',1,'F mode plane wave',status)
00108       dx1=dx/1e2        
00109       dy1=dy/1e2
00110       dz1=dz/1e2         
00111       call ftpkyj(unit,'dx',dx1,'m',status)
00112       call ftpkyj(unit,'dy',dy1,'m',status)
00113       call ftpkyj(unit,'dz',dz1,'m',status)
00114 
00115 !     write the array to the FITS file
00116       group=1
00117       fpixel=1
00118       nelements=naxes(1)*naxes(2)*naxes(3)
00119       data1=data(:,:,:,1) 
00120         write(*,*) k,":",maxval(data1)
00121       call ftppre(unit,group,fpixel,nelements,data1,status)     
00122       do k=2,6  
00123         data1=data(:,:,:,k)
00124         call ftcrhd(unit,status)
00125         call ftphpr(unit,simple,bitpix,naxis,naxes,0,1,extend,status)
00126         call ftpkyj(unit,'time',0,'restart time',status)
00127         call ftpkyj(unit,'step',step,'Step number',status)
00128         call ftpkyj(unit,'W-MODE',1,'F mode plane wave',status)
00129         dx1=dx/1e2      
00130         dy1=dy/1e2
00131         dz1=dz/1e2       
00132         call ftpkyj(unit,'dx',dx1,'m',status)
00133         call ftpkyj(unit,'dy',dy1,'m',status)
00134         call ftpkyj(unit,'dz',dz1,'m',status)
00135         write(*,*) k,":",maxval(data1)
00136 
00137         call ftppre(unit,group,fpixel,nelements,data1,status)
00138       enddo             
00139 
00140 !     close the file and free the unit number
00141       call ftclos(unit, status)
00142       call ftfiou(unit, status)
00143       if (status .ne. 0) then
00144             call ftgerr(status,errmsg)
00145             write(*,*) "err: ", errmsg  
00146             stop        
00147       endif     
00148       end subroutine fw
00149 
00150 
00151