diff --git a/src/Namelists.f90 b/src/Namelists.f90 index 84a7c3a7..b8f386ff 100644 --- a/src/Namelists.f90 +++ b/src/Namelists.f90 @@ -85,9 +85,9 @@ subroutine readNamelists(tscheme) & strat,polind,DissNb,g0,g1,g2,r_cut_model,thickStrat, & & epsS,slopeStrat,rStrat,ampStrat, & & r_LCR,nVarDiff,nVarVisc,difExp,nVarEps,interior_model, & - & nVarEntropyGrad,l_isothermal,ktopp,po,prec_angle, & - & dilution_fac,stef,tmelt,phaseDiffFac,penaltyFac, & - & epsPhase,ktopphi,kbotphi,ampForce + & nVarEntropyGrad,l_isothermal,l_vol_heat,ktopp,po, & + & prec_angle,dilution_fac,stef,tmelt,phaseDiffFac, & + & penaltyFac,epsPhase,ktopphi,kbotphi,ampForce namelist/B_external/ & & rrMP,amp_imp,expo_imp,bmax_imp,n_imp,l_imp, & @@ -990,6 +990,7 @@ subroutine writeNamelists(n_out) write(n_out,'('' radratio ='',ES14.6,'','')') radratio write(n_out,'('' l_isothermal ='',l3,'','')') l_isothermal + write(n_out,'('' l_vol_heat ='',l3,'','')') l_vol_heat write(n_out,'('' phaseDiffFac ='',ES14.6,'','')') phaseDiffFac write(n_out,'('' epsPhase ='',ES14.6,'','')') epsPhase write(n_out,'('' penaltyFac ='',ES14.6,'','')') penaltyFac @@ -1361,6 +1362,7 @@ subroutine defaultNamelists l_non_rot =.false. ! No Coriolis force ! l_anel =.false. ! Anelastic stuff ! l_isothermal =.false. ! Isothermal = 0 Grünesein ! + l_vol_heat =.false. ! No volemtric heating interior_model="None" ! Name of the interior model !---- Run time and number of threads: diff --git a/src/get_nl.f90 b/src/get_nl.f90 index ae522878..d3b030c7 100644 --- a/src/get_nl.f90 +++ b/src/get_nl.f90 @@ -30,7 +30,7 @@ module grid_space_arrays_mod use parallel_mod, only: get_openmp_blocks use constants, only: two, third, one use logic, only: l_conv_nl, l_heat_nl, l_mag_nl, l_anel, l_mag_LF, l_adv_curl, & - & l_chemical_conv, l_precession, l_centrifuge, l_phase_field, l_ehd_dep, l_ehd_die + & l_chemical_conv, l_precession, l_centrifuge, l_phase_field, l_ehd_dep, l_ehd_die, l_vol_heat implicit none @@ -460,6 +460,10 @@ subroutine get_nl(this, time, nR, nBc, lRmsCalc) this%heatTerms(:,nPhi)= & & opr * rae/rat * radratio**2/(1.0D0-radratio)**4 * or4(nR) end if + + if ( l_vol_heat ) then + this%heatTerms(:,nPhi)= one + end if end do !$omp end parallel diff --git a/src/get_td.f90 b/src/get_td.f90 index d4a53242..955b1638 100644 --- a/src/get_td.f90 +++ b/src/get_td.f90 @@ -11,7 +11,7 @@ module nonlinear_lm_mod use truncation, only: lm_max, lm_maxMag, m_min use logic, only : l_anel, l_conv_nl, l_corr, l_heat, l_anelastic_liquid, & & l_mag_nl, l_mag_kin, l_mag_LF, l_conv, l_mag, & - & l_chemical_conv, l_single_matrix, l_double_curl, l_ehd_die + & l_chemical_conv, l_single_matrix, l_double_curl, l_ehd_die, l_vol_heat use radial_functions, only: r, or2, or1, beta, epscProf, or4, temp0, orho1, l_R use physical_parameters, only: CorFac, epsc, n_r_LCR, epscXi use blocking, only: lm2l, lm2lmA, lm2lmS @@ -60,7 +60,7 @@ subroutine initialize(this,lmP_max) allocate( this%VxBrLM(lmP_max), this%VxBtLM(lmP_max), this%VxBpLM(lmP_max)) bytes_allocated = bytes_allocated + 6*lmP_max*SIZEOF_DEF_COMPLEX - if ( l_anel .or. l_ehd_die ) then + if ( l_anel .or. l_ehd_die .or. l_vol_heat ) then allocate( this%heatTermsLM(lmP_max) ) bytes_allocated = bytes_allocated+lmP_max*SIZEOF_DEF_COMPLEX end if @@ -120,7 +120,7 @@ subroutine set_zero(this) this%VStLM(lm)=zero this%VSpLM(lm)=zero end if - if ( l_anel .or. l_ehd_die ) this%heatTermsLM(lm)=zero + if ( l_anel .or. l_ehd_die .or. l_vol_heat ) this%heatTermsLM(lm)=zero if ( l_chemical_conv ) then this%VXitLM(lm)=zero this%VXipLM(lm)=zero @@ -478,7 +478,7 @@ subroutine get_dsdt(this, nR, nBc, dsdt, dVSrLM) if ( nBc == 0 ) then dsdt(1)=epsc*epscProf(nR)!+opr/epsS*divKtemp0(nR) - if ( l_anel .or. l_ehd_die ) then + if ( l_anel .or. l_ehd_die .or. l_vol_heat ) then if ( l_anelastic_liquid ) then dsdt(1)=dsdt(1)+temp0(nR)*this%heatTermsLM(1) else @@ -486,7 +486,7 @@ subroutine get_dsdt(this, nR, nBc, dsdt, dVSrLM) end if end if - if ( l_anel .or. l_ehd_die ) then + if ( l_anel .or. l_ehd_die .or. l_vol_heat ) then if ( l_anelastic_liquid ) then !$omp parallel do do lm=lm_min,lm_max diff --git a/src/logic.f90 b/src/logic.f90 index fc295ae7..c571c6d4 100644 --- a/src/logic.f90 +++ b/src/logic.f90 @@ -74,6 +74,7 @@ module logic logical :: lVerbose ! Switch for detailed information about run progress logical :: l_ehd_dep ! Switch for dilectrophoretic force logical :: l_ehd_die ! Switch for dielectric heating + logical :: l_vol_heat ! Switch for volumetric heating logical :: l_PressGraph ! Store pressure in graphic files diff --git a/src/preCalculations.f90 b/src/preCalculations.f90 index c7c53ac1..78ce9058 100644 --- a/src/preCalculations.f90 +++ b/src/preCalculations.f90 @@ -416,7 +416,7 @@ subroutine preCalc(tscheme) if ( ktops == 1 .and. kbots == 1 .and. ( tops(0,0) == 0.0_cp .and. bots(0,0) == 0.0_cp ) ) then ! Fixed entropy tops(0,0)=0.0_cp - bots(0,0)=sq4pi + !bots(0,0)=sq4pi else if ( ktops == 3 .and. kbots == 3 ) then ! Fixed temperature contrast diff --git a/src/rIter.f90 b/src/rIter.f90 index 92c8a49d..2e14aa71 100644 --- a/src/rIter.f90 +++ b/src/rIter.f90 @@ -21,7 +21,7 @@ module rIter_mod & l_precession, l_centrifuge, l_adv_curl, & & l_double_curl, l_parallel_solve, l_single_matrix,& & l_temperature_diff, l_RMS, l_phase_field, & - & l_onset, l_DTrMagSpec, l_ehd_dep, l_ehd_die + & l_onset, l_DTrMagSpec, l_ehd_dep, l_ehd_die, l_vol_heat use radial_data, only: n_r_cmb, n_r_icb, nRstart, nRstop, nRstartMag, & & nRstopMag use radial_functions, only: or2, orho1, l_R @@ -683,7 +683,7 @@ subroutine transform_to_lm_space(this, nR, lRmsCalc, dVSrLM, dVXirLM, dphidt) call spat_to_qst(this%gsa%VSr, this%gsa%VSt, this%gsa%VSp, & & dVSrLM, this%nl_lm%VStLM, this%nl_lm%VSpLM, l_R(nR)) - if ( l_anel .or. l_ehd_die ) call scal_to_SH(this%gsa%heatTerms, this%nl_lm%heatTermsLM, & + if ( l_anel .or. l_ehd_die .or. l_vol_heat ) call scal_to_SH(this%gsa%heatTerms, this%nl_lm%heatTermsLM, & & l_R(nR)) end if if ( l_chemical_conv ) then