NaN results

[zamaniM]

Hi.
I wrote this code but all results except of k_1 were NaN.Do you know what's the problem?

PROGRAM GRS_MECHANISM

IMPLICIT NONE

! The purpose Is Calculating K_1,K_2&K_4 from subroutines
!*******************************************************************VARIABLES*******************************************************
!VARIABLES OF MAIN PROGRAM

REAL :: K_1(24,6),K_2(24,6),K_4(24,6),K_5(24,6),K_6(24,6),K_7(24,6)


!VARIABLES OF NO2PHOTOLYSISRATE SUBROUTINE

REAL :: j,i,fi_radian,D_l,fi
REAL,PARAMETER :: year=2010
REAL,PARAMETER :: pi=3.141592654
REAL :: D(24,6),T_celcius(24,6),f_T(24,6),SWF(24,6),teta_s(24,6),ActinicFlux(24,6),T_kelvin(24,6),&
&K_3(24,6),time_hour(24),t_s(24),H_a(24), H_a_radian(24),teta_s_radian(24,6),cosin_tetas(24,6),&
&D_J(6),N_J_D(6),cumulative_ActinicFlux(6:24,6),Epsilon_o_b(6),Epsilon_o_b_radian(6),Delta(6),&
&Delta_radian(6),landa_e_c(6),landa_e_c_radian(6),g_M(6),g_M_radian(6),L_M(6),L_M_radian(6),&
&ozone(24,6),NO(24,6),NO_2(24,6),ox_difference(2:24,6),OX(24,6),d_cumulative_ActinicFlux(7:24,6),tangant(7:24,6)


!VARIABLES OF a_ROC2 SUBROUTINE

REAL ::Total_THC_MobileEmission(8),Daily_THC_mobileEmission(8),HC_Species(6),HC_Species_WeightingFactor(6),HC_Species_reactivityfactor(6),HC_Species_Emission(6,8),&
&Total_HC_SpeciesEmission(6,8),Daily_Benzene_Emission(99),WeakDay_Benzene_Emission(84),Weakend_Benzene_Emission(15),a(6),b(6)
REAL ::Index_WeakdayWeakend,Total_DailyBenzene,Total_WeakdayBenzene,Total_WeakendBenzene,a_ROC,SUM_DailyBenzene,SUM_WeakdayBenzene,SUM_WeakendBenzene
!***********************************************************************************************************************************

CALL NO2photolysisrate(j,i,fi_radian,D_l,fi,year,pi,D,T_celcius,f_T,SWF,teta_s,ActinicFlux,T_kelvin,K_3,time_hour,t_s,H_a,H_a_radian,teta_s_radian,cosin_tetas,&
&D_J,N_J_D,cumulative_ActinicFlux,Epsilon_o_b,Epsilon_o_b_radian,Delta,Delta_radian,landa_e_c,landa_e_c_radian,g_M,g_M_radian,L_M,L_M_radian,&
&ozone,NO,NO_2,ox_difference,ox,d_cumulative_ActinicFlux,tangant)


CALL a_ROC2(Total_THC_MobileEmission,Daily_THC_mobileEmission,HC_Species,HC_Species_WeightingFactor,HC_Species_reactivityfactor,HC_Species_Emission,&
&Total_HC_SpeciesEmission,Daily_Benzene_Emission,WeakDay_Benzene_Emission,Weakend_Benzene_Emission,a,b,i,j,Index_WeakdayWeakend,&
&Total_DailyBenzene,Total_WeakdayBenzene,Total_WeakendBenzene,a_ROC,SUM_DailyBenzene,SUM_WeakdayBenzene,SUM_WeakendBenzene)

!***********************************************************************************************************************************


OPEN(1,file='output_k_3.xml',status='old') !k_3 is calculated from subroutine(NO2photolysisrate)
DO I=1,24
READ(1,*)(K_3(I,J),J=1,6)
END DO


