CC======================================================================
CC
CC *** DMT CCN COUNTER SIMULATOR, version 13 JUN 2004
CC     MINIMUM I/O VERSION FOR "PIPE" version of model.
CC
CC     THE HYDRODYNAMIC & OTHER VARIABLES (U,V,P,C,T) ARE CALCULATED
CC     USING A FINITE VOLUME FORMULATION ON A STAGGERED C-GRID
CC 
CC *** COPYRIGHT 2004 ATHANASIOS NENES
CC
CC=======================================================================
CC
C=======================================================================
C
C *** SUBROUTINE INPUT
C *** THIS SUBROUTINE READS THE INPUT FILE FOR DEFAULT PROBLEM SETTINGS
C
C=======================================================================
C
      SUBROUTINE INPUT (PROBFIL, EOF)
      INCLUDE 'common_l.inc'
      INCLUDE 'prob.inc'
      INCLUDE 'aerosol.inc'
      CHARACTER PROBFIL*(*)
      LOGICAL EOF
C
C *** READ INPUT FROM STDIN ********************************************
C
100   READ (*,*,ERR=100,END=999) VFL, RAT, PRESF, TINS, TCOLD, THOT,
     &                           OUTFIL
C
C *** OPEN INPUT FILE **************************************************
C
      OPEN(UNIT=IOINP,FILE=PROBFIL,STATUS='OLD')
C
C *** READ HEADER OF INPUT FILE ****************************************
C
      READ(IOINP,'(//)')
C
C *** GEOMETRICAL PROBLEM PARAMETERS ***********************************
C
      READ(IOINP,'(/)')
      READ(IOINP,*) UPARAB
      READ(IOINP,*)
      READ(IOINP,*) RHINS, RHINA
C
      TINA   = TINS
      TUBLEN = 0.50                ! Tube length (fixed)
      RADIUS = 1.15e-2             ! Radius (fixed)
      RADINJ = 2.286e-3            ! Aerosol injector radius (fixed)
      IF (UPARAB) THEN
         AREA   = PI*RADIUS**2.0          ! Total section
         UMN    = VFL/AREA                ! Mean velocity
         VFA    = VFL/(1.0+RAT)           ! Flow rate inner
         TERM   = RADIUS**4.0 - RADIUS*RADIUS*VFA/PI/UMN
         RADINJ = SQRT(RADIUS*RADIUS - SQRT(TERM)) ! Radius of inject
         UINA   = UMN
         UINS   = UINA
      ELSE
         ARA    = PI*RADINJ**2.0          ! Inner section
         ARS    = PI*RADIUS**2.0 - ARA    ! Outer section
         VFA    = VFL/(1.0+RAT)           ! Flow rate inner
         VFS    = VFL - VFA               ! Flow rate out
         UINA   = VFA/ARA                 ! Inner velocity
         UINS   = VFS/ARS                 ! Outer velocity
      ENDIF
C
C *** AEROSOL PARAMETERS ***********************************************
C
      DPAVG = 0.12e-6
      SIGV  = 1.35
      ANTOT = 0.0d6
      NSEC  = 3
      PDNS  = 1760d0
      PMW   = 132d-3
      PVHF  = 3.0
C
C *** GRID PARAMETERS **************************************************
C
      READ(IOINP,'(/)')
      POLAR = .TRUE.
      READ(IOINP,*) NI, NJ
      NI = MIN(II,NI)      ! MAKE SURE NI,NJ ARE ACCEPTABLE
      NJ = MIN(JJ,NJ)
      FEXPX = 1.0
      FEXPY = 1.0
C
C ** CREATE X-GRID POINTS 
C
      CALL GRID1 (X,II,0.0,TUBLEN,2,NI-1)
      X(1)=-X(2)
      X(NI)=2.0*TUBLEN-X(NI-1)
C
C ** CREATE Y-GRID POINTS 
C
      NJAER = MAX (NJ*INT(RADINJ/RADIUS),2)
      CALL GRID1 (Y,JJ,0.0,RADINJ,2,NJAER)
      CALL GRID1 (Y,JJ,RADINJ,RADIUS-RADINJ,NJAER+1,NJ-1)
      Y(1)=-Y(2)
      Y(NJ)=2.0*RADIUS-Y(NJ-1)
C
C *** TIME PARAMETERS **************************************************
C
      STEADY = .TRUE.
      NT     = 10
      DT     = 0.5
C
C *** DEPENDENT VARIABLES **********************************************
C
      READ(IOINP,'(/)')
      READ(IOINP,*) INCALU
      INCALV=INCALU
      INCALP=INCALU
      READ(IOINP,*)
      READ(IOINP,*) INCALT, INCALC
C
C *** THERMO-PHYSICAL PROPERTIES ***************************************
C
      INPRO=.TRUE.
C
C *** TERMINATION CRITERIA FOR ITERATIONS ******************************
C
      READ(IOINP,'(/)')
      READ(IOINP,*) MAXIT
      READ(IOINP,*)
      READ(IOINP,*) SORMAX
      SNORM = 1.0
      SNORU = SNORM
      SNORV = SNORM
      SNORT = SNORM
      SNORC = SNORM
C
C *** UNDER RELAXATION *************************************************
C
      READ(IOINP,'(/)')
      READ(IOINP,*) URFU, URFP
      URFV = URFU
      READ(IOINP,*)
      READ(IOINP,*) URFT, URFC
C
C *** COORDINATES OF FIXED PRESSURE POINT ******************************
C
      IPREF = 2
      JPREF = 2
C
C *** EQUATION SOLVER SELECTION ****************************************
C
      READ(IOINP, '(/)')
      READ(IOINP,*) ISLVU, ISLVP
      ISLVV = ISLVU
      READ(IOINP,*)
      READ(IOINP,*) ISLVT, ISLVC
      ALPHA = 0.8
C
C *** TERMINATION CRITERIA FOR EQUATION SOLVER *************************
C
      LITERU = 3
      LITERV = 3
      LITERP = 5
      LITERT = 10 ! 10
      LITERC = 3 !  5
      SORU   = 1.0e-6
      SORV   = SORU
      SORP   = SORU
      SORT   = SORU
      SORC   = SORU
C
C *** INPUT-OUTPUT PARAMETERS ******************************************
C
      READ(IOINP,'(/)')
      READ(IOINP,*) SOAVGP
      IF (OUTFIL(1:3).eq.'nul' .OR. OUTFIL(1:3).eq.'NUL') THEN
         OUTPF = .FALSE.
      ELSE
         OUTPF = .TRUE.
      ENDIF
      IMON   = 5
      JMON   = 5
      NITPRI = 1500
      NSTPRI = 1
      INREAD =.FALSE.
      INWRIT =.FALSE.
      INPLOT = INCALU
      INPRES = INCALU
      INTEMP = INCALT
      INCONC = INCALC
      RSTFIL = OUTFIL
      CALL APPENDEXT (RSTFIL,'.sav',.TRUE.)
      SAVFIL = RSTFIL
C
C *** CLOSE INPUT FILE *************************************************
C
      CLOSE(UNIT=IOINP)
      RETURN
C
C *** END OF CONSOLE INPUT - TERMINATE RUNS ****************************
C
999   EOF = .TRUE.
      RETURN
      END
CC
C=======================================================================
C
C *** SUBROUTINE INIT
C *** THIS SUBROUTINE IS USED TO INITIALIZE FIELDS OF VARIABLES,
C     PHYSICAL PROPERTIES, ETC.
C
C=======================================================================
C
      SUBROUTINE INIT
      INCLUDE 'common_l.inc'
      INCLUDE 'prob.inc'
      INCLUDE 'aerosol.inc'
C
C *** WALL TEMPERATURE & WATER VAPOR CONCENTRATION *********************
C
      DTDX = ABS(THOT-TCOLD)/TUBLEN   ! Add ABS for error-proofing.
      DO 200 I=2,NI-1
         T(I,NJ) = TCOLD+DTDX*X(I)  ! temperature of wall
         C(I,NJ) = VPRES(T(I,NJ))*1e2/RGAS/T(I,NJ) ! C of water vapor
 200  CONTINUE
      T(NI,NJ) = T(NI-1,NJ)
      T(1, NJ) = T(2,   NJ)
      C(NI,NJ) = C(NI-1,NJ)
      C(1, NJ) = C(2,   NJ)
C
C *** VELOCITY AND TEMPERATURE PROFILES ********************************
C
      DO 100 I=1,NI
        DO 110 J=1,NJM1
          IF (J.LE.NJAER) THEN
             U(I,J)    = UINA
             IF (UPARAB) U(I,J)=2*UINA*(1-ABS(Y(J)/RADIUS)**2)
             T(I,J)    = TINA
             C(I,J)    = RHINA*VPRES(T(I,J))*1e2/RGAS/T(I,J)
          ELSE
             U(I,J)    = UINS
             IF (UPARAB) U(I,J)=2*UINS*(1-ABS(Y(J)/RADIUS)**2)
             T(I,J)    = TINS
             C(I,J)    = RHINS*VPRES(T(I,J))*1e2/RGAS/T(I,J)
          ENDIF
          DNCEL(I,J) = 0.0
 110    CONTINUE
 100  CONTINUE
C
C *** THERMOPHYSICAL PROPERTIES ****************************************
C
      CALL PROPS
C
C *** CALCULATE RESIDUAL SOURCES NORMALIZATION FACTORS *****************
C
      FLOWIN=0.0
      XMONIN=0.0
      DO 400 J=2,NJM1
         ARDEN=0.5*(DEN(1,J)+DEN(2,J))*R(J)*SNS(J)
         FLOWIN=FLOWIN+ARDEN*U(2,J)
         XMONIN=XMONIN+ARDEN*U(2,J)*U(2,J)
 400  CONTINUE
      SNORM=FLOWIN
      SNORU=XMONIN
      SNORV=XMONIN
      SNORT=FLOWIN*MAX(TINA, TINS)
      SNORC=FLOWIN
C
C *** END OF INIT SUBROUTINE *******************************************
C
      RETURN
      END
CC
C=======================================================================
C
C *** SUBROUTINE PROPS
C
C=======================================================================
C
      SUBROUTINE PROPS
      INCLUDE 'common_l.inc'
      INCLUDE 'prob.inc'
      INCLUDE 'aerosol.inc'
      REAL KAPAC
      LOGICAL CONSTT
C
C *** CALCULATE IN-FIELD PROPERTIES
C
      CPAIR = 1.0065e+3                           ! Cp for air (J/kg/K)
      VISCF = 1.e-5*(2.08-1.59)/(350.0-250.0)     ! Visc coef.
C
      DO 100 I=2,NIM1
      DO 100 J=2,NJM1
            DENI     = (PRESF+P(I,J))*MWA/RGAS/T(I,J)    ! Density (kg/m3)
            VISI     = 1.59e-5 + VISCF*(T(I,J)-250.0)    ! Viscosity
            KAPAC    = (4.39 + 0.071*T(I,J))*1.0D-3      ! Th.conduct.
            PR       = (P(I,J)+PRESF)/1.013e5            ! Pressure (atm)
            DIFFC    = (0.211/PR)*(T(I,J)/273d0)**1.94   ! Diffus (cm2/s)
            DIFFC    = DIFFC*1d-4                        ! cm2/s -> m2/s
         DEN(I,J) = DENI                              ! Density (kg/m3)
         VIS(I,J) = VISI                              ! Viscosity
         GAMH(I,J)= KAPAC/CPAIR
         GAMC(I,J)= DIFFC*DEN(I,J)                    ! viscosity/Sc
 100  CONTINUE
C
C *** AXIAL BOUNDARY PROPERTIES
C
      DO 110 I=2,NIM1       
         GAMH(I,1)  = GAMH(I,2)
         GAMH(I,NJ) = GAMH(I,NJM1)
         GAMC(I,1)  = GAMC(I,2)
         GAMC(I,NJ) = GAMC(I,NJM1)
         DEN(I,1)   = DEN(I,2)
         DEN(I,NJ)  = DEN(I,NJM1)
         VIS(I,1)   = VIS(I,2)
         VIS(I,NJ)  = VIS(I,NJM1)
 110  CONTINUE
C
C *** RADIAL BOUNDARY PROPERTIES
C
      DO 120 J=1,NJ
         GAMH(1,J)  = GAMH(2,J)
         GAMH(NI,J) = GAMH(NIM1,J)
         GAMC(1,J)  = GAMC(2,J)
         GAMC(NI,J) = GAMC(NIM1,J)
         DEN(1,J)   = DEN(2,J)
         DEN(NI,J)  = DEN(NIM1,J)
         VIS(1,J)   = VIS(2,J)
         VIS(NI,J)  = VIS(NIM1,J)
 120  CONTINUE
C
C *** END OF PROPS SUBROUTINE ******************************************
C
      RETURN
      END
CC
C=======================================================================
C
C *** SUBROUTINE MODU
C *** THIS SUBROUTINE PROVIDES THE MODIFICATIONS OF THE COEFFICIENTS
C     IN U-MOMENTUM EQUATION FOR THE NEAR BOUNDARY CONTROL VOLUMES
C
C=======================================================================
C
      SUBROUTINE MODU
      INCLUDE 'common_l.inc'
      INCLUDE 'prob.inc'
      INCLUDE 'aerosol.inc'
