From c5abbc2e3b7156f7b37853fbc6b529cb049ab501 Mon Sep 17 00:00:00 2001
From: Yeongdo Park <ignis@postech.ac.kr>
Date: Fri, 23 Oct 2020 19:24:01 -0700
Subject: [PATCH] sutherland diffusivity

---
 code/ysolve.f90 | 25 +++++++++++++++++++++++++
 1 file changed, 25 insertions(+)

diff --git a/code/ysolve.f90 b/code/ysolve.f90
index d79727f..3bbab83 100644
--- a/code/ysolve.f90
+++ b/code/ysolve.f90
@@ -93,6 +93,8 @@
          CALL solve_with_diffusivity (update_dm)
        ELSE IF (d_mode == 1) THEN
          CALL solve_with_diffusivity (update_dm_cd)
+       ELSE IF (d_mode == 2) THEN
+         CALL solve_with_diffusivity (update_dm_sutherland)
        ELSE
          WRITE(*,*)'ERROR, UNSUPPORTED DIFFUSIVITY ', d_mode, '. STOPPING'
          STOP
@@ -567,6 +569,15 @@
       ENDDO
      END SUBROUTINE update_dm_cd
 
+     SUBROUTINE update_dm_sutherland(r1)
+      REAL, INTENT(IN),DIMENSION(:,:) :: r1
+      INTEGER :: i
+
+      DO i=1,nx
+        dm(i) = diff * diffusivity_sutherland(r1(i,nsp)) + d_turb
+      ENDDO
+     END SUBROUTINE update_dm_sutherland
+
      SUBROUTINE fonestep(r1,f)
       REAL, INTENT(IN),DIMENSION(:,:) :: r1
       REAL, INTENT(OUT),DIMENSION(:,:) :: f
@@ -656,4 +667,18 @@
 
      END SUBROUTINE fns
 
+     REAL FUNCTION diffusivity_sutherland(c)
+         REAL, INTENT(IN) :: c
+         REAL :: theta, As, Ts, T0, T1
+         T0 = 1.0
+         T1 = (1.0 + bc)
+
+         Ts = (rvis*T1 - (T1**(3./2.))) / (T1**(3./2.) - rvis)
+         As = (T0 + Ts)
+
+         theta = (1.0 + bc * c)
+
+         diffusivity_sutherland = As * sqrt(theta) / (1. + Ts/theta)
+     END FUNCTION diffusivity_sutherland
+
     END MODULE ysolve