Adding limit on the diffusivity

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 +++++++++++++++++++!!
 
 

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
 
 

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