📄 conf2.f
字号:
* print*, 'istate after', istate
moment00(I+1)=Y(1)
moment10(I+1)=Y(2)
moment01(I+1)=Y(3)
moment11(I+1)=Y(4)
moment20(I+1)=Y(5)
moment02(I+1)=Y(6)
moment21(I+1)=Y(7)
moment12(I+1)=Y(8)
moment22(I+1)=Y(9)
moment30(I+1)=Y(10)
moment31(I+1)=Y(11)
concentration(I+1)=Y(12)
* concentration_measured(I+1)=concentration(I+1)+
* & betaconc*DRNNOF()
seed_moment10(I+1)=Y(13)
seed_moment01(I+1)=Y(14)
seed_moment11(I+1)=Y(15)
seed_moment20(I+1)=Y(16)
seed_moment02(I+1)=Y(17)
seed_moment21(I+1)=Y(18)
seed_moment12(I+1)=Y(19)
relsatn(I+1)=(concentration(I+1)-Csat(Temp(T)))/
& Csat(Temp(T))
time(I+1)=T
K = 20
DO 363 M = 1 , NU
DO 364 N = 1 , NTHETA
derivtheta(I+1,M,N)=Y(K)
F((NSTEP-1)*(Iteration-1)+Nu*(I-1)+M,N)=Y(k)
K = K + 1
364 CONTINUE
363 CONTINUE
conc_variance=conc_variance+
& (ConcData(I)-concentration(I+1))**2.0D0
IF ((int(I*Cinter/Minter/L).EQ.1) .AND. (L .LE. MSTEP)) THEN
mu00_variance=mu00_variance+
& (Mu00Data(L)-moment00(I+1))**2.0D0
mu10_variance=mu10_variance+
& (Mu10Data(L)-moment10(I+1))**2.0D0
mu01_variance=mu01_variance+
& (Mu01Data(L)-moment01(I+1))**2.0D0
mu11_variance=mu11_variance+
& (Mu11Data(L)-moment11(I+1))**2.0D0
mu20_variance=mu20_variance+
& (Mu20Data(L)-moment20(I+1))**2.0D0
mu02_variance=mu02_variance+
& (Mu02Data(L)-moment02(I+1))**2.0D0
mu21_variance=mu21_variance+
& (Mu21Data(L)-moment21(I+1))**2.0D0
mu12_variance=mu12_variance+
& (Mu12Data(L)-moment12(I+1))**2.0D0
mu22_variance=mu22_variance+
& (Mu22Data(L)-moment22(I+1))**2.0D0
mu30_variance=mu30_variance+
& (Mu30Data(L)-moment30(I+1))**2.0D0
mu31_variance=mu31_variance+
& (Mu31Data(L)-moment31(I+1))**2.0D0
L=L+1
END IF
1200 CONTINUE
787 CONTINUE
1 FORMAT(2(F13.6,3x))
2 FORMAT(I1)
3 FORMAT(2(E16.6,3x),2(F16.6,3x))
mu00_variance=mu00_variance/DFLOAT((NSTEP-1)*Nsets)
mu10_variance=mu10_variance/DFLOAT((NSTEP-1)*Nsets)
mu01_variance=mu01_variance/DFLOAT((NSTEP-1)*Nsets)
mu11_variance=mu11_variance/DFLOAT((NSTEP-1)*Nsets)
mu20_variance=mu20_variance/DFLOAT((NSTEP-1)*Nsets)
mu02_variance=mu02_variance/DFLOAT((NSTEP-1)*Nsets)
mu21_variance=mu21_variance/DFLOAT((NSTEP-1)*Nsets)
mu12_variance=mu12_variance/DFLOAT((NSTEP-1)*Nsets)
mu22_variance=mu22_variance/DFLOAT((NSTEP-1)*Nsets)
mu30_variance=mu30_variance/DFLOAT((NSTEP-1)*Nsets)
mu31_variance=mu31_variance/DFLOAT((NSTEP-1)*Nsets)
conc_variance=conc_variance/DFLOAT((NSTEP-1)*Nsets)
*
*
DO 1370 I = 1, NTHETA
DO 1375 J = 1, Nu*(NSTEP-1)*Nsets, Nu
FTV(I,J)=F(J,I)/mu00_variance
FTV(I,J+1)=F(J+1,I)/mu10_variance
FTV(I,J+2)=F(J+2,I)/mu01_variance
FTV(I,J+3)=F(J+3,I)/mu11_variance
FTV(I,J+4)=F(J+4,I)/mu20_variance
FTV(I,J+5)=F(J+5,I)/mu02_variance
FTV(I,J+6)=F(J+6,I)/mu21_variance
FTV(I,J+7)=F(J+7,I)/mu12_variance
FTV(I,J+8)=F(J+8,I)/mu22_variance
FTV(I,J+9)=F(J+9,I)/mu30_variance
FTV(I,J+10)=F(J+10,I)/mu31_variance
FTV(I,J+11)=F(J+11,I)/conc_variance
1375 CONTINUE
1370 CONTINUE
DO 1400 I = 1, NTHETA
DO 1500 J = 1, NTHETA
FTVF(I,J) = 0.0D0
DO 1600 K = 1, Nu*(NSTEP-1)*Nsets
FTVF(I,J)=FTVF(I,J)+FTV(I,K)*F(K,J)
1600 CONTINUE
1500 CONTINUE
1400 CONTINUE
*Calculate determinant of FTVF,
*which is symmetric with dimensions NTHETA by NTHETA
*
DO 200 I = 1, NTHETA
DO 300 J = 1, NTHETA
FAD(I,J) = 0.0D0
300 CONTINUE
200 CONTINUE
CALL DLFTRG (P, FTVF, NTHETA, FAD, LDFAC, IPVT)
CALL DLFDRG (P,FAD,LDFAC,IPVT,DET1,DET2)
detFTVF = DET1*10**(DET2)
PRINT*, "Determinant of FTVF = ", detFTVF
PRINT*, '-natural log (determinant of FTVF)= ',
& -DLOG(detFTVF)
PRINT*,'FTVF = '
WRITE(*,11)FTVF(1,1), FTVF(1,2),FTVF(1,3),FTVF(1,4),
&FTVF(1,5),FTVF(1,6)
WRITE(*,11)FTVF(2,1), FTVF(2,2),FTVF(2,3),FTVF(2,4),
&FTVF(2,5),FTVF(2,6)
WRITE(*,11)FTVF(3,1), FTVF(3,2),FTVF(3,3),FTVF(3,4),
&FTVF(3,5),FTVF(3,6)
WRITE(*,11)FTVF(4,1), FTVF(4,2),FTVF(4,3),FTVF(4,4),
&FTVF(4,5),FTVF(4,6)
WRITE(*,11)FTVF(5,1), FTVF(5,2),FTVF(5,3),FTVF(5,4),
&FTVF(5,5),FTVF(5,6)
WRITE(*,11)FTVF(6,1), FTVF(6,2),FTVF(6,3),FTVF(6,4),
&FTVF(6,5),FTVF(6,6)
11 FORMAT(4(E24.16,2x))
CALL DLINRG(NTHETA,FTVF,NTHETA, inv_FTVF,NTHETA)
* From chi squared distribution
chi_squared=9.49D0
g1_interval=DSQRT(chi_squared*inv_FTVF(1,1))
lnkg1_interval=DSQRT(chi_squared*inv_FTVF(2,2))
g2_interval=DSQRT(chi_squared*inv_FTVF(1,1))
lnkg2_interval=DSQRT(chi_squared*inv_FTVF(2,2))
b_interval=DSQRT(chi_squared*inv_FTVF(3,3))
lnkb_interval=DSQRT(chi_squared*inv_FTVF(4,4))
PRINT*
PRINT*,'Using chi-squared = ',chi_squared
PRINT*
PRINT*,'g1 = ', g1, ' +/- ', g1_interval
PRINT*,'ln kg1 = ', DLOG(kg1), ' +/- ', lnkg1_interval
PRINT*,'g2 = ', g2, ' +/- ', g2_interval
PRINT*,'ln kg2 = ', DLOG(kg2), ' +/- ', lnkg2_interval
PRINT*,' b = ', b, ' +/- ', b_interval
PRINT*,'ln kb = ', DLOG(kb), ' +/- ', lnkb_interval
pause
PRINT*
PRINT*, 'mu00_variance = ', mu00_variance
PRINT*, 'mu10_variance = ', mu10_variance
PRINT*, 'mu01_variance = ', mu01_variance
PRINT*, 'mu11_variance = ', mu11_variance
PRINT*, 'mu20_variance = ', mu20_variance
PRINT*, 'mu02_variance = ', mu02_variance
PRINT*, 'mu21_variance = ', mu21_variance
PRINT*, 'mu12_variance = ', mu12_variance
PRINT*, 'mu22_variance = ', mu22_variance
PRINT*, 'mu30_variance = ', mu30_variance
PRINT*, 'mu31_variance = ', mu31_variance
PRINT*, 'conc_variance = ', conc_variance
pause
PRINT*, 'mu00_sigma = ', DSQRT(mu00_variance)
PRINT*, 'mu10_sigma = ', DSQRT(mu10_variance)
PRINT*, 'mu01_sigma = ', DSQRT(mu01_variance)
PRINT*, 'mu11_sigma = ', DSQRT(mu11_variance)
PRINT*, 'mu20_sigma = ', DSQRT(mu20_variance)
PRINT*, 'mu02_sigma = ', DSQRT(mu02_variance)
PRINT*, 'mu21_sigma = ', DSQRT(mu21_variance)
PRINT*, 'mu12_sigma = ', DSQRT(mu12_variance)
PRINT*, 'mu22_sigma = ', DSQRT(mu22_variance)
PRINT*, 'mu30_sigma = ', DSQRT(mu30_variance)
PRINT*, 'mu31_sigma = ', DSQRT(mu31_variance)
PRINT*, 'conc_sigma = ', DSQRT(conc_variance)
OPEN(UNIT=50, FILE='FTVF.dat', FORM='FORMATTED',
& ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
OPEN(UNIT=60, FILE='variance.dat', FORM='FORMATTED',
& ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
OPEN(UNIT=70, FILE='interval.dat', FORM='FORMATTED',
& ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
WRITE(50,11)FTVF(1,1), FTVF(1,2),FTVF(1,3),FTVF(1,4),
& FTVF(1,5), FTVF(1,6)
WRITE(50,11)FTVF(2,1), FTVF(2,2),FTVF(2,3),FTVF(2,4),
& FTVF(2,5), FTVF(2,6)
WRITE(50,11)FTVF(3,1), FTVF(3,2),FTVF(3,3),FTVF(3,4),
& FTVF(3,5), FTVF(3,6)
WRITE(50,11)FTVF(4,1), FTVF(4,2),FTVF(4,3),FTVF(4,4),
& FTVF(4,5), FTVF(4,6)
WRITE(50,11)FTVF(5,1), FTVF(5,2),FTVF(5,3),FTVF(5,4),
& FTVF(5,5), FTVF(5,6)
WRITE(50,11)FTVF(6,1), FTVF(6,2),FTVF(6,3),FTVF(6,4),
& FTVF(6,5), FTVF(6,6)
WRITE(60,12) mu00_variance
WRITE(60,12) mu10_variance
WRITE(60,12) mu01_variance
WRITE(60,12) mu11_variance
WRITE(60,12) mu20_variance
WRITE(60,12) mu02_variance
WRITE(60,12) mu21_variance
WRITE(60,12) mu12_variance
WRITE(60,12) mu22_variance
WRITE(60,12) mu30_variance
WRITE(60,12) mu31_variance
WRITE(60,12) conc_variance
WRITE(70,12) g1_interval
WRITE(70,12) lnkg1_interval
WRITE(70,12) g2_interval
WRITE(70,12) lnkg2_interval
WRITE(70,12) b_interval
WRITE(70,12) lnkb_interval
CLOSE(UNIT=50)
CLOSE(UNIT=60)
CLOSE(UNIT=70)
CLOSE(UNIT=90)
12 FORMAT(E24.16)
print*, "done"
STOP
END
****************************************************************
SUBROUTINE FCN(NEQ,T,Y,Yprime)
INTEGER NEQ, I, J, K
INTEGER NTHETA, NU
REAL*8 T, Y(NEQ), Yprime(NEQ)
REAL*8 r0, alpha, mu00, kg1, g1, kg2, g2, kb, b
REAL*8 Csat, birth, birth2, growth1, growth2, growth3, growth4
REAL*8 Temp
REAL*8 JACOBIAN(12,12), W(12,6)
REAL*8 WPRIME(12,6), F_THETA(12,6)
COMMON /EXP_DATA/r0, alpha, mu00
COMMON /GROWTH_DATA1/kg1, g1
COMMON /GROWTH_DATA2/kg2, g2
COMMON /BIRTH_DATA/kb, b
NTHETA=6
NU=12
* print*, 'in moments'
*Moments
* 00 moment, number of particles/g solvent/minute
Yprime(1)=birth(Y(12),Csat(Temp(T)),Y(7))
* 10 moment, microns/g solvent/minute
Yprime(2)=growth1(Y(12),Csat(Temp(T)))*Y(1)
* 01 moment, microns/g solvent/minute
Yprime(3)=growth2(Y(12),Csat(Temp(T)))*Y(1)
* 11 moment, microns^2/g solvent/minute
Yprime(4)=growth1(Y(12),Csat(Temp(T)))*Y(3)
& + growth2(Y(12),Csat(Temp(T)))*Y(2)
* 20 moment, microns^2/g solvent/minute
Yprime(5)=2*growth1(Y(12),Csat(Temp(T)))*Y(2)
* 02 moment, microns^2/g solvent/minute
Yprime(6)=2*growth2(Y(12),Csat(Temp(T)))*Y(3)
* 21 moment, microns^3/g solvent/minute
Yprime(7)=2*growth1(Y(12),Csat(Temp(T)))*Y(4)+
& growth2(Y(12),Csat(Temp(T)))*Y(5)
* 12 moment, microns^3/g solvent/minute
Yprime(8)=growth1(Y(12),Csat(Temp(T)))*Y(6)+
& 2*growth2(Y(12),Csat(Temp(T)))*Y(4)
* 22 moment, microns^4/g solvent/minute
Yprime(9)=2*growth1(Y(12),Csat(Temp(T)))*Y(8)+
& 2*growth2(Y(12),Csat(Temp(T)))*Y(7)
* 30 moment, microns^4/g solvent/minute
Yprime(10)=3*growth1(Y(12),Csat(Temp(T)))*Y(5)
* 31 moment, microns^4/g solvent/minute
Yprime(11)=3*growth1(Y(12),Csat(Temp(T)))*Y(7)+
& growth2(Y(12),Csat(Temp(T)))*Y(10)
*Mass balance
* d(conc(t))/dt, g solute/g solvent/minute
Yprime(12)=-alpha*2.0D0*growth1( Y(12), Csat(Temp(T)) )
& *Y(4) - alpha*growth2( Y(12),Csat(Temp(T) ) )*Y(5)
*Moments for seed crystals only
* 10 moment, microns/g solvent/minute
Yprime(13)=growth1(Y(12),Csat(Temp(T)))*mu00
* 01 moment, microns/g solvent/minute
Yprime(14)=growth2(Y(12),Csat(Temp(T)))*mu00
* 11 moment, microns^2/g solvent/minute
Yprime(15)=growth1(Y(12),Csat(Temp(T)))*Y(14)
& + growth2(Y(12),Csat(Temp(T)))*Y(13)
* 20 moment, microns^2/g solvent/minute
Yprime(16)=2*growth1(Y(12),Csat(Temp(T)))*Y(13)
* 02 moment, microns^2/g solvent/minute
Yprime(17)=2*growth2(Y(12),Csat(Temp(T)))*Y(14)
* 21 moment, microns^3/g solvent/minute
Yprime(18)=2*growth1(Y(12),Csat(Temp(T)))*Y(15)+
& growth2(Y(12),Csat(Temp(T)))*Y(16)
* 12 moment, microns^3/g solvent/minute
Yprime(19)=growth1(Y(12),Csat(Temp(T)))*Y(17)+
& 2*growth2(Y(12),Csat(Temp(T)))*Y(15)
DO 203 I = 1 , NU
DO 204 J = 1 , NU
JACOBIAN(I,J)=0.0D0
204 CONTINUE
203 CONTINUE
* print*, 'out moments'
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -