From 59da189bb58b5b81c9ec0cebd518d924a4f0ea50 Mon Sep 17 00:00:00 2001
From: Pedro Costa <p.simoes.costa@gmail.com>
Date: Mon, 10 Feb 2025 00:06:26 +0100
Subject: [PATCH] Add scalar transport temporal stability.

---
 src/chkdt.f90 | 6 +++---
 src/main.f90  | 6 +++---
 src/param.f90 | 4 ++++
 3 files changed, 10 insertions(+), 6 deletions(-)

diff --git a/src/chkdt.f90 b/src/chkdt.f90
index 4cf70966..0efdb369 100644
--- a/src/chkdt.f90
+++ b/src/chkdt.f90
@@ -12,7 +12,7 @@ module mod_chkdt
   private
   public chkdt
   contains
-  subroutine chkdt(n,dl,dzci,dzfi,visc,u,v,w,dtmax)
+  subroutine chkdt(n,dl,dzci,dzfi,visc,alpha,u,v,w,dtmax)
     !
     ! computes maximum allowed time step
     !
@@ -20,7 +20,7 @@ subroutine chkdt(n,dl,dzci,dzfi,visc,u,v,w,dtmax)
     integer , intent(in), dimension(3) :: n
     real(rp), intent(in), dimension(3) :: dl
     real(rp), intent(in), dimension(0:) :: dzci,dzfi
-    real(rp), intent(in) :: visc
+    real(rp), intent(in) :: visc,alpha
     real(rp), intent(in), dimension(0:,0:,0:) :: u,v,w
     real(rp), intent(out) :: dtmax
     real(rp) :: dxi,dyi,dzi
@@ -72,7 +72,7 @@ subroutine chkdt(n,dl,dzci,dzfi,visc,u,v,w,dtmax)
 #if defined(_IMPDIFF) && !defined(_IMPDIFF_1D)
     dtmax = sqrt(3.)/dti
 #else
-    dtmax = min(1.65/12./visc*dlmin**2,sqrt(3.)/dti)
+    dtmax = min(1.65/12./max(visc,alpha)*dlmin**2,sqrt(3.)/dti)
 #endif
   end subroutine chkdt
 end module mod_chkdt
diff --git a/src/main.f90 b/src/main.f90
index 47d2aebe..32670735 100644
--- a/src/main.f90
+++ b/src/main.f90
@@ -51,7 +51,7 @@ program cans
   use mod_param          , only: ng,l,dl,dli, &
                                  gtype,gr, &
                                  cfl,dtmax,dt_f, &
-                                 visc, &
+                                 visc,alpha_max, &
                                  inivel,is_wallturb, &
                                  nstep,time_max,tw_max,stop_type, &
                                  restart,is_overwrite_save,nsaves_max, &
@@ -375,7 +375,7 @@ program cans
     include 'out3d.h90'
   end if
   !
-  call chkdt(n,dl,dzci,dzfi,visc,u,v,w,dt_cfl)
+  call chkdt(n,dl,dzci,dzfi,visc,alpha_max,u,v,w,dt_cfl)
   dt = merge(dt_f,min(cfl*dt_cfl,dtmax),dt_f > 0.)
   if(myid == 0) print*, 'dt_cfl = ', dt_cfl, 'dt = ', dt
   dti = 1./dt
@@ -450,7 +450,7 @@ program cans
     end if
     if(icheck > 0.and.mod(istep,max(icheck,1)) == 0) then
       if(myid == 0) print*, 'Checking stability and divergence...'
-      call chkdt(n,dl,dzci,dzfi,visc,u,v,w,dt_cfl)
+      call chkdt(n,dl,dzci,dzfi,visc,alpha_max,u,v,w,dt_cfl)
       dt = merge(dt_f,min(cfl*dt_cfl,dtmax),dt_f > 0.)
       if(myid == 0) print*, 'dt_cfl = ', dt_cfl, 'dt = ', dt
       if(dt_cfl < small) then
diff --git a/src/param.f90 b/src/param.f90
index af08a5b1..329fc36f 100644
--- a/src/param.f90
+++ b/src/param.f90
@@ -74,6 +74,7 @@ module mod_param
 real(rp), protected, allocatable, dimension(:) :: ssource
 logical , protected, allocatable, dimension(:) :: is_sforced
 real(rp), protected, allocatable, dimension(:) :: scalf
+real(rp), protected :: alphai_min,alpha_max
 !
 #if defined(_OPENACC)
 !
@@ -242,6 +243,9 @@ subroutine read_input(myid)
   else
     nscal = 0 ! negative values equivalent to nscal = 0
   end if
+  alpha_max = huge(1._rp)
+  alpha_max = minval(alphai(1:nscal))
+  alpha_max = alpha_max**(-1)
 #if defined(_BOUSSINESQ_BUOYANCY)
   if (nscal == 0) then
       if(myid == 0) print*, 'Error reading the input file: `BOUSSINESQ_BUOYANCY` requires `nscal > 0`.'