diff --git a/src/Namelists.f90 b/src/Namelists.f90 index 29501c95..84a7c3a7 100644 --- a/src/Namelists.f90 +++ b/src/Namelists.f90 @@ -77,7 +77,7 @@ subroutine readNamelists(tscheme) & mpi_transp,l_adv_curl,mpi_packing namelist/phys_param/ & - & ra,rae,raxi,pr,sc,prmag,ek,gamma,epsc0,epscxi0,radratio,Bn, & + & ra,rae,rat,raxi,pr,sc,prmag,ek,gamma,epsc0,epscxi0,radratio,Bn, & & ktops,kbots,ktopv,kbotv,ktopb,kbotb,kbotxi,ktopxi, & & s_top,s_bot,impS,sCMB,xi_top,xi_bot,impXi,xiCMB, & & nVarCond,con_DecRate,con_RadRatio,con_LambdaMatch, & @@ -322,7 +322,9 @@ subroutine readNamelists(tscheme) l_AB1 =.false. l_bridge_step=.true. l_onset =.false. - l_ehd_dep=.true. + l_ehd_dep=.false. + l_ehd_die=.false. + if ( mode == 1 ) then !-- Only convection: @@ -429,6 +431,13 @@ subroutine readNamelists(tscheme) if ( ra == 0.0_cp .and. rae == 0.0_cp ) l_heat=.false. + if ( rae > 0.0_cp ) then + l_ehd_dep = .true. + if ( rat > 0.0_cp ) then + l_ehd_die = .true. + end if + end if + if ( ek < 0.0_cp ) l_non_rot= .true. if ( l_non_rot ) then l_corr=.false. @@ -929,6 +938,7 @@ subroutine writeNamelists(n_out) write(n_out,*) "&phys_param" write(n_out,'('' ra ='',ES14.6,'','')') ra write(n_out,'('' rae ='',ES14.6,'','')') rae + write(n_out,'('' rat ='',ES14.6,'','')') rat write(n_out,'('' raxi ='',ES14.6,'','')') raxi write(n_out,'('' pr ='',ES14.6,'','')') pr write(n_out,'('' sc ='',ES14.6,'','')') sc @@ -1373,6 +1383,8 @@ subroutine defaultNamelists !----- Namelist phys_param: ra =0.0_cp + rae =0.0_cp + rat =0.0_cp raxi =0.0_cp ek =1.0e-3_cp pr =one diff --git a/src/get_nl.f90 b/src/get_nl.f90 index a5ac40e3..ae522878 100644 --- a/src/get_nl.f90 +++ b/src/get_nl.f90 @@ -23,14 +23,14 @@ module grid_space_arrays_mod use radial_functions, only: or2, orho1, beta, otemp1, visc, r, or3, & & lambda, or4, or1 use physical_parameters, only: radratio, LFfac, n_r_LCR, prec_angle, ViscHeatFac, & - & oek, po, dilution_fac, ra, rae, gamma, opr, OhmLossFac, & + & oek, po, dilution_fac, ra, rae, rat, gamma, opr, OhmLossFac, & & epsPhase, phaseDiffFac, penaltyFac, tmelt use horizontal_data, only: sinTheta, cosTheta, phi, O_sin_theta_E2, & & cosn_theta_E2, O_sin_theta 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_chemical_conv, l_precession, l_centrifuge, l_phase_field, l_ehd_dep, l_ehd_die implicit none @@ -456,6 +456,10 @@ subroutine get_nl(this, time, nR, nBc, lRmsCalc) end if ! Viscous heating and Ohmic losses ? + if ( l_ehd_die ) then + this%heatTerms(:,nPhi)= & + & opr * rae/rat * radratio**2/(1.0D0-radratio)**4 * or4(nR) + end if end do !$omp end parallel diff --git a/src/get_td.f90 b/src/get_td.f90 index 85808615..d4a53242 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_chemical_conv, l_single_matrix, l_double_curl, l_ehd_die 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) then + if ( l_anel .or. l_ehd_die ) 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 ) this%heatTermsLM(lm)=zero + if ( l_anel .or. l_ehd_die ) 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 ) then + if ( l_anel .or. l_ehd_die ) 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 ) then + if ( l_anel .or. l_ehd_die ) 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 3946658b..fc295ae7 100644 --- a/src/logic.f90 +++ b/src/logic.f90 @@ -73,6 +73,7 @@ module logic logical :: l_LCR ! Switch for zero electrical conductivity beyond r_LCR 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_PressGraph ! Store pressure in graphic files diff --git a/src/phys_param.f90 b/src/phys_param.f90 index e2142b41..94b2e0fc 100644 --- a/src/phys_param.f90 +++ b/src/phys_param.f90 @@ -38,6 +38,7 @@ module physical_parameters real(cp) :: radratio ! aspect ratio real(cp) :: ra ! Rayleigh number real(cp) :: rae ! Electric Rayleigh number + real(cp) :: rat ! Thermal heating from dielectric heat real(cp) :: raxi ! Chemical composition-based Rayleigh number real(cp) :: sc ! Schmidt number (i.e. chemical Prandtl number) real(cp) :: ek ! Ekman number diff --git a/src/rIter.f90 b/src/rIter.f90 index cb7f2452..92c8a49d 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_onset, l_DTrMagSpec, l_ehd_dep, l_ehd_die 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 ) call scal_to_SH(this%gsa%heatTerms, this%nl_lm%heatTermsLM, & + if ( l_anel .or. l_ehd_die ) call scal_to_SH(this%gsa%heatTerms, this%nl_lm%heatTermsLM, & & l_R(nR)) end if if ( l_chemical_conv ) then