! $Id$
!
!  This module provide a way for users to specify custom initial
!  conditions.
!
!  The module provides a set of standard hooks into the Pencil Code
!  and currently allows the following customizations:
!
!   Description                               | Relevant function call
!  ------------------------------------------------------------------------
!   Initial condition registration            | register_initial_condition
!     (pre parameter read)                    |
!   Initial condition initialization          | initialize_initial_condition
!     (post parameter read)                   |
!                                             |
!   Initial condition for momentum            | initial_condition_uu
!   Initial condition for density             | initial_condition_lnrho
!   Initial condition for entropy             | initial_condition_ss
!   Initial condition for magnetic potential  | initial_condition_aa
!                                             |
!   Initial condition for all in one call     | initial_condition_all
!     (called last)                           |
!
!   And a similar subroutine for each module with an "init_XXX" call.
!   The subroutines are organized IN THE SAME ORDER THAT THEY ARE CALLED.
!   First uu, then lnrho, then ss, then aa, and so on.
!
!** 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.
!
!***************************************************************
!
! HOW TO USE THIS FILE
! --------------------
!
! Change the line above to
!    linitial_condition = .true.
! to enable use of custom initial conditions.
!
! The rest of this file may be used as a template for your own initial
! conditions. Simply fill out the prototypes for the features you want
! to use.
!
! Save the file with a meaningful name, e.g. mhs_equilibrium.f90, and
! place it in the $PENCIL_HOME/src/initial_condition directory. This
! path has been created to allow users to optionally check their
! contributions in to the Pencil Code SVN repository. This may be
! useful if you are working on/using an initial condition with
! somebody else or may require some assistance from one from the main
! Pencil Code team. HOWEVER, less general initial conditions should
! not go here (see below).
!
! You can also place initial condition files directly in the run
! directory. Simply create the folder 'initial_condition' at the same
! level as the *.in files and place an initial condition file there.
! With pc_setupsrc this file is linked automatically into the local
! src directory. This is the preferred method for initial conditions
! that are not very general.
!
! To use your additional initial condition code, edit the
! Makefile.local in the src directory under the run directory in which
! you wish to use your initial condition. Add a line that says e.g.
!
!    INITIAL_CONDITION =   initial_condition/mhs_equilibrium
!
! Here mhs_equilibrium is replaced by the filename of your new file,
! not including the .f90 extension.
!
! This module is based on Tony's special module.
!
module InitialCondition
!
  use Cdata
  use General, only: keep_compiler_quiet
  use Messages
!
  implicit none
!
  include '../initial_condition.h'
  real :: l0_d, l_d1, psipn0, psipn1, psi0, psi1, lnrho0_c
  real :: m0, r0_d, rs=1.0, apara, h_d, dsteep, ngamma, &
          rho0_c, cs0_c, Tc, rpos,psi0d,pbeta,amplaa=1.0
  character (len=labellen) :: initaa_dc='nothing'
!
  namelist /initial_condition_pars/ r0_d, apara,amplaa, & 
  h_d, rho0_c, dsteep, ngamma, cs0_c, Tc, rpos, pbeta,initaa_dc
!
  real :: gamma, gamma_m1, cp1

  contains
!***********************************************************************
    subroutine register_initial_condition()
!
!  Register variables associated with this module; likely none.
!
!  07-may-09/wlad: 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
!
      use EquationOfState, only: get_gamma_etc
      real, dimension (mx,my,mz,mfarray) :: f

      real :: cp
!
      call get_gamma_etc(gamma,cp)
      gamma_m1=gamma-1.
      cp1=1./cp

      call keep_compiler_quiet(f)
!
    endsubroutine initialize_initial_condition
!***********************************************************************
    subroutine initial_condition_uu(f)
!
!  Initialize the velocity field.
!
!  07-may-09/wlad: coded

!  /mayank
      real, dimension (mx,my,mz,mfarray), intent(inout) :: f
!      
      call keep_compiler_quiet(f)
!
    endsubroutine initial_condition_uu
!***********************************************************************

    subroutine initial_condition_all(f,profiles)
!
!  Initialize logarithmic density. init_lnrho will take care of
!  converting it to linear density if you use ldensity_nolog.
!
!  07-may-09/wlad: coded

