Skip to content

Commit 92aa2cb

Browse files
AnkitBarikAnkitBarik
committed
Tidal forcing fix
1 parent b727c3e commit 92aa2cb

5 files changed

Lines changed: 67 additions & 37 deletions

File tree

‎src/nonlinear_bcs.f90‎

Lines changed: 15 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -117,7 +117,7 @@ subroutine get_b_nl_bcs(bc,br_vt_lm,br_vp_lm,b_nl_bc,aj_nl_bc)
117117
end subroutine get_b_nl_bcs
118118
!-------------------------------------------------------------------------
119119
subroutine v_rigid_boundary(nR,omega,lDeriv,vrr,vtr,vpr,cvrr,dvrdtr, &
120-
& dvrdpr,dvtdpr,dvpdpr)
120+
& dvrdpr,dvtdpr,dvpdpr,l_radial_flow_bc)
121121
!
122122
! Purpose of this subroutine is to set the velocities and their
123123
! derivatives at a fixed boundary.
@@ -129,7 +129,7 @@ subroutine v_rigid_boundary(nR,omega,lDeriv,vrr,vtr,vpr,cvrr,dvrdtr, &
129129
!-- Input of variables:
130130
integer, intent(in) :: nR ! no of radial grid point
131131
logical, intent(in) :: lDeriv ! derivatives required ?
132-
132+
logical, intent(in) :: l_radial_flow_bc ! radial flow at boundary
133133
!-- Input of boundary rotation rate
134134
real(cp), intent(in) :: omega ! boundary rotation rate
135135

@@ -153,17 +153,24 @@ subroutine v_rigid_boundary(nR,omega,lDeriv,vrr,vtr,vpr,cvrr,dvrdtr, &
153153
return
154154
end if
155155

156+
if ( .not. l_radial_flow_bc ) then
157+
vrr(:,:)=0.0_cp
158+
end if
159+
160+
vtr(:,:)=0.0_cp
161+
162+
if ( lDeriv ) then
163+
dvrdtr(:,:)=0.0_cp
164+
dvrdpr(:,:)=0.0_cp
165+
dvtdpr(:,:)=0.0_cp
166+
dvpdpr(:,:)=0.0_cp
167+
end if
168+
156169
!$omp parallel do default(shared)
157170
do nPhi=1,n_phi_max
158-
vrr(:,nPhi)=0.0_cp
159-
vtr(:,nPhi)=0.0_cp
160171
vpr(:,nPhi)=r2*rho0(nR)*sinTheta_E2(:)*omega
161172
if ( lDeriv ) then
162173
cvrr(:,nPhi) =r2*rho0(nR)*two*cosTheta(:)*omega
163-
dvrdtr(:,nPhi)=0.0_cp
164-
dvrdpr(:,nPhi)=0.0_cp
165-
dvtdpr(:,nPhi)=0.0_cp
166-
dvpdpr(:,nPhi)=0.0_cp
167174
end if
168175
end do
169176
!$omp end parallel do

‎src/rIter.f90‎

Lines changed: 27 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -568,18 +568,34 @@ subroutine transform_to_grid_space(this, nR, nBc, lViscBcCalc, lRmsCalc, &
568568
& this%gsa%dvtdpc, this%gsa%dvpdpc, l_R(nR))
569569
end if
570570
else if ( nBc == 2 ) then
571-
if ( nR == n_r_cmb .and. (.not. l_radial_flow_bc) ) then
572-
call v_rigid_boundary(nR, omega_ma, lDeriv, this%gsa%vrc, &
573-
& this%gsa%vtc, this%gsa%vpc, this%gsa%cvrc, &
574-
& this%gsa%dvrdtc, this%gsa%dvrdpc, &
575-
& this%gsa%dvtdpc,this%gsa%dvpdpc)
576-
else if ( nR == n_r_icb .and. ellip_fac_icb == 0.0_cp ) then
577-
call v_rigid_boundary(nR, omega_ic, lDeriv, this%gsa%vrc, &
578-
& this%gsa%vtc, this%gsa%vpc, &
579-
& this%gsa%cvrc, this%gsa%dvrdtc, &
580-
& this%gsa%dvrdpc, this%gsa%dvtdpc, &
581-
& this%gsa%dvpdpc)
571+
if ( nR == n_r_cmb ) then
572+
if ( l_radial_flow_bc ) then
573+
call torpol_to_spat(w_Rloc(:,nR), dw_Rloc(:,nR), z_Rloc(:,nR), &
574+
& this%gsa%vrc, this%gsa%vtc, this%gsa%vpc, l_R(nR))
575+
end if
576+
call v_rigid_boundary(nR, omega_ma, lDeriv, this%gsa%vrc, &
577+
& this%gsa%vtc, this%gsa%vpc, this%gsa%cvrc, &
578+
& this%gsa%dvrdtc, this%gsa%dvrdpc, &
579+
& this%gsa%dvtdpc, this%gsa%dvpdpc, l_radial_flow_bc)
582580
end if
581+
582+
if ( nR == n_r_icb ) then
583+
if ( ellip_fac_icb == 0.0_cp ) then
584+
call torpol_to_spat(w_Rloc(:,nR), dw_Rloc(:,nR), z_Rloc(:,nR), &
585+
& this%gsa%vrc, this%gsa%vtc, this%gsa%vpc, l_R(nR))
586+
587+
call v_rigid_boundary(nR, omega_ma, lDeriv, this%gsa%vrc, &
588+
& this%gsa%vtc, this%gsa%vpc, this%gsa%cvrc, &
589+
& this%gsa%dvrdtc, this%gsa%dvrdpc, &
590+
& this%gsa%dvtdpc, this%gsa%dvpdpc, l_radial_flow_bc)
591+
else
592+
call v_rigid_boundary(nR, omega_ma, lDeriv, this%gsa%vrc, &
593+
& this%gsa%vtc, this%gsa%vpc, this%gsa%cvrc, &
594+
& this%gsa%dvrdtc, this%gsa%dvrdpc, &
595+
& this%gsa%dvtdpc, this%gsa%dvpdpc, .false.) ! A bit hacky, should ideally separate icb and cmb BC flags, TBD
596+
end if
597+
end if
598+
583599
if ( lDeriv ) then
584600
call torpol_to_spat(dw_Rloc(:,nR), ddw_Rloc(:,nR), dz_Rloc(:,nR), &
585601
& this%gsa%dvrdrc, this%gsa%dvtdrc, &

‎src/special.f90‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -47,7 +47,7 @@ module special
4747
real(cp), public :: tide_fac20 ! Pre-factor for tidal forcing for (2,0)
4848
real(cp), public :: tide_fac22p ! Pre-factor for tidal forcing for (2,2)
4949
real(cp), public :: tide_fac22n ! Pre-factor for tidal forcing for (2,-2)
50-
logical, public :: l_radial_flow_bc ! A flag that tells whether radial flow BCs are being used
50+
logical, public :: l_radial_flow_bc ! A flag that tells whether radial flow BCs are being used
5151

5252
real(cp), public :: ampForce ! Amplitude of external body force
5353

‎src/updateS.f90‎

Lines changed: 14 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -360,6 +360,7 @@ subroutine prepareS_FD(tscheme, dsdt, phi, time)
360360

361361
!-- Local variables
362362
integer :: nR, lm_start, lm_stop, lm, l, m
363+
complex(cp):: sBC
363364

364365
!-- LU factorisation of the matrix if needed
365366
if ( .not. lSmat(0) ) then
@@ -391,10 +392,12 @@ subroutine prepareS_FD(tscheme, dsdt, phi, time)
391392
l = st_map%lm2l(lm)
392393
m = st_map%lm2m(lm)
393394
if ( ktops == 1 ) then ! Fixed temperature
394-
if (amp_tide /= 0.0_cp .and. l == 2) then
395-
call get_radial_flow_bc(time, m, tops(l,m))
395+
if (amp_tide /= 0.0_cp .and. l == 2 .and. present(time)) then ! Check if this is testLM
396+
call get_radial_flow_bc(time, m, sBC)
397+
else
398+
sBC = zero
396399
end if
397-
s_ghost(lm,nR)=tops(l,m)
400+
s_ghost(lm,nR)=tops(l,m) + sBC
398401
else ! Fixed flux
399402
!TBD
400403
s_ghost(lm,nR)=s_ghost(lm,nR)+fd_fac_top(l)*tops(l,m)
@@ -1352,13 +1355,17 @@ subroutine get_radial_flow_bc(time, m, sBC)
13521355
fac = 6.0_cp / r_cmb**2 ! l(l+1)/r_o^2 for l=2
13531356

13541357
if (m==0) then
1355-
sBC = sBC + radratio/omega_tide * fac * tide_fac20 * &
1356-
& cmplx(sin(omega_tide*time), &
1357-
& -cos(omega_tide*time), cp)
1358+
sBC = radratio/omega_tide * fac * tide_fac20 * &
1359+
& cmplx(sin(omega_tide*time), &
1360+
& -cos(omega_tide*time), cp)
1361+
end if
1362+
1363+
if (m==1) then ! Obliquity tide, yet to implement
1364+
sBC = zero
13581365
end if
13591366

13601367
if (m==2) then
1361-
sBC = sBC + radratio/omega_tide * fac * &
1368+
sBC = radratio/omega_tide * fac * &
13621369
& cmplx((-tide_fac22p + tide_fac22n) * sin(omega_tide*time), & ! B * sin(ωt)
13631370
& -(tide_fac22p + tide_fac22n) * cos(omega_tide*time), cp) ! -A * cos(ωt)
13641371
end if

‎src/updateWP.f90‎

Lines changed: 10 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -419,9 +419,9 @@ subroutine updateWP(time, s, xi, w, dw, ddw, dwdt, p, dp, dpdt, tscheme, &
419419

420420
if ( l1 == 2 .and. m1 == 0 ) then
421421
rhs1(1,2*lm-1,threadid)=rhs1(1,2*lm-1,threadid) &
422-
& + tide_fac20 * cos(omega_tide*(time))
422+
& + tide_fac20 * cos(omega_tide*time)
423423
rhs1(1,2*lm,threadid) =rhs1(1,2*lm,threadid) &
424-
& + tide_fac20 * sin(omega_tide*(time))
424+
& + tide_fac20 * sin(omega_tide*time)
425425
end if
426426

427427
if ( l1 == 2 .and. m1 == 2 ) then
@@ -443,10 +443,10 @@ subroutine updateWP(time, s, xi, w, dw, ddw, dwdt, p, dp, dpdt, tscheme, &
443443
! sign in the rotating frame (Ogilvie, 2014)
444444
rhs1(1,2*lm-1,threadid)=rhs1(1,2*lm-1,threadid) &
445445
& + (tide_fac22p + tide_fac22n) &
446-
& * cos(omega_tide*(time))
446+
& * cos(omega_tide*time)
447447
rhs1(1,2*lm,threadid) =rhs1(1,2*lm,threadid) &
448448
& + (-tide_fac22p + tide_fac22n) &
449-
& * sin(omega_tide*(time))
449+
& * sin(omega_tide*time)
450450

451451
end if
452452
end if
@@ -2690,17 +2690,17 @@ subroutine get_radial_flow_bc(time, m, nBC, wbc)
26902690

26912691
if ( nBC == n_r_cmb ) then
26922692
if ( m == 0 ) then
2693-
wbc = tide_fac20 * cmplx(cos(omega_tide*(time)), &
2694-
& sin(omega_tide*(time)), cp)
2693+
wbc = tide_fac20 * cmplx(cos(omega_tide*time), &
2694+
& sin(omega_tide*time), cp)
26952695
end if
26962696
if ( m == 2 ) then
26972697
wbc = ellip_fac_cmb * &
26982698
& cmplx(cos(omegaOsz_ma1*(time+tShift_ma1)), &
26992699
& sin(omegaOsz_ma1*(time+tShift_ma1)), cp) + &
2700-
& tide_fac22p * cmplx(cos(omega_tide*(time)), &
2701-
& -sin(omega_tide*(time)), cp) + &
2702-
& tide_fac22n * cmplx(cos(omega_tide*(time)), &
2703-
& sin(omega_tide*(time)), cp)
2700+
& tide_fac22p * cmplx(cos(omega_tide*time), &
2701+
& -sin(omega_tide*time), cp) + &
2702+
& tide_fac22n * cmplx(cos(omega_tide*time), &
2703+
& sin(omega_tide*time), cp)
27042704
end if
27052705
end if
27062706

0 commit comments

Comments
 (0)