PROGRAM MANAGE C ***** GLOBAL INPUT PARAMETERS COMMON /MANPAR/ PPROB,PYMAX,PNYSTP,PKSTEP,PDSTEP,PNBSTP,PBMIN, 1 PBMAX,PSCALE,PHASET,PHASEP,PCYCLE,PLEVEL,PSLOPE REAL PPROB, PYMAX, PNYSTP, PKSTEP, PDSTEP, PNBSTP, PBMIN, PBMAX, 1 PSCALE, PHASET, PHASEP, PCYCLE, PLEVEL, PSLOPE COMMON /MANDAT/ ISTART, IYEAR, NS, NZ, RKLO, RKHI INTEGER ISTART, IYEAR, NS, NZ REAL RKLO, RKHI PARAMETER (MAXYR = 200, MAXEST = 100) PARAMETER (MSIZE = (MAXEST * (MAXEST + 1) / 2)) REAL CATCH(0:MAXYR), SIGHT(0:MAXEST), FMATRX(0:MSIZE), 1 ZMULT(0:MAXEST), POP(0:MAXYR), G(0:MAXYR), RAWCL, CL INTEGER ISYR(0:MAXEST), IZYR(0:MAXEST), POUT C ***** LOCAL VARIABLES CHARACTER STOCK*30, FORMT*50 INTEGER IY, INITYR, I, IS, IYRCL, ILAST, N, IN, IOUT C ***** Read in the data required for application of the CLA DATA IN/2/, IOUT/3/ OPEN (IN, FILE = 'CLC.DAT') OPEN (IOUT, FILE = 'CLC.OUT') CALL RDPARS READ (IN, '(A30/)') STOCK WRITE (IOUT, '(A30/)') STOCK WRITE (*, '(A30/)') STOCK C Read the year of the first catch, the first year for which catch C limits are to be set & the phaseout option READ (IN, '(T30,I4)') INITYR WRITE (IOUT, '(A,T30,I4)') 'Year of first input data', INITYR WRITE (*, '(A,T30,I4)') 'Year of first input data', INITYR READ (IN, '(T30,I4)') IYRCL WRITE (IOUT, '(A, T30, I4)') 'Year of first catch limit', IYRCL WRITE (*, '(A, T30, I4)') 'Year of first catch limit', IYRCL READ (IN, '(T30, I4)') POUT IF (POUT .EQ. 1) THEN WRITE (IOUT, '(A, T30,A4)') 'Apply phaseout if necessary', 'YES' WRITE (*, '(A, T30,A4)') 'Apply phaseout if necessary', 'YES' ELSE WRITE (IOUT, '(A, T30,A4)') 'Apply phaseout', 'No' WRITE (*, '(A, T30,A4)') 'Apply phaseout', 'No' ENDIF C Re-scale IYRCL such that 0 is the year prior to the first input data ISTART = 0 IYEAR = IYRCL - INITYR IF (IYEAR .LE. 0 .OR. IYEAR .GT. MAXYR) STOP 'INVALID YEAR' C Initialize the catch array DO 10 I = 0, MAXYR CATCH(I) = 0 10 CONTINUE C Read in the catch data, scaling each year to the initial year READ (IN, '(// A)') FORMT WRITE (IOUT, '(/A)') 'Historic catches:' WRITE (*, '(/A)') 'Historic catches:' TOTALC = 0 DO 20 I = 0, MAXYR READ (IN, FORMT) IY, C IF (IY .LT. 0) GO TO 25 WRITE (IOUT, FORMT) IY, C WRITE (*, FORMT) IY, C IY = IY - INITYR IF (IY .LT. 0 .OR. IY .GE. IYEAR) STOP CATCH(IY) = C TOTALC = TOTALC + C 20 CONTINUE 25 IF (TOTALC .LT. 0) STOP ' ERROR: No historic catch input' C Read in the non-zero sightings estimates and information matrix READ (IN, '(//T30, I4 / A)') NS, FORMT WRITE (IOUT, '(/A)') 'Abundance estimates:' WRITE (*, '(/A)') 'Abundance estimates:' IF (NS .GT. MAXEST) STOP 'ERROR: Abundance year out of range' DO 30 N=0,NS-1 N1 = (N * (N + 1)) / 2 READ (IN, FORMT) ISYR(N), SIGHT(N), (FMATRX(N1 + J), J = 0, N) WRITE (IOUT, FORMT) ISYR(N), SIGHT(N), 1 (FMATRX (N1 + J), j = 0, N) WRITE (*, FORMT) ISYR(N), SIGHT(N), 1 (FMATRX (N1 + J), j = 0, N) ISYR(N) = ISYR(N) - INITYR IF (ISYR(N) .LT. 0 .OR. ISYR(N) .GE. IYEAR) STOP 1 'ERROR: Sight year out of range' IF (SIGHT(N) .LE. 0.0) STOP ' ERROR: Estimate not positive' 30 CONTINUE C Read in the Poisson multipliers for any zero sightings estimates READ (IN, '(/T30, I3, /A)') NZ, FORMAT IF (NZ .GT. 0) THEN WRITE (IOUT, '(/A)') 'Zero abundance estimates:' WRITE (*, '(/A)') 'Zero abundance estimates:' ENDIF IF (NZ .GT. MAXEST) STOP ' ERROR: Zero estimate array too small' DO 40 N = 0, NZ - 1 N1 = (N * (N + 1)) / 2 READ (IN, FORMT) IZYR(N), ZMULT(N) WRITE (IOUT, FORMT) IZYR(N), ZMULT(N) WRITE (*, FORMT) IZYR(N), ZMULT(N) IZYR(N) = IZYR(N) - INITYR IF (IZYR(N) .LT. 0 .OR. IZYR(N) .GE. IYEAR) STOP 1 ' Sight year out of range' IF (ZMULT(N) .LE. 0.0) STOP ' ERROR: Multiplier not positive' 40 CONTINUE WRITE (IOUT, '()') WRITE (*, '()') C Bound the range for the integration over K RKHI = 1.E7 RKLO = 0.0 C ****** C ****** Run the CLA to obtain the nominal catch limit C ****** RAWCL = CLIMIT (CATCH, SIGHT, FMATRX, ISYR, IZYR, ZMULT, POP, G) C Set the catch limits for PCYCLE years. If the catch limit may be C subject to phaseout, call POUT to apply the phaseout rule. C First set ILAST = year of the most recent abundance estimate IF (NS .GT. 0) ILAST = ISYR(NS - 1) IF (NZ .GT. 0) ILAST = MAX(IS, IZYR(NZ - 1)) C DO 100 IY = IYEAR, IYEAR + PCYCLE - 1 DO 100 IY = IYEAR, IYEAR + int(PCYCLE) - 1 IF (POUT .EQ. 1) THEN CALL PHOUT (CL, RAWCL, ILAST, IY) ELSE CL = RAWCL ENDIF WRITE (IOUT, '(A6, I5, A17, I6)') 'Year:', IY + INITYR, 1 'Catch limit:', NINT(CL) WRITE (* , '(A6, I5, A17, I6)') 'Year:', IY + INITYR, 1 'Catch limit:', NINT(CL) 100 CONTINUE STOP END C ************************************************************************************** C ************************************************************************************** C C CLC version 6 C C ************************************************************************************** C C 31 January 1994 C C ************************************************************************************** SUBROUTINE RDPARS C Read the file of input parameters COMMON /MANPAR/ PPROB,PYMAX,PNYSTP,PKSTEP,PDSTEP,PNBSTP,PBMIN, 1 PBMAX,PSCALE,PHASET,PHASEP,PCYCLE,PLEVEL,PSLOPE REAL PPROB, PYMAX, PNYSTP, PKSTEP, PDSTEP, PNBSTP, PBMIN, PBMAX, 1 PSCALE, PHASET, PHASEP, PCYCLE, PLEVEL, PSLOPE DATA IFILE /21/ C C Open the input file (only read once) C OPEN (UNIT = IFILE, FILE = 'CLC.PAR') C READ (IFILE, '(T30, F10.0)') PPROB,PYMAX, PNYSTP, PKSTEP, PDSTEP, 1 PBMIN, PBMAX, PNBSTP, PSCALE, PHASET, 1 PHASEP, PCYCLE, PLEVEL, PSLOPE CLOSE (IFILE) RETURN END C C C ************************************************************************************** C ************************************************************************************** C SUBROUTINE PHOUT (CL, RAWCL, ILAST, IY) C COMMON /MANPAR/ PPROB,PYMAX,PNYSTP,PKSTEP,PDSTEP,PBMIN,PBMAX, 1 PNBSTP,PSACLE,PHASET,PHASEP,PCYCLE,PLEVEL,PSLOPE REAL PPROB, YMAX,PNYSTP,PKSTEP,PDSTEP,PBMIN,PBMAX,PNBSTP,PSACLE, 1 PHASET, PHASEP, PCYCLE, PLEVEL, PSLOPE C PHASET Number of years without surveys before phaseout invoked C PHASEP Phaseout annual reduction proportion C REAL CL, RAWCL INTEGER ILAST, IY C C Phaseout: Reduce catch limit if there is no survey data in the C last PHASET years IF (IY .GE. ILAST + PHASET) THEN CL = RAWCL * (1.0 - PHASEP * (IY - ILAST - PHASET)) IF (CL .LT. 0.0) CL = 0.0 ELSE CL = RAWCL ENDIF C RETURN END C *************************************************************************************** REAL FUNCTION CLIMIT (CATCH,SIGHT,FMATRX,ISYR,IZYR,ZMULT,POP,G) C Run the CLA to obtain to nominal catch limit C PARAMETER (MAXSIM = 600000, MAXSTP = 500) REAL PRES(0:MAXSIM), QRES(0:MAXSIM), SS0, SS1, SS2, SS3 C COMMON /MANPAR/ PPROB,PYMAX,PNYSTP,PKSTEP,PDSTEP,PNBSTP,PBMIN, 1 PBMAX,PSCALE,PHASET,PHASEP,PCYCLE,PLEVEL,PSLOPE REAL PPROB,PYMAX,PNYSTP,PKSTEP,PDSTEP,PNBSTP,PBMIN, PBMAX, 1 PSCALE,PHASET,PHASEP,PCYCLE,PLEVEL,PSLOPE C COMMON /MANDAT/ ISTART, IYEAR, NS, NZ, RKLO, RKHI INTEGER ISTART, IYEAR, NS, NZ REAL RKLO, RKHI C REAL CATCH(0:*),SIGHT(0:*),FMATRX(0:*),ZMULT(0:*),POP(0:*),G(0:*) INTEGER ISYR(0:*), IZYR(0:*) C C Local variables: REAL SF, Y(0:MAXSTP), B(0:MAXSTP), RLGB(0:MAXSTP) INTEGER NB, NR C NB Number of bias steps C NR Number of steps for Y (the productivity parameter) C SF Deviance scale factor (SF = .5 / PSCALE**2) C CLIMIT = 0. IF (NS .LE. 0) RETURN C C Set deviance scale factor S = 1/PSCALE**2 SF = 0.5 / (PSCALE * PSCALE) C C Check the sizes of the Y and B arrays are large enough IF (PNBSTP .GT. MAXSTP .OR. PNYSTP .GT. MAXSTP) STOP 1 'Y &/or b array sizes too small' C Set sightings bias step sizes & their log values. BINC = Bias increment NB = PNBSTP BINC = (PBMAX - PBMIN) / NB DO 10 I = 0, NB - 1 B(I) = PBMIN + (I + 0.5) * BINC RLGB(I) = -ALOG(B(I)) 10 CONTINUE C C Set productivity parameter step sizes (midpoints) NR = PNYSTP YINC = PYMAX / NR DO 20 I = 0, NR - 1 Y(I) = (I + 0.5) * YINC 20 CONTINUE C PTOT = 0 N = 0 DO 50 I = 0, NR - 1 C Set R from the productivity parameter, Y 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.) 1 GOTO 40 IF (N .GE. MAXSIM) STOP 'ERROR: TOO MANY SIMULATIONS' C How much depletion covered? 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) C Calculate deviance CALL DEVIAN (SS0, SS1, SS2, SS3, SIGHT, FMATRX, ISYR, 1 IZYR, ZMULT, POP, G) C Scale likelihood and integrate over values for the bias parameter DO 35 J = 0, NB - 1 P = P + EXP(-SF * (SS0 + RLGB(J) * (SS1 + RLGB(J) * SS2) 1 + SS3 * B(J))) 35 CONTINUE C Calculate the weight for this point (& total area under likelihood) PRES(N) = P * DD PTOT = PTOT + PRES(N) C Update counter N = N + 1 C Find the next K DK = DK * PDSTEP / DD IF (DK .GT. PKSTEP) DK = PKSTEP ELSE C IF DD = 0 change the step size only DK = PKSTEP ENDIF C Set the new value of K RK = RK / (1. + DK) GOTO 30 40 CONTINUE 50 CONTINUE IF (PTOT .LE. 0.) STOP 'ERROR: PROB INTEGRATES TO ZERO' C Sort the QRES and PRES arrays in ascending order of QRES N2 = N CALL SORT (QRES, PRES, N2) C Normalize the relative likelihoods DO 60 I = 0, N - 1 PRES(I) = PRES(I) / PTOT 60 CONTINUE C Extract the desired probability level: the nominal catch limit (NCL) C is the lower PROB% of the distribution. C First calculate PRES(I), the probability that the NCL is between C QRES(I) & QRES(I + 1). P = 0 DO 70 I = 0, N - 1 P = P + PRES (I) IF (P .GT. PPROB) GOTO 80 70 CONTINUE C Interpolate to set the nominal catch limit 80 IF (I .GE. N - 1) THEN Q = QRES (N - 1) ELSE Q = (QRES (I + 1) * (PPROB - P + PRES(I)) + QRES(I) * 1 (P - PPROB)) / PRES(I) ENDIF CLIMIT = Q RETURN END C *********************************************************************************** FUNCTION CONTRL (D, Y, PLEVEL, PSLOPE) C Catch control law IF (D .LT. PLEVEL) THEN CONTRL = 0. ELSEIF (D .LT. 1.) THEN CONTRL = PSLOPE * Y * (D - PLEVEL) ELSE CONTRL = PSLOPE * Y * (1.0 - PLEVEL) ENDIF END C ********************************************************************************** SUBROUTINE DEVIAN (SS0, SS1, SS2, SS3, SIGHT, FMATRX, ISYR, 1 IZYR, ZMULT, POP, G) C Calculate deviance (-2 log likelihood) in terms of coefficients for C the bias and log bias parameters COMMON /MANDAT/ ISTART, IYEAR, NS, NZ, RKLO, RKHI INTEGER ISTART, IYAER, NS, NZ REAL RKLO, RKHI REAL SS0,SS1,SS2,SIGHT(0:*),FMATRX(0:*),ZMULT(0:*),POP(0:*),G(0:*) INTEGER ISYR(0:*), IZYR(0:*) SS0 = 0 SS1 = 0 SS2 = 0 SS3 = 0 DO 100 N = 0, NS - 1 G(N) = ALOG (SIGHT(N) / POP(ISYR(N))) K = N * (N + 1) / 2 DO 10 J = 0, N - 1 C 1st add non diagonal contributions (which are doubled up) SS0 = SS0 + 2. * G(J) * G(N) * FMATRX(K) SS1 = SS1 + 2. * (G(J) + G(N)) * FMATRX(K) SS2 = SS2 + FMATRX(K) + FMATRX(K) K = K + 1 10 CONTINUE C Now add diagnoal contribution SS0 = SS0 + G(N) * G(N) * FMATRX(K) SS1 = SS1 + 2. * G(N) * FMATRX(K) SS2 = SS2 + FMATRX(K) 100 CONTINUE C Now do the zero estimates DO 200 N = 0, NZ - 1 SS3 = SS3 + 2.0 * POP(IZYR(N)) / ZMULT(N) 200 CONTINUE RETURN END C ********************************************************************************* FUNCTION STKSIM (RK, R, POP, CATCH) C Calculate the stock trajectory with parameter RK and R: return C the current stock size C RK = notional carrying capacity; R = productivity parameter * 1.4184 COMMON /MANDAT/ ISTART, IYEAR, NS, NZ, RKLO, RKHI INTEGER ISTART, IYEAR, NS, NZ REAL RKLO, RKHI REAL CATCH(0:*), POP(0:*) POP(ISTART) = RK DO 10 I = ISTART + 1, IYEAR D = POP(I - 1) / RK POP(I) = POP (I - 1) * (1. + R * (1. - D * D)) - CATCH(I-1) IF (POP(I) .LE. 0) GOTO 20 10 CONTINUE STKSIM = POP(IYEAR) RETURN 20 STKSIM = 0. END C ********************************************************************************* SUBROUTINE SORT (ARRAY, ARRAY2, N) C SORT sorts a pair of arrays in ascending order of the first array C (using the heapsort algorithm) REAL ARRAY(0:*), ARRAY2(0:*), TEMP, TEMP2 INTEGER N, K, IR, I, J IF (N .LT. 2) RETURN K = N / 2 IR = N - 1 10 IF (K .NE. 0) THEN K = K - 1 TEMP = ARRAY(K) TEMP2 = ARRAY2(K) ELSE TEMP = ARRAY(IR) TEMP2 = ARRAY2(IR) ARRAY(IR) = ARRAY(0) ARRAY2(IR) = ARRAY2(0) IR = IR - 1 IF (IR .EQ. 0) THEN ARRAY(0) = TEMP ARRAY2(0) = TEMP2 RETURN ENDIF ENDIF I = K J = K + K + 1 20 IF (J .LE. IR) THEN IF (J .LT. IR .AND. ARRAY(J) .LT. ARRAY(J + 1)) J = J + 1 IF (TEMP .LT. ARRAY(J)) THEN ARRAY(I) = ARRAY(J) ARRAY2(I) = ARRAY2(J) I = J J = J + I + 1 ELSE J = IR + 1 ENDIF GOTO 20 ENDIF ARRAY(I) = TEMP ARRAY2(I) = TEMP2 GOTO 10 END