! $Id: newton_cooling.f90,v 1.2 2015/12/11 02:56:37 wlyra Exp $

!** AUTOMATIC CPARAM.INC GENERATION ****************************
! Declare (for generation of special_dummies.inc) the number of f array
! variables and auxiliary variables added by this module
!
! CPARAM logical, parameter :: lspecial = .true.
!
! MVAR CONTRIBUTION 0
! MAUX CONTRIBUTION 2
!
!***************************************************************

!-------------------------------------------------------------------
!
! HOW TO USE THIS FILE
! --------------------
!
! The rest of this file may be used as a template for your own
! special module.  Lines which are double commented are intended
! as examples of code.  Simply fill out the prototypes for the 
! features you want to use.
!
! Save the file with a meaningful name, eg. geo_kws.f90 and place
! it in the $PENCIL_HOME/src/special directory.  This path has
! been created to allow users ot optionally check their contributions
! in to the Pencil-Code CVS repository.  This may be useful if you
! are working on/using the additional physics with somebodyelse or 
! may require some assistance from one of the main Pencil-Code team.
!
! To use your additional physics code edit the Makefile.local in
! the src directory under the run directory in which you wish to
! use your additional physics.  Add a line with all the module 
! selections to say something like:
!
!    SPECIAL=special/nstar
!
! Where nstar it replaced by the filename of your new module
! upto and not including the .f90
!
!--------------------------------------------------------------------
!
!  This method solves for flux-limited diffusion, as per
!  
!    Levermore & Pomraning 1981, ApJ, 248, 321
!    Kley 1989, A&A, 208, 98
!  
!  11-dec-14/wlad: coded
!
module Special
!
  use Cdata
  use General, only: keep_compiler_quiet
  use Messages
!
  implicit none
!
  include '../special.h'
!
  real :: taucool_floor = impossible
  logical :: laddheatingrate=.true.
!
  namelist /special_init_pars/ taucool_floor
!
  namelist /special_run_pars/ laddheatingrate,taucool_floor
!
  type InternalPencils
     real, dimension(nx) :: tau,taucool,taucool1
  endtype InternalPencils
!
  real, dimension(nx,nz,0:nprocy-1) :: tau_column
  real, dimension(nx,ny,nz) :: dtau
!
  type (InternalPencils) :: q
!
  integer :: ikappar,itau
!
  integer :: idiag_kappam=0,idiag_kappamax=0,idiag_kappamin=0
  integer :: idiag_taum=0,idiag_taumax=0,idiag_taumin=0
  integer :: idiag_taucoolm=0,idiag_taucoolmax=0,idiag_taucoolmin=0
  integer :: idiag_dts=0
!
  real :: gamma, gamma1, gamma_m1, cp1, cv1

  contains
!***********************************************************************
    subroutine register_special
!
!  Configure pre-initialised (i.e. before parameter read) variables 
!  which should be know to be able to evaluate
!
!  14-jul-09/wlad: coded
!
      use FArrayManager, only: farray_register_auxiliary
!
      if (lroot) call svn_id( &
           "$Id: newton_cooling.f90,v 1.2 2015/12/11 02:56:37 wlyra Exp $")
!
      call farray_register_auxiliary('kappar',ikappar)
      call farray_register_auxiliary('tau',itau)
