text:
c-------------------------------------------------------------------------------
c-------------------------------------------------------------------------------
c     SUBROUTINE REG_AR1.
c
c     Date:    27-MAR-2000.
c     Author:  B.D. Santer, PCMDI, Livermore.
c     Version: Single Precision.
c
c     Performs a simple linear regression with dependent variable y and
c     independent variable x. Various other statistical quantities are
c     also computed. Also computes lag-1 autocorrelation of residuals from
c     linear regression, and uses this to adjust effective temporal degrees
c     of freedom for confidence interval about estimated regression slope.
c
c     REFERENCE:
c
c     NAGFLIB:801/715:Mk5:Nov74.
c     Routine G02CAF.
c
c     EXTERNALS:
c
c     function BETAI   - Computes incomplete Beta function.
c     function BETACF  - Computes continued fraction for BETAI.
c     function GAMMLN  - Computes Gamma function
c     function PROBF1  - Computes probability for given F-value.
c     function PROBND1 - Computes area under a normal curve.
c     subroutine AUTO1 - Computes autocovariance and autocorrelation.
c
c
      subroutine reg_ar1(nmax, n1, lag, rneff, x, y,
     &result, res, cxx, rxx)
c
c
c
      implicit real (o-z, a-h)
c
c     ** INPUT **
c     integer nmax          Max. length of x and y
c     integer n1            Actual length of x and y
c     integer lag           Max. lag for autocorrelations
c     dimension x(nmax)     Independent variable
c     dimension y(nmax)     Dependent variable
c
c     ** OUTPUT **
c     real rneff               Effective sample size
c     dimension result(31)     Array of linear regression results
c     dimension res(nmax)      Residuals from linear regression
c     dimension cxx(0 : lag)   Autocovariance function
c     dimension rxx(0 : lag)   Autocorrelation function
c
      integer nmax, n1, lag
      real rneff
      dimension x(nmax), y(nmax), result(31), res(nmax)
      dimension cxx(0 : lag), rxx(0 : lag)
c
c
c
c-------------------------------------------------------------------------------
c
c
c
      eps1  = 1.0e-8
      bogus = -99.999
c
c
c
c-------------------------------------------------------------------------------
c     ** 1.0 **
c
c     Initialize sums.
c-------------------------------------------------------------------------------
c
c
c
      sumx  = 0.0
      sumy  = 0.0
      varx  = 0.0
      vary  = 0.0
      xy    = 0.0
      ssd   = 0.0
c
      do 10 i = 1, n1
        sumx  = sumx + x(i)
        sumy  = sumy + y(i)
   10 continue
c
c
c
c-------------------------------------------------------------------------------
c     ** 1.1 **
c
c     Compute means.
c-------------------------------------------------------------------------------
c
c
c
      dn   = float(n1)
      xbar = sumx / dn
      ybar = sumy / dn
c
c
c
c-------------------------------------------------------------------------------
c     ** 1.2 **
c
c     Compute standard deviations (sdx, sdy).
c-------------------------------------------------------------------------------
c
c
c
      do 12 i = 1, n1
       x1   = x(i) - xbar
       y1   = y(i) - ybar
       xy   = xy + (x1 * y1)
       varx = varx + (x1 * x1)
       vary = vary + (y1 * y1)
   12 continue
c
      if (varx .le. eps1) then
        print *, 'In REG_AR1. Small or zero variance detected in x!'
        print *, 'varx = ', varx
        do 121 kk = 1, 31
          result(kk) = bogus
  121   continue
        goto 1000
      else
        sdx  = sqrt(varx / (dn - 1.0))
      end if
c
      if (vary .eq. 0.0) then
        print *, 'In REG_AR1. Small or zero variance detected in y!'
        print *, 'vary = ', vary
        do 122 kk = 1, 31
          result(kk) = bogus
  122   continue
        goto 1000
      else
        sdy  = sqrt(vary / (dn - 1.0))
      end if
c
c
c
c-------------------------------------------------------------------------------
c     ** 1.3 **
c
c     Compute Pearson product-moment correlation coefficient (r).
c-------------------------------------------------------------------------------
c
c
c
      r = xy / sqrt(varx * vary)
c
c
c
c-------------------------------------------------------------------------------
c     ** 1.4 **
c
c     Compute regression coefficient (b) and regression constant (a).
c-------------------------------------------------------------------------------
c
c
c
      b    = xy / varx
      a    = ybar - b * xbar
c
c
c
c-------------------------------------------------------------------------------
c     ** 1.5 **
c
c     Compute:
c     ssr       Sum of squares attributable to the regression.
c     ssd       Sum of squares of deviations about the regression.
c     sst       Total sum of squares.
c     res       Residuals from regression.
c     rxx       Autocorrelation function.
c     rneff     Effective sample size.
c-------------------------------------------------------------------------------
c
c
c
      sst = vary
      do 15 i = 1, n1
        x2     = y(i) - a - b * x(i)
        ssd    = ssd + (x2 * x2)
        res(i) = x2
   15 continue
c
      ssr = sst - ssd
c
      call auto1(nmax, n1, lag, res, cxx, rxx)
c
      term1 = 1.0 - rxx(1)
      term2 = 1.0 + rxx(1)
      term3 = n1 * 1.0
      rneff = term3 * (term1 / term2)
c
c
c
c-------------------------------------------------------------------------------
c     ** 1.6 **
c
c     Compute:
c     idfr      D.O.F. attributable to the regression.
c     idfd1     D.O.F. of deviations about the regression (unadjusted).
c     rdfd2     D.O.F. of deviations about the regression (adjusted).
c     idft      Total degrees of freedom.
c-------------------------------------------------------------------------------
c
c
c
      idft  = n1 - 1
      idfd1 = n1 - 2
      rdfd2 = rneff - 2.0
      idfr  = 1
c
c
c
c-------------------------------------------------------------------------------
c     ** 1.7 **
c
c     Compute:
c     amsr      Mean square attributable to the regression.
c     amsd1     Mean square of deviations about the regression (unadjusted).
c     amsd2     Mean square of deviations about the regression (adjusted).
c-------------------------------------------------------------------------------
c
c
c
      amsr  = ssr / float(idfr)
      amsd1 = ssd / float(idfd1)
      amsd2 = ssd / rdfd2
c
c
c
c-------------------------------------------------------------------------------
c     ** 1.8 **
c
c     Compute the F-value for the analysis of variance (f). See also D.S.
c     Wilks (Statistical Methods in the Atmospheric Sciences, Academic
c     Press, 1995, pp. 166-168).
c-------------------------------------------------------------------------------
c
c
c
      f = amsr / amsd1
