Fix until compiled.

geodynamics · Jun 25, 2024 · 09b7b67 · 09b7b67
1 parent 59514cb
commit 09b7b67
Show file tree

Hide file tree

Showing 4 changed files with 12 additions and 10 deletions.
diff --git a/src/Physics/Controls.F90 b/src/Physics/Controls.F90
@@ -79,7 +79,7 @@ Module Controls
     Namelist /Physical_Controls_Namelist/ magnetism, nonlinear, rotation, lorentz_forces, &
                 & viscous_heating, ohmic_heating, advect_reference_state, benchmark_mode, &
                 & benchmark_integration_interval, benchmark_report_interval, &
-                & momentum_advection, inertia, n_active_scalars, n_passive_scalars,
+                & momentum_advection, inertia, n_active_scalars, n_passive_scalars, &
                 & pseudo_incompressible
 
     !///////////////////////////////////////////////////////////////////////////

diff --git a/src/Physics/PDE_Coefficients.F90 b/src/Physics/PDE_Coefficients.F90
@@ -1303,7 +1303,7 @@ Subroutine Get_Custom_Reference()
         Enddo
         ref%exp_entropy = exp(ref%entropy/pressure_specific_heat)
         ref%dsdr_over_cp = ref%dsdr/pressure_specific_heat
-        Call log_deriv(ref%dsdr_over_cp(:), ref%d2s_over_cp(i), no_log=.true.)
+        Call log_deriv(ref%dsdr_over_cp(:), ref%d2s_over_cp(:), no_log=.true.)
 
     End Subroutine Get_Custom_Reference
 
@@ -2067,7 +2067,8 @@ Subroutine Compute_Diffusion_Coefs()
 
 
         If (pseudo_incompressible) Then
-            W_Diffusion_Coefs_0 = W_Diffusion_Coefs_0 - nu*(4.0d0/3.0d0)*ref%dsdr_over_cp*(dnu - ref%dlnrho -ref%dsdr_over_cp - 2.0d0/radius)
+            W_Diffusion_Coefs_0 = W_Diffusion_Coefs_0 - nu*(4.0d0/3.0d0)*ref%dsdr_over_cp*(dlnu - ref%dlnrho -ref%dsdr_over_cp)
+            W_Diffusion_Coefs_0 = W_Diffusion_Coefs_0 + nu*(8.0d0/3.0d0)*ref%dsdr_over_cp/radius
             W_Diffusion_Coefs_0  = W_Diffusion_Coefs_0  - nu*(4.0d0/3.0d0)*ref%d2s_over_cp
             W_Diffusion_Coefs_1  = W_Diffusion_Coefs_1  - nu*(7.0d0/3.0d0)*ref%dsdr_over_cp
             DW_Diffusion_Coefs_0 = DW_Diffusion_Coefs_0 + nu*ref%dsdr_over_cp/3.0d0

diff --git a/src/Physics/Sphere_Hybrid_Space.F90 b/src/Physics/Sphere_Hybrid_Space.F90
@@ -49,13 +49,14 @@ Subroutine Hybrid_Init()
         Allocate(drho_term(my_r%min:my_r%max))
         r1 = my_r%min
         r2 = my_r%max
-        over_rhor(r1:r2) = one_over_r(r1:r2)/ref%density(r1:r2)
-        over_rhorsq(r1:r2) = OneOverRSquared(r1:r2)/ref%density(r1:r2)
-        drho_term(r1:r2) = ref%dlnrho(r1:r2)+one_over_r(r1:r2)
         If (pseudo_incompressible) Then
-            over_rhosq(r1:r2) = over_rhosq(r1:r2)/ ref%exp_entropy(r1:r2)
-            over_rhor(r1:r2)  = over_rhor(r1:r2) / ref%exp_entropy(r1:r2)
-            drho_term(r1:r2)  = drho_term(r1:r2) + ref%dsdr_over_cp(r1:r2)
+            over_rhor(r1:r2)  = one_over_r(r1:r2)/(ref%density(r1:r2)*ref%exp_entropy(r1:r2))
+            over_rhorsq(r1:r2)= OneOverRSquared(r1:r2)/(ref%density(r1:r2)*ref%exp_entropy(r1:r2))
+            drho_term(r1:r2)  = ref%dlnrho(r1:r2)+one_over_r(r1:r2) + ref%dsdr_over_cp(r1:r2)
+        Else
+            over_rhor(r1:r2)  = one_over_r(r1:r2)/ref%density(r1:r2)
+            over_rhorsq(r1:r2)= OneOverRSquared(r1:r2)/ref%density(r1:r2)
+            drho_term(r1:r2)  = ref%dlnrho(r1:r2)+one_over_r(r1:r2)
         Endif
 
     End Subroutine Hybrid_Init

diff --git a/src/Physics/Sphere_Linear_Terms.F90 b/src/Physics/Sphere_Linear_Terms.F90
@@ -321,7 +321,7 @@ Subroutine Load_Linear_Coefficients()
                 ! amp = -exp(s/c_P) rho/(c_P * H_Laplacian) (ds/dr)    Pseudo-incompressible
                 If (pseudo_incompressible) Then
                     amp = -ref%exp_entropy * ref%density * ref%dsdr_over_cp / H_Laplacian
-                    add_implicit_term(weq,pvar, 0, amp, lp)
+                    Call add_implicit_term(weq,pvar, 0, amp, lp)
                 Endif