!/mayank

      use Mpicomm, only: mpibcast
      use EquationOfState, only: cs0, cs20, cs2bot, cs2top, rho0, lnrho0
      use FArrayManager,   only: farray_use_global
      use Sub,             only: get_radial_distance, location_in_proc
      use SharedVariables, only: get_shared_variable
      
      real, dimension (mx,my,mz,mfarray), optional, intent(inout):: f
      real, dimension (:,:),              optional, intent(out)  :: profiles
      real, dimension (mx, my, mz) :: rho_d1, psi, rcut, masked, xmesh, tmasked, masked3,masked4
      real, dimension (mx) :: rr_sph, rr_cyl
      real, dimension (nx) :: rho_d, lnrho_c, psipn, l_d, vphi_d, lnT_d, aphi
      real, dimension (nx) :: dzphi_d, drphi_d, pg, dzrho_d, drrho_d, kphi, pg_c, pr,beta_c
      real, pointer :: g0
      real :: m00, rpn=2.0, pg0, drphi_d0, drrho_d0
      integer :: lpos, mpos, npos, procno,procnomax

! Some constants

      call get_shared_variable('g0',g0)
      lnrho0_c=log(rho0_c)
      rs=2*g0/c_light**2
      l0_d=sqrt(g0*(r0_d*rs)**3.0)/(r0_d*rs-rs)
      l_d1=l0_d*(2.0*rs/r0_d*rs)**apara
      psipn0=-g0/(r0_d*rs-rs)
      psipn1=-g0/rs
      psi0d=psipn0+1.0/(2.0-2.0*apara)*(l0_d/(r0_d*rs))**2.0
      psi1=psipn1+1.0/(2.0-2.0*apara)*(l_d1/(2.0*rs))**2.0
!      psi0=psipn0+1.0/(2.0-2.0*apara)*(l0_d/(r0_d*rs))**2.0
!    
      m00=(g0/G_Newton)
      do n=n1,n2
        do m=1,my
          xmesh(:,m,n)=x
          call get_radial_distance(rr_sph,rr_cyl)
          rcut(:,m,n)=rr_sph          
! Specific angular momentum
!        
          l_d=l0_d*(rr_cyl(l1:l2)/r0_d)**apara
!       
! Pseudo-Newtonian potential
!        
        psipn=-g0/(rr_sph(l1:l2)*rs-rs+0.01*rs)

! Gravitational Potential for the system
        
        psi(l1:l2,m,n)=psipn+1.0/(2.0-2.0*apara)*(l_d/(rr_cyl(l1:l2)*rs+0.01*rs))**2.0
        rho_d1(l1:l2,m,n)=rho0*(1.0-gamma*(psi(l1:l2,m,n)-psi0d)/(cs0**2.0*(ngamma+1.0)))**ngamma
        enddo
      enddo
      if (location_in_proc((/rpos,xyz0(2), 0.0/),lpos,mpos,npos)) then
          psi0=psi(lpos,mpos,npos)
!          procno=iprocmask2 = abs(XX) gt 4.0
          procno=iproc
!          print*, psi0mask2 = abs(XX) gt 4.0
      else
          procno=-100 
      endif

! The following is a hack as procno is not set for other processors
!
      call find_procno(procno,procnomax)
      call mpibcast(psi0,procnomax)
      print*, iproc, procnomax, psi0
      mask: where ((-(psi-psi0) .gt. 0.0) .and. (abs(xmesh) .gt. rpos))
               masked=1.0
             elsewhere
               masked=0.0
            end where mask
!
!      mask2: where ((-(psi-psi0) .gt. 0.0) .and. (abs(xmesh) .gt. rpos))
!               tmasked=0.0
!             elsewhere
!               tmasked=1.0
!            end where mask2
      mask3: where (rcut .gt. rpn)
               masked3=1.0
            elsewhere
               masked3=0.0
            end where mask3
!
      mask4: where ((rho_d1 .gt. 0.0) .and. (abs(xmesh) .gt. rpos))
               masked4=1.0
            elsewhere
               masked4=0.0
            end where mask4
      mask2: where ((rho_d1 .gt. 0.0) .and. (abs(xmesh) .gt. rpos))
               tmasked=0.0
            elsewhere
               tmasked=1.0
            end where mask2
!
      do n=n1,n2
        do m=m1,m2
          call get_radial_distance(rr_sph,rr_cyl)
       
! Specific angular momentum
!
          l_d=l0_d*(rr_cyl(l1:l2)/r0_d)**apara
!
! Pseudo-Newtonian potential
!
        psipn=-g0/(rr_sph(l1:l2)*rs-rs+0.01*rs)
        
!  Disk Density
                  
!        rho_d=masked(l1:l2,m,n)*exp(lnrho0)*(1.0-gamma*(psi(l1:l2,m,n)-psi0)/(cs0**2*(ngamma+1.0)))**ngamma
        rho_d=masked4(l1:l2,m,n)*rho0*(1.0-gamma*(psi(l1:l2,m,n)-psi0d)/(cs0**2.0*(ngamma+1.0)))**ngamma