c
c
c
c-------------------------------------------------------------------------------
c     ** 1.9 **
c
c     Compute:
c     seb1      Standard error of regression coefficient b (unadjusted).
c     seb2      Standard error of regression coefficient b (adjusted).
c     sea1      Standard error of regression constant a (unadjusted).
c     sea2      Standard error of regression constant a (adjusted).
c-------------------------------------------------------------------------------
c
c
c
      seb1 = sqrt(amsd1 / varx)
      seb2 = sqrt(amsd2 / varx)
      s1   = (1.0 / dn) + ((xbar * xbar) / varx)
      sea1 = sqrt(amsd1 * s1)
      sea2 = sqrt(amsd2 * s1)
c
c
c
c-------------------------------------------------------------------------------
c     ** 1.10 **
c
c     Compute:
c     tb1       t-value for regression coefficient b (unadjusted).
c     tb2       t-value for regression coefficient b (adjusted).
c     ta1       t-value for regression constant a (unadjusted).
c     ta2       t-value for regression constant a (adjusted).
c-------------------------------------------------------------------------------
c
c
c
      tb1 = b / seb1
      tb2 = b / seb2
      ta1 = a / sea1
      ta2 = a / sea2
c
c
c
c-------------------------------------------------------------------------------
c     ** 1.11 **
c
c     Compute:
c     pt1       p-value for regression coefficient t-value. Effective sample
c               size adjustment for standard error (seb).
c     pt2       p-value for regression coefficient t-value. Effective sample
c               size adjustment for standard error (seb) and critical t-value.
c     pt3       p-value for regression coefficient t-value. No adjustment for
c               standard error or critical t-value.
c     pf1       p-value for regression coefficient F-value (one-tailed).
c     pf2       p-value for regression coefficient F-value (two-tailed).
c
c     The values pt1 and pt2 are use to test the null hypothesis that b = 0
c     (i.e., y is independent of x).
c
c     The values pf1 and pf2 are use to test the null hypothesis that the
c     regression is linear (goodness of linear fit). For non-replicated
c     values of y, the degrees of freedom are 1 and n - 2.
c-------------------------------------------------------------------------------
c
c
c
      ii1 = 1
      xx1 = n1 / (n1 + (tb2 * tb2))
      xx2 = rneff / (rneff + (tb2 * tb2))
      xx3 = n1 / (n1 + (tb1 * tb1))
      aa1 = n1 / 2.0
      aa2 = rneff / 2.0
      bb  = 1.0 / 2.0
c
      pt1 = betai(aa1, bb, xx1)
      pt2 = betai(aa2, bb, xx2)
      pt3 = betai(aa1, bb, xx3)
      pf1 = probf1(f, dum3, ii1, n1 - 2, 1)
      pf2 = probf1(f, dum4, ii1, n1 - 2, 2)
c
c
c
c-------------------------------------------------------------------------------
c     ** 1.12 **
c
c     Store results (result).
c-------------------------------------------------------------------------------
c
c
c
      result(1)  = xbar
      result(2)  = ybar
      result(3)  = sdx
      result(4)  = sdy
      result(5)  = r
      result(6)  = b
      result(7)  = a
      result(8)  = seb1
      result(9)  = seb2
      result(10) = sea1
      result(11) = sea2
      result(12) = tb1
      result(13) = tb2
      result(14) = ta1
      result(15) = ta2
      result(16) = ssr
      result(17) = idfr * 1.0
      result(18) = amsr
      result(19) = f
      result(20) = ssd
      result(21) = idfd1 * 1.0
      result(22) = rdfd2
      result(23) = amsd1
      result(24) = amsd2
      result(25) = sst
      result(26) = idft * 1.0
      result(27) = pt1
      result(28) = pt2
      result(29) = pt3
      result(30) = pf1
      result(31) = pf2
c
 1000 continue
      return
      end
c-------------------------------------------------------------------------------
c-------------------------------------------------------------------------------
      FUNCTION betai(a,b,x)
      REAL betai,a,b,x
CU    USES betacf,gammln
      REAL bt,betacf,gammln
      if(x.lt.0..or.x.gt.1.)pause 'bad argument x in betai'
      if(x.eq.0..or.x.eq.1.)then
        bt=0.
      else
        bt=exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.-x))
      endif
      if(x.lt.(a+1.)/(a+b+2.))then
        betai=bt*betacf(a,b,x)/a
        return
      else
        betai=1.-bt*betacf(b,a,1.-x)/b
        return
      endif
      END
c-------------------------------------------------------------------------------
c-------------------------------------------------------------------------------
      FUNCTION betacf(a,b,x)
      INTEGER MAXIT
      REAL betacf,a,b,x,EPS,FPMIN
      PARAMETER (MAXIT=100,EPS=3.e-7,FPMIN=1.e-30)
      INTEGER m,m2
      REAL aa,c,d,del,h,qab,qam,qap
      qab=a+b
      qap=a+1.
      qam=a-1.
      c=1.
      d=1.-qab*x/qap
      if(abs(d).lt.FPMIN)d=FPMIN
      d=1./d
      h=d
      do 11 m=1,MAXIT
        m2=2*m
        aa=m*(b-m)*x/((qam+m2)*(a+m2))
        d=1.+aa*d
        if(abs(d).lt.FPMIN)d=FPMIN
        c=1.+aa/c
        if(abs(c).lt.FPMIN)c=FPMIN
        d=1./d
        h=h*d*c
        aa=-(a+m)*(qab+m)*x/((a+m2)*(qap+m2))
        d=1.+aa*d
        if(abs(d).lt.FPMIN)d=FPMIN
        c=1.+aa/c
        if(abs(c).lt.FPMIN)c=FPMIN
        d=1./d
        del=d*c
        h=h*del
        if(abs(del-1.).lt.EPS)goto 1
11    continue
      pause 'a or b too big, or MAXIT too small in betacf'
1     betacf=h
      return
      END
c-------------------------------------------------------------------------------
c-------------------------------------------------------------------------------
      FUNCTION gammln(xx)
      REAL gammln,xx
      INTEGER j
      DOUBLE PRECISION ser,stp,tmp,x,y,cof(6)
      SAVE cof,stp
      DATA cof,stp/76.18009172947146d0,-86.50532032941677d0,
     *24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2,
     *-.5395239384953d-5,2.5066282746310005d0/
      x=xx
      y=x
      tmp=x+5.5d0
      tmp=(x+0.5d0)*log(tmp)-tmp
      ser=1.000000000190015d0
      do 11 j=1,6
        y=y+1.d0
        ser=ser+cof(j)/y
11    continue
      gammln=tmp+log(stp*ser/x)
      return
      END
