From 37c5e7821f010ae7c6b8cd72302fac29b60c13a4 Mon Sep 17 00:00:00 2001
From: nazmulazim <nazmul.azim@gmail.com>
Date: Mon, 20 Sep 2021 09:57:22 -0500
Subject: [PATCH] Adding limit on the diffusivity

---
 Moussa_et_al_network.f90                | 11 +++++++++++
 mainLoop_Moussa_test_repeatTimeLoop.f90 | 14 +++++++++++++-
 var_module.f90                          |  2 ++
 3 files changed, 26 insertions(+), 1 deletion(-)

diff --git a/Moussa_et_al_network.f90 b/Moussa_et_al_network.f90
index f154e27..55546e6 100644
--- a/Moussa_et_al_network.f90
+++ b/Moussa_et_al_network.f90
@@ -665,6 +665,17 @@ subroutine mesh_diffusive_backward(j,ntim,repeatInterval)
 			if (diffusivity(i,j) .lt. minDiffuLm) diffusivity(i,j) = minDiffuLm !!! Applying diffusivity lower limit
 		end do
 
+
+
+		!!! test of min_X_plus and min_T_plus
+		do i=1,ncomp-1
+            if (celerity(i,j)/diffusivity(i,j)*dx(i,j) .le. 0.1) diffusivity(1:ncomp,j) = celerity(i,j)/0.1*dx(i,j)
+            if (celerity(i,j)*celerity(i,j)/diffusivity(i,j)*dtini .le. 0.3) &
+                diffusivity(1:ncomp,j) = celerity(i,j)*celerity(i,j)/0.3*dtini
+            min_X_plus = min(min_X_plus,celerity(i,j)/diffusivity(i,j)*dx(i,j))
+            min_T_plus = min(min_T_plus,celerity(i,j)*celerity(i,j)/diffusivity(i,j)*dtini)
+        end do
+
 !!++++++++++++++++++++ Diffusive wave Backward sweep ends +++++++++++++++++++!!
 
 
diff --git a/mainLoop_Moussa_test_repeatTimeLoop.f90 b/mainLoop_Moussa_test_repeatTimeLoop.f90
index 1d9032a..76a746f 100644
--- a/mainLoop_Moussa_test_repeatTimeLoop.f90
+++ b/mainLoop_Moussa_test_repeatTimeLoop.f90
@@ -189,6 +189,11 @@ program mesh
     open(unit=1,file="D:\Project_Works\JTTI\Vermilion_CNT\input_Vermelion_CNT.txt",status='unknown')
     open(unit=1,file="D:\Project_Works\JTTI\Florence_NC\Model\input_file_737_final_nonInterpolatedSections_Moussa",status='unknown')
 
+    !!! testing of new diffusivity limit as shown in Moussa et.at (1996)
+    open(unit=1,file="D:\Project_Works\JTTI\Moussa_etal\Test_Case\compoundTrapoCS\input_Y_chn.txt",status='unknown')
+    open(unit=1,file="D:\Project_Works\JTTI\Florence_NC\Model\input_file_737_final_nonInterpolatedSections_Moussa",status='unknown')
+
+
     print*, 'Reading input file'
 
     ! read data
@@ -750,6 +755,11 @@ program mesh
 	repeatInterval = int(60.*60./dtini_given)
 
 
+    !!! test of min_X_plus and min_T_plus
+    min_X_plus = 999.
+    min_T_plus = 999.
+
+
     !! calculation loop starts !!
     do kkk = 1,totalTimeSteps-1, repeatInterval
 
@@ -874,7 +884,7 @@ program mesh
     !pause 7
 
 
-    print*, 'Finished timestep', kkk+repeatInterval-1, 'out of ', totalTimeSteps-1
+    print*, 'Finished timestep', kkk+repeatInterval-1, 'out of ', totalTimeSteps-1, 'min_X_plus',min_X_plus,'min_T_plus',min_T_plus
 
     end do ! end kkk loop
 
@@ -894,6 +904,8 @@ program mesh
 
     print*, t2-t1, 'sec'
 
+    print*, 'min_X_plus',min_X_plus,'min_T_plus',min_T_plus
+
 end program mesh
 
 
diff --git a/var_module.f90 b/var_module.f90
index cd4518b..6021420 100644
--- a/var_module.f90
+++ b/var_module.f90
@@ -23,4 +23,6 @@ module var_module
     real :: dxini,lastKnownDiffuDT, skk, minDiffuLm, maxDiffuLm!, tfin ntim
     integer :: applyNaturalSection  ! if 1, then attribute table will be activated, if 0, then rectangular channel will be applied
 
+    real :: min_X_plus, min_T_plus
+
 end module var_module