diff --git a/code/ysolve.f90 b/code/ysolve.f90 index 01032e8..1a784d1 100644 --- a/code/ysolve.f90 +++ b/code/ysolve.f90 @@ -107,10 +107,12 @@ pflame=0. DO i=1,nx pflame=pflame+y1(i,1)*hx - uxt(i)=y1(i,1) + uxt(i)=y1(i,3) ENDDO + CALL dfnonp(nx,hx,uxt,duxt,1,1) delf=1./MAXVAL(ABS(duxt)) + oldu=u(1) u=u+0.5*(hx*REAL(nx)*tar_lo-pflame)+0.5*(pflold-pflame) t_uf=0. @@ -231,29 +233,27 @@ SUBROUTINE write_pre REAL :: theta,yr,c,dy,maxdy=0.,del_f - REAL :: S_L=0.,wrate + REAL :: S_L=0.,wrate,wrate1,wrate2 + REAL :: ya,yx INTEGER :: i REAL, DIMENSION(1,nx) :: ux, dux DO i=1,nx yr=y1(i,1) + + ya=y1(i,1) + yx=y1(i,2) theta=y1(i,3) - c=1.-yr - ux(1,i)=yr - wrate=pre*yr*exp(-ac/(1.+bc*theta)) !wrate -! IF (c.le.c_cut) THEN -! wrate=min_wr -! IF (c.gt.c_ref) wrate= & -! ((exp((c-c_cut)*prof_wr)-minf)/(1.-minf)*(refwr-min_wr))+min_wr -! ENDIF - IF (theta.le.c_ref) THEN - wrate=min_wr - IF (theta.gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(theta-c_ref))+ & - min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref))) - ENDIF + ux(1,i)=theta - S_L=S_L+wrate*hx + wrate=rate_1step(yr, theta) + + wrate1=rate1_2step(ya, yx, theta) + + wrate2=rate2_2step(yx, theta) + + S_L=S_L+wrate1*hx ENDDO WRITE(*,'(a31,e14.8)') ' INTEGRAL( Wrate x dx ) => Sc :',S_L