! $Id$
!
! Initial condition for a stably tratified solar atmosphere
! from the convection zone to the corona, with a Omega shaped
! magnetic field loop.
!
!** AUTOMATIC CPARAM.INC GENERATION ****************************
! Declare (for generation of cparam.inc) the number of f array
! variables and auxiliary variables added by this module
!
! CPARAM logical, parameter :: linitial_condition = .true.
!
!***************************************************************

module InitialCondition
!
  use Cparam
  use Cdata
  use General, only: keep_compiler_quiet
  use Messages

!
  implicit none
!
  include '../initial_condition.h'
!
  ! g = surface gravity
  ! mu = mean atomic weight
  ! z_ss = photosphereic height
  ! z_tr = start of the transition region
  ! z_cr = start of the corona
  ! p_ss = pressure ont he solar surface
  ! rho_ss = density on the solar surface
  ! T_ss = temperature on the solar surface
  ! T_cr = temperature in the corona
  ! H_ss = photosphere scale height 
  ! H_cr = corona scale height
  ! xi = mu*g/R
  ! x0, z0 = center of the magnetic flux tube
  ! B0 = magnetic field strength
  ! d = diameter of the flux tube with Gaussian profile
  ! alpha = twist parameter for the magnetic field
  ! lam = exponent for Gaussian y-profile for rising flux tube

  real :: z_ss = 0, z_tr = 10, z_cr = 20
  real :: p_ss = 1, rho_ss = 1, T_ss = 1, T_cr = 150
  real :: H_ss = 1, H_cr = 150
  real :: xi = 1
  real :: x_0 = 0, z_0 = 0, B0 = 1, d = 1, alpha = 0.5, lam = 1

  namelist /initial_condition_pars/ &
          z_ss, z_tr, z_cr, &
          p_ss, rho_ss, T_ss, T_cr, &
          H_ss, H_cr, &
          xi, &
          x_0, z_0, B0, d, alpha, lam
!
  contains
!***********************************************************************
    subroutine register_initial_condition()
!
!  Register variables associated with this module; likely none.
!
!  20-may-11/simon: coded
!
      if (lroot) call svn_id( &
           "$Id$")
!
    endsubroutine register_initial_condition
!***********************************************************************
    subroutine initialize_initial_condition(f)
!
!  Initialize any module variables which are parameter dependent.
!
!  07-may-09/wlad: coded
!
      real, dimension (mx,my,mz,mfarray) :: f
!
      call keep_compiler_quiet(f)
!
    endsubroutine initialize_initial_condition
!***********************************************************************
    subroutine initial_condition_uu(f)
!
!  Initialize the velocity field.
!
!  07-may-09/wlad: coded
!
      real, dimension (mx,my,mz,mfarray), intent(inout) :: f
!
!  SAMPLE IMPLEMENTATION
!
      call keep_compiler_quiet(f)
!
    endsubroutine initial_condition_uu
!***********************************************************************
    subroutine initial_condition_lnrho(f)
!
!  Initialize logarithmic density. init_lnrho 
!  will take care of converting it to linear 
!  density if you use ldensity_nolog
!
!  07-sep-11/simon: coded
!
      use SharedVariables
      use EquationOfState
      
      real, dimension (mx,my,mz,mfarray), intent(inout) :: f
      integer :: l, m, n
      real :: r, p_ex, p_b
!
    do n = n1, n2, 1
      do m = m1, m2, 1
        do l = l1, l2, 1
          if (z(n) < z_ss) then
            f(l,m,n,ilnrho) = (1/(gamma-1))*log((1-z(n)*(gamma-1)/gamma))
            f(l,m,n,ilnTT) = log(1 - z(n)*(gamma-1)/gamma)
          elseif (z(n) < z_tr) then
            f(l,m,n,ilnrho) = -z(n)
            f(l,m,n,ilnTT) = 0
          else if (z(n) < z_cr) then
            f(l,m,n,ilnrho) = -z_tr - (z(n) - z_tr)/z_tr*log(T_cr) + z_tr/log(T_cr)*(1/T_cr**((z(n)-z_tr)/z_tr) - 1)
            f(l,m,n,ilnTT) = (z(n)-z_tr)/z_tr * log(T_cr)
          else
            f(l,m,n,ilnrho) = -z_tr - log(T_cr) + z_tr/log(T_cr)*(1/T_cr-1) - (z(n)-z_cr)/T_cr
            f(l,m,n,ilnTT) = log(T_cr)
          end if
          ! Introduce density deficit.
          r = sqrt((x(l)-x_0)**2 + (z(n)-z_0)**2)
          if (r < d) then
            p_ex = B0**2 * exp(-r**2/d**2)*((alpha*d)**2 - 2 - 2*(alpha*r)**2)/4
            ! 0.4 = (cv - cp); gamma = cp/cv; cv = 1
            p_b = 0.4*exp(f(l,m,n,ilnrho) + f(l,m,n,ilnTT))
            f(l,m,n,ilnrho) = log(exp(f(l,m,n,ilnrho)) + p_ex/p_b*exp(-y(m)**2/lam**2))
          end if
        enddo
      enddo
    enddo    

    endsubroutine initial_condition_lnrho
