! $Id$ ! ! This module takes care of entropy (initial condition ! and time advance) for a fluid consisting of gas and perfectly ! coupled pressureless dust. ! !** 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 :: lentropy = .true. ! CPARAM logical, parameter :: ltemperature = .false. ! CPARAM logical, parameter :: lthermal_energy = .false. ! ! MVAR CONTRIBUTION 1 ! MAUX CONTRIBUTION 0 ! ! PENCILS PROVIDED ss; gss(3); ee; pp; lnTT; cs2; cp1tilde; glnTT(3) ! PENCILS PROVIDED TT; TT1; Ma2; ugss; hss(3,3); hlnTT(3,3) ! PENCILS PROVIDED del2ss; del6ss; del2lnTT; cv1; fpres(3); sglnTT(3) ! !*************************************************************** module Energy ! use Cparam use Cdata use General, only: keep_compiler_quiet use EquationOfState, only: gamma, gamma_m1, cs20 use Density, only: beta_glnrho_global use Interstellar use Messages use Viscosity ! implicit none ! include 'energy.h' ! real :: ss_const=0.0 real :: T0=1.0 logical :: lupw_ss=.false. logical, pointer :: lpressuregradient_gas logical :: lviscosity_heat=.true. logical :: ladvection_energy=.true. logical :: lenergy_slope_limited=.false. character (len=labellen), dimension(ninit) :: initss='nothing' character (len=intlen) :: iinit_str integer :: pushpars2c, pushdiags2c ! should be procedure pointer (F2003) ! namelist /entropy_init_pars/ & initss, ss_const, T0, beta_glnrho_global, ladvection_energy ! namelist /entropy_run_pars/ & lupw_ss, beta_glnrho_global, ladvection_energy ! integer :: idiag_dtc=0,idiag_ethm=0,idiag_ethdivum=0,idiag_ssm=0 integer :: idiag_ugradpm=0,idiag_ethtot=0,idiag_ssmphi=0 integer :: idiag_yHm=0,idiag_yHmax=0,idiag_TTm=0,idiag_TTmax=0,idiag_TTmin=0 integer :: idiag_fconvz=0,idiag_dcoolz=0,idiag_fradz=0,idiag_fturbz=0 integer :: idiag_ssmz=0,idiag_ssmy=0,idiag_ssmx=0,idiag_TTmz=0 ! contains !*********************************************************************** subroutine register_energy ! ! Initialise variables which should know that we solve an entropy ! equation: iss, etc; increase nvar accordingly. ! ! 6-nov-01/wolf: coded ! use FArrayManager use SharedVariables, only: get_shared_variable ! integer :: ierr ! call farray_register_pde('ss',iss) ! ! Identify version number. ! if (lroot) call svn_id( & "$Id$") ! ! logical variable lpressuregradient_gas shared with hydro modules ! call get_shared_variable('lpressuregradient_gas',lpressuregradient_gas,ierr) if (ierr/=0) call fatal_error('register_energy','lpressuregradient_gas') ! endsubroutine register_energy !*********************************************************************** subroutine initialize_energy(f) ! ! Called by run.f90 after reading parameters, but before the time loop. ! ! 21-jul-2002/wolf: coded ! use EquationOfState use Gravity, only: gravz,g0 use SharedVariables, only: put_shared_variable ! real, dimension (mx,my,mz,mfarray) :: f ! integer :: i logical :: lnothing ! ! check any module dependencies. ! if (.not. leos) then call fatal_error('initialize_energy','EOS=noeos but entropy requires an EQUATION OF STATE for the fluid') endif call select_eos_variable('ss',iss) ! ! For global density gradient beta=H/r*dlnrho/dlnr, calculate actual ! gradient dlnrho/dr = beta/H. ! if (maxval(abs(beta_glnrho_global))/=0.0) then beta_glnrho_scaled=beta_glnrho_global*Omega/cs0 if (lroot) print*, 'initialize_energy: Global density gradient '// & 'with beta_glnrho_global=', beta_glnrho_global endif ! ! Turn off pressure gradient term and advection for 0-D runs. ! if (nxgrid*nygrid*nzgrid==1) then lpressuregradient_gas=.false. ladvection_energy=.false. print*, 'initialize_energy: 0-D run, turned off pressure gradient term' print*, 'initialize_energy: 0-D run, turned off advection of entropy' endif ! call put_shared_variable('lviscosity_heat',lviscosity_heat) ! call keep_compiler_quiet(f) ! endsubroutine initialize_energy !*********************************************************************** subroutine read_energy_init_pars(iostat) ! use File_io, only: parallel_unit ! integer, intent(out) :: iostat ! read(parallel_unit, NML=entropy_init_pars, IOSTAT=iostat) ! endsubroutine read_energy_init_pars !*********************************************************************** subroutine write_energy_init_pars(unit) ! integer, intent(in) :: unit ! write(unit, NML=entropy_init_pars) ! endsubroutine write_energy_init_pars !*********************************************************************** subroutine read_energy_run_pars(iostat) ! use File_io, only: parallel_unit ! integer, intent(out) :: iostat ! read(parallel_unit, NML=entropy_run_pars, IOSTAT=iostat) ! endsubroutine read_energy_run_pars !*********************************************************************** subroutine write_energy_run_pars(unit) ! integer, intent(in) :: unit ! write(unit, NML=entropy_run_pars) ! endsubroutine write_energy_run_pars !*********************************************************************** subroutine init_energy(f) ! ! Initialise entropy; called from start.f90. ! use Sub use Gravity use General, only: itoa use Initcond use InitialCondition, only: initial_condition_ss use EquationOfState, only: rho0, lnrho0, isothermal_entropy, & isothermal_lnrho_ss, eoscalc, ilnrho_pp ! real, dimension (mx,my,mz,mfarray) :: f ! integer :: j logical :: lnothing ! intent(inout) :: f ! lnothing=.true. do j=1,ninit ! if (initss(j)/='nothing') then ! lnothing=.false. iinit_str=itoa(j) ! ! select different initial conditions ! select case (initss(j)) ! case ('zero', '0'); f(:,:,:,iss) = 0. case ('const_ss'); f(:,:,:,iss)=f(:,:,:,iss)+ss_const case ('isothermal'); call isothermal_entropy(f,T0) case ('isothermal_lnrho_ss') print*, 'init_energy: Isothermal density and entropy stratification' call isothermal_lnrho_ss(f,T0,rho0) case default ! ! Catch unknown values ! write(unit=errormsg,fmt=*) 'No such value for initss(' & //trim(iinit_str)//'): ',trim(initss(j)) call fatal_error('init_energy',errormsg) endselect ! endif ! enddo ! ! Interface for user's own initial condition ! if (linitial_condition) call initial_condition_ss(f) ! endsubroutine init_energy !*********************************************************************** subroutine pencil_criteria_energy ! ! All pencils that the Entropy module depends on are specified here. ! ! 20-11-04/anders: coded ! use Density, only: beta_glnrho_scaled ! if (ldt) lpenc_requested(i_cs2)=.true. if (lpressuregradient_gas) then lpenc_requested(i_cs2)=.true. lpenc_requested(i_cp1tilde)=.true. lpenc_requested(i_glnrho)=.true. lpenc_requested(i_gss)=.true. endif if (ladvection_energy) lpenc_requested(i_ugss)=.true. if (pretend_lnTT) lpenc_requested(i_divu)=.true. if (lpressuregradient_gas) lpenc_requested(i_cp1tilde)=.true. ! if (maxval(abs(beta_glnrho_scaled))/=0.0) lpenc_requested(i_cs2)=.true. ! lpenc_diagnos2d(i_ss)=.true. ! if (idiag_ethdivum/=0) lpenc_diagnos(i_divu)=.true. if (idiag_ssm/=0 .or. idiag_ssmz/=0 .or. idiag_ssmy/=0.or.idiag_ssmx/=0) & lpenc_diagnos(i_ss)=.true. if (idiag_ethm/=0 .or. idiag_ethtot/=0 .or. idiag_ethdivum/=0) then lpenc_diagnos(i_rho)=.true. lpenc_diagnos(i_ee)=.true. endif if (idiag_fconvz/=0 .or. idiag_fturbz/=0 ) then lpenc_diagnos(i_rho)=.true. lpenc_diagnos(i_TT)=.true. !(to be replaced by enthalpy) endif if (idiag_TTm/=0 .or. idiag_TTmz/=0 .or. idiag_TTmax/=0 & .or. idiag_TTmin/=0) & lpenc_diagnos(i_TT)=.true. if (idiag_yHm/=0 .or. idiag_yHmax/=0) lpenc_diagnos(i_yH)=.true. if (idiag_dtc/=0) lpenc_diagnos(i_cs2)=.true. ! endsubroutine pencil_criteria_energy !*********************************************************************** subroutine pencil_interdep_energy(lpencil_in) ! ! Interdependency among pencils from the Entropy module is specified here. ! ! 20-11-04/anders: coded ! logical, dimension(npencils) :: lpencil_in ! if (lpencil_in(i_del2lnTT)) then lpencil_in(i_del2lnrho)=.true. lpencil_in(i_del2ss)=.true. endif if (lpencil_in(i_glnTT)) then lpencil_in(i_glnrho)=.true. lpencil_in(i_gss)=.true. endif if (lpencil_in(i_TT)) lpencil_in(i_lnTT)=.true. if (lpencil_in(i_TT1)) lpencil_in(i_lnTT)=.true. if (lpencil_in(i_ugss)) then lpencil_in(i_uu)=.true. lpencil_in(i_gss)=.true. endif if (pretend_lnTT .and. lpencil_in(i_glnTT)) lpencil_in(i_gss)=.true. if (lpencil_in(i_Ma2)) then lpencil_in(i_u2)=.true. lpencil_in(i_cs2)=.true. endif if (lpencil_in(i_hlnTT)) then if (pretend_lnTT) then lpencil_in(i_hss)=.true. else lpencil_in(i_hlnrho)=.true. lpencil_in(i_hss)=.true. endif endif ! The pencils cs2 and cp1tilde come in a bundle, so enough to request one. if (lpencil_in(i_cs2) .and. lpencil_in(i_cp1tilde)) & lpencil_in(i_cp1tilde)=.false. ! endsubroutine pencil_interdep_energy !*********************************************************************** subroutine calc_pencils_energy(f,p) ! ! Calculate Entropy pencils. ! Most basic pencils should come first, as others may depend on them. ! ! 20-11-04/anders: coded ! 16-05-12/MR: dead branch for pressure force calculation added ! use EquationOfState use Sub ! real, dimension (mx,my,mz,mfarray) :: f type (pencil_case) :: p ! real, dimension (nx,3) :: glnrhod, glnrhotot real, dimension (nx) :: eps integer :: i ! intent(in) :: f intent(inout) :: p ! ss if (lpencil(i_ss)) p%ss=f(l1:l2,m,n,iss) ! gss if (lpencil(i_gss)) call grad(f,iss,p%gss) ! pp if (lpencil(i_pp)) call eoscalc(f,nx,pp=p%pp) ! ee if (lpencil(i_ee)) call eoscalc(f,nx,ee=p%ee) ! lnTT if (lpencil(i_lnTT)) call eoscalc(f,nx,lnTT=p%lnTT) ! TT if (lpencil(i_TT)) p%TT=exp(p%lnTT) ! TT1 if (lpencil(i_TT1)) p%TT1=exp(-p%lnTT) ! cs2 and cp1tilde if (lpencil(i_cs2) .or. lpencil(i_cp1tilde)) & call pressure_gradient(f,p%cs2,p%cp1tilde) ! Ma2 if (lpencil(i_Ma2)) p%Ma2=p%u2/p%cs2 ! glnTT if (lpencil(i_glnTT)) then if (pretend_lnTT) then p%glnTT=p%gss else call temperature_gradient(f,p%glnrho,p%gss,p%glnTT) endif endif ! ugss if (lpencil(i_ugss)) & call u_dot_grad(f,iss,p%gss,p%uu,p%ugss,UPWIND=lupw_ss) !ajwm Should probably combine the following two somehow. ! hss if (lpencil(i_hss)) then call g2ij(f,iss,p%hss) endif ! del2ss if (lpencil(i_del2ss)) then call del2(f,iss,p%del2ss) endif ! del2lnTT if (lpencil(i_del2lnTT)) then call temperature_laplacian(f,p) endif ! del6ss if (lpencil(i_del6ss)) then call del6(f,iss,p%del6ss) endif ! hlnTT if (lpencil(i_hlnTT)) then if (pretend_lnTT) then p%hlnTT=p%hss else call temperature_hessian(f,p%hlnrho,p%hss,p%hlnTT) endif endif ! uud (for dust continuity equation in one fluid approximation) if (lpencil(i_uud)) p%uud(:,:,1)=p%uu ! divud (for dust continuity equation in one fluid approximation) if (lpencil(i_divud)) p%divud(:,1)=p%divu ! if (lpencil(i_fpres)) & call fatal_error('calc_pencils_energy', & 'calculation of pressure force not yet implemented'//& ' for entropy_onefluid') ! sglnTT if (lpencil(i_sglnTT)) & call fatal_error('calc_pencils_energy', & 'Pencil sglnTT not yet implemented for entropy_onefluid') ! endsubroutine calc_pencils_energy !*********************************************************************** subroutine denergy_dt(f,df,p) ! ! Calculate right hand side of entropy equation. ! ! use Conductivity, only: heat_conductivity use Diagnostics use Density, only: beta_glnrho_global, beta_glnrho_scaled use Special, only: special_calc_energy use Sub ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz,mvar) :: df type (pencil_case) :: p ! real, dimension (nx) :: Hmax integer :: j,ju ! intent(inout) :: f,p intent(out) :: df ! Hmax = 0.0 ! ! Identify module and boundary conditions. ! if (headtt.or.ldebug) print*,'denergy_dt: SOLVE denergy_dt' if (headtt) call identify_bcs('ss',iss) ! if (lhydro) then ! ! Pressure term in momentum equation (setting lpressuregradient_gas to ! .false. allows suppressing pressure term for test purposes). ! if (lpressuregradient_gas) then do j=1,3 ju=j+iuu-1 df(l1:l2,m,n,ju) = df(l1:l2,m,n,ju) - & p%cs2*(p%glnrho(:,j) + p%cp1tilde*p%gss(:,j)) enddo ! ! Add pressure force from global density gradient. ! WARNING (AJ): This may be implemented inconsistently, since we have ! here linearised rho and P independently. ! if (maxval(abs(beta_glnrho_global))/=0.0) then if (headtt) print*, 'denergy_dt: adding global pressure gradient force' do j=1,3 df(l1:l2,m,n,(iux-1)+j) = df(l1:l2,m,n,(iux-1)+j) & - p%cs2*beta_glnrho_scaled(j) enddo endif endif endif ! ! Advection term. ! if (ladvection_energy) df(l1:l2,m,n,iss) = df(l1:l2,m,n,iss) - p%ugss ! ! Calculate viscous contribution to entropy. ! if (lviscosity) call calc_viscous_heat(df,p,Hmax) ! ! Thermal conduction. ! ! if (lconductivity) call heat_conductivity(f,df,p) ! ! Entry possibility for "personal" entries. ! if (lspecial) call special_calc_energy(f,df,p) ! call calc_diagnostics_energy(f,p) ! ! ``cs2/dx^2'' for timestep ! if (lhydro.and.ldensity.and.lfirst.and.ldt) then advec_cs2=p%cs2*dxyz_2 if (headtt.or.ldebug) print*,'denergy_dt: max(advec_cs2) =',maxval(advec_cs2) endif endsubroutine denergy_dt !*********************************************************************** subroutine calc_diagnostics_energy(f,p) real, dimension (mx,my,mz,mfarray) :: f type(pencil_case) :: p ! ! Calculate entropy related diagnostics. ! call keep_compiler_quiet(f) if (ldiagnos) then call max_mn_name(p%TT,idiag_TTmax) if (idiag_TTmin/=0) call max_mn_name(-p%TT,idiag_TTmin,lneg=.true.) call sum_mn_name(p%TT,idiag_TTm) if (idiag_dtc/=0) & call max_mn_name(sqrt(advec_cs2)/cdt,idiag_dtc,l_dt=.true.) if (idiag_ethm/=0) call sum_mn_name(p%rho*p%ee,idiag_ethm) if (idiag_ethtot/=0) call integrate_mn_name(p%rho*p%ee,idiag_ethtot) if (idiag_ethdivum/=0) & call sum_mn_name(p%rho*p%ee*p%divu,idiag_ethdivum) call sum_mn_name(p%ss,idiag_ssm) endif ! ! 1D averages. Happens at every it1d timesteps, NOT at every it1 ! idiag_fradz is done in the calc_heatcond routine ! if (l1davgfirst) then if (idiag_fconvz) call xysum_mn_name_z(p%rho*p%uu(:,3)*p%TT,idiag_fconvz) call xysum_mn_name_z(p%ss,idiag_ssmz) call xzsum_mn_name_y(p%ss,idiag_ssmy) call yzsum_mn_name_x(p%ss,idiag_ssmx) call xysum_mn_name_z(p%TT,idiag_TTmz) endif ! endsubroutine calc_diagnostics_energy !*********************************************************************** subroutine energy_after_boundary(f) ! ! dummy routine ! real, dimension (mx,my,mz,mfarray) :: f intent(in) :: f ! call keep_compiler_quiet(f) ! if (lenergy_slope_limited) & call fatal_error('energy_after_boundary', & 'Slope-limited diffusion not implemented') endsubroutine energy_after_boundary !*********************************************************************** subroutine fill_farray_pressure(f) ! ! Fill f array with the pressure, to be able to calculate pressure gradient ! directly from the pressure. ! ! 18-feb-10/anders: dummy ! real, dimension (mx,my,mz,mfarray) :: f ! call keep_compiler_quiet(f) ! endsubroutine fill_farray_pressure !*********************************************************************** subroutine impose_energy_floor(f) ! ! Dummy subroutine; may not be necessary for lnTT ! real, dimension(mx,my,mz,mfarray) :: f ! call keep_compiler_quiet(f) ! endsubroutine impose_energy_floor !*********************************************************************** subroutine split_update_energy(f) ! ! Dummy subroutine ! real, dimension(mx,my,mz,mfarray), intent(inout) :: f ! call keep_compiler_quiet(f) ! endsubroutine split_update_energy !*********************************************************************** subroutine expand_shands_energy ! ! Presently dummy, for possible use ! endsubroutine expand_shands_energy !*********************************************************************** subroutine dynamical_thermal_diffusion(uc) ! ! Dummy subroutine ! real, intent(in) :: uc ! call keep_compiler_quiet(uc) ! endsubroutine dynamical_thermal_diffusion !*********************************************************************** subroutine get_slices_energy(f,slices) ! real, dimension (mx,my,mz,mfarray) :: f type (slice_data) :: slices ! call keep_compiler_quiet(f) call keep_compiler_quiet(slices%ready) ! endsubroutine get_slices_energy !*********************************************************************** subroutine rprint_energy(lreset,lwrite) ! ! Reads and registers print parameters relevant to entropy. ! ! 1-jun-02/axel: adapted from magnetic fields ! use Diagnostics use FArrayManager, only: farray_index_append ! integer :: iname,inamez,inamey,inamex,irz 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_dtc=0; idiag_ethm=0; idiag_ethdivum=0; idiag_ssm=0 idiag_ugradpm=0; idiag_ethtot=0; idiag_ssmphi=0 idiag_yHmax=0; idiag_yHm=0; idiag_TTmax=0; idiag_TTmin=0; idiag_TTm=0 idiag_fconvz=0; idiag_dcoolz=0; idiag_fradz=0; idiag_fturbz=0 idiag_ssmz=0; idiag_ssmy=0; idiag_ssmx=0; idiag_TTmz=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),'dtc',idiag_dtc) call parse_name(iname,cname(iname),cform(iname),'ethtot',idiag_ethtot) call parse_name(iname,cname(iname),cform(iname),'ethdivum',idiag_ethdivum) call parse_name(iname,cname(iname),cform(iname),'ethm',idiag_ethm) call parse_name(iname,cname(iname),cform(iname),'ssm',idiag_ssm) call parse_name(iname,cname(iname),cform(iname),'ugradpm',idiag_ugradpm) call parse_name(iname,cname(iname),cform(iname),'yHm',idiag_yHm) call parse_name(iname,cname(iname),cform(iname),'yHmax',idiag_yHmax) call parse_name(iname,cname(iname),cform(iname),'TTm',idiag_TTm) call parse_name(iname,cname(iname),cform(iname),'TTmax',idiag_TTmax) call parse_name(iname,cname(iname),cform(iname),'TTmin',idiag_TTmin) enddo ! ! Check for those quantities for which we want xy-averages. ! do inamez=1,nnamez call parse_name(inamez,cnamez(inamez),cformz(inamez),'fturbz',idiag_fturbz) call parse_name(inamez,cnamez(inamez),cformz(inamez),'fconvz',idiag_fconvz) call parse_name(inamez,cnamez(inamez),cformz(inamez),'dcoolz',idiag_dcoolz) call parse_name(inamez,cnamez(inamez),cformz(inamez),'fradz',idiag_fradz) call parse_name(inamez,cnamez(inamez),cformz(inamez),'ssmz',idiag_ssmz) call parse_name(inamez,cnamez(inamez),cformz(inamez),'TTmz',idiag_TTmz) enddo ! ! Check for those quantities for which we want xz-averages. ! do inamey=1,nnamey call parse_name(inamey,cnamey(inamey),cformy(inamey),'ssmy',idiag_ssmy) enddo ! ! Check for those quantities for which we want yz-averages. ! do inamex=1,nnamex call parse_name(inamex,cnamex(inamex),cformx(inamex),'ssmx',idiag_ssmx) enddo ! ! Check for those quantities for which we want phi-averages. ! do irz=1,nnamerz call parse_name(irz,cnamerz(irz),cformrz(irz),'ssmphi',idiag_ssmphi) enddo ! ! Write column where which entropy variable is stored. ! if (lwr) then call farray_index_append('iyH',iyH) call farray_index_append('ilnTT',ilnTT) call farray_index_append('iTT',iTT) endif ! endsubroutine rprint_energy !*********************************************************************** subroutine energy_after_timestep(f,df,dtsub) ! real, dimension(mx,my,mz,mfarray) :: f real, dimension(mx,my,mz,mvar) :: df real :: dtsub ! call keep_compiler_quiet(f,df) call keep_compiler_quiet(dtsub) ! endsubroutine energy_after_timestep !*********************************************************************** subroutine update_char_vel_energy(f) ! ! TB implemented. ! ! 25-sep-15/MR+joern: coded ! real, dimension (mx,my,mz,mfarray), intent(in) :: f call keep_compiler_quiet(f) call warning('update_char_vel_energy', & 'characteristic velocity not yet implemented for entropy_onefluid') endsubroutine update_char_vel_energy !*********************************************************************** endmodule Energy