!
      if (lroot) then
        open(4,file=trim(datadir)//'/index_special.pro',status='replace')
        write(4,*) 'tau    ',itau
        write(4,*) 'kappar ',ikappar
        close(4)
      endif
!
    endsubroutine register_special
!***********************************************************************
    subroutine initialize_special(f)
!
!  called by run.f90 after reading parameters, but before the time loop
!
!  14-jul-09/wlad: coded
!
      use EquationOfState, only: get_gamma_etc
!
      real, dimension (mx,my,mz,mfarray) :: f
      real :: cp, cv
!
      if (.not.lentropy) &
        call fatal_error('initialize_special','newton cooling requires entropy module')

      call get_gamma_etc(gamma,cp,cv)

      gamma1=1./gamma; gamma_m1=gamma-1.
      cp1=1./cp; cv1=1./cv
!
      call keep_compiler_quiet(f)
!
    endsubroutine initialize_special
!***********************************************************************
    subroutine init_special(f)
!
      real, dimension (mx,my,mz,mfarray) :: f
!
      if (lroot) print*,'calling special before boundary'
!MR: should be needed on start only if kappar is written in snapshot.
      call special_before_boundary(f)
!
    endsubroutine  init_special
!***********************************************************************
    subroutine special_before_boundary(f)
!
!  This subroutine calculates the full potential due to the turbulence.
!
!  03-oct-12/wlad: coded
!
      use EquationOfState, only: cs20,rho0
      use Sub, only: grad,dot
!
      real, dimension(mx,my,mz,mfarray), intent(inout) :: f
      real, dimension(nx) :: rho,TT
      real :: TT0,rho01,lnTT0,kappa_cgs,TTdim,rhodim
      integer :: i
!
      lnTT0=log(cs20*cp1/gamma_m1)
      TT0=exp(lnTT0)
      rho01=1./rho0
!
      do n=n1,n2; do m=m1,m2
!
        if (ldensity_nolog) then
          rho=f(l1:l2,m,n,irho)
        else
          rho=exp(f(l1:l2,m,n,ilnrho))
        endif
        TT = TT0 * (rho*rho01)**gamma_m1 * exp(f(l1:l2,m,n,iss)*cv1) !MR: better use generic eoscalc here!
!
! lnrho = log(f(i,m,n,irho))
! lnTT = lnTT0 + cv1*f(i,m,n,iss) + gamma_m1*(lnrho-lnrho0)
! TT_csg = exp(lnTT) * temperature_unit_cgs
!
        do i=1,nx 
           ! spits out opacity in cgs units
           TTdim=TT(i)*unit_temperature
           rhodim=rho(i)*unit_density
           call calc_opacity(TTdim,rhodim,kappa_cgs)
!
! kappa_code = kappa_cgs/ unit_kappa
! unit_kappa = 1/(unit_density*unit_length)
!
           f(i+nghost,m,n,ikappar) = kappa_cgs * (unit_density*unit_length)
!
        enddo
!
! Also precompute dtau for optical depth calculation. 
!
        dtau(:,m-nghost,n-nghost) = f(l1:l2,m,n,ikappar)*rho*x(l1:l2)/dy_1(m)
!
      enddo;enddo
!
!  Calculate now opacity
!
      call integrate_optical_depth(f)
!
!  Boundaries for limiter (gradient will be needed)
!
    endsubroutine special_before_boundary
!***********************************************************************
    subroutine integrate_optical_depth(f)
!
      real, dimension(mx,my,mz,mfarray), intent(inout) :: f
      real, dimension(my) :: tau_below_local,tau_above_local
      real, dimension(my) :: tau_below,tau_above
      real :: tau_above_sum,tau_below_sum
      integer :: ipy_loc,i
!
      if (lmpicomm) call calc_column_tau(f)
!
!  For every point, get stuff above and below -- assume meriodional to start with
!
      do n=n1,n2; do i=l1,l2
          !if tht < 90, gas parcel is above midplane. Sum stuff that is above is (lower thetas)
          !if tht > 90, gas parcel is below midplane. Sum stuff that is below it (higher thetas)
          !In fact, for all, take both all that is above and all that is below, take optical depth
          !  as the smaller of the two numbers. 
!
        do m=m1,m2
          tau_above_local(m) = sum(dtau(i-nghost,m1-nghost:m-nghost,n-nghost))
        enddo
        tau_above_sum=0.
        if (ipy/=0) then
          do ipy_loc=0,ipy-1
            tau_above_sum=tau_above_sum+tau_column(i-nghost,n-nghost,ipy_loc)
          enddo
        endif
        do m=m1,m2
          tau_above(m)=tau_above_sum+tau_above_local(m)
        enddo
!
        do m=m2,m1,-1
          tau_below_local(m) = sum(dtau(i-nghost,m-nghost:m2-nghost,n-nghost))
        enddo
        tau_below_sum=0.
        if (lmpicomm) then
          if (ipy/=nprocy-1) then
            do ipy_loc=ipy+1,nprocy-1
              tau_below_sum=tau_below_sum+tau_column(i-nghost,n-nghost,ipy_loc)
            enddo
          endif
        endif
        do m=m1,m2
          tau_below(m)=tau_below_sum+tau_below_local(m)
        enddo
!
        do m=m1,m2
          !if (f(i,m,n,ikappar) /= impossible) then !T >> T9
             f(i,m,n,itau)=min(tau_below(m),tau_above(m))
          !else
          !   !force it to cool! 
          !   f(i,m,n,itau)=1.
          !endif
        enddo
!
      enddo;enddo
!
    endsubroutine integrate_optical_depth
!***********************************************************************
    subroutine calc_column_tau(f)
!
!  Integrate the density vertically (meridionally): we need column 
!  density for the X-ray and cosmic ray attenuation factor.
!
      use Mpicomm,only: mpibcast_real,mpisend_real,mpirecv_real
!
      real, dimension(mx,my,mz,mfarray) :: f
      integer :: i,n,j
      real, dimension(nx,nz) :: tau_column_local,tmp
      real, dimension(nx,nz,0:nprocy-1) :: tau_column_yroot,tmp3
      integer, parameter :: y1tag=117
      integer, parameter :: y2tag=217
      integer :: partner,py,broadcaster,collector
!
      do n=1,nz; do i=1,nx
        tau_column_local(i,n)=sum(dtau(i,:,n))
      enddo; enddo
!      
!  Collect all y data
! 
      collector = ipx + ipz*nprocxy !root in y
      if (iproc == collector) then
        do py = 0,nprocy-1 
          partner = ipx + py*nprocx + ipz*nprocxy
          if (iproc == partner) then
            !data is local
            tau_column_yroot(:,:,py)=tau_column_local
          else
            call mpirecv_real(tmp,(/nx,nz/),partner,y1tag)
            tau_column_yroot(:,:,py)=tmp
          endif
        enddo
      else
        call mpisend_real(tau_column_local,(/nx,nz/),collector,y1tag)
      endif
!
!  Broadcast the data back to all y-processors
!
      broadcaster=ipx + ipz*nprocxy !root in y
      if (iproc == broadcaster) then
        do py=0,nprocy-1
          partner = ipx + py*nprocx + ipz*nprocxy
          if (iproc == partner) then 
            !data is local
            tau_column=tau_column_yroot
          else
            ! send to partner
            call mpisend_real(tau_column_yroot,(/nx,nz,nprocy/),partner,y2tag)
          endif
        enddo
      else
        call mpirecv_real(tau_column,(/nx,nz,nprocy/),broadcaster,y2tag)
      endif
!
    endsubroutine calc_column_tau
!***********************************************************************
    subroutine calc_cooling_time(f,p,taucool,taucool1)
!
      use EquationOfState, only: cs20,rho0
!
      real, dimension(mx,my,mz,mfarray) :: f
      type (pencil_case) :: p
      real, dimension(nx) :: tmp,tau_eff,Rd,OOK1
      real, dimension(nx), intent(out) :: taucool,taucool1
!
      integer :: i
!
      tmp=p%cp**1.5 * gamma1 * sqrt(gamma_m1) / (3.*sigmaSB) * p%rho * p%TT**(-2.5)
      tau_eff = 0.375*q%tau + .25*sqrt(3.) + .25/q%tau
      Rd=tmp*tau_eff
      OOK1 = (x(l1:l2)*sinth(m))**(1.5)  !inverse Keplerian frequency
      taucool = Rd*OOK1
!
      if (taucool_floor/=impossible) where (taucool < taucool_floor) taucool=taucool_floor
!
      taucool1 = 1./taucool
!
    endsubroutine calc_cooling_time
!***********************************************************************
    subroutine pencil_criteria_special
!
      lpenc_requested(i_TT)=.true.
      lpenc_requested(i_TT1)=.true.
      lpenc_requested(i_rho)=.true.
      lpenc_requested(i_cp)=.true.
      lpenc_requested(i_cv)=.true.
!
    endsubroutine pencil_criteria_special
!***********************************************************************
    subroutine calc_pencils_special(f,p)
!
!   14-jul-09/wlad: coded
!
      real, dimension(mx,my,mz,mfarray) :: f
      type (pencil_case) :: p
!
      q%tau=f(l1:l2,m,n,itau)
!
      call calc_cooling_time(f,p,q%taucool,q%taucool1)
!
    endsubroutine calc_pencils_special
!***********************************************************************
    subroutine calc_opacity(TT,rho,kk)
!
!  Piece-wise opacities from Bell et al. 1997.
!
!  01-aug-11/wlad: coded
!
      use General, only: notanumber
!
      real, intent(in) :: TT,rho
      real, intent(out) :: kk
      real :: k,a,b,logkk,logk
!
      real, save :: T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11
      logical, save :: lfirstcall=.true.
!
      if (lfirstcall) then
        T1=132. ; T2=170. ; T3=375. ; T4=390.
        T5=580. ; T6=680. ; T7=960. ; T8=1570. ; T9=3730.
        T10=1e4; T11=1e5
        lfirstcall=.false.
      endif
!
      if (TT < 0.0)        call inevitably_fatal_error("calc_opacity", "Negative temperature",.true.)
      if (rho < 0.0)       call inevitably_fatal_error("calc_opacity", "Negative density",.true.)
      if (notanumber(TT))  call inevitably_fatal_error("calc_opacity", "notanumber in temperature",.true.)
      if (notanumber(rho)) call inevitably_fatal_error("calc_opacity", "notanumber in density",.true.)
!
      if (TT <= T1) then
        k=2d-4 ; a=0 ; b= 2.1  ; kk = k * TT**b
      else if ((TT > T1) .and. (TT <= T2)) then
        k=3.   ; a=0 ; b=-0.01 ; kk = k * TT**b
      else if ((TT > T2) .and. (TT <= T3)) then
        k=0.01 ; a=0 ; b= 1.1  ; kk = k * TT**b
      else if ((TT > T3) .and. (TT <= T4)) then
        k=5d4  ; a=0 ; b=-1.5  ; kk = k * TT**b
      else if ((TT > T4) .and. (TT <= T5)) then
        k=0.1  ; a=0 ;  b= 0.7 ; kk = k * TT**b
      else if ((TT > T5) .and. (TT <= T6)) then
        k=2d15 ; a=0 ; b=-5.2  ; kk = k * TT**b
      else if ((TT > T6) .and. (TT <= T7)) then
        k=0.02 ; a=0 ; b= 0.8  ; kk = k * TT**b
      else if ((TT > T7) .and. (TT <= T8)) then
        logk=81.3010 ; a=1. ; b=-24.
        logkk = logk + a*alog10(rho) + b*alog10(TT)
        kk=10**(logkk)
      else if ((TT > T8) .and. (TT <= T9)) then
        k=1d-8 ; a=2./3 ; b=3.
        kk = k * rho**a * TT**b
      else if ((TT > T9) .and. (TT <= T10)) then
         logk=-36. ; a=1./3 ; b=10.
         logkk = logk + a*alog10(rho) + b*alog10(TT)
         kk=10**(logkk)
      else if ((TT > T10) .and. (TT <= T11)) then
         k=1.5d20 ; a=1. ; b=-2.5
         kk = k * rho**a * TT**b
      else if (TT > T11) then 
         kk=0.348
      endif
! 
    endsubroutine calc_opacity
!***********************************************************************
    subroutine read_special_init_pars(iostat)
!
      use File_io, only: parallel_unit
!
      integer, intent(out) :: iostat
!
      read(parallel_unit, NML=special_init_pars, IOSTAT=iostat)
!
    endsubroutine read_special_init_pars
!***********************************************************************
    subroutine write_special_init_pars(unit)
      integer, intent(in) :: unit

      write(unit,NML=special_init_pars)

    endsubroutine write_special_init_pars
!***********************************************************************
    subroutine read_special_run_pars(iostat)
!
      use File_io, only: parallel_unit
!
      integer, intent(out) :: iostat
!
      read(parallel_unit, NML=special_run_pars, IOSTAT=iostat)
!
    endsubroutine read_special_run_pars
!***********************************************************************
    subroutine write_special_run_pars(unit)
      integer, intent(in) :: unit
!
      write(unit,NML=special_run_pars)
!
    endsubroutine write_special_run_pars
!***********************************************************************
    subroutine rprint_special(lreset,lwrite)
!
      use Diagnostics
      use FArrayManager, only: farray_index_append
!
!  reads and registers print parameters relevant to special
!
!   14-jul-09/wlad: coded
!
      integer :: iname
      logical :: lreset,lwr
      logical, optional :: lwrite
!
      lwr = .false.
      if (present(lwrite)) lwr=lwrite
!
!  Write information to index.pro
!
      if (lreset) then
        idiag_kappam=0;idiag_kappamax=0;idiag_kappamin=0 
        idiag_taum=0;idiag_taumax=0;idiag_taumin=0 
        idiag_taucoolm=0;idiag_taucoolmax=0;idiag_taucoolmin=0 
        idiag_dts=0
      endif
!
      do iname=1,nname
        call parse_name(iname,cname(iname),cform(iname),'dts',idiag_dts)
!
        call parse_name(iname,cname(iname),cform(iname),'kappam',idiag_kappam)
        call parse_name(iname,cname(iname),cform(iname),'kappamax',idiag_kappamax)
        call parse_name(iname,cname(iname),cform(iname),'kappamin',idiag_kappamin)
!
        call parse_name(iname,cname(iname),cform(iname),'taum',idiag_taum)
        call parse_name(iname,cname(iname),cform(iname),'taumax',idiag_taumax)
        call parse_name(iname,cname(iname),cform(iname),'taumin',idiag_taumin)
!
        call parse_name(iname,cname(iname),cform(iname),'taucoolm',idiag_taucoolm)
        call parse_name(iname,cname(iname),cform(iname),'taucoolmax',idiag_taucoolmax)
        call parse_name(iname,cname(iname),cform(iname),'taucoolmin',idiag_taucoolmin)

      enddo
!
      if (lwr) then
        call farray_index_append('i_dts',idiag_dts)
!
        call farray_index_append('i_kappam',idiag_kappam)
        call farray_index_append('i_kappamax',idiag_kappamax)
        call farray_index_append('i_kappamin',idiag_kappamin)
!
        call farray_index_append('i_taum',idiag_taum)
        call farray_index_append('i_taumax',idiag_taumax)
        call farray_index_append('i_taumin',idiag_taumin)
!
        call farray_index_append('i_taucoolm',idiag_taucoolm)
        call farray_index_append('i_taucoolmax',idiag_taucoolmax)
        call farray_index_append('i_taucoolmin',idiag_taucoolmin)
      endif
!
    endsubroutine rprint_special
!***********************************************************************
    subroutine special_calc_energy(f,df,p)
!
!   calculate a additional 'special' term on the right hand side of the 
!   momentum equation.
!
!   Some precalculated pencils of data are passed in for efficiency
!   others may be calculated directly from the f array
!
      use EquationOfState, only: cs20
!      
      real, dimension (mx,my,mz,mvar+maux), intent(in) :: f
      real, dimension (mx,my,mz,mvar), intent(inout) :: df
      type (pencil_case), intent(in) :: p
      real, dimension(nx) :: rr_cyl,TT_init,heating_rate
!
!  Modified entropy equation.
!
      if (laddheatingrate) then
        rr_cyl=x(l1:l2)*sinth(m)
        TT_init = cs20/(p%cp*gamma_m1) * r_ref / rr_cyl
        heating_rate = p%cv * (p%TT-TT_init)*q%taucool1
!
        df(l1:l2,m,n,iss) = df(l1:l2,m,n,iss) - p%TT1 * heating_rate
        if (lfirst.and.ldt) dt1_max=max(dt1_max,q%taucool1/cdts)
      endif
!
      call special_calc_diagnostics(f)
!
    endsubroutine special_calc_energy
!***********************************************************************
    subroutine special_calc_diagnostics(f)

      use Diagnostics
!
      real, dimension (mx,my,mz,mvar+maux), intent(in) :: f
!
      real, dimension (nx) :: kappa
!
      if (ldiagnos) then 
      
        kappa=f(l1:l2,m,n,ikappar)
!
        call sum_mn_name(kappa,idiag_kappam)
        call max_mn_name(kappa,idiag_kappamax)
        if (idiag_kappamin/=0) call max_mn_name(-kappa,idiag_kappamin,lneg=.true.)
!
        call sum_mn_name(q%tau,idiag_taum)
        call max_mn_name(q%tau,idiag_taumax)
        if (idiag_taumin/=0) call max_mn_name(-q%tau,idiag_taumin,lneg=.true.)
!
        call sum_mn_name(q%taucool,idiag_taucoolm)
        call max_mn_name(q%taucool,idiag_taucoolmax)
        if (idiag_taucoolmin/=0) call max_mn_name(-q%taucool,idiag_taucoolmin,lneg=.true.)
!
        if (idiag_dts/=0) call max_mn_name(q%taucool1/cdts,idiag_dts,l_dt=.true.)
      endif
!
    endsubroutine special_calc_diagnostics
!***********************************************************************
!********************************************************************
!************        DO NOT DELETE THE FOLLOWING       **************
!********************************************************************
!**  This is an automatically generated include file that creates  **
!**  copies dummy routines from nospecial.f90 for any Special      **
!**  routines not implemented in this file                         **
!**                                                                **
    include '../special_dummies.inc'
!********************************************************************
endmodule Special

