Add the volumetric heating term

The volumetric heating term is enabled using the boolean l_vol_heat. This
term scales with the inverse of r^2 and is meant for the use using
buoyancy with the Rayleigh number.

Author: ygaillard

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:

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)= opr
+         end if
       end do
       !$omp end parallel
 

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,15 +478,15 @@ 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
                dsdt(1)=dsdt(1)+this%heatTermsLM(1)
             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

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
 

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