I am struggling to replicate the following fortran 77 subroutine into C# method. The main problem is that I am trying to get rid of the fortran's goto
statements.
The fortran to C# conversion looks fine, still the goto
statements remained. Is there a way all of them can be replaced with conditionals or some other ways? Thank you.
Here is the fortran subroutine:
C PROGRAM TEMP
subroutine TEMP (acl,adu,aeff,c,cair,cb,cbare,
+ cclo,count1,csum,di,ed,emcl,emsk,enbal,
+ enbal2,ere,erel,esw,eswdif,eswphy,eswpot,
+ evap,facl,fcl,fec,feff,food,h,hc,he,ht,htcl,icl,j,
+ mbody,p,po,r1,r2,rbare,rcl,
+ rclo,rclo2,rdcl,rdsk,rob,rsum,sex,sigm,sw,swf,swm,
+ ta,tbody,tcl,tcore,tmrt,tsk,v,vb,vb1,vb2,
+ vpa,vpts,wetsk,wd,wr,ws,wsum,xx)
real acl,adu,aeff,c(0:10),cair,cb,cbare,
+ cclo,csum,di,ed,emcl,emsk,enbal,
+ enbal2,ere,erel,esw,eswdif,eswphy,eswpot,
+ evap,facl,fcl,fec,feff,food,h,hc,he,ht,htcl,icl,
+ mbody,p,po,r1,r2,rbare,rcl,
+ rclo,rclo2,rdcl,rdsk,rob,rsum,sigm,sw,swf,swm,
+ ta,tbody,tcl,tcore(1:7),tmrt,tsk,v,vb,vb1,vb2,
+ vpa,vpts,wetsk,wd,wr,ws,wsum,xx
integer count1,count3,j,sex
wetsk = 0.
adu = 0.203 * mbody ** 0.425 * ht ** 0.725
hc = 2.67 + ( 6.5 * v ** 0.67)
hc = hc * (p /po) ** 0.55
feff = 0.725
C
facl = (- 2.36 + 173.51 * icl - 100.76 * icl * icl + 19.28
+ * (icl ** 3.)) / 100.
C
if (facl .gt.1.) facl = 1.
rcl = (icl/6.45)/facl
if (icl.ge.2.) y = 1.
if ((icl .gt. 0.6) .and. (icl .lt. 2.)) y = (ht - 0.2) / ht
if ((icl .le. 0.6) .and. (icl .gt. 0.3)) y = 0.5
if ((icl .le. 0.3) .and. (icl .gt. 0.)) y = 0.1
r2 = adu * (fcl - 1. + facl) / (2. * 3.14 * ht * y)
r1 = facl * adu / (2. * 3.14 * ht * y)
di = r2 - r1
C TEMPERATURE
do 90 j = 1,7
tsk = 34.
count1 = 0
tcl = (ta + tmrt + tsk) / 3.
count3 = 1
enbal2 = 0.
20 acl = adu * facl + adu * (fcl - 1.)
rclo2 = emcl*sigm *((tcl+273.2)** 4.-(tmrt+273.2)** 4.)*feff
htcl = 6.28 * ht * y * di / (rcl * alog(r2/r1) * acl)
tsk = 1. / htcl * (hc * (tcl - ta) + rclo2) + tcl
C RADIATION
aeff = adu * feff
rbare = aeff * (1.-facl) * emsk * sigm *
+ ((tmrt + 273.2) ** 4. - (tsk + 273.2) ** 4.)
rclo = feff * acl * emcl * sigm *
+ ((tmrt + 273.2) ** 4. - (tcl + 273.2) ** 4.)
rsum = rbare + rclo
C CONVECT
cbare = hc * (ta - tsk) * adu * (1. - facl)
cclo = hc * (ta - tcl ) * acl
csum = cbare + cclo
C CORE
c(0) = h + ere
c(1) = adu * rob * cb
c(2) = 18. - 0.5 * tsk
c(3) = 5.28 * adu * c(2)
c(4) = 0.0208 * c(1)
c(5) = 0.76075 * c(1)
c(6) = c(3) - c(5) - tsk * c(4)
c(7) = - c(0) * c(2) - tsk * c(3) + tsk * c(5)
c(8) = c(6) * c(6) - 4. * c(4) * c(7)
c(9) = 5.28 * adu - c(5) - c(4) * tsk
c(10) = c(9) * c(9) - 4. * c(4) *
+ (c(5) * tsk - c(0) - 5.28 * adu * tsk)
C
if (tsk.eq.36.) tsk=36.01
tcore(7) = c(0) / (5.28 * adu + c(1) * 6.3 / 3600.) + tsk
tcore(3) = c(0) / (5.28 * adu + (c(1) * 6.3 / 3600.) /
+ (1 + 0.5 * (34. -tsk))) + tsk
if (c(10) .lt. 0.) goto 22
tcore(6) = (- c(9) - c(10) ** 0.5) / (2. * c(4))
tcore(1) = (- c(9) + c(10) ** 0.5) / (2. * c(4))
22 if (c(8) .lt. 0.) goto 24
tcore(2) = (- c(6) + abs(c(8)) ** 0.5) / (2. * c(4))
tcore(5) = (- c(6) - abs(c(8)) ** 0.5) / (2. * c(4))
24 tcore(4) = c(0) / (5.28 * adu + c(1) * 1. / 40.) + tsk
C TRANSPARENCE
tbody = 0.1 * tsk + 0.9 * tcore (j)
swm = 304.94 * (tbody - 36.6) * adu / 3600000.
vpts = 6.11 * 10. ** (7.45 * tsk / (235. + tsk))
if (tbody .le. 36.6) swm = 0.
swf = 0.7 * swm
if(sex .eq. 1) sw = swm
if(sex .eq. 2) sw = swf
eswphy = - sw * evap
he = 0.633 * hc / (p * cair)
fec = 1. / (1. + 0.92 * hc * rcl)
eswpot = he * (vpa - vpts) * adu * evap * fec
wetsk = eswphy / eswpot
if (wetsk .gt. 1.) wetsk = 1.
eswdif = eswphy - eswpot
if (eswdif .le. 0.) esw = eswpot
if (eswdif .gt. 0.) esw = eswphy
if (esw .gt. 0.) esw = 0.
C DIFFERENCE
rdsk = 0.79 * 10. ** 7.
rdcl = 0.
ed = evap / (rdsk + rdcl) * adu * (1 - wetsk) * (vpa-vpts)
C VB
vb1 = 34. - tsk
vb2 = tcore(j) - 36.6
if (vb2 .lt.0.) vb2 = 0.
if (vb1 .lt.0.) vb1 = 0.
vb = (6.3 + 75. * (vb2)) / (1. + 0.5 * vb1)
C BALANCE
enbal = h + ed + ere + esw + csum + rsum + food
C COVER
if (count1 .eq.0) xx = 1.
if (count1 .eq.1) xx = 0.1
if (count1 .eq.2) xx = 0.01
if (count1 .eq.3) xx = 0.001
if (enbal .gt. 0.) tcl = tcl + xx
if (enbal .lt. 0.) tcl = tcl - xx
if ((enbal .le. 0.) .and. (enbal2 .gt. 0.)) goto 30
if ((enbal .ge. 0.) .and. (enbal2 .lt. 0.)) goto 30
enbal2 = enbal
count3 = count3 + 1
C
if (count3 .gt. 200) goto 30
goto 20
30 if ((count1 .eq.0.).or.(count1.eq.1.).or.(count1.eq.2.)) then
count1 = count1 + 1.
enbal2 = 0.
goto 20
end if
C
if (count1 .eq. 3.) then
C
if ((j .eq. 2) .or. (j .eq. 5)) goto 40
if ((j .eq. 6) .or. (j .eq. 1)) goto 50
if (j .eq. 3) goto 60
if (j .eq. 7) goto 70
if (j .eq. 4) goto 80
end if
40 if (c(8) .lt. 0.) goto 90
if ((tcore(j) .ge. 36.6) .and. (tsk .le. 34.050)) goto 80
goto 90
50 if (c(10) .lt. 0. ) goto 90
if ((tcore(j) .ge. 36.6) .and. (tsk .gt. 33.850)) goto 80
goto 90
60 if ((tcore(j) .lt. 36.6) .and. (tsk .le. 34.000)) goto 80
goto 90
70 if ((tcore(j) .lt. 36.6) .and. (tsk .gt. 34.000)) goto 80
goto 90
80 if ((j .ne. 4) .and. (vb .ge. 91.)) goto 90
if ((j. eq. 4) .and. (vb .lt. 89.)) goto 90
if (vb .gt. 90.) vb = 90.
C LOSSES
ws = sw * 3600. * 1000.
if (ws .gt.2000.) ws = 2000.
wd = ed / evap * 3600. * (-1000.)
wr = erel / evap * 3600. * (-1000.)
wsum = ws + wr + wd
goto 100
90 continue
100 return
end
And here is the C# attempt to replicate the upper fortran subroutine:
using System;
public static class GlobalMembers_TEMP
{
//C PROGRAM TEMP
// --------------------------------------------
public static void TEMP(acl, adu, aeff, c, cair, cb, cbare, cclo, count1, csum, di, ed , emcl , emsk , enbal , enbal2 , ere , erel , esw , eswdif , eswphy , eswpot , evap , facl , fcl , fec , feff , food , h , hc , he , ht , htcl , icl , j , mbody , p , po , r1 , r2 , rbare , rcl , rclo , rclo2 , rdcl , rdsk , rob , rsum , sex , sigm , sw , swf , swm , ta , tbody , tcl, tcore , tmrt , tsk , v , vb , vb1 , vb2 , vpa , vpts , wetsk , wd , wr , ws , wsum , xx)
{
float acl;
float adu;
float aeff;
float[] c = new float[11]; // ???
float cair;
float cb;
float cbare;
float cclo;
float csum;
float di;
float ed;
float emcl;
float emsk;
float enbal;
float enbal2;
float ere;
float erel;
float esw;
float eswdif;
float eswphy;
float eswpot;
float evap;
float facl;
float fcl;
float fec;
float feff;
float food;
float h;
float hc;
float he;
float ht;
float htcl;
float icl;
float mbody;
float p;
float po;
float r1;
float r2;
float rbare;
float rcl;
float rclo;
float rclo2;
float rdcl;
float rdsk;
float rob;
float rsum;
float sigm;
float sw;
float swf;
float swm;
float ta;
float tbody;
float tcl;
float[] tcore = new float[7]; // ???
float tmrt;
float tsk;
float v;
float vb;
float vb1;
float vb2;
float vpa;
float vpts;
float wetsk;
float wd;
float wr;
float ws;
float wsum;
float xx;
int count1;
int count3;
int j;
int sex;
wetsk = 0.0;
adu = 0.203 * mbody * *0.425 * ht * *0.725;
hc = 2.67 + (6.5 * v * *0.67);
hc = hc * (p / po) * *0.55;
feff = 0.725;
//C rcl = icl / 6.45
facl = (-2.36 + 173.51 * icl - 100.76 * icl * icl + 19.28 * (icl * *3.0)) / 100.0;
//C
if (facl > 1.0)
{
facl = 1.0;
}
rcl = (icl / 6.45) / facl;
if (icl >= 2.0)
{
y = 1.0;
}
if ((icl > 0.6)) && ((icl < 2.0))
{
y = (ht - 0.2) / ht;
}
if ((icl <= 0.6)) && ((icl > 0.3))
{
y = 0.5;
}
if ((icl <= 0.3)) && ((icl > 0.0))
{
y = 0.1;
}
r2 = adu * (fcl - 1.0 + facl) / (2.0 * 3.14 * ht * y);
r1 = facl * adu / (2.0 * 3.14 * ht * y);
di = r2 - r1;
//C TEMPERATURE
for (90 j = 1, 7)
{
tsk = 34.0;
count1 = 0;
tcl = (ta + tmrt + tsk) / 3.0;
count3 = 1;
enbal2 = 0.0;
g20:
acl = adu * facl + adu * (fcl - 1.0);
rclo2 = emcl * sigm * ((tcl + 273.2) * *4.0 - (tmrt + 273.2) * *4.0) * feff;
htcl = 6.28 * ht * y * di / (rcl * alog(r2 / r1) * acl);
tsk = 1.0 / htcl * (hc * (tcl - ta) + rclo2) + tcl;
//C RADIATION
aeff = adu * feff;
rbare = aeff * (1.0 - facl) * emsk * sigm * ((tmrt + 273.2) * *4.0 - (tsk + 273.2) * *4.0);
rclo = feff * acl * emcl * sigm * ((tmrt + 273.2) * *4.0 - (tcl + 273.2) * *4.0);
rsum = rbare + rclo;
//C CONVECT
cbare = hc * (ta - tsk) * adu * (1.0 - facl);
cclo = hc * (ta - tcl) * acl;
csum = cbare + cclo;
//C CORE
c[0] = h + ere;
c[1] = adu * rob * cb;
c[2] = 18.0 - 0.5 * tsk;
c[3] = 5.28 * adu * c[2];
c[4] = 0.0208 * c[1];
c[5] = 0.76075 * c[1];
c[6] = c[3] - c[5] - tsk * c[4];
c[7] = -c[0] * c[2] - tsk * c[3] + tsk * c[5];
c[8] = c[6] * c[6] - 4.0 * c[4] * c[7];
c[9] = 5.28 * adu - c[5] - c[4] * tsk;
c[10] = c[9] * c[9] - 4.0 * c[4] * (c[5] * tsk - c[0] - 5.28 * adu * tsk);
//C
if (tsk == 36.0)
{
tsk = 36.01;
}
tcore[7] = c[0] / (5.28 * adu + c[1] * 6.3 / 3600.0) + tsk;
tcore[3] = c[0] / (5.28 * adu + (c[1] * 6.3 / 3600.0) / (1 + 0.5 * (34.0 - tsk))) + tsk;
if (c[10] < 0.0)
{
goto g22;
}
tcore[6] = (-c[9] - c[10] * *0.5) / (2.0 * c[4]);
tcore[1] = (-c[9] + c[10] * *0.5) / (2.0 * c[4]);
g22:
if (c[8] < 0.0)
{
goto g24;
}
tcore[2] = (-c[6] + Math.Abs(c[8]) * *0.5) / (2.0 * c[4]);
tcore[5] = (-c[6] - Math.Abs(c[8]) * *0.5) / (2.0 * c[4]);
g24:
tcore[4] = c[0] / (5.28 * adu + c[1] * 1.0 / 40.0) + tsk;
//C TRANSPARENCE
tbody = 0.1 * tsk + 0.9 * tcore[j];
swm = 304.94 * (tbody - 36.6) * adu / 3600000.0;
vpts = 6.11 * 10.0 * *(7.45 * tsk / (235.0 + tsk));
if (tbody <= 36.6)
{
swm = 0.0;
}
swf = 0.7 * swm;
if (sex == 1)
{
sw = swm;
}
if (sex == 2)
{
sw = swf;
}
eswphy = -sw * evap;
he = 0.633 * hc / (p * cair);
fec = 1.0 / (1.0 + 0.92 * hc * rcl);
eswpot = he * (vpa - vpts) * adu * evap * fec;
wetsk = eswphy / eswpot;
if (wetsk > 1.0)
{
wetsk = 1.0;
}
eswdif = eswphy - eswpot;
if (eswdif <= 0.0)
{
esw = eswpot;
}
if (eswdif > 0.0)
{
esw = eswphy;
}
if (esw > 0.0)
{
esw = 0.0;
}
//C DIFFERENCE
rdsk = 0.79 * 10.0 * *7.0;
rdcl = 0.0;
ed = evap / (rdsk + rdcl) * adu * (1 - wetsk) * (vpa - vpts);
//C VB
vb1 = 34.0 - tsk;
vb2 = tcore[j] - 36.6;
if (vb2 < 0.0)
{
vb2 = 0.0;
}
if (vb1 < 0.0)
{
vb1 = 0.0;
}
vb = (6.3 + 75.0 * (vb2)) / (1.0 + 0.5 * vb1);
//C BALANCE
enbal = h + ed + ere + esw + csum + rsum + food;
//C COVER
if (count1 == 0)
{
xx = 1.0;
}
if (count1 == 1)
{
xx = 0.1;
}
if (count1 == 2)
{
xx = 0.01;
}
if (count1 == 3)
{
xx = 0.001;
}
if (enbal > 0.0)
{
tcl = tcl + xx;
}
if (enbal < 0.0)
{
tcl = tcl - xx;
}
if ((enbal <= 0.0)) && ((enbal2 > 0.0))
{
goto g30;
}
if ((enbal >= 0.0)) && ((enbal2 < 0.0))
{
goto g30;
}
enbal2 = enbal;
count3 = count3 + 1;
//C
if (count3 > 200)
{
goto g30;
}
goto g20;
g30:
if ((count1 == 0.0))
{
|| ((count1 == 1.0)) || ((count1 == 2.0))
{
count1 = count1 + 1.0;
enbal2 = 0.0;
goto g20;
}
}
//C
if (count1 == 3.0)
{
//C
if ((j == 2))
{
|| ((j == 5)) goto g40;
}
if ((j == 6))
{
|| ((j == 1)) goto g50;
}
if (j == 3)
{
goto g60;
}
if (j == 7)
{
goto g70;
}
if (j == 4)
{
goto g80;
}
}
g40:
if (c[8] < 0.0)
{
goto g90;
}
if ((tcore[j] >= 36.6)) && ((tsk <= 34.050))
{
goto g80;
}
goto g90;
g50:
if (c[10] < 0.0)
{
goto g90;
}
if ((tcore[j] >= 36.6)) && ((tsk > 33.850))
{
goto g80;
}
goto g90;
g60:
if ((tcore[j] < 36.6)) && ((tsk <= 34.000))
{
goto g80;
}
goto g90;
g70:
if ((tcore[j] < 36.6)) && ((tsk > 34.000))
{
goto g80;
}
goto g90;
g80:
if ((j != 4)) && ((vb >= 91.0))
{
goto g90;
}
if ((j.eq.4)) && ((vb < 89.0))
{
goto g90;
}
if (vb > 90.0)
{
vb = 90.0;
}
//C LOSSES
ws = sw * 3600.0 * 1000.0;
if (ws > 2000.0)
{
ws = 2000.0;
}
wd = ed / evap * 3600.0 * (-1000.0);
wr = erel / evap * 3600.0 * (-1000.0);
wsum = ws + wr + wd;
goto g100;
}
g100:
return;
}
}