dsrc2c.f90

河口模型 使用模拟盐水入侵、热量扩散等等 河口模型 使用模拟盐水入侵、热量扩散

!       
   90 IB1 = 1     
      IB2 = IB1+N 
      IB3 = IB2+N 
      IB4 = IB3+N 
      IB5 = IB4+N 
      IB6 = IB5+N 
      IB7 = IB6+N 
      IPARM(8) = 6*N+2*ITMAX
      IF (ISYM.NE.0) IPARM(8) = IPARM(8)+2*ITMAX
      IF (NW.GE.IPARM(8)) GO TO 110   
      IER = 42    
      IF (LEVEL.GE.0) WRITE (NOUT,100) NW,IPARM(8)
  100 FORMAT ('0','*** F A T A L     E R R O R ************'/'0',   &
     &   '    CALLED FROM ITPACK ROUTINE SSORCG '/' ',    &
     &   '    NOT ENOUGH WORKSPACE AT ',I10/' ','    SET IPARM(8) =',I10&
     &   ,' (NW)')
      GO TO 390   
!       
! ... PERMUTE TO  RED-BLACK SYSTEM IF REQUESTED 
!       
  110 NB = IPARM(9) 
      IF (NB.LT.0) GO TO 170
      N3 = 3*N    
      CALL IVFILL (N3,IWKSP,0)
      CALL PRBNDX (N,NB,IA,JA,IWKSP,IWKSP(IB2),LEVEL,NOUT,IER)      
      IF (IER.EQ.0) GO TO 130 
      IF (LEVEL.GE.0) WRITE (NOUT,120) IER,NB   
  120 FORMAT ('0','*** F A T A L     E R R O R ************'/'0',   &
     &   '    CALLED FROM ITPACK ROUTINE SSORCG '/' ',    &
     &   '    ERROR DETECTED IN SUBROUTINE  PRBNDX'/' ',  &
     &   '    WHICH COMPUTES THE RED-BLACK INDEXING'/' ','    IER = ',I5&
     &   ,' IPARM(9) = ',I5,' (NB)')  
      GO TO 390   
!       
! ... PERMUTE MATRIX AND RHS
!       
  130 IF (LEVEL.GE.2) WRITE (NOUT,140) NB       
  140 FORMAT (/10X,'ORDER OF BLACK SUBSYSTEM = ',I5,' (NB)')
      CALL PERMAT (N,IA,JA,A,IWKSP,IWKSP(IB3),ISYM,LEVEL,NOUT,IER)  
      IF (IER.EQ.0) GO TO 160 
      IF (LEVEL.GE.0) WRITE (NOUT,150) IER      
  150 FORMAT ('0','*** F A T A L     E R R O R ************'/'0',   &
     &   '    CALLED FROM ITPACK ROUTINE SSORCG '/' ',    &
     &   '    ERROR DETECTED IN SUBROUTINE  PERMAT'/' ',  &
     &   '    WHICH DOES THE RED-BLACK PERMUTATION'/' ','    IER = ',I5)
      GO TO 390   
  160 CALL PERVEC (N,RHS,IWKSP)       
      CALL PERVEC (N,U,IWKSP) 
!       
! ... SCALE LINEAR SYSTEM, U, AND RHS BY THE SQUARE ROOT OF THE     
! ... DIAGONAL ELEMENTS.    
!       
  170 CONTINUE    
      CALL VFILL (IPARM(8),WKSP,0.0D0)
      CALL SCAL (N,IA,JA,A,RHS,U,WKSP,LEVEL,NOUT,IER)     
      IF (IER.EQ.0) GO TO 190 
      IF (LEVEL.GE.0) WRITE (NOUT,180) IER      
  180 FORMAT ('0','*** F A T A L     E R R O R ************'/'0',   &
     &   '    CALLED FROM ITPACK ROUTINE SSORCG '/' ',    &
     &   '    ERROR DETECTED IN SUBROUTINE  SCAL  '/' ',  &
     &   '    WHICH SCALES THE SYSTEM   '/' ','    IER = ',I5)      
      GO TO 390   
  190 IF (LEVEL.LE.2) GO TO 220       
      WRITE (NOUT,200)      
  200 FORMAT (///1X,'IN THE FOLLOWING, RHO AND GAMMA ARE',&
     &   ' ACCELERATION PARAMETERS')  
      WRITE (NOUT,210)      
  210 FORMAT (1X,'S-PRIME IS AN INITIAL ESTIMATE FOR NEW CME')      
  220 IF (IPARM(11).NE.0) GO TO 230   
      TIMI1 = TIMER(DUMMY)  
