!*==dtrsyl.f90  processed by SPAG 7.51RB at 20:08 on  3 Mar 2022
!> \brief \b DTRSYL
!
!  =========== DOCUMENTATION ===========
!
! Online html documentation available at
!            http://www.netlib.org/lapack/explore-html/
!
!> \htmlonly
!> Download DTRSYL + dependencies
!>
!> [TGZ]
!>
!> [ZIP]
!>
!> [TXT]
!> \endhtmlonly
!
!  Definition:
!  ===========
!
!       SUBROUTINE DTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C,
!                          LDC, SCALE, INFO )
!
!       .. Scalar Arguments ..
!       CHARACTER          TRANA, TRANB
!       INTEGER            INFO, ISGN, LDA, LDB, LDC, M, N
!       DOUBLE PRECISION   SCALE
!       ..
!       .. Array Arguments ..
!       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), C( LDC, * )
!       ..
!
!
!> \par Purpose:
!  =============
!>
!> \verbatim
!>
!> DTRSYL solves the real Sylvester matrix equation:
!>
!>    op(A)*X + X*op(B) = scale*C or
!>    op(A)*X - X*op(B) = scale*C,
!>
!> where op(A) = A or A**T, and  A and B are both upper quasi-
!> triangular. A is M-by-M and B is N-by-N; the right hand side C and
!> the solution X are M-by-N; and scale is an output scale factor, set
!> <= 1 to avoid overflow in X.
!>
!> A and B must be in Schur canonical form (as returned by DHSEQR), that
!> is, block upper triangular with 1-by-1 and 2-by-2 diagonal blocks;
!> each 2-by-2 diagonal block has its diagonal elements equal and its
!> off-diagonal elements of opposite sign.
!> \endverbatim
!
!  Arguments:
!  ==========
!
!> \param[in] TRANA
!> \verbatim
!>          TRANA is CHARACTER*1
!>          Specifies the option op(A):
!>          = 'N': op(A) = A    (No transpose)
!>          = 'T': op(A) = A**T (Transpose)
!>          = 'C': op(A) = A**H (Conjugate transpose = Transpose)
!> \endverbatim
!>
!> \param[in] TRANB
!> \verbatim
!>          TRANB is CHARACTER*1
!>          Specifies the option op(B):
!>          = 'N': op(B) = B    (No transpose)
!>          = 'T': op(B) = B**T (Transpose)
!>          = 'C': op(B) = B**H (Conjugate transpose = Transpose)
!> \endverbatim
!>
!> \param[in] ISGN
!> \verbatim
!>          ISGN is INTEGER
!>          Specifies the sign in the equation:
!>          = +1: solve op(A)*X + X*op(B) = scale*C
!>          = -1: solve op(A)*X - X*op(B) = scale*C
!> \endverbatim
!>
!> \param[in] M
!> \verbatim
!>          M is INTEGER
!>          The order of the matrix A, and the number of rows in the
!>          matrices X and C. M >= 0.
!> \endverbatim
!>
!> \param[in] N
!> \verbatim
!>          N is INTEGER
!>          The order of the matrix B, and the number of columns in the
!>          matrices X and C. N >= 0.
!> \endverbatim
!>
!> \param[in] A
!> \verbatim
!>          A is DOUBLE PRECISION array, dimension (LDA,M)
!>          The upper quasi-triangular matrix A, in Schur canonical form.
!> \endverbatim
!>
!> \param[in] LDA
!> \verbatim
!>          LDA is INTEGER
!>          The leading dimension of the array A. LDA >= max(1,M).
!> \endverbatim
!>
!> \param[in] B
!> \verbatim
!>          B is DOUBLE PRECISION array, dimension (LDB,N)
!>          The upper quasi-triangular matrix B, in Schur canonical form.
!> \endverbatim
!>
!> \param[in] LDB
!> \verbatim
!>          LDB is INTEGER
!>          The leading dimension of the array B. LDB >= max(1,N).
!> \endverbatim
!>
!> \param[in,out] C
!> \verbatim
!>          C is DOUBLE PRECISION array, dimension (LDC,N)
!>          On entry, the M-by-N right hand side matrix C.
!>          On exit, C is overwritten by the solution matrix X.
!> \endverbatim
!>
!> \param[in] LDC
!> \verbatim
!>          LDC is INTEGER
!>          The leading dimension of the array C. LDC >= max(1,M)
!> \endverbatim
!>
!> \param[out] SCALE
!> \verbatim
!>          SCALE is DOUBLE PRECISION
!>          The scale factor, scale, set <= 1 to avoid overflow in X.
!> \endverbatim
!>
!> \param[out] INFO
!> \verbatim
!>          INFO is INTEGER
!>          = 0: successful exit
!>          < 0: if INFO = -i, the i-th argument had an illegal value
!>          = 1: A and B have common or very close eigenvalues; perturbed
!>               values were used to solve the equation (but the matrices
!>               A and B are unchanged).
!> \endverbatim
!
!  Authors:
!  ========
!
!> \author Univ. of Tennessee
!> \author Univ. of California Berkeley
!> \author Univ. of Colorado Denver
!> \author NAG Ltd.
!
!> \date December 2016
!
!> \ingroup doubleSYcomputational
!
!  =====================================================================
      SUBROUTINE DTRSYL(Trana,Tranb,Isgn,M,N,A,Lda,B,Ldb,C,Ldc,Scale,   &
     &                  Info)
      IMPLICIT NONE
!*--DTRSYL168
!
!  -- LAPACK computational routine (version 3.7.0) --
!  -- LAPACK is a software package provided by Univ. of Tennessee,    --
!  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
!     December 2016
!
!     .. Scalar Arguments ..
      CHARACTER Trana , Tranb
      INTEGER Info , Isgn , Lda , Ldb , Ldc , M , N
      DOUBLE PRECISION Scale
!     ..
!     .. Array Arguments ..
      DOUBLE PRECISION A(Lda,*) , B(Ldb,*) , C(Ldc,*)
!     ..
!
!  =====================================================================
!
!     .. Parameters ..
      DOUBLE PRECISION ZERO , ONE
      PARAMETER (ZERO=0.0D+0,ONE=1.0D+0)
!     ..
!     .. Local Scalars ..
      LOGICAL notrna , notrnb
      INTEGER ierr , j , k , k1 , k2 , knext , l , l1 , l2 , lnext
      DOUBLE PRECISION a11 , bignum , da11 , db , eps , scaloc , sgn ,  &
     &                 smin , smlnum , suml , sumr , xnorm
!     ..
!     .. Local Arrays ..
      DOUBLE PRECISION dum(1) , vec(2,2) , x(2,2)
!     ..
!     .. External Functions ..
      LOGICAL LSAME
      DOUBLE PRECISION DDOT , DLAMCH , DLANGE
      EXTERNAL LSAME , DDOT , DLAMCH , DLANGE
!     ..
!     .. External Subroutines ..
      EXTERNAL DLABAD , DLALN2 , DLASY2 , DSCAL , XERBLA
!     ..
!     .. Intrinsic Functions ..
      INTRINSIC ABS , DBLE , MAX , MIN
!     ..
!     .. Executable Statements ..
!
!     Decode and Test input parameters
!
      notrna = LSAME(Trana,'N')
      notrnb = LSAME(Tranb,'N')
!
      Info = 0
      IF ( .NOT.notrna .AND. .NOT.LSAME(Trana,'T') .AND.                &
     &     .NOT.LSAME(Trana,'C') ) THEN
         Info = -1
      ELSEIF ( .NOT.notrnb .AND. .NOT.LSAME(Tranb,'T') .AND.            &
     &         .NOT.LSAME(Tranb,'C') ) THEN
         Info = -2
      ELSEIF ( Isgn/=1 .AND. Isgn/=-1 ) THEN
         Info = -3
      ELSEIF ( M<0 ) THEN
         Info = -4
      ELSEIF ( N<0 ) THEN
         Info = -5
      ELSEIF ( Lda<MAX(1,M) ) THEN
         Info = -7
      ELSEIF ( Ldb<MAX(1,N) ) THEN
         Info = -9
      ELSEIF ( Ldc<MAX(1,M) ) THEN
         Info = -11
      ENDIF
      IF ( Info/=0 ) THEN
         CALL XERBLA('DTRSYL',-Info)
         RETURN
      ENDIF
!
!     Quick return if possible
!
      Scale = ONE
      IF ( M==0 .OR. N==0 ) RETURN
!
!     Set constants to control overflow
!
      eps = DLAMCH('P')
      smlnum = DLAMCH('S')
      bignum = ONE/smlnum
      CALL DLABAD(smlnum,bignum)
      smlnum = smlnum*DBLE(M*N)/eps
      bignum = ONE/smlnum
!
      smin = MAX(smlnum,eps*DLANGE('M',M,M,A,Lda,dum),                  &
     &       eps*DLANGE('M',N,N,B,Ldb,dum))
!
      sgn = Isgn
!
      IF ( notrna .AND. notrnb ) THEN
!
!        Solve    A*X + ISGN*X*B = scale*C.
!
!        The (K,L)th block of X is determined starting from
!        bottom-left corner column by column by
!
!         A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
!
!        Where
!                  M                         L-1
!        R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)].
!                I=K+1                       J=1
!
!        Start column loop (index = L)
!        L1 (L2) : column index of the first (first) row of X(K,L).
!
         lnext = 1
         DO l = 1 , N
            IF ( l>=lnext ) THEN
               IF ( l==N ) THEN
                  l1 = l
                  l2 = l
               ELSEIF ( B(l+1,l)/=ZERO ) THEN
                  l1 = l
                  l2 = l + 1
                  lnext = l + 2
               ELSE
                  l1 = l
                  l2 = l
                  lnext = l + 1
               ENDIF
!
!           Start row loop (index = K)
!           K1 (K2): row index of the first (last) row of X(K,L).
!
               knext = M
               DO k = M , 1 , -1
                  IF ( k<=knext ) THEN
                     IF ( k==1 ) THEN
                        k1 = k
                        k2 = k
                     ELSEIF ( A(k,k-1)/=ZERO ) THEN
                        k1 = k - 1
                        k2 = k
                        knext = k - 2
                     ELSE
                        k1 = k
                        k2 = k
                        knext = k - 1
                     ENDIF
!
                     IF ( l1==l2 .AND. k1==k2 ) THEN
                        suml = DDOT(M-k1,A(k1,MIN(k1+1,M)),Lda,         &
     &                         C(MIN(k1+1,M),l1),1)
                        sumr = DDOT(l1-1,C(k1,1),Ldc,B(1,l1),1)
                        vec(1,1) = C(k1,l1) - (suml+sgn*sumr)
                        scaloc = ONE
!
                        a11 = A(k1,k1) + sgn*B(l1,l1)
                        da11 = ABS(a11)
                        IF ( da11<=smin ) THEN
                           a11 = smin
                           da11 = smin
                           Info = 1
                        ENDIF
                        db = ABS(vec(1,1))
                        IF ( da11<ONE .AND. db>ONE ) THEN
                           IF ( db>bignum*da11 ) scaloc = ONE/db
                        ENDIF
                        x(1,1) = (vec(1,1)*scaloc)/a11
!
                        IF ( scaloc/=ONE ) THEN
                           DO j = 1 , N
                              CALL DSCAL(M,scaloc,C(1,j),1)
                           ENDDO
                           Scale = Scale*scaloc
                        ENDIF
                        C(k1,l1) = x(1,1)
!
                     ELSEIF ( l1==l2 .AND. k1/=k2 ) THEN
!
                        suml = DDOT(M-k2,A(k1,MIN(k2+1,M)),Lda,         &
     &                         C(MIN(k2+1,M),l1),1)
                        sumr = DDOT(l1-1,C(k1,1),Ldc,B(1,l1),1)
                        vec(1,1) = C(k1,l1) - (suml+sgn*sumr)
!
                        suml = DDOT(M-k2,A(k2,MIN(k2+1,M)),Lda,         &
     &                         C(MIN(k2+1,M),l1),1)
                        sumr = DDOT(l1-1,C(k2,1),Ldc,B(1,l1),1)
                        vec(2,1) = C(k2,l1) - (suml+sgn*sumr)
!
                        CALL DLALN2(.FALSE.,2,1,smin,ONE,A(k1,k1),Lda,  &
     &                              ONE,ONE,vec,2,-sgn*B(l1,l1),ZERO,x, &
     &                              2,scaloc,xnorm,ierr)
                        IF ( ierr/=0 ) Info = 1
!
                        IF ( scaloc/=ONE ) THEN
                           DO j = 1 , N
                              CALL DSCAL(M,scaloc,C(1,j),1)
                           ENDDO
                           Scale = Scale*scaloc
                        ENDIF
                        C(k1,l1) = x(1,1)
                        C(k2,l1) = x(2,1)
!
                     ELSEIF ( l1/=l2 .AND. k1==k2 ) THEN
!
                        suml = DDOT(M-k1,A(k1,MIN(k1+1,M)),Lda,         &
     &                         C(MIN(k1+1,M),l1),1)
                        sumr = DDOT(l1-1,C(k1,1),Ldc,B(1,l1),1)
                        vec(1,1) = sgn*(C(k1,l1)-(suml+sgn*sumr))
!
                        suml = DDOT(M-k1,A(k1,MIN(k1+1,M)),Lda,         &
     &                         C(MIN(k1+1,M),l2),1)
                        sumr = DDOT(l1-1,C(k1,1),Ldc,B(1,l2),1)
                        vec(2,1) = sgn*(C(k1,l2)-(suml+sgn*sumr))
!
                        CALL DLALN2(.TRUE.,2,1,smin,ONE,B(l1,l1),Ldb,   &
     &                              ONE,ONE,vec,2,-sgn*A(k1,k1),ZERO,x, &
     &                              2,scaloc,xnorm,ierr)
                        IF ( ierr/=0 ) Info = 1
!
                        IF ( scaloc/=ONE ) THEN
                           DO j = 1 , N
                              CALL DSCAL(M,scaloc,C(1,j),1)
                           ENDDO
                           Scale = Scale*scaloc
                        ENDIF
                        C(k1,l1) = x(1,1)
                        C(k1,l2) = x(2,1)
!
                     ELSEIF ( l1/=l2 .AND. k1/=k2 ) THEN
!
                        suml = DDOT(M-k2,A(k1,MIN(k2+1,M)),Lda,         &
     &                         C(MIN(k2+1,M),l1),1)
                        sumr = DDOT(l1-1,C(k1,1),Ldc,B(1,l1),1)
                        vec(1,1) = C(k1,l1) - (suml+sgn*sumr)