!
! Corona Density
        lnrho_c=lnrho0_c-masked3(l1:l2,m,n)*(psipn-psipn1)*gamma/cs0_c**2.0
        
        pg=exp(lnrho0)**(-1.0/ngamma)*cs0**2.0/gamma*(rho_d)**(1.0+1.0/ngamma)
        pg_c=exp(lnrho_c)*cs0_c**2.0/gamma 
        pr=pg/(pg_c)
        f(l1:l2,m,n,ilnrho) = log(rho_d+exp(lnrho_c))
       ! f(l1:l2,m,n,ilnrho) = 
        !vphi_d=l_d*x(l1:l2)/rr_cyl(l1:l2)**2*masked(l1:l2,m,n)
        vphi_d=l_d*x(l1:l2)/(rr_cyl(l1:l2)*rs)**2*masked4(l1:l2,m,n)
        f(l1:l2,m,n,iuy) = vphi_d  
        lnT_d=(log(rho_d)-lnrho0)/ngamma+log(cs0**2.0/(gamma_m1/cp1))
        f(l1:l2,m,n,ilnTT)=log(exp(lnT_d)+Tc*tmasked(l1:l2,m,n))
        enddo
      enddo
      select case (initaa_dc)
        case ('kato')
          do n=n1,n2
            do m=m1,m2
              call get_radial_distance(rr_sph,rr_cyl)
        
              drphi_d0=g0/((r0_d*rs-rs)**2+(1e-3*rs)**2)*40.0*rs/sqrt((r0_d*rs)**2.0+&
                     (1e-3*rs)**2)-l0_d**2.0/((1-apara)*sqrt((r0_d*rs)**2+(1e-3*rs)**2)**3.0)
!
              pg0=exp(lnrho0)**(-1.0/ngamma)*cs0**2.0/gamma*(exp(lnrho0))**(1.0+1.0/ngamma)
!
              drrho_d0=-gamma*drphi_d0/(cs0**2.0*(ngamma+1.0))*exp(lnrho0)*ngamma*(1.0/ &
                  (cs0**2*(ngamma+1.0)))**(ngamma-1.0)
!        
              kphi=sqrt(2.0*mu0*pg0/(pbeta*(drrho_d0)**2.0))
              f(l1:l2,m,n,iay)=amplaa*x(l1:l2)*exp(f(l1:l2,m,n,ilnrho)) !*masked(l1:l2,m,n)
            enddo
          enddo
        case ('loop')
          do n=n1,n2
            do m=m1,m2
              call get_radial_distance(rr_sph,rr_cyl)
              aphi=amplaa*rpos**2.0/(rpos**2.0+rr_sph(l1:l2)**2)**1.5*(1+15.0*rpos**2.0*&
                x(l1:l2)**2.0/(8.0*(rpos**2.0+rr_sph(l1:l2)**2)**2.0))
              f(l1:l2,m,n,iay)=aphi*x(l1:l2)  !*masked(l1:l2,m,n)
            enddo
          enddo
        case ('nothing')
!         Do nothing
        case default
          call fatal_error('initial_condition_all','no such initaa_dc: '//trim(initaa_dc))
      endselect 
!
    endsubroutine initial_condition_all
!***********************************************************************
    subroutine find_procno(procno,procnomax)
!
!  Find the absolute maximum of the velocity.
!
!  13-sep-2023/piyali: adapted from find_umax
!
      use Mpicomm, only: mpiallreduce_max, MPI_COMM_WORLD
!
      integer, intent(in) :: procno
      integer, intent(out) :: procnomax
!
      integer :: procno1
!
!  Find the maximum.
!
      procno1 = procno
      call mpiallreduce_max(procno1, procnomax, comm=MPI_COMM_WORLD)
!
    endsubroutine find_procno
!***********************************************************************
    subroutine initial_condition_ss(f)
!
!  Initialize entropy.
!
!  07-may-09/wlad: coded
! 
      real, dimension (mx,my,mz,mfarray), intent(inout) :: f
! 
         call keep_compiler_quiet(f)
! 
    endsubroutine initial_condition_ss
!***********************************************************************
    subroutine initial_condition_aa(f)
!
!  Initialize the magnetic vector potential.
!
!  07-may-09/wlad: coded
!
      real, dimension (mx,my,mz,mfarray), intent(inout) :: f

!
      call keep_compiler_quiet(f)
!
    endsubroutine initial_condition_aa
!***********************************************************************
    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