C
C *** TOP WALL *********************************************************
C
      J =NJM1
      YP=YV(NJ)-Y(NJM1)
      DO 100 I=3,NIM1
         VISCOS =0.5*(VIS(I-1,J)+VIS(I,J))
         TMULT  =VISCOS/YP
         SP(I,J)=SP(I,J)-TMULT*SEWU(I)*RV(J+1)
         AN(I,J)=0.0
 100  CONTINUE
C
C *** SYMMETRY AXIS ****************************************************
C
      J=2
      DO 200 I=2,NI
         U(I,J-1)=U(I,J)
         AS(I,J)=0.0
 200  CONTINUE
C
C *** OUTLET ***********************************************************
C
      ARDENT=0.0
      FLOW=0.0
      DO 300 J=2,NJM1
         ARDEN=0.5*(DEN(NIM1,J)+DEN(NIM1-1,J))*R(J)*SNS(J)
         ARDENT=ARDENT+ARDEN
         FLOW=FLOW+ARDEN*U(NIM1,J)
 300  CONTINUE
      UINC=(FLOWIN-FLOW)/ARDENT
      DO 310 J=2,NJM1
         U(NI,J)=U(NIM1,J)+UINC
 310  CONTINUE
C
C *** BUOYANCY *********************************************************
C
      DO 400 I=2,NIM1
         TBAR=0.0          ! Mean T
         DO 410 J=2,NJM1
            TBAR = TBAR + T(I,J)*(YV(J+1)-YV(J))*(YV(J+1)+YV(J))
 410     CONTINUE
         TBAR = TBAR/YV(NJ)**2
C
         DO 400 J=2,NJM1   ! Buoyancy source term
            VOL    = R(J)*SNS(J)*SEWU(I)
            DENB   = 0.5*(DEN(I,J)+DEN(I-1,J))
            TMPB   = 0.5*(T(I,J)+T(I-1,J))
            SU(I,J)=SU(I,J) - VOL*DENB*9.81*(TMPB/TBAR - 1.0)
 400  CONTINUE
C
C *** END OF MODU SUBROUTINE ******************************************
C
      RETURN
      END
CC
C=======================================================================
C
C *** SUBROUTINE MODV
C *** THIS SUBROUTINE PROVIDES THE MODIFICATIONS OF THE COEFFICIENTS
C     IN V-MOMENTUM EQUATION FOR THE NEAR BOUNDARY CONTROL VOLUMES
C
C=======================================================================
C
      SUBROUTINE MODV
      INCLUDE 'common_l.inc'
      INCLUDE 'prob.inc'
      INCLUDE 'aerosol.inc'
C
C *** TOP WALL *********************************************************
C
      J=NJM1
      DO 100 I=2,NIM1
         AN(I,J)=0.0
 100  CONTINUE
C
C *** SYMMETRY AXIS ****************************************************
C
      J=3
      DO 200 I=2,NIM1
         AS(I,J)=0.0
 200  CONTINUE
C
C *** OUTLET ***********************************************************
C
      I=NIM1
      DO 300 J=3,NJM1
         AE(I,J)=0.0
 300  CONTINUE
C
C *** END OF MODV SUBROUTINE *******************************************
C
      RETURN
      END
CC
C=======================================================================
C
C *** SUBROUTINE MODP
C *** THIS SUBROUTINE PROVIDES THE MODIFICATIONS OF THE COEFFICIENTS
C     IN PRESSURE-CORRECTION EQUATION FOR THE NEAR BOUNDARY CV'S
C
C=======================================================================
C
      SUBROUTINE MODP
      INCLUDE 'common_l.inc'
      INCLUDE 'prob.inc'
      INCLUDE 'aerosol.inc'
C
C ** END OF MODP SUBROUTINE *******************************************
C
      RETURN
      END
CC
C=======================================================================
C
C *** SUBROUTINE MODT
C *** THIS SUBROUTINE PROVIDES THE MODIFICATIONS OF THE COEFFICIENTS
C     IN TEMPERATURE EQUATION FOR THE NEAR BOUNDARY CONTROL VOLUMES
C
C=======================================================================
C
      SUBROUTINE MODT
      INCLUDE 'common_l.inc'
      INCLUDE 'prob.inc'
      INCLUDE 'aerosol.inc'
C
C *** SYMMETRY AXIS ***************************************************
C
      J=2
      DO 100 I=2,NIM1
         T(I,J-1)=T(I,J)
         AS(I,J)=0.0
 100  CONTINUE
C
C *** OUTLET ***********************************************************
C
      I=NIM1
      DO 200 J=2,NJM1
         T(I+1,J)=T(I,J)
         AE(I,J)=0.0
 200  CONTINUE
C
C *** TOP WALL ********************************************************
C
      J=NJM1
      DO 300 I=2,NIM1
         DN     =GAMH(I,J)*SEW(I)*RV(J+1)/(YV(J+1)-Y(J))
         SU(I,J)=SU(I,J)+DN*T(I,NJ)
         SP(I,J)=SP(I,J)-DN
         AN(I,J)=0.0
 300  CONTINUE
C
C *** ENERGY SOURCE FROM CONDENSATION ON PARTICLES ********************
C
      IF (ANTOT.GT.TINY) THEN
         DHVAP = 2.25e6*MWW                 ! Enth.of evap. (J/mole)
         CPAIR = 1.0065e+3                  ! Cp for air (J/kg/K)
         HTVAP = DHVAP/CPAIR
         DO 400 J=2,NJM1
         DO 400 I=2,NIM1
            SU(I,J) = SU(I,J) + DNCEL(I,J)*HTVAP
  400    CONTINUE
      ENDIF
C
C *** END OF MODT SUBROUTINE ******************************************
C
      RETURN
      END
CC
C=======================================================================
C
C *** SUBROUTINE MODC
C *** THIS SUBROUTINE PROVIDES THE MODIFICATIONS OF THE COEFFICIENTS
C     IN CONCENTRATION EQUATION FOR THE NEAR BOUNDARY CONTROL VOLUMES
C
C=======================================================================
C
      SUBROUTINE MODC
      INCLUDE 'common_l.inc'
      INCLUDE 'prob.inc'
      INCLUDE 'aerosol.inc'
C
C *** SYMMETRY AXIS ***************************************************
C
      J=2
      DO 100 I=2,NIM1
         C(I,J-1)=C(I,J)
         AS(I,J) =0.0
 100  CONTINUE
C
C *** OUTLET ***********************************************************
C
      I=NIM1
      DO 200 J=2,NJM1
         C(I+1,J)=C(I,J)
         AE(I,J) =0.0
 200  CONTINUE
C
C *** TOP WALL ********************************************************
C
      J=NJM1
      DO 300 I=2,NIM1
         CWALL   = VPRES(T(I,NJ))*1e2/RGAS/T(I,NJ)
         DN      = GAMC(I,J)*SEW(I)*RV(J+1)/(YV(J+1)-Y(J))
         SU(I,J) = SU(I,J)+DN*CWALL
         SP(I,J) = SP(I,J)-DN
         AN(I,J) = 0.0
 300  CONTINUE
C
C *** WATER VAPOR SINK FROM CONDENSATION ON PARTICLES ****************
C
      DO 400 J=2,NJM1
      DO 410 I=2,NIM1
         DNCEL(I,J) = ZERO                          ! dNwater in cell
         SU(I,J)    = SU(I,J)-DNCEL(I,J)*DEN(I,J)   ! Water condens. rate
 410  CONTINUE
 400  CONTINUE
C
C *** END OF MODC SUBROUTINE ******************************************
C
      RETURN
      END
CC
C=======================================================================
C
C *** SUBROUTINE USROUT
C *** THIS SUBROUTINE IS USED TO PROVIDE SPECIAL PRINTOUTS.
C
C=======================================================================
C
      SUBROUTINE USROUT (ICODE)
      INCLUDE 'common_l.inc'
      INCLUDE 'prob.inc'
      INCLUDE 'aerosol.inc'
      CHARACTER HDR(JJ)*19, REC*32000
      REAL SPR(II,JJ)
C
C *** INITIAL OUTPUT - HEADER ******************************************
C
      IF (ICODE.EQ.1) THEN
C 
C *** CALLED AFTER EVERY ITERATION **************************************
C
      ELSE IF (ICODE.EQ.2) THEN
C
C *** CALLED AFTER EVERY TIME STEP **************************************
C
      ELSE IF (ICODE.EQ.3) THEN
C
C *** FINAL OUTPUT ******************************************************
C
      ELSE IF (ICODE.EQ.4) THEN
C
      CALL CHRBLN (OUTFIL, IOE)        ! Find last non-blank character
      NIX=INDEX(OUTFIL,'.')-1          ! Find '.' (if exists)
      IF (NIX.GT.8) NIX=8              ! Last character before '.'
      IF (NIX.LE.0) NIX=MIN(IOE,8)     ! If no dot found, use nonblank chars
C
C *** CALCULATE SUPERSATURATION FIELD
C
      IF (INCONC) THEN
      DO 355 J=2,NJM1   
      DO 355 I=2,NIM1
         SPR(I,J)=(C(I,J)*RGAS*T(I,J)/VPRES(T(I,J))/1e2 - 1e0)*1e2
  355 CONTINUE
