diff --git a/code/ysolve.f90 b/code/ysolve.f90 index be1c76d..887a1e9 100644 --- a/code/ysolve.f90 +++ b/code/ysolve.f90 @@ -159,14 +159,14 @@ INTEGER :: i,j,k,nd REAL, DIMENSION(2,nx) :: ux,dux,d2ux REAL, DIMENSION(10,nx) :: sav - + sav=0. ! refwr=pre*1.*exp(-ac/(1.+bc*c_cut)) ! minf=exp((c_ref-c_cut)*prof_wr) refwr=pre*1.*exp(-ac/(1.+bc*c_ref)) - - + + WRITE(500,*) 'VARIABLES = "X","Yr","C","U","Wrate","|DEL(Y)|","Sdr"' WRITE(500,*) ' "Sdd","Sd","(1/C)/(dC/dx)","DIV(rho*Dmu*Gra(C))"' @@ -175,7 +175,7 @@ ux(2,i)=y1(i,2) ! T IF (ux(1,i).gt.1.) ux(1,i)=1. ENDDO - + nd=2 CALL dfnonp(nx,hx,ux(:,:),dux(:,:),nd,1) CALL d2fnonp(nx,hx,ux(:,:),d2ux(:,:),nd,1) @@ -252,25 +252,41 @@ REAL :: ya,yx INTEGER :: i REAL, DIMENSION(1,nx) :: ux, dux - + + if ( reaction_type == "onestep" ) then + DO i=1,nx yr=y1(i,1) + theta=y1(i,2) + ux(1,i)=theta + + wrate=rate_1step(yr, theta) + + S_L=S_L+wrate*hx + ENDDO + + else if ( reaction_type == "twostep" ) then + + DO i=1,nx ya=y1(i,1) yx=y1(i,2) theta=y1(i,3) ux(1,i)=theta - wrate=rate_1step(yr, theta) - wrate1=rate1_2step(ya, yx, theta) wrate2=rate2_2step(yx, theta) S_L=S_L+wrate1*hx ENDDO - + + else + WRITE(*,*) 'ERROR, UNDEFINED REACTION TYPE ', reaction_type + stop + end if + WRITE(*,'(a31,e14.8)') ' INTEGRAL( Wrate x dx ) => Sc :',S_L CALL dfnonp(nx,hx,ux(1,:),dux(1,:),1,1) @@ -280,7 +296,7 @@ maxdy=MAX(ABS(dy),maxdy) ENDDO del_f=1./maxdy - + WRITE(*,'(a13,e14.8,a25,e14.8)') ' Grid size : ',hx,' / Laminar flame speed : ',u(1) WRITE(*,'(a19,e14.8,a3,f9.5,a20)') & ' Flame thickness : ',del_f,' / ',del_f/hx,' grids in the flame.'