c-------------------------------------------------------------------------------
c-------------------------------------------------------------------------------
c     FUNCTION PROBF1.
c
c     Date:    6-MAR-1992.
c     Author:  B.D. Santer, MPI, Hamburg.
c     Version: Single Precision.
c
c     The output is either the one- or two-tailed test area: i.e., the
c     area under an F-curve (with N1 and N2 degrees of freedom) to the
c     right of X if X exceeds 1.0 (one-tailed test) or twice this area
c     (two-tailed test).
c
c     Note: if X is less than 1.0, this function gives the area to the
c     right of 1/X with reversed order for the degrees of freedom. This
c     ensures the accuracy of the numerical algorithm.
c
c     REFERENCE:
c
c     M. Abramowitz and I.A. Stegun.
c     Handbook of Mathematical Functions.
c     Dover, 1970, page 947 (26.6.15).
c
c     EXTERNALS:
c
c     function PROBND1 - Calculates the area under a normal curve.
c
c
c
      function probf1(y, x, n1, n2, id)
c
c
c
      implicit real (o-z, a-h)
c
c     ** INPUT **
c     real y            Calculated F-value (unchanged on output)
c     real x            Inverse of Y if Y is less than 1.0
c     integer n1, n2    Degrees of freedom
c     integer id        Identifier for one- or two-tailed test
c
c     ** OUTPUT **
c     real probf1       Significance level (p-value) for F-value
c
c
c
      integer n1, n2, id
      real y, x, probf1
c
c
c
c-------------------------------------------------------------------------------
c
c
c
      x = y
      if (x .gt. 1.0) goto 1
      x  = 1.0 / x
      n  = n1
      n1 = n2
      n2 = n
c
    1 continue
      dn1 = float(n1)
      dn2 = float(n2)
      term1 = 2.0 / (9.0 * dn1)
      term2 = 2.0 / (9.0 * dn2)
      term3 = ((x ** (1.0 / 3.0)) * (1.0 - term2)) - (1.0 - term1)
      term4 = sqrt(term1 + ((x ** (2.0 / 3.0)) * term2))
      term5 = term3 / term4
      if (id .eq. 2) probf1 = 2.0 * probnd1(term5)
c
c
c
c-------------------------------------------------------------------------------
c     The numerical algorithm can have problems when the F-value is
c     close to 1.0 and the degrees of freedom are small. Therefore,
c     insure that the probabilities returned cannot exceed 1.0.
c-------------------------------------------------------------------------------
c
c
c
      if (id .eq. 1) probf1 = probnd1(term5)
      if (probf1 .gt. 1.0) probf1 = 1.0
      return
      end
c-------------------------------------------------------------------------------
c-------------------------------------------------------------------------------
c     FUNCTION PROBND1.
c
c     Date:    6-MAR-1992.
c     Author:  B.D. Santer, MPI, Hamburg.
c     Version: Single Precision.
c
c     Calculates the area under a normal curve (mean=0.0, variance=1.0)
c     to the right of x. The accuracy is better than 7.5 * 10.**-8.
c
c     REFERENCE:
c
c     M. Abramowitz and I.A. Stegun.
c     Handbook of Mathematical Functions.
c     Dover, 1970, pages 931-932 (26.2.1 and 26.2.17).
c
c
c
      function probnd1(x)
c
c
c
      implicit real (o-z, a-h)
c
c     ** INPUT **
c     real x         Input z score
c
c     ** OUTPUT **
c     real probnd1   Area under normal curve to right of x
c
c
c
      real x, probnd1
c
c
c
c-------------------------------------------------------------------------------
c
c
c
      b1 =  0.319381530
      b2 = -0.356563782
      b3 =  1.781477937
      b4 = -1.821255978
      b5 =  1.330274429
      p  =  0.2316419
      pi =  4.0 * atan(1.0)
c
      t = 1.0 / (1.0 + (p * x))
      term1 = ((((b1 * t) + (b2 * (t ** 2))) + (b3 * (t ** 3))) + (b4 *
     &(t ** 4))) + (b5 * (t ** 5))
c
c
c
c-------------------------------------------------------------------------------
c     If x is large, set probnd1 to zero to avoid error in computing
c     z (underflow).
c-------------------------------------------------------------------------------
c
c
c
      probnd1 = 0.0
      if (x .gt. 7.0) goto 1
      z = (1.0 / sqrt(2.0 * pi)) * exp(- ((x * x) / 2.0))
      probnd1 = z * term1
    1 continue
      return
      end
c-------------------------------------------------------------------------------
c-------------------------------------------------------------------------------
c     SUBROUTINE AUTO1.
c
c     Date:    27-OCT-1994.
c     Author:  B.D. Santer, LLNL, Livermore.
c     Version: Single precision.
c
c     Computes autocovariance (cxx) and autocorrelation (rxx) functions
c     of an input data series x(t). User specifies the number of lags
c     computed.
c
c     REFERENCE:
c
c     G.M. Jenkins and D.G. Watts.
c     Spectral Analysis and its Applications.
c     Holden-Day, 1968, pages 180-182, (5.3.25) and (5.3.33).
c
c
c
      subroutine auto1(n, n1, lag, x, cxx, rxx)
c
c
c
      implicit real (o-z, a-h)
c
c     ** INPUT **
c     integer n           Maximum dimension of x
c     integer n1          Actual dimension of x
c     integer lag         Number of lags computed
c     dimension x(n)      Input data series
c
c     ** OUTPUT **
c     dimension cxx(0 : lag)    Autocovariance function
c     dimension rxx(0 : lag)    Autocovariance function
c
      integer n, n1
      dimension x(n), cxx(0 : lag), rxx(0 : lag)
c
c
c
c-------------------------------------------------------------------------------
c
c
c
c-------------------------------------------------------------------------------
c     ** 1.0 **
c
c     Compute sample mean of x.
c-------------------------------------------------------------------------------
c
c
c
      sum1 = 0.0
      sum2 = 0.0
c
      do 100 it = 1, n1
        sum1 = sum1 + x(it)
  100 continue
      xbar = sum1 / float(n1)
c
c
c
c-------------------------------------------------------------------------------
c     ** 1.1 **
c
c     Compute total variance of x.
c-------------------------------------------------------------------------------
c
c
c
      do 110 it = 1, n1
        sum2 = sum2 + ((x(it) - xbar) ** 2)
  110 continue
      c0 = sum2 / float(n1)
c
c
c
c-------------------------------------------------------------------------------
c     ** 1.2 **
c
c     Compute autocovariance and autocorrelation for each lag (0 to lag).
c-------------------------------------------------------------------------------
c
c
c
      do 121 k = 0, lag
        sum3 = 0.0
        n2   = n1 - k
        do 120 it = 1, n2
          sum3 = sum3 + ((x(it) - xbar) * (x(it + k) - xbar))
  120   continue