!       
! ... SPECIAL PROCEDURE FOR FULLY ADAPTIVE CASE.
!       
  230 CONTINUE    
      IF (.NOT.ADAPT) GO TO 250       
      IF (.NOT.BETADT) GO TO 240      
      CALL VFILL (N,WKSP(IB1),1.D0)   
      BETNEW = PBETA(N,IA,JA,A,WKSP(IB1),WKSP(IB2),WKSP(IB3))/      &
     &   DBLE(FLOAT(N))     
      BETAB = DMAX1(BETAB,.25D0,BETNEW) 
  240 CALL OMEG (0.D0,1)    
      IS = 0      
!       
! ... INITIALIZE FORWARD PSEUDO-RESIDUAL
!       
  250 CALL DCOPY (N,RHS,1,WKSP(IB1),1)
      CALL DCOPY (N,U,1,WKSP(IB2),1)  
      CALL PFSOR (N,IA,JA,A,WKSP(IB2),WKSP(IB1))
      CALL VEVMW (N,WKSP(IB2),U)      
!       
! ... ITERATION SEQUENCE    
!       
      ITMAX1 = ITMAX+1      
      DO 270 LOOP = 1,ITMAX1
         IN = LOOP-1
         IF (MOD(IN,2).EQ.1) GO TO 260
!       
! ... CODE FOR THE EVEN ITERATIONS.   
!       
!     U           = U(IN)       WKSP(IB2) = C(IN) 
!     WKSP(IB1)   = U(IN-1)     WKSP(IB3) = C(IN-1)       
!       
         CALL ITSRCG (N,IA,JA,A,RHS,U,WKSP(IB1),WKSP(IB2),WKSP(IB3),&
     &      WKSP(IB4),WKSP(IB5),WKSP(IB6),WKSP(IB7))      
!       
         IF (HALT) GO TO 300
         GO TO 270
!       
! ... CODE FOR THE ODD ITERATIONS.    
!       
!     U           = U(IN-1)     WKSP(IB2) = C(IN-1)       
!     WKSP(IB1)   = U(IN)       WKSP(IB3) =C(IN)
!       
  260    CALL ITSRCG (N,IA,JA,A,RHS,WKSP(IB1),U,WKSP(IB3),WKSP(IB2),&
     &      WKSP(IB4),WKSP(IB5),WKSP(IB6),WKSP(IB7))      
!       
         IF (HALT) GO TO 300
  270 CONTINUE    
!       
! ... ITMAX HAS BEEN REACHED
!       
      IF (IPARM(11).NE.0) GO TO 280   
      TIMI2 = TIMER(DUMMY)  
      TIME1 = DBLE(TIMI2-TIMI1)       
  280 IF (LEVEL.GE.1) WRITE (NOUT,290) ITMAX    
  290 FORMAT ('0','*** W A R N I N G ************'/'0',   &
     &   '    IN ITPACK ROUTINE SSORCG'/' ','    FAILURE TO CONVERGE IN'&
     &   ,I5,' ITERATIONS') 
      IER = 43    
      IF (IPARM(3).EQ.0) RPARM(1) = STPTST      
      GO TO 330   
