PROGRAM MANAGE C ***** GLOBAL INPUT PARAMETERS COMMON /MANPAR/ PPROB,PYMAX,PNYSTP,PKSTEP,PDSTEP,PNBSTP,PBMIN, 1 PBMAX,PSCALE,PHASET,PHASEP,PCYCLE,PLEVEL,PSLOPE,PYMIN REAL PPROB, PYMAX, PNYSTP, PKSTEP, PDSTEP, PNBSTP, PBMIN, PBMAX, 1 PSCALE, PHASET, PHASEP, PCYCLE, PLEVEL, PSLOPE,PYMIN 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 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 READ (IN, '(T30,I4)') IYRCL WRITE (IOUT, '(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' ELSE WRITE (IOUT, '(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:' TOTALC = 0 DO 20 I = 0, MAXYR READ (IN, FORMT) IY, C IF (IY .LT. 0) GO TO 25 WRITE (IOUT, 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:' 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) 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) WRITE (IOUT, '(/A)') 'Zero abundance estimates:' 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) 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, '()') 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)) DO 100 IY = IYEAR, IYEAR + 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) 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,PYMIN REAL PPROB, PYMAX, PNYSTP, PKSTEP, PDSTEP, PNBSTP, PBMIN, PBMAX, 1 PSCALE, PHASET, PHASEP, PCYCLE, PLEVEL, PSLOPE,PYMIN 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,PYMIN,PNYSTP,PKSTEP, 1 PDSTEP, 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,PYMIN REAL PPROB, YMAX,PNYSTP,PKSTEP,PDSTEP,PBMIN,PBMAX,PNBSTP,PSACLE, 1 PHASET, PHASEP, PCYCLE, PLEVEL, PSLOPE,PYMIN 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,PYMIN REAL PPROB,PYMAX,PNYSTP,PKSTEP,PDSTEP,PNBSTP,PBMIN, PBMAX, 1 PSCALE,PHASET,PHASEP,PCYCLE,PLEVEL,PSLOPE,PYMIN 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 - PYMIN) / NR DO 20 I = 0, NR - 1 Y(I) = (I + 0.5) * YINC + PYMIN 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