!
                        suml = DDOT(M-k2,A(k1,MIN(k2+1,M)),Lda,         &
     &                         C(MIN(k2+1,M),l2),1)
                        sumr = DDOT(l1-1,C(k1,1),Ldc,B(1,l2),1)
                        vec(1,2) = C(k1,l2) - (suml+sgn*sumr)
!
                        suml = DDOT(M-k2,A(k2,MIN(k2+1,M)),Lda,         &
     &                         C(MIN(k2+1,M),l1),1)
                        sumr = DDOT(l1-1,C(k2,1),Ldc,B(1,l1),1)
                        vec(2,1) = C(k2,l1) - (suml+sgn*sumr)
!
                        suml = DDOT(M-k2,A(k2,MIN(k2+1,M)),Lda,         &
     &                         C(MIN(k2+1,M),l2),1)
                        sumr = DDOT(l1-1,C(k2,1),Ldc,B(1,l2),1)
                        vec(2,2) = C(k2,l2) - (suml+sgn*sumr)
!
                        CALL DLASY2(.FALSE.,.FALSE.,Isgn,2,2,A(k1,k1),  &
     &                              Lda,B(l1,l1),Ldb,vec,2,scaloc,x,2,  &
     &                              xnorm,ierr)
                        IF ( ierr/=0 ) Info = 1
!
                        IF ( scaloc/=ONE ) THEN
                           DO j = 1 , N
                              CALL DSCAL(M,scaloc,C(1,j),1)
                           ENDDO
                           Scale = Scale*scaloc
                        ENDIF
                        C(k1,l1) = x(1,1)
                        C(k1,l2) = x(1,2)
                        C(k2,l1) = x(2,1)
                        C(k2,l2) = x(2,2)
                     ENDIF
                  ENDIF
!
               ENDDO
            ENDIF
!
         ENDDO
!
      ELSEIF ( .NOT.notrna .AND. notrnb ) THEN
!
!        Solve    A**T *X + ISGN*X*B = scale*C.
!
!        The (K,L)th block of X is determined starting from
!        upper-left corner column by column by
!
!          A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
!
!        Where
!                   K-1        T                    L-1
!          R(K,L) = SUM [A(I,K)**T*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)]
!                   I=1                          J=1
!
!        Start column loop (index = L)
!        L1 (L2): column index of the first (last) row of X(K,L)
!
         lnext = 1
         DO l = 1 , N
            IF ( l>=lnext ) THEN
               IF ( l==N ) THEN
                  l1 = l
                  l2 = l
               ELSEIF ( B(l+1,l)/=ZERO ) THEN
                  l1 = l
                  l2 = l + 1
                  lnext = l + 2
               ELSE
                  l1 = l
                  l2 = l
                  lnext = l + 1
               ENDIF
!
!           Start row loop (index = K)
!           K1 (K2): row index of the first (last) row of X(K,L)
!
               knext = 1
               DO k = 1 , M
                  IF ( k>=knext ) THEN
                     IF ( k==M ) THEN
                        k1 = k
                        k2 = k
                     ELSEIF ( A(k+1,k)/=ZERO ) THEN
                        k1 = k
                        k2 = k + 1
                        knext = k + 2
                     ELSE
                        k1 = k
                        k2 = k
                        knext = k + 1
                     ENDIF
!
                     IF ( l1==l2 .AND. k1==k2 ) THEN
                        suml = DDOT(k1-1,A(1,k1),1,C(1,l1),1)
                        sumr = DDOT(l1-1,C(k1,1),Ldc,B(1,l1),1)
                        vec(1,1) = C(k1,l1) - (suml+sgn*sumr)
                        scaloc = ONE
!
                        a11 = A(k1,k1) + sgn*B(l1,l1)
                        da11 = ABS(a11)
                        IF ( da11<=smin ) THEN
                           a11 = smin
                           da11 = smin
                           Info = 1
                        ENDIF
                        db = ABS(vec(1,1))
                        IF ( da11<ONE .AND. db>ONE ) THEN
                           IF ( db>bignum*da11 ) scaloc = ONE/db
                        ENDIF
                        x(1,1) = (vec(1,1)*scaloc)/a11
!
                        IF ( scaloc/=ONE ) THEN
                           DO j = 1 , N
                              CALL DSCAL(M,scaloc,C(1,j),1)
                           ENDDO
                           Scale = Scale*scaloc
                        ENDIF
                        C(k1,l1) = x(1,1)
!
                     ELSEIF ( l1==l2 .AND. k1/=k2 ) THEN
!
                        suml = DDOT(k1-1,A(1,k1),1,C(1,l1),1)
                        sumr = DDOT(l1-1,C(k1,1),Ldc,B(1,l1),1)
                        vec(1,1) = C(k1,l1) - (suml+sgn*sumr)
!
                        suml = DDOT(k1-1,A(1,k2),1,C(1,l1),1)
                        sumr = DDOT(l1-1,C(k2,1),Ldc,B(1,l1),1)
                        vec(2,1) = C(k2,l1) - (suml+sgn*sumr)
!
                        CALL DLALN2(.TRUE.,2,1,smin,ONE,A(k1,k1),Lda,   &
     &                              ONE,ONE,vec,2,-sgn*B(l1,l1),ZERO,x, &
     &                              2,scaloc,xnorm,ierr)
                        IF ( ierr/=0 ) Info = 1
!
                        IF ( scaloc/=ONE ) THEN
                           DO j = 1 , N
                              CALL DSCAL(M,scaloc,C(1,j),1)
                           ENDDO
                           Scale = Scale*scaloc
                        ENDIF
                        C(k1,l1) = x(1,1)
                        C(k2,l1) = x(2,1)
!
                     ELSEIF ( l1/=l2 .AND. k1==k2 ) THEN
!
                        suml = DDOT(k1-1,A(1,k1),1,C(1,l1),1)
                        sumr = DDOT(l1-1,C(k1,1),Ldc,B(1,l1),1)
                        vec(1,1) = sgn*(C(k1,l1)-(suml+sgn*sumr))
!
                        suml = DDOT(k1-1,A(1,k1),1,C(1,l2),1)
                        sumr = DDOT(l1-1,C(k1,1),Ldc,B(1,l2),1)
                        vec(2,1) = sgn*(C(k1,l2)-(suml+sgn*sumr))
!
                        CALL DLALN2(.TRUE.,2,1,smin,ONE,B(l1,l1),Ldb,   &
     &                              ONE,ONE,vec,2,-sgn*A(k1,k1),ZERO,x, &
     &                              2,scaloc,xnorm,ierr)
                        IF ( ierr/=0 ) Info = 1
!
                        IF ( scaloc/=ONE ) THEN
                           DO j = 1 , N
                              CALL DSCAL(M,scaloc,C(1,j),1)
                           ENDDO
                           Scale = Scale*scaloc
                        ENDIF
                        C(k1,l1) = x(1,1)
                        C(k1,l2) = x(2,1)
!
                     ELSEIF ( l1/=l2 .AND. k1/=k2 ) THEN
!
                        suml = DDOT(k1-1,A(1,k1),1,C(1,l1),1)
                        sumr = DDOT(l1-1,C(k1,1),Ldc,B(1,l1),1)
                        vec(1,1) = C(k1,l1) - (suml+sgn*sumr)
!
                        suml = DDOT(k1-1,A(1,k1),1,C(1,l2),1)
                        sumr = DDOT(l1-1,C(k1,1),Ldc,B(1,l2),1)
                        vec(1,2) = C(k1,l2) - (suml+sgn*sumr)
!
                        suml = DDOT(k1-1,A(1,k2),1,C(1,l1),1)
                        sumr = DDOT(l1-1,C(k2,1),Ldc,B(1,l1),1)
                        vec(2,1) = C(k2,l1) - (suml+sgn*sumr)
!
                        suml = DDOT(k1-1,A(1,k2),1,C(1,l2),1)
                        sumr = DDOT(l1-1,C(k2,1),Ldc,B(1,l2),1)
                        vec(2,2) = C(k2,l2) - (suml+sgn*sumr)
!
                        CALL DLASY2(.TRUE.,.FALSE.,Isgn,2,2,A(k1,k1),   &
     &                              Lda,B(l1,l1),Ldb,vec,2,scaloc,x,2,  &
     &                              xnorm,ierr)
                        IF ( ierr/=0 ) Info = 1
!
                        IF ( scaloc/=ONE ) THEN
                           DO j = 1 , N
                              CALL DSCAL(M,scaloc,C(1,j),1)
                           ENDDO
                           Scale = Scale*scaloc
                        ENDIF
                        C(k1,l1) = x(1,1)
                        C(k1,l2) = x(1,2)
                        C(k2,l1) = x(2,1)
                        C(k2,l2) = x(2,2)
                     ENDIF
                  ENDIF
!
               ENDDO
            ENDIF
         ENDDO
!
      ELSEIF ( .NOT.notrna .AND. .NOT.notrnb ) THEN
!
!        Solve    A**T*X + ISGN*X*B**T = scale*C.
!
!        The (K,L)th block of X is determined starting from
!        top-right corner column by column by
!
!           A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L)
!
!        Where
!                     K-1                            N
!            R(K,L) = SUM [A(I,K)**T*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T].
!                     I=1                          J=L+1
!
!        Start column loop (index = L)
!        L1 (L2): column index of the first (last) row of X(K,L)
!
         lnext = N
         DO l = N , 1 , -1
            IF ( l<=lnext ) THEN
               IF ( l==1 ) THEN
                  l1 = l
                  l2 = l
               ELSEIF ( B(l,l-1)/=ZERO ) THEN
                  l1 = l - 1
                  l2 = l
                  lnext = l - 2
               ELSE
                  l1 = l
                  l2 = l
                  lnext = l - 1
               ENDIF
!
!           Start row loop (index = K)
!           K1 (K2): row index of the first (last) row of X(K,L)
!
               knext = 1
               DO k = 1 , M
                  IF ( k>=knext ) THEN
                     IF ( k==M ) THEN
                        k1 = k
                        k2 = k
                     ELSEIF ( A(k+1,k)/=ZERO ) THEN
                        k1 = k
                        k2 = k + 1
                        knext = k + 2
                     ELSE
                        k1 = k
                        k2 = k
                        knext = k + 1
                     ENDIF
!
                     IF ( l1==l2 .AND. k1==k2 ) THEN
                        suml = DDOT(k1-1,A(1,k1),1,C(1,l1),1)
                        sumr = DDOT(N-l1,C(k1,MIN(l1+1,N)),Ldc,         &
     &                         B(l1,MIN(l1+1,N)),Ldb)
                        vec(1,1) = C(k1,l1) - (suml+sgn*sumr)
                        scaloc = ONE
!
                        a11 = A(k1,k1) + sgn*B(l1,l1)
                        da11 = ABS(a11)
                        IF ( da11<=smin ) THEN
                           a11 = smin
                           da11 = smin
                           Info = 1
                        ENDIF
                        db = ABS(vec(1,1))
                        IF ( da11<ONE .AND. db>ONE ) THEN
                           IF ( db>bignum*da11 ) scaloc = ONE/db
                        ENDIF
                        x(1,1) = (vec(1,1)*scaloc)/a11
!
                        IF ( scaloc/=ONE ) THEN
                           DO j = 1 , N
                              CALL DSCAL(M,scaloc,C(1,j),1)
                           ENDDO
                           Scale = Scale*scaloc
                        ENDIF
                        C(k1,l1) = x(1,1)
!
                     ELSEIF ( l1==l2 .AND. k1/=k2 ) THEN
!
                        suml = DDOT(k1-1,A(1,k1),1,C(1,l1),1)
                        sumr = DDOT(N-l2,C(k1,MIN(l2+1,N)),Ldc,         &
     &                         B(l1,MIN(l2+1,N)),Ldb)
                        vec(1,1) = C(k1,l1) - (suml+sgn*sumr)
!
                        suml = DDOT(k1-1,A(1,k2),1,C(1,l1),1)
                        sumr = DDOT(N-l2,C(k2,MIN(l2+1,N)),Ldc,         &
     &                         B(l1,MIN(l2+1,N)),Ldb)
                        vec(2,1) = C(k2,l1) - (suml+sgn*sumr)
!
                        CALL DLALN2(.TRUE.,2,1,smin,ONE,A(k1,k1),Lda,   &
     &                              ONE,ONE,vec,2,-sgn*B(l1,l1),ZERO,x, &
     &                              2,scaloc,xnorm,ierr)
                        IF ( ierr/=0 ) Info = 1
!
                        IF ( scaloc/=ONE ) THEN
                           DO j = 1 , N
                              CALL DSCAL(M,scaloc,C(1,j),1)
                           ENDDO
                           Scale = Scale*scaloc
                        ENDIF
                        C(k1,l1) = x(1,1)
                        C(k2,l1) = x(2,1)
!
                     ELSEIF ( l1/=l2 .AND. k1==k2 ) THEN
!
                        suml = DDOT(k1-1,A(1,k1),1,C(1,l1),1)
                        sumr = DDOT(N-l2,C(k1,MIN(l2+1,N)),Ldc,         &
     &                         B(l1,MIN(l2+1,N)),Ldb)
                        vec(1,1) = sgn*(C(k1,l1)-(suml+sgn*sumr))
!
                        suml = DDOT(k1-1,A(1,k1),1,C(1,l2),1)
                        sumr = DDOT(N-l2,C(k1,MIN(l2+1,N)),Ldc,         &
     &                         B(l2,MIN(l2+1,N)),Ldb)
                        vec(2,1) = sgn*(C(k1,l2)-(suml+sgn*sumr))
!
                        CALL DLALN2(.FALSE.,2,1,smin,ONE,B(l1,l1),Ldb,  &
     &                              ONE,ONE,vec,2,-sgn*A(k1,k1),ZERO,x, &
     &                              2,scaloc,xnorm,ierr)
                        IF ( ierr/=0 ) Info = 1
!
                        IF ( scaloc/=ONE ) THEN
                           DO j = 1 , N
                              CALL DSCAL(M,scaloc,C(1,j),1)
                           ENDDO
                           Scale = Scale*scaloc
                        ENDIF
                        C(k1,l1) = x(1,1)
                        C(k1,l2) = x(2,1)
!
                     ELSEIF ( l1/=l2 .AND. k1/=k2 ) THEN
!
                        suml = DDOT(k1-1,A(1,k1),1,C(1,l1),1)
                        sumr = DDOT(N-l2,C(k1,MIN(l2+1,N)),Ldc,         &
     &                         B(l1,MIN(l2+1,N)),Ldb)
                        vec(1,1) = C(k1,l1) - (suml+sgn*sumr)
