
c********************************************************

      subroutine mhdfldsslab(t1)
 
      implicit double precision (a-h,o-z)
 
      include 'maxdef.h'

      double precision t1
      double precision bx1d(NZMAX+2), by1d(NZMAX+2), bz1d 
      common /fieldslab/ bx1d, by1d, bz1d

c     
c      Box lengths in terms of the turbulence correlation length
c
c
      integer NX, NY, NZ
      double precision lenx,leny,lenz
      common /box/ lenx, leny, lenz, NX, NY, NZ

      double precision delx,dely,delz
      common /grids/ delx,dely,delz
c
      logical debug
      common /debug/ debug

      double precision  b0, dbob
      common /bfield/  b0, dbob

      integer nkmin,nkmax
      double precision alpha,lambdax,lambday
      double precision eslab,fracx
      common /slabmodel/ eslab,fracx,lambdax,lambday,alpha,
     &                    nkmin,nkmax

      double precision bkx2,bky2

      integer ikind
      common /ikind/ ikind

      integer printb
      common /printb/ printb
      
      integer iseed_1,iseed_2,iseed_3
      common /seed/ iseed_1,iseed_2,iseed_3

      character*80 filename
c
      double precision TWOPI
      parameter ( TWOPI = 2.0d0*3.14159265358979323846d0 )
c
      integer ndim(1)
      integer i,nkmaxp,nktot,iseed1,iseed2,nk,kk

      double precision ee,efinal,enorm,bxmag,bymag
      double precision dbx2,dby2,db2,phi1,phi2

      double precision efinalx, efinaly,eex,eey,enormx,enormy      

      logical zloadedslab,zloaded2d
      common /zloaded/ zloadedslab,zloaded2d

c      data zloadedslab /.FALSE./
c      save zloadedslab      

      double precision kn2, skhalfx,skhalfy,k,boxlength
      

      if ( ikind .eq. 0 ) then
       if ( zloadedslab ) return
      endif
c      write(*,*)'lambda',lambda 
c      write(*,*)'lambdax',lambdax 
c      write(*,*)'lambday',lambday 
      ndim(1) = NZ
      bx1d(1) = 0.0d0
      bx1d(2) = 0.0d0
      by1d(1) = 0.0d0
      by1d(2) = 0.0d0
      boxlength = lenz*lambdax
c      boxlength = lenz
      
c
c     The physical maximum wavenumber set by array constraints. 
c     WARNING: NZ *must* be a power of 2

      nktot = NZ/2 
c      

c     efinal is the magnetic field energy that we want
c 	efinal = 0.0
      efinal = (dbob*dbob*b0*b0)*eslab


c     fracx is dbx^2/db^2 (the value is  between 0 and 1) 

      bxmag = dsqrt(fracx)
      bymag = dsqrt(1-fracx)
      efinalx = efinal*(fracx)
      efinaly = efinal*(1-fracx)
     
c
      if ( .not. zloadedslab ) then
         zloadedslab = .TRUE.

c     set bxk and byk at k=0
       
        do nk = 1, nktot

          i = 2*nk+1
          k = TWOPI*dble(nk)/boxlength
          bx1d(i)   = 0.0d0
          bx1d(i+1) = 0.0d0
          by1d(i)   = 0.0d0
          by1d(i+1) = 0.0d0
          skhalfx = 0.0d0
          skhalfy = 0.0d0

c  
c         This is spectrum for a model 
c                |b(k)|^2 \approx 1/[1+(k*lambda)^2]^(alpha/2) 
c
          if(nk .ge. nkmin .AND. nk .le. nkmax) then
             skhalfx = (1+(k*lambdax)**2)**(-alpha/4.0)
             skhalfy = (1+(k*lambdax)**2)**(-alpha/4.0)
c	      skhalfx = k**(-5/6)
c	      skhalfy = k**(-5/6)
             phi1 = TWOPI*dble(ran1(iseed_1))
             phi2 = TWOPI*dble(ran1(iseed_1))

