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