!********************************************************************
    subroutine initial_condition_aa(f)
!
!  Initialize entropy.
!
!  07-sep-11/simon: coded
!
!       use SharedVariables
!       use EquationOfState
      use Poisson
      use Sub

      real, dimension (mx,my,mz,mfarray) :: f      
      integer :: l, j, ju
      real :: r, theta
      real, dimension (nx,ny,nz,3) :: jj, tmpJ  ! This is phi for poisson.f90
!
      do n = n1, n2, 1
        do m = m1, m2, 1
          do l = l1, l2, 1
            theta = atan2((z(n)-z_0), (x(l)-x_0))
            r = sqrt((x(l)-x_0)**2 + (z(n)-z_0)**2)            
            f(l,m,n,iay) = B0*exp(-r**2/d**2)
            f(l,m,n,iax) = -alpha*r*f(l,m,n,iay) * sin(theta)
            f(l,m,n,iaz) =  alpha*r*f(l,m,n,iay) * cos(theta)
          enddo
        enddo
      enddo
!
!  Compute curl(B) = J for the Poisson solver
      do m=m1,m2
         do n=n1,n2
            call curl(f,iaa,jj(:,m-nghost,n-nghost,:))
         enddo
      enddo
      tmpJ = -jj
!  Use the Poisson solver to solve \nabla^2 A = -J for A
      do j=1,3
        call inverse_laplacian(tmpJ(:,:,:,j))
      enddo
      
!  Overwrite the f-array with the correct vector potential A
      do j=1,3
          ju=iaa-1+j
          f(l1:l2,m1:m2,n1:n2,ju) = tmpJ(:,:,:,j)
      enddo
!       
    endsubroutine initial_condition_aa
!***********************************************************************
    subroutine initial_condition_ss(f)
!
!  Initialize entropy.
!
!  07-sep-11/simon: coded
!
      use SharedVariables
      use EquationOfState

      real, dimension (mx,my,mz,mfarray), intent(inout) :: f      
      integer :: l, m, n
      real :: log_T_b, log_T_0
      real :: r, z_tilde, log_T_ex  ! working variables for smoothing
!
    ! initialize the temperature of the medium and the bubble
    do n = n1, n2, 1
      do m = m1, m2, 1
        do l = l1, l2, 1
          if (z(n) < z_ss) then
            f(l,m,n,ilnTT) = log(1 - z(n)*(gamma-1)/gamma)
          elseif (z(n) < z_tr) then
            f(l,m,n,ilnTT) = 0
          else if (z(n) < z_cr) then
            f(l,m,n,ilnTT) = (z(n)-z_tr)/z_tr * log(T_cr)
          else
            f(l,m,n,ilnTT) = log(T_cr)
          end if
        enddo
      enddo
    enddo
!      
    endsubroutine initial_condition_ss
!***********************************************************************
    subroutine initial_condition_lncc(f)
!
!  Initialize entropy.
!
!  07-sep-11/simon: coded
!
      real, dimension (mx,my,mz,mfarray), intent(inout) :: f      
      integer :: l, m, n
!
!     ! initialize the temperature of the medium and the bubble
!     if (passive_scalar == 1) then
!         do n = n1, n2, 1
!             do m = m1, m2, 1
!                 do l = l1, l2, 1
!                     ! check if this point lies in the bubble
!                     if (((x(l) - x_b)**2 + (y(m) - y_b)**2 + (z(n) - z_b)**2) .le. r_b**2) then
!                         f(l,m,n,ilncc) = 1.
!                     else
!                         f(l,m,n,ilncc) = 0.
!                     endif
!                 enddo
!             enddo
!         enddo
!     endif
!      
    endsubroutine initial_condition_lncc
!***********************************************************************
    subroutine read_initial_condition_pars(iostat)
!
      use File_io, only: parallel_unit
!
      integer, intent(out) :: iostat
!
      read(parallel_unit, NML=initial_condition_pars, IOSTAT=iostat)
!
    endsubroutine read_initial_condition_pars
!***********************************************************************
    subroutine write_initial_condition_pars(unit)
!
      integer, intent(in) :: unit
!
      write(unit, NML=initial_condition_pars)
!
    endsubroutine write_initial_condition_pars
!***********************************************************************
!
!********************************************************************
!************        DO NOT DELETE THE FOLLOWING       **************
!********************************************************************
!**  This is an automatically generated include file that creates  **
!**  copies dummy routines from noinitial_condition.f90 for any    **
!**  InitialCondition routines not implemented in this file        **
!**                                                                **
    include '../initial_condition_dummies.inc'
!********************************************************************
endmodule InitialCondition