! $Id$
!
! This module serves as a sample for a special_XXX module that
! introduces additional primitive variables. Use this as a basis for your
! own special_ module if you need one.
!
! To ensure it is kept up to date, we also use it for production stuff.
! This sample modules solves the dynamical alpha quenching equation,
! involving mean-field theory, which explains the presence of a number
! of non-generic routines
!
!
!** 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 2
! MAUX CONTRIBUTION 0
!
!***************************************************************
!
module Special
!
use Cparam
use Cdata
use General, only: keep_compiler_quiet
use Messages
!
implicit none
!
include '../special.h'
!
character (len=labellen) :: initalpm='zero',Omega_profile='nothing', initetam='constant'
!
! input parameters
real :: amplalpm=.1, ampletat=1.
real :: kx_alpm=1.,ky_alpm=1.,kz_alpm=1.
real :: Omega_ampl=.0
!
namelist /special_init_pars/ &
initalpm,initetam,amplalpm,ampletat,kx_alpm,ky_alpm,kz_alpm
!
! run parameters
real :: kf_alpm=1., alpmdiff=0.
logical :: ladvect_alpm=.false.
!
namelist /special_run_pars/ &
initetam,kf_alpm,ladvect_alpm,alpmdiff, &
Omega_profile,Omega_ampl
!
! Declare any index variables necessary
!
! other variables (needs to be consistent with reset list below)
integer :: idiag_etatm=0,idiag_etmax=0,idiag_etrms=0
integer :: idiag_alpmm=0,idiag_ammax=0,idiag_amrms=0,idiag_alpmmz=0
!
logical, pointer :: lmeanfield_theory
real, pointer :: meanfield_etat,eta
!
contains
!
!***********************************************************************
subroutine register_special()
!
! Initialise variables which should know that we solve for passive
! scalar: ialpm; increase nvar accordingly
!
! 6-jul-02/axel: coded
!
use FArrayManager
!
! register ialpm in the f-array and set lalpm=.false.
!
call farray_register_pde('alpm',ialpm)
call farray_register_pde('etat',ietat)
lalpm=.true.
! letat=.true.
!
! Identify version number.
!
if (lroot) call svn_id( &
"$Id$")
!
! Writing files for use with IDL
!
if (lroot) then
if (maux == 0) then
if (nvar < mvar) write(4,*) ',alpm $'
if (nvar == mvar) write(4,*) ',alpm'
else
write(4,*) ',alpm $'
endif
write(15,*) 'alpm = fltarr(mx,my,mz,2)*one'
endif
!
endsubroutine register_special
!***********************************************************************
subroutine initialize_special(f)
!
! Perform any necessary post-parameter read initialization
! Dummy routine
!
! 24-nov-02/tony: coded
!
real, dimension (mx,my,mz,mvar+maux) :: f
!
! set to zero and then call the same initial condition
! that was used in start.csh
!
! set the magnetic Reynold number :
! Rm_alpm=etat_alpm/eta
! write(*,*) 'Dhruba', Rm_alpm,etat_alpm
!
call keep_compiler_quiet(f)
!
endsubroutine initialize_special
!***********************************************************************
subroutine init_special(f)
!
! initialise passive scalar field; called from start.f90
!
! 6-jul-2001/axel: coded
!
use Mpicomm
! use Density
use Sub
use Initcond
!
real, dimension (mx,my,mz,mvar+maux) :: f
!
! inititialize alpm and etat. Not that the f-array in ietat only contains
! the part not already included in meanfield_etat
!
select case (initalpm)
case ('zero'); f(:,:,:,ialpm)=0.; f(:,:,:,ietat)=0.
case ('constant'); f(:,:,:,ialpm)=amplalpm; f(:,:,:,ietat)=ampletat
case default; call stop_it('init_alpm: bad initalpm='//trim(initalpm))
endselect
!
endsubroutine init_special
!***********************************************************************
subroutine calc_pencils_special(f,p)
!
! Calculate Hydro pencils.
! Most basic pencils should come first, as others may depend on them.
!
! 24-nov-04/tony: coded
!
real, dimension (mx,my,mz,mvar+maux) :: f
type (pencil_case) :: p
!
intent(in) :: f
intent(inout) :: p
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(p)
!
endsubroutine calc_pencils_special
!***********************************************************************
subroutine dspecial_dt(f,df,p)
!
! dynamical alpha quenching equation
! dalpm/dt=-2*etat*kf2*(EMF*BB/Beq2+alpm/Rm)
! detat/dt=-(1/3)*(EMF*JJ-kf*EMF*BB)/(eta*kf2+kf*sqrt(EMF*JJ-kf*EMF*BB))
!
! 18-nov-04/axel: coded
! 17-sep-09/axel+koen: added etat evolution
!
use Sub
use Diagnostics
use SharedVariables, only : get_shared_variable
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (mx,my,mz,mvar) :: df
real, dimension (nx,3) :: galpm
real, dimension (nx) :: alpm,etat,ugalpm,EMFdotB,EMFdotJ
real, dimension (nx) :: divflux,del2alpm,EJ_kfEB
type (pencil_case) :: p
integer :: ierr
!
intent(in) :: f
intent(out) :: df
!
! identify module and boundary conditions
!
if (headtt) call identify_bcs('alpm',ialpm)
if (headtt) call identify_bcs('etat',ietat)
!
! get meanfield_etat and eta. Leave df(l1:l2,m,n,ialpm) unchanged
! if lmeanfield_theory is false.
!
call get_shared_variable('lmeanfield_theory',lmeanfield_theory,ierr)
if (ierr/=0) &
call fatal_error("dspecial_dt: ", &
"cannot get shared var lmeanfield_theory")
if (lmeanfield_theory) then
call get_shared_variable('meanfield_etat',meanfield_etat,ierr)
if (ierr/=0) &
call fatal_error("dspecial_dt: ", &
"cannot get shared var meanfield_etat")
call get_shared_variable('eta',eta,ierr)
if (ierr/=0) &
call fatal_error("dspecial_dt: ", "cannot get shared var eta")
!
! Abbreviations
!
alpm=f(l1:l2,m,n,ialpm)
etat=f(l1:l2,m,n,ietat)+meanfield_etat
!
! dynamical quenching equation
! with advection flux proportional to uu
!
if (lmagnetic) then
call dot_mn(p%mf_EMF,p%bb,EMFdotB)
call divflux_from_Omega_effect(p,divflux)
df(l1:l2,m,n,ialpm)=df(l1:l2,m,n,ialpm)&
-2*kf_alpm**2*(etat*EMFdotB+eta*alpm)-etat*divflux
if (ladvect_alpm) then
call grad(f,ialpm,galpm)
call dot_mn(p%uu,galpm,ugalpm)
df(l1:l2,m,n,ialpm)=df(l1:l2,m,n,ialpm)-ugalpm
endif
if (alpmdiff/=0) then
call del2(f,ialpm,del2alpm)
df(l1:l2,m,n,ialpm)=df(l1:l2,m,n,ialpm)+alpmdiff*del2alpm
endif
!
! etat evolution
!
select case (initetam)
case ('evolving');
!
! compute d_t eta= tau/3 *d_t <u^2> with 1/tau=kf^2(eta+etat)
! and d_t <u^2>= -2*J.E_EMF+2kfB.E_EMF
!
call dot_mn(p%mf_EMF,p%jj,EMFdotJ)
EJ_kfEB=EMFdotJ-kf_alpm*EMFdotB
df(l1:l2,m,n,ietat)=df(l1:l2,m,n,ietat)&
-(2./3.)*EJ_kfEB/kf_alpm**2/(eta+etat)
case ('constant'); df(:,:,:,ietat)=0.
endselect
endif
endif
!
! diagnostics
!
if (ldiagnos) then
if (idiag_etatm/=0) call sum_mn_name(etat,idiag_etatm)
if (idiag_etmax/=0) call max_mn_name(alpm,idiag_etmax)
if (idiag_etrms/=0) call sum_mn_name(alpm**2,idiag_etrms,lsqrt=.true.)
if (idiag_alpmm/=0) call sum_mn_name(alpm,idiag_alpmm)
if (idiag_ammax/=0) call max_mn_name(alpm,idiag_ammax)
if (idiag_amrms/=0) call sum_mn_name(alpm**2,idiag_amrms,lsqrt=.true.)
endif
!
endsubroutine dspecial_dt
!***********************************************************************
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)
!
! reads and registers print parameters relevant for magnetic fields
!
! 6-jul-02/axel: coded
!
use Diagnostics
use FArrayManager, only: farray_index_append
!
integer :: iname,inamez
logical :: lreset,lwr
logical, optional :: lwrite
!
lwr = .false.
if (present(lwrite)) lwr=lwrite
!
! reset everything in case of reset
! (this needs to be consistent with what is defined above!)
!
if (lreset) then
idiag_etatm=0; idiag_etmax=0; idiag_etrms=0
idiag_alpmm=0; idiag_ammax=0; idiag_amrms=0; idiag_alpmmz=0
endif
!
! check for those quantities that we want to evaluate online
!
do iname=1,nname
call parse_name(iname,cname(iname),cform(iname),'etatm',idiag_etatm)
call parse_name(iname,cname(iname),cform(iname),'etmax',idiag_etmax)
call parse_name(iname,cname(iname),cform(iname),'etrms',idiag_etrms)
call parse_name(iname,cname(iname),cform(iname),'alpmm',idiag_alpmm)
call parse_name(iname,cname(iname),cform(iname),'ammax',idiag_ammax)
call parse_name(iname,cname(iname),cform(iname),'amrms',idiag_amrms)
enddo
!
! check for those quantities for which we want xy-averages
!
do inamez=1,nnamez
call parse_name(inamez,cnamez(inamez),cformz(inamez),&
'alpmmz',idiag_alpmmz)
enddo
!
! write column where which magnetic variable is stored
!
if (lwr) then
call farray_index_append('i_etatm',idiag_etatm)
call farray_index_append('i_etmax',idiag_etmax)
call farray_index_append('i_etrms',idiag_etrms)
call farray_index_append('i_alpmm',idiag_alpmm)
call farray_index_append('i_ammax',idiag_ammax)
call farray_index_append('i_amrms',idiag_amrms)
call farray_index_append('ispecial',ialpm)
endif
!
endsubroutine rprint_special
!***********************************************************************
subroutine divflux_from_Omega_effect(p,divflux)
!
! Omega effect coded (normally used in context of mean field theory)
! Can do uniform shear (0,Sx,0), and the cosx*cosz profile (solar CZ).
!
! 30-apr-05/axel: coded
!
real, dimension (nx) :: divflux
type (pencil_case) :: p
!
intent(in) :: p
intent(out) :: divflux
!
! Fi = a*eps_ijl Slk BjBk
!
select case (Omega_profile)
case ('nothing'); if (headtt) print*,'Omega_profile=nothing'
divflux=0 ! or we will be using uninitialized memory...
case ('(0,Sx,0)')
if (headtt) print*,'divflux: uniform shear, S=',Omega_ampl
divflux=Omega_ampl*(p%bb(:,1)*p%bij(:,1,3)-p%bb(:,2)*p%bij(:,2,3))
case ('(0,cosx*cosz,0)')
if (headtt) print*,'divflux: solar shear, S=',Omega_ampl
divflux=Omega_ampl*((p%bb(:,2)*p%bij(:,2,1)- p%bb(:,3)*p%bij(:,3,1)&
+.5*p%bb(:,3)*p%bij(:,1,3)+.5*p%bb(:,1)*p%bij(:,3,3))&
*cos(x(l1:l2))*sin(z(n))&
-(p%bb(:,2)*p%bij(:,2,3)- p%bb(:,1)*p%bij(:,1,3)&
+.5*p%bb(:,3)*p%bij(:,1,1)+.5*p%bb(:,1)*p%bij(:,3,1))&
*sin(x(l1:l2))*cos(z(n))&
+(p%bb(:,1)**2-p%bb(:,3)**2)&
*sin(x(l1:l2))*sin(z(n)))
case default; print*,'Omega_profile=unknown'
divflux=0 ! or we will be using uninitialized memory...
endselect
!
endsubroutine divflux_from_Omega_effect
!***********************************************************************
!
!********************************************************************
!************ 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