program initB implicit none integer(kind=kind(1000)) :: i, j, k integer(kind=kind(1000)), parameter :: nx = 100, ny= 160, nz = 80 real(kind=kind(0.d0)) :: xi, xf, Lx, yi, yf, Ly, zi, zf, Lz real(kind=kind(0.d0)), allocatable :: x(:),y(:),z(:) real(kind=kind(0.d0)), allocatable :: bx(:,:,:),by(:,:,:),bz(:,:,:),ro(:,:,:),pr(:,:,:) real(kind=kind(0.d0)), parameter :: f = 4.d-1, z0 = 1.d+1, c0 = 5.d-1, c1 = 5.d-1, & r1 = 2.2d+1, L1 = 3.d0, rn = 3.d+1, ry = 8.d+1,& eps = 1.d0, Bc = 1.d0 write(*,*) ' ' write(*,*) '**************************************' write(*,*) '* *' write(*,*) '* Fp7 European Network Soteria *' write(*,*) '* *' write(*,*) '* Flare/C.M.E. Initiation Challenge *' write(*,*) '* *' write(*,*) '* Birn, J et al ApJ 645 (2006) *' write(*,*) '* *' write(*,*) '* Initial conditions *' write(*,*) '* *' write(*,*) '* Suggested compiler: g95 *' write(*,*) '* *' write(*,*) '**************************************' write(*,*) ' ' call alloc(1) call gridB call alloc(2) call fieldsB call writeB call dealloc write(*,*) ' ' write(*,*) '**************************************' write(*,*) '* *' write(*,*) '* Fp7 European Network Soteria *' write(*,*) '* *' write(*,*) '* Flare/C.M.E. Initiation Challenge *' write(*,*) '* *' write(*,*) '* Birn, J et al ApJ 645 (2006) *' write(*,*) '* *' write(*,*) '* End of the program *' write(*,*) '* *' write(*,*) '**************************************' write(*,*) ' ' contains subroutine alloc(nf) integer(kind=kind(10)) :: nf select case(nf) case(1) ! grid write(*,*) ' allocating grid variables ' write(*,*) ' ' allocate(x(nx),y(ny),z(nz)) x(:) = 0.d0 y(:) = 0.d0 z(:) = 0.d0 case(2) ! fields write(*,*) ' allocating field quantites ' write(*,*) ' ' allocate(bx(nx,ny,nz),by(nx,ny,nz),bz(nx,ny,nz),ro(nx,ny,nz),pr(nx,ny,nz)) bx(:,:,:) = 0.d0 by(:,:,:) = 0.d0 bz(:,:,:) = 0.d0 ro(:,:,:) = 0.d0 pr(:,:,:) = 0.d0 end select return end subroutine alloc subroutine dealloc write(*,*) ' deallocating grid variables ' write(*,*) ' ' deallocate(x,y,z) deallocate(bx,by,bz,ro,pr) return end subroutine dealloc subroutine gridB Lx = 5.d+1 Ly = 8.d+1 Lz = 4.d+1 xi = -2.5d+1 yi = -4.0d+1 zi = 0.d0 xf = xi+Lx yf = yi+Ly zf = zi+Lz do i = 1, nx x(i) = xi+Lx/(nx-1)*(i-1) end do do j = 1, ny y(j) = yi+Ly/(ny-1)*(j-1) end do do k = 1, nz z(k) = zi+Lz/(nz-1)*(k-1) end do return end subroutine gridB subroutine fieldsB real(kind=kind(1.d0)) :: r, rd, sphi, cphi, B0, & B1, B0p, Bhat, Bhatp, xi, Br, Bphi do i = 1, nx do j = 1, ny do k = 1, nz r = y(j)*y(j)+(z(k)+z0)*(z(k)+z0) rd = r+1.d-3 sphi = y(j)/r cphi = (z(k)+z0)/r B0 = Bc*sqrt(c0+(1.-(r-z0)/(ry-z0))*(1.-(r-z0)/(ry-z0))+c1/(1.+((r-r1)/L1)*((r-r1)/L1))) if(r > ry) B0 = Bc*sqrt(c0+c1/(1.+((r-r1)/L1*(r-r1)/L1))) B1 = Bc*sqrt(c0+(1.-(rd-z0)/(ry-z0))*(1.-(rd-z0)/(ry-z0))+c1/(1.+(((rd-r1)/L1)*((rd-r1)/L1)))) if(rd > ry) B1 = Bc*sqrt(c0+c1/(1.+((rd-r1)/L1)*((rd-r1)/L1))) B0p = (B1-B0)/1.d-3 Bhat = B0*sqrt(f+(1.-f)*(rn/r)*(rn/r)); Bhatp = -B0*(1-f)*rn*rn/(r**3*sqrt(f+(1.-f)*(rn/r)*(rn/r))); Bhatp = Bhatp+Bhat*B0p; xi = r*Bhat*x(i)/(eps*rn*rn); Br = -Bhat*tanh(xi); Bphi = sqrt(1.-f)*(rn/r)*B0/cosh(xi); bx(i,j,k) = eps*(rn/r)*(xi*tanh(xi)*(rn*Bhatp/Bhat+rn/r)-B0p/B0); by(i,j,k) = cphi*Bphi+sphi*Br; bz(i,j,k) = -sphi*Bphi+cphi*Br; ro(i,j,k) = 1.d0+f*1.d0 pr(i,j,k) = f*B0*B0/(cosh(xi)*cosh(xi)) end do end do end do return end subroutine fieldsB subroutine writeB character(len=2) :: snx, snz character(len=3) :: sny write(*,*) ' writing data for Paraview ' write(*,*) ' ' write(*,*) ' check the "endianity"! ' write(*,*) ' ' write(snx,'(i2)') nx-1 write(sny,'(i3)') ny-1 write(snz,'(i2)') nz-1 open(100,file='init_birn.vtr') write(100,'(a21)') '' write(100,'(a72)') '' write(100,'(a47)') '' write(100,'(a32)') '' write(100,'(a40)') '' write(100,'(a51)') '' write(100,1000) bx write(100,'(a12)') '' write(100,'(a51)') '' write(100,1000) by write(100,'(a12)') '' write(100,'(a51)') '' write(100,1000) bz write(100,'(a12)') '' write(100,'(a51)') '' write(100,1000) ro write(100,'(a12)') '' write(100,'(a51)') '' write(100,1001) pr write(100,'(a12)') '' write(100,'(a12)'),'' write(100,'(a13)'),'' write(100,'(a50)'),'' do i = 1, nx write(100,1002),x(i) end do write(100,'(a12)'),'' write(100,'(a50)'),'' do j = 1, ny write(100,1002),y(j) end do write(100,'(a12)'),'' write(100,'(a50)'),'' do k = 1, nz write(100,1002),z(k) end do write(100,'(a12)') '' write(100,'(a14)') '' write(100,'(a8)') '' write(100,'(a18)') '' write(100,'(a10)') '' close(100) 1000 format(5(1pe12.4)) 1001 format(f5.2) 1002 format(f7.2) return end subroutine writeB end program initB