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

00001 # 1 "makecoeff.f90"
00002 # 1 "<built-in>"
00003 # 1 "<command line>"
00004 # 1 "makecoeff.f90"
00005 program makecoeff
00006  !#define threed 1
00007   ! program to make coefficients for helioseismology runs
00008   implicit none
00009   include 'comments.f90'
00010   include 'param.f90'
00011   real B0(nx,ny,nz_max,3),U0(nx,ny,nz_max,3),bztemp(nx,ny,nz_max,3)
00012   real pressure(nx,ny,nz_max),gamma(nx,ny,nz_max),rho(nx,ny,nz_max)
00013   real g_3d(nx,ny,nz_max)
00014   real vtemp(nx,ny,nz_max,3)
00015   real coeff(nx,ny,nz_max,ncoeff)
00016   ! real*4 coeffout(nx,ny,nz,ncoeff)
00017   real temp(nx,ny,nz_max),fac,data(nx,ny,nz_max)
00018   real global(nx,ny,nz_ghostglobal+2)
00019   real global_out(nx,ny,nz_ghostglobal)
00020 
00021   character*80 fn,fnb
00022   integer k,k1,i,j,n_rec
00023 
00024   interface
00025      function grad(datain)
00026        include 'param.f90'
00027        real datain(nx,ny,nz_max),grad(nx,ny,nz_max,3)
00028      end function grad
00029   end interface
00030   interface
00031      function ave(datain)
00032        include 'param.f90'
00033        real datain(nx,ny,nz_max),ave(nx,ny,nz_max,3)
00034      end function ave
00035   end interface
00036 
00037 !INITIALIZE MPI stuff
00038  call mpi_init(ierr)
00039 
00040 ! Give each node its number and the total number of nodes.
00041 ! Node 0 will handle all io.
00042 ! Also mark the top and bottom boxes.
00043   !write(*,*) "Main"
00044   call mpi_comm_size (MPI_COMM_WORLD, nprocs, ierr)
00045   call mpi_comm_rank (MPI_COMM_WORLD, my_rank, ierr)
00046 
00047 
00048 !first read in the files
00049 
00050   if(my_rank .eq. 0) then
00051      fn="spot/bx1_wavesim.fits"
00052      call fr(fn,global,nx,ny,nz_ghostglobal+2)
00053   endif
00054   call distributeall(global,temp)
00055 
00056   b0(:,:,:,1)=temp
00057   if(my_rank .eq. 0) then
00058      fn="spot/by1_wavesim.fits"
00059      call fr(fn,global,nx,ny,nz_ghostglobal+2)
00060   endif
00061   call distributeall(global,temp)
00062   b0(:,:,:,2)=temp
00063 
00064   if(my_rank .eq. 0) then
00065      fn="spot/bz1_wavesim.fits"
00066      call fr(fn,global,nx,ny,nz_ghostglobal+2)
00067   endif
00068   call distributeall(global,temp)
00069   b0(:,:,:,3)=temp
00070 
00071   if(my_rank .eq. 0) then
00072      fn="spot/pressure_wavesim.fits"
00073      call fr(fn,global,nx,ny,nz_ghostglobal+2)
00074   endif
00075   call distributeall(global,temp)
00076   pressure=temp
00077 
00078   if(my_rank .eq. 0) then
00079      fn="spot/rho_wavesim.fits"
00080      call fr(fn,global,nx,ny,nz_ghostglobal+2)
00081   endif
00082   call distributeall(global,temp)
00083   rho=temp
00084 
00085 
00086  if(my_rank .eq. 0) then
00087      fn="spot/gamma1_wavesim.fits"
00088      call fr(fn,global,nx,ny,nz_ghostglobal+2)
00089   endif
00090   call distributeall(global,temp)
00091   gamma=temp
00092 
00093   if(my_rank .eq. 0) then
00094      fn="spot/ux1_wavesim.fits"
00095      call fr(fn,global,nx,ny,nz_ghostglobal+2)
00096   endif
00097   call distributeall(global,temp)
00098   u0(:,:,:,1)=temp
00099 
00100   if(my_rank .eq. 0) then
00101      fn="spot/uy1_wavesim.fits"
00102      call fr(fn,global,nx,ny,nz_ghostglobal+2)
00103   endif
00104   call distributeall(global,temp)
00105   u0(:,:,:,2)=temp
00106 
00107   if(my_rank .eq. 0) then
00108      fn="spot/uz1_wavesim.fits"
00109      call fr(fn,global,nx,ny,nz_ghostglobal+2)
00110   endif
00111   call distributeall(global,temp)
00112   u0(:,:,:,3)=temp
00113 
00114   if(my_rank .eq. 0) then
00115      fn="spot/g_wavesim.fits"
00116      call fr(fn,global,nx,ny,nz_ghostglobal+2)
00117   endif
00118   call distributeall(global,temp)
00119   g_3d(:,:,:)=-temp
00120   !g_3d=27400
00121 
00122 
00123 
00124   write(*,*) "DISTRIBUTED"
00125 
00126 
00127 
00128 
00129 
00130 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00131  ! B0=0
00132 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00133   ! next calculate the coefficient array
00134   coeff(:,:,:,1)=rho
00135  ! coeff(:,:,:,2:4)=grad(rho)
00136   vtemp=grad(alog(rho))
00137   coeff(:,:,:,2)=rho*vtemp(:,:,:,1)
00138   coeff(:,:,:,3)=rho*vtemp(:,:,:,2)
00139   coeff(:,:,:,4)=rho*vtemp(:,:,:,3)
00140   coeff(:,:,:,5)= Gamma * pressure / rho
00141 
00142 ! coeff(:,:,:,9:11)=Grad(pressure)
00143   vtemp=Grad(alog(pressure))
00144   coeff(:,:,:,9)=pressure*vtemp(:,:,:,1)
00145   coeff(:,:,:,10)=pressure*vtemp(:,:,:,2)
00146   coeff(:,:,:,11)=pressure*vtemp(:,:,:,3)
00147 
00148 
00149   ! coeff(:,:,:,9:10)=0
00150   ! coeff(:,:,:,11)=-rho *g
00151   coeff(:,:,:,6)= Gamma*pressure*coeff(:,:,:,2)/rho -coeff(:,:,:,9)
00152   coeff(:,:,:,7)= Gamma*pressure*coeff(:,:,:,3)/rho -coeff(:,:,:,10)
00153   coeff(:,:,:,8)= Gamma*pressure*coeff(:,:,:,4)/rho -coeff(:,:,:,11)
00154   ! coeff(:,:,:,6)= -coeff(:,:,:,9)
00155   ! coeff(:,:,:,7)= -coeff(:,:,:,10)
00156   ! coeff(:,:,:,8)= -coeff(:,:,:,11)
00157   coeff(:,:,:,12)= 1. / rho
00158   coeff(:,:,:,13)= g_3d / rho
00159   coeff(:,:,:,14)= g_3d
00160   !magnetic field stuff
00161   coeff(:,:,:,15)= B0(:,:,:,1)
00162   coeff(:,:,:,16)= B0(:,:,:,2)
00163   coeff(:,:,:,17)= B0(:,:,:,3)
00164   vtemp=Grad(B0(:,:,:,1))
00165   coeff(:,:,:,19)= vtemp(:,:,:,3)
00166   coeff(:,:,:,20)= -vtemp(:,:,:,2)
00167   vtemp=Grad(B0(:,:,:,2))
00168   coeff(:,:,:,18)= -vtemp(:,:,:,3)
00169   coeff(:,:,:,20)= coeff(:,:,:,20)+ vtemp(:,:,:,1)
00170   vtemp=Grad(B0(:,:,:,3))
00171   coeff(:,:,:,18)= coeff(:,:,:,18)+vtemp(:,:,:,2)
00172   coeff(:,:,:,19)= coeff(:,:,:,19)-vtemp(:,:,:,1)
00173   !velocity field stuff
00174   coeff(:,:,:,21)= U0(:,:,:,1)
00175   coeff(:,:,:,22)= U0(:,:,:,2)
00176   coeff(:,:,:,23)= U0(:,:,:,3)
00177   vtemp=Grad(U0(:,:,:,1))
00178   coeff(:,:,:,24)= vtemp(:,:,:,1)
00179   coeff(:,:,:,27)= vtemp(:,:,:,2)
00180   coeff(:,:,:,30)= vtemp(:,:,:,3)
00181   vtemp=Grad(U0(:,:,:,2))
00182   coeff(:,:,:,25)= vtemp(:,:,:,1)
00183   coeff(:,:,:,28)= vtemp(:,:,:,2)
00184   coeff(:,:,:,31)= vtemp(:,:,:,3)
00185   vtemp=Grad(U0(:,:,:,3))
00186   coeff(:,:,:,26)= vtemp(:,:,:,1)
00187   coeff(:,:,:,29)= vtemp(:,:,:,2)
00188   coeff(:,:,:,32)= vtemp(:,:,:,3)
00189 
00190   coeff(:,:,:,33)= U0(:,:,:,1)/rho
00191   coeff(:,:,:,34)= U0(:,:,:,2)/rho
00192   coeff(:,:,:,35)= U0(:,:,:,3)/rho
00193   coeff(:,:,:,36)=Gamma*pressure
00194   !---------------------
00195   Bztemp=b0(:,:,:,:)
00196   where (abs(Bztemp) .lt. 1e-9) Bztemp=1e-9
00197   !---------------------
00198   coeff(:,:,:,37)=B0(:,:,:,1)/Bztemp(:,:,:,3)
00199   coeff(:,:,:,38)=B0(:,:,:,2)/Bztemp(:,:,:,3)
00200   vtemp(:,:,:,1)=pressure+&
00201        & (b0(:,:,:,1)**2+b0(:,:,:,2)**2-b0(:,:,:,3)**2)/(8.*pi)
00202   coeff(:,:,:,39:41)=Grad(vtemp(:,:,:,1))
00203 
00204   !---------------------
00205   Bztemp=ave(b0(:,:,:,3))
00206   where (abs(Bztemp) .lt. 1e-6) Bztemp=1e-6
00207   !---------------------
00208 
00209   vtemp(:,:,:,1)=B0(:,:,:,1)*B0(:,:,:,3)
00210   coeff(:,:,:,42:44)=Grad(vtemp(:,:,:,1))
00211   coeff(:,:,:,42)=coeff(:,:,:,42)/bztemp(:,:,:,1)
00212   coeff(:,:,:,43)=coeff(:,:,:,43)/bztemp(:,:,:,2)
00213   coeff(:,:,:,44)=coeff(:,:,:,44)/bztemp(:,:,:,3)
00214   vtemp(:,:,:,1)=B0(:,:,:,2)*B0(:,:,:,3)
00215   coeff(:,:,:,45:47)=Grad(vtemp(:,:,:,1))
00216   coeff(:,:,:,45)=coeff(:,:,:,45)/bztemp(:,:,:,1)
00217   coeff(:,:,:,46)=coeff(:,:,:,46)/bztemp(:,:,:,2)
00218   coeff(:,:,:,47)=coeff(:,:,:,47)/bztemp(:,:,:,3)
00219   coeff(:,:,:,48)=1./bztemp(:,:,:,3)
00220 
00221   vtemp(:,:,:,1)=pressure+&
00222        & (b0(:,:,:,1)**2+b0(:,:,:,2)**2-b0(:,:,:,3)**2)/(8.*pi)*0
00223 
00224            coeff(:,:,:,39:41)=Grad(vtemp(:,:,:,1))
00225 !do the output
00226   fnb="coeff_wavesim.fits"
00227   fn=fnb
00228   data=coeff(:,:,:,1)
00229   call collectall(global,data)
00230 
00231   global_out=global(:,:,2:nz_ghostglobal+1)
00232   if(my_rank .eq. 0) then
00233      call fw(fn,global_out,nx,ny,nz_ghostglobal,0)
00234   endif
00235   do k=2,ncoeff
00236  n_rec=k-1
00237  if(n_rec .le. 10) write(fn,'(I1)') n_rec-1
00238  if(n_rec .gt. 10) write(fn,'(I2)') n_rec-1
00239         fn=trim(fnb)
00240         data=coeff(:,:,:,k)
00241         call collectall(global,data)
00242 ! write(*,*) global(1,1,nz_ghostglobal+1),data(1,1,nz_ghostglobal+1)
00243         global_out=global(:,:,2:nz_ghostglobal+1)
00244         if(my_rank .eq. 0) then
00245            call fw(fn,global_out,nx,ny,nz_ghostglobal,k-1)
00246         endif
00247   enddo
00248    call mpi_finalize(ierr)
00249 
00250 end program makecoeff
00251 !!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00252 function grad(datain)
00253   implicit none
00254   include 'param.f90'
00255   real datain(nx,ny,nz_max)
00256   real grad( nx,ny,nz_max,3)
00257   grad(2:nx-1,:,:,1)=(datain(3:nx,:,:)-datain(1:nx-2,:,:))/(2.*dx)
00258 
00259 
00260 
00261   grad(:,:,2:nz_max-1,3)=(datain(:,:,3:nz_max)&
00262        & -datain(:,:,1:nz_max-2))/(2.*dz)
00263 ! grad(:,:,:,3)=0
00264 
00265 ! x,y and z directions are periodic
00266   grad( 1,:,:,1)=(datain(2,:,:)-datain(nx ,:,:))/(2.*dx)
00267   grad(nx,:,:,1)=(datain(1,:,:)-datain(nx-1,:,:))/(2.*dx)
00268 
00269 
00270 
00271 
00272 
00273   grad(:,:, 1,3)=grad(:,:,2,3)
00274   grad(:,:,nz_max,3)=grad(:,:,nz_max-1,3)
00275 
00276 
00277 end function grad
00278 !!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00279 function ave(datain)
00280   implicit none
00281   include 'param.f90'
00282   real datain(nx,ny,nz_max)
00283   real ave( nx,ny,nz_max,3)
00284   ave(2:nx-1,:,:,1)=(datain(3:nx,:,:)+datain(1:nx-2,:,:))/(2.)
00285 
00286 
00287 
00288   ave(:,:,2:nz+1,3)=(datain(:,:,3:nz+2)+datain(:,:,1:nz))/(2.)
00289 
00290 ! x,y and z directions are periodic
00291   ave( 1,:,:,1)=(datain(2,:,:)+datain(nx ,:,:))/(2.)
00292   ave(nx,:,:,1)=(datain(1,:,:)+datain(nx-1,:,:))/(2.)
00293 
00294 
00295 
00296 
00297 
00298   ave(:,:, 1,3)=ave(:,:,2,3)
00299   ave(:,:,nz+2,3)=ave(:,:,nz+1,3)
00300 end function ave
00301 !========================================================================================================================
00302       subroutine fr(fn,data,nx,ny,nz)
00303       implicit none
00304       integer nx,ny,nz,blocksize
00305       character*80 fn
00306       real data(nx,ny,nz)
00307       integer readwrite,group,fpixel,nelements,status,unit
00308       logical er
00309       character*30 errmsg
00310 
00311       readwrite=0
00312       blocksize=1
00313       status=0
00314       call ftgiou(unit,status)
00315       call ftopen(unit,fn,readwrite,blocksize,status)
00316       group=1
00317       fpixel=1
00318       nelements=nx*ny*nz
00319       call ftgpve(unit,group,fpixel,nelements,0.,data,er,status)
00320       call ftclos(unit, status)
00321       call ftfiou(unit, status)
00322 
00323       if(status .ne. 0) then
00324         call ftgerr(status,errmsg)
00325         write(*,*) "reading: ",fn
00326  write(*,*) "readin error: ",status,errmsg
00327  stop
00328       endif
00329       end subroutine fr
00330 !========================================================================================================================
00331       subroutine fw(fn,data,nx,ny,nz,n_rec)
00332       implicit none
00333       character*80 fn,record,infn
00334       character*30 errmsg
00335       real data(nx,ny,nz),alpha,b_0
00336       integer status,unit,blocksize,bitpix,naxis,naxes(3)
00337       integer i,j,group,fpixel,nelements,nx,ny,nz
00338       integer n,n_rec,readwrite,nkeys,nspace,iunit
00339       logical simple,extend
00340 
00341       status=0
00342       write(*,*) n_rec,minval(data),maxval(data),fn
00343 ! Get an unused Logical Unit Number to use to create the FITS file
00344       call ftgiou(unit,status)
00345       call ftgiou(iunit,status)
00346       blocksize=1
00347       readwrite=1
00348       if(n_rec .eq. 0) then
00349  call ftinit(unit,fn,blocksize,status)
00350       else
00351  call ftnopn(unit,fn,blocksize,status)
00352  call ftcrhd(unit,status)
00353       endif
00354       simple=.true.
00355       bitpix=-32
00356       naxis=3
00357       naxes(1)=nx
00358       naxes(2)=ny
00359       naxes(3)=nz
00360       extend=.true.
00361       call ftphpr(unit,simple,bitpix,naxis,naxes,0,1,extend,status)
00362 
00363       if (status .ne. 0) then
00364  call ftgerr(status,errmsg)
00365  write(*,*) "error:",errmsg
00366  stop
00367       endif
00368    !OK now we need to open a data file to retrieve the header information
00369       call ftgiou(iunit,status)
00370       readwrite=0
00371       infn= "spot/bx1_wavesim.fits"
00372       call ftopen(iunit,infn,readwrite,blocksize,status)
00373    !how many key words do we have?
00374       call ftghsp(iunit,nkeys,nspace,status)
00375    ! Read each 80-character keyword record, and print it out.
00376       do i = 9, nkeys
00377           call ftgrec(iunit,i,record,status)
00378           call ftprec(unit,record,status)
00379       enddo
00380       call ftclos(iunit, status)
00381       call ftfiou(iunit, status)
00382 !OK HEADER DONE
00383 ! write the array to the FITS file
00384       group=1
00385       fpixel=1
00386       nelements=nx*ny*nz
00387       call ftppre(unit,group,fpixel,nelements,data,status)
00388 !make next hdu
00389       call ftclos(unit, status)
00390       call ftfiou(unit, status)
00391       end subroutine fw
00392 
00393 !------------------------------------------------------------
00394 subroutine collectall(rdinfull,rdin)
00395   implicit none
00396   include 'param.f90'
00397 
00398   real rdinfull(nx,ny,nz_ghostglobal+2)
00399   real rdin(nx,ny,nz_max)
00400   real temp(nx,ny,nz_max)
00401   integer*8 i
00402 
00403 !Note uses different breakup to main program.
00404   write(*,*),"nprocs",nprocs
00405 
00406   dn=(nz_global+2)/(nprocs)
00407 !excludes ghosts
00408   do i=1,nprocs
00409      start(i)=(i-1)*dn+1+2
00410      ends(i)=(i)*dn+2
00411   enddo
00412   start(1)=3
00413   ends(nprocs)=nz_global+4
00414 
00415   nz=(ends(my_rank+1)-start(my_rank+1)+1)+4
00416   write(*,*) nz,nz_ghostglobal+2
00417 
00418   if(my_rank .ne. 0) then
00419      call mpi_send(rdin,nx*ny*nz_max,MPI_REAL,&
00420           & 0,2,MPI_COMM_WORLD,ierr)
00421   endif
00422   if(my_rank .eq. 0) then
00423      !first deal with node 0
00424      rdinfull(:,:,1:nz)=rdin(:,:,1:nz)
00425      do i=1,nprocs-1
00426         call mpi_recv(temp,nx*ny*nz_max,MPI_REAL,&
00427              & i,2,MPI_COMM_WORLD,istatus,ierr)
00428         rdinfull(:,:,start(i+1):ends(i+1)+2)=temp(:,:,3:ends(i+1)-start(i+1)+5)
00429      enddo
00430   endif
00431 end subroutine collectall
00432 
00433 
00434 !====================================================
00435 
00436 
00437 
00438 
00439 
00440 !===============================================================
00441 subroutine distributeall(rdinfull,rdin)
00442   implicit none
00443   include 'param.f90'
00444   real rdinfull(nx,ny,nz_ghostglobal+2)
00445   real rdin(nx,ny,nz_max)
00446   real temp(nx,ny,nz_max)
00447   integer*8 i
00448 !Note uses different breakup to main program.
00449   dn=(nz_global)/(nprocs)
00450 !excludes ghosts
00451   do i=1,nprocs
00452      start(i)=(i-1)*dn+1+2
00453      ends(i)=(i)*dn+2
00454   enddo
00455   start(1)=3
00456   ends(nprocs)=nz_global+4
00457 
00458 
00459 
00460   if(my_rank .eq. 0) then
00461      !first deal with node 0
00462      rdin(:,:,1:nz_max)=rdinfull(:,:,1:nz_max)
00463      do i=1,nprocs-1
00464         temp(:,:,1:ends(i+1)-start(i+1)+5)=rdinfull(:,:,start(i+1)-2:ends(i+1)+2)
00465         call mpi_send(temp,nx*ny*nz_max,MPI_REAL,&
00466              & i,9,MPI_COMM_WORLD,ierr)
00467      enddo
00468   endif
00469   if(my_rank .ne. 0) then
00470      call mpi_recv(rdin,nx*ny*nz_max,MPI_REAL,&
00471           & 0,9,MPI_COMM_WORLD,istatus,ierr)
00472   endif
00473 end subroutine distributeall