OPEN(2,file='output_f_T.xml',status='old') !f_T is calculated from subroutine(NO2photolysisrate)
DO I=1,24
READ(2,*)(f_T(I,J),J=1,6)
END DO

OPEN(3,file='a_ROC.txt',status='old') !a_ROC is calculated from subroutine(a_ROC2)
READ(3,*)(a_ROC)

OPEN(4,file='T_kelvin.txt',status='old') !T_kelvin is calculated from subroutine(NO2photolysisrate)
DO I=1,24
READ(4,*)(T_kelvin(I,J),J=1,6)
END DO

DO I=1,24
DO J=1,6
K_1(I,J)=a_ROC*K_3(I,J)*f_T(I,J)
K_2(I,J)=(3.79*(10**-12)*EXP(242/T_kelvin(I,J)))*(5482.0/(3.79*(10**-12)*60.0))
K_4(I,J)=1.7*(10**-12)*EXP(-1370/T_kelvin(I,J))*(5482.0/(3.79*(10**-12)*60.0))
K_5(I,J)=6.76*(10**-12)*(5482.0/(3.79*(10**-12)*60.0))
K_6(I,J)=10**-13*(5482.0/(3.79*(10**-12)*60.0))
K_7(I,J)=K_6(I,J)
END DO
END DO


OPEN(5,file='K_1.xml')
DO I=1,24
WRITE(5,"(6(1x,1pe14.7))") (K_1(I,J),J=1,6)
END DO

OPEN(6,file='K_2.xml')
DO I=1,24
WRITE(6,"(6(1x,1pe14.7))") (K_2(I,J),J=1,6)
END DO

OPEN(7,file='K_4.xml')
DO I=1,24
WRITE(7,"(6(1x,1pe14.7))") (K_4(I,J),J=1,6)
END DO

OPEN(8,file='K_5.xml')
DO I=1,24
WRITE(8,"(6(1x,1pe14.7))") (K_5(I,J),J=1,6)
END DO

OPEN(9,file='K_6.xml')
DO I=1,24
WRITE(9,"(6(1x,1pe14.7))") (K_6(I,J),J=1,6)
END DO



END PROGRAM GRS_MECHANISM

SUBROUTINE NO2photolysisrate(j,i,fi_radian,D_l,fi,year,pi,D,T_celcius,f_T,SWF,teta_s,ActinicFlux,T_kelvin,K_3,time_hour,t_s,H_a,H_a_radian,teta_s_radian,cosin_tetas,&
&D_J,N_J_D,cumulative_ActinicFlux,Epsilon_o_b,Epsilon_o_b_radian,Delta,Delta_radian,landa_e_c,landa_e_c_radian,g_M,g_M_radian,L_M,L_M_radian,&
&ozone,NO,NO_2,ox_difference,ox,d_cumulative_ActinicFlux,tangant)
END SUBROUTINE NO2photolysisrate

SUBROUTINE a_ROC2(Total_THC_MobileEmission,Daily_THC_mobileEmission,HC_Species,HC_Species_WeightingFactor,HC_Species_reactivityfactor,HC_Species_Emission,&
&Total_HC_SpeciesEmission,Daily_Benzene_Emission,WeakDay_Benzene_Emission,Weakend_Benzene_Emission,a,b,i,j,Index_WeakdayWeakend,&
&Total_DailyBenzene,Total_WeakdayBenzene,Total_WeakendBenzene,a_ROC,SUM_DailyBenzene,SUM_WeakdayBenzene,SUM_WeakendBenzene)

END SUBROUTINE a_ROC2

 

[gummibaer]

Occurs to me you could do a lot in simplifying your numerical formulae here. Most of it seems to be constant1 * exp (variable)
Maybe you could get a clue on the magnitude of the values and see if you have overflows or underflows.

K5 and K6 seem to be constant. Why compute the same thing over and over again ?

And then it might be a good idea to debug your prog and see the values of T_kelvin as they go into your program. If they are very small or very big this EXP function could get out of hand.

Norbert