Add scalar transport temporal stability.

src/chkdt.f90

@@ -12,15 +12,15 @@ 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
     !
     implicit none
     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

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

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`.'