diff --git a/code/ysolve.f90 b/code/ysolve.f90 index 8d835ff..0f5d671 100644 --- a/code/ysolve.f90 +++ b/code/ysolve.f90 @@ -116,7 +116,6 @@ REAL :: c,yr,wrate,dely,sdr,sdd,sd,uu,onelw,dd INTEGER :: i,j,k,nd REAL, DIMENSION(1,nx) :: ux,dux,d2ux - REAL, DIMENSION(1,nx) :: duy REAL, DIMENSION(10,nx) :: sav sav=0. @@ -134,15 +133,11 @@ IF (ux(1,i).gt.1.) ux(1,i)=1. ENDDO - dm = (/ (diff * diffusivity_T076(1.-y1(i,1,1)), i=1,nx) /) - nd=1 CALL dfnonp(nx,hx,ux(1,:),dux(1,:),nd,1) - CALL dfnonp(nx,hx,dm,duy(1,:),nd,1) nd=1 CALL d2fnonp(nx,hx,ux(1,:),d2ux(1,:),nd,1) - DO i=1,nx yr=ux(1,i) c=1.-yr @@ -164,14 +159,14 @@ dely=ABS(dux(1,i)) sdr=wrate/dely - sdd=-(dm(i)*d2ux(1,i)+duy(1,i)*dux(1,i))/dely + sdd=-dm(i)*d2ux(1,i)/dely IF (dely.eq.0.) THEN sdr=0.; sdd=0. ENDIF sd=sdr+sdd uu=u(1,1,1) onelw=(-dux(1,i))/c - dd=-(dm(i)*d2ux(1,i)+duy(1,i)*dux(1,i)) + dd=-dm(i)*d2ux(1,i) if (c.eq.0.) onelw=0. sav(1,i)=sav(1,i)+yr @@ -419,7 +414,6 @@ REAL, INTENT(IN),DIMENSION(:,:,:) :: r1 REAL, INTENT(OUT),DIMENSION(:,:,:) :: f REAL, DIMENSION(3,nx) :: ux,dux,d2ux - REAL, DIMENSION(1,nx) :: duy INTEGER :: i,j,k REAL :: wrate,Ly,Dy @@ -432,10 +426,7 @@ ux(3,i)=r1(i,j,k) ! Y ENDDO - dm = (/ (diff * diffusivity_T076(1.-r1(i,1,1)), i=1,nx) /) - CALL dfnonp(nx,hx,ux(1:3,:),dux(1:3,:),3,1) - CALL dfnonp(nx,hx,dm,duy(1,:),1,1) CALL d2fnonp(nx,hx,ux(3:3,:),d2ux(1,:),1,1) DO i=1,nx @@ -454,7 +445,7 @@ ! -0.5*( d(u*Y)/dx + u*dY/dx + Y*du/dx ) + D*d2Y/d2x f(i,j,k)=-0.5*( dux(1,i) + ux(2,i)*dux(3,i) + & ux(3,i)*dux(2,i) ) & - + dm(i)*d2ux(1,i) + duy(1,i)*dux(3,i) - wrate + + dm(i)*d2ux(1,i) - wrate IF (i.eq.nx) THEN Ly=ux(2,nx)*dux(3,nx) ! Ly = u*dYr/dx f(nx,1,1)=-wrate