c         Use gasdev() from Numerical Recipe to generate random number
c         This routine returns a normally distributed deviate with 
c         zero mean and unit variance, using ran1(idum) as a source of 
c         uniform deviates
           
             bx1d(i)   = skhalfx*dcos(phi1)
             bx1d(i+1) = skhalfx*dsin(phi1)
             by1d(i)   = skhalfy*dcos(phi2)
             by1d(i+1) = skhalfy*dsin(phi2)
          endif 
          

        enddo
      endif 
c      if(printb .eq. 1) then
c         filename = 'bkx1d_b.dat'
c         call writetofile(filename,bx1d,NZ+2)
c         filename = 'bky1d_b.dat'
c         call writetofile(filename,by1d,NZ+2)
c      endif
c    
c     Compute raw energy
c     
      ee = 0.0
      eex = 0.0
      eey = 0.0
      bkx2 = 0.0
      bky2 = 0.0

      do i=1, NZ+2
         eex = eex + bx1d(i)*bx1d(i)
      enddo

      do i=1, NZ+2
         eey = eey + by1d(i)*by1d(i)
      enddo

      enormx = dsqrt(efinalx/(2*eex))
      enormy = dsqrt(efinaly/(2*eey))
c      write(*,*) 'enorm in slab part =',enormx,enormy

C     Normalize magnetic field energy = efinal

      do i=1,NZ+2

         bx1d(i) = bx1d(i)*enormx
         by1d(i) = by1d(i)*enormy

      enddo

      do i=1, NZ+2
         bkx2 = bkx2 + bx1d(i)*bx1d(i)
         bky2 = bky2 + by1d(i)*by1d(i)
         
      enddo

c      write(*,*) 'bkx2=',bkx2*2
c      write(*,*) 'bky2=',bky2*2
c      write(*,*) 'bk2=',bky2*2+bkx2*2

      if(printb .eq. 1) then
         filename = 'bkx1d_a.dat'
         call writetofile(filename,bx1d,NZ+2)
         filename = 'bky1d_a.dat'
         call writetofile(filename,by1d,NZ+2)
      endif

   
     
c
c     Transform to z space      
c

      call four2(bx1d,ndim,1,-1,-1)
      call four2(by1d,ndim,1,-1,-1)

c     Factor 2 needs to complete Inverse Fourier Transform
    
      do i=1,NZ+2
         bx1d(i) = bx1d(i)*2.0
         by1d(i) = by1d(i)*2.0
      enddo
      
      bz1d = b0

c
c     Check normalization in z space
c       
      dbx2  = 0.0
      dby2  = 0.0
      db2   = 0.0
      do i=1,NZ
         dbx2 = dbx2 + bx1d(i)*bx1d(i)
         dby2 = dby2 + by1d(i)*by1d(i)
         db2  = db2  + bx1d(i)*bx1d(i)+by1d(i)*by1d(i)
      enddo
      dbx2 = dbx2/dble(NZ)
      dby2 = dby2/dble(NZ) 
      db2  = db2/dble(NZ)
      
c      if(printb .eq. 1) then
c         filename = 'bx1d.dat'
c         call writetofile(filename,bx1d,NZ+2)
c         filename = 'by1d.dat'
c         call writetofile(filename,by1d,NZ+2)
c      endif

c      Output file from c_file_write give wrong name
c      filename = 'bx1d'
c      call c_file_write( 1, NZ, filename, bx1d)
c      filename = 'by1d'
c      call c_file_write( 1, NZ, filename, by1d)


c      write(*,*) 'end of mhdfldsslab()'
      return
      end






c********************************************************
      subroutine mhdflds2d(t)
      implicit double precision (a-h,o-z)
      double precision t 

      include 'maxdef.h'
      
      
      dimension a(NNMAX)
      integer n(2),nhalf(2)
      double precision bx,by
      common /field2d/ bx(NNMAX),by(NNMAX) 

      double precision efinal2d,phi1,phi2

      double precision lenx,leny,lenz
      common /box/ lenx,leny,lenz,NX,NY,NZ
