さあ！諸君！捕鯨問題だ！

IWCサイトの捕獲枠計算プログラムは、163頁が壊れてるから貼っとくわ。
http://www.iwcoffice.org/_documents/sci_com/workshops/RIWC-44-pp153-167-AnnexI.pdf
==================
REP. INT. WHAL, COMMN 44, 1991 　　 163


PTOT = 0.
N = 0
DO 50 I =0, NR -1

C Set R from the productivity parameter, Y [Documentation Eqn 5]
R = 1.4184*Y(I)

C Step size for K
DK = PKSTEP
D = 1.
RK = RKHI

C Use function STKSIM to set up the Nth population trajectory
C i. e. set up the POP array
30 IF (RK .LE. RKLO .OR. STKSIM(RK, R,POP,CATCH) .LE. 0.) GOTO 40
IF (N .GE. MAXSIM) STOP 'ERROR: TOO MANY SIMULATIONS'

C How much depletion covered? [Documentation Eqns 6 and 12]
DD = D - POP(IYEAR)/RK
D = POP (IYEAR) /RK
P = 0.0

IF (DD. GT. 0.) THEN

C Compute the internal catch limit corresponding to D and Y(I)
QRES(N) = CONTRL(D, Y(I),PLEVEL,PSLOPE) *POP(IYEAR)

Calculate deviance [Documentation Eqn 9]
CALL DEVIAN (SS0,SS1,SS2,SS3, SIGHT, FMATRX,ISYR,IZYR,ZMULT,
+ POP, G)

C Scale likelihood and integrate over values for the bias parameter
C [Documentation Eqns 9, 10 & 11]
DO 35 J = 0, NB - 1
P = P + EXP(-SF* (SS0 + RLGB(J)*(SS1 + RLGB(J)*SS2)
+ 　 + SS3*B(J) ) )
35 CONTINUE

C Calculate the weight for this point (& total area under likelihood)
C [Documentation Eqn 11]
PRES (N) = P*DD
PTOT = PTOT + PRES (N)

C Update counter
N = N+1

C Find the next K [Documentation Eqn 13]
DK = DK*PDSTEP/DD
IF (DK .GT. PKSTEP) DK = PKSTEP

ELSE
C If DD=0 change the step size onl-y [Documentation Eqn 13]
　 DK = PKSTEP
ENDIF

C Set the new value of K [Documentation Eqn 14]
RK = RK/(1. + DK)

GOTO 30
40 CONTINUE
50 CONTINUE

== =163頁おわり ========== =