C
      IF (OUTPF) THEN !  SUPERSATURATION PROFILES
         OPEN(UNIT=IOAUX,FILE=OUTFIL(1:NIX)//'.SPR',STATUS='UNKNOWN')!,RECL=8000)
         DO 358 I=2,NI-1
            WRITE (HDR(I),'(A,1PE11.5,A)') '"S(Y=',Y(I),')" '
  358    CONTINUE
         WRITE (REC,*) '"X" ',(HDR(I),I=2,NJ-1)
         CALL CHRBLN (REC,IRE)
         WRITE (IOAUX,*) REC(1:IRE)
CCCCC         WRITE (IOAUX,*) '"X" ',(HDR(I),I=2,NJ-1)
         DO 356 I=2,NIM1
            WRITE (REC,*) X(I),(SPR(I,J),J=2,NJM1)
            CALL CHRBLN (REC,IRE)
            WRITE (IOAUX,*) REC(1:IRE)
CCCCC            WRITE(IOAUX,*) X(I),(SPR(I,J),J=2,NJM1)
  356    CONTINUE
         CLOSE(UNIT=IOAUX,STATUS='KEEP')
      ENDIF
C
      ENDIF
C
      IF (INCONC .AND. OUTPF) THEN  ! SUPERSATURATION CONTOURS
         OPEN(UNIT=IOAUX,FILE=OUTFIL(1:NIX)//'.SAT',STATUS='UNKNOWN')
         DO 350 J=2,NJM1
         DO 350 I=2,NIM1                         
            WRITE(IOAUX,*) X(I),Y(J),SPR(I,J)
  350    CONTINUE
         CLOSE(UNIT=IOAUX,STATUS='KEEP')
      ENDIF
C
C *** CALIBRATION BINARY FILE
C
      IF (INCONC .AND. OUTPF) THEN
         OPEN(UNIT=IOAUX,FILE=OUTFIL(1:NIX)//'.CLB',STATUS='UNKNOWN',
     &                   FORM='UNFORMATTED')
         WRITE(IOAUX) NI,NJAER
         WRITE(IOAUX) (XU(I),I=1,NI)
         WRITE(IOAUX) (YV(J),J=1,NJAER+1)
         WRITE(IOAUX) ((U(I,J),I=1,NI),J=1,NJAER)
         WRITE(IOAUX) ((V(I,J),I=1,NI),J=1,NJAER)
         WRITE(IOAUX) ((P(I,J),I=1,NI),J=1,NJAER)
         WRITE(IOAUX) ((T(I,J),I=1,NI),J=1,NJAER)
         WRITE(IOAUX) ((C(I,J),I=1,NI),J=1,NJAER)
         WRITE(IOAUX) PRESF,RHINA,TINA
         WRITE(IOAUX) PDNS,PMW,PVHF
         CLOSE(UNIT=IOAUX,STATUS='KEEP')
      ENDIF
C
C *** PRANDL NUMBER CONTOURS
C
      IF (.FALSE.) THEN    ! CC      IF (INCONC) THEN
         OPEN(UNIT=IOAUX,FILE=OUTFIL(1:NIX)//'.PR',STATUS='UNKNOWN',
     &                   RECL=8000)
         DO 380 J=2,NJM1
         DO 380 I=2,NIM1
            WRITE(IOAUX,*) X(I),Y(J),VIS(I,J)/GAMH(I,J)
  380    CONTINUE
         CLOSE(UNIT=IOAUX,STATUS='KEEP')
      ENDIF
C
C *** SCHMIDT NUMBER CONTOURS
C
      IF (.FALSE.) THEN    ! CC      IF (INCONC) THEN
         OPEN(UNIT=IOAUX,FILE=OUTFIL(1:NIX)//'.SC',STATUS='UNKNOWN',
     &                   RECL=8000)
         DO 390 J=2,NJM1
         DO 390 I=2,NIM1
            WRITE(IOAUX,*) X(I),Y(J),VIS(I,J)/GAMC(I,J)
  390    CONTINUE
         CLOSE(UNIT=IOAUX,STATUS='KEEP')
      ENDIF
C
C *** DENSITY CONTOURS
C
      IF (.FALSE.) THEN   ! CC      IF (INCONC) THEN
         OPEN(UNIT=IOAUX,FILE=OUTFIL(1:NIX)//'.DEN',STATUS='UNKNOWN',
     &                   RECL=8000)
         DO 400 J=2,NJM1
         DO 400 I=2,NIM1
            WRITE(IOAUX,*) X(I),Y(J),DEN(I,J)
  400    CONTINUE
         CLOSE(UNIT=IOAUX,STATUS='KEEP')
      ENDIF
C
C *** VISCOSITY CONTOURS
C
      IF (.FALSE.) THEN   ! CC      IF (INCONC) THEN
         OPEN(UNIT=IOAUX,FILE=OUTFIL(1:NIX)//'.VIS',STATUS='UNKNOWN',
     &                   RECL=8000)
         DO 410 J=2,NJM1
         DO 410 I=2,NIM1
            WRITE(IOAUX,*) X(I),Y(J),VIS(I,J)
  410    CONTINUE
         CLOSE(UNIT=IOAUX,STATUS='KEEP')
      ENDIF
C
C *** VELOCITIES *****************************************************
C
      IF (INPLOT .AND. OUTPF) THEN
         OPEN(UNIT=IOAUX,FILE=OUTFIL(1:NIX)//'.VEL',STATUS='UNKNOWN')
         DO 920 J=2,NJM1
         FY = (Y(J)-Y(J-1))/(Y(J+1)-Y(J-1))
         DO 920 I=2,NIM1
            FX  = (X(I)-X(I-1))/(X(I+1)-X(I-1))
            UMN = FX*U(I+1,J) + (1.0-FX)*U(I,J)
            VMN = FY*V(I,J+1) + (1.0-FY)*V(I,J)
            WRITE (IOAUX,*) X(I), Y(J), UMN, VMN
  920    CONTINUE
         CLOSE(UNIT=IOAUX,STATUS='KEEP')
      ENDIF
C
C *** PRESURE CONTOURS ************************************************
C
      IF (INPRES .AND. OUTPF) THEN
         OPEN(UNIT=IOAUX,FILE=OUTFIL(1:NIX)//'.PRS',STATUS='UNKNOWN')
         DO 930 J=2,NJM1
         DO 930 I=2,NIM1
            WRITE(IOAUX,*) X(I),Y(J),P(I,J)
  930    CONTINUE
         CLOSE(UNIT=IOAUX,STATUS='KEEP')
      ENDIF
C
C *** TEMPERATURE CONTOURS ********************************************
C
      IF (INTEMP .AND. OUTPF) THEN
         OPEN(UNIT=IOAUX,FILE=OUTFIL(1:NIX)//'.TMP',STATUS='UNKNOWN')
         DO 940 J=2,NJM1
         DO 940 I=2,NIM1
            WRITE(IOAUX,*) X(I),Y(J),T(I,J)
  940    CONTINUE
         CLOSE(UNIT=IOAUX,STATUS='KEEP')
      ENDIF
C
C *** CONCENTRATION CONTOURS ******************************************
C
      IF (INCONC .AND. OUTPF) THEN
         OPEN(UNIT=IOAUX,FILE=OUTFIL(1:NIX)//'.CON',STATUS='UNKNOWN')
         DO 950 J=2,NJM1
         DO 950 I=2,NIM1
            WRITE(IOAUX,*) X(I),Y(J),C(I,J)
  950    CONTINUE
         CLOSE(UNIT=IOAUX,STATUS='KEEP')
      ENDIF
C
C *** OUTPUT AVG SUPERSATURATION TO STDOUT ***************************
C
      NISMP = MAX(NI*INT(SOAVGP/100.),0) ! Smooth over last SOAVGP% of instr.
      AVGS  = ZERO                       ! Zero average counter
      J     = 2                          ! Centerline supersaturation
      DO 960 I= NIM1, NIM1-NISMP, -1
         AVGS = AVGS + SPR(I,J)
  960 CONTINUE
      AVGS = AVGS/FLOAT(NISMP+1)
      WRITE (6,*) AVGS                ! Write exit (average) S to STDOUT
      CALL FLUSH                      ! Flush output
C     Call a C function to flush stdout
      CALL trvflush()
C
      ENDIF
C
C *** END OF USROUT SUBROUTINE *****************************************
C
      RETURN
      END



C=======================================================================
C
C *** BLOCK DATA SUBPROGRAM
C *** THIS SUBROUTINE PROVIDES INITIAL (DEFAULT) VALUES TO PROGRAM
C     PARAMETERS VIA DATA STATEMENTS
C
C=======================================================================
C
      BLOCK DATA
      INCLUDE 'common_l.inc'
C
C *** GRID VARIABLES ***************************************************
C
      DATA POLAR/.FALSE./
      DATA NI,NJ/3,3/
      DATA X(1),Y(1),X(2),Y(2),X(3),Y(3)/2*-0.5,2*0.5,2*1.5/
C
C *** TIME VARIABLES ***************************************************
C
      DATA STEADY/.TRUE./
      DATA NT,IT,DT,TIME/1,0,1.0,0.0/
C
C *** DEPENDENT VARIABLES **********************************************
C
      DATA INCALU,INCALV,INCALP,INCALT,INCALC/5*.FALSE./
C
C *** THERMO-PHYSICAL PROPERTIES ***************************************
C
      DATA INPRO/.FALSE./
      DATA VISCOS,DENSIT,PRANDT,SCHMID/4*1.0/
C
C *** TERMINATION CRITERIA FOR ITERATIONS ******************************
C
      DATA MAXIT,NITER,SORMAX/1,0,1.0E-6/
      DATA SNORM,SNORU,SNORV,SNORT,SNORC/5*1.0/
      DATA RESORU,RESORV,RESORM,RESORT,RESORC/5*0.0/
C
C *** UNDER-RELAXATION *************************************************
C
      DATA URFU,URFV,URFP,URFT,URFC/5*1.0/
C
C *** FIXED PRESSURE POINT *********************************************
C
      DATA IPREF,JPREF/2,2/
C
C *** EQUATION SOLVER SELECTION ****************************************
C
      DATA ISLVU,ISLVV,ISLVP,ISLVT,ISLVC/5*1/
      DATA ALPHA/0.8/
C
C *** TERMINATION CRITERIA FOR EQUATION SOLVER *************************
C
      DATA LITERU,LITERV,LITERP,LITERT,LITERC/1,1,1,1,1/
      DATA SORU,SORV,SORP,SORT,SORC/5*1.0E-6/
C
C *** INPUT-OUTPUT PARAMETERS ******************************************
C
      DATA IOINP,IOUSR,IOFIL,IOAUX,IOCON/1,2,3,4,6/
      DATA OUTFIL,RSTFIL,SAVFIL/'RESULT.','PHIDAT.','PHIDAT.'/
      DATA HEDU,HEDV,HEDP,HEDT,HEDC/
     &     '*-*- U- VELOCITY -*-*','*-*- V- VELOCITY -*-*',
     &     '*-*-*- PRESSURE *-*-*','*-*- TEMPERATURE -*-*',
     &     '*-* CONCENTRATION *-*'/
      DATA INREAD,INWRIT,INPLOT,INPRES,INTEMP,INCONC/6*.FALSE./
      DATA IMON,JMON,NITPRI,NSTPRI/2,2,1,1/
C
C *** GENERAL VARIABLES ************************************************
C
      DATA GREAT,TINY,PI/1.0E30,1.0E-30,3.1415927/
C
C *** END OF BLOCK DATA SUBPROGRAM *************************************
C
      END
C=======================================================================
C
C *** BLOCK DATA BLKPAR
C *** THIS SUBROUTINE PROVIDES INITIAL (DEFAULT) VALUES TO PROGRAM
C     PARAMETERS VIA DATA STATEMENTS
C
C=======================================================================
C
      BLOCK DATA BLKPAR
      INCLUDE 'aerosol.inc'
C
C *** Model parameters
C
      DATA RGAS /8.314/               ! Gas constant (J/mol/K)
      DATA GRAV /9.81/                ! Gravitational constant (m/s2)
      DATA MWA  /29.e-3/              ! Mol.weight of air (kg/mol)
      DATA MWW  /18.e-3/              ! Mol.weight of water (kg/mol)
C
      DATA NSEC /10/                  ! Default # of sections
      DATA ZERO /0.0e0/               ! Double Precision zero
C
C *** LSODE, VODE parameters
C
      DATA ISTATE /1/                 ! Integer flag (1=start integration)  
      DATA ITOL   /1/                 ! Atol (below) is a scalar
      DATA RTOL   /1.0d-5/            ! Relative tolerance parameter 
      DATA ATOL   /1.0d-25/           ! Absolute tolerance parameter 
      DATA ITASK  /1/                 ! Output values of q at t = tout.
      DATA IOPT   /0/                 ! No optional inputs used.
      DATA MF     /22/                ! Integration method
C
C *** END OF BLOCK DATA SUBPROGRAM *************************************
C
      END

C=======================================================================
C
C *** PROGRAM MAIN
C *** THIS IS THE MAIN ROUTINE OF THE PROGRAM.
C     IT IS USED TO PROVIDE INFORMATIONS SPECIFIC TO EACH NEW PROBLEM.
C     SUBROUTINE CONTRO IS CALLED TO SOLVE THE PROBLEM.
C
C=======================================================================
C
      PROGRAM MAIN
      INCLUDE 'common_l.inc'
C
      CHARACTER PROBFIL*256
      LOGICAL   EOF, FRST
C
C *** READ INPUT FILE PATH *********************************************
C
      READ (*,'(A)') PROBFIL
      EOF  = .FALSE.          ! EOF flag for STDIN input
      FRST = .TRUE.           ! First run follows
C
C *** MODEL "LOOP" *****************************************************
C
      DO WHILE (.TRUE.)
         CALL INPUT (PROBFIL, EOF)  ! Input parameters
         IF (EOF) GOTO 100          ! If EOF encountered, terminate
C
         IF (FRST) THEN             ! First time around things
           CALL GEOM           ! Setup geometry
           DO 110 I=1,NI       ! Initial field values
           DO 110 J=1,NJ
              U(I,J)=0.0
              V(I,J)=0.0
              P(I,J)=0.0
              T(I,J)=0.0
              C(I,J)=0.0
              PP(I,J)=0.0
              DEN(I,J)=DENSIT
              VIS(I,J)=VISCOS
              GAMH(I,J)=VISCOS/PRANDT
              GAMC(I,J)=VISCOS/SCHMID
 110       CONTINUE
           FRST = .FALSE.       ! Set "first" time flag to .false.
         ENDIF
C
         CALL CONTRO                ! Solve problem
      ENDDO
C
C *** END OF MAIN PROGRAM **********************************************
C
100   STOP
      END
CC
C=======================================================================
C
C *** SUBROUTINE CONTRO
C *** THIS IS THE PROGRAM CONTROL SUBROUTINE
C
C=======================================================================
C
      SUBROUTINE CONTRO
      INCLUDE 'common_l.inc'
C
      CALL INIT               ! Initialize "stuff"
      NITER = 0
C
C *** READ IN OLD FLOW FIELD *******************************************
C
      IF (INREAD) THEN
         OPEN(UNIT=IOAUX,FILE=RSTFIL,STATUS='OLD',FORM='UNFORMATTED',
     &        ERR=120)
         IF (STEADY) THEN
           READ(IOAUX) NI,NJ,NITER
           MAXIT=NITER+MAXIT
         ELSE
           READ(IOAUX) NI,NJ,IT,TIME
           NT=IT+NT
         ENDIF
         READ(IOAUX) ((U(I,J),I=1,NI),J=1,NJ)
         READ(IOAUX) ((V(I,J),I=1,NI),J=1,NJ)
         READ(IOAUX) ((P(I,J),I=1,NI),J=1,NJ)
         READ(IOAUX) ((T(I,J),I=1,NI),J=1,NJ)
         READ(IOAUX) ((C(I,J),I=1,NI),J=1,NJ)
         READ(IOAUX) ((VIS(I,J),I=1,NI),J=1,NJ)
         READ(IOAUX) ((DEN(I,J),I=1,NI),J=1,NJ)
         READ(IOAUX) ((GAMH(I,J),I=1,NI),J=1,NJ)
         READ(IOAUX) ((GAMC(I,J),I=1,NI),J=1,NJ)
         CLOSE(UNIT=IOAUX,STATUS='KEEP')
         GOTO 121
  120    CONTINUE
  121    CONTINUE
      ENDIF
C
C *** UPDATE FLUID PROPERTIES ******************************************
C
      CALL PROPS
C
C *** INITIAL OUTPUT - FIELDS ******************************************
C
      CALL USROUT(1)
C
C-----------------------------------------------------------------------
CHAPTER 2 2 2 2 2 2 2 2 2 2 2  SOLUTION LOOP  2 2 2 2 2 2 2 2 2 2 2 2 2
C-----------------------------------------------------------------------
C
C *** TIME STEP LOOP ***************************************************
C
      IF (STEADY) GOTO 215
 200  IT=IT+1
      TIME=TIME+DT
      DO 210 I=1,NI
      DO 210 J=1,NJ
         UO(I,J)=U(I,J)
         VO(I,J)=V(I,J)
         TO(I,J)=T(I,J)
         CO(I,J)=C(I,J)
         DENO(I,J)=DEN(I,J)
 210  CONTINUE
      NITER=0
C
C *** ITERATION LOOP ***************************************************
C
 215  CONTINUE
 220  NITER=NITER+1
C
C *** SOLVE EQUATIONS **************************************************
C
      IF (INCALU) CALL CALCU
      IF (INCALV) CALL CALCV
      IF (INCALP) CALL CALCP
      IF (INCALT) CALL CALCT
      IF (INCALC) CALL CALCC
C
C *** UPDATE FLUID PROPERITIES *****************************************
C
      IF (INPRO) CALL PROPS
C
C *** CALCULATE NORMALIZED RESIDUALS ***********************************
C
      RESORM=RESORM/SNORM
      RESORU=RESORU/SNORU
      RESORV=RESORV/SNORV
      RESORT=RESORT/SNORT
      RESORC=RESORC/SNORC
C
C *** CONVERGENCE TESTS ************************************************
C
      SORCE=AMAX1(RESORM,RESORU,RESORV,RESORT,RESORC)
      IF (SORCE.LE.SORMAX) THEN
        ICONV=2
      ELSE IF (NITER.EQ.MAXIT) THEN
        ICONV=1
      ELSE
        ICONV=0
      ENDIF
C
C *** PRINTOUT SPECIFIED BY ITERATION PRINTING INDEX *******************
C
      CALL USROUT(2)
C
C *** END OF ITERATION LOOP ********************************************
C
      IF (ICONV.EQ.0) GOTO 220
C
C *** PRINTOUT SPECIFIED BY TIME PRINTING INDEX ************************
C
      CALL USROUT(3)
C
C *** END OF TIME STEP LOOP ********************************************
C
      IF (.NOT.STEADY.AND.IT.LT.NT) GOTO 200
C
C-----------------------------------------------------------------------
CHAPTER 3 3 3 3 3 3 3 3 3 3  FINAL OPERATIONS  3 3 3 3 3 3 3 3 3 3 3 3 3
C-----------------------------------------------------------------------
C
 300  CONTINUE
      CALL USROUT(4)
C
C *** COPY OUT FLOW FIELD **********************************************
C
      IF (INWRIT) THEN
        OPEN(IOAUX,FILE=SAVFIL,STATUS='UNKNOWN',FORM='UNFORMATTED')
        IF (STEADY) THEN
          WRITE(IOAUX) NI,NJ,NITER
        ELSE
          WRITE(IOAUX) NI,NJ,IT,TIME
        ENDIF
        WRITE(IOAUX) ((U(I,J),I=1,NI),J=1,NJ)
        WRITE(IOAUX) ((V(I,J),I=1,NI),J=1,NJ)
        WRITE(IOAUX) ((P(I,J),I=1,NI),J=1,NJ)
        WRITE(IOAUX) ((T(I,J),I=1,NI),J=1,NJ)
        WRITE(IOAUX) ((C(I,J),I=1,NI),J=1,NJ)
        WRITE(IOAUX) ((VIS(I,J),I=1,NI),J=1,NJ)
        WRITE(IOAUX) ((DEN(I,J),I=1,NI),J=1,NJ)
        WRITE(IOAUX) ((GAMH(I,J),I=1,NI),J=1,NJ)
        WRITE(IOAUX) ((GAMC(I,J),I=1,NI),J=1,NJ)
        CLOSE(UNIT=IOAUX,STATUS='KEEP')
      ENDIF
C
C *** END OF SUBROUTINE CONTRO *****************************************
C
      RETURN
      END
CC
C=======================================================================
C
C *** SUBROUTINE GEOM
C *** THIS SUBROUTINE CALCULATES GEOMETRICAL QUANTITIES
C
C=======================================================================
C
      SUBROUTINE GEOM
      INCLUDE 'common_l.inc'
C
C-----------------------------------------------------------------------
CHAPTER 1 1 1 1 1 1 1 1 1 1 1  SCALAR CELL  1 1 1 1 1 1 1 1 1 1 1 1 1 1
C-----------------------------------------------------------------------
C
      NIM1=NI-1
      NJM1=NJ-1
C
      DO 100 J=1,NJ
         R(J)=1.0
         IF (POLAR) R(J)=Y(J)
  100 CONTINUE
C
      DXPW(1)=0.0
      DXEP(NI)=0.0
      DO 110 I=1,NIM1
         DXEP(I)=X(I+1)-X(I)
         DXPW(I+1)=DXEP(I)
  110 CONTINUE
C
      DYPS(1)=0.0
      DYNP(NJ)=0.0
      DO 120 J=1,NJM1
         DYNP(J)=Y(J+1)-Y(J)
         DYPS(J+1)=DYNP(J)
  120 CONTINUE
C
      SEW(1)=0.0
      SEW(NI)=0.0
      DO 130 I=2,NIM1
         SEW(I)=0.5*(DXEP(I)+DXPW(I))
  130 CONTINUE
C
      SNS(1)=0.0
      SNS(NJ)=0.0
      DO 140 J=2,NJM1
         SNS(J)=0.5*(DYNP(J)+DYPS(J))
  140 CONTINUE
C
C-----------------------------------------------------------------------
CHAPTER 2 2 2 2 2 2 2 2 2 2 2 2  U - CELL  2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
C-----------------------------------------------------------------------
C
      XU(1)=0.0
      DO 200 I=2,NI
         XU(I)=0.5*(X(I)+X(I-1))
  200 CONTINUE
C
      DXPWU(1)=0.0
      DXPWU(2)=0.0
      DXEPU(1)=0.0
      DXEPU(NI)=0.0
      DO 210 I=2,NIM1
         DXEPU(I)=XU(I+1)-XU(I)
         DXPWU(I+1)=DXEPU(I)
  210 CONTINUE
C
      SEWU(1)=0.0
      SEWU(2)=0.0
      DO 220 I=3,NIM1
         SEWU(I)=0.5*(DXEPU(I)+DXPWU(I))
  220 CONTINUE
C
C-----------------------------------------------------------------------
CHAPTER 3 3 3 3 3 3 3 3 3 3 3 3  V - CELL  3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
C-----------------------------------------------------------------------
C
      YV(1)=0.0
      RV(1)=0.0
      DO 300 J=2,NJ
         RV(J)=0.5*(R(J)+R(J-1))
         RCV(J)=0.5*(RV(J)+RV(J-1))
         YV(J)=0.5*(Y(J)+Y(J-1))
  300 CONTINUE
C
      DYPSV(1)=0.0
      DYPSV(2)=0.0
      DYNPV(NJ)=0.0
      DO 310 J=2,NJM1
         DYNPV(J)=YV(J+1)-YV(J)
         DYPSV(J+1)=DYNPV(J)
  310 CONTINUE
C
      SNSV(1)=0.0
      SNSV(2)=0.0
      SNSV(NJ)=0.0
      DO 320 J=3,NJM1
         SNSV(J)=0.5*(DYNPV(J)+DYPSV(J))
  320 CONTINUE
C
C *** END OF GEOM SUBROUTINE *******************************************
C
      RETURN
      END
CC
C=======================================================================
C
C *** SUBROUTINE CALCU
C *** THIS SUBROUTINE SOLVES THE U-MOMENTUM EQUATION
C
C=======================================================================
C
      SUBROUTINE CALCU
      INCLUDE 'common_l.inc'
C
C-----------------------------------------------------------------------
CHAPTER 1 1 1 1 1 1 1 1  ASSEMBLY OF COEFFICIENTS  1 1 1 1 1 1 1 1 1 1 1
C-----------------------------------------------------------------------
C
      CT=0.0
      CTO=0.0
      DO 100 I=3,NIM1
        DO 110 J=2,NJM1
C
C *** COMPUTE AREAS AND VOLUME *****************************************
C
          AREAN=RV(J+1)*SEWU(I)
          AREAS=RV(J)*SEWU(I)
          AREAEW=R(J)*SNS(J)
          VOL=R(J)*SEWU(I)*SNS(J)
C
C *** CALCULATE CONVECTION COEFFICIENTS ********************************
C
          GN=0.5*(DEN(I,J+1)+DEN(I,J))*V(I,J+1)
          GNW=0.5*(DEN(I-1,J)+DEN(I-1,J+1))*V(I-1,J+1)
          GS=0.5*(DEN(I,J-1)+DEN(I,J))*V(I,J)
          GSW=0.5*(DEN(I-1,J)+DEN(I-1,J-1))*V(I-1,J)
          GE=0.5*(DEN(I+1,J)+DEN(I,J))*U(I+1,J)
          GP=0.5*(DEN(I,J)+DEN(I-1,J))*U(I,J)
          GW=0.5*(DEN(I-1,J)+DEN(I-2,J))*U(I-1,J)
          CN=0.5*(GN+GNW)*AREAN
          CS=0.5*(GS+GSW)*AREAS
          CE=0.5*(GE+GP)*AREAEW
          CW=0.5*(GP+GW)*AREAEW
C
C *** CALCULATE DIFFUSION COEFFICIENTS *********************************
C
          VISN=0.25*(VIS(I,J)+VIS(I,J+1)+VIS(I-1,J)+VIS(I-1,J+1))
          VISS=0.25*(VIS(I,J)+VIS(I,J-1)+VIS(I-1,J)+VIS(I-1,J-1))
          DN=VISN*AREAN/DYNP(J)
          DS=VISS*AREAS/DYPS(J)
          DE=VIS(I,J)*AREAEW/DXEPU(I)
          DW=VIS(I-1,J)*AREAEW/DXPWU(I)
C
C *** ASSEMBLE MAIN COEFFICIENTS ***************************************
C
          AN(I,J)=AMAX1(ABS(0.5*CN),DN)-0.5*CN
          AS(I,J)=AMAX1(ABS(0.5*CS),DS)+0.5*CS
          AE(I,J)=AMAX1(ABS(0.5*CE),DE)-0.5*CE
          AW(I,J)=AMAX1(ABS(0.5*CW),DW)+0.5*CW
C
C *** CALCULATE UNSTEADY TERM ******************************************
C
          IF (.NOT.STEADY) THEN
            CT=0.5*(DEN(I,J)+DEN(I-1,J))*VOL/DT
            CTO=0.5*(DENO(I,J)+DENO(I-1,J))*VOL/DT
          ENDIF
          SU(I,J)=CTO*UO(I,J)
          SP(I,J)=-CTO
C
C *** ADD CONTRIBUTION TO SOURCE TERMS DUE TO THE FACT *****************
C *** THAT CONTINUITY IS NOT EXACTLY SATISFIED *************************
C
          SMP=CT-CTO+CN-CS+CE-CW
          CP=AMAX1(0.0,SMP)
          SU(I,J)=SU(I,J)+CP*U(I,J)
          SP(I,J)=SP(I,J)-CP
C
C *** ASSEMBLE LINEARIZED SOURCE TERMS *********************************
C
          DU(I,J)=AREAEW
          SU(I,J)=SU(I,J)+DU(I,J)*(P(I-1,J)-P(I,J))
          DUDXP=(U(I+1,J)-U(I,J))/SEW(I)
          DUDXM=(U(I,J)-U(I-1,J))/SEW(I-1)
          SU(I,J)=SU(I,J)+(VIS(I,J)*DUDXP-VIS(I-1,J)*DUDXM)/SEWU(I)*VOL
          GAMP=0.25*(VIS(I,J)+VIS(I-1,J)+VIS(I,J+1)+VIS(I-1,J+1))
          DVDXP=RV(J+1)*(V(I,J+1)-V(I-1,J+1))/DXEP(I)
          GAMM=0.25*(VIS(I,J)+VIS(I-1,J)+VIS(I,J-1)+VIS(I-1,J-1))
          DVDXM=RV(J)*(V(I,J)-V(I-1,J))/DXEP(I)
          SU(I,J)=SU(I,J)+(GAMP*DVDXP-GAMM*DVDXM)/SNS(J)/R(J)*VOL
C
  110   CONTINUE
  100 CONTINUE
C
C-----------------------------------------------------------------------
CHAPTER 2 2 2 2 2 2 2 2  PROBLEM MODIFICATIONS  2 2 2 2 2 2 2 2 2 2 2 2
C-----------------------------------------------------------------------
C
      CALL MODU
C
C-----------------------------------------------------------------------
CHAPTER 3 3 3 3 3 3 3  FINAL ASSEMPLY OF COEFFICIENT  3 3 3 3 3 3 3 3 3
C-----------------------------------------------------------------------
C
      RESORU=0.0
      DO 300 I=3,NIM1
        DO 310 J=2,NJM1
C
C *** CALCULATE AP COEFFICIENTS ****************************************
C
          AP(I,J)=AN(I,J)+AS(I,J)+AE(I,J)+AW(I,J)-SP(I,J)
C
C *** CALCULATE COEFFICIENTS OF VELOCITY CORRECTION ********************
C
          DU(I,J)=DU(I,J)/AP(I,J)
C
C *** RAITHBY'S MODIFICATION *******************************************
C
C          DU(I,J)=DU(I,J)/(AP(I,J)-URFU*(AP(I,J)+SP(I,J))+1.0E-20)
C
C *** CALCULATE RESIDUALS **********************************************
C
          RESOR=AN(I,J)*U(I,J+1)+AS(I,J)*U(I,J-1)+AE(I,J)*U(I+1,J)+
     &          AW(I,J)*U(I-1,J)-AP(I,J)*U(I,J)+SU(I,J)
          VOL=R(J)*SEWU(I)*SNS(J)
          SORVOL=GREAT*VOL
          IF (-SP(I,J).GT.0.5*SORVOL) RESOR=RESOR/SORVOL
          RESORU=RESORU+ABS(RESOR)
C
C *** UNDER-RELAXATION *************************************************
C
          AP(I,J)=AP(I,J)/URFU
          SU(I,J)=SU(I,J)+(1.0-URFU)*AP(I,J)*U(I,J)
          DU(I,J)=DU(I,J)*URFU
C
  310   CONTINUE
  300 CONTINUE
C
C-----------------------------------------------------------------------
CHAPTER 4 4 4 4 4 4 4  SOLUTION OF DIFFERENCE EQUATIONS  4 4 4 4 4 4 4 4
C-----------------------------------------------------------------------
C
      IF (ISLVU.EQ.1) THEN
        CALL ADISOL(3,2,NIM1,NJM1,U,LITERU,SORU)
      ELSE IF(ISLVU.EQ.2) THEN
        CALL SIMSOL(3,2,NIM1,NJM1,U,LITERU)
      ENDIF
C
C *** END OF CALCU SUBROUTINE ******************************************
C
      RETURN
      END
CC
C=======================================================================
C
C *** SUBROUTINE CALCV
C *** THIS SUBROUTINE SOLVES THE V-MOMENTUM EQUATION
C
C=======================================================================
C
      SUBROUTINE CALCV
      INCLUDE 'common_l.inc'
C
C-----------------------------------------------------------------------
CHAPTER 1 1 1 1 1 1 1 1  ASSEMBLY OF COEFFICIENTS  1 1 1 1 1 1 1 1 1 1 1
C-----------------------------------------------------------------------
C
      CT=0.0
      CTO=0.0
      DO 100 I=2,NIM1
        DO 110 J=3,NJM1
C
C *** COMPUTE AREAS AND VOLUME *****************************************
C
          AREAN=RCV(J+1)*SEW(I)
          AREAS=RCV(J)*SEW(I)
          AREAEW=RV(J)*SNSV(J)
          VOL=RV(J)*SEW(I)*SNSV(J)
C
C *** CALCULATE CONVECTION COEFFICIENTS ********************************
C
          GN=0.5*(DEN(I,J+1)+DEN(I,J))*V(I,J+1)
          GP=0.5*(DEN(I,J)+DEN(I,J-1))*V(I,J)
          GS=0.5*(DEN(I,J-1)+DEN(I,J-2))*V(I,J-1)
          GE=0.5*(DEN(I+1,J)+DEN(I,J))*U(I+1,J)
          GSE=0.5*(DEN(I,J-1)+DEN(I+1,J-1))*U(I+1,J-1)
          GW=0.5*(DEN(I,J)+DEN(I-1,J))*U(I,J)
          GSW=0.5*(DEN(I,J-1)+DEN(I-1,J-1))*U(I,J-1)
          CN=0.5*(GN+GP)*AREAN
          CS=0.5*(GP+GS)*AREAS
          CE=0.5*(GE+GSE)*AREAEW
          CW=0.5*(GW+GSW)*AREAEW
C
C *** CALCULATE DIFUSION COEFFICIENTS **********************************
C
          VISE=0.25*(VIS(I,J)+VIS(I+1,J)+VIS(I,J-1)+VIS(I+1,J-1))
          VISW=0.25*(VIS(I,J)+VIS(I-1,J)+VIS(I,J-1)+VIS(I-1,J-1))
          DN=VIS(I,J)*AREAN/DYNPV(J)
          DS=VIS(I,J-1)*AREAS/DYPSV(J)
          DE=VISE*AREAEW/DXEP(I)
          DW=VISW*AREAEW/DXPW(I)
C
C *** ASSEMBLE MAIN COEFFICIENTS ***************************************
C
          AN(I,J)=AMAX1(ABS(0.5*CN),DN)-0.5*CN
          AS(I,J)=AMAX1(ABS(0.5*CS),DS)+0.5*CS
          AE(I,J)=AMAX1(ABS(0.5*CE),DE)-0.5*CE
          AW(I,J)=AMAX1(ABS(0.5*CW),DW)+0.5*CW
C
C *** CALCULATE UNSTEADY TERM ******************************************
C
          IF (.NOT.STEADY) THEN
            CT=0.5*(DEN(I,J)+DEN(I,J-1))*VOL/DT
            CTO=0.5*(DENO(I,J)+DENO(I,J-1))*VOL/DT
          ENDIF
          SU(I,J)=CTO*VO(I,J)
          SP(I,J)=-CTO
C
C *** ADD CONTRIBUTION TO SOURCE TERMS DUE TO THE FACT *****************
C *** THAT CONTINUITY IS NOT EXACTLY SATISFIED *************************
C
          SMP=CT-CTO+CN-CS+CE-CW
          CP=AMAX1(0.0,SMP)
          SU(I,J)=SU(I,J)+CP*V(I,J)
          SP(I,J)=SP(I,J)-CP
C
C *** ASSEMBLE LINEARIZED SOURCE TERMS *********************************
C
          DV(I,J)=0.5*(AREAN+AREAS)
          SU(I,J)=SU(I,J)+DV(I,J)*(P(I,J-1)-P(I,J))
          IF (POLAR) SP(I,J)=SP(I,J)-VIS(I,J)*VOL/RV(J)**2
          DUDYP=(U(I+1,J)-U(I+1,J-1))/DYPS(J)
          GAMP=0.25*(VIS(I,J)+VIS(I+1,J)+VIS(I,J-1)+VIS(I+1,J-1))
          GAMM=0.25*(VIS(I,J)+VIS(I-1,J)+VIS(I,J-1)+VIS(I-1,J-1))
          DUDYM=(U(I,J)-U(I,J-1))/DYPS(J)
          SU(I,J)=SU(I,J)+(GAMP*DUDYP-GAMM*DUDYM)/SEW(I)*VOL
          DVDYP=(V(I,J+1)-V(I,J))/SNS(J)
          RGAMP=VIS(I,J)*R(J)
          DVDYM=(V(I,J)-V(I,J-1))/SNS(J-1)
          RGAMM=VIS(I,J-1)*R(J-1)
          SU(I,J)=SU(I,J)+(RGAMP*DVDYP-RGAMM*DVDYM)/(R(J)*SNS(J))*VOL
C
  110   CONTINUE
  100 CONTINUE
C
C-----------------------------------------------------------------------
CHAPTER 2 2 2 2 2 2 2 2  PROBLEM MODIFICATIONS  2 2 2 2 2 2 2 2 2 2 2 2
C-----------------------------------------------------------------------
C
      CALL MODV
C
C-----------------------------------------------------------------------
CHAPTER 3 3 3 3 3 3 3  FINAL ASSEMPLY OF COEFFICIENT  3 3 3 3 3 3 3 3 3
C-----------------------------------------------------------------------
C
      RESORV=0.0
      DO 300 I=2,NIM1
        DO 310 J=3,NJM1
C
C *** CALCULATE AP COEFFICIENTS ****************************************
C
          AP(I,J)=AN(I,J)+AS(I,J)+AE(I,J)+AW(I,J)-SP(I,J)
C
C *** CALCULATE COEFFICIENTS OF VELOCITY CORRECTION ********************
C
          DV(I,J)=DV(I,J)/AP(I,J)
C
C *** RAITHBY'S MODIFICATION *******************************************
C
C          DV(I,J)=DV(I,J)/(AP(I,J)-URFV*(AP(I,J)+SP(I,J))+1.0E-20)
C
C *** CALCULATE RESIDUALS **********************************************
C
          RESOR=AN(I,J)*V(I,J+1)+AS(I,J)*V(I,J-1)+AE(I,J)*V(I+1,J)+
     &          AW(I,J)*V(I-1,J)-AP(I,J)*V(I,J)+SU(I,J)
          VOL=RV(J)*SEW(I)*SNSV(J)
          SORVOL=GREAT*VOL
          IF (-SP(I,J).GT.0.5*SORVOL) RESOR=RESOR/SORVOL
          RESORV=RESORV+ABS(RESOR)
C
C *** UNDER-RELAXATION *************************************************
C
          AP(I,J)=AP(I,J)/URFV
          SU(I,J)=SU(I,J)+(1.0-URFV)*AP(I,J)*V(I,J)
          DV(I,J)=DV(I,J)*URFV
C
  310   CONTINUE
  300 CONTINUE
C
C-----------------------------------------------------------------------
CHAPTER 4 4 4 4 4 4 4  SOLUTION OF DIFFERENCE EQUATIONS  4 4 4 4 4 4 4 4
C-----------------------------------------------------------------------
C
      IF (ISLVV.EQ.1) THEN
        CALL ADISOL(2,3,NIM1,NJM1,V,LITERV,SORV)
      ELSE IF (ISLVV.EQ.2) THEN
        CALL SIMSOL(2,3,NIM1,NJM1,V,LITERV)
      ENDIF
C
C *** END OF CALCV SUBROUTINE ******************************************
C
      RETURN
      END
CC
C=======================================================================
C
C *** SUBROUTINE CALCP
C *** THIS SUBROUTINE SOLVES THE PRESSURE CORRECTION EQUATION
C
C=======================================================================
C
      SUBROUTINE CALCP
      INCLUDE 'common_l.inc'
C
C-----------------------------------------------------------------------
CHAPTER 1 1 1 1 1 1 1 1  ASSEMBLY OF COEFFICIENTS  1 1 1 1 1 1 1 1 1 1 1
C-----------------------------------------------------------------------
C
      CT=0.0
      CTO=0.0
      RESORM=0.0
      DO 100 I=2,NIM1
        DO 110 J=2,NJM1
C
C *** COMPUTE AREAS AND VOLUME *****************************************
C
          AREAN=RV(J+1)*SEW(I)
          AREAS=RV(J)*SEW(I)
          AREAEW=R(J)*SNS(J)
          VOL=R(J)*SNS(J)*SEW(I)
C
C *** CALCULATE COEFFICIENTS *******************************************
C
          DENN=0.5*(DEN(I,J)+DEN(I,J+1))
          DENS=0.5*(DEN(I,J)+DEN(I,J-1))
          DENE=0.5*(DEN(I,J)+DEN(I+1,J))
          DENW=0.5*(DEN(I,J)+DEN(I-1,J))
          AN(I,J)=DENN*AREAN*DV(I,J+1)
          AS(I,J)=DENS*AREAS*DV(I,J)
          AE(I,J)=DENE*AREAEW*DU(I+1,J)
          AW(I,J)=DENW*AREAEW*DU(I,J)
C
C *** CALCULATE SOURCE TERMS *******************************************
C
          CN=DENN*V(I,J+1)*AREAN
          CS=DENS*V(I,J)*AREAS
          CE=DENE*U(I+1,J)*AREAEW
          CW=DENW*U(I,J)*AREAEW
          IF (.NOT.STEADY) THEN
            CT=DEN(I,J)*VOL/DT
            CTO=DENO(I,J)*VOL/DT
          ENDIF
          SMP=CT-CTO+CN-CS+CE-CW
          SP(I,J)=0.0
          SU(I,J)=-SMP
C
C *** COMPUTE SUM OF ABSOLUTE MASS SOURCES *****************************
C
          RESORM=RESORM+ABS(SMP)
C
  110   CONTINUE
  100 CONTINUE
C
C-----------------------------------------------------------------------
CHAPTER 2 2 2 2 2 2 2 2  PROBLEM MODIFICATIONS  2 2 2 2 2 2 2 2 2 2 2 2
C-----------------------------------------------------------------------
C
      CALL MODP
C
C-----------------------------------------------------------------------
CHAPTER 3 3 3 3 3 3 3  FINAL ASSEMPLY OF COEFFICIENT  3 3 3 3 3 3 3 3 3
C-----------------------------------------------------------------------
C
      DO 300 I=2,NIM1
        DO 310 J=2,NJM1
          AP(I,J)=AN(I,J)+AS(I,J)+AE(I,J)+AW(I,J)-SP(I,J)
 310    CONTINUE
 300  CONTINUE
C
C-----------------------------------------------------------------------
CHAPTER 4 4 4 4 4 4 4  SOLUTION OF DIFFERENCE EQUATIONS  4 4 4 4 4 4 4 4
C-----------------------------------------------------------------------
C
      IF (ISLVP.EQ.1) THEN
        CALL ADISOL(2,2,NIM1,NJM1,PP,LITERP,SORP)
      ELSE IF(ISLVP.EQ.2) THEN
        CALL SIMSOL(2,2,NIM1,NJM1,PP,LITERP)
      ENDIF
C
C-----------------------------------------------------------------------
CHAPTER 5 5 5 5 5 5 5  CALCULATE VELOCITIES AND PRESSURE  5 5 5 5 5 5 5
C-----------------------------------------------------------------------
C
C *** VELOCITIES *******************************************************
C
      DO 500 I=2,NIM1
        DO 510 J=2,NJM1
          IF (I.NE.2) U(I,J)=U(I,J)+DU(I,J)*(PP(I-1,J)-PP(I,J))
          IF (J.NE.2) V(I,J)=V(I,J)+DV(I,J)*(PP(I,J-1)-PP(I,J))
  510   CONTINUE
  500 CONTINUE
C
C *** PRESSURES (WITH PROVISION FOR UNDER-RELAXATION) ******************
C
      PPREF=PP(IPREF,JPREF)
      DO 530 I=2,NIM1
        DO 540 J=2,NJM1
          P(I,J)=P(I,J)+URFP*(PP(I,J)-PPREF)
          PP(I,J)=0.0
  540   CONTINUE
  530 CONTINUE
C
C *** END OF CALCP SUBROUTINE ******************************************
C
      RETURN
      END
CC
C=======================================================================
C
C *** SUBROUTINE CALCT
C *** THIS SUBROUTINE SOLVES THE TEMPERATURE EQUATION
C
C=======================================================================
C

      SUBROUTINE CALCT
      INCLUDE 'common_l.inc'
C
C-----------------------------------------------------------------------
CHAPTER 1 1 1 1 1 1 1 1  ASSEMBLY OF COEFFICIENTS  1 1 1 1 1 1 1 1 1 1 1
C-----------------------------------------------------------------------
C
      CT=0.0
      CTO=0.0
      DO 100 I=2,NIM1
        DO 110 J=2,NJM1
C
C *** COMPUTE AREAS AND VOLUME *****************************************
C
          AREAN=RV(J+1)*SEW(I)
          AREAS=RV(J)*SEW(I)
          AREAEW=R(J)*SNS(J)
          VOL=R(J)*SNS(J)*SEW(I)
C
C *** CALCULATE CONVECTION COEFFICIENTS ********************************
C
          GN=0.5*(DEN(I,J)+DEN(I,J+1))*V(I,J+1)
          GS=0.5*(DEN(I,J)+DEN(I,J-1))*V(I,J)
          GE=0.5*(DEN(I,J)+DEN(I+1,J))*U(I+1,J)
          GW=0.5*(DEN(I,J)+DEN(I-1,J))*U(I,J)
          CN=GN*AREAN
          CS=GS*AREAS
          CE=GE*AREAEW
          CW=GW*AREAEW
C
C *** CALCULATE DIFFUSION COEFFICIENTS *********************************
C
          GAMN=0.5*(GAMH(I,J)+GAMH(I,J+1))
          GAMS=0.5*(GAMH(I,J)+GAMH(I,J-1))
          GAME=0.5*(GAMH(I,J)+GAMH(I+1,J))
          GAMW=0.5*(GAMH(I,J)+GAMH(I-1,J))
          DN=GAMN*AREAN/DYNP(J)
          DS=GAMS*AREAS/DYPS(J)
          DE=GAME*AREAEW/DXEP(I)
          DW=GAMW*AREAEW/DXPW(I)
C
C *** ASSEMBLE MAIN COEFFICIENTS ***************************************
C
          AN(I,J)=AMAX1(ABS(0.5*CN),DN)-0.5*CN
          AS(I,J)=AMAX1(ABS(0.5*CS),DS)+0.5*CS
          AE(I,J)=AMAX1(ABS(0.5*CE),DE)-0.5*CE
          AW(I,J)=AMAX1(ABS(0.5*CW),DW)+0.5*CW
C
C *** CALCULATE UNSTEADY TERM ******************************************
C
          IF (.NOT.STEADY) THEN
            CT=DEN(I,J)*VOL/DT
            CTO=DENO(I,J)*VOL/DT
          ENDIF
          SU(I,J)=CTO*TO(I,J)
          SP(I,J)=-CTO
C
C *** ADD CONTRIBUTION TO SOURCE TERMS DUE TO THE FACT *****************
C *** THAT CONTINUITY IS NOT EXACTLY SATISFIED *************************
C
          SMP=CT-CTO+CN-CS+CE-CW
          CP=AMAX1(0.0,SMP)
          SU(I,J)=SU(I,J)+CP*T(I,J)
          SP(I,J)=SP(I,J)-CP
C
  110   CONTINUE
  100 CONTINUE
C
C-----------------------------------------------------------------------
CHAPTER 2 2 2 2 2 2 2 2  PROBLEM MODIFICATIONS  2 2 2 2 2 2 2 2 2 2 2 2
C-----------------------------------------------------------------------
C
      CALL MODT
C
C-----------------------------------------------------------------------
CHAPTER 3 3 3 3 3 3 3  FINAL ASSEMPLY OF COEFFICIENT  3 3 3 3 3 3 3 3 3
C-----------------------------------------------------------------------
C
      RESORT=0.0
      DO 300 I=2,NIM1
        DO 310 J=2,NJM1
C
C *** CALCULATE AP COEFFICIENTS ****************************************
C
          AP(I,J)=AN(I,J)+AS(I,J)+AE(I,J)+AW(I,J)-SP(I,J)
C
C *** CALCULATE RESIDUALS **********************************************
C
          RESOR=AN(I,J)*T(I,J+1)+AS(I,J)*T(I,J-1)+AE(I,J)*T(I+1,J)
     &         +AW(I,J)*T(I-1,J)-AP(I,J)*T(I,J)+SU(I,J)
          VOL=R(J)*SEW(I)*SNS(J)
          SORVOL=GREAT*VOL
          IF (-SP(I,J).GT.0.5*SORVOL) RESOR=RESOR/SORVOL
          RESORT=RESORT+ABS(RESOR)
C
C *** UNDER-RELAXATION *************************************************
C
          AP(I,J)=AP(I,J)/URFT
          SU(I,J)=SU(I,J)+(1.0-URFT)*AP(I,J)*T(I,J)
C
  310   CONTINUE
  300 CONTINUE
C
C-----------------------------------------------------------------------
CHAPTER 4 4 4 4 4 4 4  SOLUTION OF DIFFERENCE EQUATIONS  4 4 4 4 4 4 4 4
C-----------------------------------------------------------------------
C
      IF (ISLVT.EQ.1) THEN
        CALL ADISOL(2,2,NIM1,NJM1,T,LITERT,SORT)
      ELSE IF (ISLVT.EQ.2) THEN
        CALL SIMSOL(2,2,NIM1,NJM1,T,LITERT)
      ENDIF
C
C *** END OF CALCT SUBROUTINE ******************************************
C
      RETURN
      END
CC
C=======================================================================
C
C *** SUBROUTINE CALCC
C *** THIS SUBROUTINE SOLVES THE CONCENTRATION EQUATION
C
C=======================================================================
C
      SUBROUTINE CALCC
      INCLUDE 'common_l.inc'
C
C-----------------------------------------------------------------------
CHAPTER 1 1 1 1 1 1 1 1  ASSEMBLY OF COEFFICIENTS  1 1 1 1 1 1 1 1 1 1 1
C-----------------------------------------------------------------------
C
      CT=0.0
      CTO=0.0
      DO 100 I=2,NIM1
        DO 110 J=2,NJM1
C
C *** COMPUTE AREAS AND VOLUME *****************************************
C
          AREAN=RV(J+1)*SEW(I)
          AREAS=RV(J)*SEW(I)
          AREAEW=R(J)*SNS(J)
          VOL=R(J)*SNS(J)*SEW(I)
C
C *** CALCULATE CONVECTION COEFFICIENTS ********************************
C
          GN=0.5*(DEN(I,J)+DEN(I,J+1))*V(I,J+1)
          GS=0.5*(DEN(I,J)+DEN(I,J-1))*V(I,J)
          GE=0.5*(DEN(I,J)+DEN(I+1,J))*U(I+1,J)
          GW=0.5*(DEN(I,J)+DEN(I-1,J))*U(I,J)
          CN=GN*AREAN
          CS=GS*AREAS
          CE=GE*AREAEW
          CW=GW*AREAEW
C
C *** CALCULATE DIFFUSION COEFFICIENTS *********************************
C
          GAMN=0.5*(GAMC(I,J)+GAMC(I,J+1))
          GAMS=0.5*(GAMC(I,J)+GAMC(I,J-1))
          GAME=0.5*(GAMC(I,J)+GAMC(I+1,J))
          GAMW=0.5*(GAMC(I,J)+GAMC(I-1,J))
          DN=GAMN*AREAN/DYNP(J)
          DS=GAMS*AREAS/DYPS(J)
          DE=GAME*AREAEW/DXEP(I)
          DW=GAMW*AREAEW/DXPW(I)
C
C *** ASSEMBLE MAIN COEFFICIENTS ***************************************
C
          AN(I,J)=AMAX1(ABS(0.5*CN),DN)-0.5*CN
          AS(I,J)=AMAX1(ABS(0.5*CS),DS)+0.5*CS
          AE(I,J)=AMAX1(ABS(0.5*CE),DE)-0.5*CE
          AW(I,J)=AMAX1(ABS(0.5*CW),DW)+0.5*CW
C
C *** CALCULATE UNSTEADY TERM ******************************************
C
          IF (.NOT.STEADY) THEN
            CT=DEN(I,J)*VOL/DT
            CTO=DENO(I,J)*VOL/DT
          ENDIF
          SU(I,J)=CTO*CO(I,J)
          SP(I,J)=-CTO
C
C *** ADD CONTRIBUTION TO SOURCE TERMS DUE TO THE FACT *****************
C *** THAT CONTINUITY IS NOT EXACTLY SATISFIED *************************
C
          SMP=CT-CTO+CN-CS+CE-CW
          CP=AMAX1(0.0,SMP)
          SU(I,J)=SU(I,J)+CP*C(I,J)
          SP(I,J)=SP(I,J)-CP
C
  110   CONTINUE
  100 CONTINUE
C
C-----------------------------------------------------------------------
CHAPTER 2 2 2 2 2 2 2 2  PROBLEM MODIFICATIONS  2 2 2 2 2 2 2 2 2 2 2 2
C-----------------------------------------------------------------------
C
      CALL MODC
C
C-----------------------------------------------------------------------
CHAPTER 3 3 3 3 3 3 3  FINAL ASSEMPLY OF COEFFICIENT  3 3 3 3 3 3 3 3 3
C-----------------------------------------------------------------------
C
      RESORC=0.0
      DO 300 I=2,NIM1
        DO 310 J=2,NJM1
C
C *** CALCULATE AP COEFFICIENTS ****************************************
C
          AP(I,J)=AN(I,J)+AS(I,J)+AE(I,J)+AW(I,J)-SP(I,J)
C
C *** CALCULATE RESIDUALS **********************************************
C
          RESOR=AN(I,J)*C(I,J+1)+AS(I,J)*C(I,J-1)+AE(I,J)*C(I+1,J)
     &         +AW(I,J)*C(I-1,J)-AP(I,J)*C(I,J)+SU(I,J)
          VOL=R(J)*SEW(I)*SNS(J)
          SORVOL=GREAT*VOL
          IF (-SP(I,J).GT.0.5*SORVOL) RESOR=RESOR/SORVOL
          RESORC=RESORC+ABS(RESOR)
C
C *** UNDER-RELAXATION *************************************************
C
          AP(I,J)=AP(I,J)/URFC
          SU(I,J)=SU(I,J)+(1.0-URFC)*AP(I,J)*C(I,J)
C
  310   CONTINUE
  300 CONTINUE
C
C-----------------------------------------------------------------------
CHAPTER 4 4 4 4 4 4 4  SOLUTION OF DIFFERENCE EQUATIONS  4 4 4 4 4 4 4 4
C-----------------------------------------------------------------------
C
      IF (ISLVC.EQ.1) THEN
        CALL ADISOL(2,2,NIM1,NJM1,C,LITERC,SORC)
      ELSE IF (ISLVT.EQ.2) THEN
        CALL SIMSOL(2,2,NIM1,NJM1,C,LITERC)
      ENDIF
C
C *** END OF CALCC SUBROUTINE ******************************************
C
      RETURN
      END
CC
C=======================================================================
C
C *** SUBROUTINE ADISOL
C *** THIS SUBROUTINE SOLVES THE SYSTEM OF DISCRETIZED EQUATIONS
C     BY EMPLOYING THE ADI (ALTERNATE DIRECTION IMPLICIT) ALGORITHM
C
C=======================================================================
C
      SUBROUTINE ADISOL(ISTART,JSTART,IEND,JEND,PHI,LITER,SOR)
      INCLUDE 'common_l.inc'
      DIMENSION PHI(II,JJ)
C
C *** START ITERATIONS *************************************************
C
      L=0
 100  L=L+1
C
C *** SOLVE EQUATIONS **************************************************
C
      CALL LISOLI(ISTART,JSTART,IEND,JEND,PHI)
      CALL LISOLJ(ISTART,JSTART,IEND,JEND,PHI)
C
C *** CALCULATE RESIDUALS **********************************************
C
      SUM=0.0
      DO 120 I=ISTART,IEND
        DO 110 J=JSTART,JEND
          RES=AN(I,J)*PHI(I,J+1)+AS(I,J)*PHI(I,J-1)+AE(I,J)*PHI(I+1,J)
     &       +AW(I,J)*PHI(I-1,J)-AP(I,J)*PHI(I,J)+SU(I,J)
          SUM=SUM+ABS(RES)
 110    CONTINUE
 120  CONTINUE
C
C *** TEST FOR CONVERGENCE *********************************************
C
      IF (SUM.GT.SOR.AND.L.LT.LITER) GOTO 100
C
C *** END OF SUBROUTINE ADISOL *****************************************
C
      RETURN
      END
CC
C=======================================================================
C
C *** SUBROUTINE LISOLI
C *** THIS SUBROUTINE IS A LINEAR EQUATION SOLVER USING A LINE BY LINE
C     TRI-DIAGONAL MATRIX ALGORITHM SWEEPING IN THE I-TH DIRECTION
C
C=======================================================================
C
      SUBROUTINE LISOLI(ISTART,JSTART,IEND,JEND,PHI)
      INCLUDE 'common_l.inc'
      DIMENSION PHI(II,JJ),AA(JJ),BB(JJ),CC(JJ),DD(JJ)
C
C *** COMMENCE W-E SWEEP ***********************************************
C
      JSTM1=JSTART-1
      AA(JSTM1)=0.0
      DO 100 I=ISTART,IEND
C
C *** COMMENCE S-N TRAVERSE ********************************************
C
        CC(JSTM1)=PHI(I,JSTM1)
        DO 110 J=JSTART,JEND
C
C *** ASSEMBLE TDMA COEFFICIENTS ***************************************
C
          AA(J)=AN(I,J)
          BB(J)=AS(I,J)
          CC(J)=AE(I,J)*PHI(I+1,J)+AW(I,J)*PHI(I-1,J)+SU(I,J)
          DD(J)=AP(I,J)
C
C *** CALCULATE COEFFICIENTS OF RECURRENCE FORMULA *********************
C
          TERM=1.0/(DD(J)-BB(J)*AA(J-1))
          AA(J)=AA(J)*TERM
          CC(J)=(CC(J)+BB(J)*CC(J-1))*TERM
  110   CONTINUE
C
C *** OBTAIN NEW PHI'S *************************************************
C
        DO 120 J=JEND,JSTART,-1
          PHI(I,J)=AA(J)*PHI(I,J+1)+CC(J)
  120   CONTINUE
C
C *** VISIT NEXT J-TH LINE *********************************************
C
  100 CONTINUE
C
C *** END OF SUBROUTINE LISOLI *****************************************
C
      RETURN
      END
CC
C=======================================================================
C
C *** SUBROUTINE LISOLJ
C *** THIS SUBROUTINE IS A LINEAR EQUATION SOLVER USING A LINE BY LINE
C     TRI-DIAGONAL MATRIX ALGORITHM SWEEPING IN THE J-TH DIRECTION
C
C=======================================================================
C
      SUBROUTINE LISOLJ(ISTART,JSTART,IEND,JEND,PHI)
      INCLUDE 'common_l.inc'
      DIMENSION PHI(II,JJ),AA(II),BB(II),CC(II),DD(II)
C
C *** COMMENCE S-N SWEEP ***********************************************
C
      ISTM1=ISTART-1
      AA(ISTM1)=0.0
      DO 100 J=JSTART,JEND
C
C *** COMMENCE W-E TRAVERSE ********************************************
C
        CC(ISTM1)=PHI(ISTM1,J)
        DO 110 I=ISTART,IEND
C
C *** ASSEMBLE TDMA COEFFICIENTS ***************************************
C
          AA(I)=AE(I,J)
          BB(I)=AW(I,J)
          CC(I)=AS(I,J)*PHI(I,J-1)+AN(I,J)*PHI(I,J+1)+SU(I,J)
          DD(I)=AP(I,J)
C
C *** CALCULATE COEFFICIENTS OF RECURRENCE FORMULA *********************
C
          TERM=1.0/(DD(I)-BB(I)*AA(I-1))
          AA(I)=AA(I)*TERM
          CC(I)=(CC(I)+BB(I)*CC(I-1))*TERM
  110   CONTINUE
C
C *** OBTAIN NEW PHI'S *************************************************
C
        DO 120 I=IEND,ISTART,-1
          PHI(I,J)=AA(I)*PHI(I+1,J)+CC(I)
  120   CONTINUE
C
C *** VISIT NEXT I-TH LINE *********************************************
C
  100 CONTINUE
C
C *** END OF SUBROUTINE LISOLJ *****************************************
C
      RETURN
      END
CC
C=======================================================================
C
C *** SUBROUTINE SIMSOL
C *** THIS SUBROUTINE SOLVES THE SYSTEM OF DISCRETIZED EQUATIONS
C     BY EMPLOYING THE SIM (STRONGLY IMPLICIT METHOD) ALGORITHM
C
C=======================================================================
C
      SUBROUTINE SIMSOL(ISTART,JSTART,IEND,JEND,PHI,ITER)
      INCLUDE 'common_l.inc'
      DIMENSION PHI(II,JJ),ANSX(II,JJ),AEX(II,JJ),SPX(II,JJ),
     &          PSIX(II,JJ),SUX(II,JJ),F(II,JJ)
      SAVE ANSX,AEX,SPX,PSIX,SUX,F
C
C-----------------------------------------------------------------------
CHAPTER 1 1 1 1 1 1 1 1 1 1 1  SET UP CONSTANTS  1 1 1 1 1 1 1 1 1 1 1 1
C-----------------------------------------------------------------------
C
      SMALL1=10E-35
      ALARG1=10E+10
      DO 100 I=1,NI
      DO 100 J=1,NJ
        SUX(I,J)=0.0
        ANSX(I,J)=0.0
        AEX(I,J)=0.0
        SPX(I,J)=0.0
        PSIX(I,J)=0.0
  100 CONTINUE
      IFLIP=0
C
C-----------------------------------------------------------------------
CHAPTER 2 2 2 2 2 2 2 2 2 2 2  SIM METHOD  2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
C-----------------------------------------------------------------------
C
      DO 200 L=1,ITER
C
C *** ALTERNATION OF SWEEP DIRECTION ( S->N , W->E ) *******************
C
  108 CONTINUE
      IF (IFLIP.EQ.0 ) THEN
      IFLIP=1
      DO 110 JK=JSTART,JEND
      DO 110 IJ=ISTART,IEND
C
C *** SUX USED FOR STORING THE RESIDUUM ********************************
C
        SUX(IJ,JK)=SU(IJ,JK)+
     &             AS(IJ,JK)*PHI(IJ,JK-1)+AW(IJ,JK)*PHI(IJ-1,JK)+
     &             AE(IJ,JK)*PHI(IJ+1,JK)+AN(IJ,JK)*PHI(IJ,JK+1)-
     &             AP(IJ,JK)*PHI(IJ,JK)
        IF (ABS(SUX(IJ,JK)).LT.SMALL1) SUX(IJ,JK)=SMALL1
        IF (ABS(SUX(IJ,JK)).GT.ALARG1) SUX(IJ,JK)=ALARG1
C
C *** A-COEFF. USED FOR STORING THE COEFF. OF THE TRIANGLE MATRICES ****
C
        ASX=-AS(IJ,JK)/(1.0+ALPHA*AEX(IJ,JK-1))
        AWX=-AW(IJ,JK)/(1.0+ALPHA*ANSX(IJ-1,JK))
        APX=AP(IJ,JK)+ASX*(ALPHA*AEX(IJ,JK-1)-ANSX(IJ,JK-1))
     &     +AWX*(ALPHA*ANSX(IJ-1,JK)-AEX(IJ-1,JK))
        RAP=1.0/APX
        AEX(IJ,JK) =(-AE(IJ,JK)-ALPHA*ASX*AEX(IJ,JK-1))*RAP
        ANSX(IJ,JK)=(-AN(IJ,JK)-ALPHA*AWX*ANSX(IJ-1,JK))*RAP
C
C *** SP USED FOR STORING UPPER TRIANGLE MATRIX TIMES SOLUTION VECTOR **
C
        SPX(IJ,JK)=(SUX(IJ,JK)-ASX*SPX(IJ,JK-1)-AWX*SPX(IJ-1,JK))*RAP
C
  110 CONTINUE
C
C *** CALCULATION OF SOLUTION VECTOR ***********************************
C
      DO 120 JK=JEND,JSTART,-1
      DO 120 IJ=IEND,ISTART,-1
        PSIX(IJ,JK)=SPX(IJ,JK)-AEX(IJ,JK)*PSIX(IJ+1,JK)-ANSX(IJ,JK)*
     &              PSIX(IJ,JK+1)
  120 CONTINUE
C
      ELSE
C
C *** ALTERNATION OF SWEEP DIRECTION ( N->S , W->E ) *******************
C
      IFLIP=0
      DO 150 JK=JEND,JSTART,-1
      DO 150 IJ=ISTART,IEND
C
C *** SUX USED FOR STORING THE RESIDUUM ********************************
C
        SUX(IJ,JK)=SU(IJ,JK)+
     &            AS(IJ,JK)*PHI(IJ,JK-1)+AW(IJ,JK)*PHI(IJ-1,JK)+
     &            AE(IJ,JK)*PHI(IJ+1,JK)+AN(IJ,JK)*PHI(IJ,JK+1)-
     &            AP(IJ,JK)*PHI(IJ,JK)
        IF (ABS(SUX(IJ,JK)).LT.SMALL1) SUX(IJ,JK)=SMALL1
        IF (ABS(SUX(IJ,JK)).GT.ALARG1) SUX(IJ,JK)=ALARG1
C
C *** A-COEFF. USED FOR STORING THE COEFF. OF THE TRIANGLE *************
C
        ANX=-AN(IJ,JK)/(1.0+ALPHA*AEX(IJ,JK+1))
        AWX=-AW(IJ,JK)/(1.0+ALPHA*ANSX(IJ-1,JK))
        APX=AP(IJ,JK)+ANX*(ALPHA*AEX(IJ,JK+1)-ANSX(IJ,JK+1))
     &      +AWX*(ALPHA*ANSX(IJ-1,JK)-AEX(IJ-1,JK))
        RAP=1.0/APX
        AEX(IJ,JK) =(-AE(IJ,JK)-ALPHA*ANX*AEX(IJ,JK+1))*RAP
        ANSX(IJ,JK)=(-AS(IJ,JK)-ALPHA*AWX*ANSX(IJ-1,JK))*RAP
C
C *** SP USED FOR STORING UPPER TRIANGLE MATRIX TIMES SOLUTION VECTOR **
C
        SPX(IJ,JK)=(SUX(IJ,JK)-ANX*SPX(IJ,JK+1)-AWX*SPX(IJ-1,JK))*RAP
C
  150 CONTINUE
C
C *** CALCULATION OF SOLUTION VECTOR ***********************************
C
      DO 160 JK=JSTART,JEND
      DO 160 IJ=IEND,ISTART,-1
        PSIX(IJ,JK)=SPX(IJ,JK)-AEX(IJ,JK)*PSIX(IJ+1,JK)-ANSX(IJ,JK)*
     &             PSIX(IJ,JK-1)
  160 CONTINUE
C
      END IF
C
C-----------------------------------------------------------------------
CHAPTER 3 3 3 3 3 3 3 3  RELAXATION FACTOR AND ALPHA  3 3 3 3 3 3 3 3 3
C-----------------------------------------------------------------------
C
C *** CALCULATE OPTIMUM RELAXATION FACTOR ******************************
C
      DO 170 JK=JSTART,JEND
      DO 170 IJ=ISTART,IEND
        F(IJ,JK)=-AS(IJ,JK)*PSIX(IJ,JK-1)-AW(IJ,JK)*PSIX(IJ-1,JK)-
     &            AE(IJ,JK)*PSIX(IJ+1,JK)-AN(IJ,JK)*PSIX(IJ,JK+1)+
     &            AP(IJ,JK)*PSIX(IJ,JK)
  170 CONTINUE
      SUM1=0.0
      SUM2=0.0
      DO 180 JK=JSTART,JEND
      DO 180 IJ=ISTART,IEND
        SUM1=SUM1+SUX(IJ,JK)*F(IJ,JK)
        SUM2=SUM2+F(IJ,JK)*F(IJ,JK)
  180 CONTINUE
      TEST=SIGN(1.0,-SUM2)
      TEST=AMAX1(0.0,TEST)
      OME=(SUM1+TEST)/(SUM2+TEST)
      TEST=SIGN(1.0,OME)
      TEST=AMAX1(0.0,TEST)
      OME=OME*TEST
      ALPHA=ALPHA*TEST
      OME=AMIN1(1.0,OME)
C
C *** CALCULATE NEW PHI ************************************************
C
      DO 190 JK=1,JEND
      DO 190 IJ=1,IEND
        PHI(IJ,JK)=PHI(IJ,JK)+OME*PSIX(IJ,JK)
  190 CONTINUE
C
C *** ALTER SWEEP DIRECTION OR RETURN **********************************
C
  200 CONTINUE
C
C *** END OF SIM SUBROUTINE ********************************************
C
      RETURN
      END


CC
C=======================================================================
C
C *** SUBROUTINE GRID1
C *** THIS SUBROUTINE CREATES A UNIFORM COMPUTATIONAL GRID
C
C=======================================================================
C
      SUBROUTINE GRID1(X,II,XSTART,XLENGTH,ISTART,IEND)
C
      DIMENSION X(II)
      NI=IEND-ISTART+1
      DX=XLENGTH/FLOAT(NI)
      X(ISTART)=XSTART+0.5*DX
      DO 100 I=ISTART+1,IEND
         X(I)=X(I-1)+DX
 100  CONTINUE
      RETURN
      END
CC
C=======================================================================
C
C *** SUBROUTINE GRID2
C *** THIS SUBROUTINE CREATES A NON-UNIFORM COMPUTATIONAL GRID
C *** GEOMETRICAL PROGRESION FORMULA (CONDENSES NEAR A WALL)
C
C=======================================================================
C
      SUBROUTINE GRID2(X,II,XSTART,XLENGTH,ISTART,IEND,FEXPX)
C
      DIMENSION X(II)
C
C *** RANGES FOR PARAMETER FEXPX
C
      IF (FEXPX .LT. 0.1) FEXPX=0.1
      IF (FEXPX .GT. 5.0) FEXPX=5.0
C
C *** CALCULATE GRID POSITIONS
C
      NI=IEND-ISTART+1
      DX=XLENGTH/FLOAT(NI)
      IF (FEXPX.NE.1.0)
     &DX=2*XLENGTH*(1.0-FEXPX)/(1.0+FEXPX-FEXPX**NI-FEXPX**(NI+1))
      X(ISTART)=XSTART+0.5*DX
      DX=FEXPX*DX
      DO 100 I=ISTART+1,IEND
         X(I)=X(I-1)+DX
         DX=FEXPX*DX
 100  CONTINUE
      RETURN
      END
CC
C=======================================================================
C
C *** SUBROUTINE GRID3
C *** THIS SUBROUTINE CREATES A NON-UNIFORM COMPUTATIONAL GRID
C *** ANDERSON FORMULA for a=0.5 
C *** CONDENSES NEAR XSTART AND XSTART+XLENGTH
C
C=======================================================================
C
      SUBROUTINE GRID3(X,II,XSTART,XLENGTH,ISTART,IEND,BETA)
C
      DIMENSION X(II)
C
C *** RANGES FOR PARAMETER BETA
C
      IF (BETA .LT. 1.000001) BETA=1.000001
      IF (BETA .GT. 10.00000) BETA=10.00000
C
C *** CALCULATE GRID POSITIONS
C
      NI=IEND-ISTART+1
      DO 100 I=1,NI
         XBAR=FLOAT(I)/FLOAT(NI+1)    ! NI-2 points internal
         EKV=2.0*(XBAR-0.5)
         BP1BM1=((BETA+1)/(BETA-1))**EKV
         DX=0.5*XLENGTH*((BETA+1)*BP1BM1-BETA+1)/(1+BP1BM1)
         X(ISTART+I-1)=XSTART+DX
  100 CONTINUE
      RETURN
      END
CC
C=======================================================================
C
C *** SUBROUTINE GRID4
C *** THIS SUBROUTINE CREATES A NON-UNIFORM COMPUTATIONAL GRID
C *** ANDERSON FORMULA - CONDENSES NEAR INTERNAL POINT XINT
C *** [ XINT ] IS THE DISTANCE OF INTERNAL POINT FROM [ XSTART ]
C
C=======================================================================
C
      SUBROUTINE GRID4(X,II,XSTART,XLENGTH,XINT,ISTART,IEND,TAU)
C
      DIMENSION X(II)
C
C *** RANGES FOR PARAMETER TAU
C
      IF (TAU .LT. 0.000001) TAU=0.000001
      IF (TAU .GT. 25.00000) TAU=25.00000
C
C *** CALCULATE GRID POSITIONS
C
      NI=IEND-ISTART+1
      XCH=XINT/XLENGTH
      BETA=0.5*ALOG((1+(EXP(TAU)-1)*XCH)/(1+(EXP(-TAU)-1)*XCH))/TAU
      DO 100 I=1,NI
         XBAR=FLOAT(I)/FLOAT(NI+1)    ! NI-2 points are internal
         DX=XINT*(1+SINH(TAU*(XBAR-BETA))/SINH(TAU*BETA))
         X(ISTART+I-1)=XSTART+DX
  100 CONTINUE
      RETURN
      END

CC*************************************************************************
CC
CC  TOOLBOX LIBRARY v.1.0 (May 1995)
CC
CC  Program unit   : SUBROUTINE CHRBLN
CC  Purpose        : Position of last non-blank character in a string
CC  Author         : Athanasios Nenes
CC
CC  ======================= ARGUMENTS / USAGE =============================
CC
CC  STR        is the CHARACTER variable containing the string examined
CC  IBLK       is a INTEGER variable containing the position of last non 
CC             blank character. If string is all spaces (ie '   '), then
CC             the value returned is 1.
CC
CC  EXAMPLE:
CC             STR = 'TEST1.DAT     '
CC             CALL CHRBLN (STR, IBLK)
CC          
CC  after execution of this code segment, "IBLK" has the value "9", which
CC  is the position of the last non-blank character of "STR".
CC
CC***********************************************************************
CC
      SUBROUTINE CHRBLN (STR, IBLK)
CC
CC***********************************************************************
      CHARACTER*(*) STR
C
      IBLK = 1                       ! Substring pointer (default=1)
      ILEN = LEN(STR)                ! Length of string 
      DO 10 i=ILEN,1,-1              
         IF (STR(i:i).NE.' ' .AND. STR(i:i).NE.CHAR(0)) THEN
            IBLK = i
            RETURN
         ENDIF
10    CONTINUE
      RETURN
C
      END




CC*************************************************************************
CC
CC  TOOLBOX LIBRARY v.1.0 (May 1995)
CC
CC  Program unit   : SUBROUTINE SHFTRGHT
CC  Purpose        : RIGHT-JUSTIFICATION FUNCTION ON A STRING
CC  Author         : Athanasios Nenes
CC
CC  ======================= ARGUMENTS / USAGE =============================
CC
CC  STRING     is the CHARACTER variable with the string to be justified
CC
CC  EXAMPLE:
CC             STRING    = 'AAAA    '
CC             CALL SHFTRGHT (STRING)
CC          
CC  after execution of this code segment, STRING contains the value
CC  '    AAAA'.
CC
CC*************************************************************************
CC
      SUBROUTINE SHFTRGHT (CHR)
CC
CC***********************************************************************
      CHARACTER CHR*(*)
C
      I1  = LEN(CHR)             ! Total length of string
      CALL CHRBLN(CHR,I2)        ! Position of last non-blank character
      IF (I2.EQ.I1) RETURN
C
      DO 10 I=I2,1,-1            ! Shift characters
         CHR(I1+I-I2:I1+I-I2) = CHR(I:I)
         CHR(I:I) = ' '
10    CONTINUE
      RETURN
C
      END




CC*************************************************************************
CC
CC  TOOLBOX LIBRARY v.1.0 (May 1995)
CC
CC  Program unit   : SUBROUTINE RPLSTR
CC  Purpose        : REPLACE CHARACTERS OCCURING IN A STRING
CC  Author         : Athanasios Nenes
CC
CC  ======================= ARGUMENTS / USAGE =============================
CC
CC  STRING     is the CHARACTER variable with the string to be edited
CC  OLD        is the old character which is to be replaced
CC  NEW        is the new character which OLD is to be replaced with
CC  IERR       is 0 if everything went well, is 1 if 'NEW' contains 'OLD'.
CC             In this case, this is invalid, and no change is done.
CC
CC  EXAMPLE:
CC             STRING    = 'AAAA'
CC             OLD       = 'A'
CC             NEW       = 'B' 
CC             CALL RPLSTR (STRING, OLD, NEW)
CC          
CC  after execution of this code segment, STRING contains the value
CC  'BBBB'.
CC
CC*************************************************************************
CC
      SUBROUTINE RPLSTR (STRING, OLD, NEW, IERR)
CC
CC***********************************************************************
      CHARACTER STRING*(*), OLD*(*), NEW*(*)
C
C *** INITIALIZE ********************************************************
C
      ILO = LEN(OLD)
C
C *** CHECK AND SEE IF 'NEW' CONTAINS 'OLD', WHICH CANNOT ***************
C      
      IP = INDEX(NEW,OLD)
      IF (IP.NE.0) THEN
         IERR = 1
         RETURN
      ELSE
         IERR = 0
      ENDIF
C
C *** PROCEED WITH REPLACING *******************************************
C      
10    IP = INDEX(STRING,OLD)      ! SEE IF 'OLD' EXISTS IN 'STRING'
      IF (IP.EQ.0) RETURN         ! 'OLD' DOES NOT EXIST ; RETURN
      STRING(IP:IP+ILO-1) = NEW   ! REPLACE SUBSTRING 'OLD' WITH 'NEW'
      GOTO 10                     ! GO FOR NEW OCCURANCE OF 'OLD'
C
      END
        

CC*************************************************************************
CC
CC  TOOLBOX LIBRARY v.1.0 (May 1995)
CC
CC  Program unit   : SUBROUTINE Pushend 
CC  Purpose        : Positions the pointer of a sequential file at its end
CC                   Simulates the ACCESS='APPEND' clause of a F77L OPEN
CC                   statement with Standard Fortran commands.
CC
CC  ======================= ARGUMENTS / USAGE =============================
CC
CC  Iunit      is a INTEGER variable, the file unit which the file is 
CC             connected to.
CC
CC  EXAMPLE:
CC             CALL PUSHEND (10)
CC          
CC  after execution of this code segment, the pointer of unit 10 is 
CC  pushed to its end.
CC
CC***********************************************************************
CC
      SUBROUTINE Pushend (Iunit)
CC
CC***********************************************************************
C
      LOGICAL OPNED
C
C *** INQUIRE IF Iunit CONNECTED TO FILE ********************************
C
      INQUIRE (UNIT=Iunit, OPENED=OPNED)
      IF (.NOT.OPNED) GOTO 25
C
C *** Iunit CONNECTED, PUSH POINTER TO END ******************************
C
10    READ (Iunit,'()', ERR=20)
      GOTO 10
C
C *** RETURN POINT ******************************************************
C
20    BACKSPACE (Iunit)
25    RETURN
      END



CC*************************************************************************
CC
CC  TOOLBOX LIBRARY v.1.0 (May 1995)
CC
CC  Program unit   : SUBROUTINE APPENDEXT
CC  Purpose        : Fix extension in file name string
CC
CC  ======================= ARGUMENTS / USAGE =============================
CC
CC  Filename   is the CHARACTER variable with the file name
CC  Defext     is the CHARACTER variable with extension (including '.',
CC             ex. '.DAT')
CC  Overwrite  is a LOGICAL value, .TRUE. overwrites any existing extension
CC             in "Filename" with "Defext", .FALSE. puts "Defext" only if 
CC             there is no extension in "Filename".
CC
CC  EXAMPLE:
CC             FILENAME1 = 'TEST.DAT'
CC             FILENAME2 = 'TEST.DAT'
CC             CALL APPENDEXT (FILENAME1, '.TXT', .FALSE.)
CC             CALL APPENDEXT (FILENAME2, '.TXT', .TRUE. )
CC          
CC  after execution of this code segment, "FILENAME1" has the value 
CC  'TEST.DAT', while "FILENAME2" has the value 'TEST.TXT'
CC
CC***********************************************************************
CC
      SUBROUTINE Appendext (Filename, Defext, Overwrite)
CC
CC***********************************************************************
      CHARACTER*(*) Filename, Defext
      LOGICAL       Overwrite
C
      CALL CHRBLN (Filename, Iend)
      IF (Filename(1:1).EQ.' ' .AND. Iend.EQ.1) RETURN  ! Filename empty
      Idot = INDEX (Filename, '.')                      ! Append extension ?
      IF (Idot.EQ.0) Filename = Filename(1:Iend)//Defext
      IF (Overwrite .AND. Idot.NE.0)
     &              Filename = Filename(:Idot-1)//Defext
      RETURN
      END





C=======================================================================
C
C *** FUNCTION VPRES
C *** THIS FUNCTION CALCULATES SATURATED WATER VAPOUR PRESSURE AS A 
C     FUNCTION OF TEMPERATURE. VALID FOR TEMPERATURES BETWEEN -50 AND 
C     50 C.
C
C ======================== ARGUMENTS / USAGE ===========================
C
C  INPUT:
C     [T] 
C     REAL variable.
C     Ambient temperature expressed in Kelvin. 
C
C  OUTPUT:
C     [VPRES] 
C     REAL variable.  
C     Saturated vapor pressure expressed in mbar.
C
C=======================================================================
C
      REAL FUNCTION VPRES (T)
      REAL A(0:6), T
      DATA A/6.107799610E+0, 4.436518521E-1, 1.428945805E-2,
     &       2.650648471E-4, 3.031240396E-6, 2.034080948E-8,
     &       6.136820929E-11/
C
C Calculate polynomial (without exponentiation).
C
CCC      TTEMP = MAX(MIN(T-273.0,50.0), -50.0) ! Assure correct range
      TTEMP = T-273.0
      VPRES = A(6)*TTEMP
      DO I=5,1,-1
         VPRES = (VPRES + A(I))*TTEMP
      ENDDO
      VPRES = VPRES + A(0)
C
ccc      VPRES = 18.3036d0 - 3816.44d0/(DBLE(T)-46.13d0)
ccc      VPRES = EXP(VPRES)
ccc      VPRES = VPRES*(10.0d0/7.6d0)
C
C End of FUNCTION VPRES
C
      RETURN
      END

