fork download 
  1.       PROGRAM MANAGE
  2.  
  3. C ***** GLOBAL INPUT PARAMETERS
  4.  
  5.       COMMON /MANPAR/ PPROB,PYMAX,PNYSTP,PKSTEP,PDSTEP,PNBSTP,PBMIN,
  6.      1             PBMAX,PSCALE,PHASET,PHASEP,PCYCLE,PLEVEL,PSLOPE
  7.       REAL PPROB, PYMAX, PNYSTP, PKSTEP, PDSTEP, PNBSTP, PBMIN, PBMAX,
  8.      1	   PSCALE, PHASET, PHASEP, PCYCLE, PLEVEL, PSLOPE
  9.  
  10.       COMMON /MANDAT/ ISTART, IYEAR, NS, NZ, RKLO, RKHI
  11.       INTEGER ISTART, IYEAR, NS, NZ
  12.       REAL RKLO, RKHI
  13.  
  14.       PARAMETER (MAXYR = 200, MAXEST = 100)
  15.       PARAMETER (MSIZE = (MAXEST * (MAXEST + 1) / 2))
  16.       REAL CATCH(0:MAXYR), SIGHT(0:MAXEST), FMATRX(0:MSIZE),
  17.      1     ZMULT(0:MAXEST), POP(0:MAXYR), G(0:MAXYR), RAWCL, CL
  18.       INTEGER ISYR(0:MAXEST), IZYR(0:MAXEST), POUT
  19.  
  20. C ***** LOCAL VARIABLES
  21.  
  22.       CHARACTER STOCK*30, FORMT*50
  23.       INTEGER IY, INITYR, I, IS, IYRCL, ILAST, N, IN, IOUT
  24.  
  25. C ***** Read in the data required for application of the CLA
  26.  
  27.       DATA IN/2/, IOUT/3/
  28.       OPEN (IN, FILE = 'CLC.DAT')
  29.       OPEN (IOUT, FILE = 'CLC.OUT')
  30.  
  31.       CALL RDPARS
  32.       READ (IN, '(A30/)') STOCK
  33.       WRITE (IOUT, '(A30/)') STOCK
  34.       WRITE (*,    '(A30/)') STOCK
  35.  
  36. C 	Read the year of the first catch, the first year for which catch
  37. C	limits are to be set & the phaseout option
  38.  
  39.       READ (IN, '(T30,I4)') INITYR
  40.       WRITE (IOUT, '(A,T30,I4)') 'Year of first input data', INITYR
  41.       WRITE (*,    '(A,T30,I4)') 'Year of first input data', INITYR
  42.       READ (IN, '(T30,I4)') IYRCL
  43.       WRITE (IOUT, '(A, T30, I4)') 'Year of first catch limit', IYRCL
  44.       WRITE (*,    '(A, T30, I4)') 'Year of first catch limit', IYRCL
  45.       READ (IN, '(T30, I4)') POUT
  46.       IF (POUT .EQ. 1) THEN
  47.         WRITE (IOUT, '(A, T30,A4)') 'Apply phaseout if necessary', 'YES'
  48.         WRITE (*,    '(A, T30,A4)') 'Apply phaseout if necessary', 'YES'
  49.       ELSE
  50.         WRITE (IOUT, '(A, T30,A4)') 'Apply phaseout', 'No'
  51.         WRITE (*,    '(A, T30,A4)') 'Apply phaseout', 'No'
  52.       ENDIF
  53.  
  54. C	Re-scale IYRCL such that 0 is the year prior to the first input data
  55.  
  56.       ISTART = 0
  57.       IYEAR = IYRCL - INITYR
  58.       IF (IYEAR .LE. 0 .OR. IYEAR .GT. MAXYR) STOP 'INVALID YEAR'
  59.  
  60. C	Initialize the catch array
  61.  
  62.       DO 10 I = 0, MAXYR
  63.         CATCH(I) = 0
  64.  10   CONTINUE
  65.  
  66. C	Read in the catch data, scaling each year to the initial year
  67.  
  68.       READ (IN, '(// A)') FORMT
  69.       WRITE (IOUT, '(/A)') 'Historic catches:'
  70.       WRITE (*,    '(/A)') 'Historic catches:'
  71.       TOTALC = 0
  72.       DO 20 I = 0, MAXYR
  73.         READ (IN, FORMT) IY, C 
  74.         IF (IY .LT. 0) GO TO 25
  75.         WRITE (IOUT, FORMT) IY, C
  76.         WRITE (*,    FORMT) IY, C
  77.         IY = IY - INITYR
  78.         IF (IY .LT. 0 .OR. IY .GE. IYEAR) STOP 
  79.         CATCH(IY) = C
  80.         TOTALC = TOTALC + C
  81.  20   CONTINUE
  82.  
  83.  25   IF (TOTALC .LT. 0) STOP ' ERROR: No historic catch input'
  84.  
  85. C	Read in the non-zero sightings estimates and information matrix
  86.  
  87.       READ (IN, '(//T30, I4 / A)') NS, FORMT
  88.       WRITE (IOUT, '(/A)') 'Abundance estimates:'
  89.       WRITE (*,    '(/A)') 'Abundance estimates:'
  90.       IF (NS .GT. MAXEST) STOP 'ERROR: Abundance year out of range'
  91.       DO 30 N=0,NS-1
  92.         N1 = (N * (N + 1)) / 2
  93.         READ (IN, FORMT) ISYR(N), SIGHT(N), (FMATRX(N1 + J), J = 0, N)
  94.         WRITE (IOUT, FORMT) ISYR(N), SIGHT(N), 
  95.      1                           (FMATRX (N1 + J), j = 0, N)
  96.         WRITE (*,    FORMT) ISYR(N), SIGHT(N),
  97.      1                           (FMATRX (N1 + J), j = 0, N)
  98.         ISYR(N) = ISYR(N) - INITYR
  99.         IF (ISYR(N) .LT. 0 .OR. ISYR(N) .GE. IYEAR) STOP 
  100.      1                         'ERROR: Sight year out of range'
  101.         IF (SIGHT(N) .LE. 0.0) STOP ' ERROR: Estimate not positive'
  102.  30   CONTINUE
  103.  
  104. C	Read in the Poisson multipliers for any zero sightings estimates
  105.  
  106.       READ (IN, '(/T30, I3, /A)') NZ, FORMAT
  107.       IF (NZ .GT. 0) THEN
  108.         WRITE (IOUT, '(/A)') 'Zero abundance estimates:'
  109.         WRITE (*,    '(/A)') 'Zero abundance estimates:'
  110.       ENDIF
  111.       IF (NZ .GT. MAXEST) STOP ' ERROR: Zero estimate array too small'
  112.       DO 40 N = 0, NZ - 1
  113.         N1 = (N * (N + 1)) / 2
  114.         READ (IN, FORMT) IZYR(N), ZMULT(N)
  115.         WRITE (IOUT, FORMT) IZYR(N), ZMULT(N)
  116.         WRITE (*,    FORMT) IZYR(N), ZMULT(N)
  117.         IZYR(N) = IZYR(N) - INITYR
  118.         IF (IZYR(N) .LT. 0 .OR. IZYR(N) .GE. IYEAR) STOP 
  119.      1                                    ' Sight year out of range'
  120.         IF (ZMULT(N) .LE. 0.0) STOP ' ERROR: Multiplier not positive'
  121.  40   CONTINUE
  122.       WRITE (IOUT, '()')
  123.       WRITE (*,    '()')
  124.  
  125. C	Bound the range for the integration over K
  126.  
  127.       RKHI = 1.E7
  128.       RKLO = 0.0
  129.  
  130. C ******
  131. C ****** Run the CLA to obtain the nominal catch limit
  132. C ******
  133.  
  134.       RAWCL = CLIMIT (CATCH, SIGHT, FMATRX, ISYR, IZYR, ZMULT, POP, G)
  135. C	Set the catch limits for PCYCLE years. If the catch limit may be
  136. C	subject to phaseout, call POUT to apply the phaseout rule.
  137.  
  138. C	First set ILAST = year of the most recent abundance estimate
  139.  
  140.       IF (NS .GT. 0) ILAST = ISYR(NS - 1)
  141.       IF (NZ .GT. 0) ILAST = MAX(IS, IZYR(NZ - 1))
  142.  
  143. C     DO 100 IY = IYEAR, IYEAR + PCYCLE - 1
  144.       DO 100 IY = IYEAR, IYEAR + int(PCYCLE) - 1
  145.         IF (POUT .EQ. 1) THEN
  146.           CALL PHOUT (CL, RAWCL, ILAST, IY)
  147.         ELSE
  148.           CL = RAWCL
  149.         ENDIF
  150.         WRITE (IOUT, '(A6, I5, A17, I6)') 'Year:', IY + INITYR,
  151.      1                                         'Catch limit:', NINT(CL)
  152.         WRITE (*   , '(A6, I5, A17, I6)') 'Year:', IY + INITYR,
  153.      1                                         'Catch limit:', NINT(CL)
  154.  100  CONTINUE
  155.  
  156.       STOP
  157.       END
  158.  
  159. C **************************************************************************************
  160. C **************************************************************************************
  161. C
  162. C	CLC version 6
  163. C
  164. C **************************************************************************************
  165. C
  166. C	31 January 1994
  167. C
  168. C **************************************************************************************
  169.  
  170.       SUBROUTINE RDPARS
  171.  
  172. C	Read the file of input parameters
  173.  
  174.       COMMON /MANPAR/ PPROB,PYMAX,PNYSTP,PKSTEP,PDSTEP,PNBSTP,PBMIN,
  175.      1		     PBMAX,PSCALE,PHASET,PHASEP,PCYCLE,PLEVEL,PSLOPE
  176.       REAL PPROB, PYMAX, PNYSTP, PKSTEP, PDSTEP, PNBSTP, PBMIN, PBMAX,
  177.      1		PSCALE, PHASET, PHASEP, PCYCLE, PLEVEL, PSLOPE
  178.       DATA IFILE /21/
  179. C
  180. C	Open the input file (only read once)
  181. C
  182.       OPEN (UNIT = IFILE, FILE = 'CLC.PAR')
  183. C
  184.       READ (IFILE, '(T30, F10.0)') PPROB,PYMAX, PNYSTP, PKSTEP, PDSTEP,
  185.      1			     PBMIN,  PBMAX,  PNBSTP, PSCALE, PHASET,
  186.      1			     PHASEP, PCYCLE, PLEVEL, PSLOPE
  187.       CLOSE (IFILE)
  188.       RETURN
  189.       END
  190. C
  191. C
  192. C **************************************************************************************
  193. C **************************************************************************************
  194. C
  195.       SUBROUTINE PHOUT (CL, RAWCL, ILAST, IY)
  196. C
  197.       COMMON /MANPAR/ PPROB,PYMAX,PNYSTP,PKSTEP,PDSTEP,PBMIN,PBMAX,
  198.      1	     PNBSTP,PSACLE,PHASET,PHASEP,PCYCLE,PLEVEL,PSLOPE
  199.       REAL PPROB, YMAX,PNYSTP,PKSTEP,PDSTEP,PBMIN,PBMAX,PNBSTP,PSACLE,
  200.      1     PHASET, PHASEP, PCYCLE, PLEVEL, PSLOPE
  201. C	PHASET Number of years without surveys before phaseout invoked
  202. C	PHASEP Phaseout annual reduction proportion
  203. C
  204.       REAL CL, RAWCL
  205.       INTEGER ILAST, IY
  206. C
  207. C	Phaseout: Reduce catch limit if there is no survey data in the
  208. C	last PHASET years
  209.       IF (IY .GE. ILAST + PHASET) THEN
  210.         CL = RAWCL * (1.0 - PHASEP * (IY - ILAST - PHASET))
  211.         IF (CL .LT. 0.0) CL = 0.0
  212.       ELSE
  213.         CL = RAWCL
  214.       ENDIF
  215. C
  216.       RETURN
  217.       END
  218.  
  219. C ***************************************************************************************
  220.  
  221.       REAL FUNCTION CLIMIT (CATCH,SIGHT,FMATRX,ISYR,IZYR,ZMULT,POP,G)
  222.  
  223. C	Run the CLA to obtain to nominal catch limit
  224. C
  225.       PARAMETER (MAXSIM = 600000, MAXSTP = 500)
  226.       REAL PRES(0:MAXSIM), QRES(0:MAXSIM), SS0, SS1, SS2, SS3
  227. C
  228.       COMMON /MANPAR/ PPROB,PYMAX,PNYSTP,PKSTEP,PDSTEP,PNBSTP,PBMIN,
  229.      1	         PBMAX,PSCALE,PHASET,PHASEP,PCYCLE,PLEVEL,PSLOPE
  230. 	REAL PPROB,PYMAX,PNYSTP,PKSTEP,PDSTEP,PNBSTP,PBMIN, PBMAX,
  231.      1		  PSCALE,PHASET,PHASEP,PCYCLE,PLEVEL,PSLOPE
  232. C
  233.       COMMON /MANDAT/ ISTART, IYEAR, NS, NZ, RKLO, RKHI
  234.       INTEGER ISTART, IYEAR, NS, NZ
  235.       REAL    RKLO, RKHI
  236. C
  237.       REAL CATCH(0:*),SIGHT(0:*),FMATRX(0:*),ZMULT(0:*),POP(0:*),G(0:*)
  238.       INTEGER ISYR(0:*), IZYR(0:*)
  239. C
  240. C	Local variables:
  241.  
  242.       REAL SF, Y(0:MAXSTP), B(0:MAXSTP), RLGB(0:MAXSTP)
  243.       INTEGER NB, NR
  244.  
  245. C	NB	Number of bias steps
  246. C	NR	Number of steps for Y (the productivity parameter)
  247. C	SF	Deviance scale factor (SF = .5 / PSCALE**2)
  248. C
  249.       CLIMIT = 0.
  250.       IF (NS .LE. 0) RETURN
  251. C
  252. C	Set deviance scale factor S = 1/PSCALE**2
  253.       SF = 0.5 / (PSCALE * PSCALE)
  254. C
  255. C	Check the sizes of the Y and B arrays are large enough
  256.       IF (PNBSTP .GT. MAXSTP .OR. PNYSTP .GT. MAXSTP) STOP
  257.      1                             'Y &/or b array sizes too small'
  258.  
  259. C	Set sightings bias step sizes & their log values. BINC = Bias increment
  260.       NB = PNBSTP
  261.       BINC = (PBMAX - PBMIN) / NB
  262.       DO 10 I = 0, NB - 1
  263.         B(I) = PBMIN + (I + 0.5) * BINC
  264.         RLGB(I) = -ALOG(B(I))
  265.  10   CONTINUE
  266. C
  267. C	Set productivity parameter step sizes (midpoints)
  268.       NR = PNYSTP
  269.       YINC = PYMAX / NR
  270.       DO 20 I = 0, NR - 1
  271.         Y(I) = (I + 0.5) * YINC
  272.  20   CONTINUE
  273. C
  274.  
  275.       PTOT = 0
  276.       N = 0
  277.       DO 50 I = 0, NR - 1
  278.  
  279. C	  Set R from the productivity parameter, Y
  280.         R = 1.4184 * Y(I)
  281.  
  282. C	  Step size for K
  283.         DK = PKSTEP
  284.         D = 1.
  285.         RK = RKHI
  286.  
  287. C	  Use function STKSIM to set up the Nth population trajectory
  288. C	  i.e. set up the pop array
  289.  30     IF (RK .LE. RKLO .OR. STKSIM (RK, R, POP, CATCH) .LE. 0.)
  290.      1                                                     GOTO 40
  291. 	  IF (N .GE. MAXSIM) STOP 'ERROR: TOO MANY SIMULATIONS'
  292.  
  293. C	    How much depletion covered?
  294.           DD = D - POP(IYEAR) / RK
  295.           D = POP(IYEAR) / RK
  296.           P = 0.0
  297.  
  298.           IF (DD .GT. 0.) THEN
  299.  
  300. C	      Compute the internal catch limit corresponding to D and Y(I)
  301.             QRES(N) = CONTRL (D, Y(I), PLEVEL, PSLOPE) * POP(IYEAR)
  302.  
  303. C	      Calculate deviance
  304.             CALL DEVIAN (SS0, SS1, SS2, SS3, SIGHT, FMATRX, ISYR,
  305.      1                                            IZYR, ZMULT, POP, G)
  306.  
  307. C	      Scale likelihood and integrate over values for the bias parameter
  308.             DO 35 J = 0, NB - 1
  309.               P = P + EXP(-SF * (SS0 + RLGB(J) * (SS1 + RLGB(J) * SS2)
  310.      1                                              + SS3 * B(J)))
  311.  
  312.  35         CONTINUE
  313.  
  314.  
  315. C	      Calculate the weight for this point (& total area under likelihood)
  316.             PRES(N) = P * DD
  317.             PTOT = PTOT + PRES(N)
  318.  
  319. C	      Update counter
  320.             N = N + 1
  321.  
  322. C	      Find the next K
  323.             DK = DK * PDSTEP / DD
  324.             IF (DK .GT. PKSTEP) DK = PKSTEP
  325.  
  326.             ELSE
  327. C	    IF DD = 0 change the step size only 
  328.               DK = PKSTEP
  329.             ENDIF
  330.  
  331. C	    Set the new value of K
  332.           RK = RK / (1. + DK)
  333.  
  334.         GOTO 30
  335.  40     CONTINUE
  336.  50   CONTINUE
  337.  
  338.  
  339.       IF (PTOT .LE. 0.) STOP 'ERROR: PROB INTEGRATES TO ZERO'
  340.  
  341. C	Sort the QRES and PRES arrays in ascending order of QRES
  342.       N2 = N
  343.       CALL SORT (QRES, PRES, N2)
  344.  
  345. C	Normalize the relative likelihoods
  346.       DO 60 I = 0, N - 1
  347.         PRES(I) = PRES(I) / PTOT
  348.  60   CONTINUE
  349.  
  350. C	Extract the desired probability level: the nominal catch limit (NCL)
  351. C	is the lower PROB% of the distribution.
  352. C	First calculate PRES(I), the probability that the NCL is between
  353. C	   QRES(I) & QRES(I + 1).
  354.       P = 0
  355.       DO 70 I = 0, N - 1
  356.         P = P + PRES (I)
  357.         IF (P .GT. PPROB) GOTO 80
  358.  70   CONTINUE
  359.  
  360. C	Interpolate to set the nominal catch limit
  361.  80   IF (I .GE. N - 1) THEN
  362.         Q = QRES (N - 1)
  363.       ELSE
  364.         Q = (QRES (I + 1) * (PPROB - P + PRES(I)) + QRES(I) *
  365.      1                                    (P - PPROB)) / PRES(I)
  366.       ENDIF
  367.  
  368.       CLIMIT = Q
  369.  
  370.       RETURN
  371.       END
  372.  
  373. C ***********************************************************************************
  374.  
  375.       FUNCTION CONTRL (D, Y, PLEVEL, PSLOPE)
  376.  
  377. C	Catch control law
  378.  
  379.       IF (D .LT. PLEVEL) THEN
  380.         CONTRL = 0.
  381.       ELSEIF (D .LT. 1.) THEN
  382.         CONTRL = PSLOPE * Y * (D - PLEVEL)
  383.       ELSE
  384.         CONTRL = PSLOPE * Y * (1.0 - PLEVEL)
  385.       ENDIF
  386.  
  387.       END
  388.  
  389. C **********************************************************************************
  390.  
  391.       SUBROUTINE DEVIAN (SS0, SS1, SS2, SS3, SIGHT, FMATRX, ISYR,
  392.      1                                         IZYR, ZMULT, POP, G)
  393.  
  394. C	Calculate deviance (-2 log likelihood) in terms of coefficients for 
  395. C	the bias and log bias parameters
  396.  
  397.       COMMON /MANDAT/ ISTART, IYEAR, NS, NZ, RKLO, RKHI
  398.       INTEGER ISTART, IYAER, NS, NZ
  399.       REAL    RKLO, RKHI
  400.  
  401.       REAL SS0,SS1,SS2,SIGHT(0:*),FMATRX(0:*),ZMULT(0:*),POP(0:*),G(0:*)
  402.       INTEGER ISYR(0:*), IZYR(0:*)
  403.  
  404.       SS0 = 0
  405.       SS1 = 0
  406.       SS2 = 0
  407.       SS3 = 0
  408.  
  409.       DO 100 N = 0, NS - 1
  410.  
  411.         G(N) = ALOG (SIGHT(N) / POP(ISYR(N)))
  412.         K = N * (N + 1) / 2
  413.         DO 10 J = 0, N - 1
  414. C	    1st add non diagonal contributions (which are doubled up)
  415.           SS0 = SS0 + 2. * G(J) * G(N) * FMATRX(K)
  416.           SS1 = SS1 + 2. * (G(J) + G(N)) * FMATRX(K)
  417.           SS2 = SS2 + FMATRX(K) + FMATRX(K)
  418.           K = K + 1
  419.  10     CONTINUE
  420. C	  Now add diagnoal contribution
  421.         SS0 = SS0 + G(N) * G(N) * FMATRX(K)
  422.         SS1 = SS1 + 2. * G(N) * FMATRX(K)
  423.         SS2 = SS2 + FMATRX(K)
  424.  100  CONTINUE
  425.  
  426. C	Now do the zero estimates
  427.       DO 200 N = 0, NZ - 1
  428.         SS3 = SS3 + 2.0 * POP(IZYR(N)) / ZMULT(N)
  429.  200  CONTINUE
  430.  
  431.       RETURN
  432.       END
  433.  
  434. C *********************************************************************************
  435.  
  436.       FUNCTION STKSIM (RK, R, POP, CATCH)
  437.  
  438. C	Calculate the stock trajectory with parameter RK and R: return
  439. C	the current stock size
  440. C	RK = notional carrying capacity; R = productivity parameter * 1.4184
  441.  
  442.       COMMON /MANDAT/ ISTART, IYEAR, NS, NZ, RKLO, RKHI
  443.       INTEGER ISTART, IYEAR, NS, NZ
  444.       REAL    RKLO, RKHI
  445.  
  446.       REAL CATCH(0:*), POP(0:*)
  447.  
  448.       POP(ISTART) = RK
  449.       DO 10 I = ISTART + 1, IYEAR
  450.         D = POP(I - 1) / RK
  451.         POP(I) = POP (I - 1) * (1. + R * (1. - D * D)) - CATCH(I-1)
  452.         IF (POP(I) .LE. 0) GOTO 20
  453.  10   CONTINUE
  454.       STKSIM = POP(IYEAR)
  455.       RETURN
  456.  
  457.  20   STKSIM = 0.
  458.  
  459.       END
  460.  
  461. C *********************************************************************************
  462.  
  463.       SUBROUTINE SORT (ARRAY, ARRAY2, N)
  464.  
  465. C	SORT sorts a pair of arrays in ascending order of the first array
  466. C	(using the heapsort algorithm)
  467.  
  468.       REAL ARRAY(0:*), ARRAY2(0:*), TEMP, TEMP2
  469.       INTEGER N, K, IR, I, J
  470.       IF (N .LT. 2) RETURN
  471.       K = N / 2
  472.       IR = N - 1
  473.  10   IF (K .NE. 0) THEN
  474.         K = K - 1
  475.         TEMP = ARRAY(K)
  476.         TEMP2 = ARRAY2(K)
  477.       ELSE
  478. 	TEMP = ARRAY(IR)
  479. 	TEMP2 = ARRAY2(IR)
  480.         ARRAY(IR) = ARRAY(0)
  481.         ARRAY2(IR) = ARRAY2(0)
  482.         IR = IR - 1
  483.         IF (IR .EQ. 0) THEN
  484.           ARRAY(0) = TEMP
  485.           ARRAY2(0) = TEMP2
  486.           RETURN
  487.         ENDIF
  488.       ENDIF
  489.       I = K
  490.       J = K + K + 1
  491.  20   IF (J .LE. IR) THEN
  492.         IF (J .LT. IR .AND. ARRAY(J) .LT. ARRAY(J + 1)) J = J + 1
  493.         IF (TEMP .LT. ARRAY(J)) THEN
  494.           ARRAY(I) = ARRAY(J)
  495.           ARRAY2(I) = ARRAY2(J)
  496.           I = J
  497.           J = J + I + 1
  498.         ELSE
  499.           J = IR + 1
  500.         ENDIF
  501.         GOTO 20
  502.       ENDIF
  503.       ARRAY(I) = TEMP
  504.       ARRAY2(I) = TEMP2
  505.       GOTO 10
  506.  
  507.       END
Not running #stdin #stdout 0s 0KB
comments (?)
stdin 
Standard input is empty
stdout 
Standard output is empty
https://ideone.com/1HJaWg
language:
Text (text 6.10)
created:
5 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!