!
                        suml = DDOT(k1-1,A(1,k1),1,C(1,l2),1)
                        sumr = DDOT(N-l2,C(k1,MIN(l2+1,N)),Ldc,         &
     &                         B(l2,MIN(l2+1,N)),Ldb)
                        vec(1,2) = C(k1,l2) - (suml+sgn*sumr)
!
                        suml = DDOT(k1-1,A(1,k2),1,C(1,l1),1)
                        sumr = DDOT(N-l2,C(k2,MIN(l2+1,N)),Ldc,         &
     &                         B(l1,MIN(l2+1,N)),Ldb)
                        vec(2,1) = C(k2,l1) - (suml+sgn*sumr)
!
                        suml = DDOT(k1-1,A(1,k2),1,C(1,l2),1)
                        sumr = DDOT(N-l2,C(k2,MIN(l2+1,N)),Ldc,         &
     &                         B(l2,MIN(l2+1,N)),Ldb)
                        vec(2,2) = C(k2,l2) - (suml+sgn*sumr)
!
                        CALL DLASY2(.TRUE.,.TRUE.,Isgn,2,2,A(k1,k1),Lda,&
     &                              B(l1,l1),Ldb,vec,2,scaloc,x,2,xnorm,&
     &                              ierr)
                        IF ( ierr/=0 ) Info = 1
!
                        IF ( scaloc/=ONE ) THEN
                           DO j = 1 , N
                              CALL DSCAL(M,scaloc,C(1,j),1)
                           ENDDO
                           Scale = Scale*scaloc
                        ENDIF
                        C(k1,l1) = x(1,1)
                        C(k1,l2) = x(1,2)
                        C(k2,l1) = x(2,1)
                        C(k2,l2) = x(2,2)
                     ENDIF
                  ENDIF
!
               ENDDO
            ENDIF
         ENDDO
!
      ELSEIF ( notrna .AND. .NOT.notrnb ) THEN
!
!        Solve    A*X + ISGN*X*B**T = scale*C.
!
!        The (K,L)th block of X is determined starting from
!        bottom-right corner column by column by
!
!            A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L)
!
!        Where
!                      M                          N
!            R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T].
!                    I=K+1                      J=L+1
!
!        Start column loop (index = L)
!        L1 (L2): column index of the first (last) row of X(K,L)
!
         lnext = N
         DO l = N , 1 , -1
            IF ( l<=lnext ) THEN
               IF ( l==1 ) THEN
                  l1 = l
                  l2 = l
               ELSEIF ( B(l,l-1)/=ZERO ) THEN
                  l1 = l - 1
                  l2 = l
                  lnext = l - 2
               ELSE
                  l1 = l
                  l2 = l
                  lnext = l - 1
               ENDIF
!
!           Start row loop (index = K)
!           K1 (K2): row index of the first (last) row of X(K,L)
!
               knext = M
               DO k = M , 1 , -1
                  IF ( k<=knext ) THEN
                     IF ( k==1 ) THEN
                        k1 = k
                        k2 = k
                     ELSEIF ( A(k,k-1)/=ZERO ) THEN
                        k1 = k - 1
                        k2 = k
                        knext = k - 2
                     ELSE
                        k1 = k
                        k2 = k
                        knext = k - 1
                     ENDIF
!
                     IF ( l1==l2 .AND. k1==k2 ) THEN
                        suml = DDOT(M-k1,A(k1,MIN(k1+1,M)),Lda,         &
     &                         C(MIN(k1+1,M),l1),1)
                        sumr = DDOT(N-l1,C(k1,MIN(l1+1,N)),Ldc,         &
     &                         B(l1,MIN(l1+1,N)),Ldb)
                        vec(1,1) = C(k1,l1) - (suml+sgn*sumr)
                        scaloc = ONE
!
                        a11 = A(k1,k1) + sgn*B(l1,l1)
                        da11 = ABS(a11)
                        IF ( da11<=smin ) THEN
                           a11 = smin
                           da11 = smin
                           Info = 1
                        ENDIF
                        db = ABS(vec(1,1))
                        IF ( da11<ONE .AND. db>ONE ) THEN
                           IF ( db>bignum*da11 ) scaloc = ONE/db
                        ENDIF
                        x(1,1) = (vec(1,1)*scaloc)/a11
!
                        IF ( scaloc/=ONE ) THEN
                           DO j = 1 , N
                              CALL DSCAL(M,scaloc,C(1,j),1)
                           ENDDO
                           Scale = Scale*scaloc
                        ENDIF
                        C(k1,l1) = x(1,1)
!
                     ELSEIF ( l1==l2 .AND. k1/=k2 ) THEN
!
                        suml = DDOT(M-k2,A(k1,MIN(k2+1,M)),Lda,         &
     &                         C(MIN(k2+1,M),l1),1)
                        sumr = DDOT(N-l2,C(k1,MIN(l2+1,N)),Ldc,         &
     &                         B(l1,MIN(l2+1,N)),Ldb)
                        vec(1,1) = C(k1,l1) - (suml+sgn*sumr)
!
                        suml = DDOT(M-k2,A(k2,MIN(k2+1,M)),Lda,         &
     &                         C(MIN(k2+1,M),l1),1)
                        sumr = DDOT(N-l2,C(k2,MIN(l2+1,N)),Ldc,         &
     &                         B(l1,MIN(l2+1,N)),Ldb)
                        vec(2,1) = C(k2,l1) - (suml+sgn*sumr)
!
                        CALL DLALN2(.FALSE.,2,1,smin,ONE,A(k1,k1),Lda,  &
     &                              ONE,ONE,vec,2,-sgn*B(l1,l1),ZERO,x, &
     &                              2,scaloc,xnorm,ierr)
                        IF ( ierr/=0 ) Info = 1
!
                        IF ( scaloc/=ONE ) THEN
                           DO j = 1 , N
                              CALL DSCAL(M,scaloc,C(1,j),1)
                           ENDDO
                           Scale = Scale*scaloc
                        ENDIF
                        C(k1,l1) = x(1,1)
                        C(k2,l1) = x(2,1)
!
                     ELSEIF ( l1/=l2 .AND. k1==k2 ) THEN
!
                        suml = DDOT(M-k1,A(k1,MIN(k1+1,M)),Lda,         &
     &                         C(MIN(k1+1,M),l1),1)
                        sumr = DDOT(N-l2,C(k1,MIN(l2+1,N)),Ldc,         &
     &                         B(l1,MIN(l2+1,N)),Ldb)
                        vec(1,1) = sgn*(C(k1,l1)-(suml+sgn*sumr))
!
                        suml = DDOT(M-k1,A(k1,MIN(k1+1,M)),Lda,         &
     &                         C(MIN(k1+1,M),l2),1)
                        sumr = DDOT(N-l2,C(k1,MIN(l2+1,N)),Ldc,         &
     &                         B(l2,MIN(l2+1,N)),Ldb)
                        vec(2,1) = sgn*(C(k1,l2)-(suml+sgn*sumr))
!
                        CALL DLALN2(.FALSE.,2,1,smin,ONE,B(l1,l1),Ldb,  &
     &                              ONE,ONE,vec,2,-sgn*A(k1,k1),ZERO,x, &
     &                              2,scaloc,xnorm,ierr)
                        IF ( ierr/=0 ) Info = 1
!
                        IF ( scaloc/=ONE ) THEN
                           DO j = 1 , N
                              CALL DSCAL(M,scaloc,C(1,j),1)
                           ENDDO
                           Scale = Scale*scaloc
                        ENDIF
                        C(k1,l1) = x(1,1)
                        C(k1,l2) = x(2,1)
!
                     ELSEIF ( l1/=l2 .AND. k1/=k2 ) THEN
!
                        suml = DDOT(M-k2,A(k1,MIN(k2+1,M)),Lda,         &
     &                         C(MIN(k2+1,M),l1),1)
                        sumr = DDOT(N-l2,C(k1,MIN(l2+1,N)),Ldc,         &
     &                         B(l1,MIN(l2+1,N)),Ldb)
                        vec(1,1) = C(k1,l1) - (suml+sgn*sumr)
!
                        suml = DDOT(M-k2,A(k1,MIN(k2+1,M)),Lda,         &
     &                         C(MIN(k2+1,M),l2),1)
                        sumr = DDOT(N-l2,C(k1,MIN(l2+1,N)),Ldc,         &
     &                         B(l2,MIN(l2+1,N)),Ldb)
                        vec(1,2) = C(k1,l2) - (suml+sgn*sumr)
!
                        suml = DDOT(M-k2,A(k2,MIN(k2+1,M)),Lda,         &
     &                         C(MIN(k2+1,M),l1),1)
                        sumr = DDOT(N-l2,C(k2,MIN(l2+1,N)),Ldc,         &
     &                         B(l1,MIN(l2+1,N)),Ldb)
                        vec(2,1) = C(k2,l1) - (suml+sgn*sumr)
!
                        suml = DDOT(M-k2,A(k2,MIN(k2+1,M)),Lda,         &
     &                         C(MIN(k2+1,M),l2),1)
                        sumr = DDOT(N-l2,C(k2,MIN(l2+1,N)),Ldc,         &
     &                         B(l2,MIN(l2+1,N)),Ldb)
                        vec(2,2) = C(k2,l2) - (suml+sgn*sumr)
!
                        CALL DLASY2(.FALSE.,.TRUE.,Isgn,2,2,A(k1,k1),   &
     &                              Lda,B(l1,l1),Ldb,vec,2,scaloc,x,2,  &
     &                              xnorm,ierr)
                        IF ( ierr/=0 ) Info = 1
!
                        IF ( scaloc/=ONE ) THEN
                           DO j = 1 , N
                              CALL DSCAL(M,scaloc,C(1,j),1)
                           ENDDO
                           Scale = Scale*scaloc
                        ENDIF
                        C(k1,l1) = x(1,1)
                        C(k1,l2) = x(1,2)
                        C(k2,l1) = x(2,1)
                        C(k2,l2) = x(2,2)
                     ENDIF
                  ENDIF
!
               ENDDO
            ENDIF
         ENDDO
!
      ENDIF
!
!
!     End of DTRSYL
!
      END SUBROUTINE DTRSYL


text:
      SUBROUTINE hqr(a,n,np,wr,wi)
      INTEGER n,np
      REAL a(np,np),wi(np),wr(np)
      INTEGER i,its,j,k,l,m,nn
      REAL anorm,p,q,r,s,t,u,v,w,x,y,z
      anorm=0.
      do 12 i=1,n
        do 11 j=max(i-1,1),n
          anorm=anorm+abs(a(i,j))
11      continue
12    continue
      nn=n
      t=0.
1     if(nn.ge.1)then
        its=0
2       do 13 l=nn,2,-1
          s=abs(a(l-1,l-1))+abs(a(l,l))
          if(s.eq.0.)s=anorm
          if(abs(a(l,l-1))+s.eq.s)goto 3
13      continue
        l=1
3       x=a(nn,nn)
        if(l.eq.nn)then
          wr(nn)=x+t
          wi(nn)=0.
          nn=nn-1
        else
          y=a(nn-1,nn-1)
          w=a(nn,nn-1)*a(nn-1,nn)
          if(l.eq.nn-1)then
            p=0.5*(y-x)
            q=p**2+w
            z=sqrt(abs(q))
            x=x+t
            if(q.ge.0.)then
              z=p+sign(z,p)
              wr(nn)=x+z
              wr(nn-1)=wr(nn)
              if(z.ne.0.)wr(nn)=x-w/z
              wi(nn)=0.
              wi(nn-1)=0.
            else
              wr(nn)=x+p
              wr(nn-1)=wr(nn)
              wi(nn)=z
              wi(nn-1)=-z
            endif
            nn=nn-2
          else
            if(its.eq.30)pause 'too many iterations in hqr'
            if(its.eq.10.or.its.eq.20)then
              t=t+x
              do 14 i=1,nn
                a(i,i)=a(i,i)-x
14            continue
              s=abs(a(nn,nn-1))+abs(a(nn-1,nn-2))
              x=0.75*s
              y=x
              w=-0.4375*s**2
            endif
            its=its+1
            do 15 m=nn-2,l,-1
              z=a(m,m)
              r=x-z
              s=y-z
              p=(r*s-w)/a(m+1,m)+a(m,m+1)
              q=a(m+1,m+1)-z-r-s
              r=a(m+2,m+1)
              s=abs(p)+abs(q)+abs(r)
              p=p/s
              q=q/s
              r=r/s
              if(m.eq.l)goto 4
              u=abs(a(m,m-1))*(abs(q)+abs(r))
              v=abs(p)*(abs(a(m-1,m-1))+abs(z)+abs(a(m+1,m+1)))
              if(u+v.eq.v)goto 4
15          continue
4           do 16 i=m+2,nn
              a(i,i-2)=0.
              if (i.ne.m+2) a(i,i-3)=0.
16          continue
            do 19 k=m,nn-1
              if(k.ne.m)then
                p=a(k,k-1)
                q=a(k+1,k-1)
                r=0.
                if(k.ne.nn-1)r=a(k+2,k-1)
                x=abs(p)+abs(q)+abs(r)
                if(x.ne.0.)then
                  p=p/x
                  q=q/x
                  r=r/x
                endif
              endif
              s=sign(sqrt(p**2+q**2+r**2),p)
              if(s.ne.0.)then
                if(k.eq.m)then
                  if(l.ne.m)a(k,k-1)=-a(k,k-1)
                else
                  a(k,k-1)=-s*x
                endif
                p=p+s
                x=p/s
                y=q/s
                z=r/s
                q=q/p
                r=r/p
                do 17 j=k,nn
                  p=a(k,j)+q*a(k+1,j)
                  if(k.ne.nn-1)then
                    p=p+r*a(k+2,j)
                    a(k+2,j)=a(k+2,j)-p*z
                  endif
                  a(k+1,j)=a(k+1,j)-p*y
                  a(k,j)=a(k,j)-p*x
17              continue
                do 18 i=l,min(nn,k+3)
                  p=x*a(i,k)+y*a(i,k+1)
                  if(k.ne.nn-1)then
                    p=p+z*a(i,k+2)
                    a(i,k+2)=a(i,k+2)-p*r
                  endif
                  a(i,k+1)=a(i,k+1)-p*q
                  a(i,k)=a(i,k)-p
18              continue
              endif
19          continue
            goto 2
          endif
        endif
      goto 1
      endif
      return
      END

ext: for

max_stars_repo_path: NR-Functions/Numerical Recipes- Example & Functions/Functions/hqr.for


text:
MODULE gold_seaice_data

  USE gold_seaice_lib
  USE gold_seaice_netcdf
  IMPLICIT NONE

