Help regarding Fortran code

Help
1

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
2

Welcome to this forum. You posted the code as ordinary text, which means it suffers from all manner of markup. Could you repost it as code - use the </> symbol.

Also, indicate where in the code you would like to see a write statement. It is unclear to me where this should go.

And finally, can you indicate what you tried to write out the variable of interest, as it is rather basic and you have already some write statements as examples.

4 



MODULE GLOBAL  
IMPLICIT NONE  
  
!INTEGER,PARAMETER :: NI=202 
!INTEGER,PARAMETER :: NJ=202   
!INTEGER,PARAMETER :: NIJ=202  

INTEGER,PARAMETER:: J=152

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=2.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
`
5

Thanks for your reply sir,
i want to write “TE” for every “TC”

6

If you are only interested into the values of TE (because TC seems to regularly increment by 1.0), redirect the output of the terminal into a permanent record and extract them from there it with a regular expression. As an example (instance of Linux Debian 13)

gfortran test.f90 -o executable
./executable > report.txt
awk '{if (/TE/) print $6}' report.txt > my_TE_data.txt

In case you work with Windows which isn’t shipped with AWK, first install git followed by tortoise git as a graphical interface to git. Among the icons the later adds to the pull-down menu of the file manager is one reading Git Bash. See here:

cropped

A click on this opens a new shell which provides a couple of elementary Linux tools* including Gnu AWK 5.0.0 (at time of writing). If AWK is new to you, Benjamin Porter compiled Awk: Hack the planet['s text]! including test data. Else for instance the gray first, and the 2023 second edition of the AWK book by the language’s authors.

* It is constrained compared to a (parallel) installation of Linux, but lighter in footprint than e.g. a Linux WSL on Windows.

7

I’d amend the code so that the output is handled in a separate procedure. For example:

SUBROUTINE COPCALC
    USE GLOBAL  
    use iso_fortran_env, only: output_unit, error_unit
 ...
        COP(I)=(Q_EVAP-Q_CO-Q_LTTCSAT)/Q_TOTIN
        SCE(I)=Q_EVAP-Q_CO-Q_LTTCSAT
        call output_iteration(.true.)
    END DO
...
    contains
    
    subroutine output_iteration(debugging)
        logical, intent(in) :: debugging
        if (debugging) then
            write(error_unit,'(f0.1,a,f0.6)') TC,' ',TE
        else
            write(output_unit,*)'TB=',TB,'I=',I,'TCVAR=',TCVAR(I),'XA=',XA
            write(output_unit,*)'TA=',TA,'Q_SH=',Q_SH,'Q_DE_E=',Q_DE_E
            write(output_unit,*)'TD=',TD,'XC=',XC,'TE=',TE
            write(output_unit,*)'TC=',TC,'Q_SC=',Q_SC
            write(output_unit,*)'Q_DE=',Q_DE,'Q_AD=',Q_AD
            write(output_unit,*)'Q_REG=',Q_REG
            write(output_unit,*)'Q_TOTIN=',Q_TOTIN ,'Q_EVAP=',Q_EVAP,'Q_CO=',Q_CO
            write(output_unit,*)'Q_LTTCSAT=',Q_LTTCSAT,'COP',COP(I),'SCE',SCE(I)
            write(output_unit,*)'=========================================================='
        endif
    end subroutine output_iteration
END SUBROUTINE COPCALC

You can then use the shell redirection to choose the output you want. For example
gfortran 1.f90
./a.out 2> deb_data

You can plot that file directly in gnuplot with
plot ‘deb_data’ with lines