time dependent chemical parameter, and sl 0.21

code/m_chemistry.f90

@@ -6,6 +6,10 @@ implicit none
 
 real, private :: coef(10)
 real, private :: lambda_onestep
+real, private :: lambda1_twostep
+real, private :: lambda2_twostep
+real, private :: beta1_twostep
+real, private :: hrp_twostep
 
 contains
 
@@ -13,12 +17,44 @@ subroutine init_chemistry
     if ( reaction_type == "onestep" ) then
         lambda_onestep = pre * exp ( - beta / hrp )
     else if ( reaction_type == "twostep" ) then
+        lambda1_twostep = lambda1
+        lambda2_twostep = lambda2
+        beta1_twostep = beta1
+        hrp_twostep = hrp
     else
         stop
     end if
 
 end subroutine init_chemistry
 
+
+subroutine update_chemistry (t)
+
+    real :: t
+
+    real :: factor
+
+    real :: relax_duration = 60.
+
+    if ( reaction_type == "onestep" ) then
+        lambda_onestep = pre * exp ( - beta / hrp )
+    else if ( reaction_type == "twostep" ) then
+        if (t < relax_duration) then
+            factor = (relax_duration + t) / relax_duration / 2.
+        else
+            factor = 1.
+        end if
+
+        lambda1_twostep = factor * lambda1
+        lambda2_twostep = factor * lambda2
+
+    else
+        stop
+    end if
+
+end subroutine update_chemistry
+
+
 real function rate_1step (yr, theta)
 
     real, intent(in)  :: yr
@@ -76,8 +112,8 @@ real function rate1_2step (ya, yx, theta)
     if(theta.lt.0.) t_reduce=0.
     if(theta.gt.1.) t_reduce=1.
 
-    rate1_2step = lambda1 * y1 * y2 * & 
-                  exp (-(beta1*(1. - t_reduce))/(1. - hrp*(1. - t_reduce)))
+    rate1_2step = lambda1_twostep * y1 * y2 * & 
+                  exp (-(beta1_twostep*(1. - t_reduce))/(1. - hrp_twostep*(1. - t_reduce)))
 
 end function rate1_2step
 
@@ -98,7 +134,7 @@ real function rate2_2step (yx, theta)
     if(theta.lt.0.) t_reduce=0.
     if(theta.gt.1.) t_reduce=1.
 
-    rate2_2step = lambda2 * yx * yx
+    rate2_2step = lambda2_twostep * yx * yx
 
 end function rate2_2step
 

code/m_parameters.f90

@@ -46,6 +46,10 @@ module m_parameters
 
      CHARACTER(100) :: reaction_type
 
+     real :: rate_relaxation (2,100)
+
+     integer :: table_size
+
 contains
 
      SUBROUTINE SET_CHEMISTRY
@@ -125,6 +129,14 @@ contains
         READ(itape,*) cdum,le_x
         WRITE(otape,*) cdum,le_x
 
+        ! READ(itape,*) cdum,table_size
+        ! WRITE(otape,*) cdum,table_size
+
+        ! do i = 1:table_size
+        ! READ(itape,*) cdum,rate_relaxation(1,i),rate_relaxation(2,i)
+        ! WRITE(otape,*) cdum,rate_relaxation(1,i),rate_relaxation(2,i)
+        ! end do
+
        IF (.not.read_stdin) THEN
          CLOSE(itape)
        END IF

code/sample/4pi-IC1/itape

@@ -17,9 +17,9 @@
     min_wr     0.        ! 5.0e-14
    prof_wr     1.0       
      lewis     1.0
-       lambda1   100000.
-       lambda2   63.6E+5
-       beta1   9.1
-       hrp   0.8136765
+       lambda1    6.25E+4       7.7E+4       5E+4      
+       lambda2    318E+4   445.2E+4   489.72E+4   
+       beta1   5
+       hrp   0.75
        le_a  1.
-       le_x  0.2
+       le_x  0.20

code/ysolve.f90

@@ -96,6 +96,8 @@
 
        yold = y1
 
+       CALL update_chemistry(t_now)
+
        CALL SET_BC
 
        CALL RK4(fns)
@@ -331,6 +333,8 @@
 
        CALL READ_INTRO
 
+       CALL init_chemistry
+
        ALLOCATE(   u(nx),STAT=ierr)    ; u=0.
        ALLOCATE(  y1(nx,3),STAT=ierr)  ; y1=0.
        ALLOCATE(  y2(nx,3),STAT=ierr)  ; y2=0.