CONTAINS
  ! Read in netCDF restart files for GOLDSTEIN sea-ice
  SUBROUTINE inm_netcdf_sic
    USE genie_global, ONLY: write_status, gui_restart, istep_sic
    IMPLICIT NONE

    REAL, DIMENSION(maxi,maxj) :: hght_read, frac_read, temp_read, albd_read

    INTEGER :: iday, ifail, ncid
    LOGICAL :: lexist
    CHARACTER(LEN=200) :: fnamein
    REAL :: timestep
    INTEGER :: i, j, icell
    REAL :: tmp_val(4)

    timestep = 24.0 * 60.0 * 60.0 * yearlen / REAL(nyear)
    IF (gui_restart) THEN
       PRINT *, 'READING SEAICE GUI RESTART FILE: gui_restart_goldsteinseaice.nc'
       fnamein = 'gui_restart_goldsteinseaice.nc'
    ELSE
       fnamein = TRIM(filenetin)
    END IF
    ifail = 0
    INQUIRE(FILE=TRIM(fnamein), EXIST=lexist)
    IF (.NOT. lexist) THEN
       PRINT *, ' Missing file ', TRIM(fnamein)
       PRINT *, ' Correct error and try again '
       CALL write_status('ERRORED')
    END IF

    PRINT *, ' goldstein sea-ice: Opening restart file for read: ', &
         & TRIM(filenetin)

    CALL open_file_nc(TRIM(fnamein), ncid)
    IF (gui_restart) &
         & CALL getscalari_data_nc(ncid, 'istep_sic', istep_sic, ifail)
    CALL getscalari_data_nc(ncid,'ioffset', ioffset_rest, ifail)
    CALL getscalari_data_nc(ncid,'iyear',   iyear_rest,   ifail)
    CALL getscalari_data_nc(ncid,'imonth',  imonth_rest,  ifail)
    CALL getscalari_data_nc(ncid,'iday',    iday,         ifail)
    CALL get2d_data_nc(ncid,'sic_height', maxi, maxj, hght_read, ifail)
    CALL get2d_data_nc(ncid,'sic_cover',  maxi, maxj, frac_read, ifail)
    CALL get2d_data_nc(ncid,'sic_temp',   maxi, maxj, temp_read, ifail)
    CALL get2d_data_nc(ncid,'sic_albedo', maxi, maxj, albd_read, ifail)
    CALL close_file_nc(TRIM(filenetin),ncid)

    day_rest = iday
    ioffset_rest = MOD(ioffset_rest, NINT(yearlen))

    day_rest = day_rest + timestep / (24 * 60 * 60.0)
    ! This bit so that we don't get too far out in our count....
    ! Anchor to a day if we start drifting.
    ! Means timestep can never be less than 1/1000 of a day!!!!
    IF (ABS(iday - day_rest) <= 1.0E-3) day_rest = iday
    IF (day_rest >= 31) THEN
       day_rest = day_rest - 30
       imonth_rest = imonth_rest + 1
       IF (imonth_rest == 13) THEN
          imonth_rest = 1
          iyear_rest = iyear_rest + 1
       END IF
    END IF
    iday = NINT(day_rest)
    IF (debug_init) PRINT *, 'day in goldstein sea-ice restart is now', day_rest
    IF (debug_init) PRINT *, 'iday in goldstein sea-ice restart is now', iday

    varice(1,1:maxi,1:maxj) = hght_read(:,:)
    varice(2,1:maxi,1:maxj) = frac_read(:,:)
    tice(1:maxi,1:maxj) = temp_read(:,:)
    albice(1:maxi,1:maxj) = albd_read(:,:)

    ! Write out averages for restart checks
    IF (debug_init) &
         & WRITE (*,320) 'Avg height','Avg area', 'Avg T', 'Avg albedo'
    tmp_val = 0
    icell = 0
    DO j = 1, maxj
       DO i = 1, maxi
          IF (k1(i,j) <= maxk .AND. varice(2,i,j) > 0.0) THEN
             icell = icell + 1
             tmp_val(1:2) = tmp_val(1:2) + varice(:,i,j)
             tmp_val(3) = tmp_val(3) + tice(i,j)
             tmp_val(4) = tmp_val(4) + albice(i,j)
          END IF
       END DO
    END DO
    IF (icell > 0) THEN
       IF (debug_init) WRITE (*,310) tmp_val(1) / icell, tmp_val(2) / icell, &
            & tmp_val(3) / icell, tmp_val(4) / icell
    ELSE
       IF (debug_init) WRITE (*,310) 0.0, 0.0, 0.0, 0.0
    END IF

310 FORMAT(4f13.9)
320 FORMAT(4a13)
  END SUBROUTINE inm_netcdf_sic

  SUBROUTINE outm_netcdf_sic(istep)
    USE netcdf
    USE genie_global, ONLY: writing_gui_restarts
    IMPLICIT NONE

    INTEGER istep

    REAL, DIMENSION(maxi,maxj) :: hght_write, frac_write, temp_write, albd_write
    REAL :: lons1(maxi), lats1(maxj)
    INTEGER :: landmask(maxi,maxj)
    INTEGER :: i, j, nhghtid, nfracid, ntempid, nalbdid
    INTEGER :: nlon1id, nlongit1id, nlat1id, nlatit1id, nrecsid(1), ioffsetid
    INTEGER :: istepsicid
    INTEGER :: dim1pass(2)
    CHARACTER(LEN=200) :: fname

    INTEGER :: iday, iyearid, imonthid, idayid
    CHARACTER(LEN=10) :: yearstring
    CHARACTER(LEN=2) :: monthstring, daystring
    CHARACTER(LEN=7) :: datestring

    INTEGER :: ncid
    REAL :: timestep

    INTEGER :: icell
    REAL :: tmp_val(4)

    timestep = yearlen / REAL(nyear)

    ! output file format is yyyy_mm_dd

    iday = NINT(day_rest)
    IF (MOD(istep,iwstp) == 0 .OR. writing_gui_restarts) THEN
       ! WRITE A RESTART.....
       lons1 = nclon1
       lats1 = nclat1
       landmask = 0
       ! Next line contains an explicit reference to the values
       ! used in the topography grid (i.e. < 90) - not ideal,
       ! but the routine doesn't know about maxk
       WHERE (k1(1:maxi,1:maxj) < 90) landmask = 1

       hght_write = REAL(varice(1,:,:) * landmask)
       frac_write = REAL(varice(2,:,:) * landmask)
       temp_write = REAL(tice)
       albd_write = REAL(albice)

       WRITE (datestring,'(i7.7)') istep
       WRITE (yearstring,'(i10)') iyear_rest
       WRITE (monthstring,'(i2.2)') imonth_rest
       WRITE (daystring,'(i2.2)') iday

       IF (writing_gui_restarts) THEN
          fname = 'gui_restart_goldsteinseaice.nc'
       ELSE
          fname = TRIM(dirnetout) // '/goldsic_restart_' // &
               & TRIM(ADJUSTL(yearstring)) // '_' // monthstring // '_' // &
               & daystring // '.nc'
       END IF
       IF (debug_end) PRINT *, ' Opening netcdf restart file for write: ', TRIM(fname)
       CALL check_err(NF90_CREATE(TRIM(fname), NF90_CLOBBER, ncid))
       CALL check_err(NF90_DEF_DIM(ncid, 'nrecs', 1, nrecsid(1)))
       CALL check_err(NF90_DEF_DIM(ncid, 'longitude', maxi, nlon1id))
       CALL check_err(NF90_DEF_DIM(ncid, 'latitude', maxj, nlat1id))

       CALL check_err(NF90_DEF_VAR(ncid, 'longitude', NF90_REAL, &
            & (/ nlon1id /), nlongit1id))
       CALL check_err(NF90_DEF_VAR(ncid, 'latitude', NF90_REAL, &
            & (/ nlat1id /), nlatit1id))
       dim1pass(1) = nlon1id
       dim1pass(2) = nlat1id
       CALL check_err(NF90_DEF_VAR(ncid, 'ioffset', NF90_INT, &
            & nrecsid, ioffsetid))
       CALL check_err(NF90_DEF_VAR(ncid, 'iyear', NF90_INT, nrecsid, iyearid))
       CALL check_err(NF90_DEF_VAR(ncid, 'imonth', NF90_INT, nrecsid, imonthid))
       CALL check_err(NF90_DEF_VAR(ncid, 'iday', NF90_INT, nrecsid, idayid))
       IF (writing_gui_restarts) &
            & CALL check_err(NF90_DEF_VAR(ncid, 'istep_sic', NF90_INT, &
            &                             nrecsid, istepsicid))

       CALL check_err(NF90_DEF_VAR(ncid, 'sic_height', NF90_DOUBLE, &
            & dim1pass, nhghtid))
       CALL check_err(NF90_DEF_VAR(ncid, 'sic_cover', NF90_DOUBLE, &
            & dim1pass, nfracid))
       CALL check_err(NF90_DEF_VAR(ncid, 'sic_temp', NF90_DOUBLE, &
            & dim1pass, ntempid))
       CALL check_err(NF90_DEF_VAR(ncid, 'sic_albedo', NF90_DOUBLE, &
            & dim1pass, nalbdid))
       CALL check_err(NF90_ENDDEF(ncid))

       CALL check_err(NF90_PUT_VAR(ncid, iyearid, iyear_rest))
       CALL check_err(NF90_PUT_VAR(ncid, imonthid, imonth_rest))
       CALL check_err(NF90_PUT_VAR(ncid, idayid, iday))
       CALL check_err(NF90_PUT_VAR(ncid, ioffsetid, ioffset_rest))
       IF (writing_gui_restarts) &
            & CALL check_err(NF90_PUT_VAR(ncid, istepsicid, istep))

       CALL check_err(NF90_PUT_VAR(ncid, nlongit1id, lons1))
       CALL check_err(NF90_PUT_VAR(ncid, nlatit1id, lats1))
       CALL check_err(NF90_PUT_VAR(ncid, nhghtid, hght_write))
       CALL check_err(NF90_PUT_VAR(ncid, nfracid, frac_write))
       CALL check_err(NF90_PUT_VAR(ncid, ntempid, temp_write))
       CALL check_err(NF90_PUT_VAR(ncid, nalbdid, albd_write))

       CALL check_err(NF90_CLOSE(ncid))

       IF (debug_end) WRITE (*,320) 'Avg height','Avg area', 'Avg T', 'Avg albedo'
       tmp_val = 0
       icell = 0
       DO j = 1, maxj
          DO i = 1, maxi
             IF (k1(i,j) <= maxk .AND. varice(2,i,j) > 0.0) THEN
                icell = icell + 1
                tmp_val(1:2) = tmp_val(1:2) + varice(:,i,j)
                tmp_val(3) = tmp_val(3) + tice(i,j)
                tmp_val(4) = tmp_val(4) + albice(i,j)
             END IF
          END DO
       END DO
       IF (debug_end) then
          IF (icell > 0) THEN
             WRITE (*,310) tmp_val(1) / icell, tmp_val(2) / icell, &
                  & tmp_val(3) / icell, tmp_val(4) / icell
          ELSE
             WRITE (*,310) 0.0, 0.0, 0.0, 0.0
       END IF
       end if