!       
! ... METHOD HAS CONVERGED  
!       
  300 IF (IPARM(11).NE.0) GO TO 310   
      TIMI2 = TIMER(DUMMY)  
      TIME1 = DBLE(TIMI2-TIMI1)       
  310 IF (LEVEL.GE.1) WRITE (NOUT,320) IN       
  320 FORMAT (/1X,'SSORCG  HAS CONVERGED IN ',I5,' ITERATIONS')     
!       
! ... PUT SOLUTION INTO U IF NOT ALREADY THERE. 
!       
  330 CONTINUE    
      IF (MOD(IN,2).EQ.1) CALL DCOPY (N,WKSP(IB1),1,U,1)  
!       
! ... UNSCALE THE MATRIX, SOLUTION, AND RHS VECTORS.      
!       
      CALL UNSCAL (N,IA,JA,A,RHS,U,WKSP)
!       
! ... UN-PERMUTE MATRIX,RHS, AND SOLUTION       
!       
      IF (IPARM(9).LT.0) GO TO 360    
      CALL PERMAT (N,IA,JA,A,IWKSP(IB2),IWKSP(IB3),ISYM,LEVEL,NOUT, &
     &   IERPER)  
      IF (IERPER.EQ.0) GO TO 350      
      IF (LEVEL.GE.0) WRITE (NOUT,340) IERPER   
  340 FORMAT ('0','*** F A T A L     E R R O R ************'/'0',   &
     &   '    CALLED FROM ITPACK ROUTINE SSORCG '/' ',    &
     &   '    ERROR DETECTED IN SUBROUTINE  PERMAT'/' ',  &
     &   '    WHICH UNDOES THE RED-BLACK PERMUTATION   '/' ',       &
     &   '    IER = ',I5)   
      IF (IER.EQ.0) IER = IERPER      
      GO TO 390   
  350 CALL PERVEC (N,RHS,IWKSP(IB2))  
      CALL PERVEC (N,U,IWKSP(IB2))    
!       
! ... OPTIONAL ERROR ANALYSIS 
!       
  360 IDGTS = IPARM(12)     
      IF (IDGTS.LT.0) GO TO 370       
      IF (IPARM(2).LE.0) IDGTS = 0    
      CALL PERROR5 (N,IA,JA,A,RHS,U,WKSP,DIGIT1,DIGIT2,IDGTS)
!       
! ... SET RETURN PARAMETERS IN IPARM AND RPARM  
!       
  370 IF (IPARM(11).NE.0) GO TO 380   
      TIMJ2 = TIMER(DUMMY)  
      TIME2 = DBLE(TIMJ2-TIMJ1)       
  380 IPARM(8) = IPARM(8)-2*(ITMAX-IN)
      IF (ISYM.NE.0) IPARM(8) = IPARM(8)-2*(ITMAX-IN)     
      IF (IPARM(3).NE.0) GO TO 390    
      IPARM(1) = IN 
      IPARM(9) = NB 
      RPARM(2) = CME
      RPARM(3) = SME
      RPARM(5) = OMEGA      
      RPARM(6) = SPECR      
      RPARM(7) = BETAB      
      RPARM(9) = TIME1      
      RPARM(10) = TIME2     
      RPARM(11) = DIGIT1    
      RPARM(12) = DIGIT2    
!       
  390 CONTINUE    
      IERR = IER  
      IF (LEVEL.GE.3) CALL ECHALL (N,IA,JA,A,RHS,IPARM,RPARM,2)     
!       
      RETURN      
      END 
      SUBROUTINE SSORSI (NN,IA,JA,A,RHS,U,IWKSP,NW,WKSP,IPARM,RPARM,&
     &   IERR)    