c     real lenx,leny,lenz
      
      double precision delx,dely,delz
      common /grids/ delx,dely,delz
c     
      double precision  b0, dbob
      common /bfield/  b0, dbob

      
      integer nkmin,nkmax
      double precision alpha,lambda_parx,lambda_pary
      double precision eslab,fracx
      common /slabmodel/ eslab,fracx,lambda_parx,lambda_pary,alpha,
     &                    nkmin,nkmax
      double precision kmax,kmin
      integer nxhalf,nyhalf
      double precision alpha2,lambda_per
      double precision e2d
      double precision xi
      common /twodmodel/ e2d,lambda_per,alpha2,
     &            kmin,kmax,nxhalf,nyhalf,xi


      integer ikind
      common /ikind/ ikind

      integer printb
      common /printb/ printb

      integer nn	
      
      logical zloadedslab,zloaded2d
      common /zloaded/ zloadedslab,zloaded2d
      
      double precision boxlengthx,boxlengthy

      double precision zk1,zk2,rzk,zk,deltakx,deltaky,skhalf
      double precision b2tot,a2tot,b2totx,b2totxn,a2totx,btot2
      logical debug
      common /debug/debug
 
      character*80 filename
        
      parameter ( TWOPI = 2.0d0*3.14159265358979323846d0 )

      integer iseed_1,iseed_2,iseed_3
      common /seed/ iseed_1,iseed_2,iseed_3
        
c     integer strlen
        



        nn = (NX+2)*NY
c        write(*,*)'nn=',nn

        nnhalf = (nxhalf+2)*nyhalf
        
	
	if ( ikind .eq. 0 ) then
           if ( zloaded2d ) return
        endif
c        write(*,*) 'nxhalf=',nxhalf
c        write(*,*) 'nyhalf=',nyhalf
c        write(*,*) 'kmin=',kmin
c        write(*,*) 'kmax=',kmax
c        write(*,*) 'xi',xi

        if(.not. zloaded2d) zloaded2d = .true.
        
        n(1) = NX
        n(2) = NY
        
        nhalf(1) = nxhalf
        nhalf(2) = nyhalf
        
        
c     target energy in the 2D modes is E2D
        
        efinal2d = e2d*(dbob*b0)**2 
c        write(*,*) 'efinal2d=',efinal2d

c        write(*,*) 'lambda_parx',lambda_parx
 
        boxlengthx = lenx*lambda_parx
        boxlengthy = leny*lambda_parx
   
c        boxlengthx = lenx
c        boxlengthy = leny

        deltakx = TWOPI/boxlengthx
        deltaky = TWOPI/boxlengthy


c        write(*,*)'nkmin2d,nkmax2d',nkmin2d,nkmax2d
c        write(*,*)'kmin,kmax',kmin,kmax
        
        
        zk1= kmin !1.0e-10   lower inclusive k-limit for spectrum
        zk2= kmax !600.0d0   upper k-limit)
        

c        write(*,*) 'zk1=',zk1,'zk2=',zk2        
C
	do  i=1,nn
           a(i)=0.d0
	enddo
        
        nexcited = 0
           write(*,*)'seed2d', iseed_2

        do  j = 4,nnhalf,2
c           a(j) = 0.d0
c           a(j-1) = 0.d0
           loc= j 
           call kloc(kx,ky,nhalf,loc)

c          write(*,*) 'kx=',kx,'  ky=',ky
           zkx  = (sqrt(xi))*deltakx*dble(kx)
	   zky  = (1./sqrt(xi))*deltaky*dble(ky)
           if (j .eq. 6) write(*,*) zkx,zky
 
           zksq = zkx*zkx + zky*zky
           zk = dsqrt(zksq)
           if(zk .ge. zk1 .and. zk .le. zk2) then 