310    FORMAT(4f13.9)
320    FORMAT(4a13)
    END IF

    day_rest = day_rest + timestep
    ! This bit so that we don't get too far out in our count....
    ! Anchor to a day if we start drifting.
    ! Means timestep can never be less than 1/1000 of a day!!!!
    IF (ABS(iday - day_rest) <= 1.0E-3) THEN
       PRINT *, 'CORRECTING TIME-LAG! in outm_netcdf in genie-goldstein', &
            & iday, day_rest
       day_rest = iday
    END IF
    IF (day_rest >= 31) THEN
       day_rest = day_rest - 30
       imonth_rest = imonth_rest + 1
       IF (imonth_rest == 13) THEN
          imonth_rest = 1
          iyear_rest = iyear_rest + 1
       END IF
    END IF
    iday = NINT(day_rest)
  END SUBROUTINE outm_netcdf_sic

  SUBROUTINE diagsic(istep)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: istep

    REAL, DIMENSION(2) :: hamax, hamin, haaav
    REAL :: ticemax, ticemin, ticeaav, albicemax, albicemin, albiceaav, area
    INTEGER :: i, j, l, icepts

    DO l = 1, 2
       hamax(l) = MAXVAL(varice(l,:,:))
       hamin(l) = MINVAL(varice(l,:,:))
    END DO
    ticemax = MAXVAL(tice)
    ticemin = MINVAL(tice)
    albicemax = MAXVAL(albice)
    albicemin = MINVAL(albice)
    haaav = 0.0
    ticeaav = 0.0
    albiceaav = 0.0
    icepts = 0
    area = 0.0

    DO j = 1, maxj
       DO i = 1, maxi
          IF (varice(2,i,j) > 0.0) THEN
             haaav(1) = haaav(1) + varice(1,i,j) * ds(j)
             haaav(2) = haaav(2) + varice(2,i,j) * ds(j)
             ticeaav = ticeaav + tice(i,j)
             albiceaav = albiceaav + albice(i,j)
             icepts = icepts + 1
             area = area + ds(j)
          END IF
       END DO
    END DO
    IF (icepts > 0) THEN
       haaav = haaav / area
       ticeaav = ticeaav / area
       albiceaav = albiceaav / area
    END IF
    PRINT *, 'Sea-ice properties at time ', istep
    PRINT *, 'Ice height : max = ', hamax(1), '; min = ', hamin(1), &
         & '; avg = ', haaav(1)
    PRINT *, 'Ice cover  : max = ', hamax(2), '; min = ', hamin(2), &
         & '; avg = ', haaav(2)
    PRINT *, 'Ice temp.  : max = ', ticemax, '; min = ', ticemin, &
         & '; avg = ', ticeaav
    PRINT *, 'Ice albd.  : max = ', albicemax, '; min = ', albicemin, &
         & '; avg = ', albiceaav
  END SUBROUTINE diagsic

  ! Extra diagnostic routine for c-goldstein v2 with seasonal cycle
  ! calculate average over nyear timesteps.
  SUBROUTINE diagosc_sic(istep, iout, ext, fx_delta, fw_delta)
    USE genie_util, ONLY: check_unit, check_iostat
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: istep, iout
    CHARACTER(LEN=3), INTENT(IN) :: ext
    REAL, DIMENSION(:,:), INTENT(IN) :: fx_delta, fw_delta

    INTEGER :: ios, l
    REAL :: rnyear, work((maxi+1) * (maxj+1))

    rnyear = 1.0 / nyear
    DO l = 1, 2
       haavg(l,:,:) = SUM(varice(l,:,:) * rnyear)
       dthaavg(l,:,:) = SUM(dtha(l,:,:) * rnyear)
    END DO
    ticeavg = SUM(tice * rnyear)
    albiceavg = SUM(albice * rnyear)
    fxdelavg = SUM(fx_delta * rnyear)
    fwdelavg = SUM(fw_delta * rnyear)

    IF (iout == 1) THEN
       PRINT *, 'Sea-ice : writing averaged data at istep ', istep

       ! write averaged data (a near-copy of outm.f) not a restart
       ! as such, therefore can write less accurate, more economical output
       CALL check_unit(2, __LINE__, __FILE__)
       OPEN(2, FILE=outdir_name(1:lenout) // lout // '.osc.' // ext, IOSTAT=ios)
       CALL check_iostat(ios, __LINE__, __FILE__)
       WRITE (2,10,IOSTAT=ios) haavg
       CALL check_iostat(ios, __LINE__, __FILE__)
       WRITE (2,10,IOSTAT=ios) ticeavg
       CALL check_iostat(ios, __LINE__, __FILE__)
       WRITE (2,10,IOSTAT=ios) albiceavg
       CALL check_iostat(ios, __LINE__, __FILE__)
       WRITE (2,10,IOSTAT=ios) dthaavg(1,:,:)
       CALL check_iostat(ios, __LINE__, __FILE__)
       WRITE (2,10,IOSTAT=ios) dthaavg(2,:,:)
       CALL check_iostat(ios, __LINE__, __FILE__)
       WRITE (2,10,IOSTAT=ios) fxdelavg
       CALL check_iostat(ios, __LINE__, __FILE__)
       WRITE (2,10,IOSTAT=ios) fwdelavg
       CALL check_iostat(ios, __LINE__, __FILE__)
       CLOSE(2,IOSTAT=ios)
       CALL check_iostat(ios, __LINE__, __FILE__)

       PRINT *, 'Writing sea-ice mean annual netCDF file at time', istep
       CALL ini_netcdf_sic(istep, 2)
       CALL write_netcdf_sic(k1, haavg, ticeavg, albiceavg, dthaavg, &
            & fxdelavg, fwdelavg, work, maxi, maxj, 2)
       CALL end_netcdf_sic(2)
       PRINT *

       ! increment average counter
       iav = iav + 1

       ! perform diagnostics on averaged data, either by rewriting other diag
       ! routines to accept data as argument, or by simply copying code,
       ! otherwise diagnose by integrating one (short) step from .avg file.
       PRINT *, 'Sea-ice : resetting averaged data arrays at step', istep
       haavg = 0.0
       ticeavg = 0.0
       albiceavg = 0.0
       dthaavg = 0.0
       dthaavg = 0.0
       fxdelavg = 0.0
       fwdelavg = 0.0
       PRINT *
    END IF

10  FORMAT(e14.7)
  END SUBROUTINE diagosc_sic

END MODULE gold_seaice_data

max_stars_repo_name: derpycode/cgenie

max_issues_repo_name: derpycode/cgenie

max_forks_repo_name: derpycode/cgenie


text:
      SUBROUTINE DT_SHMAKI(Na,Nca,Nb,Ncb,Ijp,Ppn,Mode)

C***********************************************************************
C Initialisation of Glauber formalism. This subroutine has to be       *
C called once (in case of target emulsions as often as many different  *
C target nuclei are considered) before events are sampled.             *
C         NA / NCA   mass number/charge of projectile nucleus          *
C         NB / NCB   mass number/charge of target     nucleus          *
C         IJP        identity of projectile (hadrons/leptons/photons)  *
C         PPN        projectile momentum (for projectile nuclei:       *
C                    momentum per nucleon) in target rest system       *
C         MODE = 0   Glauber formalism invoked                         *
C              = 1   fitted results are loaded from data-file          *
C              = 99  NTARG is forced to be 1                           *
C                    (used in connection with GLAUBERI-card only)      *
C This version dated 22.03.96 is based on the original SHMAKI-routine  *
C and revised by S. Roesler.                                           *
C***********************************************************************

      IMPLICIT NONE
      DOUBLE PRECISION aecmhi , aecmlo , amp , amp2 , amt , amt2 ,
     &                 aq2hi , aq2lo , daecm , daq2 , ecm , elab , ONE ,
     &                 plab , Ppn , THREE , TINY10 , TWO , ZERO
      INTEGER i , iback , ibin , icout , iestep , Ijp , iqstep ,
     &        iveout , j , Mode , Na , Nb , Nca , Ncb , nidx , ntarg
      SAVE

      INCLUDE 'inc/dtflka'

      PARAMETER (ZERO=0.0D0,TINY10=1.0D-10,ONE=1.0D0,TWO=2.0D0,
     &           THREE=3.0D0)

C emulsion treatment
      INCLUDE 'inc/dtcomp'
C Glauber formalism: parameters
      INCLUDE 'inc/dtglam'
C Lorentz-parameters of the current interaction
      INCLUDE 'inc/dtltra'
C properties of photon/lepton projectiles
      INCLUDE 'inc/dtgpro'
C kinematical cuts for lepton-nucleus interactions
      INCLUDE 'inc/dtlcut'
C Glauber formalism: cross sections
      INCLUDE 'inc/dtglxs'
C cuts for variable energy runs
      INCLUDE 'inc/dtvare'
C nucleon-nucleon event-generator
      INCLUDE 'inc/dtmodl'
C Glauber formalism: flags and parameters for statistics
      INCLUDE 'inc/dtglgp'

      DATA ntarg , icout , iveout/0 , 0 , 0/

C     CALL DT_HISHAD
C     STOP

      ntarg = ntarg + 1
      IF ( Mode.EQ.99 ) ntarg = 1
      nidx = -ntarg
      IF ( Mode.EQ.-1 ) nidx = ntarg

      IF ( (icout.LT.15) .AND. (MCGene.NE.4) ) icout = icout + 1

      IF ( LPRi.GT.4 .AND. icout.EQ.1 ) WRITE (LOUt,99010)
99010 FORMAT (//,1X,'SHMAKI:    Glauber formalism (Shmakov et. al) -',
     &        ' initialization',/,12X,'--------------------------',
     &        '-------------------------',/)

      IF ( Mode.EQ.2 ) THEN
         CALL DT_XSGLAU(Na,Nb,Ijp,ZERO,VIRt,UMO,1,1,nidx)
         CALL DT_SHFAST(Mode,Ppn,iback)
         STOP ' Glauber pre-initialization done'
      END IF
      IF ( Mode.EQ.1 ) THEN
         CALL DT_PROFBI(Na,Nb,Ppn,ntarg)
      ELSE
         iback = 1
         IF ( Mode.EQ.3 ) CALL DT_SHFAST(Mode,Ppn,iback)
         IF ( iback.EQ.1 ) THEN
C lepton-nucleus (variable energy runs)
            IF ( (Ijp.EQ.3) .OR. (Ijp.EQ.4) .OR. (Ijp.EQ.10) .OR.
     &           (Ijp.EQ.11) ) THEN

               IF ( (icout.LT.15) .AND. (MCGene.NE.4) .AND. LPRi.GT.4 )
     &              WRITE (LOUt,99020) Nb , Ncb
99020          FORMAT (1X,'variable energy run:     projectile-id:  7',
     &                 '    target A/Z: ',I3,' /',I3,/,/,8X,
     &                 'E_cm (GeV)    Q^2 (GeV^2)',
     &                 '    Sigma_tot (mb)     Sigma_in (mb)',/,7X,
     &                 '--------------------------------',
     &                 '------------------------------')
               aecmlo = LOG10(MIN(UMO,ECMli))
               aecmhi = LOG10(MIN(UMO,ECMhi))
               iestep = NEB - 1
               daecm = (aecmhi-aecmlo)/DBLE(iestep)
               IF ( aecmlo.EQ.aecmhi ) iestep = 0
               DO i = 1 , iestep + 1
                  ecm = 10.0D0**(aecmlo+DBLE(i-1)*daecm)
                  IF ( Q2Hi.GT.0.1D0 ) THEN
                     IF ( Q2Li.LT.0.01D0 ) THEN
                        CALL DT_XSGLAU(Na,Nb,7,ZERO,ZERO,ecm,i,1,nidx)

                        IF ( (icout.LT.15) .AND. (MCGene.NE.4) .AND.
     &                       LPRi.GT.4 ) WRITE (LOUt,99070) ECMnn(i) ,
     &                       ZERO , XSTot(i,1,ntarg) , XSPro(i,1,ntarg)
                        Q2Li = 0.01D0
                        ibin = 2
                     ELSE
                        ibin = 1
                     END IF
                     iqstep = NQB - ibin
                     aq2lo = LOG10(Q2Li)
                     aq2hi = LOG10(Q2Hi)
                     daq2 = (aq2hi-aq2lo)/MAX(DBLE(iqstep),ONE)
                     DO j = ibin , iqstep + ibin
                        Q2 = 10.0D0**(aq2lo+DBLE(j-ibin)*daq2)
                        CALL DT_XSGLAU(Na,Nb,7,ZERO,Q2,ecm,i,j,nidx)

                        IF ( (icout.LT.15) .AND. (MCGene.NE.4) .AND.
     &                       LPRi.GT.4 ) WRITE (LOUt,99070) ECMnn(i) ,
     &                       Q2G(j) , XSTot(i,j,ntarg) ,
     &                       XSPro(i,j,ntarg)
                     END DO
                  ELSE
                     CALL DT_XSGLAU(Na,Nb,7,ZERO,ZERO,ecm,i,1,nidx)

                     IF ( (icout.LT.15) .AND. (MCGene.NE.4) .AND.
     &                    LPRi.GT.4 ) WRITE (LOUt,99070) ECMnn(i) ,
     &                    ZERO , XSTot(i,1,ntarg) , XSPro(i,1,ntarg)
                  END IF
               END DO
               iveout = 1
C hadron/photon/nucleus-nucleus
            ELSE IF ( (ABS(VARehi).GT.ZERO) .AND.
     &                (ABS(VARehi).GT.ABS(VARelo)) ) THEN
               IF ( (icout.LT.15) .AND. (MCGene.NE.4) ) THEN

                  IF ( LPRi.GT.4 ) WRITE (LOUt,99030) Na , Nb , Ncb
99030             FORMAT (1X,'variable energy run:    projectile-id:',
     &                    I3,'    target A/Z: ',I3,' /',I3,/)

                  IF ( LPRi.GT.4 ) WRITE (LOUt,99040)
99040             FORMAT ('  E_cm (GeV)  E_Lab (GeV)  sig_tot^pp (mb)',
     &                    '  Sigma_tot (mb)  Sigma_prod (mb)',/,
     &                    ' -------------------------------------',
     &                    '--------------------------------------')
               END IF
               aecmlo = LOG10(VARclo)
               aecmhi = LOG10(VARchi)
               iestep = NEB - 1
               daecm = (aecmhi-aecmlo)/DBLE(iestep)
               IF ( aecmlo.EQ.aecmhi ) iestep = 0
               DO i = 1 , iestep + 1
                  ecm = 10.0D0**(aecmlo+DBLE(i-1)*daecm)
                  amp = 0.938D0
                  amt = 0.938D0
                  amp2 = amp**2
                  amt2 = amt**2
                  elab = (ecm**2-amp2-amt2)/(TWO*amt)
                  plab = SQRT((elab+amp)*(elab-amp))
                  CALL DT_XSGLAU(Na,Nb,Ijp,ZERO,VIRt,ecm,i,1,nidx)

                  IF ( (icout.LT.15) .AND. (MCGene.NE.4) .AND.
     &                 LPRi.GT.4 ) WRITE (LOUt,99050) ecm , plab ,
     &                 SIGsh , XSTot(i,1,ntarg) , XSPro(i,1,ntarg)
99050             FORMAT (1X,F9.1,1X,E11.3,1X,F12.2,8X,F10.3,8X,F8.3)
               END DO
               iveout = 1
            ELSE
               CALL DT_XSGLAU(Na,Nb,Ijp,ZERO,VIRt,UMO,1,1,nidx)
            END IF
         END IF
      END IF

      IF ( (icout.LT.15) .AND. (iveout.EQ.0) .AND. (MCGene.NE.4) .AND.
     &     (IOGlb.NE.100) ) THEN

         IF ( LPRi.GT.4 ) WRITE (LOUt,99060) Na , Nca , Nb , Ncb ,
     &        ECMnn(1) , SIGsh*10.0D0 , ROSh , BSLope , NSIteb ,
     &        NSTatb , XSPro(1,1,ntarg)
99060    FORMAT (38X,'projectile','      target',/,1X,
     &           'Mass number / charge',17X,I3,' /',I3,6X,I3,' /',I3,/,
     &           /,1X,'Nucleon-nucleon c.m. energy',9X,F10.2,' GeV',/,/,
     &           1X,'Parameters of elastic scattering amplitude:',/,5X,
     &           'sigma =',F7.2,' mb',6X,'rho = ',F9.4,6X,'slope = ',
     &           F4.1,' GeV^-2',/,/,1X,'Number of b-steps',4X,I3,8X,
     &           'statistics at each b-step',4X,I5,/,/,1X,
     &           'Prod. cross section  ',5X,F10.4,' mb',/)
      END IF
99070 FORMAT (9X,F6.1,9X,F6.2,8X,F8.3,11X,F8.3)

      END SUBROUTINE

ext: f

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


! { dg-do compile }
! { dg-options "-fcoarray=single" }
!
! PR fortran/51632
!
! Was rejected before as __def_init and __copy were
! resolved and coarray components aren't valid in this
! context
!
module periodic_2nd_order_module
  implicit none

  type periodic_2nd_order
    real, allocatable :: global_f(:)[:]
  contains
    procedure :: output
  end type

contains
  subroutine output (this)
    class(periodic_2nd_order), intent(in) :: this
  end subroutine
end module


text:
test2/test.cpp
test2/tb.cpp
test2/main.cpp


text:
C***BEGIN PROLOGUE  DWEB
C***REFER TO  DDASPK
C***DATE WRITTEN   950914   (YYMMDD)
C
C***AUTHORS  A. C. Hindmarsh, P. N. Brown
C            Lawrence Livermore National Laboratory
C            L-316, P.O. Box 808
C            Livermore, CA 94551, USA
C
C***DESCRIPTION
C
C-----------------------------------------------------------------------
C Example program for DDASPK.
C DAE system derived from ns-species interaction PDE in 2 dimensions.
C
C This is the double precision version.
C-----------------------------------------------------------------------
C
C This program solves a DAE system that arises from a system
C of partial differential equations.  The PDE system is a food web
C population model, with predator-prey interaction and diffusion on
C the unit square in two dimensions.  The dependent variable vector is
C
C         1   2        ns
C   c = (c , c , ..., c  )
C
C and the PDEs are as follows..
C
C     i               i      i
C   dc /dt  =  d(i)*(c    + c   )  +  R (x,y,c)  (i=1,...,ns/2)
C                     xx     yy        i
C
C                     i      i
C   0       =  d(i)*(c    + c   )  +  R (x,y,c)  (i=(ns/2)+1,...,ns)
C                     xx     yy        i
C
C where
C                  i          ns         j
C   R (x,y,c)  =  c *(b(i) + sum a(i,j)*c )
C    i                       j=1
C
C The number of species is ns = 2*np, with the first np being prey and
C the last np being predators.  The coefficients a(i,j), b(i), d(i) are
C
C   a(i,i) = -a  (all i)
C   a(i,j) = -g  (i .le. np, j .gt. np)
C   a(i,j) =  e  (i .gt. np, j .le. np)
C   b(i) =  b*(1 + alpha*x*y + beta*sin(4*pi*x)*sin(4*pi*y))  (i .le. np)
C   b(i) = -b*(1 + alpha*x*y + beta*sin(4*pi*x)*sin(4*pi*y))  (i .gt. np)
C   d(i) = dprey  (i .le. np)
C   d(i) = dpred  (i .gt. np)
C
C The various scalar parameters are set in subroutine setpar.
C
C The boundary conditions are.. normal derivative = 0.
C A polynomial in x and y is used to set the initial conditions.
C
C The PDEs are discretized by central differencing on a MX by MY mesh.
C
C The DAE system is solved by DDASPK with three different method options:
C (1) direct band method for the linear systems (internal Jacobian),
C (2) preconditioned Krylov method for the linear systems, without
C     block-grouping in the reaction-based factor, and
C (3) preconditioned Krylov method for the linear systems, with
C     block-grouping in the reaction-based factor.
C
C In the Krylov cases, the preconditioner is the product of two factors:
C (a) The spatial factor uses a fixed number of Gauss-Seidel iterations
C based on the diffusion terms only.
C (b) The reaction-based factor is a block-diagonal matrix based on
C the partial derivatives of the interaction terms R only.
C With block-grouping, only a subset of the ns by ns blocks are computed.
C An integer flag, JPRE, is set in the main program to specify whether
C the preconditioner is to use only one of the two factors or both,
C and in which order.
C
C The reaction-based preconditioner factor is set up and solved in
C seven subroutines --
C   DMSET2, DRBDJA, DRBDPS  in the case of no block-grouping, and
C   DGSET2, GSET1, DRBGJA, DRBGPS  in the case of block-grouping.
C These routines are provided separately for general use on problems
C arising from reaction-transport systems.
C
C Two output files are written.. one with the problem description and
C performance statistics on unit LOUT = 9, and one with solution
C profiles at selected output times on unit LCOUT = 10.
C The solution file is written only in the case of the direct method.
C-----------------------------------------------------------------------
C Note.. in addition to the main program and subroutines given below,
C this program requires the BLAS routine DAXPY.
C-----------------------------------------------------------------------
C References
C [1] Peter N. Brown and Alan C. Hindmarsh,
C     Reduced Storage Matrix Methods in Stiff ODE Systems,
C     J. Appl. Math. & Comp., 31 (1989), pp. 40-91.
C [2] Peter N. Brown, Alan C. Hindmarsh, and Linda R. Petzold,
C     Using Krylov Methods in the Solution of Large-Scale Differential-
C     Algebraic Systems, SIAM J. Sci. Comput., 15 (1994), pp. 1467-1488.
C [3] Peter N. Brown, Alan C. Hindmarsh, and Linda R. Petzold,
C     Consistent Initial Condition Calculation for Differential-
C     Algebraic Systems, LLNL Report UCRL-JC-122175, August 1995;
C     submitted to SIAM J. Sci. Comp.
C-----------------------------------------------------------------------
C***ROUTINES CALLED
C   SETPAR, DGSET2, CINIT, DDASPK, OUTWEB
C
C***END PROLOGUE  DWEB
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      EXTERNAL RESWEB, JACRS, PSOLRS
C
C Set output unit numbers for main output and tabulated solution.
      DATA LOUT/9/, LCOUT/10/
C
C Dimension solution arrays and work arrays.
C
C When INFO(12) = 0, with INFO(5) = 0, INFO(6) = 1:
C    The length required for RWORK is
C              50 + (2*ML+MU+11)*NEQ + 2*(NEQ/(ML+MU+1) + 1) .
C    For MX = MY = (even number) and ML = MU = NS*MX, this length is
C              50 + (3*NS*MX + 11)*NEQ + MY .
C    The length required for IWORK is  40 + 2*NEQ .
C
C When INFO(12) = 1:
C    The length required for RWORK is
C              91 + 19*NEQ + LENWP = 91 + 19*NEQ + NS*NS*NGRP .
C    The length required for IWORK is
C              40 + NEQ + LENIWP = 40 + NEQ + NS*NGRP .
C
C The dimensions for the various arrays are set below using parameters
C   MAXN    which must be .ge. NEQ = NS*MX*MY,
C   MAXS    which must be .ge. NS,
C   MAXM    which must be .ge. MAX(MX,MY).
C
      PARAMETER (MAXS = 2, MAXM = 20, MAXN = 800,
     1           LRW = 50 + (3*MAXS*MAXM + 11)*MAXN + MAXM,
     2           LIW = 40 + 2*MAXN)
C
      DIMENSION CC(MAXN), CCPRIME(MAXN), RWORK(LRW), IWORK(LIW),
     1          INFO(20), RPAR(MAXN), IPAR(2)
C
C The COMMON blocks /PPAR1/ and /PPAR2/ contain problem parameters.
C
      COMMON /PPAR1/ AA, EE, GG, BB, DPREY, DPRED
      COMMON /PPAR2/ NP, NS, AX, AY, ACOEF(MAXS,MAXS), BCOEF(MAXS),
     1               DX, DY, ALPH, BETA, FPI, DIFF(MAXS),
     2               COX(MAXS), COY(MAXS), MX, MY, MXNS
C
C Open output files.
      OPEN(UNIT=LOUT,FILE='wdout',STATUS='unknown')
      OPEN(UNIT=LCOUT,FILE='wccout',STATUS='unknown')
C
C Call SETPAR to set basic problem parameters.
      CALL SETPAR
C
C Set remaining problem parameters.
      NEQ = NS*MX*MY
      MXNS = MX*NS
      DX = AX/REAL(MX-1)
      DY = AY/REAL(MY-1)
      DO 10 I = 1,NS
        COX(I) = DIFF(I)/DX**2
 10     COY(I) = DIFF(I)/DY**2
C
      WRITE(LOUT,20)NS
 20   FORMAT(//' Example program for DDASPK package'//
     1   ' Food web problem with NS species, NS =',I4/
     2   ' Predator-prey interaction and diffusion on a 2-D square'/)
      WRITE(LOUT,30) AA,EE,GG,BB,DPREY,DPRED, ALPH,BETA
 30   FORMAT(' Matrix parameters..  a =',E12.4,'   e =',E12.4,
     1   '   g =',E12.4/21x,' b parameter =',E12.4//
     2   ' Diffusion coefficients.. dprey =',E12.4,'   dpred =',E12.4//
     3   ' Rate parameters alpha =',E12.4,' and beta =',E12.4/)
      WRITE(LOUT,40) MX,MY,NEQ
 40   FORMAT(' Mesh dimensions (MX,MY) =',2I4,
     1   5x,' Total system size is NEQ =',I7/)
C
C Here set the flat initial guess for the predators.
      PREDIC = 1.0D5
C
C Set remaining method parameters for DDASPK.
C These include the INFO array and tolerances.
C
      DO 50 I = 1,20
 50     INFO(I) = 0
C
C Here set INFO(11) = 1, indicating I.C. calculation requested.
      INFO(11) = 1
C
C Here set INFO(14) = 1 to get the computed initial values.
      INFO(14) = 1
C
C Here set INFO(15) = 1 to signal that a preconditioner setup routine
C is to be called in the Krylov case.
      INFO(15) = 1
C
C Here set INFO(16) = 1 to get alternative error test (on the
C differential variables only).
      INFO(16) = 1
C
C Here set the tolerances.
      RTOL = 1.0D-5
      ATOL = RTOL
C
      WRITE(LOUT,70)RTOL,ATOL,INFO(11),PREDIC,INFO(16)
 70   FORMAT(' Tolerance parameters.. RTOL =',E10.2,'   ATOL =',E10.2//
     1   ' Internal I.C. calculation flag INFO(11) =',I2,
     2   '   (0 = off, 1 = on)'/
     3   ' Predator I.C. guess =',E10.2//
     4   ' Alternate error test flag INFO(16) =',I2,
     5        '  (0 = off, 1 = on)')
C
C Set NOUT = number of output times.
      NOUT = 18
C
C Loop over method options:
C METH = 0 means use INFO(12) = 0 (direct)
C METH = 1 means use INFO(12) = 1 (Krylov) without block-grouping in
C          the reaction-based factor in the preconditioner.
C METH = 2 means use INFO(12) = 1 (Krylov) with block-grouping in
C          the reaction-based factor in the preconditioner.
C A block-grouping flag JBG, communicated through IPAR, is set to
C 0 (no block-grouping) or 1 (use block-grouping) with METH = 1 or 2.
C Reset INFO(1) = 0 and INFO(11) = 1.
C
      DO 300 METH = 0,2
      INFO(12) = MIN(METH,1)
      INFO(1) = 0
      INFO(11) = 1
      JBG = METH - 1
      IPAR(2) = JBG
C
      WRITE(LOUT,80)INFO(12)
 80   FORMAT(//80('.')//
     1   ' Linear solver method flag INFO(12) =',I2,
     2   '   (0 = direct, 1 = Krylov)'/)
C
C In the case of the direct method, set INFO(6) = 1 to signal a banded
C Jacobian, set IWORK(1) = IWORK(2) = MX*NS, the half-bandwidth, and
C call SETID to set the IWORK segment ID indicating the differential
C and algebraic components.
      IF (INFO(12) .EQ. 0) THEN
        INFO(6) = 1
        IWORK(1) = MXNS
        IWORK(2) = MXNS
        CALL SETID (MX, MY, NS, NP, 40, IWORK)
        WRITE(LOUT,90)MXNS
 90     FORMAT(' Difference-quotient banded Jacobian,'
     1         ' half-bandwidths =',I4)
        ENDIF
C
C In the case of the Krylov method, set and print various
C preconditioner parameters.
      IF (INFO(12) .EQ. 1) THEN
C First set the preconditioner choice JPRE.
C  JPRE = 1 means reaction-only (block-diagonal) factor A_R
C  JPRE = 2 means spatial factor (Gauss-Seidel) A_S
C  JPRE = 3 means A_S * A_R
C  JPRE = 4 means A_R * A_S
C Use IPAR to communicate JPRE to the preconditioner solve routine.
        JPRE = 3
        IPAR(1) = JPRE
        WRITE(LOUT,100)JPRE
 100    FORMAT(' Preconditioner flag is JPRE =',I3/
     1      '  (1 = reaction factor A_R, 2 = spatial factor A_S,',
     2      ' 3 = A_S*A_R, 4 = A_R*A_S )'/)
C Here call DMSET2 if JBG = 0, or DGSET2 if JBG = 1, to set the 2-D
C mesh parameters and block-grouping data, and the IWORK segment ID
C indicating the differential and algebraic components.
        IF (JBG .EQ. 0) THEN
          CALL DMSET2 (MX, MY, NS, NP, 40, IWORK)
          WRITE(LOUT,110)
 110      FORMAT(' No block-grouping in reaction factor')
          ENDIF
        IF (JBG .EQ. 1) THEN
          NXG = 5
          NYG = 5
          NG = NXG*NYG
          CALL DGSET2 (MX, MY, NS, NP, NXG, NYG, 40, IWORK)
          WRITE(LOUT,120)NG,NXG,NYG
 120      FORMAT(' Block-grouping in reaction factor'/
     1           ' Number of groups =',I5,
     2           '   (NGX by NGY, NGX =',I3,',  NGY =',I3,')')
          ENDIF
        ENDIF
C
C Set the initial T and TOUT, and call CINIT to set initial values.
      T = 0.0D0
      TOUT = 1.0D-8
      CALL CINIT (CC, CCPRIME, PREDIC, RPAR)
C
      NLI = 0
      NNI = 0
C
      WRITE(LOUT,140)
 140  FORMAT(//'   t',7X,'NSTEP  NRE  NNI  NLI  NPE  NQ',4X,
     1   'H',10X,'AVLIN')
C
C Loop over output times, call DDASPK, and print performance data.
C The first call, with IOUT = 0, is to calculate initial values only.
C After the first call, reset INFO(11) = 0 and the initial TOUT.
C
      DO 200 IOUT = 0,NOUT
C
       CALL DDASPK (RESWEB, NEQ, T, CC, CCPRIME, TOUT, INFO, RTOL, ATOL,
     1            IDID, RWORK,LRW, IWORK,LIW, RPAR, IPAR, JACRS, PSOLRS)
C
        NST = IWORK(11)
        NRE = IWORK(12)
        NPE = IWORK(13)
        NNIDIF = IWORK(19) - NNI
        NNI = IWORK(19)
        NLIDIF = IWORK(20) - NLI
        NLI = IWORK(20)
        NQU = IWORK(8)
        HU = RWORK(7)
        AVLIN = 0.0D0
        IF (NNIDIF .GT. 0) AVLIN = REAL(NLIDIF)/REAL(NNIDIF)
C
        IF (METH .EQ. 0) THEN
          IMOD3 = IOUT - 3*(IOUT/3)
          IF (IMOD3 .EQ. 0) CALL OUTWEB (T, CC, NS, MX, MY, LCOUT)
          ENDIF
C
        WRITE(LOUT,150)T,NST,NRE,NNI,NLI,NPE,NQU,HU,AVLIN
 150    FORMAT(E10.2,I5,I6,3I5,I4,E11.2,F9.4)
C
        IF (IDID .LT. 0) THEN
          WRITE(LOUT,160)T
 160      FORMAT(//' Final time reached =',E12.4//)
          GO TO 210
          ENDIF
C
        IF (TOUT .GT. 0.9D0) TOUT = TOUT + 1.0D0
        IF (TOUT .LT. 0.9D0) TOUT = TOUT*10.0D0
        IF (IOUT .EQ. 0) THEN
          INFO(11) = 0
          TOUT = 1.0D-8
          NLI = 0
          NNI = 0
          ENDIF
 200    CONTINUE
C
 210  CONTINUE
      LENRW = IWORK(18)
      LENIW = IWORK(17)
      NST = IWORK(11)
      NRE = IWORK(12)
      NPE = IWORK(13)
      NNI = IWORK(19)
      NLI = IWORK(20)
      NPS = IWORK(21)
      IF (NNI .GT. 0) AVLIN = REAL(NLI)/REAL(NNI)
      NCFN = IWORK(15)
      NCFL = IWORK(16)
      WRITE(LOUT,220) LENRW,LENIW,NST,NRE,NPE,NPS,NNI,NLI,AVLIN,
     1               NCFN,NCFL
 220   FORMAT(//' Final statistics for this run..'/
     1   ' RWORK size =',I8,'   IWORK size =',I6/
     2   ' Number of time steps            =',I5/
     3   ' Number of residual evaluations  =',I5/
     4   ' Number of Jac. or prec. evals.  =',I5/
     5   ' Number of preconditioner solves =',I5/
     6   ' Number of nonlinear iterations  =',I5/
     7   ' Number of linear iterations     =',I5/
     8   ' Average Krylov subspace dimension =',F8.4/
     9   I3,' nonlinear conv. failures,',i5,' linear conv. failures')
C
 300   CONTINUE

      STOP
C------  End of main program for DWEB example program ------------------
      END

      SUBROUTINE SETPAR
C-----------------------------------------------------------------------
C This routine sets the basic problem parameters, namely
C AX, AY, NS, MX, MY,  problem coefficients ACOEF, BCOEF, DIFF,
C ALPH, BETA, using parameters NP, AA, EE, GG, BB, DPREY, DPRED.
C-----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MAXS = 2)
      COMMON /PPAR1/ AA, EE, GG, BB, DPREY, DPRED
      COMMON /PPAR2/ NP, NS, AX, AY, ACOEF(MAXS,MAXS), BCOEF(MAXS),
     1               DX, DY, ALPH, BETA, FPI, DIFF(MAXS),
     2               COX(MAXS), COY(MAXS), MX, MY, MXNS
C
      AX = 1.0D0
      AY = 1.0D0
      NP = 1
      MX = 20
      MY = 20
      AA = 1.0D0
      EE = 1.0D4
      GG = 0.5D-6
      BB = 1.0D0
      DPREY = 1.0D0
      DPRED = 0.05D0
      ALPH = 50.0D0
      BETA = 100.0D0
      NS = 2*NP
      DO 20 J = 1,NP
        DO 10 I = 1,NP
          ACOEF(NP+I,J) = EE
          ACOEF(I,NP+J) = -GG
 10       CONTINUE
        ACOEF(J,J) = -AA
        ACOEF(NP+J,NP+J) = -AA
        BCOEF(J) = BB
        BCOEF(NP+J) = -BB
        DIFF(J) = DPREY
        DIFF(NP+J) = DPRED
 20     CONTINUE
      PI = 3.141592653589793D0
      FPI = 4.0D0*PI
C
      RETURN
C------------  End of Subroutine SETPAR  -------------------------------
      END

      SUBROUTINE SETID (MX, MY, NS, NSD, LID, IWORK)
C-----------------------------------------------------------------------
C This routine sets the ID array in IWORK, indicating which components
C are differential and which are algebraic.
C-----------------------------------------------------------------------
      DIMENSION IWORK(*)
C
      NSDP1 = NSD + 1
      DO 40 JY = 1,MY
        I00 = MX*NS*(JY-1) + LID
        DO 30 JX = 1,MX
          I0 = I00 + NS*(JX-1)
          DO 10 I = 1,NSD
 10         IWORK(I0+I) = 1
          DO 20 I = NSDP1,NS
 20         IWORK(I0+I) = -1
 30       CONTINUE
 40     CONTINUE
C
      RETURN
C------------  End of Subroutine SETID  --------------------------------
      END

      SUBROUTINE CINIT (CC, CCPRIME, PREDIC, RPAR)
C-----------------------------------------------------------------------
C This routine computes and loads the vectors of initial values.
C-----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION CC(*), CCPRIME(*), RPAR(*)
      PARAMETER (MAXS = 2)
      COMMON /PPAR2/ NP, NS, AX, AY, ACOEF(MAXS,MAXS), BCOEF(MAXS),
     1               DX, DY, ALPH, BETA, FPI, DIFF(MAXS),
     2               COX(MAXS), COY(MAXS), MX, MY, MXNS
C
C Load CC.
      NPP1 = NP + 1
      DO 30 JY = 1,MY
        Y = REAL(JY-1)*DY
        ARGY = 16.0D0*Y*Y*(AY-Y)*(AY-Y)
        IYOFF = MXNS*(JY-1)
        DO 20 JX = 1,MX
          X = REAL(JX-1)*DX
          ARGX = 16.0D0*X*X*(AX-X)*(AX-X)
          IOFF = IYOFF + NS*(JX-1)
          FAC = 1.0D0 + ALPH*X*Y + BETA*SIN(FPI*X)*SIN(FPI*Y)
          DO 10 I = 1,NP
 10         CC(IOFF + I) = 10.0D0 + REAL(I)*ARGX*ARGY
          DO 15 I = NPP1,NS
 15         CC(IOFF + I) = PREDIC
 20       CONTINUE
 30     CONTINUE
C
C Load CCPRIME.
      T = 0.0D0
      CALL FWEB (T, CC, CCPRIME, RPAR)
      DO 60 JY = 1,MY
        IYOFF = MXNS*(JY-1)
        DO 50 JX = 1,MX
          IOFF = IYOFF + NS*(JX-1)
          DO 40 I = NPP1,NS
 40         CCPRIME(IOFF+I) = 0.0D0
 50     CONTINUE
 60   CONTINUE
C
      RETURN
C------------  End of Subroutine CINIT  --------------------------------
      END

      SUBROUTINE OUTWEB (T, C, NS, MX, MY, LUN)
C-----------------------------------------------------------------------
C This routine prints the values of the individual species densities
C at the current time T, to logical unit LUN.
C-----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION C(NS,MX,MY)
C
      WRITE(LUN,10) T
 10   FORMAT(/1X,79('-')/30X,'At time t = ',E16.8/1X,79('-') )
C
      DO 40 I = 1,NS
        WRITE(LUN,20) I
 20     FORMAT(' the species c(',i2,') values are ')
        DO 30 JY = MY,1,-1
          WRITE(LUN,25) (C(I,JX,JY),JX=1,MX)
 25       FORMAT(6(1X,G12.6))
 30       CONTINUE
        WRITE(LUN,35)
 35     FORMAT(1X,79('-'),/)
 40     CONTINUE
C
      RETURN
C------------  End of Subroutine OUTWEB  -------------------------------
      END

      SUBROUTINE RESWEB (T, U, UPRIME, CJ, DELTA, IRES, RPAR, IPAR)
C-----------------------------------------------------------------------
C This routine computes the residual vector, using Subroutine FWEB
C for the right-hand sides.
C-----------------------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION U(*),UPRIME(*),DELTA(*),RPAR(*),IPAR(*)
      PARAMETER (MAXS = 2)
      COMMON /PPAR2/ NP, NS, AX, AY, ACOEF(MAXS,MAXS), BCOEF(MAXS),
     1               DX, DY, ALPH, BETA, FPI, DIFF(MAXS),
     2               COX(MAXS), COY(MAXS), MX, MY, MXNS
C
      CALL FWEB (T, U, DELTA, RPAR)
C
      DO 30 JY = 1,MY
         IYOFF = MXNS*(JY-1)
         DO 20 JX = 1,MX
            IC0 = IYOFF + NS*(JX-1)
            DO 10 I = 1,NS
               ICI = IC0 + I
               IF (I .GT. NP) THEN
                  DELTA(ICI) = -DELTA(ICI)
               ELSE
                  DELTA(ICI) = UPRIME(ICI) - DELTA(ICI)
               ENDIF
 10         CONTINUE
 20      CONTINUE
 30   CONTINUE
C
      RETURN
C------------  End of Subroutine RESWEB  -------------------------------
      END

      SUBROUTINE FWEB (T, CC, CRATE, RPAR)
C-----------------------------------------------------------------------
C This routine computes the right-hand sides of all the equations
C and returns them in the array CRATE.
C The interaction rates are computed by calls to WEBR, and these are
C saved in RPAR(1),...,RPAR(NEQ) for use later in preconditioning.
C-----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION CC(*), CRATE(*), RPAR(*)
      PARAMETER (MAXS = 2)
      COMMON /PPAR2/ NP, NS, AX, AY, ACOEF(MAXS,MAXS), BCOEF(MAXS),
     1               DX, DY, ALPH, BETA, FPI, DIFF(MAXS),
     2               COX(MAXS), COY(MAXS), MX, MY, MXNS
C
      DO 60 JY = 1,MY
        IYOFF = MXNS*(JY-1)
        IDYU = MXNS
        IF (JY .EQ. MY) IDYU = -MXNS
        IDYL = MXNS
        IF (JY .EQ. 1) IDYL = -MXNS
        DO 40 JX = 1,MX
          IC = IYOFF + NS*(JX-1) + 1
C Get interaction rates at one point (X,Y).
          CALL WEBR (T, JX, JY, CC(IC), RPAR(IC))
          IDXU = NS
          IF (JX .EQ. MX) IDXU = -NS
          IDXL = NS
          IF (JX .EQ. 1) IDXL = -NS
          DO 20 I = 1,NS
            ICI = IC + I - 1
C Do differencing in Y.
            DCYLI = CC(ICI) - CC(ICI-IDYL)
            DCYUI = CC(ICI+IDYU) - CC(ICI)
C Do differencing in X.
            DCXLI = CC(ICI) - CC(ICI-IDXL)
            DCXUI = CC(ICI+IDXU) - CC(ICI)
C Collect terms and load CRATE elements.
            CRATE(ICI) = COY(I)*(DCYUI - DCYLI) + COX(I)*(DCXUI - DCXLI)
     1                  + RPAR(ICI)
 20         CONTINUE
 40       CONTINUE
 60    CONTINUE
      RETURN
C------------  End of Subroutine FWEB  ---------------------------------
      END

      SUBROUTINE WEBR (T, JX, JY, C, CRATE)
C-----------------------------------------------------------------------
C This routine computes one block of the interaction term R of the
C system, namely block (JX,JY), for use in preconditioning.
C-----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION C(*), CRATE(*)
      PARAMETER (MAXS = 2)
      COMMON /PPAR2/ NP, NS, AX, AY, ACOEF(MAXS,MAXS), BCOEF(MAXS),
     1               DX, DY, ALPH, BETA, FPI, DIFF(MAXS),
     2               COX(MAXS), COY(MAXS), MX, MY, MXNS
C
      Y = REAL(JY-1)*DY
      X = REAL(JX-1)*DX
      DO 10 I = 1,NS
 10     CRATE(I) = 0.0D0
      DO 15 J = 1,NS
        CALL DAXPY (NS, C(J), ACOEF(1,J), 1, CRATE, 1)
 15     CONTINUE
      FAC = 1.0D0 + ALPH*X*Y + BETA*SIN(FPI*X)*SIN(FPI*Y)
      DO 20 I = 1,NS
 20     CRATE(I) = C(I)*(BCOEF(I)*FAC + CRATE(I))
      RETURN
C------------  End of Subroutine WEBR  ---------------------------------
      END

      SUBROUTINE JACRS(RES, IRES, NEQ, T, CC, CCPRIME, REWT, SAVR, WK,
     1                  H, CJ, WP, IWP, IER, RPAR, IPAR)
C-----------------------------------------------------------------------
C This routine interfaces to Subroutine DRBGJA or Subroutine DRBGJA,
C depending on the flag JBG = IPAR(2), to generate and preprocess the
C block-diagonal Jacobian corresponding to the reaction term R.
C If JBG = 0, we call DRBDJA, with no block-grouping.
C If JBG = 1, we call DRBGJA, and use block-grouping.
C Array RPAR, containing the current R vector, is passed to DRBDJA and
C DRBGJA as argument R0, consistent with the loading of RPAR in FWEB.
C The external name WEBR is passed, as the routine which computes the
C individual blocks of the R vector.
C-----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      EXTERNAL WEBR
      DIMENSION CC(*), CCPRIME(*), REWT(*), SAVR(*), WK(*), WP(*),
     1          IWP(*), RPAR(*), IPAR(*)
C
      JBG = IPAR(2)
      IF (JBG .EQ. 0) THEN
          CALL DRBDJA (T, CC, RPAR, WEBR, WK, REWT, CJ, WP, IWP, IER)
        ELSE
          CALL DRBGJA (T, CC, RPAR, WEBR, WK, REWT, CJ, WP, IWP, IER)
        ENDIF
C
      RETURN
C------------  End of Subroutine JACRS  --------------------------------
      END

      SUBROUTINE PSOLRS (NEQ, T, CC, CCPRIME, SAVR, WK, CJ, WT, WP, IWP,
     1                   B, EPLIN, IER, RPAR, IPAR)
C-----------------------------------------------------------------------
C This routine applies the inverse of a product preconditioner matrix
C to the vector in the array B.  Depending on the flag JPRE, this
C involves a call to GS, for the inverse of the spatial factor, and/or
C a call to DRBDPS or DRBGPS for the inverse of the reaction-based
C factor (CJ*I_d - dR/dy).  The latter factor uses block-grouping
C (with a call to DRBGPS) if JBG = 1, and does not (with a call to
C DRBDPS) if JBG = 0.  JBG is communicated as IPAR(2).
C The array B is overwritten with the solution.
C-----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION CC(*), CCPRIME(*), SAVR(*), WK(*), WP(*), IWP(*), B(*),
     1          RPAR(*), IPAR(*)
C
      JPRE = IPAR(1)
      IER = 0
      HL0 = 1.0D0/CJ
C
      JBG = IPAR(2)
C
      IF (JPRE .EQ. 2 .OR. JPRE .EQ. 3) CALL GS (NEQ, HL0, B, WK)
C
      IF (JPRE .NE. 2 ) THEN
        IF (JBG .EQ. 0) CALL DRBDPS (B, WP, IWP)
        IF (JBG .EQ. 1) CALL DRBGPS (B, WP, IWP)
        ENDIF
C
      IF (JPRE .EQ. 4) CALL GS (NEQ, HL0, B, WK)
C
      RETURN
C------------  End of Subroutine PSOLRS  -------------------------------
      END

      SUBROUTINE GS (N, HL0, Z, X)
C-----------------------------------------------------------------------
C This routine provides the inverse of the spatial factor for a
C product preconditoner in an NS-species reaction-diffusion problem.
C It performs ITMAX = 5 Gauss-Seidel iterations to compute an
C approximation to (A_S)-inverse * Z, where A_S = I - hl0*Jd, and Jd
C represents the diffusion contributions to the Jacobian.
C The solution is vector returned in Z.
C-----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION Z(*), X(*)
      PARAMETER (MAXS = 2)
      COMMON /PPAR2/ NP, NS, AX, AY, ACOEF(MAXS,MAXS), BCOEF(MAXS),
     1               DX, DY, ALPH, BETA, FPI, DIFF(MAXS),
     2               COX(MAXS), COY(MAXS), MX, MY, MXNS
      DATA ITMAX/5/
C
      DIMENSION BETA1(MAXS), GAMMA1(MAXS), BETA2(MAXS), GAMMA2(MAXS),
     1          DINV(MAXS)
C
C-----------------------------------------------------------------------
C Write matrix as A = D - L - U.
C Load local arrays BETA, BETA2, GAMMA1, GAMMA2, and DINV.
C-----------------------------------------------------------------------
      DO 10 I = 1,NS
        ELAMDA = 1.0D0/(1.0D0 + 2.0D0*HL0*(COX(I) + COY(I)))
        BETA1(I) = HL0*COX(I)*ELAMDA
        BETA2(I) = 2.0D0*BETA1(I)
        GAMMA1(I) = HL0*COY(I)*ELAMDA
        GAMMA2(I) = 2.0D0*GAMMA1(I)
        DINV(I) = ELAMDA
 10     CONTINUE
C-----------------------------------------------------------------------
C Begin iteration loop.
C Load array X with (D-inverse)*Z for first iteration.
C-----------------------------------------------------------------------
      ITER = 1
C Zero X in all its components, since X is added to Z at the end.
      DO 15 II = 1,N
 15      X(II) = 0.0d0
C
      DO 50 JY = 1,MY
        IYOFF = MXNS*(JY-1)
        DO 40 JX = 1,MX
          IC = IYOFF + NS*(JX-1)
          DO 30 I = 1,NS
            ICI = IC + I
            X(ICI) = DINV(I)*Z(ICI)
            Z(ICI) = 0.0D0
 30         CONTINUE
 40       CONTINUE
 50     CONTINUE
      GO TO 160
C-----------------------------------------------------------------------
C Calculate (D-inverse)*U*X.
C-----------------------------------------------------------------------
 70   CONTINUE
      ITER = ITER + 1
      JY = 1
        JX = 1
        IC = NS*(JX-1)
        DO 75 I = 1,NS
          ICI = IC + I
 75       X(ICI) = BETA2(I)*X(ICI+NS) + GAMMA2(I)*X(ICI+MXNS)
        DO 85 JX = 2,MX-1
          IC = NS*(JX-1)
          DO 80 I = 1,NS
            ICI = IC + I
 80         X(ICI) = BETA1(I)*X(ICI+NS) + GAMMA2(I)*X(ICI+MXNS)
 85       CONTINUE
        JX = MX
        IC = NS*(JX-1)
        DO 90 I = 1,NS
          ICI = IC + I
 90       X(ICI) = GAMMA2(I)*X(ICI+MXNS)
      DO 115 JY = 2,MY-1
        IYOFF = MXNS*(JY-1)
          JX = 1
          IC = IYOFF
          DO 95 I = 1,NS
            ICI = IC + I
 95         X(ICI) = BETA2(I)*X(ICI+NS) + GAMMA1(I)*X(ICI+MXNS)
          DO 105 JX = 2,MX-1
            IC = IYOFF + NS*(JX-1)
            DO 100 I = 1,NS
              ICI = IC + I
 100          X(ICI) = BETA1(I)*X(ICI+NS) + GAMMA1(I)*X(ICI+MXNS)
 105        CONTINUE
          JX = MX
          IC = IYOFF + NS*(JX-1)
          DO 110 I = 1,NS
            ICI = IC + I
 110        X(ICI) = GAMMA1(I)*X(ICI+MXNS)
 115      CONTINUE
      JY = MY
      IYOFF = MXNS*(JY-1)
        JX = 1
        IC = IYOFF
        DO 120 I = 1,NS
          ICI = IC + I
 120      X(ICI) = BETA2(I)*X(ICI+NS)
        DO 130 JX = 2,MX-1
          IC = IYOFF + NS*(JX-1)
          DO 125 I = 1,NS
            ICI = IC + I
 125      X(ICI) = BETA1(I)*X(ICI+NS)
 130      CONTINUE
        JX = MX
        IC = IYOFF + NS*(JX-1)
        DO 135 I = 1,NS
          ICI = IC + I
 135      X(ICI) = 0.0D0
C-----------------------------------------------------------------------
C Calculate (I - (D-inverse)*L)*X.
C-----------------------------------------------------------------------
 160  CONTINUE
      JY = 1
        DO 175 JX = 2,MX-1
          IC = NS*(JX-1)
          DO 170 I = 1,NS
            ICI = IC + I
 170        X(ICI) = X(ICI) + BETA1(I)*X(ICI-NS)
 175      CONTINUE
        JX = MX
        IC = NS*(JX-1)
        DO 180 I = 1,NS
          ICI = IC + I
 180      X(ICI) = X(ICI) + BETA2(I)*X(ICI-NS)
      DO 210 JY = 2,MY-1
        IYOFF = MXNS*(JY-1)
          JX = 1
          IC = IYOFF
          DO 185 I = 1,NS
            ICI = IC + I
 185        X(ICI) = X(ICI) + GAMMA1(I)*X(ICI-MXNS)
          DO 200 JX = 2,MX-1
            IC = IYOFF + NS*(JX-1)
            DO 195 I = 1,NS
              ICI = IC + I
              X(ICI) = (X(ICI) + BETA1(I)*X(ICI-NS))
     1             + GAMMA1(I)*X(ICI-MXNS)
 195          CONTINUE
 200        CONTINUE
            JX = MX
            IC = IYOFF + NS*(JX-1)
            DO 205 I = 1,NS
              ICI = IC + I
              X(ICI) = (X(ICI) + BETA2(I)*X(ICI-NS))
     1             + GAMMA1(I)*X(ICI-MXNS)
 205          CONTINUE
 210        CONTINUE
      JY = MY
      IYOFF = MXNS*(JY-1)
        JX = 1
        IC = IYOFF
        DO 215 I = 1,NS
          ICI = IC + I
 215      X(ICI) = X(ICI) + GAMMA2(I)*X(ICI-MXNS)
        DO 225 JX = 2,MX-1
          IC = IYOFF + NS*(JX-1)
          DO 220 I = 1,NS
            ICI = IC + I
            X(ICI) = (X(ICI) + BETA1(I)*X(ICI-NS))
     1           + GAMMA2(I)*X(ICI-MXNS)
 220        CONTINUE
 225      CONTINUE
        JX = MX
        IC = IYOFF + NS*(JX-1)
        DO 230 I = 1,NS
          ICI = IC + I
          X(ICI) = (X(ICI) + BETA2(I)*X(ICI-NS))
     1         + GAMMA2(I)*X(ICI-MXNS)
 230      CONTINUE
C-----------------------------------------------------------------------
C Add increment X to Z.
C-----------------------------------------------------------------------
      DO 300 I = 1,N
 300    Z(I) = Z(I) + X(I)
C
      IF (ITER .LT. ITMAX) GO TO 70
      RETURN
C------------  End of Subroutine GS  -----------------------------------
      END

max_stars_repo_name: faribas/PyDAS

max_issues_repo_name: faribas/PyDAS

max_forks_repo_name: faribas/PyDAS


text:
        SUBROUTINE YCLMCR (RESULT, REFRAT, RETARD, NX, MAXLEN, NUPDAT,
     *                     ISTAT)
*
*  Module number:
*
*  Module name: YCLMCR
*
*  Keyphrase:
*  ----------
*       Polarimetry processing
*
*  Description:
*  ------------
*       Inverts mirror reflection
*
*  FORTRAN name: YCLMCR
*
*  Keywords of accessed files and tables:
*  --------------------------------------
*       none
*
*  Subroutines Called:
*  -------------------
*  CDBS:
*       yclcrt, yclcpl, yclcst, ymsput
*
*  SDAS:
*
*
*  History:
*  --------
*  Version      Date    Author		  Description
*       1       Mar97   M.D. De La Pena   Original Implementation for CALFOS
*                                         based upon code from R. Allen.
*       1.1     Jul97   M.D. De La Pena   Mods for new parameter
*-------------------------------------------------------------------------------
*
* INPUTS:
*       RESULT - Array containing Stokes parameters, linear and circular
*                polarization, polarization angle, and error estimates
*       REFRAT - Array of the reflection ratio of the mirrors, rp/rs
*       RETARD - Array of the phase of a retarder
*       NX     - Number of channels
*       MAXLEN - Maximum number of channels
*       NUPDAT - Number of quantities being updated
*       ISTAT  - Return status
*
* OUTPUT:
*       RESULT - Modified input array
*       ISTAT  - Modified return status
*----------------------------------------------------------------------------
C
C UMSPUT DESTINATIONS -- CB, DAO, 4-SEP-87
C
      INTEGER STDOUT
      PARAMETER (STDOUT = 1)
      INTEGER STDERR
      PARAMETER (STDERR = 2)
C
C Passed Parameters
C
      INTEGER  NX, MAXLEN, NUPDAT, ISTAT
      REAL     RESULT(MAXLEN,NUPDAT)
      DOUBLE PRECISION RETARD(MAXLEN), REFRAT(MAXLEN)
C
C Local Variables
C
      CHARACTER*80 CONTXT
      INTEGER      CHAN, I, J, K
      DOUBLE PRECISION CHIRAD, DELTRAD, PSIRAD
      DOUBLE PRECISION SUM, RP, RS, DPI
      DOUBLE PRECISION EQUIV(4,4), RETMAT(4,4), NLPMAT(4,4)
C
C CHIRAD is the angle of the plane of incidence = 0.0 in the COSTAR frame
C PSIRAD is the polarizing axis which is perpendicular to the plane of
C   incidence == 90 degrees in the COSTAR frame
C DELTRAD is the retardance
C RS is the reflectivity perpendicular to the plane of incidence
C RP is the reflectivity parallel to the plane of incidence
C
      PARAMETER (CHIRAD = 0.0D0)
      PARAMETER (RS     = 1.0D0)
      PARAMETER (DPI    = 3.1415926536D0)
C
C----------------------------------------------------------------------------
C
C Loop over all the channels, DELTRAD and RP are a function of channel
C
      PSIRAD = DPI / 2.0D0
C
      DO 40 CHAN = 1, NX
          DELTRAD = -RETARD(CHAN)
          RP      = 1.0D0 / REFRAT(CHAN)
C
C Generate the Mueller Matrix for the retarder
C
          CALL YCLCRT (CHIRAD, DELTRAD, RETMAT, ISTAT)
          IF (ISTAT .NE. 0) THEN
              CONTXT = 'ERROR in computing the Mueller Matrix for' //
     *                 ' the retarder.'
              GO TO 999
          END IF
C
C Generate the Mueller Matrix for the polarizer
C
          CALL YCLCPL (RS, RP, PSIRAD, NLPMAT, ISTAT)
          IF (ISTAT .NE. 0) THEN
              CONTXT = 'ERROR in computing the Mueller Matrix for' //
     *                 ' the polarizer.'
              GO TO 999
          END IF
C
          DO 30 I = 1, 4
              DO 20 J =1, 4
                  SUM = 0.0D0
                  DO 10 K = 1, 4
                      SUM = SUM + NLPMAT(I,K) * RETMAT(K,J)
   10             CONTINUE
                  EQUIV(I,J) = SUM
   20         CONTINUE
   30     CONTINUE
C
          CALL YCLCST(RESULT, MAXLEN, NUPDAT, CHAN, EQUIV, ISTAT)
          IF (ISTAT .NE. 0) THEN
              CONTXT = 'ERROR in computing the normalized Stokes' //
     *                 ' parameters.'
              GO TO 999
          END IF
C
   40 CONTINUE
C
      ISTAT = 0
      GO TO 1000
999   CALL YMSPUT(CONTXT,STDOUT+STDERR,0,ISTAT)
      ISTAT = 1
1000  RETURN
      END

max_stars_repo_name: iraf-community/stsdas


text:
c do continuum subtraction by subtracting a low order Legendre polynomial
c fit to the good data points.
c I added a sigma-clipping routine before the fit, which may not
c work great since the continuum varies, avsigclip or something
c would be better.

c do the subtraction in place, or not
      subroutine contsubleg(wave,spec,serr,nsp,iord)

c      parameter (IORDER=15)

      real wave(nsp),spec(nsp),serr(nsp)
      real wscale(nsp),specclip(nsp)
      real wfit(nsp),ffit(nsp),efit(nsp)
      real a(iord),uu(nsp,iord),vv(iord,iord),ww(iord)
      real avals(iord)

      external fleg

      include 'pcredshift.h'

      ncliprad=100
      clipsig=2.0
      call localsigclip(nsp,spec,specclip,ncliprad,clipsig)

      wrange = wave(nsp) - wave(1)
      do i=1,nsp
c rescale wfit to run between 0 and 1 for improved polynomial fitting
         wscale(i) = (wave(i) - wave(1)) / wrange
      end do

      npoly=iord
      nfit=0
      do i=1,nsp
         if(specclip(i) .gt. bad .and. serr(i) .gt. bad
     $        .and. serr(i) .lt. 1.0e6) then
            nfit=nfit+1
            wfit(nfit)=wscale(i)
            ffit(nfit)=specclip(i)
            efit(nfit)=serr(i)
         end if
      end do
      if (nfit .le. npoly) return

c do fit
      call bigsvdfit(wfit,ffit,efit,nfit,a,npoly,uu,vv,ww,nsp,iord,
     $     chisq,fleg)

c evaluate and subtract fit
      do i=1,nsp
c         call fleg(wave(i),avals,npoly)
         call fleg(wscale(i),avals,npoly)
         fitval = 0.
         do j=1,npoly
            fitval = fitval + a(j)*avals(j)
         end do
         if (spec(i) .gt. bad) then
            spec(i) = spec(i) - fitval
         end if
      end do

      return
      end

cccccccccccccccccccccccccccccc
c from Numerical Recipes, generates Legendre polynomials
c via recurrence relation

      SUBROUTINE fleg(x,pl,nl)
      INTEGER nl
      REAL x,pl(nl)
      INTEGER j
      REAL d,f1,f2,twox
      pl(1)=1.
      pl(2)=x
      if(nl.gt.2) then
        twox=2.*x
        f2=x
        d=1.
        do 11 j=3,nl
          f1=d
          f2=f2+twox
          d=d+1.
          pl(j)=(f2*pl(j-1)-f1*pl(j-2))/d
11      continue
      endif
      return
      END