c
        cxx(k) = sum3 / float(n1)
        if (c0 .eq. 0.0) then
          print *, 'AUTO1. zero variance for c0! '
          rxx(k) = 1.0
        else
          rxx(k) = cxx(k) / c0
        end if
c
c       if (k .eq. 0) write(*, 990)
c       write(*,991) k, cxx(k), rxx(k)
  121 continue
c
  990 format(t2, 'Lag', t10, 'Autocovariance', t30,
     &'Autocorrelation' )
  991 format(i3, t10, f12.5, t30, f12.5)
c
      return
      end

ext: f

max_stars_repo_name: xylar/cdat

max_issues_repo_name: xylar/cdat

max_forks_repo_name: CDAT/uvcdat


text:
! { dg-do run }
!
! PR 45420: [OOP] polymorphic TBP call in a CLASS DEFAULT clause
!
! Contributed by Salvatore Filippone <salvatore.filippone@uniroma2.it>

module base_mat_mod

 type  :: base_sparse_mat
 contains
   procedure, pass(a) :: get_fmt => base_get_fmt
 end type base_sparse_mat

contains

 function base_get_fmt(a) result(res)
   implicit none
   class(base_sparse_mat), intent(in) :: a
   character(len=5) :: res
   res = 'NULL'
 end function base_get_fmt

end module base_mat_mod

module d_base_mat_mod

 use base_mat_mod

 type, extends(base_sparse_mat) :: d_base_sparse_mat
 contains
   procedure, pass(a) :: get_fmt => d_base_get_fmt
 end type d_base_sparse_mat

 type, extends(d_base_sparse_mat) :: x_base_sparse_mat
 contains
   procedure, pass(a) :: get_fmt => x_base_get_fmt
 end type x_base_sparse_mat

contains

 function d_base_get_fmt(a) result(res)
   implicit none
   class(d_base_sparse_mat), intent(in) :: a
   character(len=5) :: res
   res = 'DBASE'
 end function d_base_get_fmt

 function x_base_get_fmt(a) result(res)
   implicit none
   class(x_base_sparse_mat), intent(in) :: a
   character(len=5) :: res
   res = 'XBASE'
 end function x_base_get_fmt

end module d_base_mat_mod

program bug20
  use d_base_mat_mod
  class(d_base_sparse_mat), allocatable  :: a

  allocate(x_base_sparse_mat :: a)
  if (a%get_fmt()/="XBASE") call abort()

  select type(a)
  type is (d_base_sparse_mat)
    call abort()
  class default
    if (a%get_fmt()/="XBASE") call abort()
  end select

end program bug20

! { dg-final { cleanup-modules "base_mat_mod d_base_mat_mod" } }

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


PROGRAM insertion_sort
!
!  Purpose:
!    To read a series of integer values from an input data file
!    and sort them using an insertion sort.  After the values
!    are sorted, they will be written back to the standard
!    output device.
!
!  Record of revisions:
!      Date       Programmer          Description of change
!      ====       ==========          =====================
!    12/23/06    S. J. Chapman        Original code
!
IMPLICIT NONE

! Derived data type to store integer values in
TYPE :: int_value
   INTEGER :: value
   TYPE (int_value), POINTER :: next_value
END TYPE

! Data dictionary: declare variable types & definitions
TYPE (int_value), POINTER :: head   ! Pointer to head of list
CHARACTER(len=20) :: filename       ! Input data file name
INTEGER :: istat                    ! Status: 0 for success
INTEGER :: nvals = 0                ! Number of data read
TYPE (int_value), POINTER :: ptr    ! Ptr to new value
TYPE (int_value), POINTER :: ptr1   ! Temp ptr for search
TYPE (int_value), POINTER :: ptr2   ! Temp ptr for search
TYPE (int_value), POINTER :: tail   ! Pointer to tail of list
INTEGER :: temp                     ! Temporary variable

! Get the name of the file containing the input data.
WRITE (*,*) 'Enter the file name with the data to be sorted: '
READ (*,'(A20)') filename

! Open input data file.
OPEN ( UNIT=9, FILE=filename, STATUS='OLD', ACTION='READ', &
       IOSTAT=istat )

! Was the OPEN successful?
fileopen: IF ( istat == 0 ) THEN         ! Open successful

   ! The file was opened successfully, so read the data value
   ! to sort, allocate a variable for it, and locate the proper
   ! point to insert the new value into the list.
   input: DO
      READ (9, *, IOSTAT=istat) temp      ! Get value
      IF ( istat /= 0 ) EXIT input        ! Exit on end of data
      nvals = nvals + 1                   ! Bump count

      ALLOCATE (ptr,STAT=istat)           ! Allocate space
      ptr%value = temp                    ! Store number

      ! Now find out where to put it in the list.
      new: IF (.NOT. ASSOCIATED(head)) THEN ! No values in list
         head => ptr                      ! Place at front
         tail => head                     ! Tail pts to new value
         NULLIFY (ptr%next_value)         ! Nullify next ptr
      ELSE
         ! Values already in list.  Check for location.
         front: IF ( ptr%value < head%value ) THEN
            ! Add at front of list
            ptr%next_value => head
            head => ptr
         ELSE IF ( ptr%value >= tail%value ) THEN
            ! Add at end of list
            tail%next_value => ptr
            tail => ptr
            NULLIFY ( tail%next_value )
         ELSE
            ! Find place to add value
            ptr1 => head
            ptr2 => ptr1%next_value
            search: DO
               IF ( (ptr%value >= ptr1%value) .AND. &
                    (ptr%value < ptr2%value) ) THEN
                  ! Insert value here
                  ptr%next_value => ptr2
                  ptr1%next_value => ptr
                  EXIT search
               END IF
               ptr1 => ptr2
               ptr2 => ptr2%next_value
            END DO search
         END IF front
      END IF new
   END DO input

   ! Now, write out the data.
   ptr => head
   output: DO
      IF ( .NOT. ASSOCIATED(ptr) ) EXIT   ! Pointer valid?
      WRITE (*,'(1X,I10)') ptr%value      ! Yes: Write value
      ptr => ptr%next_value               ! Get next pointer
   END DO output

ELSE fileopen

   ! Else file open failed.  Tell user.
   WRITE (*,'(1X,A,I6)') 'File open failed--status = ', istat

END IF fileopen

END PROGRAM insertion_sort


MODULE saprc99_mosaic_8bin_vbs2_Model

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!  Completely defines the model saprc99_mosaic_8bin_vbs2
!    by using all the associated modules
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  USE saprc99_mosaic_8bin_vbs2_Precision
  USE saprc99_mosaic_8bin_vbs2_Parameters
  USE saprc99_mosaic_8bin_vbs2_Global
  USE saprc99_mosaic_8bin_vbs2_Function
  USE saprc99_mosaic_8bin_vbs2_Integrator
  USE saprc99_mosaic_8bin_vbs2_Rates
  USE saprc99_mosaic_8bin_vbs2_Jacobian
  USE saprc99_mosaic_8bin_vbs2_LinearAlgebra
  USE saprc99_mosaic_8bin_vbs2_Monitor
  USE saprc99_mosaic_8bin_vbs2_Util

