program initFG implicit none integer(kind=kind(1000)) :: i, j, k integer(kind=kind(1000)), parameter :: nx = 30, ny= 20, nz = 25 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 :: rg = 3.75d-1, a=1.d-1, q = -1.d0, B0 = 1.d0, & Bt = 9.d0*B0, B1=0.3, eps=1.d-10 write(*,*) ' ' write(*,*) '*****************************************' write(*,*) '* *' write(*,*) '* Fp7 European Network Soteria *' write(*,*) '* *' write(*,*) '* Flare/C.M.E. Initiation Challenge *' write(*,*) '* *' write(*,*) '* Fan, Y. & Gibson, S.E. ApJ 609 (2004) *' write(*,*) '* *' write(*,*) '* Initial conditions *' write(*,*) '* *' write(*,*) '* Suggested compiler: g95 *' write(*,*) '* *' write(*,*) '*****************************************' write(*,*) ' ' call alloc(1) call gridFG call alloc(2) call fieldsFG call writeFG call dealloc write(*,*) ' ' write(*,*) '*****************************************' write(*,*) '* *' write(*,*) '* Fp7 European Network Soteria *' write(*,*) '* *' write(*,*) '* Flare/C.M.E. Initiation Challenge *' write(*,*) '* *' write(*,*) '* Fan, Y. & Gibson, S.E. ApJ 609 (2004) *' write(*,*) '* *' write(*,*) '* End of the program *' write(*,*) '* *' write(*,*) '* *' 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 gridFG Lx = 1.5d0 Ly = 1.d0 Lz = 1.25d0 xi = 0.d0 yi = 0.d0 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 gridFG subroutine fieldsFG real(kind=kind(1.d0)) :: xc, yc, zc, & r, cth, sth, cfi, sfi, omega, & dA, p0mpth, p0mpr, bfi, & br, bth, kappa, elle, alpha xc = Lx/2 yc = Ly/2 zc = 0.d0 write(*,*) ' force-free arcade i.c. ' write(*,*) ' ' do i = 1, nx do j = 1, ny do k = 1, nz r = sqrt((x(i)-xc)*(x(i)-xc)+(y(j)-yc)*(y(j)-yc)+(z(k)-zc)*(z(k)-zc)) cth = (y(j)-yc)/(r+eps) sth = sqrt(1.d0-cth*cth) cfi = (x(i)-xc)/(r*sth+eps) sfi = (z(k)-zc)/(r*sth+eps) omega = sqrt(r*r+rg*rg-2.d0*r*rg*sth*sth) dA = -q*Bt*exp(-omega*omega/(a*a))*omega p0mpth = -2.d0/omega*r*rg*sth*cth p0mpr = 1.d0/omega*(r-rg*sth*sth) bfi = a*Bt/(r*sth+eps)*exp(-omega*omega/(a*a)) br = dA*p0mpth/(r*r*sth+eps) bth = -dA*p0mpr/(r*sth+eps) ro(i,j,k) = 1.d0 pr(i,j,k) = 1.d0 kappa = acos(-1.d0)/Ly elle = 2.d0/Lz if(elle <= kappa) then alpha = sqrt(kappa*kappa-elle*elle) else elle = kappa alpha = 0.d0 endif bx(i,j,k) = -bfi*sfi+cth*cfi*bth+sth*cfi*br+alpha/kappa*B1*cos(kappa*(y(j)-yc))*exp(-elle*(z(k)-zc)); by(i,j,k) = -sth*bth+cth*br-exp(-(z(k)-zc)/(a*a)/16.0)+elle/kappa*B1*cos(kappa*(y(j)-yc))*exp(-elle*(z(k)-zc)); bz(i,j,k) = cfi*bfi+cth*sfi*bth+sth*sfi*br-B1*sin(kappa*(y(j)-yc))*exp(-elle*(z(k)-zc)); end do end do end do return end subroutine fieldsFG subroutine writeFG character(len=2) :: snx, sny, snz write(*,*) ' writing data for Paraview ' write(*,*) ' ' write(*,*) ' check the "endianity"! ' write(*,*) ' ' write(snx,'(i2)') nx-1 write(sny,'(i2)') 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 writeFG end program initFG