c             skhalf = (1+(zk*lambda_per)**2)**(-alpha2/4.0)
c             rzk = (1.d0/zk)**1.5d0
              phi = TWOPI*dble(ran1(iseed_2))
              
              call lloc(kx,ky,n,loc2)

c             a(loc2) =   rzk*skhalf*dsin(phi)
c             a(loc2-1) = rzk*skhalf*dcos(phi)
             a(loc2) = dsin(phi)*(1+(zk*lambda_per)**2)**(-7.0/6.0)
             a(loc2-1) = dcos(phi)*(1+(zk*lambda_per)**2)**(-7.0/6.0)
c              write(*,*) 'zkx=',zkx,'zky=',zky
c              write(*,*) 'a(loc2)=',a(loc2),'a(loc2-1)',a(loc2-1)

           endif
               

       enddo       


c... enforce reality condition

        call realf(a,n)
       

c-- normalize energy

        b2tot = 0.d0
        do j = 1,nn
           loc = j
           call kloc(kx,ky,n,loc)
           if(kx.gt.0.or.ky.gt.0) then
              zkx = deltakx*dble(kx)
              zky = deltaky*dble(ky)
              b2tot =b2tot+dble(zkx*zkx+zky*zky)*dble(a(j)*a(j))
           endif
        enddo
c.....normalization factors
        fact = dsqrt(efinal2d/(2.0d0*b2tot))  
        do j = 1,nn
           a(j) = fact* a(j)
        enddo
c-------normalization complete

c--    print |a(k)|^2
 
        if (printb .eq. 1) then

           open(13,file='ak2.dat')
              do j=2,nn,2
                 loc = j
                 call kloc(kx,ky,n,loc)
                 write(13,*) kx,ky,a(loc)*a(loc)+a(loc-1)*a(loc-1)
              enddo
           close(13)
           
        endif

                                    
c     calculate  B(k)

       do j=2,nn,2
         loc = j
         call kloc(kx,ky,n,loc)
	 zkx     = deltakx*dble(kx)
	 zky     = deltaky*dble(ky)
c        write(*,*)'zky=',zky  
	 bx(j-1) = zky*a(j  )
	 bx(j  ) = - zky*a(j-1)
	 by(j-1) = - zkx*a(j  )
	 by(j  ) = zkx*a(j-1)
         
c           
       enddo

       call realf(bx,n)
       call realf(by,n)



        b2tot = 0.d0
        a2tot = 0.d0
        btot2 = 0.0d0
        do j = 1,nn

          loc = j
          call kloc(kx,ky,n,loc)
          if(kx.gt.0.or.ky.gt.0) then
            a2tot =a2tot+dble(a(j)*a(j))
	    zkx   = deltakx*dble(kx)
	    zky   = deltaky*dble(ky)
            b2tot =b2tot+dble(zkx*zkx+zky*zky)*dble(a(j)*a(j))
            btot2 = btot2 + bx(j)*bx(j)+by(j)*by(j)
          endif
        enddo

c        write(*,*)'b2tot before FFT',b2tot
c        write(*,*)'btot2 before FFT',btot2


c     have Bx(k) --> Bx(x), By(k) --> By(k)

       
c       write(*,*)'n',n

       call four2(bx,n,2,-1,-1)
       call four2(by,n,2,-1,-1)
       call four2(a,n,2,-1,-1)

       do l=1,nn
	  bx(l) = bx(l) * 2.d0
	  by(l) = by(l) * 2.d0
	  a(l) = a(l) * 2.d0
       enddo

c <b**2> in x-space
      b2totxn=0.0d0
      do i=1,nn
	loc=i
	call kloc(kx,ky,n,loc)
	if(kx.gt.0.or.ky.gt.0)then	
	   b2totxn=b2totxn+dble(bx(i)*bx(i)+by(i)*by(i))
	endif
      enddo

		
      b2totxn=b2totxn/dble(NX*NY)		      
