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