END MODULE saprc99_mosaic_8bin_vbs2_Model


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

!
!  *********************************************************************
!  *         SUBROUTINE TO CALCULATE GRAZING RATE CONSTANTS            *
!  *                  MULTIPLE SPECIES                                 *
!  *********************************************************************
!
      SUBROUTINE CONSB2(X,GRAMOR,GRADET,ZOOD,ZOOC,ITNUM,LGRAZ,T,
     1                  PERIOD,GRAMX,LCOUPL,ZMAX)

!**** Name   Type Size   I/O Description
!**** ------ ---- ------ --- -------------------------------------------------
!**** AA     R*8  NUNUCO*I   BLOOM algae stoichiometry matrix excluding
!****             NUSPEC     carbon
!**** ALGFLX R*8  NUSPEC Total algal flux (mg C/l.period) [POSITIVE = DECREASE ]
!**** ALGSTC R*8  NTONUT*    Algae stoichiometry matrix including carbon
!****             NUSPEC
!**** ALTFLX R*8  NTONUT Total algal flux (mg/l.d) [POSITIVE = ALGAE DECREASE ]
!**** BOTFLX R*8  NTONUT Total bottom detritus flux (mg/l.period)
!**** CDETR  R*8  NTONUT     detritus concentration for the nutrients
!****                        including carbon (g/m3)
!**** CTODRY R*8  NUSPEC I   dry weigth to carbon ratio of the algae species
!**** DETRIT R*8  NUNUCO I   BLOOM detritus concentration for the
!****                        nutrients excluding carbon (mg/m3)
!**** DETFLX R*8  NTONUT     Flux of nutrients from detritus
!**** DISFLX R*8  NTONUT Total dissolved nutrients flux (mg/l.period)
!**** GDETPR R*8  NUGRAZ I   Preference for detritus
!**** GDETFF R*8  NUGRAZ I   Feacal fraction for detritus
!**** GDETO  R*8  NUNUCO     First order decay rate of detritus due to grazing
!****                        for the previous iteration step
!**** GFECFR R*8  NUSPEC*I   Feacal fraction for algae types
!****             NUGRAZ
!**** GMORO  R*8  NUSPEC     First order decay rate of algae due to grazing
!****                        for the previous iteration step
!**** GRADET R*8  NUNUCO O   First order decay rate of detritus due to grazing
!**** GRAMOR R*8  NUSPEC O   First order decay rate of algae due to grazing
!**** GRAMX  R*8  NUSPEC I   Maximum sustainable grazing rate
!**** GRZBRU R*8         Bruto growth rate grazer (mg C/l.d)
!**** GRZFIL R*8             Realized filtration rate
!**** GRZFOO R*8         Available food (mg C/l)
!**** GRZFM  R*8  NUGRAZ I   Maximum clearance rate grazer (l/mg C.d)
!**** GRZGM  R*8  NUGRAZ I   Maximum relative growth grazer (1/d)
!**** GRZGRZ R*8         Grazing rate (1/d)
!**** GRZMET R*8         Respiration grazer (1/d)
!**** GRZMM  R*8  NUGRAZ I   Maximum relative mortality grazer (1/d)
!**** GRZMO  R*8  NUGRAZ I   Monod term filtr.r. in rel. to foodconc. (mg C/l)
!**** GRZNEW R*8  NUGRAZ     biomass of the grazers provided for the current
!****                        timestep (gC/m3) - can be adapted
!**** GRZOLD R*8  NUGRAZ     biomass of the grazers during the previous
!****                        timestep (gC/m3)
!**** GRZPMX R*8         Max. sustainable bruto growth rate (mg C/l.d)
!**** GRZRM  R*8  NUGRAZ I   Maximum daily ration grazer (mg C/mg C.d)
!**** GRZRAT R*8             Realized feeding rate
!**** GRZRE  R*8  NUGRAZ     Routine respiration coefficient (fraction)
!**** GRZSE  R*8  NUGRAZ I   Standard respiration coefficient (1/d)
!**** GRZST  R*8  NTONUT*    Grazer stoichiometry matrix including carbon
!****             NUGRAZ
!**** GRZUPT R*8  NTONUT Total intake per element (mg/l.d)
!**** GTMPFM R*8  NUGRAZ I   Temperature coefficients filtration rate
!**** GTMPGM R*8  NUGRAZ I   Temperature coefficients maximum growth rate
!**** GTMPMM R*8  NUGRAZ I   Temperature coefficients maximum mortality rate
!**** GTMPRE R*8  NUGRAZ I   Temperature coefficients routine metabolism
!**** GTMPRM R*8  NUGRAZ I   Temperature coefficients feeding rate
!**** GTMPSE R*8  NUGRAZ I   Temperature coefficients standard metabolism
!**** GTODET R*8  NUGRAZ I   Fraction of excretion of grazers to the water colum
!**** IFILSP I*4             Counter in the loop over the grazer types
!**** IPERM  I*4         I   Maximum number of steps in grazing iteration loop
!**** ITNUM  I*4         I   Step number in grazing iteration loop
!**** LCOUPL I*4         I   Switch for stand alone (0) or coupled (1)
!**** LGRAZ  I*4         IO  Switch for start of next iteration step (1)
!**** LZERO  I*4             Switch for zero biomass of food
!**** NTONUT I*4             Number of nutrients for CONSBL, NUNUCO + 1
!**** NUGRAZ I*4         I   Number of grazers
!**** NUNUCO I*4         I   Number of nutrients for BLOOM, excluding carbon
!**** NUSPEC I*4         I   Total number of algae species
!**** PERIOD R*8         I   Duration of period (d)
!**** T      R*8         I   Temperature (degrees Celsius)
!**** TOTGRA R*8             Total grazing rate (1/d)
!**** X      R*8  MX+1   I   A.o. algae biomass (mgC/m3)
!**** ZOOC   R*8  NUGRAZ I   biomass of the grazers during the previous
!****                        timestep (gC/m3)
!**** ZOOD   R*8  NUGRAZ I   biomass of the grazers provided for the current
!****                        timestep (gC/m3)
!**** ZOONUT R*8  NUNUCO*I   BLOOM grazer stoichiometry matrix excluding
!****             NUGRAZ     carbon
!**** ZOOPR  R*8  NUSPEC*I   Preference for algae types
!****             NUGRAZ
      IMPLICIT REAL*8 (A-H,O-Z)
      INCLUDE 'blmdim.inc'
      INCLUDE 'size.inc'
      INCLUDE 'phyt2.inc'
      INCLUDE 'phyt1.inc'
      INCLUDE 'graas.inc'
      INCLUDE 'ioblck.inc'
      DIMENSION X(MX),GRAMOR(MT),ZOOD(0:MG),ZOOC(MG),GRZOLD(MG),
     1          GRZNEW(MG),DISFLX(MN+1),DETFLX(MN+1),BOTFLX(MN+1),
     2          ALGFLX(MT),GRZUPT(MN+1),ALGSTC(MN+1,MT),GRZST(MN+1,MG),
     3          DFLUX(MN+1),ALTFLX(MN+1),GRZFLX(MN+1),TOTFLX(MN+1),
     4          CDETR(MN+1),GRADET(MN),GMORO(MT),GDETO(MN)
      SAVE GMORO,GDETO,LZERO

