Dear experts, I have a FORTRAN code, this code is evaluating a variable in between but i am not able to write that result, Please help. following is code, i want to fetch data of “TE” for each iteration of “TC”.
MODULE GLOBAL
IMPLICIT NONE
!INTEGER,PARAMETER :: NI=202
!INTEGER,PARAMETER :: NJ=202
!INTEGER,PARAMETER :: NIJ=202
INTEGER,PARAMETER:: J=351
DOUBLE PRECISION,SAVE :: TC,TA,XC,XA,TB,TD,XB,XD,TE,XE,XF,Q_SH,Q_SC,Q_DE,Q_AD
DOUBLE PRECISION,SAVE :: GAMATB,GAMATD
DOUBLE PRECISION,SAVE :: TAVG,PCSAT,PESAT,PMEAN,TCSAT,TESAT,TMEAN_SAT
DOUBLE PRECISION, PARAMETER :: K=10.226
DOUBLE PRECISION, PARAMETER :: N=1.99
DOUBLE PRECISION, PARAMETER :: R=0.488
DOUBLE PRECISION, PARAMETER :: C=2823
DOUBLE PRECISION, PARAMETER :: X0=0.318
DOUBLE PRECISION, PARAMETER :: CAC_REF=0.8485
DOUBLE PRECISION, PARAMETER :: CVR=2.344
DOUBLE PRECISION, PARAMETER :: CPR=4.9
DOUBLE PRECISION, PARAMETER :: CHTF=1.8024
DOUBLE PRECISION, PARAMETER :: HEVAP=1368.0
DOUBLE PRECISION,SAVE :: RM,AGAMA,XGAMA,IGAMA
DOUBLE PRECISION,SAVE,DIMENSION(J) :: TCVAR,COP,SCE
END MODULE GLOBAL
!======================================================================
!===> PROGRAM MAIN ***********************************************
!======================================================================
PROGRAM MAIN
USE GLOBAL
IMPLICIT NONE
WRITE(*,*)"PROGRAM RUNNING MASS RECOVERY"
CALL COPCALC
! CALL GEOMET
! CALL START
!MAINLOOP: DO WHILE(.NOT.LSTOP)
! !CALL DENSE
! CALL BOUND
! CALL OLDVAL
! CALL COEFF
! END DO MAINLOOP
WRITE(*,*)"PROGRAM FINISHED SOLIDIFICATION"
END PROGRAM MAIN
SUBROUTINE COPCALC
USE GLOBAL
IMPLICIT NONE
INTEGER :: I
DOUBLE PRECISION :: TCSTART
DOUBLE PRECISION :: GAMATC,GAMATE,Q_DE_E,Q_REG
DOUBLE PRECISION :: GAMATA,GAMATCSAT,XTCSAT
DOUBLE PRECISION :: Q_TOTIN,Q_EVAP,Q_CO,Q_LTTCSAT,COPDIF,COPMAX,TCOPT
LOGICAL :: FINDMAX
CHARACTER (LEN=80) :: FILENAME
!======================================================================
!********************** INPUT DATA HERE *******************************
TA=300.0 ; TCSAT=300.0 ; TESAT=273.0;
RM=0.0
!======================================================================
TCSTART=350
DO I=1,J
TCVAR(I)=TCSTART+(I-1)*1.0
TC=TCVAR(I)
XC=X0*EXP(-K*(TC/TCSAT-1)**N)
XA=X0*EXP(-K*(TA/TESAT-1)**N)
TB=((LOG(X0/XA)/K)**(1.0/N)+1.0)*TCSAT
TD=((LOG(X0/XC)/K)**(1.0/N)+1.0)*TESAT
CALL GRID
Q_SH=0.175*(TB-TA)+2.245E-3*0.5*(TB**2-TA**2)
Q_SH=Q_SH+RM*CAC_REF*(TB-TA)+CVR*XA*(TB-TA)
Q_SC=0.175*(TC-TD)+2.245E-3*0.5*(TC**2-TD**2)
Q_SC=Q_SC+RM*CAC_REF*(TC-TD)+CVR*XC*(TC-TD)
Q_DE=0.175*(TC-TB)+2.245E-3*0.5*(TC**2-TB**2)
Q_DE=Q_DE+RM*CAC_REF*(TC-TB)
AGAMA=1/N
XGAMA=K*((TB/TCSAT-1.0)**N)
CALL INGAMA2ND
GAMATB=IGAMA
AGAMA=1/N
XGAMA=K*((TC/TCSAT-1.0)**N)
CALL INGAMA2ND
GAMATC=IGAMA
Q_DE=Q_DE+CPR*X0*TCSAT*(GAMATB-GAMATC)/(N*(K**(1/N)))
XB=XA
AGAMA=(1.0+N)/N
XGAMA=LOG(X0/XC)
CALL INGAMA2ND
GAMATC=IGAMA
XGAMA=LOG(X0/XB)
CALL INGAMA2ND
GAMATB=IGAMA
Q_DE=Q_DE-C*R*(XC-XB)
Q_DE=Q_DE-C*R*X0*(GAMATC-GAMATB)/(K**(1/N))
Q_AD=0.175*(TD-TA)+2.245E-3*0.5*(TD**2-TA**2)
Q_AD=Q_AD+RM*CAC_REF*(TD-TA)
AGAMA=1/N
XGAMA=K*((TD/TESAT-1.0)**N)
CALL INGAMA2ND
GAMATD=IGAMA
AGAMA=1/N
XGAMA=K*((TA/TESAT-1.0)**N)
CALL INGAMA2ND
GAMATA=IGAMA
Q_AD=Q_AD+CPR*X0*TESAT*(GAMATA-GAMATD)/(N*(K**(1/N)))
XD=XC
AGAMA=(1.0+N)/N
XGAMA=LOG(X0/XA)
CALL INGAMA2ND
GAMATA=IGAMA
XGAMA=LOG(X0/XD)
CALL INGAMA2ND
GAMATD=IGAMA
Q_AD=Q_AD-C*R*(XD-XA)
Q_AD=Q_AD-C*R*X0*(GAMATD-GAMATA)/(K**(1/N))
Q_AD=Q_AD+CPR*TESAT*X0*(GAMATD-GAMATA)/(K**(1/N))
!CALCULATION OF Q_REG
IF(TB.GE.TD) THEN
Q_DE_E=0.0
Q_REG=(0.175+RM*CAC_REF+CVR*XC)*(TC-TE)+(2.245E-3)*0.5*(TC**2-TE**2)
ELSE
Q_DE_E=0.175*(TE-TB)+2.245E-3*0.5*(TE**2-TB**2)
Q_DE_E=Q_DE_E+RM*CAC_REF*(TE-TB)
AGAMA=1/N
XGAMA=K*((TB/TCSAT-1.0)**N)
CALL INGAMA2ND
GAMATB=IGAMA
AGAMA=1/N
XGAMA=K*((TE/TCSAT-1.0)**N)
CALL INGAMA2ND
GAMATE=IGAMA
Q_DE_E=Q_DE_E+CPR*X0*TCSAT*(GAMATB-GAMATE)/(N*(K**(1/N)))
XB=XA
XF=X0*EXP(-K*(TE/TCSAT-1)**N)
AGAMA=(1.0+N)/N
XGAMA=LOG(X0/XF)
CALL INGAMA2ND
GAMATE=IGAMA
XGAMA=LOG(X0/XB)
CALL INGAMA2ND
GAMATB=IGAMA
Q_DE_E=Q_DE_E-C*R*(XF-XB)
Q_DE_E=Q_DE_E-C*R*X0*(GAMATE-GAMATB)/(K**(1/N))
Q_REG=Q_SH+Q_DE_E
END IF
Q_TOTIN=Q_SH+Q_DE-Q_REG
Q_EVAP=(XB-XC)*HEVAP
Q_CO=4.043*(XB-XC)*(TCSAT-TESAT)-4.1E-3*0.5*(XB-XC)*(TCSAT**2-TESAT**2)
IF(TA.GT.TCSAT) THEN
Q_LTTCSAT=0.0
ELSE
Q_LTTCSAT=0.175*(TCSAT-TA)+2.245E-3*0.5*(TCSAT**2-TA**2)
Q_LTTCSAT=Q_LTTCSAT+RM*CAC_REF*(TCSAT-TA)
AGAMA=1/N
XGAMA=K*((TCSAT/TESAT-1.0)**N)
CALL INGAMA2ND
GAMATCSAT=IGAMA
AGAMA=1/N
XGAMA=K*((TA/TESAT-1.0)**N)
CALL INGAMA2ND
GAMATA=IGAMA
Q_LTTCSAT=Q_LTTCSAT+CPR*X0*TESAT*(GAMATA-GAMATCSAT)/(N*(K**(1/N)))
XTCSAT=X0*EXP(-K*(TCSAT/TESAT-1)**N)
AGAMA=(1.0+N)/N
XGAMA=LOG(X0/XA)
CALL INGAMA2ND
GAMATA=IGAMA
XGAMA=LOG(X0/XTCSAT)
CALL INGAMA2ND
GAMATCSAT=IGAMA
Q_LTTCSAT=Q_LTTCSAT-C*R*(XTCSAT-XA)
Q_LTTCSAT=Q_LTTCSAT-C*R*X0*(GAMATCSAT-GAMATA)/(K**(1/N))
Q_LTTCSAT=Q_LTTCSAT+CPR*TESAT*X0*(GAMATCSAT-GAMATA)/(K**(1/N))
END IF
COP(I)=(Q_EVAP-Q_CO-Q_LTTCSAT)/Q_TOTIN
SCE(I)=Q_EVAP-Q_CO-Q_LTTCSAT
WRITE(*,*)'TB=',TB,'I=',I,'TCVAR=',TCVAR(I),'XA=',XA
WRITE(*,*)'TA=',TA,'Q_SH=',Q_SH,'Q_DE_E=',Q_DE_E
WRITE(*,*)'TD=',TD,'XC=',XC,'TE=',TE
WRITE(*,*)'TC=',TC,'Q_SC=',Q_SC
WRITE(*,*)'Q_DE=',Q_DE,'Q_AD=',Q_AD
WRITE(*,*)'Q_REG=',Q_REG
WRITE(*,*)'Q_TOTIN=',Q_TOTIN ,'Q_EVAP=',Q_EVAP,'Q_CO=',Q_CO
WRITE(*,*)'Q_LTTCSAT=',Q_LTTCSAT,'COP',COP(I),'SCE',SCE(I)
WRITE(*,*)'=========================================================='
END DO
FINDMAX=.FALSE.
I=1
DO WHILE(.NOT.FINDMAX)
I=I+1
COPDIF=COP(I+1)-COP(I)
IF(COPDIF.LT.0.0) THEN
FINDMAX=.TRUE.
COPMAX=COP(I)
TCOPT=TCVAR(I)
END IF
END DO
WRITE(FILENAME,FMT='(A6,I2.2,A2)') "output", INT(RM*10), ".m"
OPEN(UNIT=109,FILE=FILENAME)
WRITE(109,*)"close all"
WRITE(109,*)"clear all"
WRITE(109,*)"TCOPT=",TCOPT
WRITE(109,*)"COPMAX=",COPMAX
WRITE(109,*) "TC=["
DO I=1,J
WRITE(109,*) TCVAR(I)
END DO
WRITE(109,*)"];"
WRITE(109,*) "COP=["
DO I=1,J
WRITE(109,*) COP(I)
END DO
WRITE(109,*)"];"
WRITE(109,*)"figure(1);"
WRITE(109,*)"plot(TC,COP)"
WRITE(109,*)"xlabel('TC (K)')"
WRITE(109,*)"ylabel('COP')"
WRITE(109,*) "SCE=["
DO I=1,J
WRITE(109,*) SCE(I)
END DO
WRITE(109,*)"];"
WRITE(109,*)"figure(2);"
WRITE(109,*)"plot(TC,SCE)"
WRITE(109,*)"xlabel('TC (K)')"
WRITE(109,*)"ylabel('SCE (kJ/kg)')"
WRITE(*,*)'I AM HERE'
CLOSE(109)
END SUBROUTINE COPCALC
SUBROUTINE GRID
USE GLOBAL
IMPLICIT NONE
INTEGER :: I
DOUBLE PRECISION :: TE1,TE2,TEMP1,TEMP2,DIF_TE,GAMATE1,GAMATE2
!DOUBLE PRECISION :: TA_PRIME1,TA_PRIME2,DIF_TA_PRIME,DIF_TC_PRIME
!DOUBLE PRECISION :: XC_PRIME1,XC_PRIME2,IGAMAC,IGAMAA
!DOUBLE PRECISION :: TC_PRIMEO,TA_PRIMEO,XC_PRIMEO,XA_PRIMEO,PMEANO,TMEAN_SATO
!DOUBLE PRECISION :: DIF_TC_P,DIF_TA_P,DIF_XC_P,DIF_XA_P,DIF_PMEAN
DOUBLE PRECISION,PARAMETER :: ERT=1.0E-8
!DOUBLE PRECISION :: ERRT,ERRX,ERRP,REL
LOGICAL :: LCONV1
!======================================================================
!********************** INPUT DATA HERE *******************************
LCONV1=.FALSE.
I=1
DO WHILE(.NOT.LCONV1)
IF(TB.GE.TD) THEN
IF(I.EQ.1) THEN
TE1=TB
TE2=TD
ELSE
TE1=TE2
TE2=TE
END IF
TEMP1=(2.245E-3)*(TE1)**2+2*(0.175+RM*CAC_REF)*TE1+CVR*(XA+XC)*TE1
TEMP1=TEMP1-((2.245E-3)*0.5*(TA)**2+(0.175+RM*CAC_REF)*TA+CVR*XA*TA)
TEMP1=TEMP1-((2.245E-3)*0.5*(TC)**2+(0.175+RM*CAC_REF)*TC+CVR*XC*TC)
TEMP2=(2.245E-3)*(TE2)**2+2*(0.175+RM*CAC_REF)*TE2+CVR*(XA+XC)*TE2
TEMP2=TEMP2-((2.245E-3)*0.5*(TA)**2+(0.175+RM*CAC_REF)*TA+CVR*XA*TA)
TEMP2=TEMP2-((2.245E-3)*0.5*(TC)**2+(0.175+RM*CAC_REF)*TC+CVR*XC*TC)
!CALCULATION OF TE USING SECANT METHOD
TE=TE2-TEMP2*(TE2-TE1)/(TEMP2-TEMP1)
ELSE
IF(I.EQ.1) THEN
TE1=TB
TE2=TD
ELSE
TE1=TE2
TE2=TE
END IF
!CALCULATION OF F1 STARTS HERE
Q_SH=0.175*(TB-TA)+2.245E-3*0.5*(TB**2-TA**2)
Q_SH=Q_SH+RM*CAC_REF*(TB-TA)+CVR*XA*(TB-TA)
Q_SC=0.175*(TC-TD)+2.245E-3*0.5*(TC**2-TD**2)
Q_SC=Q_SC+RM*CAC_REF*(TC-TD)+CVR*XC*(TC-TD)
XE=X0*EXP(-K*(TE1/TESAT-1)**N)
XF=X0*EXP(-K*(TE1/TCSAT-1)**N)
Q_DE=0.175*(TE1-TB)+2.245E-3*0.5*(TE1**2-TB**2)
Q_DE=Q_DE+RM*CAC_REF*(TE1-TB)
AGAMA=1/N
XGAMA=K*((TB/TCSAT-1.0)**N)
CALL INGAMA2ND
GAMATB=IGAMA
AGAMA=1/N
XGAMA=K*((TE1/TCSAT-1.0)**N)
CALL INGAMA2ND
GAMATE1=IGAMA
Q_DE=Q_DE+CPR*X0*TCSAT*(GAMATB-GAMATE1)/(N*(K**(1/N)))
XB=XA
AGAMA=(1.0+N)/N
XGAMA=LOG(X0/XF)
CALL INGAMA2ND
GAMATE1=IGAMA
XGAMA=LOG(X0/XB)
CALL INGAMA2ND
GAMATB=IGAMA
Q_DE=Q_DE-C*R*(XF-XB)
Q_DE=Q_DE-C*R*X0*(GAMATE1-GAMATB)/(K**(1/N))
Q_AD=0.175*(TD-TE1)+2.245E-3*0.5*(TD**2-TE1**2)
Q_AD=Q_AD+RM*CAC_REF*(TD-TE1)
AGAMA=1/N
XGAMA=K*((TD/TESAT-1.0)**N)
CALL INGAMA2ND
GAMATD=IGAMA
AGAMA=1/N
XGAMA=K*((TE1/TESAT-1.0)**N)
CALL INGAMA2ND
GAMATE1=IGAMA
Q_AD=Q_AD+CPR*X0*TESAT*(GAMATE1-GAMATD)/(N*(K**(1/N)))
XD=XC
AGAMA=(1.0+N)/N
XGAMA=LOG(X0/XE)
CALL INGAMA2ND
GAMATE1=IGAMA
XGAMA=LOG(X0/XD)
CALL INGAMA2ND
GAMATD=IGAMA
Q_AD=Q_AD-C*R*(XD-XE)
Q_AD=Q_AD-C*R*X0*(GAMATD-GAMATE1)/(K**(1/N))
Q_AD=Q_AD+CPR*TESAT*X0*(GAMATD-GAMATE1)/(K**(1/N))
TEMP1=Q_SH+Q_DE-Q_SC-Q_AD
!CALCULATION OF F2 STARTS HERE
Q_SH=0.175*(TB-TA)+2.245E-3*0.5*(TB**2-TA**2)
Q_SH=Q_SH+RM*CAC_REF*(TB-TA)+CVR*XA*(TB-TA)
Q_SC=0.175*(TC-TD)+2.245E-3*0.5*(TC**2-TD**2)
Q_SC=Q_SC+RM*CAC_REF*(TC-TD)+CVR*XC*(TC-TD)
XE=X0*EXP(-K*(TE2/TESAT-1)**N)
XF=X0*EXP(-K*(TE2/TCSAT-1)**N)
Q_DE=0.175*(TE2-TB)+2.245E-3*0.5*(TE2**2-TB**2)
Q_DE=Q_DE+RM*CAC_REF*(TE2-TB)
AGAMA=1/N
XGAMA=K*((TB/TCSAT-1.0)**N)
CALL INGAMA2ND
GAMATB=IGAMA
AGAMA=1/N
XGAMA=K*((TE2/TCSAT-1.0)**N)
CALL INGAMA2ND
GAMATE2=IGAMA
Q_DE=Q_DE+CPR*X0*TCSAT*(GAMATB-GAMATE2)/(N*(K**(1/N)))
XB=XA
AGAMA=(1.0+N)/N
XGAMA=LOG(X0/XF)
CALL INGAMA2ND
GAMATE2=IGAMA
XGAMA=LOG(X0/XB)
CALL INGAMA2ND
GAMATB=IGAMA
Q_DE=Q_DE-C*R*(XF-XB)
Q_DE=Q_DE-C*R*X0*(GAMATE2-GAMATB)/(K**(1/N))
Q_AD=0.175*(TD-TE2)+2.245E-3*0.5*(TD**2-TE2**2)
Q_AD=Q_AD+RM*CAC_REF*(TD-TE2)
AGAMA=1/N
XGAMA=K*((TD/TESAT-1.0)**N)
CALL INGAMA2ND
GAMATD=IGAMA
AGAMA=1/N
XGAMA=K*((TE2/TESAT-1.0)**N)
CALL INGAMA2ND
GAMATE2=IGAMA
Q_AD=Q_AD+CPR*X0*TESAT*(GAMATE2-GAMATD)/(N*(K**(1/N)))
XD=XC
AGAMA=(1.0+N)/N
XGAMA=LOG(X0/XE)
CALL INGAMA2ND
GAMATE2=IGAMA
XGAMA=LOG(X0/XD)
CALL INGAMA2ND
GAMATD=IGAMA
Q_AD=Q_AD-C*R*(XD-XE)
Q_AD=Q_AD-C*R*X0*(GAMATD-GAMATE2)/(K**(1/N))
Q_AD=Q_AD+CPR*TESAT*X0*(GAMATD-GAMATE2)/(K**(1/N))
TEMP2=Q_SH+Q_DE-Q_SC-Q_AD
!CALCULATION OF TE USING SECANT METHOD
TE=TE2-TEMP2*(TE2-TE1)/(TEMP2-TEMP1)
END IF
DIF_TE=ABS(TE-TE2)
! WRITE(*,*)I,DIF_TE
IF(DIF_TE.LT.ERT) LCONV1=.TRUE.
I=I+1
END DO
END SUBROUTINE GRID
SUBROUTINE INGAMA2ND
USE GLOBAL
IMPLICIT NONE
INTEGER :: ICOUNT
DOUBLE PRECISION, PARAMETER :: ACU=1.0E-15
DOUBLE PRECISION :: DELXGAMA,IGAMA1,IGAMA2,IGAMAOLD,IGAMADIF
LOGICAL :: LCONV
DELXGAMA=XGAMA/10.0
LCONV=.FALSE.
IGAMA=0.0
IGAMA1=(XGAMA**(AGAMA-1.0))*EXP(-XGAMA)
ICOUNT=1
IGAMAOLD=0.0
DO WHILE(.NOT.LCONV)
IGAMA2=((XGAMA+ICOUNT*DELXGAMA)**(AGAMA-1.0))*EXP(-(XGAMA+ICOUNT*DELXGAMA))
IGAMA=IGAMA+0.5*(IGAMA1+IGAMA2)*DELXGAMA
ICOUNT=ICOUNT+1
IGAMA1=IGAMA2
IGAMADIF=ABS(IGAMA-IGAMAOLD)
IF(IGAMADIF.LT.ACU) LCONV=.TRUE.
IGAMAOLD=IGAMA
END DO
!WRITE(*,*)'IGAMA=',IGAMA,'ICOUNT=',ICOUNT
!STOP
END SUBROUTINE INGAMA2ND