forums.silverfrost.com

[Alvaro]

Hi all,
I have a problem, while running a old fortran 77 code in my Plato 3 the program it gives me the following error:
"attempt to call a subroutine as if were a real (kind=1) function"
EQ_OF_TIME1: at line 423[+0044]
TLinke: at line 220[+03de]
main: at line 169 [+159d]

I'm calling the soubrutine like:
"call TLinke(doy,tAira,gsra1,gsra2,gsra3,solelv,TL,xId)" ) (line 169)
While the subrutine line 220 says:
"eot=eq_of_time1(doy)" (doy is double precision)
While eq_of_time1 is a function, where in line 440 there is an expresion like:
"eq_of_time1=c1*dcos(ed) + c2*dcos(2.d0*ed)+..."

Do you know what am I doing wrong? I think the problem it is related with some kind of double precision numbers, but I am not sure at all.

Thanks anyway

[PaulLaidler]

You will need to provide more information.
Is eq_of_time1 defined as a function or a subroutine or
perhaps even a statement function?

Perhaps the most likely possibility is that it is defined as a function
and called first in your program as a subroutine and then as a function.

The first call will tell the compiler what to expect (a subroutine) and the second call will indicate an error.

[Alvaro]

Hi Paul.
The subroutine is:
"subroutine TLinke(doy,tAira,gsra1,gsra2,gsra3,
> solelv,TL,xId)

DOUBLE PRECISION doy,time,lat,lon,decl,omega,ed, r_cor
DOUBLE PRECISION, DIMENSION(24):: solelv,tt
PARAMETER (xI00=1368.)
PARAMETER (GMR=-3.41632e-5)
PARAMETER (H=148.)
PARAMETER (Lon = 6.23d0)
PARAMETER (Lat=37.43d0)
REAL xI0,tmp
REAL, DIMENSION(24):: tAira,gsra1,gsra2,gsra3
REAL, DIMENSION(24):: TL,p,xId,snh
parameter (pi=3.141592653589793d0)
parameter (rtd=180.d0/pi)
parameter (dtr=1.d0/rtd)

ed=doy*2.d0*pi/365.d0
r_cor=1+0.03344*cos(ed)
xI0=xI00*r_cor
write(*,*) "doy", doy," xI0 =",xI0
do i=1,24
eot=eq_of_time1(doy)
write(*,*) " eot = ", eot
if(i.EQ.1) then
tt(1)=dble(i)-1.5
else
tt(i)=dble(i)-1.5
endif
decl=declination(doy,tt(i))
write(*,*) "decl ", decl*rtd
omega=h_angle(tt(i),eot,lon)
write(*,*) "omega", omega
solelv(i)=sol_el(decl,omega,lat)
write(*,*) " solelv ", solelv(i)*rtd
snh(i)=sin(solelv(i))
write(*,*) "sinh = ", snh(i)
tmp=gsra2(i)/(xI0*snh(i)*0.84)
write(*,*) " gsra2 tmp = ", gsra2(i),tmp
if(solelv(i)*rtd.GT.10.) then

TL(i)=(-snh(i)/0.027)*LOG(gsra2(i)/(xI0*snh(i)*0.84))
else
TL(i)=0.
endif
write(*,*) " TL ", TL(i)
p(i)=EXP(GMR*H/tAira(i))
if(TL(i).GT.0.) then

xId(i)= xI0*EXP(-TL(i)/(0.9+9.4*snh(i)))
else
xId(i)=0.
endif
enddo
open(10,file='alv1.txt')
do i=1,24
write(10,50) i,tt(i),solelv(i)*rtd,snh(i),gsra2(i),TL(i),
> p(i),tAira(i),xId(i)
enddo
50 format(i3,1x,8f12.6,1x,f16.6)
end"

And the function:
"function eq_of_time(doy,time)
*.......................................................................
* compute the equation of time in radians, convert to ????? .
*.......................................................................
implicit double precision (a-h,o-z)
parameter (pi=3.141592653589793d0)
parameter (c0=0.000075d0)
parameter (c1=0.001868d0)
parameter (c2=-0.014615d0)
parameter (c4=-0.032077d0)
parameter (c5=-0.040849d0)

dn=doy-1.d0
ed=(dn+((time-12.d0)/24.d0))*2.d0*pi/365.d0
eq_of_time = c0+
> c1*dcos(ed)+c2*dcos(2.d0*ed)+
> c4*dsin(ed)+c5*dsin(2.d0*ed)
c eq_of_time=eq_of_time*(180.d0/pi)*24.d0*60.d0/360.d0
return
end"
Thanks in advance