!       
!     ITPACK 2C MAIN SUBROUTINE  SSORSI  (SYMMETRIC SUCCESSIVE RELAX- 
!                                         ATION SEMI-ITERATION)     
!     EACH OF THE MAIN SUBROUTINES:   
!           JCG, JSI, SOR, SSORCG, SSORSI, RSCG, RSSI     
!     CAN BE USED INDEPENDENTLY OF THE OTHERS   
!       
! ... FUNCTION:   
!       
!          THIS SUBROUTINE, SSORSI, DRIVES THE  SYMMETRIC SOR-SI    
!          ALGORITHM.       
!       
! ... PARAMETER LIST:       
!       
!          N      INPUT INTEGER.  DIMENSION OF THE MATRIX. (= NN)   
!          IA,JA  INPUT INTEGER VECTORS.  THE TWO INTEGER ARRAYS OF 
!                 THE SPARSE MATRIX REPRESENTATION.       
!          A      INPUT D.P. VECTOR.  THE D.P. ARRAY OF THE SPARSE  
!                 MATRIX REPRESENTATION.
!          RHS    INPUT D.P. VECTOR.  CONTAINS THE RIGHT HAND SIDE  
!                 OF THE MATRIX PROBLEM.
!          U      INPUT/OUTPUT D.P. VECTOR.  ON INPUT, U CONTAINS THE 
!                 INITIAL GUESS TO THE SOLUTION. ON OUTPUT, IT CONTAINS 
!                 THE LATEST ESTIMATE TO THE SOLUTION.    
!          IWKSP  INTEGER VECTOR WORKSPACE OF LENGTH 3*N  
!          NW     INPUT INTEGER.  LENGTH OF AVAILABLE WKSP.  ON OUTPUT, 
!                 IPARM(8) IS AMOUNT USED.      
!          WKSP   D.P. VECTOR USED FOR WORKING SPACE.  SSORSI       
!                 NEEDS THIS TO BE IN LENGTH AT LEAST  5*N
!          IPARM  INTEGER VECTOR OF LENGTH 12.  ALLOWS USER TO SPECIFY
!                 SOME INTEGER PARAMETERS WHICH AFFECT THE METHOD.  IF
!          RPARM  D.P. VECTOR OF LENGTH 12. ALLOWS USER TO SPECIFY SOME 
!                 D.P. PARAMETERS WHICH AFFECT THE METHOD.
!          IER    OUTPUT INTEGER.  ERROR FLAG. (= IERR)   
!       
! ... SSORSI SUBPROGRAM REFERENCES:   
!       
!          FROM ITPACK    BISRCH, CHEBY, CHGSI, DFAULT, ECHALL, ECHOUT, 
!                         ITERM, TIMER, ITSRSI, IVFILL, OMEG,       
!                         OMGSTR, PARSI, PBETA, PERMAT, PERROR5,     
!                         PERVEC, PFSOR, PMULT, PRBNDX, PSSOR1,     
!                         PSTOP, PVTBV, QSORT, SBELM, SCAL, DCOPY,  
!                         DDOT, SUM3, TSTCHG, UNSCAL, VEVPW, VFILL, 
!                         VOUT, WEVMW 
!          SYSTEM         DABS, DLOG, DLOG10, DBLE(AMAX0), DMAX1, DBLE(F
!                         MOD, DSQRT  
!       
!     VERSION:  ITPACK 2C (MARCH 1982)
!       
!     CODE WRITTEN BY:  DAVID KINCAID, ROGER GRIMES, JOHN RESPESS   
!                       CENTER FOR NUMERICAL ANALYSIS     
!                       UNIVERSITY OF TEXAS     
!                       AUSTIN, TX  78712       
!                       [phone]
!       
!     FOR ADDITIONAL DETAILS ON THE   
!          (A) SUBROUTINE SEE TOMS ARTICLE 1982 
!          (B) ALGORITHM  SEE CNA REPORT 150    
!       
!     BASED ON THEORY BY:  DAVID YOUNG, DAVID KINCAID, LOU HAGEMAN  
!       
!     REFERENCE THE BOOK:  APPLIED ITERATIVE METHODS      
!                          L. HAGEMAN, D. YOUNG 
!                          ACADEMIC PRESS, 1981 
!       
!     **************************************************  
!     *               IMPORTANT NOTE                   *  
!     *                                                *  
!     *      WHEN INSTALLING ITPACK ROUTINES ON A      *  
!     *  DIFFERENT COMPUTER, RESET SOME OF THE VALUES  *  
!     *  IN  SUBROUTNE DFAULT.   MOST IMPORTANT ARE    *  
!     *                                                *  
!     *   DRELPR      MACHINE RELATIVE PRECISION       *  
!     *   RPARM(1)    STOPPING CRITERION               *  
!     *                                                *  
!     *   ALSO CHANGE SYSTEM-DEPENDENT ROUTINE         *  
!     *   SECOND USED IN TIMER                         *  
!     *                                                *  
!     **************************************************  
!       
!     SPECIFICATIONS FOR ARGUMENTS    
!       
      INTEGER IA(1),JA(1),IWKSP(1),IPARM(12),NN,NW,IERR   
      DOUBLE PRECISION A(1),RHS(NN),U(NN),WKSP(NW),RPARM(12)
!       
!     SPECIFICATIONS FOR LOCAL VARIABLES
!       
      INTEGER IB1,IB2,IB3,IB4,IB5,IDGTS,IER,IERPER,ITMAX1,LOOP,N,NB,N3
      DOUBLE PRECISION BETNEW,DIGIT1,DIGIT2,PBETA,TEMP,TIME1,TIME2,TOL
!       
! *** BEGIN: ITPACK COMMON  
!       
      INTEGER IN,IS,ISYM,ITMAX,LEVEL,NOUT       
      COMMON /ITCOM1/ IN,IS,ISYM,ITMAX,LEVEL,NOUT 
!       
      LOGICAL ADAPT,BETADT,CASEII,HALT,PARTAD   
      COMMON /ITCOM2/ ADAPT,BETADT,CASEII,HALT,PARTAD     
!       
      DOUBLE PRECISION BDELNM,BETAB,CME,DELNNM,DELSNM,FF,GAMMA,OMEGA,QA,&
     &   QT,RHO,RRR,SIGE,SME,SPECR,SPR,DRELPR,STPTST,UDNM,ZETA      
      COMMON /ITCOM3/ BDELNM,BETAB,CME,DELNNM,DELSNM,FF,GAMMA,OMEGA,QA, &
     &   QT,RHO,RRR,SIGE,SME,SPECR,SPR,DRELPR,STPTST,UDNM,ZETA      
!       
! *** END  : ITPACK COMMON  
!       
! ... VARIABLES IN COMMON BLOCK - ITCOM1
!       
!     IN     - ITERATION NUMBER       
!     ISYM   - SYMMETRIC/NONSYMMETRIC STORAGE FORMAT SWITCH 
!     IS     - ITERATION NUMBER WHEN PARAMETERS LAST CHANGED
!     ITMAX  - MAXIMUM NUMBER OF ITERATIONS ALLOWED       
!     LEVEL  - LEVEL OF OUTPUT CONTROL SWITCH   
!     NOUT   - OUTPUT UNIT NUMBER     
!       
! ... VARIABLES IN COMMON BLOCK - ITCOM2
!       
!     ADAPT  - FULLY ADAPTIVE PROCEDURE SWITCH  
!     BETADT - SWITCH FOR ADAPTIVE DETERMINATION OF BETA  
!     CASEII - ADAPTIVE PROCEDURE CASE SWITCH   
!     HALT   - STOPPING TEST SWITCH   
!     PARTAD - PARTIALLY ADAPTIVE PROCEDURE SWITCH
!       
! ... VARIABLES IN COMMON BLOCK - ITCOM3
!       
!     BDELNM - TWO NORM OF B TIMES DELTA-SUPER-N
!     BETAB  - ESTIMATE FOR THE SPECTRAL RADIUS OF LU MATRIX
!     CME    - ESTIMATE OF LARGEST EIGENVALUE   
!     DELNNM - INNER PRODUCT OF PSEU