diff --git a/src/Namelists.f90 b/src/Namelists.f90 index 8dd3a898..297f21f8 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,raxi,pr,sc,prmag,ek,epsc0,epscxi0,radratio,Bn, & + & ra,rae,raxi,pr,sc,prmag,ek,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, & @@ -323,6 +323,7 @@ subroutine readNamelists(tscheme) l_AB1 =.false. l_bridge_step=.true. l_onset =.false. + l_ehd_dep=.true. if ( mode == 1 ) then !-- Only convection: @@ -427,7 +428,7 @@ subroutine readNamelists(tscheme) l_chemical_conv = .false. end if - if ( ra == 0.0_cp ) l_heat=.false. + if ( ra == 0.0_cp .and. rae == 0.0_cp ) l_heat=.false. if ( ek < 0.0_cp ) l_non_rot= .true. if ( l_non_rot ) then @@ -928,6 +929,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,'('' raxi ='',ES14.6,'','')') raxi write(n_out,'('' pr ='',ES14.6,'','')') pr write(n_out,'('' sc ='',ES14.6,'','')') sc diff --git a/src/get_nl.f90 b/src/get_nl.f90 index d0a156ce..9db315f4 100644 --- a/src/get_nl.f90 +++ b/src/get_nl.f90 @@ -22,15 +22,15 @@ module grid_space_arrays_mod use truncation, only: n_phi_max, nlat_padded use radial_functions, only: or2, orho1, beta, otemp1, visc, r, or3, & & lambda, or4, or1 - use physical_parameters, only: LFfac, n_r_LCR, prec_angle, ViscHeatFac, & - & oek, po, dilution_fac, ra, opr, OhmLossFac, & + use physical_parameters, only: radratio, LFfac, n_r_LCR, prec_angle, ViscHeatFac, & + & oek, po, dilution_fac, ra, rae, 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_chemical_conv, l_precession, l_centrifuge, l_phase_field, l_ehd_dep implicit none @@ -46,6 +46,7 @@ module grid_space_arrays_mod real(cp), allocatable :: VSr(:,:), VSt(:,:), VSp(:,:) real(cp), allocatable :: VXir(:,:), VXit(:,:), VXip(:,:) real(cp), allocatable :: heatTerms(:,:), phiTerms(:,:) + real(cp), allocatable :: DEPFr(:,:) !----- Fields calculated from these help arrays by legtf: real(cp), allocatable :: vrc(:,:), vtc(:,:), vpc(:,:) @@ -182,6 +183,12 @@ subroutine initialize(this) bytes_allocated=bytes_allocated+2*n_phi_max*nlat_padded*SIZEOF_DEF_REAL end if + if ( l_ehd_dep ) then + allocate( this%DEPFr(nlat_padded,n_phi_max) ) + this%DEPFr(:,:)=0.0_cp + bytes_allocated=bytes_allocated + 1*n_phi_max*nlat_padded*SIZEOF_DEF_REAL + end if + end subroutine initialize !---------------------------------------------------------------------------- subroutine finalize(this) @@ -195,6 +202,7 @@ subroutine finalize(this) deallocate( this%VxBr, this%VxBt, this%VxBp, this%VSr, this%VSt, this%VSp ) if ( l_chemical_conv ) deallocate( this%VXir, this%VXit, this%VXip ) if ( l_precession ) deallocate( this%PCr, this%PCt, this%PCp ) + if ( l_ehd_dep ) deallocate( this%DEPFr ) if ( l_centrifuge ) deallocate( this%CAr, this%CAt ) if ( l_adv_curl ) deallocate( this%cvtc, this%cvpc ) if ( l_phase_field ) deallocate( this%phic, this%phiTerms ) @@ -257,6 +265,12 @@ subroutine get_nl(this, time, nR, nBc, lRmsCalc) & this%cbtc(:,nPhi)*this%brc(:,nPhi) ) end if ! Lorentz force required ? + if ( l_ehd_dep .and. (nBc == 0 .or. lRmsCalc) .and. nR>n_r_LCR ) then + !------ Get the dielectrophoretic force: + !---- r**2* RaE * eta**2 / (1-eta)**4 * Sc / r**5 + this%DEPFr(:,nPhi)= rae * opr * radratio**2/(1.0D0-radratio)**4 * this%sc(:,nPhi) * or3(nR) + end if ! DEP force required ? + if ( l_conv_nl .and. (nBc == 0 .or. lRmsCalc) ) then if ( l_adv_curl ) then ! Advection is \curl{u} \times u diff --git a/src/logic.f90 b/src/logic.f90 index 0458213a..3946658b 100644 --- a/src/logic.f90 +++ b/src/logic.f90 @@ -72,6 +72,7 @@ module logic logical :: l_perpPar ! Switch for calculation of of kinetic energy perpendicular+parallel to the rotation axis 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_PressGraph ! Store pressure in graphic files diff --git a/src/phys_param.f90 b/src/phys_param.f90 index 88f3694f..1d3452ab 100644 --- a/src/phys_param.f90 +++ b/src/phys_param.f90 @@ -37,6 +37,7 @@ module physical_parameters !-- Dimensionless parameters: real(cp) :: radratio ! aspect ratio real(cp) :: ra ! Rayleigh number + real(cp) :: rae ! Electric Rayleigh number 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 576c3f47..cb7f2452 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_onset, l_DTrMagSpec, l_ehd_dep use radial_data, only: n_r_cmb, n_r_icb, nRstart, nRstop, nRstartMag, & & nRstopMag use radial_functions, only: or2, orho1, l_R @@ -667,6 +667,10 @@ subroutine transform_to_lm_space(this, nR, lRmsCalc, dVSrLM, dVXirLM, dphidt) this%gsa%Advr(:, nPhi)=this%gsa%Advr(:,nPhi) + this%gsa%CAr(:,nPhi) this%gsa%Advt(:, nPhi)=this%gsa%Advt(:,nPhi) + this%gsa%CAt(:,nPhi) end if + + if ( l_ehd_dep ) then + this%gsa%Advr(:, nPhi)=this%gsa%Advr(:,nPhi) + this%gsa%DEPFr(:,nPhi) + end if end do !$omp end parallel