[PaulLaidler]

These look OK. It is the way that the function is called in the main program that is important. It is presumably called in at least two places, the first time it is called as a subroutine when in fact it is a function.

If this does not explain it, we will need to see the main program in detail.

[Alvaro]

Hi Paul

I've been searching, but I don't find that the factor is define in two places. If you don't mind I place here the whole code. It is a code to Calculate the Direct Solar Radiation in a given place.

here it goes, it's a bit long:

"PROGRAM radCal3
c IMPLICIT NONE
******************************************************************************
* This program is developed to calculate the Direct Normal Solar Radiation *
******************************************************************************


**************************************************************************
* *
* Define the variables: *
* *
**************************************************************************

REAL, DIMENSION(721) :: rh,tAir,wd,dsr,gsr1,gsr2,gsr3!,ws
REAL, DIMENSION(24):: rha,tAira,wda,dsra,gsra1,gsra2,gsra3!,wsa
INTEGER, DIMENSION(721) :: mm,dd,hh,mn,ss
INTEGER, DIMENSION(13) :: nm
INTEGER :: i,j,k,nday,ierror,nl_max,dl!,nd,nmonth,nl
INTEGER :: n_avg
REAL :: T,wdd,GSR,RHH,tmp,p,xId,Lat,Lon
CHARACTER(110) :: ofileR,ofileF,ifileR,ifileF
CHARACTER(110) :: rootiR,rootiF,rootoR,rootoF
CHARACTER(140) :: input_fileR,output_fileR,input_fileF,
> output_fileF
CHARACTER(50) :: ch,ichars
c CHARACTER :: ch1*1
DIMENSION :: ch(1000)
CHARACTER(:: name
PARAMETER (Lon = 6.23)
PARAMETER (Lat=37.43)
PARAMETER (nl_max=10000)
PARAMETER (n_avg=30)
PARAMETER (dl=24)
DIMENSION hour(dl)
DATA nm / 0,31,59,90,120,151,181,212,243,273,304,334,365 /
DIMENSION p(dl),xId(dl),TL(dl)
DOUBLE PRECISION doy
DOUBLE PRECISION, DIMENSION(24):: solelv

**************************************************************************
* *
* INPUT AND OUTPUT FILES (Defined for September(09) 2008) *
* *
* P = day number of the month *
* ifile = defined as HeatingAndTemperature_2008-09- *
**************************************************************************


! The root of the input file is:

rootiR = "C:\Documents and Settings\MAAMR\Mis documentos\Datos Ob
>servacion\July\E. Met. PS10 01-07-08.xml"

c rootoR = "C:\Documents and Settings\MAAMR\Mis documentos\Datos Ob
c >servacion\Output\July\"

c !ifile = 'E. Met. PS10 01-07-08.xml'

input_fileR = rootiR

c !output_fileR = ofile

c PRINT *, input_fileR
OPEN (UNIT=10,FILE=input_fileR,STATUS="OLD",
>ACTION="READ",POSITION="REWIND",IOSTAT=ierror)
IF (ierror /=0) then
print *,ierror, "NOT ABLE TO OPEN INPUT FILE R"
STOP 123456
END IF

! Read the file, and create a 1 dimensional array with the variables from Goblar Solar Radiation thrid sensor

DO i=1,64
READ(10,*)
END DO
write(*,*) "Read header"

DO k=1,nl_max
c write(*,*) "k = ",k
read(10, 312) name
c write(*,*) "name ", name
if (name.NE."</Table>") then
READ(10,310) mm(k),dd(k),hh(k),mn(k),ss(k)
c write(*,*

[Alvaro]

CONTINUE FROM:( if (name.NE."</Table>") then)
READ(10,310) mm(k),dd(k),hh(k),mn(k),ss(k)
c write(*,*) "month ",mm(k)," day ",dd(k)
READ(10,309) ichars
rh(k)=RHH(ichars)
c write(*,*) "k:", k, " RH:",rh(k)
READ(10,309) ichars
tAir(k)=T(ichars)
c write(*,*) "k:", k, " T",tAir(k)
READ (10,309) ichars
wd(k)=wdd(ichars)
c write(*,*) "k:", k, " wd",wd(k)
read(10,*) !nchars,irchars
!ws(k)=wss(ncharts,ichars)
read(10,309) dsr(k)
c write(*,*) "k:", k, " dsr",dsr(k)
read(10,309) ichars
gsr1(k)=gsr(ichars)
c write(*,*) "k:", k, " gsr1",gsr1(k)
read(10,309) ichars
gsr2(k)=gsr(ichars)
c write(*,*) "k:", k, " gsr2",gsr2(k)
read(10,311) gsr3(k)
c write(*,*) "k:", k, " gsr3",gsr3(k)
read(10,*)
c PRINT *, tAir(k)
else
nl=k-1
goto 1200
endif
enddo


309 FORMAT(50x,(A))
310 FORMAT(57x,i2,1x,i2,1x,i2,1x,i2,1x,i2)
311 FORMAT(50x,f6.0)
312 FORMAT(2x,(a))
1200 CONTINUE
doy=dble(nm(mm(1))+dd(1))
write(*,*) " DOY = ",doy
write(*,*) "Input read correctly!"

do j=0,23
t1=0.
t2=0.
t3=0.
c t4=0.
t5=0.
t6=0.
t7=0.
t8=0.
do i=j*30+1,j*30+n_avg
t1=t1+rh(i)
t2=t2+tAir(i)
t3=t3+wd(i)
C t4=t4+ws(i)
t5=t5+dsr(i)
t6=t6+gsr1(i)
t7=t7+gsr2(i)
t8=t8+gsr3(i)
c write(*,*) "i j tmp", i,j,tmp
end do
rha(j+1)=t1/real(n_avg)
tAira(j+1)=t2/n_avg
wda(j+1)=t3/n_avg
c wsa(j+1)=t4/n_avg
dsra(j+1)=t5/n_avg
gsra1(j+1)=t6/n_avg
gsra2(j+1)=t7/n_avg
gsra3(j+1)=t8/n_avg
c write(*,*) "j+1 ",j+1, "GRSA ", gsra1(j+1)
enddo


PRINT *, "The hourly average values for global radiation are"
write(*,*) "h rh T wd dsr gsr1 gsr2 gsr3"
do i=1,24
write(*,8 i,rha(i),TAir(i),wd(i),dsr(i),gsra1(i),gsra2(i),
>gsra3(i)
enddo
88 format(i2,1x,f3.0,1x,f5.1,1x,f5.0,1x,f5.0,3(f6.0,1x))

call TLinke(doy,tAira,gsra1,gsra2,gsra3,solelv,TL,
>xId)



close(Unit=10)

end program radCal3




subroutine TLinke(doy,tAira,gsra1,gsra2,gsra3,
> solelv,TL,xId)
************************************************************************************
* LINKE'S TURBIDITY FACTOR *
* The Linke's turbidity factor is an extintion factor, is the number of pure *
* dry air masses which results in the same extinction than the real atmosphere. *
* The extintion factor would be expressed in terms of Surface Global Radiation *
* and the solar position. *
* The current code will calculate the Turbidity Factor for each hour, and later *
* does a daily mean average for each Turbidity Factor. *
* *
* The actual code calculates the Linke's factor from the observations, but it *
* has to be obtain from the forecasted values of SGR. *
************************************************************************************
DOUBLE PRECISION doy,time,lat,lon,decl,omega,ed, r_cor
DOUBLE PRECISION, DIMENSION(24):: solelv,tt
c PARAMETER (PI=3.14159265358979323)
PARAMETER (xI00=1368.)
PARAMETER (GMR=-3.41632e-5)
PARAMETER (H=148.)
PARAMETER (Lon = 6.23d0)
PARAMETER (Lat=37.43d0)
REAL xI0,tmp
REAL, DIMENSION(24):: tAira,gsra1,gsra2,gsra3
REA

[Alvaro]

CONTINUE: (REAL, DIMENSION(24):: tAira,gsra1,gsra2,gsra3)
REAL, DIMENSION(24):: TL,p,xId,snh
parameter (pi=3.141592653589793d0)
parameter (rtd=180.d0/pi)
parameter (dtr=1.d0/rtd)

ed=doy*2.d0*pi/365.d0
r_cor=1+0.03344*cos(ed)
xI0=xI00*r_cor
write(*,*) "doy", doy," xI0 =",xI0
do i=1,24
eot=eq_of_time1(doy)
write(*,*) " eot = ", eot
if(i.EQ.1) then
tt(1)=dble(i)-1.5
else
tt(i)=dble(i)-1.5
endif
decl=declination(doy,tt(i))
write(*,*) "decl ", decl*rtd
omega=h_angle(tt(i),eot,lon)
write(*,*) "omega", omega
solelv(i)=sol_el(decl,omega,lat)
write(*,*) " solelv ", solelv(i)*rtd
snh(i)=sin(solelv(i))
write(*,*) "sinh = ", snh(i)
tmp=gsra2(i)/(xI0*snh(i)*0.84)
write(*,*) " gsra2 tmp = ", gsra2(i),tmp
if(solelv(i)*rtd.GT.10.) then

TL(i)=(-snh(i)/0.027)*LOG(gsra2(i)/(xI0*snh(i)*0.84))
else
TL(i)=0.
endif
write(*,*) " TL ", TL(i)
p(i)=EXP(GMR*H/tAira(i))
if(TL(i).GT.0.) then
xId(i)= xI0*EXP(-TL(i)/(0.9+9.4*snh(i)))
else
xId(i)=0.
endif
enddo
open(10,file='alv1.txt')
do i=1,24
write(10,50) i,tt(i),solelv(i)*rtd,snh(i),gsra2(i),TL(i),
> p(i),tAira(i),xId(i)
enddo
50 format(i3,1x,8f12.6,1x,f16.6)
end



function T(ichars)

INTEGER n
CHARACTER ichars*50, ch1*1,ch2*1,ch3*1,ch4*1,ch5*1

read(ichars(1:1),'(a)') ch1
read(ichars(2:2),'(a)') ch2
read(ichars(3:3),'(a)') ch3
read(ichars(4:4),'(a)') ch4
read(ichars(5:5),'(a)') ch5

if(ch1.EQ.'-') then
if(ch3.EQ.'.') then
read(ichars(1:4),'(f4.1)') T
elseif(ch3.EQ.'<') then
read(ichars(1:2),'(i2)') n
T=real(n)
elseif(ch4.eq.'.') then
read(ichars(1:5),'(f5.1)') T
elseif(ch4.EQ.'<') then
read(ichars(1:3),'(i3)') n
T=real(n)
endif
elseif(ch1.NE.'-') then
if(ch2.EQ.'.') then
read(ichars(1:3),'(f3.1)') T
elseif(ch2.EQ.'<') then
read(ichars(1:1),'(i1)') n
T=real(n)
elseif(ch3.EQ.'.') then
read(ichars(1:4),'(f4.1)') T
elseif(ch3.EQ.'<') then
read(ichars(1:2),'(i2)') n
T=real(n)
endif
endif
return
end

function WDD(ichars)
INTEGER n
CHARACTER ichars*50, ch1*1,ch2*1,ch3*1,ch4*1,ch5*1

read(ichars(1:1),'(a)') ch1
read(ichars(2:2),'(a)') ch2
read(ichars(3:3),'(a)') ch3
read(ichars(4:4),'(a)') ch4
read(ichars(5:5),'(a)') ch5

if(ch2.EQ.'<') then
read(ichars(1:1),'(i1)') n
wdd=real(n)
elseif(ch2.EQ.'.') then
read(ichars(1:3),'(f3.1)') wdd
endif
if(ch3.EQ.'<') then
read(ichars(1:2),'(i2)') n
wdd=real(n)
elseif(ch3.EQ.'.') then
read(ichars(1:4),'(f4.1)') wdd
endif
if(ch4.EQ.'<') then
read(ichars(1:3),'(i3)') n
wdd=real(n)
elseif(ch4.EQ.'.') then
read(ichars(1:5),'(f5.1)') wdd
endif
return
end

function GSR(ichars)
INTEGER n
CHARACTER ichars*50, ch1*1,ch2*1,ch3*1,ch4*1,ch5*1

read(ichars(1:1),'(a)') ch1
read(ichars(2:2),'(a)') ch2
read(ichars(3:3),'(a)') ch3
read(ichars(4:4),'(a)') ch4
read(ichars(5:5),'(a)') ch5

if(ch2.EQ.'<') then
read(ichars(1:1),'(i1)') n
gsr=real(n)
elseif(ch2.EQ.'.') then
read(ichars(1:3),'(f3.1)') gsr
endif
if(ch3.EQ.'<') then
read(ichars(1:2),'(i2)') n
gsr=real(n)
elseif(ch3.EQ.'.') then
read(ichars(1:4),'(f4.1)') gsr
endif
if(ch4.EQ.'<') then
read(ichars(1:3),'(i3)') n
gsr=real(n)
elseif(ch4.EQ.'.') then
read(ichars(1:5),'(f5.1)') gsr
endif
if(ch5.EQ.'<') then
read(ichars(1:

[Alvaro]

The are several fuctions which are similar defined in the code, it continue later like:
function declination(doy,tt)
*.......................................................................
* declination in radians .
* enter day number and compute solar declination .
*.......................................................................
implicit double precision (a-h,o-z)
parameter (pi=3.141592653589793d0)
parameter (rtd=180.d0/pi)
parameter (dtr=1.d0/rtd)
parameter (c1=0.006918d0)
parameter (c2=-0.399912d0)
parameter (c3=-0.006758d0)
parameter (c4=-0.002697d0)
parameter (c5=0.070257d0)
parameter (c6=0.000907d0)
parameter (c7=0.001480d0)

dn=doy-1.d0
ed=(dn+((tt-12.d0)/24.d0))*2.d0*pi/365.d0
declination = c1 +
> c2*dcos(ed) + c3*dcos(2.d0*ed) + c4*dcos(3.d0*ed) +
> c5*dsin(ed) + c6*dsin(2.d0*ed) + c7*dsin(3.d0*ed)
d_degr = rtd*declination
frac = dabs(d_degr) - dint( dabs(d_degr) )
c minutes = idnint(frac*60.d0)
return
end


function eq_of_time1(doy)
*.......................................................................
* compute the equation of time in hours!! convert to ????? .
*.......................................................................
implicit double precision (a-h,o-z)
parameter (pi=3.141592653589793d0)
parameter (c1=0.0072d0)
parameter (c2=-0.0528d0)
parameter (c3=-0.0012d0)
parameter (c4=-0.1229d0)
parameter (c5=-0.1565d0)
parameter (c6=-0.0041d0)

ed=doy*2.d0*pi/366.d0
eq_of_time1=c1*dcos(ed) + c2*dcos(2.d0*ed) + c3*dcos(3.d0*ed) +
> c4*dsin(ed) + c5*dsin(2.d0*ed) + c6*dsin(3.d0*ed)
c eq_of_time1=eq_of_time1*60.d0
return
end

function h_angle(tt,e,xlambda)
*.......................................................................
* Omega = the hour angle, defined as TST + 12 h, .
* time = GMT, (for MEZT a correction of -2.d0 is needed), .
* xlambda/15 = correction for the longitude of the observer, .
* e = equation of time, -12.0 h comes from omega = TST + 12.0 h .
* .
* Omega is computed in radians .
*.......................................................................
implicit double precision (a-h,o-z)
parameter (pi=3.141592653589793d0)
parameter (rtd=180.d0/pi)
parameter (dtr=1.d0/rtd)

h_angle=(360.d0/24.d0)*dtr*(tt-(xlambda/15.d0)+e-12.d0)
return
end




function sol_el(decl,omega,xlatitude)
*.......................................................................
* The equation to obtain the solar elevation (solel) is from .
* A.C. Velds 6.1.3, p 132 .
* The solar elevation sol_el is compted in radians .
*.......................................................................
implicit double precision (a-h,o-z)
parameter (pi =3.141592653589793d0 )
parameter (rtd=180.d0/pi)
parameter (dtr=1.d0/rtd)

xlat=xlatitude*dtr
sinh=dsin(decl)*dsin(xlat)+dcos(decl)*dcos(xlat)*dcos(omega)
sol_el=dasin(sinh)
return
end

[Alvaro]

and finally:
function sol_az(tt,decl,omega,solel,xlatitude)
*.......................................................................
* Computation af solar azimuth in radians .
*.......................................................................
implicit double precision (a-h,o-z)
parameter (pi=3.141592653589793d0 )
parameter (rtd=180.d0/pi)
parameter (dtr=1.d0/rtd)

xlat=xlatitude*dtr

cosa=(-dsin(decl)*dcos(xlat) +
> dcos(decl)*sin(xlat)*dcos(omega))/dcos(solel)
azc=dacos(cosa)
azc_degr=azc*rtd

sina=dcos(decl)*dsin(omega)/dcos(solel)
azs=dasin(sinA)
azs_degr=azs*rtd

if(tt.lt.12.d0) then
sol_az = 180.d0 + azs_degr
else
sol_az = 180.d0 - azs_degr
endif
return
end


Thanks for your help

[PaulLaidler]

I am not sure about the details but I think EQ_OF_TIME1 is called as single precision and defined as double precision.

Also you have 365 and 366 instead of 360 for radian conversion.

[JohnHorspool]

Aren't 365 and 366 used as days in a year and a leap year?

[Alvaro]

So, how can I call a function in double precision?

[PaulLaidler]

Declare the function name as double precision in the calling program.

[Alvaro]

Thanks, I declared the functions as reals and now it works properly.
Thanks!

[PaulLaidler]

You will also need to fix the incorrect values for conversion to radians.

Angle in radians = 2*pi/360*(Angle in degrees)