! $Id$ ! !** 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 0 ! !*************************************************************** ! module Special ! use Cparam use Cdata use General, only: keep_compiler_quiet use Messages, only: svn_id, fatal_error ! implicit none ! include '../special.h' ! real :: Kpara=0.,Kperp=0.,Ksat=0.,Kc=0.,cool_RTV=0.,Kchrom=0. real :: exp_RTV=0.,cubic_RTV=0.,tanh_RTV=0. real :: tau_inv_newton=0.,exp_newton=0. real :: tanh_newton=0.,cubic_newton=0. real :: width_newton=0.,width_RTV=0.,hyper3_chi=0. real :: init_time=0.,lnTT0_chrom=0.,width_lnTT_chrom=0. real :: hcond_grad_iso=0.,chi_aniso=0. real :: tau_inv_spitzer=0. real :: bullets_x0=0.,bullets_dx=0. real :: bullets_t0=0.,bullets_dt=0. real :: bullets_h0=0.,hyper3_diffrho=0. logical :: lfilter_farray=.false.,lrad_loss=.false. real, dimension(mvar) :: filter_strength=0.01 real :: Ltot,R_hyper3=0. real, dimension (nx) :: diffus_chi ! character (len=labellen), dimension(3) :: iheattype='nothing' character (len=labellen) :: loop_frac='full' real, dimension(2) :: heat_par_exp=(/0.,1./) real, dimension(2) :: heat_par_exp2=(/0.,1./) real, dimension(3) :: heat_par_gauss=(/0.,1.,0./) real, dimension(3) :: heat_par_vandoors=(/1.,0.1,50./) ! namelist /special_run_pars/ & Kpara,Kperp,cool_RTV,tau_inv_newton,exp_newton,init_time, & iheattype,heat_par_exp,heat_par_exp2,heat_par_gauss, & width_newton,tanh_newton,cubic_newton,Kchrom, & lnTT0_chrom,width_lnTT_chrom,width_RTV, & exp_RTV,cubic_RTV,tanh_RTV,hcond_grad_iso,Ksat,Kc, & tau_inv_spitzer,hyper3_chi,bullets_x0,bullets_dx, & bullets_t0,bullets_dt,bullets_h0,hyper3_diffrho, & lfilter_farray,filter_strength,loop_frac,& heat_par_vandoors,R_hyper3,lrad_loss,chi_aniso ! ! variables for print.in ! integer :: idiag_dtchi2=0 ! DIAG_DOC: heatconduction ! DIAG DOC: in special module integer :: idiag_dthyper3=0 integer :: idiag_dtrad=0 ! DIAG_DOC: radiative loss from RTV integer :: idiag_dtspitzer=0 ! DIAG_DOC: Spitzer heat conduction ! DIAG_DOC: time step integer :: idiag_dtnewt=0 integer :: idiag_qmax=0 ! DIAG_DOC: max of heat flux vector integer :: idiag_qrms=0 ! DIAG_DOC: rms of heat flux vector ! ! variables for video.in ! integer :: ivid_logQ=0; ivid_rtv=0; ivid_spitzer=0 real, target, dimension(:,:), allocatable :: rtv_xy, rtv_xy2, rtv_xy3, rtv_xy4 real, target, dimension(:,:), allocatable :: rtv_xz, rtv_yz, rtv_xz2 real, target, dimension(:,:,:,:,:), allocatable :: rtv_r real, target, dimension(:,:), allocatable :: logQ_xy, logQ_xy2, logQ_xy3, logQ_xy4 real, target, dimension(:,:), allocatable :: logQ_xz, logQ_xz2, logQ_yz real, target, dimension(:,:,:,:,:), allocatable :: logQ_r real, target, dimension(:,:), allocatable :: spitzer_xy, spitzer_xy2, spitzer_xy3, spitzer_xy4 real, target, dimension(:,:), allocatable :: spitzer_xz, spitzer_xz2, spitzer_yz real, target, dimension(:,:,:,:,:), allocatable :: spitzer_r ! ! miscellaneous variables ! real, save, dimension (mx) :: lnTT_init_prof,lnrho_init_prof integer, save, dimension(mseed) :: nano_seed ! real :: Kspitzer_para_SI = 2e-11, Kspitzer_para=0. real :: Ksaturation_SI = 7e7,Ksaturation=0.,ln_unit_TT=0. ! contains !*********************************************************************** subroutine register_special() ! if (lroot) call svn_id( & "$Id$") ! endsubroutine register_special !*********************************************************************** subroutine initialize_special(f) ! ! Called after reading parameters, but before the time loop. ! ! 13-sep-10/bing: coded ! use Slices_methods, only: alloc_slice_buffers real, dimension (mx,my,mz,mfarray), intent(inout) :: f real, dimension (mx) :: xtmp character (len=*), parameter :: filename='/strat.dat' integer :: lend,unit=12 real :: dummy=1. ! ln_unit_TT = alog(real(unit_temperature)) ! Kspitzer_para = Kspitzer_para_SI /unit_density/unit_velocity**3./ & unit_length*unit_temperature**(3.5) ! Ksaturation = Ksaturation_SI /unit_velocity**3.*unit_temperature**1.5 ! call keep_compiler_quiet(f) ! inquire(IOLENGTH=lend) dummy ! if (lrun .and. (tau_inv_newton/=0.)) then open(unit,file=trim(directory_snap)//filename, & form='unformatted',status='unknown',recl=lend*mx) read(unit) xtmp read(unit) lnrho_init_prof read(unit) lnTT_init_prof close(unit) endif ! select case (loop_frac) case ('full') Ltot = Lxyz(1) + 2*xyz0(1) case ('half') Ltot = 2*Lxyz(1) + 2*xyz0(1) case default call fatal_error('initialize_initial_condition','Wrong loop_frac') Ltot = 0. endselect ! if (ivid_rtv/=0) & call alloc_slice_buffers(rtv_xy,rtv_xz,rtv_yz,rtv_xy2,rtv_xy3,rtv_xy4,rtv_xz2,rtv_r) if (ivid_logQ/=0) & call alloc_slice_buffers(logQ_xy,logQ_xz,logQ_yz,logQ_xy2,logQ_xy3,logQ_xy4,logQ_xz2,logQ_r) if (ivid_spitzer/=0) & call alloc_slice_buffers(spitzer_xy,spitzer_xz,spitzer_yz, & spitzer_xy2,spitzer_xy3,spitzer_xy4,spitzer_xz2,spitzer_r) ! endsubroutine initialize_special !*********************************************************************** 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 pencil_criteria_special() ! ! All pencils that this special module depends on are specified here. ! ! 04-sep-10/bing: coded ! if (tau_inv_spitzer/=0.) then lpenc_requested(i_uglnrho)=.true. lpenc_requested(i_divu)=.true. lpenc_requested(i_cp1)=.true. lpenc_requested(i_lnrho)=.true. lpenc_requested(i_glnrho)=.true. lpenc_requested(i_lnTT)=.true. lpenc_requested(i_glnTT)=.true. endif ! if (chi_aniso/=0.) then lpenc_requested(i_bij)=.true. lpenc_requested(i_bunit)=.true. lpenc_requested(i_glnTT)=.true. lpenc_requested(i_hlnTT)=.true. lpenc_requested(i_glnrho)=.true. endif ! if (Kpara/=0.) then lpenc_requested(i_cp1)=.true. lpenc_requested(i_lnrho)=.true. lpenc_requested(i_lnTT)=.true. lpenc_requested(i_glnTT)=.true. lpenc_requested(i_del2lnTT)=.true. endif ! if (hcond_grad_iso/=0.) then lpenc_requested(i_glnTT)=.true. lpenc_requested(i_hlnTT)=.true. lpenc_requested(i_del2lnTT)=.true. lpenc_requested(i_glnrho)=.true. lpenc_requested(i_cp1)=.true. endif ! if (Kchrom/=0.) then lpenc_requested(i_cp1)=.true. lpenc_requested(i_rho1)=.true. lpenc_requested(i_lnTT)=.true. lpenc_requested(i_glnTT)=.true. lpenc_requested(i_hlnTT)=.true. lpenc_requested(i_del2lnTT)=.true. endif ! if (cool_RTV/=0.) then lpenc_requested(i_cp1)=.true. lpenc_requested(i_rho)=.true. lpenc_requested(i_lnTT)=.true. lpenc_requested(i_lnrho)=.true. endif ! if ((tau_inv_newton/=0.).or.(bullets_h0/=0.)) then lpenc_requested(i_lnTT)=.true. lpenc_requested(i_lnrho)=.true. endif ! if (iheattype(1)/='nothing') then lpenc_requested(i_TT1)=.true. lpenc_requested(i_rho1)=.true. lpenc_requested(i_cp1)=.true. lpenc_requested(i_glnrho)=.true. endif ! if (Ksat /= 0.) then lpenc_requested(i_TT)=.true. lpenc_requested(i_glnTT)=.true. lpenc_requested(i_hlnTT)=.true. lpenc_requested(i_glnrho)=.true. lpenc_requested(i_cp1)=.true. endif ! if (R_hyper3 /= 0.) lpenc_requested(i_u2)=.true. ! endsubroutine pencil_criteria_special !*********************************************************************** subroutine pencil_interdep_special(lpencil_in) ! logical, dimension (npencils) :: lpencil_in ! call keep_compiler_quiet(lpencil_in) ! endsubroutine pencil_interdep_special !*********************************************************************** subroutine calc_pencils_special(f,p) ! real, dimension (mx,my,mz,mfarray) :: f type (pencil_case) :: p ! call keep_compiler_quiet(f) call keep_compiler_quiet(p) ! endsubroutine calc_pencils_special !*********************************************************************** subroutine rprint_special(lreset,lwrite) ! ! reads and registers print parameters relevant to special ! ! 04-sep-10/bing: coded ! use Diagnostics, only: parse_name use FArrayManager, only: farray_index_append ! integer :: iname 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_dtchi2=0. idiag_dthyper3=0. idiag_dtrad=0. idiag_dtnewt=0 idiag_dtspitzer=0 idiag_qmax=0 idiag_qrms=0 ivid_logQ=0; ivid_rtv=0; ivid_spitzer=0 endif ! ! iname runs through all possible names that may be listed in print.in ! do iname=1,nname call parse_name(iname,cname(iname),cform(iname),'dtchi2',idiag_dtchi2) call parse_name(iname,cname(iname),cform(iname),'dthyper3',idiag_dthyper3) call parse_name(iname,cname(iname),cform(iname),'dtrad',idiag_dtrad) call parse_name(iname,cname(iname),cform(iname),'dtnewt',idiag_dtnewt) call parse_name(iname,cname(iname),cform(iname),'dtspitzer',idiag_dtspitzer) call parse_name(iname,cname(iname),cform(iname),'qmax',idiag_qmax) call parse_name(iname,cname(iname),cform(iname),'qrms',idiag_qrms) enddo ! ! check for those quantities for which we want video slices ! do iname = 1,nnamev call parse_name(iname,cnamev(iname),cformv(iname),'logQ',ivid_logQ) call parse_name(iname,cnamev(iname),cformv(iname),'rtv',ivid_rtv) call parse_name(iname,cnamev(iname),cformv(iname),'spitzer',ivid_spitzer) enddo ! ! write column where which variable is stored ! if (lwr) then call farray_index_append('i_dtchi2',idiag_dtchi2) call farray_index_append('i_dthyper3',idiag_dthyper3) call farray_index_append('i_dtrad',idiag_dtrad) call farray_index_append('i_dtnewt',idiag_dtnewt) call farray_index_append('i_dtspitzer',idiag_dtspitzer) call farray_index_append('i_qmax',idiag_qmax) call farray_index_append('i_qrms',idiag_qrms) endif ! endsubroutine rprint_special !*********************************************************************** subroutine dspecial_dt(f,df,p) ! ! 12-may-11/bingert: coded ! use EquationOfState, only: gamma use Deriv, only: der6 use Diagnostics, only : max_mn_name, sum_mn_name use Sub, only: identify_bcs, multsv, dot, del6 ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz,mvar) :: df type (pencil_case) :: p real, dimension (nx) :: rhs,tmp,hyper3_coeff,hc real, dimension (nx,3) :: K1 integer :: i,j ! intent(in) :: p intent(inout) :: df ! ! Identify module and boundary conditions. ! if (headtt.or.ldebug) print*,'dspecial_dt: SOLVE dspecial_dt' ! if (lfilter_farray) call filter_farray(f,df) ! ! Add hyper diffusion with fixed Reynoldsnmumber ! if (R_hyper3 /= 0.) then hyper3_coeff = sqrt(p%u2)/dxmax_pencil/R_hyper3 call del6(f, ilnTT, hc, IGNOREDX=.true.) df(l1:l2,m,n,ilnTT) = df(l1:l2,m,n,ilnTT) + hyper3_coeff * hc call del6(f, ilnrho, hc, IGNOREDX=.true.) df(l1:l2,m,n,ilnrho) = df(l1:l2,m,n,ilnrho) + hyper3_coeff * hc call del6(f, iux, hc, IGNOREDX=.true.) df(l1:l2,m,n,iux) = df(l1:l2,m,n,iux) + hyper3_coeff * hc call del6(f, iuy, hc, IGNOREDX=.true.) df(l1:l2,m,n,iuy) = df(l1:l2,m,n,iuy) + hyper3_coeff * hc call del6(f, iuz, hc, IGNOREDX=.true.) df(l1:l2,m,n,iuz) = df(l1:l2,m,n,iuz) + hyper3_coeff * hc ! if (lfirst.and.ldt) dt1_max=max(dt1_max,tmp/0.01) endif ! endsubroutine dspecial_dt !*********************************************************************** subroutine get_slices_special(f,slices) ! ! Write slices for animation of special variables. ! ! 26-jun-06/tony: dummy ! use Slices_methods, only: assign_slices_scal ! real, dimension (mx,my,mz,mfarray) :: f type (slice_data) :: slices ! ! Loop over slices ! select case (trim(slices%name)) ! case ('rtv'); call assign_slices_scal(slices,rtv_xy,rtv_xz,rtv_yz,rtv_xy2, & rtv_xy3,rtv_xy4,rtv_xz2,rtv_r) ! case ('logQ'); call assign_slices_scal(slices,logQ_xy,logQ_xz,logQ_yz,logQ_xy2, & logQ_xy3,logQ_xy4,logQ_xz2,logQ_r) ! case ('spitzer') call assign_slices_scal(slices,spitzer_xy,spitzer_xz,spitzer_yz,spitzer_xy2, & spitzer_xy3,spitzer_xy4,spitzer_xz2,spitzer_r) endselect ! call keep_compiler_quiet(f) ! endsubroutine get_slices_special !*********************************************************************** subroutine special_after_timestep(f,df,dt_,llast) ! ! 10-oct-12/bing: coded ! use EquationOfState, only: gamma ! logical, intent(in) :: llast real, dimension(mx,my,mz,mfarray), intent(inout) :: f real, dimension(mx,my,mz,mvar), intent(inout) :: df real, intent(in) :: dt_ integer,save :: j=35 real :: lnTT_SI,dt_step,dt_rad,ln_coeff,one_m_alpha,lnTT_res integer :: l logical :: notdone ! real, parameter, dimension (37) :: intlnT = (/ & 8.74982, 8.86495, 8.98008, 9.09521, 9.21034, 9.44060, 9.67086 & , 9.90112, 10.1314, 10.2465, 10.3616, 10.5919, 10.8221, 11.0524 & , 11.2827, 11.5129, 11.7432, 11.9734, 12.2037, 12.4340, 12.6642 & , 12.8945, 13.1247, 13.3550, 13.5853, 13.8155, 14.0458, 14.2760 & , 14.5063, 14.6214, 14.7365, 14.8517, 14.9668, 15.1971, 15.4273 & , 15.6576, 69.0776 /) real, parameter, dimension (36) :: alpha_rad = (/ & 67.4287, 40.0384, 21.3332, 20.0001, 9.33333, 4.66647 & , 2.33345,-0.566705,-0.199849, 0.199849, 0.899692, 1.33405 & , 1.66611, 0.833237,-0.166805, -1.49977, 0.00000, 0.566656 & , 0.233155,-0.966534,-0.633508, 0.800159, 0.433348, -0.0998807 & ,-0.499988, -1.00000, -1.96697, -3.06643, -3.60040, -5.80186 & , -1.86549, -2.66724,-0.166734, 0.566900, 0.666504, 6.23231 /) real, parameter, dimension (36) :: lnchi_rad = (/ & -690.935, -448.121, -280.147, -268.022, -169.777, -125.719 & , -103.157, -74.4422, -78.1590, -82.2545, -89.5060, -94.1066 & , -97.7002, -88.4950, -77.2118, -61.8654, -79.4776, -86.2624 & , -82.1925, -67.2755, -71.4930, -89.9794, -85.1652, -78.0439 & , -72.6083, -65.7003, -52.1185, -36.4226, -28.6767, 3.51167 & , -54.4966, -42.5892, -80.0138, -91.1629, -92.6996, -179.847 /) ! if (llast .and. lrad_loss) then ! do l=l1,l2 do m=m1,m2 do n=n1,n2 notdone=.true. dt_rad = dt*unit_time lnTT_SI = f(l,m,n,ilnTT) + ln_unit_TT do while (notdone) ! if (lnTT_SI > intlnT(j) .and. & lnTT_SI <= intlnT(j+1) ) then ! we are the propper intervall, compute the factors ! one_m_alpha = (1-alpha_rad(j)) ! ln_coeff = f(l,m,n,ilnrho) + alog(real(1.5**2./8.7*(gamma-1.)/0.6/m_p/k_B))+ lnchi_rad(j) ln_coeff = ln_coeff + alog(real(unit_temperature/unit_velocity**2./unit_density/unit_length**6.)) ! ! compute resulting temperature in SI units ! ! first check if we would get negative temperatues if (-one_m_alpha*dt_rad*exp(ln_coeff)+exp(one_m_alpha*lnTT_SI) < 0) then lnTT_res = intlnT(j)*0.9 else lnTT_res =1./one_m_alpha*alog( & -one_m_alpha*dt_rad*exp(ln_coeff) & +exp(one_m_alpha*lnTT_SI)) endif ! if (lnTT_res >= intlnT(j)) then ! everything is fine f(l,m,n,ilnTT) = lnTT_res - ln_unit_TT notdone = .false. else ! timestep to large, we need to repeat for another intervall dt_step = exp(-ln_coeff)/one_m_alpha* & (exp(one_m_alpha*lnTT_SI) - exp(one_m_alpha*intlnT(j))) dt_rad = dt_rad - dt_step f(l,m,n,ilnTT) = intlnT(j)-ln_unit_TT lnTT_SI = f(l,m,n,ilnTT) + ln_unit_TT j = j-1 if (j<=0) notdone=.false. endif elseif (lnTT_SI == intlnT(1)) then notdone=.false. else j = j + sign(1.,lnTT_SI-intlnT(j)) if (j <= 0) then j=1 notdone=.false. elseif (j >= 37) then j=36 notdone=.false. endif endif enddo ! while loop enddo enddo enddo endif ! if (llast) ! call keep_compiler_quiet(df) call keep_compiler_quiet(dt_) call keep_compiler_quiet(llast) ! endsubroutine special_after_timestep !*********************************************************************** subroutine special_calc_energy(f,df,p) ! ! Calculate an additional 'special' term on the right hand side of the ! entropy equation. ! ! Some precalculated pencils of data are passed in for efficiency ! others may be calculated directly from the f array ! ! 06-oct-03/tony: coded ! use Sub, only: del6 ! real, dimension (mx,my,mz,mfarray), intent(in) :: f real, dimension (mx,my,mz,mvar), intent(inout) :: df type (pencil_case), intent(in) :: p ! real, dimension (nx) :: hc integer :: itmp ! diffus_chi=0. ! if (hcond_grad_iso/=0.) call calc_heatcond_glnTT_iso(df,p) if (Kpara/=0.) call calc_heatcond_spitzer(f,df,p) if (cool_RTV/=0.) call calc_heat_cool_RTV(df,p) if (tau_inv_newton/=0.) call calc_heat_cool_newton(df,p) if (iheattype(1)/='nothing') call calc_artif_heating(df,p) if (Kchrom/=0.) call calc_heatcond_kchrom(df,p) if (bullets_h0/=0.) call calc_bullets_heating(df,p) if (chi_aniso/=0.) call calc_heatcond_constchi(df,p) ! if (ltemperature) then itmp = ilnTT else if (lentropy) then itmp = iss else call fatal_error('hyper3_chi special','only for ltemperature') endif ! if (hyper3_chi /= 0.) then call del6(f,itmp,hc,IGNOREDX=.true.) df(l1:l2,m,n,itmp) = df(l1:l2,m,n,itmp) + hyper3_chi*hc ! ! due to ignoredx hyper3_chi has [1/s] ! if (lfirst.and.ldt) dt1_max=max(dt1_max,hyper3_chi/0.01) endif ! if (lfirst.and.ldt) then maxdiffus=max(maxdiffus,diffus_chi) endif ! endsubroutine special_calc_energy !*********************************************************************** subroutine special_calc_density(f,df,p) ! use Sub, only: del6 ! real, dimension (mx,my,mz,mfarray), intent(in) :: f real, dimension (mx,my,mz,mvar), intent(inout) :: df type (pencil_case), intent(in) :: p ! real, dimension (nx) :: hc ! if (hyper3_diffrho /= 0.) then if (ldensity_nolog.and.(ilnrho /= irho)) & call fatal_error('hyper3_diffrho special','please check') ! call del6(f,ilnrho,hc,IGNOREDX=.true.) df(l1:l2,m,n,ilnrho) = df(l1:l2,m,n,ilnrho) + hyper3_diffrho*hc ! ! due to ignoredx hyper3_diffrho has [1/s] ! if (lfirst.and.ldt) dt1_max=max(dt1_max,hyper3_diffrho/0.01) endif ! call keep_compiler_quiet(p) ! endsubroutine special_calc_density !*********************************************************************** subroutine calc_heatcond_spitzer(f,df,p) ! ! Computes Spitzer heat conduction along the 1D loop ! rhs = Div K T^2.5 Grad(T) ! ! 10-oct-04/bing: coded ! use Diagnostics, only: max_mn_name use EquationOfState, only: gamma use Sub, only: dot2,dot,cubic_step use Slices_methods, only: store_slices ! real, dimension (mx,my,mz,mfarray), intent(in) :: f real, dimension (mx,my,mz,mvar), intent(inout) :: df type (pencil_case), intent(in) :: p real, dimension (nx) :: chi,glnTT2,rhs,u_spitzer real, dimension (nx) :: chi_sat,chi_c,gKc,gKpara,gKsat !,glnTT2_upwind real, dimension (nx,3) :: tmpv2,tmpv !,glnTT_upwind real :: Ksatb,Kcb integer :: i,j ! chi=Kpara*exp(p%lnTT*2.5-p%lnrho)*p%cp1 ! if (init_time /= 0.) chi=chi*cubic_step(real(t),init_time,init_time) ! ! do i=1,3 ! call der_upwind(f,-p%glnTT,ilnTT,glnTT_upwind(:,i),i) ! enddo ! call dot2(glnTT_upwind,glnTT2_upwind) call dot2(p%glnTT,glnTT2) ! gKpara = 3.5 * glnTT2 !_upwind ! if (Ksat /= 0.) then Ksatb = Ksat*7.28d7 /unit_velocity**3. * unit_temperature**1.5 chi_sat = Ksatb * sqrt(p%TT/max(tini,glnTT2))*p%cp1 tmpv(:,:)=0. do i=1,3 do j=1,3 tmpv(:,i)=tmpv(:,i)+p%glnTT(:,j)*p%hlnTT(:,j,i) enddo enddo do i=1,3 tmpv2(:,i) = p%glnrho(:,1) + 1.5*p%glnTT(:,1) - tmpv(:,1)/max(tini,glnTT2) enddo call dot(tmpv2,p%glnTT,gKsat) where (chi > chi_sat) chi = chi_sat gKpara = gKsat endwhere endif ! ! limit the diffusion speed to the speed of light ! if (Kc /= 0.) then Kcb = Kc * dxmin * c_light * cdtv chi_c = Kcb call dot(p%glnTT+p%glnrho,p%glnTT,gKc) where (chi > chi_c) chi = chi_c gKpara = gKc endwhere endif ! rhs = gamma * chi * (gKpara + p%del2lnTT) ! ! Add to energy equation. if (ltemperature) then if (ltemperature_nolog) then call fatal_error('calc_heatcond_spitzer', & 'not implented for ltemperature_nolog') else df(l1:l2,m,n,ilnTT)=df(l1:l2,m,n,ilnTT) + rhs endif else if (lentropy) then if (pretend_lnTT) then df(l1:l2,m,n,ilnTT)=df(l1:l2,m,n,ilnTT) + rhs else df(l1:l2,m,n,iss)=df(l1:l2,m,n,iss) + rhs/(gamma*p%cp1) endif else call fatal_error('calc_heatcond_spitzer', & 'not implented for lentropy and not pretend_lnTT') endif ! if (lfirst.and.ldt) then diffus_chi=diffus_chi+gamma*chi*dxyz_2 ! u_spitzer = 7./2.*gamma*chi*sum(abs(p%glnTT)*dline_1,2) ! ! MR: dt1_max should not be used directly here ! dt1_max=max(dt1_max,u_spitzer/cdt) ! if (ldiagnos.and.idiag_dtspitzer/=0) then call max_mn_name(max(diffus_chi/cdtv,u_spitzer/cdt), & idiag_dtspitzer,l_dt=.true.) endif endif ! if (lvideo) then ! ! slices ! if (ivid_spitzer/=0) & call store_slices(1./chi,spitzer_xy,spitzer_xz,spitzer_yz, & spitzer_xy2,spitzer_xy3,spitzer_xy4,spitzer_xz2,spitzer_r) endif ! endsubroutine calc_heatcond_spitzer !*********************************************************************** subroutine calc_heatcond_kchrom(df,p) ! use Diagnostics, only: max_mn_name use EquationOfState, only: gamma use Sub, only: cubic_step, cubic_der_step, dot2,cross,dot ! real, dimension (mx,my,mz,mvar), intent(inout) :: df type (pencil_case), intent(in) :: p real, dimension(nx) :: K_chrom,glnTT_2 real, dimension(nx) :: gK_chrom,rhs,chix,tmp real, dimension(nx,3) :: tmpv,tmpv2 integer :: i,j ! call dot2(p%glnTT,glnTT_2) ! K_chrom=Kchrom*(1.-cubic_step(glnTT_2,lnTT0_chrom,width_lnTT_chrom)) chix=p%rho1*K_chrom*p%cp1 ! tmpv(:,1) = p%hlnTT(:,2,3)-p%hlnTT(:,3,2) tmpv(:,2) = p%hlnTT(:,3,1)-p%hlnTT(:,1,3) tmpv(:,3) = p%hlnTT(:,1,2)-p%hlnTT(:,2,1) ! call cross(tmpv,p%glnTT,tmpv2) ! tmpv(:,:)=0. do i=1,3 do j=1,3 tmpv(:,i) = tmpv(:,i) + p%glnTT(:,j)*p%hlnTT(:,i,j) enddo tmpv(:,i)=tmpv(:,i) + tmpv2(:,i) enddo ! call dot(2*tmpv,p%glnTT,tmp) ! gK_chrom= - Kchrom*cubic_der_step(glnTT_2,lnTT0_chrom,width_lnTT_chrom) ! rhs = tmp*gK_chrom + glnTT_2*K_chrom + K_chrom*p%del2lnTT ! rhs = p%rho1*p%cp1*rhs ! df(l1:l2,m,n,ilnTT)=df(l1:l2,m,n,ilnTT) + rhs ! ! Check maximum diffusion from thermal diffusion. ! if (lfirst.and.ldt) then diffus_chi=diffus_chi+gamma*chix*dxyz_2 if (ldiagnos.and.idiag_dtchi2/=0.) then call max_mn_name(diffus_chi/cdtv,idiag_dtchi2,l_dt=.true.) endif endif ! endsubroutine calc_heatcond_kchrom !*********************************************************************** subroutine calc_heat_cool_RTV(df,p) ! ! Computes the radiative loss in the optical thin corona. ! Zero order estimation for the electron number densities by ! computing ionization of hydrogen using Saha-equation. ! ! Electron Temperature should be used for the radiative loss ! L = n_e * n_H * Q(T_e) ! ! 04-sep-10/bing: coded ! use EquationOfState, only: gamma use Diagnostics, only: max_mn_name use Messages, only: warning use Sub, only: cubic_step use Slices_methods, only: store_slices ! real, dimension (mx,my,mz,mvar), intent(inout) :: df type (pencil_case), intent(in) :: p ! real, dimension (nx) :: lnQ,rtv_cool,lnTT_SI,lnneni,delta_lnTT real :: unit_lnQ ! unit_lnQ=3*alog(real(unit_velocity))+& 5*alog(real(unit_length))+alog(real(unit_density)) lnTT_SI = p%lnTT + alog(real(unit_temperature)) ! ! calculate ln(ne*ni) : ! ln(ne*ni) = ln( 1.17*rho^2/(1.34*mp)^2) ! lnneni = 2*p%lnrho + alog(1.17) - 2*alog(1.34)-2.*alog(real(m_p)) ! lnneni = 2.*(p%lnrho+61.4412 +alog(real(unit_mass))) ! call getlnQ(lnTT_SI,lnQ,delta_lnTT) ! rtv_cool = lnQ-unit_lnQ+lnneni-p%lnTT-p%lnrho rtv_cool = gamma*p%cp1*exp(rtv_cool) ! if (exp_RTV/=0.) then call warning('cool_RTV','exp_RTV not yet implemented') elseif (tanh_RTV/=0.) then rtv_cool=rtv_cool*cool_RTV* & 0.5*(1.-tanh(width_RTV*(p%lnrho-tanh_RTV))) ! elseif (cubic_RTV/=0.) then rtv_cool=rtv_cool*cool_RTV* & (1.-cubic_step(p%lnrho,cubic_RTV,width_RTV)) ! endif ! if (init_time /= 0.) rtv_cool = rtv_cool * cubic_step(real(t),init_time,init_time) ! ! add to temperature equation ! if (ltemperature) then if (ltemperature_nolog) then call fatal_error('calc_heat_cool_RTV', & 'not implemented for ltemperature_nolog') else df(l1:l2,m,n,ilnTT) = df(l1:l2,m,n,ilnTT)-rtv_cool endif else if (lentropy) then if (pretend_lnTT) then df(l1:l2,m,n,ilnTT) = df(l1:l2,m,n,ilnTT)-rtv_cool else df(l1:l2,m,n,iss) = df(l1:l2,m,n,iss)-rtv_cool/(gamma*p%cp1) endif else call fatal_error('calc_heat_cool_RTV','not implemented for lentropy') endif ! if (lfirst.and.ldt) then dt1_max=max(dt1_max,rtv_cool/cdts) dt1_max=max(dt1_max,abs(rtv_cool/max(tini,delta_lnTT))) if (ldiagnos.and.idiag_dtrad/=0) then itype_name(idiag_dtrad)=ilabel_max_dt call max_mn_name(rtv_cool/cdts,idiag_dtrad,l_dt=.true.) call max_mn_name(abs(rtv_cool/max(tini,delta_lnTT)),idiag_dtrad,l_dt=.true.) endif endif ! ! fill video slices ! if (lvideo) then if (ivid_rtv/=0) & call store_slices(rtv_cool,rtv_xy,rtv_xz,rtv_yz,rtv_xy2,rtv_xy3,rtv_xy4,rtv_xz2,rtv_r) if (ivid_logQ/=0) & call store_slices(lnQ*0.43429448,logQ_xy,logQ_xz,logQ_yz, & logQ_xy2,logQ_xy3,logQ_xy4,logQ_xz2,logQ_r) endif ! endsubroutine calc_heat_cool_RTV !*********************************************************************** subroutine getlnQ(lnTT,get_lnQ,delta_lnTT) ! ! input: lnTT in SI units ! output: lnP [p]= W * m^3 ! real, parameter, dimension (36) :: intlnT = (/ & 8.98008, 9.09521, 9.21034, 9.44060, 9.67086, 9.90112, & 10.1314, 10.2465, 10.3616, 10.5919, 10.8221, 11.0524, & 11.2827, 11.5129, 11.7432, 11.9734, 12.2037, 12.4340,& 12.6642, 12.8945, 13.1247, 13.3550, 13.5853, 13.8155, & 14.0458, 14.2760, 14.5063, 14.6214, 14.7365, 14.8517, & 14.9668, 15.1971, 15.4273, 15.6576, 15.8878, 16.1181 /) real, parameter, dimension (36) :: intlnQ = (/ & -83.9292, -82.8931, -82.4172, -81.2275, -80.5291, -80.0532, & -80.1837, -80.2067, -80.1837, -79.9765, -79.6694, -79.2857,& -79.0938, -79.1322, -79.4776, -79.4776, -79.3471, -79.2934,& -79.5159, -79.6618, -79.4776, -79.3778, -79.4008, -79.5159,& -79.7462, -80.1990, -80.9052, -81.3196, -81.9874, -82.2023,& -82.5093, -82.5477, -82.4172, -82.2637, -82.1793 , -82.2023 /) ! real, dimension (nx), intent(in) :: lnTT real, dimension (nx), intent (out) :: get_lnQ,delta_lnTT real :: slope,ordinate integer :: i,j=18 logical :: notdone ! get_lnQ=-1000. delta_lnTT = 100. ! do i=1,nx notdone=.true. ! do while (notdone) ! if (lnTT(i) >= intlnT(j) .and. lnTT(i) < intlnT(j+1)) then ! ! define slope and ordinate for linear interpolation ! slope=(intlnQ(j+1)-intlnQ(j))/(intlnT(j+1)-intlnT(j)) ordinate=intlnQ(j) - slope*intlnT(j) ! get_lnQ(i) = slope*lnTT(i) + ordinate delta_lnTT(i) = intlnT(j+1) - intlnT(j) notdone = .false. else j = j + sign(1.,lnTT(i)-intlnT(j)) if (j <= 0) then j=1 notdone=.false. elseif (j >= 35) then j = 35 !call fatal_error('get_lnQ','lnTT to large') notdone=.false. endif endif enddo enddo ! endsubroutine getlnQ !*********************************************************************** subroutine calc_heat_cool_newton(df,p) ! ! newton cooling dT/dt = -1/tau * (T-T0) ! dlnT/dt = 1/tau *(1-T0/T-1) ! where ! tau = (rho0/rho)^alpha ! ! 13-sep-10/bing: coded ! use Diagnostics, only: max_mn_name use EquationOfState, only: lnrho0 use Sub, only: cubic_step ! real, dimension (mx,my,mz,mvar), intent(inout) :: df type (pencil_case), intent(in) :: p ! real, dimension (nx) :: newton,tau_inv_tmp ! if (headtt) print*,'special_calc_energy: newton cooling',tau_inv_newton ! ! Get reference temperature newton = exp(lnTT_init_prof(l1:l2)-p%lnTT)-1. ! ! Multiply by density dependend time scale if (exp_newton/=0.) then tau_inv_tmp = tau_inv_newton * & exp(-exp_newton*(lnrho0-p%lnrho)) ! elseif (tanh_newton/=0.) then tau_inv_tmp = tau_inv_newton * & 0.5*(1+tanh(width_newton*(p%lnrho-tanh_newton))) ! elseif (cubic_newton/=0.) then tau_inv_tmp = tau_inv_newton * & cubic_step(p%lnrho,cubic_newton,width_newton) ! endif ! ! Adjust time scale by the initialization time if (init_time /= 0.) tau_inv_tmp = tau_inv_tmp * cubic_step(real(t),init_time,init_time) ! newton = newton * tau_inv_tmp ! ! Add newton cooling term to entropy ! df(l1:l2,m,n,ilnTT) = df(l1:l2,m,n,ilnTT) + newton ! if (lfirst.and.ldt) then dt1_max=max(dt1_max,tau_inv_tmp/cdts) if (ldiagnos.and.idiag_dtnewt/=0.) then call max_mn_name(tau_inv_tmp/cdts,idiag_dtnewt,l_dt=.true.) endif endif ! endsubroutine calc_heat_cool_newton !*********************************************************************** subroutine calc_artif_heating(df,p) ! ! Subroutine to calculate intrinsic heating. ! Activated by setting iheattype = exp, exp2 and/or gauss ! Maximum of 3 different possibible heating types ! Also set the heating parameters for the chosen heating functions. ! ! 22-sept-10/Tijmen: coded ! use EquationOfState, only: gamma use Diagnostics, only: max_mn_name use General, only: random_number_wrapper,random_seed_wrapper, & normal_deviate,notanumber use Sub, only: cubic_step, dot2 ! real, dimension (mx,my,mz,mvar) :: df real, dimension (nx) :: heatinput,doors_heat,doors_a real, dimension (nx) :: x_Mm,heat_nano,rhs,height real, dimension (nx) :: heat_event,heat_event1D,Hlength real, dimension (nx) :: Htemp real, dimension (15,2) :: a_arr integer, dimension(mseed) :: global_rstate real :: heat_unit,heat_flux real :: nano_sigma_t,nano_time,nano_start,nano_sigma_z real :: nano_flare_energy,nano_pos_x,nano_pos_z,nano_pos_y real :: nano_amplitude real :: slope real, dimension(2) :: event_pos type (pencil_case) :: p integer :: i,j,k logical :: notdone ! ! the unit W/m^3 is heat_unit= unit_density*unit_velocity**3/unit_length ! ! coordinate along the loop in [Mm] x_Mm = x(l1:l2)*unit_length*1e-6 ! ! The height in [Mm] height = Ltot/pi*sin(x(l1:l2)/Ltot*pi)*unit_length*1e-6 ! heatinput = 0. heat_flux = 0. ! do i=1,3 if (headtt) print*,'iheattype:',iheattype(i) select case(iheattype(i)) case ('nothing') ! do nothing case ('one-sided') heatinput=heatinput + & heat_par_exp(1)*exp(-x(l1:l2)/heat_par_exp(2))/heat_unit ! heatinput=heatinput + & heat_par_exp2(1)*exp(-x(l1:l2)/heat_par_exp2(2))/heat_unit ! heat_flux=heat_flux - heat_par_exp(1)*heat_par_exp(2)*1e6 ! if (headtt) print*,'Flux of exp heating: ',heat_flux,' [ W m^(-2)]' ! case ('exp') heatinput=heatinput + & heat_par_exp(1)/heat_unit*exp(-height/heat_par_exp(2)) ! ! add the flux at a ref height of 4 Mm ! if (headtt) then heat_flux = heat_par_exp(1)*heat_par_exp(2)*1e6*exp(-4./heat_par_exp(2)) write (*,*) 'Exp heating, scale height: ',heat_par_exp(2),' [ Mm ]' write (*,*) 'Exp heating, amplitude',heat_par_exp(1)/heat_unit,' [ Wm^-3 ]' write (*,*) 'Exp heating, flux at z=4 Mm: ', heat_flux write (*,*) 'Exp heating, flux at z=0 Mm: ', heat_par_exp(1)*heat_par_exp(2)*1e6 endif ! case ('exp2') ! A second exponential function ! For Sven et al. 2010 set: ! heat_par_exp= (10 , 0.2 ) ! heat_par_exp2= (1e-4 , 4.) ! heatinput=heatinput + & heat_par_exp2(1)/heat_unit*exp(-height/heat_par_exp2(2)) ! ! add the flux at a ref height of 4 Mm ! if (headtt) then heat_flux = heat_par_exp2(1)*heat_par_exp2(2)*1e6*exp(-4./heat_par_exp2(2)) write (*,*) 'Exp heating, scale height: ',heat_par_exp(2),' [ Mm ]' write (*,*) 'Exp heating, amplitude',heat_par_exp2(1)/heat_unit,' [ Wm^-3 ]' write (*,*) 'Exp heating, flux at z=4 Mm: ', heat_flux write (*,*) 'Exp heating, flux at z=0 Mm: ', heat_par_exp2(1)*heat_par_exp2(2)*1e6 endif ! case ('gauss') ! heat_par_gauss(1) is Center (z in Mm) ! heat_par_gauss(2) is Width (sigma) ! heat_par_gauss(3) is the amplitude (Flux) ! heatinput=heatinput + & heat_par_gauss(3)*exp(-((x_Mm-heat_par_gauss(1))**2/ & (2*heat_par_gauss(2)**2)))/heat_unit ! case ('nanof') ! simulate nanoflare heating =) ! height dependend amplitude und duration ! idea: call random numbers , make condition when to flare, when ! flare get position ! then over time release a flare in shape of gaussian. Also ! distribution is gaussian around a point. ! gaussian distribution is gained via the dierfc function (double ! prec, inverser erf function) ! we draw for new flares as long as nano_time is /= 0. when done we ! reset nano_time to 0. again. ! ! Later we will implement more options for nanoflaring! =D ! ! SAVE GLOBAL SEED ! LOAD NANO_SEED call random_seed_wrapper(GET=global_rstate) call random_seed_wrapper(PUT=nano_seed) ! if (nano_start .eq. 0.) then ! If 0 roll to see if a nanoflare occurs call random_number_wrapper(nano_start) ! if (nano_start .gt. 0.95 ) then ! 5% chance for a nanoflare to occur, then get the location. call normal_deviate(nano_pos_z) nano_pos_z=nano_pos_z*lz call random_number_wrapper(nano_pos_y) call random_number_wrapper(nano_pos_x) nano_time=60. else ! else no nanoflare, reset nano_start to 0. for the next ! timestep nano_start=0. endif endif ! if (nano_start .ne. 0.) then ! if nano_start is not 0. then there is a nanoflare! ! 2nd assumption , nanoflare takes 60 seconds =) nano_flare_energy=10.d17 !joules nano_sigma_z=0.5 nano_sigma_t=2.5 ! nano_amplitude=nano_flare_energy/(pi/2*nano_sigma_t*nano_sigma_z*1.d6 ) ! heat_nano=nano_amplitude*exp(-((nano_time-5.))**2/( 2*2.**2))* & exp(-((x_Mm-nano_pos_z)**2/ (2*0.5**2))) nano_time=nano_time-dt*unit_time ! heatinput=heatinput + heat_nano/heat_unit ! if (nano_time .le. 0.) nano_start=0. end if ! !SAVE NANO_SEED !RESTORE GLOBAL SEED call random_seed_wrapper(GET=nano_seed) call random_seed_wrapper(PUT=global_rstate) ! ! amp_t = exp(-((t-nano_time)**2/(2.*nano_dur**2))) !-> no idea how to properly implement ! spread = exp(-((x_Mm-nano_pos)**2/(2.*nano_spread**2))) ! ! issue! How to put in timelike guassian for a random event starting ! at a random time? ! case ('event') ! one small point heating event (gaussian to prevent too strong gradients) ! one point in time, one point in space! if (t*unit_time .gt. 150. .AND. t*unit_time .lt. 1000.) then event_pos(1)=7.5 event_pos(2)=15. heat_event=10.*exp(-((250.-t*unit_time))**2/(2*(20.*unit_time)**2))* & exp(-((x_Mm-event_pos(1))**2/ (2*0.2**2))) * & exp(-((x_Mm-event_pos(2))**2/ (2*0.2**2))) heatinput=heatinput + heat_event/heat_unit endif ! case ('event1D') ! one small point heating event (gaussian to prevent too strong gradients) ! one point in time, one point in space! if (t*unit_time .gt. 300. .AND. t*unit_time .lt. 10000.) then if (t*unit_time .gt. 300. .AND. t*unit_time .lt. 301.) & print*,'EVENTTTT!!!!!' event_pos(1)=10. heat_event1D=10.*exp(-((400.-t))**2/( 2*50.**2))* & exp(-((x_Mm-event_pos(1))**2/ (2*0.2**2))) heatinput=heatinput + heat_event1D/heat_unit endif ! !--------------------------------------------------------------------- case ('vandoors') !heating based on "Van Doorsselaere et al. 2007" !pi defined? ! heat_par_vandoors(1)=1. !amplitude ! heat_par_vandoors(2)=0.1 !(temp) some amplitude to add a physical heating ! heat_par_vandoors(3)=50. Mm, !density scale height, later dynamic? ! if (headtt) print*,heat_par_vandoors,pi ! ! --find a -- !better in the header of the subroutine? a_arr(:,1)=(/0.5,0.6,0.7,0.8,0.9,1.,1.2,1.4,1.6,2.,2.4,2.8,3.2,3.6,1000. /) !L/pi/H a_arr(:,2)=(/0.03426,0.04169,0.04933,0.05717,0.06522,0.07349,0.09068,0.10877,& 0.12778,0.16853,0.21287,0.26046,0.31068,0.36260,125./) !a ! if (notanumber(p%glnrho)) then ! print*,p%glnrho call fatal_error('p%glnrho','NaN found') print*,'-------------------------------' endif ! call dot2(p%glnrho,Htemp,FAST_SQRT=.true.) ! Htemp=sqrt(Htemp**2.) Hlength=(Htemp+tiny(0D0))**(-1.)/Ltot/pi ! j=6 ! do k=1,nx notdone=.true. ! do while (notdone) if (Hlength(k) >= a_arr(j,1) .and. Hlength(k)= 14) then j = 14 notdone=.false. endif endif enddo slope=(a_arr(j+1,2)-a_arr(j,2))/(a_arr(j+1,1)-a_arr(j,1)) doors_a(k)=slope*(Hlength(k)-a_arr(j,1))+a_arr(j,2) enddo ! if (notanumber(doors_a)) then print*,Hlength call fatal_error('update points','NaN found') print*,'-------------------------------' endif ! where(doors_a < 0) doors_a=0D0 endwhere where(doors_a > 0.5) doors_a=0.5 endwhere ! doors_heat= heat_par_vandoors(1)* (cos(pi* height/Ltot)**2.+ & 6.*doors_a*cos(pi*height/Ltot)*cos(3.*pi*height/Ltot)) where(doors_heat < 0) doors_heat=0 endwhere ! do k=0,nx if (doors_heat(k) >1. ) then print*,doors_heat(k),k call fatal_error('doors_heat','toobig!') endif enddo ! heatinput=heatinput + doors_heat!/heat_unit !print*,max(transpose(heat_par_vandoors(1)*cos(pi* height/Ltot)**2.) ) ! ! print*,max((heat_par_vandoors(1)*cos(pi* height/ L)**2.+6.*heat_par_vandoors(2)*cos(pi*height/L)*cos(3.*pi*height/L))/& ! heat_unit) !--------------------------------------------------------------------- case default if (headtt) call fatal_error('calc_artif_heating', & 'Please provide correct iheattype') endselect enddo ! if (headtt) print*,'Total flux for all types:',heat_flux ! ! Add to energy equation ! rhs = p%TT1*p%rho1*gamma*p%cp1*heatinput if (init_time /=0.) rhs=rhs*cubic_step(real(t),init_time,init_time) ! if (ltemperature .and. (.not. ltemperature_nolog)) then df(l1:l2,m,n,ilnTT) = df(l1:l2,m,n,ilnTT) + rhs else if (lentropy .and. (.not. pretend_lnTT)) then df(l1:l2,m,n,iss) = df(l1:l2,m,n,iss) + rhs/gamma/p%cp1 else call fatal_error('calc_artif_heating', & 'not for current set of thermodynamic variables') endif ! if (lfirst.and.ldt) then if (ldiagnos.and.idiag_dtnewt/=0) then ! call max_mn_name(rhs/cdts,idiag_dtnewt,l_dt=.true.) endif dt1_max=max(dt1_max,rhs/cdts) endif ! endsubroutine calc_artif_heating !*********************************************************************** subroutine calc_heatcond_glnTT_iso(df,p) ! ! L = Div( Grad(lnT)^2 Grad(T)) ! use Diagnostics, only : max_mn_name use Sub, only : dot2,dot,multsv,multmv use EquationOfState, only : gamma ! real, dimension (mx,my,mz,mvar) :: df type (pencil_case) :: p real, dimension (nx,3) :: tmpv real, dimension (nx) :: glnT2,glnT_glnr real, dimension (nx) :: tmp,rhs,chi integer :: i ! intent(in) :: p intent(out) :: df ! call dot2(p%glnTT,glnT2) call dot(p%glnTT,p%glnrho,glnT_glnr) ! do i=1,3 tmpv(:,i) = p%glnTT(:,1)*p%hlnTT(:,1,i) + & p%glnTT(:,2)*p%hlnTT(:,2,i) + & p%glnTT(:,3)*p%hlnTT(:,3,i) enddo call dot(p%glnTT,tmpv,tmp) ! chi = glnT2*hcond_grad_iso ! rhs = 2*tmp+glnT2*(glnT2+p%del2lnTT+glnT_glnr) ! if (ltemperature .and. (.not. ltemperature_nolog)) then df(l1:l2,m,n,ilnTT)=df(l1:l2,m,n,ilnTT) + p%cp1*rhs*gamma*hcond_grad_iso else if (lentropy .and. (.not. pretend_lnTT)) then df(l1:l2,m,n,iss)=df(l1:l2,m,n,iss) + rhs*hcond_grad_iso else call fatal_error('calc_heatcond_glnTT_iso','only for ltemperature') endif ! if (lfirst.and.ldt) then diffus_chi=diffus_chi+gamma*chi*dxyz_2 if (ldiagnos.and.idiag_dtchi2/=0.) then call max_mn_name(diffus_chi/cdtv,idiag_dtchi2,l_dt=.true.) endif endif ! endsubroutine calc_heatcond_glnTT_iso !*********************************************************************** subroutine calc_heatcond_constchi(df,p) ! ! field-aligned heat conduction that is proportional to rho ! L = Div (b K rho b*Grad(T)) ! K = chi [m^2/s] * cV [J/kg/K] = hcond1 [m^4/s^3/K] ! Define chi = K_0/rho ! use Diagnostics, only: max_mn_name use EquationOfState, only: gamma use Sub ! real, dimension(mx,my,mz,mvar), intent(inout) :: df type (pencil_case), intent(in) :: p ! real, dimension(nx) :: rhs, glnrho_b, fdiff real, dimension(nx) :: glnTT_H,tmp,hlnTT_Bij,glnTT_b real, dimension(nx,3) :: hhh,tmpv integer :: i,j,k ! call dot(p%glnrho,p%bunit,glnrho_b) call dot(p%glnTT,p%bunit,glnTT_b) ! do i = 1,3 hhh(:,i) = 0. do j = 1,3 tmp(:) = 0. do k = 1,3 tmp(:) = tmp(:)-2.*p%bunit(:,k)*p%bij(:,k,j) enddo hhh(:,i) = hhh(:,i)+p%bunit(:,j)*(p%bij(:,i,j)+p%bunit(:,i)*tmp(:)) enddo enddo call dot(hhh,p%glnTT,glnTT_H) ! call multmv(p%hlnTT,p%bunit,tmpv) call dot(tmpv,p%bunit,hlnTT_Bij) ! rhs = (glnTT_H + hlnTT_Bij + (glnrho_b + glnTT_b)*glnTT_b)*chi_aniso*gamma ! if (ltemperature) then df(l1:l2,m,n,ilnTT) = df(l1:l2,m,n,ilnTT) + rhs else call fatal_error('calc_heatcond_constchi special','only for ltemperature') endif ! ! ANISOTROPIC PART NOT INCLUDED IN THE TIMESTEP CALCULATION if (lfirst .and. ldt) then advec_cs2 = max(advec_cs2,chi_aniso*maxval(dxyz_2)) fdiff = gamma * chi_aniso * dxyz_2 diffus_chi = diffus_chi+fdiff if (ldiagnos .and. (idiag_dtchi2 /= 0)) then call max_mn_name(fdiff/cdtv,idiag_dtchi2,l_dt=.true.) endif endif ! endsubroutine calc_heatcond_constchi !*********************************************************************** subroutine der_upwind(f,uu,k,df,j) ! ! High order upwind derivative of variable. ! Useful for advecting. ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (nx,3) :: uu real, dimension (nx) :: df,fac integer :: j,k,l ! intent(in) :: f,uu,k,j intent(out) :: df ! if (nghost /= 3) call fatal_error('der_upwind','only for nghost==3') if (lspherical_coords.or.lcylindrical_coords) & call fatal_error('der_upwind','NOT IMPLEMENTED for non-cartesian grid') ! if (j == 1) then if (nxgrid /= 1) then fac = 1./6.*dx_1(l1:l2) do l=1,nx if (uu(l,1) > 0.) then df(l) = 11.*f(nghost+l ,m,n,k)-18.*f(nghost+l-1,m,n,k) & + 9.*f(nghost+l-2,m,n,k)- 2.*f(nghost+l-3,m,n,k) else df(l) = 2.*f(nghost+l+3,m,n,k)- 9.*f(nghost+l+2,m,n,k) & +18.*f(nghost+l+1,m,n,k)-11.*f(nghost+l ,m,n,k) endif enddo df = fac*df else df=0. if (ip<=5) print*, 'der_upwind: Degenerate case in x-direction' endif elseif (j == 2) then if (nygrid /= 1) then fac = 1./6.*dy_1(m) do l=1,nx if (uu(l,2) > 0.) then df(l) = 11.*f(nghost+l,m ,n,k)-18.*f(nghost+l,m-1,n,k) & + 9.*f(nghost+l,m-2,n,k)- 2.*f(nghost+l,m-3,n,k) else df(l) = 2.*f(nghost+l,m+3,n,k)- 9.*f(nghost+l,m+2,n,k) & +18.*f(nghost+l,m+1,n,k)-11.*f(nghost+l,m ,n,k) endif enddo df = fac*df else df=0. if (ip<=5) print*, 'der_upwind: Degenerate case in y-direction' endif elseif (j == 3) then if (nzgrid /= 1) then fac = 1./6.*dz_1(n) do l=1,nx if (uu(l,3) > 0.) then df(l) = 11.*f(nghost+l,m,n ,k)-18.*f(nghost+l,m,n-1,k) & + 9.*f(nghost+l,m,n-2,k)- 2.*f(nghost+l,m,n-3,k) else df(l) = 2.*f(nghost+l,m,n+3,k)- 9.*f(nghost+l,m,n+2,k) & +18.*f(nghost+l,m,n+1,k)-11.*f(nghost+l,m,n ,k) endif enddo df = fac*df else df=0. if (ip<=5) print*, 'der_upwind: Degenerate case in z-direction' endif endif ! endsubroutine der_upwind !*********************************************************************** subroutine calc_bullets_heating(df,p) ! use EquationOfState, only : gamma ! real, dimension (mx,my,mz,mvar) :: df type (pencil_case) :: p ! df(l1:l2,m,n,ilnTT) = df(l1:l2,m,n,ilnTT) + bullets_h0* & exp( - ((x(l1:l2)-bullets_x0)/bullets_dx)**2.)* & exp( - ((t -bullets_t0)/bullets_dt)**2.)* & p%cp1*gamma*exp(-p%lnrho -p%lnTT) ! endsubroutine calc_bullets_heating !*********************************************************************** subroutine filter_farray(f,df) ! ! Reduce noise of farray using mesh independend hyper diffusion. ! Filter strength is equal for all variables up to now and has ! to be smaller than 1/64 for numerical stability. ! ! 13-jan-12/bing: coded ! use Sub, only: del6 ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz,mvar) :: df real, dimension (nx) :: del6_fj integer :: j ! if (dt/=0.) then do j=1,mvar if (filter_strength(j)/=0.) then call del6(f,j,del6_fj,IGNOREDX=.true.) df(l1:l2,m,n,j) = df(l1:l2,m,n,j) + & filter_strength(j)*del6_fj/dt_beta_ts(itsub) endif enddo endif ! endsubroutine filter_farray !*********************************************************************** !******************************************************************** !************ 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