!****
!*  1 Make a copy of the algal and grazers stochiometry
!****
      ZMAX = 1.D7
      IF (ITNUM.LE.1) THEN
        LZERO = 0
      ELSEIF (LZERO.EQ.1) THEN
        LZERO = 0
        LGRAZ = 0
        RETURN
      ENDIF

      NTONUT = NUNUCO + 1
      DO 280 I=1,NUNUCO
        GRADET(I) = 0.0D0
  280 CONTINUE

      IF ((ITNUM.EQ.1).AND.(ITNUM.LT.IPERM)) THEN
        LGRAZ = 1
      ELSE
        LGRAZ = 0
      ENDIF
      DO 210 J = 1, NUSPEC
          GRAMOR(J) = 0.0D0
          ALGSTC(1,J)=1.0
          DO 220 I = 2, NTONUT
              ALGSTC(I,J) = AA(I-1, J) * CTODRY ( J)
  220     CONTINUE
  210 CONTINUE

      DO 230 J = 1, NUGRAZ
          GRZST(1,J)=1.0
          DO 240 I = 2, NTONUT
              GRZST(I,J) = ZOONUT(I-1, J)
  240     CONTINUE
  230 CONTINUE

!     The amount of carbon in detritus is not known in BLOOM
!     Here we estimate the amount from the amount of nitrogen
!     in detritus. The ratio is taken is the average ratio in the
!     current algae species composition

!     Calcuate the average C/N ratio in the algae
      K = NUROWS
      ALGC = 0.0D0
      ALGN = 0.0D0
      DO 260 J=1,NUSPEC
        K = K + 1
        ALGC = ALGC + X(K)
        ALGN = ALGN + X(K) * ALGSTC(2,J)
  260 CONTINUE

!     Convert DET-N to DET-C and fill the rest of the CDETR array
      IF (ALGN.GT.1.0D-6) THEN
        CDETR(1) = 1.D-3 * DETRIT(1) * ALGC / ALGN
      ELSEIF (CDETR(2).GT.1.0D-6) THEN
        CDETR(1) = 1.D-3 * DETRIT(1) * CDETR(1) / CDETR(2)
      ELSEIF (DETRIT(1).LT.1.0D-6) THEN
        CDETR(1) = 0.0
      ELSE
        CDETR(1) = 10.0D-3 * DETRIT(1)
        WRITE(*,*) 'CONSBL: Could not establish C/N ratio detritus!',
     1             ' Used 10.'
      ENDIF
      DO 270 I=2,NTONUT
        CDETR(I) = 1.D-3 * DETRIT(I-1)
  270 CONTINUE

!
!     The routine does not produce if biomass in GRZNE1 (segfun)
!     evers becomes 0.0, so safeguard added 940524 AH
!
      DO 10 J=1,NUGRAZ
        IF (ZOOC(J).LT.1.0D-20) THEN
          GRZOLD(J) = 1.0D-20
        ELSE
          GRZOLD(J) = ZOOC(J)
        ENDIF
        IF (ZOOD(J).LT.1.0D-20) THEN
          GRZNEW(J) = 1.0D-20
        ELSE
          GRZNEW(J) = ZOOD(J)
        ENDIF
   10 CONTINUE

!     IF ( DEBUG ) THEN
!        WRITE ( IDBG, 9000) ( I, GRZOLD(I), GRZNEW(I), I=1,NUGRAZ)
!     ENDIF

!****
!*  2 Loop over the grazers
!****
      DO 200 IFILSP = 1, NUGRAZ

!****
!*  3 Initialize output variables
!****

          DO 40 I = 1, NTONUT
              DISFLX(I) = 0.0D0
              DETFLX(I) = 0.0D0
              BOTFLX(I) = 0.0D0
   40     CONTINUE
          DO 50 I = 1, NUSPEC
              ALGFLX(I) = 0.0D0
   50     CONTINUE

!****
!*  4 Check if grazer growth or mortality exceeds constraints
!****

          IF ((GRZOLD(IFILSP) .GT. 0.0).AND.(LCOUPL.NE.0)) THEN
              IF ((GRZNEW(IFILSP) - GRZOLD(IFILSP)) .GE. 0.0) THEN
                  IF (GRZNEW(IFILSP) .GT. GRZOLD(IFILSP) *
     &            (1.0 + GRZGM(IFILSP) *
     &            EXP(GTMPGM(IFILSP) *
     &            (T - 20.0)) * PERIOD)) THEN
                      GRZNEW(IFILSP) = GRZOLD(IFILSP) *
     &                (1.0 + GRZGM(IFILSP) *
     &                EXP(GTMPGM(IFILSP) *
     &                (T - 20.0)) * PERIOD)
                  ENDIF
              ELSE
                  IF (GRZNEW(IFILSP) .LT. GRZOLD(IFILSP) *
     &            (1.0 - GRZMM(IFILSP) *
     &            EXP(GTMPMM(IFILSP) *
     &            (T - 20.0)) * PERIOD)) THEN
                      GRZNEW(IFILSP) = GRZOLD(IFILSP) *
     &                (1.0 - GRZMM(IFILSP) *
     &                EXP(GTMPMM(IFILSP) *
     &                (T - 20.0)) * PERIOD)
                  ENDIF
              ENDIF
          ENDIF

!****
!*  5 Calculate total available amount of food (mg C/l)
!****
          GRZFOO = CDETR(1) * GDETPR(IFILSP)
          K = NUROWS
          DO 60 I = 1, NUSPEC
              K = K + 1
              GRZFOO = GRZFOO + X(K) * ZOOPR(I,IFILSP) * 1.D-3
   60     CONTINUE

!****
!*  6 Calculate grazing (1/d) rate and optionally maximum grazing rate
!****
          IF (GRZFOO .GT. 0.0) THEN
              IF (GRZFOO .LT. ((EXP(GTMPRM(IFILSP) *
     &        (T - 20.0)) * GRZRM(IFILSP)) /
     &        (EXP(GTMPFM(IFILSP) * (T - 20.0)) *
     &        GRZFM(IFILSP) *
     &        (GRZFOO / (GRZFOO + GRZMO(IFILSP)))))) THEN
                  GRZGRZ = ((GRZOLD(IFILSP) *
     &            EXP(GTMPFM(IFILSP) * (T - 20.0)) *
     &            GRZFM(IFILSP)) *
     &            (GRZFOO / (GRZFOO + GRZMO(IFILSP))))
              ELSE
                  GRZGRZ = ((GRZOLD(IFILSP) *
     &            EXP(GTMPRM(IFILSP) * (T - 20.0)) *
     &            GRZRM(IFILSP)) / GRZFOO)
              ENDIF
          ELSE
              GRZGRZ = 0.0
              LZERO = 2
          ENDIF
!         IF (GRZGRZ.GT.GRAMX) THEN
!           GRZGRZ=GRAMX
!           LZERO = 1
!         ENDIF

!****
!*  7 Calculate realized feeding (mg C/mg C.d) and filtration rate (l/mg C.d)
!*    GRZRAT and GRZFIL are output variables, not used at the moment
!****
          IF (GRZOLD(IFILSP) .GT. 0.0) THEN
              GRZRAT = GRZFOO * GRZGRZ / GRZOLD(IFILSP)
              GRZFIL = GRZGRZ / GRZOLD(IFILSP)
          ELSE
              GRZRAT = 0.0
              GRZFIL = 0.0
          ENDIF

!****
!*  8 Calculate flux and biomass for detritus and phytoplankton groups
!****
          DO 70 I = 1, NTONUT
              DFLUX(I) = -(CDETR(I) * GRZGRZ * GDETPR(IFILSP))
   70     CONTINUE
          K = NUROWS
          DO 80 I = 1, NUSPEC
              K = K + 1
              ALGFLX(I) = -(X(K) * GRZGRZ * ZOOPR(I,IFILSP)) * 1.D-3
              GRAMOR(I) = GRAMOR(I) + GRZGRZ * ZOOPR(I,IFILSP)
   80     CONTINUE

!****
!*  9 Calculate total intake per element (mg/d)
!****
          DO 100 J = 1, NTONUT
              GRZUPT(J) = -(DFLUX(J) * (1.0 - GDETFF(IFILSP)))
              DETFLX(J) = DFLUX(J) - DFLUX(J) *
     &            GDETFF(IFILSP) * GTODET(IFILSP)
              BOTFLX(J) = BOTFLX(J) - DFLUX(J) *
     &            GDETFF(IFILSP) * (1.0-GTODET(IFILSP))
              DO 90 I = 1, NUSPEC
                  GRZUPT(J) = GRZUPT(J) - ALGFLX(I) * ALGSTC(J,I) *
     1            (1.0 - GFECFR(I,IFILSP))
                  DETFLX(J) = DETFLX(J) - ALGFLX(I) * ALGSTC(J,I) *
     1                GFECFR(I,IFILSP) * GTODET(IFILSP)
                  BOTFLX(J) = BOTFLX(J) - ALGFLX(I) * ALGSTC(J,I) *
     1                GFECFR(I,IFILSP) * (1.0 - GTODET(IFILSP))
   90         CONTINUE
  100     CONTINUE

!****
!* 10 Calculate limiting element for growth
!****
          GRZPMX = 10.0E20
          DO 110 I = 1, NTONUT
              IF (GRZST(I,IFILSP) .GT. 0.0) THEN
                  GRZPRT = GRZUPT(I) /
     &            GRZST(I,IFILSP)
                  IF (GRZPRT .LT. GRZPMX) THEN
                      GRZPMX = GRZPRT
                  ENDIF
              ENDIF
  110     CONTINUE

!****
!* 11 Calculate routine respiration (mgC/l.d)
!****
          DO 120 I = 1, NTONUT
              DISFLX(I) = (GRZPMX * EXP(GTMPRE(IFILSP) *
     &        (T - 20.0)) * GRZRE(IFILSP) *
     &        GRZST(I,IFILSP))
              GRZUPT(I) = GRZUPT(I) - DISFLX(I)
  120     CONTINUE
          GRZPMX = GRZPMX * (1.0 - EXP(GTMPRE(IFILSP) *
     &    (T - 20.0)) * GRZRE(IFILSP))

!****
!* 12 Calculate the standard respiration (1/d)
!****
          GRZMET = EXP(GTMPSE(IFILSP) * (T - 20.0)) *
     &             GRZSE(IFILSP)

!****
!* 13 Correct for length of period (d)
!****
          DO 130 I = 1, NUSPEC
              ALGFLX(I) = ALGFLX(I) * PERIOD
  130     CONTINUE
          DO 140 I = 1, NTONUT
              GRZUPT(I) = GRZUPT(I) * PERIOD
              DETFLX(I) = DETFLX(I) * PERIOD
              BOTFLX(I) = BOTFLX(I) * PERIOD
              DISFLX(I) = DISFLX(I) * PERIOD
  140     CONTINUE
          GRZPMX = GRZPMX * PERIOD
          GRZMET = GRZMET * PERIOD

!****
!* 14 Calculate bruto growth (mg C/period)
!****
          GRZBRU = GRZNEW(IFILSP) - GRZOLD(IFILSP) * (1.0 - GRZMET)

!****
!* 15 Correct bruto growth if intake can not sustain respiration and growth
!****
          IF ((GRZBRU .GT. GRZPMX).AND.(LCOUPL.NE.0)) THEN
              GRZBRU = GRZPMX
              GRZNEW(IFILSP) = GRZOLD(IFILSP) * (1.0 - GRZMET) + GRZPMX
          ENDIF

!****
!* 16 Add respiration to dissolved nutrients
!****
          DO 150 I = 1, NTONUT
              DISFLX(I) = DISFLX(I) + GRZOLD(IFILSP) * GRZMET *
     &        GRZST(I,IFILSP)

!****
!* 17 If there is bruto growth, subtract nutrients from the intake
!****
              IF (GRZBRU .GE. 0.0) THEN
                  GRZUPT(I) = GRZUPT(I) - GRZBRU *
     &            GRZST(I,IFILSP)
!****
!* 18 If their is mortality, add the nutrients to the detritus pool
!****
              ELSE
                  DETFLX(I) = DETFLX(I) - GRZBRU *
     &                GRZST(I,IFILSP) * GTODET(IFILSP)
                  BOTFLX(I) = BOTFLX(I) - GRZBRU *
     &                GRZST(I,IFILSP) * (1.0 - GTODET(IFILSP))
              ENDIF