c      write(*,*)'b2totxn=',b2totxn
    
      if(printb .eq. 1) then
         filename = 'bx2d.dat'
         call writetofile(filename,bx,nn)
         filename = 'by2d.dat'
         call writetofile(filename,by,nn)
         filename = 'axy.dat'
         call writetofile(filename,a,nn)
      endif
 

ccccc  Output file from c_file_write give wrong name
c      filename = 'bx2d'
c      call c_file_write( 1, NX*NY,filename, bx)
c      filename = 'by2d'
c      call c_file_write( 1, NX*NY,filename, by)
c      write(*,*) 'end of mhdflds2d()' 

c      write(*,*)'end mhsflds2d'
      RETURN   
      END

c********************************************************
      subroutine mhdflds2d2(t)
      implicit double precision (a-h,o-z)
      double precision t 

      include 'maxdef.h'
      
      
      dimension a(NNMAX)
      integer n(2),nhalf(2)
      double precision bx,by
      common /field2d/ bx(NNMAX),by(NNMAX) 

      double precision efinal2d,phi1,phi2

      double precision lenx,leny,lenz
      common /box/ lenx,leny,lenz,NX,NY,NZ
c     real lenx,leny,lenz
      
      double precision delx,dely,delz
      common /grids/ delx,dely,delz
c     
      double precision  b0, dbob
      common /bfield/  b0, dbob

      
      integer nkmin,nkmax
      double precision alpha,lambda_parx,lambda_pary
      double precision eslab,fracx
      common /slabmodel/ eslab,fracx,lambda_parx,lambda_pary,alpha,
     &                    nkmin,nkmax
      double precision kmax,kmin
      integer nxhalf,nyhalf
      double precision alpha2,lambda_per
      double precision e2d
      double precision xi
      common /twodmodel/ e2d,lambda_per,alpha2,
     &            kmin,kmax,nxhalf,nyhalf,xi


      integer ikind
      common /ikind/ ikind

      integer printb
      common /printb/ printb

      integer nn	
      
      logical zloadedslab,zloaded2d
      common /zloaded/ zloadedslab,zloaded2d
      
      double precision boxlengthx,boxlengthy

      double precision zk1,zk2,rzk,zk,deltakx,deltaky,skhalf
      double precision b2tot,a2tot,b2totx,b2totxn,a2totx,btot2
      logical debug
      common /debug/debug
 
      character*80 filename
        
      parameter ( TWOPI = 2.0d0*3.14159265358979323846d0 )

      integer iseed_1,iseed_2,iseed_3
      common /seed/ iseed_1,iseed_2,iseed_3
        
c     integer strlen
        

	if ( ikind .eq. 0 ) then
           if ( zloaded2d ) return
        endif
c        write(*,*) 'nxhalf=',nxhalf
c        write(*,*) 'nyhalf=',nyhalf
c        write(*,*) 'kmin=',kmin
c        write(*,*) 'kmax=',kmax
c        write(*,*) 'xi',xi

        if(.not. zloaded2d) zloaded2d = .true.
        
	open(11,FILE='bx2d.dat')
	do i=1,NX*NY
	   read(11,*) bx(i)
	enddo
	close(11)

	open(12,FILE='by2d.dat')
	do i=1,NX*NY
	   read(12,*) by(i)
	enddo
	close(12)

      RETURN   
      END

        SUBROUTINE REALF(OM,N)
	implicit double precision (a-h,o-z)
        DIMENSION OM(1),N(2)
        NY = N(2)
        DO 1 J = 1,NY/2+1
        KY = J - 1
        CALL LLOC(0,KY,N,LOC)
        CALL LLOC(0,-KY,N,NLOC)
        OM(NLOC) = -1.d0 * OM(LOC)
        OM(NLOC-1) = OM(LOC-1)
 1      CONTINUE
        RETURN
        END


      SUBROUTINE writetofile(flnm,a,n)
     
      character*80 flnm
      double precision a(n)
      integer n,i

      open(1,file='flnm')
      do i = 1,n
         write(1,500) a(i)
      enddo
 500  format(F25.15)
      close(1)
      return
      end

