From 2ce464838aeb4e1c9d5a9c87aa35c1b7f0797a8b Mon Sep 17 00:00:00 2001
From: Yann Gaillard <yann@gaillard.onl>
Date: Tue, 11 Feb 2025 15:04:42 +0100
Subject: [PATCH] Centrifugal Buoyancy Term

Add the centrifugal buoyancy via the terms CAr andCAt. This is only
possible whenn Ek>0 annd the buoyancy coefficient gamma also is
above zero.
---
 src/Namelists.f90  |  7 ++++---
 src/get_nl.f90     | 11 +++++++++--
 src/phys_param.f90 |  1 +
 3 files changed, 14 insertions(+), 5 deletions(-)

diff --git a/src/Namelists.f90 b/src/Namelists.f90
index fa741969..29501c95 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,epsc0,epscxi0,radratio,Bn,      &
+      &    ra,rae,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,     &
@@ -551,7 +551,7 @@ subroutine readNamelists(tscheme)
       end if
 
       !-- If dilution factor is not zero, then centrifugal force on
-      if (dilution_fac == 0.0_cp) then
+      if (dilution_fac == 0.0_cp .and. (l_non_rot .or. gamma == 0.0_cp)) then
          l_centrifuge = .false.
       else
          l_centrifuge = .true.
@@ -574,7 +574,7 @@ subroutine readNamelists(tscheme)
       end if
 
       if ( l_centrifuge .and. .not.  &
-      &    (l_anel .and. .not. l_isothermal .and. (index(interior_model, "NONE")/=0)) )  then
+      &    (l_anel .and. .not. l_isothermal .and. (index(interior_model, "NONE")/=0)) .and. gamma == 0.0_cp )  then
          call abortRun("This case is not implemented.")
          ! centrifugal acceleration implemented for anelastic polytropic background so far
       end if
@@ -934,6 +934,7 @@ subroutine writeNamelists(n_out)
       write(n_out,'(''  sc              ='',ES14.6,'','')') sc
       write(n_out,'(''  prmag           ='',ES14.6,'','')') prmag
       write(n_out,'(''  ek              ='',ES14.6,'','')') ek
+      write(n_out,'(''  gamma           ='',ES14.6,'','')') gamma
       write(n_out,'(''  po              ='',ES14.6,'','')') po
       write(n_out,'(''  stef            ='',ES14.6,'','')') stef
       write(n_out,'(''  tmelt           ='',ES14.6,'','')') tmelt
diff --git a/src/get_nl.f90 b/src/get_nl.f90
index 30f5b990..a5ac40e3 100644
--- a/src/get_nl.f90
+++ b/src/get_nl.f90
@@ -23,7 +23,7 @@ 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, opr, OhmLossFac, &
+        &                         oek, po, dilution_fac, ra, rae, gamma, opr, OhmLossFac, &
         &                         epsPhase, phaseDiffFac, penaltyFac, tmelt
    use horizontal_data, only: sinTheta, cosTheta, phi, O_sin_theta_E2, &
        &                      cosn_theta_E2, O_sin_theta
@@ -369,7 +369,7 @@ subroutine get_nl(this, time, nR, nBc, lRmsCalc)
             &                 this%vrc(:,nPhi)*cosTheta(:))
          end if ! precession term required ?
 
-         if ( l_centrifuge .and. nBc ==0 ) then
+         if ( l_centrifuge .and. nBc == 0 .and. gamma == 0.0_cp) then
             !if ( l_anel ) then
             !   this%CAr(:,nPhi) = dilution_fac*r(nR)*sinTheta(:)**4* &
             !   &       ( -ra*opr*this%sc(:,nPhi) )
@@ -387,6 +387,13 @@ subroutine get_nl(this, time, nR, nBc, lRmsCalc)
             !end if
          end if ! centrifuge
 
+         if (l_centrifuge .and. gamma > 0) then
+            ! ---- r^2 *          (-gamma) / Ek^2 * sin(theta)^2            * T
+            this%CAr(:,nPhi) =   -gamma * oek**2 *  r(nR)**3 * sinTheta(:)**2 * this%sc(:,nPhi)
+            ! ---- sin(theta)/r * (-gamma) / Ek^2 * sin(theta) * cos(theta) * T
+            this%CAt(:,nPhi) =   -gamma * oek**2             * sinTheta(:)**2 * cosTheta(:) *  this%sc(:,nPhi)
+         end if
+
          if ( l_mag_nl ) then
 
             if ( nBc == 0 .and. nR>n_r_LCR ) then
diff --git a/src/phys_param.f90 b/src/phys_param.f90
index 88b68d8b..e2142b41 100644
--- a/src/phys_param.f90
+++ b/src/phys_param.f90
@@ -54,6 +54,7 @@ module physical_parameters
    real(cp) :: ViscHeatFac    ! Prefactor for viscous heating: :math:`Di\,Pr/Ra`
    real(cp) :: OhmLossFac     ! Prefactor for Ohmic heating: :math:`Di\,Pr/(Ra\,E\,Pm^2)`
    real(cp) :: DissNb         ! Dissipation number
+   real(cp) :: gamma          ! buoyant thermal expasion coefficient
    real(cp) :: ThExpNb        ! Thermal expansion * temperature :math:`\alpha_0 T_0`
    real(cp) :: GrunNb         ! Grüneisen paramater :math:`\Gamma=(\gamma-1)/\alpha T`
    real(cp) :: epsS           ! Deviation from the adiabat