!****
!* 19 Convert the intake not used for bruto growth to the detritus pool
!****
              DETFLX(I) = DETFLX(I) + GRZUPT(I) * GTODET(IFILSP)
              BOTFLX(I) = BOTFLX(I) + GRZUPT(I) * (1.0 - GTODET(IFILSP))
  150     CONTINUE

!****
!* 20 Check massbalans
!****
          DO 170 J=1,NTONUT
              ALTFLX(J) = 0.0E0
              DO 160 I=1,NUSPEC
                  ALTFLX(J) = ALTFLX(J) + ALGFLX(I) * ALGSTC(J,I)
  160         CONTINUE
              GRZFLX(J) = (GRZNEW(IFILSP) - GRZOLD(IFILSP)) *
     &        GRZST(J,IFILSP)
              TOTFLX(J) = DISFLX(J) + DETFLX(J) + BOTFLX(J) +
     &        ALTFLX(J) + GRZFLX(J)
!             IF ( DEBUG ) WRITE ( IDBG), 9030) J, TOTFLX(J)
!
  170     CONTINUE

!****
!* 21 Update total fluxes detritus
!****
          DO 290 I=1,NUNUCO
            IF (CDETR(1).GT.1.D-12) THEN
!             GRADET will be used for calculation of the change in
!             detritus concentration according to
!             DETR = DETR * EXP(-GRADET*DT)
!             uit de flux  DD=(DETFLX + BOTFLX) kan GRADET afgeleid worden
!             GRADET = - LOG(1+DD/DETR) / DT
              IF ((DETFLX(I+1)+BOTFLX(I+1))/CDETR(I+1).LE.-1.D0) THEN
                WRITE(*,*) 'CONSBL: Grazing of detritus becomes ',
     1                     'unstable!'
                WRITE(*,*) '        This is caused by an unrealistic ',
     1                     'high grazing pressure.'
                WRITE(*,*) '        Please adapt the grazing pressure ',
     1                     'in the input.'
                STOP 1
              ENDIF
              GRADET(I) = GRADET(I) - DLOG(1.0D0 + (DETFLX(I+1)+
     1                    BOTFLX(I+1))/(CDETR(I+1))) / PERIOD
!           ELSEIF (ABS(DETFLX(1)+BOTFLX(1)).GT.1.D-12) THEN
!             WRITE(*,*) 'CONSBL: No detritus to add feacal fraction ',
!    1                   'grazers to!'
!             WRITE(*,*) '        Nutients will become available ',
!    1                   'immediately as dissolved.'
            ENDIF
  290     CONTINUE
!         WRITE(*,*) 'CDETR  = ',CDETR(1)
!         WRITE(*,*) 'DFLUX  = ',DFLUX(1)
!         WRITE(*,*) 'DETFLX = ',DETFLX(1)
!         WRITE(*,*) 'BOTFLX = ',BOTFLX(1)
!         WRITE(*,*) 'ALGFLX = ',ALGFLX(1)

!****
!* 22 End loop over filter feeders
!****
  200 CONTINUE

!****
!* 23 Save biomass grazers for next time step
!****
      TOTGRA = 0.0
      DO 250 J=1,NUGRAZ
          TOTGRA=TOTGRA+GRAMOR(J)
  250 CONTINUE

      IF ((ITNUM.GT.1).AND.(ITNUM.LT.IPERM).AND.(LGRAZ.NE.2)) THEN
        DO 300 I=1,NUNUCO
          IF (ABS(GRADET(I)-GDETO(I)).GT.1.0D-06) LGRAZ=1
  300   CONTINUE
        IF (LGRAZ.EQ.0) THEN
          DO 310 J=1,NUSPEC
            IF (ABS(GRAMOR(J)-GMORO(J)).GT.1.0D-06) LGRAZ=1
  310     CONTINUE
        ENDIF
      ENDIF

      IF ((ITNUM.LT.IPERM).AND.(LGRAZ.EQ.1)) THEN
        DO 320 I=1,NUNUCO
          GDETO(I)=GRADET(I)
  320   CONTINUE
        DO 330 J=1,NUSPEC
          GMORO(J)=GRAMOR(J)
  330   CONTINUE
      ENDIF

      IF (LZERO.EQ.2) THEN
        LZERO = 0
        LGRAZ = 0
      ENDIF

!     WRITE(244,'(I2,5(1PE12.4))') ITNUM,GRZOLD(1),GRZNEW(1),GRZFOO,
!    1                             TOTGRA,GRAMX

      RETURN

 9000 FORMAT ( ' CONSBL: ' / ( ' Biomass grazer ', I2, ' old ', E12.5,
     1         '   wanted ', E12.5))
 9010 FORMAT ( '   FLUX nut ',I2, ' dis, det, bot   ', 3E12.5)
 9020 FORMAT ( '   FLUX alg ',I2, '    ', E12.5)
 9030 FORMAT ( '   TOTFLX   ',I2, '    ', E12.5)
 9040 FORMAT ( '   TOT alg  ',I2, '    ', E12.5)

      END

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
	SUBROUTINE GSMODE  ( mode, iret )
C************************************************************************
C* GSMODE								*
C* 									*
C* This subroutine sets the ploting mode for map/graph coordinate	*
C* plotting.  It may be used to change mode within programs.		*
C*									*
C* GSMODE  ( MODE, IRET )						*
C*									*
C* Input parameters:							*
C*	MODE		INTEGER		Plotting mode 			*
C*					   0 = no change		*
C* 					   1 = map coordinates		*
C* 					   2 = graph coordinates	*
C* 									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* G. Chatters/RDS	 8/85						*
C* M. desJardins/GSFC	 4/85	GEMPLT 3.1				*
C* M. desJardins/NMC	07/91	Check for existence of device		*
C* K. Brill/NMC		08/91	Added DEVCHR.CMN to check ddev		*
C* S. Jacobs/NMC	 1/95	Added multi-window common block		*
C* M. Linda/GSC		 2/97	Removed GFLUSH				*
C* M. Li/GSC		 6/00	Removed erroneous setting of nmode	*
C************************************************************************
	INCLUDE		'XYDEF.CMN'
	INCLUDE		'ERROR.PRM'
C------------------------------------------------------------------------
	iret = NORMAL
C
C*	Set mode flag in common.  Update common areas if there is a
C*	change.
C
	IF  ( ( mode .eq. 1 ) .or. ( mode .eq. 2 ) )  THEN
	    IF  ( igmode .ne. mode )  THEN
		igmode = mode
		CALL UPDPXY
	    END IF
C
C*	    Check for invalid mode.
C
	  ELSE IF  ( mode .ne. 0 )  THEN
	    iret = NIMODE
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text: #redirect Bicycle Parking

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search