```c   This program determines appropriate kiln
c   schedules given the specific gravity of the species.
c
c
c   Steve Verrill, 6/10/96
c
c
c

implicit double precision (a-h,o-z)

integer temp(40),depr(40),mc(40)
integer number(40)

integer harsh(40)

integer iorigord(40)

integer thick

real*8 gb(40),sdgb(40)

real*8 logsdgb(40),logprob(40),sdgb22(40)

real*8 loglik(40),liksort(40)

real temp0,depr0,mc0
integer temp1,depr1,mc1,harsh1
integer temp2,depr2

c      character*20 results

character*80 argv,intwrite
character*1 t,mcsched,alpha(6)

data alpha/"A","B","C","D","E","F"/

t = "T"

temp( 1) =      1
temp( 2) =      2
temp( 3) =      2
temp( 4) =      2
temp( 5) =      2
temp( 6) =      3
temp( 7) =      3
temp( 8) =      3
temp( 9) =      3
temp(10) =      4
temp(11) =      5
temp(12) =      5
temp(13) =      5
temp(14) =      6
temp(15) =      6
temp(16) =      6
temp(17) =      6
temp(18) =      6
temp(19) =      6
temp(20) =      6
temp(21) =      7
temp(22) =      7
temp(23) =      8
temp(24) =      8
temp(25) =      8
temp(26) =      8
temp(27) =      8
temp(28) =      8
temp(29) =     10
temp(30) =     10
temp(31) =     10
temp(32) =     10
temp(33) =     10
temp(34) =     10
temp(35) =     11
temp(36) =     12
temp(37) =     12
temp(38) =     12
temp(39) =     13
temp(40) =     14

depr( 1) =      1
depr( 2) =      2
depr( 3) =      2
depr( 4) =      2
depr( 5) =      4
depr( 6) =      1
depr( 7) =      1
depr( 8) =      2
depr( 9) =      2
depr(10) =      2
depr(11) =      2
depr(12) =      2
depr(13) =      3
depr(14) =      2
depr(15) =      3
depr(16) =      3
depr(17) =      3
depr(18) =      3
depr(19) =      4
depr(20) =      4
depr(21) =      3
depr(22) =      6
depr(23) =      2
depr(24) =      3
depr(25) =      3
depr(26) =      4
depr(27) =      4
depr(28) =      4
depr(29) =      4
depr(30) =      4
depr(31) =      4
depr(32) =      4
depr(33) =      5
depr(34) =      5
depr(35) =      4
depr(36) =      5
depr(37) =      5
depr(38) =      7
depr(39) =      4
depr(40) =      6

mc( 1) =      2
mc( 2) =      2
mc( 3) =      3
mc( 4) =      4
mc( 5) =      4
mc( 6) =      3
mc( 7) =      4
mc( 8) =      3
mc( 9) =      4
mc(10) =      2
mc(11) =      2
mc(12) =      4
mc(13) =      3
mc(14) =      4
mc(15) =      1
mc(16) =      2
mc(17) =      3
mc(18) =      6
mc(19) =      1
mc(20) =      4
mc(21) =      2
mc(22) =      2
mc(23) =      3
mc(24) =      2
mc(25) =      4
mc(26) =      2
mc(27) =      3
mc(28) =      4
mc(29) =      3
mc(30) =      4
mc(31) =      5
mc(32) =      6
mc(33) =      4
mc(34) =      6
mc(35) =      4
mc(36) =      4
mc(37) =      5
mc(38) =      5
mc(39) =      3
mc(40) =      3

number( 1) =      2
number( 2) =      1
number( 3) =     25
number( 4) =      1
number( 5) =     19
number( 6) =      2
number( 7) =      2
number( 8) =     40
number( 9) =      9
number(10) =      1
number(11) =      1
number(12) =      1
number(13) =      1
number(14) =     63
number(15) =      1
number(16) =      1
number(17) =      2
number(18) =      1
number(19) =      1
number(20) =     18
number(21) =      2
number(22) =      1
number(23) =      1
number(24) =      3
number(25) =      1
number(26) =      2
number(27) =      3
number(28) =      4
number(29) =      1
number(30) =     20
number(31) =      3
number(32) =      2
number(33) =     16
number(34) =      2
number(35) =      1
number(36) =      1
number(37) =      1
number(38) =      2
number(39) =      9
number(40) =      1

gb( 1) =    0.91450
gb( 2) =    0.65900
gb( 3) =    0.75412
gb( 4) =    0.53200
gb( 5) =    0.67453
gb( 6) =    0.86250
gb( 7) =    0.61950
gb( 8) =    0.64773
gb( 9) =    0.68100
gb(10) =    0.58000
gb(11) =    0.74000
gb(12) =    0.53000
gb(13) =    0.41000
gb(14) =    0.57454
gb(15) =    0.66000
gb(16) =    0.40000
gb(17) =    0.62500
gb(18) =    0.40000
gb(19) =    0.51000
gb(20) =    0.50856
gb(21) =    0.75000
gb(22) =    0.36000
gb(23) =    0.56000
gb(24) =    0.61000
gb(25) =    0.60000
gb(26) =    0.51000
gb(27) =    0.49333
gb(28) =    0.47000
gb(29) =    0.48000
gb(30) =    0.44905
gb(31) =    0.39800
gb(32) =    0.34500
gb(33) =    0.44806
gb(34) =    0.34000
gb(35) =    0.40000
gb(36) =    0.22600
gb(37) =    0.46000
gb(38) =    0.33500
gb(39) =    0.35011
gb(40) =    0.33400

sdgb( 1) =    0.08273
sdgb( 2) =    0.06921
sdgb( 3) =    0.10557
sdgb( 4) =    0.06921
sdgb( 5) =    0.17456
sdgb( 6) =    0.06718
sdgb( 7) =    0.07000
sdgb( 8) =    0.11637
sdgb( 9) =    0.07114
sdgb(10) =    0.06921
sdgb(11) =    0.06921
sdgb(12) =    0.06921
sdgb(13) =    0.06921
sdgb(14) =    0.10735
sdgb(15) =    0.06921
sdgb(16) =    0.06921
sdgb(17) =    0.02121
sdgb(18) =    0.06921
sdgb(19) =    0.06921
sdgb(20) =    0.10029
sdgb(21) =    0.04243
sdgb(22) =    0.06921
sdgb(23) =    0.06921
sdgb(24) =    0.14000
sdgb(25) =    0.06921
sdgb(26) =    0.05657
sdgb(27) =    0.05508
sdgb(28) =    0.06782
sdgb(29) =    0.06921
sdgb(30) =    0.11489
sdgb(31) =    0.03704
sdgb(32) =    0.02121
sdgb(33) =    0.08762
sdgb(34) =    0.04243
sdgb(35) =    0.06921
sdgb(36) =    0.06921
sdgb(37) =    0.06921
sdgb(38) =    0.02121
sdgb(39) =    0.13640
sdgb(40) =    0.06921

harsh( 1) =     1
harsh( 2) =     1
harsh( 3) =     1
harsh( 4) =     1
harsh( 5) =     1
harsh( 6) =     2
harsh( 7) =     2
harsh( 8) =     2
harsh( 9) =     2
harsh(10) =     2
harsh(11) =     3
harsh(12) =     3
harsh(13) =     3
harsh(14) =     3
harsh(15) =     3
harsh(16) =     3
harsh(17) =     3
harsh(18) =     3
harsh(19) =     3
harsh(20) =     3
harsh(21) =     4
harsh(22) =     4
harsh(23) =     4
harsh(24) =     4
harsh(25) =     4
harsh(26) =     4
harsh(27) =     4
harsh(28) =     4
harsh(29) =     5
harsh(30) =     5
harsh(31) =     5
harsh(32) =     5
harsh(33) =     5
harsh(34) =     5
harsh(35) =     6
harsh(36) =     7
harsh(37) =     7
harsh(38) =     8
harsh(39) =     9
harsh(40) =     10

c 4    write(6,5)
c 5    format(/,1x,'What name do you want for your results file?',/)
c 7    format(a80)

c 9    open(unit=7,file=results,status="new",err=1990)
c      go to 2000

c 1990 write(6,1992) results
c 1992 format(/,1x,'It was not possible to open file',1x,a20,/,1x,
c     x'for writing.  Please use a different file.',/)
c      go to 4

2000 do 10 i = 1,40

logprob(i) = log(number(i)/268.0)
logsdgb(i) = log(sdgb(i))
sdgb22(i) = 2.0*sdgb(i)*sdgb(i)

10   continue

c      write(6,20)
c      write(7,20)
c 20   format(/,1x,'What is specific gravity of the species in',
c     x' which you are interested?',/)
c      write(7,*) sg0

call getarg(1,argv)
write(intwrite,7) argv
7    format(a80)

call getarg(2,argv)
write(intwrite,7) argv

do 40 i = 1,40

loglik(i) = -logsdgb(i) - (sg0 - gb(i))*
x          (sg0 - gb(i))/sdgb22(i)
x          + logprob(i)

40   continue

c
c   sort the loglik's
c

call rsort(loglik,liksort,iorigord,40)

c
c   obtain the estimated schedules from the regression approach
c
temp0 = 13.65 - 13.62*sg0
depr0 =  5.20 -  3.95*sg0
mc0   =  4.51 -  1.56*sg0

temp1 = nint(temp0)
depr1 = nint(depr0)
mc1 = nint(mc0)

write(6,44) temp1,mc1,depr1
c      write(7,44) temp1,depr1,mc1

44   format('We give results for both the regression',
x' and classification approaches below.  If the two',
x' approaches yield schedules that differ significantly,',
x' it might be wise to choose the mildest of the',
x' schedules.  In this case we would also appreciate',
x' hearing from you (there might be a bug in the program)',
x'',
x'.',
x'Regression Approach',
x'Rounded to the nearest integer, the temp,'
x' mc, and deprpredicted',
x' by the regression equations are:',
x1x,3(i3,2x))

if (thick .eq. 50) then

if (temp1 .lt. 5) then

depr1 = depr1 - 1

else if (temp1 .eq. 5) then

temp1 = temp1 - 2
depr1 = depr1 - 1

else if (temp1 .lt. 9) then

temp1 = temp1 - 3
depr1 = depr1 - 1

else

temp1 = temp1 - 2
depr1 = depr1 - 1

endif

write(6,41) temp1,mc1,depr1
41      format('Because the thickness is 50 mm,',
x   ' the temp, mc, and depr values are changed',
x   ' to:',
x   1x,3(i3,2x))

endif

mcsched = alpha(mc1)

write(6,42)
42   format('')

if (temp1 .lt. 10) then

if (temp1 .eq. 0) temp1 = 1
if (depr1 .eq. 0) depr1 = 1

write(6,1045) t,temp1,mcsched,depr1
1045    format('')
write(6,45) t,temp1,mcsched,depr1
45      format('This corresponds to the schedule',
x   2x,'.(Click on the box',
x   ' to obtain the full schedule.)')

else

write(6,1046) t,temp1,mcsched,depr1
1046    format('')
write(6,46) t,temp1,mcsched,depr1
46      format('This corresponds to the schedule',
x   2x,'.(Click on the box',
x   ' to obtain the full schedule.)')

endif

write(6,48)
48   format('')

c 44   format(1x,'Rounded to the nearest integer, the temp,'
c     x' depr, and mc',/,1x,'predicted',
c     x' by the regression equations are',//,
c     x1x,3(i3,2x),//)

c23456789012345678901234567890123456789012345678901234567890123456789012

if (thick .eq. 2538) then

write(6,49)
49      format('Classification Approach',
x   'In the table below we list the schedules with',
x   ' the 5 largest relative likelihoods.  These are the',
x   ' schedules',
x   ' (among the 40 currently in use) that are most likely to be',
x   ' appropriate for the specific gravity under consideration.',
x   ' ')

else

write(6,1049)
1049    format('Classification Approach',
x   'In the table below we list the schedules with',
x   ' the 5 largest relative likelihoods.  These are',
x   ' the schedules',
x   ' (among the 40 currently in use) that are most likely to be',
x   ' appropriate for the specific gravity under consideration.',
x   '  They are then modified',
x   ' to reflect the fact that the thickness',
x   ' is 50 mm.')

endif

if (thick .eq. 2538) then

write(6,50)
50      format('The temperature, moisture content,',
x   ' wet bulb depression, harshness, ',
x   'relative likelihood, and schedule values',
x   ' are (click on a box to obtain a detailed schedule)',
x   ':')

write(6,8010)
8010    format("")

write(6,8012)
8012    format(" ")

write(6,8014)

8014    format(" temp  mc  depr  harsh  rel likeli",
x   "  schedule")

do 70 i = 40,36,-1

id = iorigord(i)
diff = liksort(i) - liksort(40)
if (diff .lt. -100.0d0) then
ratio = 0.0d0
else
ratio = dexp(liksort(i) - liksort(40))
endif

temp1 = temp(id)
mc1 = mc(id)
depr1 = depr(id)
harsh1 = harsh(id)
mcsched = alpha(mc1)

if (temp1 .lt. 10) then

write(6,65) temp1,mc1,depr1,harsh1,ratio,
x                     t,temp1,mcsched,depr1,
x                     t,temp1,mcsched,depr1
65            format('',i2,' ',i1,' ',i1,
x                ' ',i2,' ',f5.3,' ',
x                '',
x'',
x'')

else

write(6,66) temp1,mc1,depr1,harsh1,ratio,
x                     t,temp1,mcsched,depr1,
x                     t,temp1,mcsched,depr1
66            format('',i2,' ',i1,' ',i1,
x                ' ',i2,' ',f5.3,' ',
x                '',
x'',
x'')

endif

70      continue

write(6,8040)
8040    format("")

else

write(6,1050)
1050    format('The temperature, moisture content,',
x   ' wet bulb depression, harshness, ',
x   'relative likelihood, schedule, and modified schedule',
x   ' values',
x   ' are (click on a box to obtain a detailed modified',
x   ' schedule)',
x   ':')

write(6,8010)

write(6,8012)

write(6,9014)
9014    format(" temp  mc  depr  harsh  rel likeli",
x   "  schedule  modified schedule")

do 1070 i = 40,36,-1

id = iorigord(i)
diff = liksort(i) - liksort(40)
if (diff .lt. -100.0d0) then
ratio = 0.0d0
else
ratio = dexp(liksort(i) - liksort(40))
endif

temp1 = temp(id)
mc1 = mc(id)
depr1 = depr(id)
harsh1 = harsh(id)
mcsched = alpha(mc1)

if (temp1 .lt. 5) then

temp2 = temp1
depr2 = depr1 - 1

else if (temp1 .eq. 5) then

temp2 = temp1 - 2
depr2 = depr1 - 1

else if (temp1 .lt. 9) then

temp2 = temp1 - 3
depr2 = depr1 - 1

else

temp2 = temp1 - 2
depr2 = depr1 - 1

endif

if (temp2 .lt. 1) temp2 = 1
if (depr2 .lt. 1) depr2 = 1

if (temp1 .lt. 10) then

write(6,1065) temp1,mc1,depr1,harsh1,ratio,
x                     t,temp1,mcsched,depr1,
x                     t,temp2,mcsched,depr2,
x                     t,temp2,mcsched,depr2
1065          format('',i2,' ',i1,' ',i1,
x                ' ',i2,' ',f5.3,' ',
x                a1,i1,'-',a1,i1,' ',
x                '',
x'',
x'')

else if (temp2 .lt. 10) then

write(6,1066) temp1,mc1,depr1,harsh1,ratio,
x                     t,temp1,mcsched,depr1,
x                     t,temp2,mcsched,depr2,
x                     t,temp2,mcsched,depr2
1066          format('',i2,' ',i1,' ',i1,
x                ' ',i2,' ',f5.3,' ',
x                a1,i2,'-',a1,i1,' ',
x                '',
x'',
x'')

else

write(6,1067) temp1,mc1,depr1,harsh1,ratio,
x                     t,temp1,mcsched,depr1,
x                     t,temp2,mcsched,depr2,
x                     t,temp2,mcsched,depr2
1067          format('',i2,' ',i1,' ',i1,
x                ' ',i2,' ',f5.3,' ',
x                a1,i2,'-',a1,i1,' ',
x                '',
x'',
x'')

endif

1070    continue

write(6,8040)

endif

stop

end

c   This subroutine attempts to perform an n*log(n) sort
c   rather than an n*n sort.
c   It sorts reals.
c   It also keeps track of the original order of the x's.
c
c
c   Steve Verrill, 5/2/89
c
c
c
subroutine rsort(x,xsort,iorigord,n)

implicit double precision (a-h,o-z)

real*8 x(n),xsort(n)
integer iorigord(n)

do 10 i = 1,n

iorigord(i) = i

10    continue

num = 2

if (x(1) .lt. x(2)) then

xsort(1) = x(1)
iorigord(1) = 1
xsort(2) = x(2)
iorigord(2) = 2

else

xsort(1) = x(2)
iorigord(1) = 2
xsort(2) = x(1)
iorigord(2) = 1

endif

do 1000 j = 3,n

c
c   insert the next value
c

if (x(j) .ge. xsort(num)) go to 200
if (x(j) .le. xsort(1)) go to 160

il = 1
iu = num
145   numb = iu - il
if (numb .eq. 1) go to 170
numbd2 = numb/2
itest = il + numbd2
if (x(j) .le. xsort(itest)) go to 150
il = itest
go to 145
150   iu = itest
go to 145

160   il = 0
170   num = num + 1

do 180 k = num,il + 2,-1

km1 = k - 1
xsort(k) = xsort(km1)
iorigord(k) = iorigord(km1)

180   continue

190   ilp1 = il + 1
xsort(ilp1) = x(j)
iorigord(ilp1) = j
go to 1000

200   num = num + 1
xsort(num) = x(j)
iorigord(num) = j

1000  continue

return

end
```