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

00001 program makecoeff
00002  !#define threed 1
00003   ! program to make coefficients for helioseismology runs
00004   implicit none
00005   include 'comments.f90'
00006   include 'param.f90'
00007   real B0(nx,ny,nz_ghostglobal+2,3),U0(nx,ny,nz_ghostglobal+2,3),bztemp(nx,ny,nz_ghostglobal+2,3)
00008   real pressure(nx,ny,nz_ghostglobal+2),gamma(nx,ny,nz_ghostglobal+2),rho(nx,ny,nz_ghostglobal+2)
00009   real vtemp(nx,ny,nz_ghostglobal+2,3)
00010   real coeff(nx,ny,nz_ghostglobal+2,ncoeff)
00011   !  real*4 coeffout(nx,ny,nz,ncoeff)
00012   real temp(nx,ny,nz_ghostglobal+2),fac,data(nx,ny,nz_ghostglobal)
00013   character*80 fn,fnb
00014   integer k,k1,i,j,n_rec
00015   interface
00016      function grad(datain) 
00017        include 'param.f90'
00018        real datain(nx,ny,nz_ghostglobal+2),grad(nx,ny,nz_ghostglobal+2,3)
00019      end function grad
00020   end interface
00021   interface
00022      function ave(datain)
00023        include 'param.f90'
00024        real datain(nx,ny,nz_ghostglobal+2),ave(nx,ny,nz_ghostglobal+2,3)
00025      end function ave
00026   end interface
00027 
00028 !first read in the files
00029 
00030         
00031   fn="spot/bx1_wavesim.fits"
00032   call fr(fn,temp,nx,ny,nz_ghostglobal+2)
00033   b0(:,:,:,1)=temp
00034   fn="spot/by1_wavesim.fits"
00035   call fr(fn,temp,nx,ny,nz_ghostglobal+2)
00036   b0(:,:,:,2)=temp
00037   fn="spot/bz1_wavesim.fits"
00038   call fr(fn,temp,nx,ny,nz_ghostglobal+2)
00039   b0(:,:,:,3)=temp
00040 
00041   fn="spot/pressure_wavesim.fits"
00042   call fr(fn,pressure,nx,ny,nz_ghostglobal+2)
00043 
00044   fn="spot/rho_wavesim.fits"
00045   call fr(fn,rho,nx,ny,nz_ghostglobal+2)
00046 
00047   fn="spot/gamma1_wavesim.fits"
00048   call fr(fn,gamma,nx,ny,nz_ghostglobal+2)
00049 
00050   fn="spot/ux1_wavesim.fits"
00051   call fr(fn,temp,nx,ny,nz_ghostglobal+2)
00052   u0(:,:,:,1)=temp
00053   fn="spot/uy1_wavesim.fits"
00054   call fr(fn,temp,nx,ny,nz_ghostglobal+2)
00055   u0(:,:,:,2)=temp
00056   fn="spot/uz1_wavesim.fits"
00057   call fr(fn,temp,nx,ny,nz_ghostglobal+2)
00058   u0(:,:,:,3)=temp
00059 
00060 
00061 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
00062  ! B0=0
00063 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00064   ! next calculate the coefficient array
00065   coeff(:,:,:,1)=rho 
00066  ! coeff(:,:,:,2:4)=grad(rho)
00067   vtemp=grad(alog(rho))
00068   coeff(:,:,:,2)=rho*vtemp(:,:,:,1)
00069   coeff(:,:,:,3)=rho*vtemp(:,:,:,2)
00070   coeff(:,:,:,4)=rho*vtemp(:,:,:,3)
00071   coeff(:,:,:,5)= Gamma * pressure / rho
00072 
00073 !  coeff(:,:,:,9:11)=Grad(pressure)
00074   vtemp=Grad(alog(pressure)) 
00075   coeff(:,:,:,9)=pressure*vtemp(:,:,:,1)
00076   coeff(:,:,:,10)=pressure*vtemp(:,:,:,2)
00077   coeff(:,:,:,11)=pressure*vtemp(:,:,:,3)
00078 
00079 
00080   !  coeff(:,:,:,9:10)=0
00081   !  coeff(:,:,:,11)=-rho *g 
00082   coeff(:,:,:,6)= Gamma*pressure*coeff(:,:,:,2)/rho -coeff(:,:,:,9)
00083   coeff(:,:,:,7)= Gamma*pressure*coeff(:,:,:,3)/rho -coeff(:,:,:,10)
00084   coeff(:,:,:,8)= Gamma*pressure*coeff(:,:,:,4)/rho -coeff(:,:,:,11)
00085   !  coeff(:,:,:,6)= -coeff(:,:,:,9)
00086   !  coeff(:,:,:,7)= -coeff(:,:,:,10)
00087   !  coeff(:,:,:,8)= -coeff(:,:,:,11)
00088   coeff(:,:,:,12)= 1. / rho
00089   coeff(:,:,:,13)= g / rho
00090   coeff(:,:,:,14)= g
00091   !magnetic field stuff
00092   coeff(:,:,:,15)= B0(:,:,:,1)
00093   coeff(:,:,:,16)= B0(:,:,:,2)
00094   coeff(:,:,:,17)= B0(:,:,:,3)
00095   vtemp=Grad(B0(:,:,:,1))
00096   coeff(:,:,:,19)=  vtemp(:,:,:,3)
00097   coeff(:,:,:,20)= -vtemp(:,:,:,2)
00098   vtemp=Grad(B0(:,:,:,2))
00099   coeff(:,:,:,18)= -vtemp(:,:,:,3)
00100   coeff(:,:,:,20)= coeff(:,:,:,20)+ vtemp(:,:,:,1)
00101   vtemp=Grad(B0(:,:,:,3))
00102   coeff(:,:,:,18)= coeff(:,:,:,18)+vtemp(:,:,:,2)
00103   coeff(:,:,:,19)= coeff(:,:,:,19)-vtemp(:,:,:,1)
00104   !velocity field stuff
00105   coeff(:,:,:,21)= U0(:,:,:,1)
00106   coeff(:,:,:,22)= U0(:,:,:,2)
00107   coeff(:,:,:,23)= U0(:,:,:,3)
00108   vtemp=Grad(U0(:,:,:,1))
00109   coeff(:,:,:,24)= vtemp(:,:,:,1)
00110   coeff(:,:,:,27)= vtemp(:,:,:,2)
00111   coeff(:,:,:,30)= vtemp(:,:,:,3)
00112   vtemp=Grad(U0(:,:,:,2))
00113   coeff(:,:,:,25)= vtemp(:,:,:,1)
00114   coeff(:,:,:,28)= vtemp(:,:,:,2)
00115   coeff(:,:,:,31)= vtemp(:,:,:,3)
00116   vtemp=Grad(U0(:,:,:,3))
00117   coeff(:,:,:,26)= vtemp(:,:,:,1)
00118   coeff(:,:,:,29)= vtemp(:,:,:,2)
00119   coeff(:,:,:,32)= vtemp(:,:,:,3)
00120 
00121   coeff(:,:,:,33)= U0(:,:,:,1)/rho
00122   coeff(:,:,:,34)= U0(:,:,:,2)/rho
00123   coeff(:,:,:,35)= U0(:,:,:,3)/rho
00124   coeff(:,:,:,36)=Gamma*pressure
00125   !---------------------
00126   Bztemp=b0(:,:,:,:)
00127   where (abs(Bztemp) .lt. 1e-9) Bztemp=1e-9
00128   !---------------------
00129   coeff(:,:,:,37)=B0(:,:,:,1)/Bztemp(:,:,:,3)
00130   coeff(:,:,:,38)=B0(:,:,:,2)/Bztemp(:,:,:,3)
00131   vtemp(:,:,:,1)=pressure+&
00132        &     (b0(:,:,:,1)**2+b0(:,:,:,2)**2-b0(:,:,:,3)**2)/(8.*pi)
00133   coeff(:,:,:,39:41)=Grad(vtemp(:,:,:,1))
00134 
00135   !---------------------
00136   Bztemp=ave(b0(:,:,:,3))
00137   where (abs(Bztemp) .lt. 1e-6) Bztemp=1e-6
00138   !---------------------
00139 
00140   vtemp(:,:,:,1)=B0(:,:,:,1)*B0(:,:,:,3)
00141   coeff(:,:,:,42:44)=Grad(vtemp(:,:,:,1))
00142   coeff(:,:,:,42)=coeff(:,:,:,42)/bztemp(:,:,:,1)
00143   coeff(:,:,:,43)=coeff(:,:,:,43)/bztemp(:,:,:,2)
00144   coeff(:,:,:,44)=coeff(:,:,:,44)/bztemp(:,:,:,3)
00145   vtemp(:,:,:,1)=B0(:,:,:,2)*B0(:,:,:,3)
00146   coeff(:,:,:,45:47)=Grad(vtemp(:,:,:,1))
00147   coeff(:,:,:,45)=coeff(:,:,:,45)/bztemp(:,:,:,1)
00148   coeff(:,:,:,46)=coeff(:,:,:,46)/bztemp(:,:,:,2)
00149   coeff(:,:,:,47)=coeff(:,:,:,47)/bztemp(:,:,:,3)
00150   coeff(:,:,:,48)=1./bztemp(:,:,:,3)
00151 
00152   !reduce Lorentz force in upper domain.
00153   vtemp(:,:,:,1)=pressure+&
00154        &     (b0(:,:,:,1)**2+b0(:,:,:,2)**2-b0(:,:,:,3)**2)/(8.*pi)*0
00155 
00156         do k=nz_ghostglobal+1-50,nz_ghostglobal+1
00157            vtemp(:,:,k,1)=pressure(:,:,k)+&
00158        &     (b0(:,:,k,1)**2+b0(:,:,k,2)**2-b0(:,:,k,3)**2)/(8.*pi)* &
00159        &    (coeff(nx,ny,k,1)/coeff(nx,ny,nz_ghostglobal+1-51,1))**4*0
00160 !           coeff(i,j,k,48)=coeff(i,j,k,48)*(fac+(1-fac)&
00161 !                &*(coeff(nx,ny,k,1)/coeff(nx,ny,nz_ghostglobal+1-51,1)))
00162         enddo
00163            coeff(:,:,:,39:41)=Grad(vtemp(:,:,:,1))
00164 !do the output
00165   fnb="coeff_wavesim.fits"
00166   fn=fnb
00167   data=coeff(:,:,2:nz_ghostglobal+1,1)
00168   write(*,*) maxval(data)
00169   call fw(fn,data,nx,ny,nz_ghostglobal,0) 
00170   do k=2,ncoeff
00171         n_rec=k-1
00172         data=coeff(:,:,2:nz_ghostglobal+1,k)
00173         if(n_rec .le. 10) write(fn,'(I1)') n_rec-1
00174         if(n_rec .gt. 10) write(fn,'(I2)') n_rec-1
00175         fn=trim(fnb)//'+'//trim(fn)
00176         call  fw(fn,data,nx,ny,nz_ghostglobal,k-1)  
00177   enddo
00178   write(*,*) "Success"
00179 
00180 end program makecoeff
00181 !!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00182 function grad(datain)
00183   implicit none
00184   include 'param.f90'
00185   real datain(nx,ny,nz_ghostglobal+2)
00186   real grad(  nx,ny,nz_ghostglobal+2,3)
00187   grad(2:nx-1,:,:,1)=(datain(3:nx,:,:)-datain(1:nx-2,:,:))/(2.*dx)
00188 #if threed 
00189   grad(:,2:ny-1,:,2)=(datain(:,3:ny,:)-datain(:,1:ny-2,:))/(2.*dy)
00190 #endif
00191   grad(:,:,2:nz_ghostglobal+1,3)=(datain(:,:,3:nz_ghostglobal+2)&
00192        &                          -datain(:,:,1:nz_ghostglobal))/(2.*dz)
00193 !  grad(:,:,:,3)=0
00194 
00195 ! x,y and z directions are periodic
00196   grad( 1,:,:,1)=(datain(2,:,:)-datain(nx  ,:,:))/(2.*dx)
00197   grad(nx,:,:,1)=(datain(1,:,:)-datain(nx-1,:,:))/(2.*dx)
00198 #if threed 
00199   grad(:, 1,:,2)=(datain(:,2,:)-datain(:,ny,:  ))/(2.*dy)
00200   grad(:,ny,:,2)=(datain(:,1,:)-datain(:,ny-1,:))/(2.*dy)
00201 #endif
00202 
00203   grad(:,:, 1,3)=grad(:,:,2,3)
00204   grad(:,:,nz_ghostglobal+2,3)=grad(:,:,nz_ghostglobal+1,3)
00205 
00206 
00207 end function grad
00208 !!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00209 function ave(datain)
00210   implicit none
00211   include 'param.f90'
00212   real datain(nx,ny,nz_ghostglobal+2)
00213   real ave(  nx,ny,nz_ghostglobal+2,3)
00214   ave(2:nx-1,:,:,1)=(datain(3:nx,:,:)+datain(1:nx-2,:,:))/(2.)
00215 #if threed 
00216   ave(:,2:ny-1,:,2)=(datain(:,3:ny,:)+datain(:,1:ny-2,:))/(2.)
00217 #endif
00218   ave(:,:,2:nz+1,3)=(datain(:,:,3:nz+2)+datain(:,:,1:nz))/(2.)
00219 
00220 ! x,y and z directions are periodic
00221   ave( 1,:,:,1)=(datain(2,:,:)+datain(nx  ,:,:))/(2.)
00222   ave(nx,:,:,1)=(datain(1,:,:)+datain(nx-1,:,:))/(2.)
00223 #if threed 
00224   ave(:, 1,:,2)=(datain(:,2,:)+datain(:,ny,:  ))/(2.)
00225   ave(:,ny,:,2)=(datain(:,1,:)+datain(:,ny-1,:))/(2.)
00226 #endif
00227 
00228   ave(:,:, 1,3)=ave(:,:,2,3)
00229   ave(:,:,nz+2,3)=ave(:,:,nz+1,3)
00230 end function ave
00231 !========================================================================================================================
00232       subroutine fr(fn,data,nx,ny,nz)
00233       implicit none     
00234       integer nx,ny,nz,blocksize        
00235       character*80 fn
00236       real data(nx,ny,nz)
00237       integer readwrite,group,fpixel,nelements,status,unit
00238       logical er
00239       character*30 errmsg       
00240 
00241       readwrite=0       
00242       blocksize=1
00243       status=0  
00244       call ftgiou(unit,status)  
00245       call ftopen(unit,fn,readwrite,blocksize,status)
00246       group=1
00247       fpixel=1
00248       nelements=nx*ny*nz
00249       call ftgpve(unit,group,fpixel,nelements,0.,data,er,status)
00250       call ftclos(unit, status)
00251       call ftfiou(unit, status)
00252 
00253       if(status .ne. 0) then
00254         call ftgerr(status,errmsg)
00255         write(*,*) "reading: ",fn
00256         write(*,*) "readin error: ",status,errmsg
00257         stop
00258       endif                     
00259       end subroutine fr  
00260 !========================================================================================================================
00261       subroutine fw(fn,data,nx,ny,nz,n_rec)
00262       implicit none     
00263       character*80 fn,record,infn
00264       character*30 errmsg
00265       real data(nx,ny,nz),alpha,b_0     
00266       integer status,unit,blocksize,bitpix,naxis,naxes(3)
00267       integer i,j,group,fpixel,nelements,nx,ny,nz
00268       integer n,n_rec,readwrite,nkeys,nspace,iunit
00269       logical simple,extend
00270 
00271       status=0
00272       write(*,*) n_rec,minval(data),maxval(data)        
00273 !     Get an unused Logical Unit Number to use to create the FITS file
00274       call ftgiou(unit,status)
00275       call  ftgiou(iunit,status)
00276       blocksize=1
00277       readwrite=1  
00278       if(n_rec .eq. 0)  then
00279         call ftinit(unit,fn,blocksize,status)
00280       else 
00281         call ftnopn(unit,fn,blocksize,status)
00282         call ftcrhd(unit,status)
00283       endif             
00284       simple=.true.
00285       bitpix=-32
00286       naxis=3
00287       naxes(1)=nx
00288       naxes(2)=ny
00289       naxes(3)=nz       
00290       extend=.true.
00291       call ftphpr(unit,simple,bitpix,naxis,naxes,0,1,extend,status)
00292 
00293       if (status .ne. 0) then
00294         call ftgerr(status,errmsg) 
00295         write(*,*) "error:",errmsg
00296         stop
00297       endif     
00298    !OK now we need to open a data file to retrieve the header information
00299       call  ftgiou(iunit,status)
00300       readwrite=0
00301       infn= "spot/bx1_wavesim.fits"
00302       call ftopen(iunit,infn,readwrite,blocksize,status)
00303    !how many key words do we have?
00304       call ftghsp(iunit,nkeys,nspace,status)
00305    ! Read each 80-character keyword record, and print it out.
00306       do i = 9, nkeys
00307           call ftgrec(iunit,i,record,status)
00308           call ftprec(unit,record,status)
00309       enddo
00310       call ftclos(iunit, status)
00311       call ftfiou(iunit, status)
00312 !OK HEADER DONE
00313 !     write the array to the FITS file
00314       group=1
00315       fpixel=1
00316       nelements=nx*ny*nz
00317       call ftppre(unit,group,fpixel,nelements,data,status)      
00318 !make next hdu
00319       call ftclos(unit, status)
00320       call ftfiou(unit, status) 
00321       end subroutine fw
00322 
00323 
00324