      PARAMETER( MXX=29999, MCOL=171 )
      IMPLICIT REAL*4(A-H,O-Z)
C ***************************************************************C
C                                                                C
C   PROGRAM:  MAP5                                               C
C                                                                C
C   DATE:  8-24-90                                               C
C   PROGRAMMER:  J.  FASTOOK                                     C
C                                                                C
C   FUNCTION:                                                    C
C           THIS IS A PROGRAM TO MODEL THE FLOW OF A GLACIER     C
C           WITH MATERIAL PROPERTIES READ FOR EACH NODAL POINT   C
C           AND THE AVERAGE USED FOR THE ELEMENT.                C
C           IT CALCULATES AN ELEMENT MATRIX ALA BECKER, 2-D      C
C           PROGRAM WITH LINEAR SHAPE FUNCTIONS.                 C
C           DOES TIME DEPENDENT CASE USING A LUMPED CAPACITANCE  C
C           MATRIX AND A BACKWARD DIFFERENCE SCHEME. ALLOWING    C
C           FOR VERY QUICK OPERATION, HOWEVER, MATERIAL PROPS    C
C           ARE HELD CONSTANT THRUOUT THE TIME ITERATIONS        C
C           SO THAT THE FINAL EQUILIBRIUM STATE DIFFERS SLIGHLTY C
C           FROM ONE WHICH RECALCULATED THE STIFFNESS MAATRIX AT C
C           EACH TIME STEP                                       C
C           THIS IS HANDLED BY OUTPUTING A COMPLETE DATA SET AT  C
C           THE END OF THE TIME STEP, (INDEED AT THE END OF      C
C           EVERY PROGRAM RUN SO THAT BY TIME STEPPING ONE STEP  C
C           AT A TIME AND USING THE NEW DATA SET IT IS AS IF THE C
C           COMPLETE MATRIX WERE REDEVELOPED AT EACH STEP.       C
C                                                                C
C ***************************************************************C
C     PROGRAM FOR STEADY AND UNSTEADY STATE FLOW ANALYSIS
C     USING FINITE ELEMENTS
C RUN BY MAP5 EXEC A:
C       FI * CLEAR
C       FI 30 DISK INPUT&1 HEAD B
C       FI 31 DISK INPUT&1 GRID B
C       FI 32 DISK INPUT&1 DIFF B
C       FI 33 DISK INPUT&1 TIME B
C       FI 34 DISK OUT&2 TIME B
C       FI 10 DISK OUTLINE DATA B (LRECL 80 RECFM F
C       FI 7 DISK OUT1&2 DATA B (LRECL 130 RECFM F
C       FI 12 DISK OUT2&2 DATA B (LRECL 130 RECFM F
C       FI 11 DISK OUT3&2 DATA B (LRECL 130
C       FI 21 DISK OUT4&2 DATA B (LRECL 130
C       FI 20 DISK LINEM DATA B (LRECL 80
C       FI 26 DISK OUTPUT&2 DATA B (LRECL 130
C       FI 13 DISK MAT&2 DATA B
C       FI 18 DISK VOL&2 DATA B
C       FI 17 DISK VOLT&2 DATA B
C       GL TXT  VSF2FORT CMSLIB IGLLIB ESSL
C       LOAD MAP5 (START NOMAP
      CHARACTER*80 HED,SCRTCH
      COMMON /SYMARG/ NUMNP,MBAND,A(MCOL,MXX),Q(MXX)
      DIMENSION AMASS(9), B(MXX), X(MXX), Y(MXX), D(MXX),
     &          KODE(MXX), CONST(MXX), LM(5), IX(3), E(3,3),
     &          KX(MXX,4), P(5), S(5,5), DD(5), IDT(MXX),
     &          ADOT(MXX), BDROCK(MXX), FLOWA(MXX), SLDGB(MXX),
     &          PSURF(MXX), PPSURF(MXX), FRACT(MXX),
     &          CNEW(MXX), QHOLD(MXX), HTICE(MXX), THICK(MXX),
     &          HFIT(MXX), IBFLUX(MXX,2), BFLUX(MXX), DEPB(MXX),
     &          UNDEPB(MXX), SLOPE(MXX), SLOPN(MXX), KZ(MXX)
      DIMENSION TTIME(2000),VVOL(2000)
      DIMENSION NTYPE(MXX), NNODE(MXX),
     &          DPSIX(9), DPSIY(9), DXDS(2,2), DSDX(2,2),
     &          PSI(4), DPSI(4,2), CNST(MXX),
     &          XY(2,4), XI(2,9), ETA(2,9), W(2,9),
     &          AADOT(MXX), AFRACT(MXX), AFLOWA(MXX),
     &          ABDRCK(MXX), ASLDGB(MXX)
      REAL*4 T(MXX)
      DATA HFIT /MXX*0.0/
C     CALL XUFLOW(0)
C FOLLOWING IS DEFAULT SNOWLINE ELEVATION AT POLE AND GRADIENT
      AMASS(7)=-3000.
      AMASS(8)=.56E-3
      AMASS(8)=.001
      AMASS(9)=-14.
C **** FOLLOWING IS SEALEVEL REFERENCED TO PRESENT=0.
      NTSTEP=0
      SEALEV=0.
C **** FOLLOWING SETS RATE OF CONVERGENCE, UP TO 5 WORKS WELL
      CONV=1.
      TIME=0.0
      RHOW=1.092
      PG = 0.089866
      NUMCOL=0
      NUMLEV=0
      NCOL=MCOL
C
C     INITIALIZE INTEGRATION POINTS AND WEIGHTS
C     GAUSSIAN QUADRATURE OF ORDER THREE QUADRILATERALS
      CALL GAUSINIT(XI,ETA,W)
C
      DO 10 I=1,4
        LM(I)=0
10    CONTINUE
C
C **** FOLLOWING READN FOR NEW SPLIT DATA SETS
      CALL READN(MXX, IDEP, HED, NUMNP, NUMEL, NUMGBC, NDT,
     &           INTER, DT, KODE, X, Y, HTICE, ADOT, FRACT,
     &           PSURF, RHOI, BDROCK, UNDEPB, FLOWA, SLDGB,
     &           THICK, KX, CONST, IBFLUX, BFLUX, QHOLD,
     &           NTYPE, NNODE, NCOL, DEPB, AADOT, AFRACT,
     &           ABDRCK, PPSURF, AFLOWA, ASLDGB, IDT, AMASS)
      DO 90 I=1,NUMNP
        Q(I)=HTICE(I)
        T(I)=HTICE(I)
90    CONTINUE
C
C SET LINEARIZATION CONSTANT USING INITIAL CONFIGURATION
      CALL NCONST(MXX, IDEP, X, Y, KX, NTYPE, NUMNP, NUMEL,
     &            AFRACT, ASLDGB, LM, AFLOWA, BDROCK, DEPB,
     &            UNDEPB, KODE, PG, Q, CNEW, SLOPE, RHOI)
C
      WRITE(*,*) 'TIME STEP=',DT
      WRITE(*,*) 'INPUT 1 TO CALL ADJUST, 0 TO BYPASS'
      READ(*,*) IADJ
      IF(IADJ.EQ.1) THEN
C
C CALCULATE SLOPES IN CASE NEEDED BY ADJUST
        CALL NODESL(MXX, NUMNP, NUMEL, KX, SLOPE, SLOPN)
C
C ENTER INTERACTIVE DATA SET MANIPULATOR
        CALL ADJUST(HED, NUMNP, NUMEL, X, Y, HTICE, ADOT, FRACT,
     &         PSURF, BDROCK, DEPB, FLOWA, SLDGB, THICK, KX, CONST,
     &              NNODE, KODE, HFIT, NUMCOL, NUMLEV, NUMGBC, NDT,
     &              INTER, DT, IBFLUX, BFLUX, MXX, IDT, SLOPN, AMASS,
     &              TIME, NTSTEP,TTIME,VVOL)
C
      ENDIF
      DO 91 I=1,NUMNP
        Q(I)=HTICE(I)
91    CONTINUE
C
C LOAD NODAL PROPERTIES INTO ELEMENT PROPERTIES
      CALL ELPROP(MXX, NUMEL, NTYPE, KX, ADOT, AADOT, FRACT, AFRACT,
     &            BDROCK, ABDRCK, PSURF, PPSURF, FLOWA, AFLOWA,
     &            SLDGB, ASLDGB)
C
C
C
C
C
C                      *** MAIN LOOP ***
C
C
C
215   CONTINUE
C
C       FOLLOWING CALCULATES NODE SLOPE FROM ELEMENT SLOPE
        CALL NODESL(MXX, NUMNP, NUMEL, KX, SLOPE, SLOPN)
C
C       FOLLOWING ADJUSTS ADOT FOR FITTEED ACCUMULATION
        DO 101 I=1,NUMNP
          IF(IDT(I).GT.0) ADOT(I)=AFUNCT(TIME, IDT(I), AMASS,
     &                                 REAL(Q(I)), BDROCK(I),
     &                                 SLOPN(I),X(I),Y(I))
101     CONTINUE
C
C LOADS NEW NODAL MATERIAL PROPERTIES INTO ELEMENT MATERIAL PROPERTIES
        CALL ELPROP(MXX, NUMEL, NTYPE, KX, ADOT, AADOT, FRACT, AFRACT,
     &            BDROCK, ABDRCK, PSURF, PPSURF, FLOWA, AFLOWA,
     &            SLDGB, ASLDGB)
C
C CALCULATES VOLUMES (FLOTATION AND TOTAL) AND AREAL EXTENT
        CALL VOLUME(MXX, TIME, NUMNP, NUMEL, X, Y, KX, Q, BDROCK,
     &             DEPB, SEALEV, RHOI, RHOW, VOL, AMASS)
C
      NTSTEP=NTSTEP+1
      TTIME(NTSTEP)=TIME
      VVOL(NTSTEP)=VOL*1.E-15
C
C       FORM STIFFNESS,CAPACITANCE AND LOAD
      CALL FORMA(MXX,MCOL,X,Y,KX,NTYPE,NUMNP,NUMEL,NCOL,ETA,XI,W,CONST,
     &     ADOT,FRACT,BDROCK,PSURF,RHOI,FLOWA,SLDGB,
     &     T,KODE,NUMGBC,IBFLUX,BFLUX,
     &        MBAND,LM,AADOT,AFRACT,ABDRCK,AFLOWA,ASLDGB,D,B,A)
C
C       FORM EFFECTIVE CONDUCTIVITY MATRIX FOR TIME INCREMENT
          DT2=1.0/DT
          DO 320 N=1,NUMNP
            IF (KODE(N).EQ.0) THEN
              IF (D(N).NE.0.) THEN
                D(N)=DT2*D(N)
                A(MBAND+1,N)=A(MBAND+1,N)+D(N)
              ENDIF
            ENDIF
320       CONTINUE
          LL=0
C
C
C
C         LOOP ON NUMBER OF TIME STEPS *******************************
          CALL GRSTRT(500,1)
          CALL WINDOW(0.,100.,0.,100.)
          CALL DRAWOUT
          DO 450 L=1,NDT
C
C           CALCULATE EFFECTIVE LOAD MATRIX
            DO 360 I=1,NUMNP
              Q(I)=B(I)
              IF (KODE(I).EQ.0) Q(I)=B(I)+D(I)*T(I)
360         CONTINUE
C
            CALL ASYMSL(1)
            CALL ASYMSL(2)
            TIME=TIME+DT
            WRITE(*,*) 'TIME=',TIME
            LL=LL+1
            HMAX=-1.E30
            DIFF=0.
            NDIFF=0
            DO 5511 JK=1,NUMNP
              IF(KODE(JK).NE.1) THEN
                DIFF=DIFF+(Q(JK)-PSURF(JK))**2
                NDIFF=NDIFF+1
              ENDIF
              IF(Q(JK).GT.HMAX) THEN
                HMAX=Q(JK)
                NMAX=JK
              ENDIF
              IF (Q(JK).LE.UNDEPB(JK)) Q(JK)=UNDEPB(JK)
C
C ****        TO DEAL WITH BED BELOW SEA LEVEL AND
C ****        SURFACE BELOW FLOTATION HEIGHT
              IF(UNDEPB(JK).LE.0.) THEN
                FLOT=(1.-RHOW/RHOI)*UNDEPB(JK)
C               IF(Q(JK).LE.FLOT) Q(JK)=FLOT
                IF(Q(JK).LE.FLOT) Q(JK)=0.
                IF(BDROCK(JK).LE.-9999.) Q(JK)=0.
              ENDIF
C ****        END FLOTATION PART ****
5511        CONTINUE
C
C THIS IS SPECIAL TO DEAL WITH TERRYS TEMP DEP STUFF, PLEASE REMOVE
C       DO 5512 JK=1,NUMNP
C         IF(Q(JK).GT.1250.) THEN
C           FRACT(JK)=FRACT(JK)+.05
C           IF(FRACT(JK).GT.1.) FRACT(JK)=1.
C         ELSE
C           FRACT(JK)=FRACT(JK)-.05
C           IF(FRACT(JK).LT.0.) FRACT(JK)=0.
C         ENDIF
C 5512  CONTINUE
C       CALL ELPROP(MXX, NUMEL, NTYPE, KX, ADOT, AADOT, FRACT, AFRACT,
C      &            BDROCK, ABDRCK, PSURF, PPSURF, FLOWA, AFLOWA,
C      &            SLDGB, ASLDGB)
C END TERRYS SPECIAL STUFF
C
C OUTPUT TIME STEP INFO TO SCREEN
            WRITE(*,*) 'MXX SURF=',HMAX,' AT NODE',NMAX,' DIFF=',
     &              SQRT(DIFF/REAL(NDIFF))
C
C OBTAIN NEW LINEARIZATION CONSTANT FROM LATEST SOLUTION
            CALL NCONST(MXX, IDEP, X, Y, KX, NTYPE, NUMNP, NUMEL,
     &             AFRACT, ASLDGB, LM, AFLOWA, BDROCK, DEPB, UNDEPB,
     &             KODE, PG, Q, CNEW, SLOPE, RHOI)
C
C DERIVE SLOPE IN CASE NEEDED BY MASS BALANCE PARAMETERIZATION
            CALL NODESL(MXX, NUMNP, NUMEL, KX, SLOPE, SLOPN)
            DO 102 I=1,NUMNP
              IF(IDT(I).GT.0) ADOT(I)=AFUNCT(TIME, IDT(I), AMASS,
     &                            REAL(Q(I)), BDROCK(I), SLOPN(I),
     &                            X(I),Y(I))
102         CONTINUE
C
C LOAD NEW NODE PROPERTIES INTO ELEMENT PROPERTIES
            CALL ELPROP(MXX, NUMEL, NTYPE, KX, ADOT, AADOT, FRACT,
     &             AFRACT, BDROCK, ABDRCK, PSURF, PPSURF, FLOWA,
     &             AFLOWA, SLDGB, ASLDGB)
C
C CALCULATE VOLUMES (FLOTATION AND TOTAL) AND AREA
            CALL VOLUME(MXX, TIME, NUMNP, NUMEL, X, Y, KX, Q, BDROCK,
     &            DEPB, SEALEV, RHOI, RHOW, VOL, AMASS)
C
      NTSTEP=NTSTEP+1
      TTIME(NTSTEP)=TIME
      VVOL(NTSTEP)=VOL*1.E-15
            DO 380 I=1,NUMNP
              T(I)=Q(I)
380         CONTINUE
C
C LOAD NEW LINEARIZATION CONSTANT INTO OLD
            DO 381 I=1,NUMEL
C             CONST(I)=CNEW(I)
              CONST(I) = (CONV*CONST(I) + CNEW(I))/(CONV+1)
381         CONTINUE
C
C           FORM STIFF,CAP,LOAD FOR NEXT PASS
      CALL FORMA(MXX,MCOL,X,Y,KX,NTYPE,NUMNP,NUMEL,NCOL,ETA,XI,W,CONST,
     &     ADOT,FRACT,BDROCK,PSURF,RHOI,FLOWA,SLDGB,
     &     T,KODE,NUMGBC,IBFLUX,BFLUX,
     &        MBAND,LM,AADOT,AFRACT,ABDRCK,AFLOWA,ASLDGB,D,B,A)
C
            DO 3201 N=1,NUMNP
              IF (KODE(N).EQ.0) THEN
                IF (D(N).NE.0.0) THEN
                  D(N)=DT2*D(N)
                  A(MBAND+1,N)=A(MBAND+1,N)+D(N)
                ENDIF
              ENDIF
3201        CONTINUE
C
C PRINT OUT SOLUTION FOR APPROPRIATE TIME STEPS
            IF(LL.GE.INTER) THEN
              DO 731 N=1,NUMNP
                HTICE(N)=Q(N)
731           CONTINUE
              WRITE(SCRTCH,*) 'TIME=',TIME
              WRITE(34) SCRTCH
              WRITE(34) (HTICE(I),I=1,NUMNP)
              WRITE(34) (ADOT(I),I=1,NUMNP)
              WRITE(34) (DEPB(I),I=1,NUMNP)
              WRITE(34) (CONST(I),I=1,NUMEL)
              LL=0
            ENDIF
          CALL PLOTSOL(NUMNP,X,Y,HTICE,DEPB,KODE)
450       CONTINUE
          CALL GRSTOP1
C
C     END OF TIME STEP SECTION ***************************************
C
C
C
      HMAX=-1.E30
      DIFF=0.
      NDIFF=0
      DO 730 N=1,NUMNP
        IF(KODE(N).NE.1) THEN
          DIFF=DIFF+(Q(N)-PSURF(N))**2
          NDIFF=NDIFF+1
        ENDIF
        HTICE(N)=Q(N)
        IF(HTICE(N).GT.HMAX) THEN
          HMAX=HTICE(N)
          NMAX=N
        ENDIF
        IF (HTICE(N).LE.UNDEPB(N)) HTICE(N)=UNDEPB(N)
C ****  TO DEAL WITH BED BELOW SEA LEVEL
C ****  AND SURFACE BELOW FLOTATION HEIGHT
        IF(UNDEPB(N).LE.0.) THEN
          FLOT=(1.-RHOW/RHOI)*UNDEPB(N)
C         IF(HTICE(N).LE.FLOT) HTICE(N)=FLOT
          IF(HTICE(N).LE.FLOT) HTICE(N)=0.
          IF(BDROCK(N).LE.-9999.) HTICE(N)=0.
        ENDIF
C ****  END FLOTATION PART ****
730   CONTINUE
C
C OUTPUT TIME STEP INFO TO SCREEN
      WRITE(*,*) 'MXX SURF=', HMAX, ' AT NODE',NMAX,' DIFF=',
     &            SQRT(DIFF/REAL(NDIFF))
C
C     FORMAT STATEMENTS
C
1001  FORMAT(I6,I4,1P2E12.5,0PF10.2,F7.2,F9.2,F10.3,F10.1,F10.5,2F10.5,
     &       I5,F10.3)
1002  FORMAT(5I5,1PE17.10)
1003  FORMAT(I10,4F10.0)
1007  FORMAT(2I5,E13.6)
2002  FORMAT(2I2,1PE10.3,2E10.3,/,10X,6E10.3)
2003  FORMAT('  N  I  J  K  L',/,
     &10X,'        COND-  SPECIFIC   HEAT      DENSITY     VERT HORI'
     &,/,
     &10X,'      UCTIVITY   HEAT   GENERATED             VELOCITY',/,
     &10X,'       CAL/M/D/YR       CAL/M**3/YR             M/YR',/,
     &10X,'               CAL/KG/D             KG/M**3')
2004  FORMAT(5I3,7E14.4)
2005  FORMAT(' TIME=', E12.5,/,(I6,E14.6,I6,E14.6,I6,E14.6,I6,E14.6,
     & I6,E14.6,I6,E14.6))
2006  FORMAT('    I    J              H    TEMPERATURE')
2007  FORMAT(2I5,2E15.6)
2009  FORMAT('  M ',14X,'K',14X,'C',14X,'D',14X,'Q',/,(I4,4E15.6))
2010  FORMAT(' AXISYMMETRIC SOLID BODY')
2011  FORMAT('TWO DIMENSIONAL PLANE BODY ')
2020  FORMAT(' CARD NO. ',I4, ' OUT OF ORDER')
2021  FORMAT(' BAD CARD NO. ',I4)
3005  FORMAT(2F10.0,G13.6)
C
C
      WRITE(*,8989)
8989  FORMAT (' END OF RUN')
C     END OF RUN, ENTER ADJUST TO GO ROUND AGAIN.
      WRITE(*,*) 'INPUT 1 TO CALL ADJUST -9 TO SKIP AND END'
      READ(*,*,END=999) IADJ
      IF(IADJ.EQ.-9) GOTO 999
C
C ENTER INTERACTIVE DATA SET MANIPULATOR
      CALL ADJUST(HED, NUMNP, NUMEL, X, Y, HTICE, ADOT, FRACT, PSURF,
     &        BDROCK, DEPB, FLOWA, SLDGB, THICK, KX, CONST, NNODE, KODE,
     &            HFIT, NUMCOL, NUMLEV, NUMGBC, NDT, INTER, DT,
     &            IBFLUX, BFLUX, MXX, IDT, SLOPN, AMASS, TIME,
     &            NTSTEP, TTIME, VVOL)
C
      DO 92 I=1,NUMNP
        Q(I)=HTICE(I)
92    CONTINUE
C
C LOAD NODAL PROPERTIES INTO ELEMENT PROPERTIES
      CALL ELPROP(MXX, NUMEL, NTYPE, KX, ADOT, AADOT, FRACT, AFRACT,
     &            BDROCK, ABDRCK, PSURF, PPSURF, FLOWA, AFLOWA, SLDGB,
     &            ASLDGB)
C
      WRITE(*,*) 'INPUT 1 TO CONTINUE WITH NEW SET '
      READ(*,*,END=999) IADJ
      ITER=1
C     GOTO BEGINNING OF MAIN LOOP
C
C
C
      GOTO 215
C
C
C
999   CONTINUE
C
C VERBOSE WRITER OFF
C     CALL WRITER(HED, NUMNP, NUMEL, X, Y, HTICE, ADOT, FRACT, PSURF,
C    &            RHOI, BDROCK, FLOWA, SLDGB, THICK, KX, CONST, NNODE,
C    &            KODE, HFIT, NUMCOL, NUMLEV, NUMGBC, NDT, INTER, DT,
C    &            IBFLUX, BFLUX, AADOT, AFRACT, AFLOWA, ABDRCK,
C    &            ASLDGB)
C
      WRITE(18,2000) -99999.,2.,0,HED
2000  FORMAT(10X,G13.6,2X,G13.6,I13,/,A80)
      STOP
      END
C
      SUBROUTINE SHAPE(NTYPE,XI,ET,PSI,DPSI)
C ELEMENT SHAPE FUNCTIONS AND DERIVATIVES AT LOCAL COORDINATES (XI,ET)
      IMPLICIT REAL*4(A-H,O-Z)
      DIMENSION PSI(4),DPSI(4,2)
      IF(NTYPE.EQ.2) GOTO 100
      PSI(1)=.25*(1.-XI)*(1.-ET)
      PSI(2)=.25*(1.+XI)*(1.-ET)
      PSI(3)=.25*(1.+XI)*(1.+ET)
      PSI(4)=.25*(1.-XI)*(1.+ET)
      DPSI(1,2)=-.25*(1.-XI)
      DPSI(2,2)=-.25*(1.+XI)
      DPSI(3,2)=.25*(1.+XI)
      DPSI(4,2)=.25*(1.-XI)
      DPSI(1,1)=-.25*(1.-ET)
      DPSI(2,1)=.25*(1.-ET)
      DPSI(3,1)=.25*(1.+ET)
      DPSI(4,1)=-.25*(1.+ET)
      RETURN
100   CONTINUE
      PSI(1)=1.-XI-ET
      PSI(2)=XI
      PSI(3)=ET
      DPSI(1,2)=-1.
      DPSI(2,2)=0.
      DPSI(3,2)=1.
      DPSI(1,1)=-1.
      DPSI(2,1)=1.
      DPSI(3,1)=0.
      RETURN
      END
C
      SUBROUTINE NCONST(MXX,IDEP,X,Y,KX,NTYPE,NUMNP,NUMEL,AFRACT,ASLDGB,
     &    LM,AFLOWA,BDROCK,DEPB,UNDEPB,KODE,PG,Q,CNEW,SLOPE,RHOI)
      IMPLICIT REAL*4(A-H,O-Z)
C CALCULATES LINEARIZATION CONSTANT FROM CURRENT SOLUTION
      DIMENSION KODE(MXX),AFRACT(MXX),ASLDGB(MXX),AFLOWA(MXX)
      DIMENSION BDROCK(MXX),DEPB(MXX),UNDEPB(MXX),SLOPE(MXX)
      DIMENSION KX(MXX,4),X(MXX),Y(MXX),NTYPE(MXX)
      DIMENSION LM(5),CNEW(MXX)
      DIMENSION DPSIX(9),DPSIY(9),DXDS(2,2),DSDX(2,2)
      DIMENSION PSI(4),DPSI(4,2)
      DIMENSION XY(2,4)
      REAL*4 Q(MXX)
      RHOR=4.0
C
C IDEP=1 DOES UNDEPRESSION AND DEPRESSION, IDEP=0 RIGID BED
C
      IF(IDEP.EQ.1) THEN
        FDEP1=RHOR/(RHOR-RHOI)
        FDEP2=RHOI/(RHOI-RHOR)
      ELSE
        FDEP1=1.
        FDEP2=0.
      ENDIF
      DO 600 J = 1,NUMEL
        IF(NTYPE(J).EQ.1) THEN
          NODEN=4
          CENTX=0.0
          CENTY=0.0
        ELSE
          NODEN=3
          CENTX=1.D0/3.D0
          CENTY=1.D0/3.D0
        ENDIF
        SUMHH=0.
        SUMX=0.
        SUMY=0.
        DO 560 I = 1,NODEN
          LM(I) = KX(J,I)
560     CONTINUE
        I=LM(1)
        JJ=LM(2)
        K=LM(3)
        L=LM(4)
        XY(1,1)=X(I)
        XY(1,2)=X(JJ)
        XY(1,3)=X(K)
        IF(NTYPE(J).EQ.1) XY(1,4)=X(L)
        XY(2,1)=Y(I)
        XY(2,2)=Y(JJ)
        XY(2,3)=Y(K)
        IF(NTYPE(J).EQ.1) XY(2,4)=Y(L)
        CALL SHAPE(NTYPE(J),CENTX,CENTY,PSI,DPSI)
C
C CALCULATE DXDS...EQUATION (5.3.6)
C
        DO 565 I=1,2
        DO 565 L=1,2
          DXDS(I,L)=0.0
          DO 565 K=1,NODEN
            DXDS(I,L)=DXDS(I,L)+DPSI(K,L)*XY(I,K)
565     CONTINUE
C
C CALCULATE DSDX...EQUATION (5.2.7)
C
        DETJ=(DXDS(1,1)*DXDS(2,2)-DXDS(1,2)*DXDS(2,1))
        IF (DETJ.LE.0.0) THEN
          WRITE(12,5544) J,DETJ
          WRITE(*,5544) J,DETJ
          WRITE(12,5545) (JJ,XY(1,JJ),XY(2,JJ),JJ=1,4)
          WRITE(*,5545) (JJ,XY(1,JJ),XY(2,JJ),JJ=1,4)
5545      FORMAT(1X,I5,1X,1PE10.3,E10.3)
          STOP
5544      FORMAT(' BAD JACOBIAN',I5,1PE10.3,/,1X,8E10.3)
        ENDIF
        DSDX(1,1)=DXDS(2,2)/DETJ
        DSDX(2,2)=DXDS(1,1)/DETJ
        DSDX(1,2)=-DXDS(1,2)/DETJ
        DSDX(2,1)=-DXDS(2,1)/DETJ
C
C CALCULATE D(PSI)/DX...EQUATION (5.3.5)
C
        DO 570 I=1,NODEN
          DPSIX(I)=DPSI(I,1)*DSDX(1,1)+DPSI(I,2)*DSDX(2,1)
          DPSIY(I)=DPSI(I,1)*DSDX(1,2)+DPSI(I,2)*DSDX(2,2)
570     CONTINUE
        DO 580 I = 1,NODEN
          SUMX = SUMX + Q(LM(I))*DPSIX(I)
          SUMY = SUMY + Q(LM(I))*DPSIY(I)
          IF(BDROCK(LM(I)).LE.-9999.) THEN
            DEPB(LM(I))=0.
          ELSE
            DEPB(LM(I))=FDEP2*Q(LM(I))+FDEP1*UNDEPB(LM(I))
          ENDIF
          THIK=(Q(LM(I))-UNDEPB(LM(I)))*FDEP1
          IF(THIK.GT.0.) SUMHH=SUMHH+THIK
580     CONTINUE
C
        DELH = SUMX**2 + SUMY**2
        DELH = SQRT(DELH)
        SLOPE(J)=DELH
        HH = SUMHH/FLOAT(NODEN)
        TERM1 = AFRACT(J)*((PG/ASLDGB(J))**2)*(HH**3)*DELH
        TERM2=(1.-AFRACT(J))*(.2)*((PG/AFLOWA(J))**3)*(HH**5)*(DELH**2)
        CNEW(J) = TERM1 + TERM2
  600 CONTINUE
      RETURN
      END
C
      SUBROUTINE ELPROP(MXX,NUMEL,NTYPE,KX,ADOT,AADOT,FRACT,AFRACT,
     &BDROCK,ABDRCK,PSURF,PPSURF,FLOWA,AFLOWA,SLDGB,ASLDGB)
      IMPLICIT REAL*4(A-H,O-Z)
C LOAD NODAL PROPERTIES INTO ELEMENT PROPERTIES
      DIMENSION NTYPE(MXX),ADOT(MXX),AADOT(MXX),FRACT(MXX),AFRACT(MXX),
     &BDROCK(MXX),ABDRCK(MXX),PSURF(MXX),PPSURF(MXX),FLOWA(MXX),
     &AFLOWA(MXX),SLDGB(MXX),ASLDGB(MXX),LM(4),KX(MXX,4)
      DO 100 N = 1,NUMEL
        IF(NTYPE(N).EQ.1) THEN
          NODEN=4
          NINT=9
        ELSE
          NODEN=3
          NINT=4
        ENDIF
        AAADOT=0.
        AAFRCT=0.
        AABDRK=0.
        APSURF=0.
        AAFLOW=0.
        AASLDG=0.
        ADNSTY=0.
        DO 90 I=1,4
          LM(I)=KX(N,I)
          AAADOT=AAADOT+ADOT(LM(I))
          AAFRCT=AAFRCT+FRACT(LM(I))
          AABDRK=AABDRK+BDROCK(LM(I))
          APSURF=APSURF+PSURF(LM(I))
          AAFLOW=AAFLOW+FLOWA(LM(I))
          AASLDG=AASLDG+SLDGB(LM(I))
90      CONTINUE
        DENOM=1./FLOAT(NODEN)
        AAADOT=AAADOT*DENOM
        AADOT(N)=AAADOT
        AAFRCT=AAFRCT*DENOM
        AFRACT(N)=AAFRCT
        AABDRK=AABDRK*DENOM
        ABDRCK(N)=AABDRK
        ADNSTY=ADNSTY*DENOM
        PPSURF(N)=APSURF
        AAFLOW=AAFLOW*DENOM
        AFLOWA(N)=AAFLOW
        AASLDG=AASLDG*DENOM
        ASLDGB(N)=AASLDG
100   CONTINUE
      RETURN
      END
C
C
      SUBROUTINE FORMA(MXX,MCOL,X,Y,KX,NTYPE,NUMNP,NUMEL,NCOL,ETA,XI,W,
     &     CONST,ADOT,FRACT,BDROCK,PSURF,RHOI,FLOWA,SLDGB,
     &     T,KODE,NUMGBC,IBFLUX,BFLUX,
     &        MBAND,LM,AADOT,AFRACT,ABDRCK,AFLOWA,ASLDGB,D,B,A)
      IMPLICIT REAL*4(A-H,O-Z)
      DIMENSION AADOT(MXX),AFRACT(MXX),AFLOWA(MXX)
      DIMENSION ABDRCK(MXX),ASLDGB(MXX),T(MXX),KODE(MXX),IBFLUX(MXX,2)
      DIMENSION BFLUX(MXX)
      DIMENSION ADOT(MXX),FRACT(MXX),PSURF(MXX),FLOWA(MXX)
      DIMENSION BDROCK(MXX),SLDGB(MXX)
      DIMENSION KX(MXX,4),X(MXX),Y(MXX),NTYPE(MXX),D(MXX),B(MXX)
      DIMENSION A(MCOL,MXX),LM(MXX),CONST(MXX)
      DIMENSION P(5),S(5,5),DD(5)
      DIMENSION DPSIX(9),DPSIY(9),DXDS(2,2),DSDX(2,2)
      DIMENSION PSI(4),DPSI(4,2)
      DIMENSION XY(2,4),XI(2,9),ETA(2,9),W(2,9)
C **********************************************************************
C     FORM CONDUCTIVITY MATRIX FOR COMPLETE BODY
C **********************************************************************
      DO 110 I=1,NUMNP
      D(I)=0.0
      B(I)=0.0
      DO 110 J=1,NCOL
  110 A(J,I)=0.0
      MBAND=(NCOL-1)/2
C     WRITE(12,2003)
C
      DO 160 N=1,NUMEL
      IF(NTYPE(N).EQ.1) THEN
        NODEN=4
        NINT=9
      ELSE
        NODEN=3
        NINT=4
      ENDIF
      DO 125 I=1,4
      LM(I)=KX(N,I)
125   CONTINUE
C
C
C     2. FORM ELEMENT CONDUCTIVITY MATRIX
C
      DO 130 I=1,4
      DD(I)=0.0
      P(I)=0.0
      DO 130 J=1,4
  130 S(I,J)=0.0
C
      I=LM(1)
      J=LM(2)
      K=LM(3)
      L=LM(4)
      XY(1,1)=X(I)
      XY(1,2)=X(J)
      XY(1,3)=X(K)
C     IF(NTYPE(N).EQ.1) XY(1,4)=X(L)
      XY(2,1)=Y(I)
      XY(2,2)=Y(J)
      XY(2,3)=Y(K)
C     IF(NTYPE(N).EQ.1) XY(2,4)=Y(L)
      IF(NODEN.EQ.4) THEN
        XY(1,4)=X(L)
        XY(2,4)=Y(L)
      ENDIF
C
C FORM ELEMENT MATRIX AND VECTORS
C
C BEGIN INTEGRATION POINT LOOP
      DO 139 L=1,NINT
      CALL SHAPE(NTYPE(N),XI(NTYPE(N),L),ETA(NTYPE(N),L),PSI,DPSI)
C CALCULATE DXDS...EQUATION (5.3.6)
      DO 133 I=1,2
      DO 133 J=1,2
      DXDS(I,J)=0.0
      DO 133 K=1,4
 133   DXDS(I,J)=DXDS(I,J)+DPSI(K,J)*XY(I,K)
C CALCULATE DSDX...EQUATION (5.2.7)
      DETJ=(DXDS(1,1)*DXDS(2,2)-DXDS(1,2)*DXDS(2,1))
      IF(DETJ.LE.0.0) GOTO 141
      DENOM=1./DETJ
      DSDX(1,1)=DXDS(2,2)*DENOM
      DSDX(2,2)=DXDS(1,1)*DENOM
      DSDX(1,2)=-DXDS(1,2)*DENOM
      DSDX(2,1)=-DXDS(2,1)*DENOM
C CALCULATE D(PSI)/DX...EQUATION (5.3.5)
      DO 136 I=1,4
      DPSIX(I)=DPSI(I,1)*DSDX(1,1)+DPSI(I,2)*DSDX(2,1)
      DPSIY(I)=DPSI(I,1)*DSDX(1,2)+DPSI(I,2)*DSDX(2,2)
 136  CONTINUE
C ACCUMULATE INTEGRATION POINT VALUES OF INTEGRALS
      FAC=DETJ*W(NTYPE(N),L)
      DO 138 I=1,4
      DD(I)=DD(I)+PSI(I)*FAC
C THIS IS LUMPED CAPACITANCE MATRIX
      P(I)=P(I)+AADOT(N)*PSI(I)*FAC
      DO 138 J=1,4
C
      IF(CONST(N).GT.1.E-30) THEN
        TERM1=CONST(N)*(DPSIX(I)*DPSIX(J)+DPSIY(I)*DPSIY(J))
        S(I,J)=S(I,J)+TERM1*FAC
      ENDIF
 138  CONTINUE
 139  CONTINUE
      GOTO 142
141   WRITE(12,1100) DETJ,((XY(MM,NN),NN=1,4),MM=1,2)
1100  FORMAT(' BAD JACOBIAN at 161',G13.6,/,4G13.6,/,4G13.6)
      STOP
142   CONTINUE
C
C     3. ADD ELEMENT CONDUCTIVITY TO COMPLETE CONDUCTIVITY MATRIX
C
      DO 155 L=1,4
      I=LM(L)
      D(I)=D(I)+DD(L)
C THIS (D) IS LUMPED CAPACITANCE MATRIX
      B(I)=B(I)+P(L)
      DO 155 M=1,4
      J=LM(M)-I+1
C OPTIONAL CHECK OF BANDWIDTH, CAN BE REMOVED FOR MORE SPEED
C     IF(J.GT.NCOL) THEN
C     WRITE(12,3000)
C     STOP
C     ENDIF
3000  FORMAT('INCREASE BANDWIDTH AND DIMENSION OF A()')
      A(J+MBAND,I)=A(J+MBAND,I)+S(L,M)
      ZZZ=J+MBAND
  155 CONTINUE
  160 CONTINUE
C
C
C     BOUNDARY CONDITIONS
C
C
      IF (NUMGBC) 2102,2102,1901
1901  CONTINUE
      DO 2101 N=1,NUMGBC
      I=IBFLUX(N,1)
      J=IBFLUX(N,2)
      TEMP=BFLUX(N)
      XL=SQRT((X(J)-X(I))**2+(Y(J)-Y(I))**2)
      TEMP=XL*TEMP*.5
      B(I)=B(I)+TEMP
      B(J)=B(J)+TEMP
2101  CONTINUE
2102  CONTINUE
C REPEAT FROM ADJUST HERE
C BRACNCH 215 BACK HERE, NOW MOVED
C
C     2. TEMPERATURE BOUNDARY CONDITIONS
C
C  *** THIS IS THE LOOP THAT CAUSES PROBLEMS IN VECTOR COMPILE
      DO 300 N=1,NUMNP
      IF(KODE(N)) 225,300,225
225   CONTINUE
      DO 275 M=MAX0(N-MBAND,1),MIN0(N+MBAND,NUMNP)
      K=N+MBAND-M+1
      B(M)=B(M)-A(K,M)*T(N)
275   CONTINUE
  300 CONTINUE
      DO 1300 N=1,NUMNP
      IF(KODE(N)) 1225,1300,1225
1225  CONTINUE
      DO 1250 M=1,MBAND
      A(M,N)=0.0
1250  A(MBAND+1+M,N)=0.0
      DO 1275 M=MAX0(N-MBAND,1),MIN0(N+MBAND,NUMNP)
      K=N+MBAND-M+1
1275  A(K,M)=0.0
1300  CONTINUE
      DO 1400 N=1,NUMNP
      IF(KODE(N)) 1325,1400,1325
1325  A(MBAND+1,N)=1.0
      B(N)=T(N)
1400  CONTINUE
      RETURN
      END
      SUBROUTINE VOLUME(MXX,TIME,NUMNP,NUMEL,X,Y,KX,Q,BDROCK,DEPB,
     &SEALEV,RHOI,RHOW,VOL,AMASS)
      IMPLICIT REAL*4(A-H,O-Z)
C CALCULATES VOLUMES (FLOTATION AND TOTAL) AND AREA
      DIMENSION BDROCK(MXX),KX(MXX,4),X(MXX),Y(MXX),LM(4)
      DIMENSION DEPB(MXX)
      REAL*4 Q(MXX)
      DIMENSION AMASS(9)
      VOL=0.0
      VOL1=0.
      AREATOT=0.0
      DO 100 I=1,NUMEL
        SUMH=0.
        SUMT=0.
        DO 90 J=1,4
        LM(J)=KX(I,J)
        IF(BDROCK(LM(J)).LE.-9999.) GOTO 100
        IF(DEPB(LM(J)).LT.SEALEV) THEN
          FLOT=(1.-RHOW/RHOI)*(DEPB(LM(J))-SEALEV)
        ELSE
          FLOT=DEPB(LM(J))
        ENDIF
        HEIGHT=Q(LM(J))-FLOT
        THICK=Q(LM(J))-DEPB(LM(J))
        IF(HEIGHT.GT.1.) SUMH=SUMH+HEIGHT
        IF(THICK.GT.1.) SUMT=SUMT+THICK
90    CONTINUE
      HEIGHT=SUMH*.25
      THICK=SUMT*.25
      AREA=0.5*((X(LM(2))-X(LM(1)))*(Y(LM(3))-Y(LM(2)))-
     &          (X(LM(3))-X(LM(2)))*(Y(LM(2))-Y(LM(1)))+
     &          (X(LM(4))-X(LM(3)))*(Y(LM(1))-Y(LM(4)))-
     &          (X(LM(1))-X(LM(4)))*(Y(LM(4))-Y(LM(3))))
C IF THICKNESS LT 1 METER, NEGLECT
      IF(HEIGHT.GT.1.) THEN
        VOL=VOL+AREA*HEIGHT
        VOL1=VOL1+AREA*THICK
        AREATOT=AREATOT+AREA
      ENDIF
100   CONTINUE
      IF(AREATOT.GT.0.) THEN
        AVGHGT=VOL/AREATOT
      ELSE
        AVGHGT=0.
      ENDIF
      WRITE(*,1000) VOL*1.E-15,AREATOT*1.E-12,AVGHGT
      WRITE(*,1000) VOL1*1.E-15,AMASS(7),AMASS(9)
C     WRITE(18,2000) TIME,VOL*1.E-15
      WRITE(18,2000) TIME,VOL1*1.E-15
2000  FORMAT(10X,G13.6,2X,G13.6)
1000  FORMAT(' V=',G13.6,' MEG KM**3, A=',G13.6,' MEG KM**2, <H>=',G13.6
     &)
      WRITE(17,1004) TIME,VOL1*1E-15,VOL*1E-15,AREATOT*1E-12,AMASS(7),
     &              AMASS(9),AVGHGT
1004  FORMAT(1X,1P,G13.6,0P,3F7.3,F6.0,F6.1,F7.1)               
1001  FORMAT(' TIME=',G13.6)
1002  FORMAT(' TOTAL VOLUME=',G13.6,' MEG KM**3')
      RETURN
      END
C
      SUBROUTINE NODESL(MXX,NUMNP,NUMEL,KX,SLOPE,SLOPN)
C CALCULATES SLOES FOR ACCUMULATION PARAMETERIZATIONS
      PARAMETER(MMXX=29999)
      IMPLICIT REAL*4(A-H,O-Z)
      DIMENSION KX(MXX,4),SLOPE(MXX),SLOPN(MXX),ICOUNT(MMXX)
      DO 100 I=1,NUMNP
        SLOPN(I)=0.
        ICOUNT(I)=0
100   CONTINUE
      DO 200 I=1,NUMEL
        DO 200 J=1,4
          SLOPN(KX(I,J))=SLOPN(KX(I,J))+SLOPE(I)
          ICOUNT(KX(I,J))=ICOUNT(KX(I,J))+1
200   CONTINUE
      DO 300 I=1,NUMNP
        IF(ICOUNT(I).GT.0) THEN
          SLOPN(I)=SLOPN(I)/REAL(ICOUNT(I))
        ENDIF
300   CONTINUE
      RETURN
      END
      SUBROUTINE GAUSINIT(XI,ETA,W)
      IMPLICIT REAL*4(A-H,O-Z)
      DIMENSION XI(2,9), ETA(2,9), W(2,9)
C     GAUSSIAN QUADRATURE OF ORDER THREE QUADRILATERALS
      XI(1,1)=-SQRT(3./5.)
      XI(1,2)=0.
      XI(1,3)=-XI(1,1)
      XI(1,4)=XI(1,1)
      XI(1,5)=0.
      XI(1,6)=XI(1,3)
      XI(1,7)=XI(1,1)
      XI(1,8)=0.
      XI(1,9)=XI(1,3)
      ETA(1,1)=XI(1,1)
      ETA(1,2)=XI(1,1)
      ETA(1,3)=XI(1,1)
      ETA(1,4)=0.
      ETA(1,5)=0.
      ETA(1,6)=0.
      ETA(1,7)=XI(1,3)
      ETA(1,8)=XI(1,3)
      ETA(1,9)=XI(1,3)
      W(1,1)=25./81.
      W(1,2)=40./81.
      W(1,3)=W(1,1)
      W(1,4)=W(1,2)
      W(1,5)=64./81.
      W(1,6)=W(1,2)
      W(1,7)=W(1,1)
      W(1,8)=W(1,2)
      W(1,9)=W(1,1)
C   GAUSSIAN QUADRATURE OF ORDER THREE TRIANGLES
      XI(2,1)=1./3.
      XI(2,2)=2./15.
      XI(2,3)=XI(2,2)
      XI(2,4)=11./15.
      ETA(2,1)=XI(2,1)
      ETA(2,2)=XI(2,4)
      ETA(2,3)=XI(2,2)
      ETA(2,4)=ETA(2,3)
      W(2,1)=-27./96.
      W(2,2)=25./96.
      W(2,3)=W(2,2)
      W(2,4)=W(2,2)
      END
      SUBROUTINE ASYMSL(KKK)
      IMPLICIT REAL*4(A-H,O-Z)
      PARAMETER(MXX=29999,MCOL=171)
      COMMON /SYMARG/ NEQ,MBAND,A(MCOL,MXX),Q(MXX)
      DATA TOLER /1.E-20/
      NBAND=MBAND+1
      NCOLS=2*NBAND-1
      GO TO (1000,2000),KKK
1000  CONTINUE
      KMIN=NBAND+1
      DO 50 N=1,NEQ
      IF(A(NBAND,N).EQ.0.0) GOTO 60
      IF(A(NBAND,N).EQ.1.0) GOTO 20
      C=1./A(NBAND,N)
      DO 10 K=KMIN,NCOLS
      IF(A(K,N).EQ.0.0) GOTO 10
      A(K,N)=C*A(K,N)
10    CONTINUE
20    CONTINUE
      DO 40 L=2,NBAND
      JJ=NBAND-L+1
      I=N+L-1
      IF(I.GT.NEQ) GOTO 40
      IF(A(JJ,I).EQ.0.0) GOTO 40
      KI=NBAND+2-L
      KF=NCOLS+1-L
      J=NBAND
      DO 30 K=KI,KF
      J=J+1
C     IF(ABS(A(JJ,I)).LT.TOLER .OR. ABS(A(J,N)).LT.TOLER) GOTO 30
      A(K,I)=A(K,I)-A(JJ,I)*A(J,N)
C     PRINT *,A(JJ,I),A(J,N)
C     IF(ABS(A(JJ,I)).LT.1.E-20) WRITE(7,*) JJ,I,A(JJ,I)
C     IF(ABS(A(J,N)).LT.1.E-20) WRITE(7,*) J,N,A(J,N)
30    CONTINUE
40    CONTINUE
50    CONTINUE
      RETURN
60    CONTINUE
C     WRITE(12,70) N,A(NBAND,N)
      WRITE(*,70) N,A(NBAND,N)
70    FORMAT(1H1,5X,' SET OF EQUATIONS MAY BE SINGULAR',//,5X,
     &'DIAGONAL TERM OF EQUATION',I5,13H IS EQUAL TO ,E15.8)
      STOP
2000  CONTINUE
      DO 130 N=1,NEQ
      IF(A(NBAND,N).EQ.0.0) GOTO 160
      IF(A(NBAND,N).EQ.1.0) GOTO 110
      Q(N)=Q(N)/A(NBAND,N)
110   CONTINUE
      DO 120 L=2,NBAND
      JJ=NBAND-L+1
      I=N+L-1
      IF(I.GT.NEQ) GOTO 120
      IF(A(JJ,I).EQ.0.0) GOTO 120
      Q(I)=Q(I)-A(JJ,I)*Q(N)
120   CONTINUE
130   CONTINUE
C
C BACK SUBSTITUTION
C
      LL=NBAND+1
      DO 150 M=1,NEQ
      N=NEQ+1-M
      DO 140 L=LL,NCOLS
      IF(A(L,N).EQ.0.0) GOTO 140
      K=N+L-NBAND
      Q(N)=Q(N)-A(L,N)*Q(K)
140   CONTINUE
150   CONTINUE
      RETURN
160   CONTINUE
C     WRITE(12,170) N,A(NBAND,N)
170   FORMAT(1H1,5X,31H SET OF EQUATIONS ARE SINGULAR ,///,5X,25HDIAGONA
     &L TERM OF EQUATION,I5,13H IS EQUAL TO ,E15.8)
      STOP
      END
C
C
      SUBROUTINE SCALE3(AR, RR, NN, II)
      IMPLICIT REAL*4(A-H,O-Z)
      DIMENSION AR(NN)
      AMIN=1.E30
      AMAX=-1.E30
      DO 100 I=1,NN
        AMAX=MAX(AR(I),AMAX)
        AMIN=MIN(AR(I),AMIN)
100   CONTINUE
      IF(AMIN.EQ.AMAX) AMAX=AMIN+1.
      AR(NN+1)=AMIN
      AR(NN+2)=(AMAX-AMIN)/RR
      AMAX=AMAX+50.*AR(NN+2)
      AMIN=AMIN-50.*AR(NN+2)
      AR(NN+1)=AMIN
      AR(NN+2)=(AMAX-AMIN)/RR
      RETURN
      END
C
      SUBROUTINE SCALE2(AR, RR, NN, II)
      IMPLICIT REAL*4(A-H,O-Z)
      DIMENSION AR(NN)
      AMIN=1.E30
      AMAX=-1.E30
      DO 100 I=1,NN
        AMAX=MAX(AR(I),AMAX)
        AMIN=MIN(AR(I),AMIN)
100   CONTINUE
      IF(AMIN.EQ.AMAX) AMAX=AMIN+1.
      AR(NN+1)=AMIN
      AR(NN+2)=(AMAX-AMIN)/RR
      RETURN
      END
C
      SUBROUTINE SORTIX(N,IX)
      DIMENSION IX(N)
      LOGICAL ANYEXC
      ANYEXC=.TRUE.
C WHILE ANY EXCHANGES
30    IF(ANYEXC) THEN
        ANYEXC=.FALSE.
        DO 40 I=1,N-1
          IF(IX(I).LT.IX(I+1)) THEN
            CALL SWAP(IX(I),IX(I+1))
            ANYEXC=.TRUE.
          ENDIF
40      CONTINUE
        GOTO 30
      ENDIF
      RETURN
      END
C
      SUBROUTINE SWAP(IX1,IX2)
      ITEMP=IX1
      IX1=IX2
      IX2=ITEMP
      RETURN
      END
C
      SUBROUTINE ELMDUP(NN,IX)
      DIMENSION IX(NN)
      NNOUT=NN
      DO 100 N=2,NN
80      IF(IX(N-1).EQ.IX(N) .AND. (IX(N).NE.0)) THEN
          NNOUT=NNOUT-1
          DO 90 J=N,NN-1
            IX(J)=IX(J+1)
90        CONTINUE
          IX(NN)=0
          GOTO 80
        ENDIF
100   CONTINUE
      NN=NNOUT
      RETURN
      END
      INTEGER FUNCTION IFIND(X, Y, XC, YC)
      IMPLICIT REAL*4(A-H,O-Z)
      DIMENSION XC(5),YC(5),TEST(4)
      XC(5)=XC(1)
      YC(5)=YC(1)
      DO 100 I=1,4
        AX=X-XC(I)
        AY=Y-YC(I)
        BX=XC(I+1)-XC(I)
        BY=YC(I+1)-YC(I)
        TEST(I)=-(AX*BY-AY*BX)
100   CONTINUE
      IF((TEST(1).GE.0.) .AND.
     &   (TEST(2).GE.0.) .AND.
     &   (TEST(3).GE.0.) .AND.
     &   (TEST(4).GE.0.))THEN
        IFIND=1
      ELSE
        IFIND=0
      ENDIF
      RETURN
      END
      FUNCTION AFUNCT(TIME, IDOT, AMASS, HT, BD, SL, X, Y)
      IMPLICIT REAL*4(A-H,O-Z)
      COMMON /SNO/ SNOLIN,SNO(20)
      DIMENSION AMASS(9)
      FPC(S,E)=(-7.044897460E+02)*S**2+( 3.286502080E+01)*S+
     1  (-8.363957700E-08)*E**2+( 2.803513780E-04)*E+
     1  (-8.280271290E-02)
      FPX(S,E)=(-1.011654790E+03)*S**2+( 5.757759090E+01)*S+
     1  (-1.945462600E-07)*E**2+( 6.485106420E-04)*E+
     1  (-1.714566350E-01)
      FSX(S,E)=( 1.613718720E+02)*S**2+(-3.064633370E+00)*S+
     1  (-2.274655910E-07)*E**2+( 9.221490470E-04)*E+
     1  (-3.658012150E-01)
      FTM(S,E)=(-2.475505130E+03)*S**2+( 7.388871770E+01)*S+
     1  (-7.129281700E-07)*E**2+( 3.180227710E-03)*E+
     1  (-6.778311130E-01)
      FSM(S,E)=(-1.272729690E+04)*S**2+( 1.028782810E+02)*S+
     1  ( 1.790857030E-07)*E**2+(-1.209567770E-03)*E+
     1  ( 2.244830130E+00)
C
C FOLLOWING SECTION IS QUITE VARIABLE. IT CONTAINS CALLS TO COEF AND
C COEF1 TO VARY VARIOUS MASS BALANCE PARAMETERS IN A SINUSOIDAL FASHION,
C AND WILL BE DIFFERENT FOR DIFFERENT EXPERIMENTS.
C
C     SNOLIN=0.
C       AMASS(7)=COEF(TIME,0.,10500.,-900.,-2100.)
C     IF(TIME.LT.70000.) THEN
C       AMASS(7)=COEF1(TIME,0.,20000.,0.,-2300.)
C     ELSE
C       AMASS(7)=COEF1(TIME,70000.,90000.,-2300.,0.)
C     ENDIF
C     DELSNNP=COEF(TIME,0.,20500.,0.,-2300.)
C     AMASS(7)=DELSNNP
C     DELSNEQ=COEF(TIME,-25000.,-14500.,5600.,3800.)
C     AMASS(8)=(DELSNEQ-DELSNNP)*1E-7
C     IF(TIME.LE.20500.) THEN
C       AMASS(7)=COEF1(TIME,0.,20500.,0.,-2300.)
C     ELSEIF(TIME.LE.61500.) THEN
C       AMASS(7)=COEF(TIME,20500.,41000.,-2300.,-1000.)
C     ELSE
C       AMASS(7)=COEF1(TIME,61500.,82000.,-2300.,0.)
C     ENDIF
C     IF(TIME.LE.20500.) THEN
C       AMASS(7)=COEF1(TIME,0.,20500.,-2000.,0.)
C     ELSEIF(TIME.LE.61500.) THEN
C       AMASS(7)=COEF(TIME,20500.,41000.,0.,-3000.)
C     ELSE
C       AMASS(7)=COEF1(TIME,61500.,82000.,0.,-2000.)
C     ENDIF
C       AMASS(7)=COEF(TIME,-3600.,6400.,0.,-2000.)
C       AMASS(7)=COEF(TIME,-3000.,7000.,-1000.,-4000.)
C END VARIABLE SECTION
C
C BRANCH FOR VARIOUS MASS BALANCE RELATIONSHIPS
      GOTO(10,20,30,40,50,60,70,200,210,220,230,240,
     1     300,310,320,330,340,345,347,348),IDOT
C
10    CONTINUE
      AMARG=AMASS(1)
      APEAK=AMASS(2)
      ADOME=AMASS(3)
      HE=AMASS(4)
      HP=AMASS(5)
      HM=AMASS(6)
      GOTO 100
C
20    CONTINUE
C **** IDOT=2 POLAR - A
      AMARG=0.0
      APEAK=.5
      ADOME=.05
      HE=0.
      HP=1250.
      HM=2500.
      SNOLIN=AMASS(7)+AMASS(8)*SQRT(X**2+Y**2)
      GOTO 100
C
30    CONTINUE
C **** IDOT=3 POLAR - B
      AMARG=0.0
      APEAK=1.
      ADOME=.10
      HE=0.
      HP=1250.
      HM=2500.
      SNOLIN=AMASS(7)+AMASS(8)*SQRT(X**2+Y**2)
      GOTO 100
C
40    CONTINUE
C **** IDOT=4 POLAR - C
      AMARG=-1.0
      APEAK=1.5
      ADOME=.1
      HE=100.
      HP=1250.
      HM=2500.
      SNOLIN=AMASS(7)+AMASS(8)*SQRT(X**2+Y**2)
      GOTO 100
C
50    CONTINUE
C **** IDOT=5 MARITIME - A
      AMARG=-2.0
      APEAK=2.
      ADOME=.15
      HE=1500.
      HP=1750.
      HM=2000.
      SNOLIN=AMASS(7)+AMASS(8)*SQRT(X**2+Y**2)
      GOTO 100
C
60    CONTINUE
C **** IDOT=6 MARITIME - B
      AMARG=-2.0
      APEAK=2.
      ADOME=.15
      HE=1000.
      HP=1500.
      HM=2000.
      SNOLIN=AMASS(7)+AMASS(8)*SQRT(X**2+Y**2)
      GOTO 100
C
70    CONTINUE
C **** IDOT=7 MARITIME - C
      AMARG=-3.0
      APEAK=2.
      ADOME=.15
      HE=1000.
      HP=1500.
      HM=2000.
      SNOLIN=AMASS(7)+AMASS(8)*SQRT(X**2+Y**2)
      GOTO 100
C
100   CONTINUE
      AMASS(1)=AMARG
      AMASS(2)=APEAK
      AMASS(3)=ADOME
      AMASS(4)=HE
      AMASS(5)=HP
      AMASS(6)=HM
      HTT=HT-SNOLIN
      IF(HTT.GT.HM) THEN
        AFUNCT=ADOME
      ELSE
        IF(HTT.LE.0.) THEN
          AFUNCT=AMARG
        ELSE
          IF(HTT.GT.0. .AND. HTT.LE.HE) THEN
            A1=AMARG
            H1=0.
            A2=0.
            H2=HE
          ENDIF
          IF(HTT.GT.HE .AND. HTT.LE.HP) THEN
            A1=0.
            H1=HE
            A2=APEAK
            H2=HP
          ENDIF
          IF(HTT.GT.HP .AND. HTT.LE.HM) THEN
            A1=APEAK
            H1=HP
            A2=ADOME
            H2=HM
          ENDIF
          AFUNCT=(A1*(H2-HTT)-A2*(H1-HTT))/(H2-H1)
        ENDIF
      ENDIF
      IF(BD.LE.-9999.) AFUNCT=-5.
      RETURN
C
200   CONTINUE
C **** IDOT=8 SM (SUBPOLAR MARITIME) SLOPE DEPENDENT
      ADOME=.10
      SMAX=.01
      B0=-5.
      HEQ=700.
      B0DOT=B0/HEQ
      BNEQ=FSM(SMAX,HEQ)
      HADD=-BNEQ/B0DOT
      IF(HT.LE.HEQ) THEN
        BN=B0-B0DOT*HT
      ELSE IF(HT.LE.(HEQ+HADD)) THEN
        BN=FSM(SL,HT)+B0DOT*(HADD-HT+HEQ)
        IF(BN.LE.0.) BN=ADOME
      ELSE
        BN=FSM(SL,HT)
        IF(BN.LE.0.) BN=ADOME
      ENDIF
      AFUNCT=BN
      GOTO 1000
C
210   CONTINUE
C **** IDOT=9 TM (TEMPERATE MARITIME) SLOPE DEPENDENT
      ADOME=.10
      SMAX=.03
      B0=-8.
      HEQ=1125.
      B0DOT=B0/HEQ
      BNEQ=FTM(SMAX,HEQ)
      HADD=-BNEQ/B0DOT
      IF(HT.LE.HEQ) THEN
        BN=B0-B0DOT*HT
      ELSE IF(HT.LE.(HEQ+HADD)) THEN
        BN=FTM(SL,HT)+B0DOT*(HADD-HT+HEQ)
        IF(BN.LE.0.) BN=ADOME
      ELSE
        BN=FTM(SL,HT)
        IF(BN.LE.0.) BN=ADOME
      ENDIF
      AFUNCT=BN
      GOTO 1000
C
220   CONTINUE
C **** IDOT=10 SX (SUBPOLAR MIX) SLOPE DEPENDENT
      ADOME=.10
      SMAX=.07
      B0=-5.
      HEQ=1250.
      B0DOT=B0/HEQ
      BNEQ=FSX(SMAX,HEQ)
      HADD=-BNEQ/B0DOT
      IF(HT.LE.HEQ) THEN
        BN=B0-B0DOT*HT
      ELSE IF(HT.LE.(HEQ+HADD)) THEN
        BN=FSX(SL,HT)+B0DOT*(HADD-HT+HEQ)
        IF(BN.LE.0.) BN=ADOME
      ELSE
        BN=FSX(SL,HT)
        IF(BN.LE.0.) BN=ADOME
      ENDIF
      AFUNCT=BN
      GOTO 1000
C
230   CONTINUE
C **** IDOT=11 PX (POLAR MIX) SLOPE DEPENDENT
      ADOME=.05
      SMAX=.04
      B0=-2.5
      HEQ=600.
      B0DOT=B0/HEQ
      BNEQ=FPX(SMAX,HEQ)
      HADD=-BNEQ/B0DOT
      IF(HT.LE.HEQ) THEN
        BN=B0-B0DOT*HT
      ELSE IF(HT.LE.(HEQ+HADD)) THEN
        BN=FPX(SL,HT)+B0DOT*(HADD-HT+HEQ)
        IF(BN.LE.0.) BN=ADOME
      ELSE
        BN=FPX(SL,HT)
        IF(BN.LE.0.) BN=ADOME
      ENDIF
      AFUNCT=BN
      GOTO 1000
C
240   CONTINUE
C **** IDOT=12 PC (POLAR CONTINENTAL) SLOPE DEPENDENT
      ADOME=.05
      SMAX=.01
      B0=-1.
      HEQ=300.
      B0DOT=B0/HEQ
      BNEQ=FPC(SMAX,HEQ)
      HADD=-BNEQ/B0DOT
      IF(HT.LE.HEQ) THEN
        BN=B0-B0DOT*HT
      ELSE IF(HT.LE.(HEQ+HADD)) THEN
        BN=FPC(SL,HT)+B0DOT*(HADD-HT+HEQ)
        IF(BN.LE.0.) BN=ADOME
      ELSE
        BN=FPC(SL,HT)
        IF(BN.LE.0.) BN=ADOME
      ENDIF
      AFUNCT=BN
      GOTO 1000
C
300   CONTINUE
C **** IDOT=13 SM EXPONENTIAL
      A1=-7.57655
      A2=2.56891
      RL1=2.773327E-6
      RL2=2.916233E-7
      SNOLIN=SNO(13)
      SNOLIN=AMASS(7)+AMASS(8)*SQRT(X**2+Y**2)
      GOTO 350
C
310   CONTINUE
C **** IDOT=14 TM EXPONENTIAL
      A1=-11.6877
      A2=3.67407
      RL1=1.083791E-6
      RL2=3.641119E-8
      SNOLIN=SNO(14)
      SNOLIN=AMASS(7)+AMASS(8)*SQRT(X**2+Y**2)
      GOTO 350
C
320   CONTINUE
C **** IDOT=15 SX EXPONENTIAL
      A1=-5.685629
      A2=0.857527
      RL1=1.322827E-6
      RL2=9.461793E-8
      SNOLIN=SNO(15)
      SNOLIN=AMASS(7)+AMASS(8)*SQRT(X**2+Y**2)
      GOTO 350
C
330   CONTINUE
C **** IDOT=16 PX EXPONENTIAL
      A1=-2.84904
      A2=.837780
      RL1=8.534590E-6
      RL2=3.405352E-8
      SNOLIN=SNO(16)
      SNOLIN=AMASS(7)+AMASS(8)*SQRT(X**2+Y**2)
      GOTO 350
C
340   CONTINUE
C **** IDOT=17 PC EXPONENTIAL
      A1=-1.29363
      A2=.299798
      RL1=1.502293E-5
      RL2=3.661944E-8
      SNOLIN=SNO(17)
      SNOLIN=AMASS(7)+AMASS(8)*SQRT(X**2+Y**2)
      GOTO 350
C
345   CONTINUE
C **** IDOT=18 ANTARCTIC VALUUES
C **** FOLLOWING IS CATCHALL FOR DIFFERENT MASS BALANCE RELATIONSHIPS
C **** AND MAY HAVE DIFFERENT VALUES FOR DIFFERENT EXPERIMENTS
C **** DOUBLE OF 17
      A1=-2.*1.29363
      A2=2.*.299798
      RL1=1.502293E-5
      RL2=3.661944E-8
C **** CTRL VALUES
C     A1=-.838547
C     A2=.887967
C     RL1=3.471435E-7
C     RL2=3.918019E-8
C **** REGSST VALUES
C     A1=-1.22491
C     A2=1.21169
C     RL1=3.736239E-7
C     RL2=8.897420E-8
C **** WARM MIN VALUES
C     A1=-15.8931
C     A2=15.3500
C     RL1=3.919554E-7
C     RL2=2.362680E-7
      SNOLIN=SNO(18)
      SNOLIN=AMASS(7)+AMASS(8)*SQRT(X**2+Y**2)
      GOTO 350
C
347   CONTINUE
C **** IDOT=19 GROSSWALD'S DIPPING SNOWLINE SCHEME
      DIST=SQRT(X**2+Y**2)
C FOLLOWING FOR 100M/100KM WITH SNOWLINE AT -1000 M AT POLE
      SNSLOPE=.001
C FOLLOWING IS FOR SIBERIAN SNOWLINE DATA, 2E-03
      SNSLOPE=2.E-03
      SNSLOPE=AMASS(8)
      SNZERO=AMASS(7)
C FOLLOWING FOR -1 M AT S.L., +1 M AT 2000 M ELEVATION, OR 2M/2000M
      BALGRAD=.001
      BALGRAD=AMASS(9)
      SNLINE=SNZERO+SNSLOPE*DIST
      ACC=(HT-SNLINE)*BALGRAD
      IF(ACC.GT.2.) THEN
        ACC=ACC-2.*(ACC-2.)
        IF(ACC.LT..05) ACC=.05
      ELSEIF(ACC.LT.-2.) THEN
        ACC=-2.
      ENDIF
      AFUNCT=ACC
      RETURN
348   CONTINUE
C **** IDOT=20 MIKE'S METEROLOGICAL MODEL
      TNSL=AMASS(9)
c variable stuff, to modify TNSL
c     TNSL=COEF(TIME,-9550.,10950.,0.,-20.)
c     TNSL=COEF(TIME,0.,2500.,0.,-20.)
c end variable stuff of TNSL
      AMASS(9)=TNSL
      AFUNCT=ACCUM(X*.001,Y*.001,HT,SL,0.,TNSL)
      RETURN
350   CONTINUE
C
C AGAIN A VARIABLE SECTION FOR DIFFERENT EXPERIMENTS, HERE VARYING
C SPECIFIC PARAMETERS IN THE DIFFERENT MASS BALANCE SCHEMES. THIS CAN BE
C USED TO SMOOTHLY SWITCH FROM ONE ZONE TO ANOTHER.
C FOLLOWING FOR TIME VARYING ACCUM CURVES
C COEF IS FUNCTION COEF(TIME,TINIT,TFINAL,VINIT,VFINAL)
C COEF DOES SINUSOIDAL VARIATION REGARDLES OF TIME
C COEF1 FOR TIME<TINIT = VINIT, TIME>TFINAL = VFINAL, IN BETWEEN IS
C COSINE INTERPOLATION
C FOLLOWING ARE INITIAL TIME AND FINAL TIME (HALF CYCLE)
C     TINIT=0.
C     TFINAL=5000.
C FOLLOWING 4 FOR CONSTANT
      A11=A1
      A22=A2
      RL11=RL1
      RL22=RL2
C FOLLOWING 4 FOR TIME VARYING SINUSOIDAL
C     A11=COEF(TIME,TINIT,TFINAL,A1,A1)
C     A22=COEF(TIME,TINIT,TFINAL,A2,2.*A2)
C     RL11=COEF(TIME,TINIT,TFINAL,RL1,RL1*.01)
C     RL22=COEF(TIME,TINIT,TFINAL,RL2,RL2)
C FOLLOWING 4 FOR COSINE INTERPOLATED STEP
C BEGINNING AND END OF TIME STEP
C     TINIT=1000.
C     TFINAL=22000.
C     A11=COEF1(TIME,TINIT,TFINAL,A1,-14.074)
C     A22=COEF1(TIME,TINIT,TFINAL,A2,4.83112)
C     RL11=COEF1(TIME,TINIT,TFINAL,RL1,1.308677E-06)
C     RL22=COEF1(TIME,TINIT,TFINAL,RL2,1.61273E-7)
C     A11=COEF1(TIME,TINIT,TFINAL,A1,A1*18.95)
C     A22=COEF1(TIME,TINIT,TFINAL,A2,17.287*A2)
C     RL11=COEF1(TIME,TINIT,TFINAL,RL1,RL1*1.129)
C     RL22=COEF1(TIME,TINIT,TFINAL,RL2,6.03*RL2)
C     A11=COEF1(TIME,TINIT,TFINAL,A1,A1*1.461)
C     A22=COEF1(TIME,TINIT,TFINAL,A2,1.365*A2)
C     RL11=COEF1(TIME,TINIT,TFINAL,RL1,RL1*1.076)
C     RL22=COEF1(TIME,TINIT,TFINAL,RL2,2.271*RL2)
C END OF VARIABLE SECTION
C
      HTT=HT-SNOLIN
      IF(HTT.LT.0.) HTT=0.
      ARG=-RL11*HTT**2
      IF(ARG.GT.-180.) THEN
        TERM1=A11*EXP(ARG)
      ELSE
        TERM1=0.
      ENDIF
      ARG=-RL22*HTT**2
      IF(ARG.GT.-180.) THEN
        TERM2=A22*EXP(ARG)
      ELSE
        TERM2=0.
      ENDIF
      AFUNCT=TERM1+TERM2
1000  CONTINUE
      RETURN
      END
      FUNCTION COEF(TIME, TINIT, TFINAL, VINIT, VFINAL)
      IMPLICIT REAL*4(A-H,O-Z)
C STARTING AT TIME=TINIT SINUSOIDALLY VARY FROM VINIT TO VFINAL WITH
C PERIOD = TFINAL-TINIT. CONTINUE PERIODICITY.
      IF(TIME.LT.TINIT) THEN
        COEF=VINIT
      ELSE
        ZEROL=0.5*(VINIT+VFINAL)
        AMPL=0.5*(VINIT-VFINAL)
        PHASE=(TIME-TINIT)/(TFINAL-TINIT)*3.1415927
        COEF=ZEROL+AMPL*COS(PHASE)
      ENDIF
      END
C
      FUNCTION COEF1(TIME, TINIT, TFINAL, VINIT, VFINAL)
      IMPLICIT REAL*4(A-H,O-Z)
C STARTING AT TIME TINIT VARY FROM VINIT TO VFINAL SUNUSOIDALLY. FOR
C TIME < TINIT VALUE=VINIT, FOR TIME > VFINAL VALUE= VFINAL.
      IF(TIME.LE.TINIT) THEN
        COEF1=VINIT
      ELSE IF(TIME.GE.TFINAL) THEN
        COEF1=VFINAL
      ELSE
        ZEROL=0.5*(VINIT+VFINAL)
        AMPL=0.5*(VINIT-VFINAL)
        PHASE=(TIME-TINIT)/(TFINAL-TINIT)*3.1415927
        COEF1=ZEROL+AMPL*COS(PHASE)
      ENDIF
      END
      SUBROUTINE READN(MAX,IDEP,HED,NUMNP,NUMEL,NUMGBC,NDT,INTER,DT,
     &KODE,X,Y,HTICE,ADOT,FRACT,PSURF,RHOI,BDROCK,UNDEPB,FLOWA,SLDGB,
     &THICK,KX,CONST,IBFLUX,BFLUX,QHOLD,NTYPE,NNODE,NCOL,DEPB,
     &AADOT,AFRACT,ABDRCK,PPSURF,AFLOWA,ASLDGB,IDT,AMASS)
      IMPLICIT REAL*4(A-H,O-Z)
C READ FROM SPLIT DATA SETS, PRODUCED BY UNPACK, REASSEMBLED BY PACK
      COMMON /SNO/ SNOLIN,SNO(20)
      CHARACTER HED*80
      DIMENSION AMASS(9)
      DIMENSION KODE(MAX),X(MAX),Y(MAX),HTICE(MAX),ADOT(MAX),FRACT(MAX),
     &PSURF(MAX),BDROCK(MAX),FLOWA(MAX),SLDGB(MAX),THICK(MAX),
     &KX(MAX,4),CONST(MAX),IBFLUX(MAX,2),BFLUX(MAX),QHOLD(MAX),IDT(MAX),
     &NNODE(MAX),NTYPE(MAX),AADOT(MAX),AFRACT(MAX),ABDRCK(MAX),
     &PPSURF(MAX),AFLOWA(MAX),ASLDGB(MAX),UNDEPB(MAX),LM(4),
     &CFIT(1000),DEPB(MAX)
      DO 17 I=1,20
        SNO(I)=0.
17    CONTINUE
      IDEP=1
C
C **** SPECIAL REMOVE FOR DEPRESSED BED
C     IDEP=0
C
      RHOI=.917
      RHOR=4.0
C    MAKE SURE RHOR IS ALSO 4.0 IN OTHER ROUTINES, IT IS NOT PASSED
C
C READ INPUT HEADER
      READ(30,1000,END=999) HED,NUMNP,NUMEL,NUMCOL,NUMLEV,NUMGBC,NDT,
     &                        INTER,DT
C
C READ INPUT GRID, THINGS THAT NEVER CHANGE
      READ(31) HED
      READ(31) (KODE(I),I=1,NUMNP)
      READ(31) (X(I),I=1,NUMNP)
      READ(31) (Y(I),I=1,NUMNP)
      READ(31) (PSURF(I),I=1,NUMNP)
      READ(31) (BDROCK(I),I=1,NUMNP)
      READ(31) (KX(I,1),KX(I,2),KX(I,3),KX(I,4),I=1,NUMEL)
      READ(31) (IBFLUX(I,1),IBFLUX(I,2),BFLUX(I),I=1,NUMGBC)
C
C READ INPUT DIFFER, THINGS THAT DIFFER FROM ONE GRID TO THE NEXT
      READ(32) HED
      READ(32) (ADOT(I),I=1,NUMNP)
      READ(32) (FRACT(I),I=1,NUMNP)
      READ(32) (FLOWA(I),I=1,NUMNP)
      READ(32) (SLDGB(I),I=1,NUMNP)
      DO 123 I=1,NUMNP
        IF(FLOWA(I).LE.0.) FLOWA(I)=1E-6
        IF(SLDGB(I).LE.0.) SLDGB(I)=1E-9
123   CONTINUE
C
C READ INPUT TIME, THINGS THAT CHANGE WITH TIME
      READ(33) HED
      READ(33) (HTICE(I),I=1,NUMNP)
      READ(33) (AJUNK,I=1,NUMNP)
      READ(33) (AJUNK,I=1,NUMNP)
      READ(33) (CONST(I),I=1,NUMEL)
1000  FORMAT (A80,/,7I5,F8.0)
c tag on vriable section
      print *,'no variation'
c     print *,'vary tnsl from -14 to -18 then up to -10'
      IF(NUMNP.GT.MAX) THEN
        PRINT *,'NUMNP=',NUMNP,' MAX=',MAX,' INCREASE MAX'
        STOP
      ENDIF
      DO 100 N=1,NUMNP
        IDT(N)=0
        IF(ADOT(N).LT.-100.) THEN
         IDT(N)=-(100+ADOT(N))
         ADOT(N)=AFUNCT(0.,IDT(N),AMASS,HTICE(N),BDROCK(N),0.,X(N),Y(N))
        ENDIF
        IF(PSURF(N).LT.BDROCK(N)) PSURF(N)=BDROCK(N)
        THICK(N)=HTICE(N)-BDROCK(N)
        IF(IDEP.EQ.1 .AND. BDROCK(N).GT.-9999.) THEN
          UNDEPB(N)=RHOI*(PSURF(N)-BDROCK(N))/RHOR+BDROCK(N)
        ELSE
          UNDEPB(N)=BDROCK(N)
        ENDIF
100   CONTINUE
      NELMAX=-10000
      NCLMAX=-10000
      DO 90 N=1,NUMEL
        NTYPE(N)=1
        IF(KX(N,4).EQ.0) NTYPE(N)=2
        NODEN=4
        NNODE(N)=4
        IF(NTYPE(N).EQ.2) NODEN=3
        IF(NTYPE(N).EQ.2) NNODE(N)=3
        IMAX=-10000
        IMIN=10000
        DO 87 I=1,NODEN
          IF(KX(N,I).GT.IMAX) IMAX=KX(N,I)
          IF(KX(N,I).LT.IMIN) IMIN=KX(N,I)
87      CONTINUE
C **** FOR MAP1 WITH ASYMSL ****
        NCOL1=2*(IMAX-IMIN)+1
C *** FOR MAP2 WITH ESSL ROUTINE SGBF/SGBS
        NCOL11=3*(IMAX-IMIN)+16
        IF(NCOL1.GT.NCLMAX) THEN
          NCLMAX=NCOL1
          NELMAX=N
        ENDIF
90    CONTINUE
      WRITE(*,*) 'MAX NCOL=',NCLMAX,' IN ELEMENT',NELMAX
      IF(NCLMAX.GT.NCOL) THEN
        WRITE(*,*) 'CURRENTLY USING NCOL=',NCOL,' INCREASE'
        STOP
      ENDIF
C
      SUMC = 0.0
      DO 95 N = 1,NUMEL
        QHOLD(N) = 0.0
        IF (CONST(N).EQ.0.0 .AND. DT.EQ.0.0) CONST(N)=1.92E7
        SUMC = SUMC + CONST(N)
        IF(NTYPE(N).EQ.1) THEN
          NODEN=4
          NINT=9
        ELSE
          NODEN=3
          NINT=4
        ENDIF
        AAADOT=0.
        AAFRCT=0.
        AABDRK=0.
        APSURF=0.
        AAFLOW=0.
        AASLDG=0.
        DO 125 I=1,4
          LM(I)=KX(N,I)
          AAADOT=AAADOT+ADOT(LM(I))
          AAFRCT=AAFRCT+FRACT(LM(I))
          AABDRK=AABDRK+BDROCK(LM(I))
          APSURF=APSURF+PSURF(LM(I))
          AAFLOW=AAFLOW+FLOWA(LM(I))
          AASLDG=AASLDG+SLDGB(LM(I))
125     CONTINUE
        DENOM=1./FLOAT(NODEN)
        AAADOT=AAADOT*DENOM
        AADOT(N)=AAADOT
        AAFRCT=AAFRCT*DENOM
        AFRACT(N)=AAFRCT
        AABDRK=AABDRK*DENOM
        ABDRCK(N)=AABDRK
        APSURF=APSURF*DENOM
        PPSURF(N)=APSURF
        AAFLOW=AAFLOW*DENOM
        AFLOWA(N)=AAFLOW
        AASLDG=AASLDG*DENOM
        ASLDGB(N)=AASLDG
95    CONTINUE
C
      SUMC = SUMC/NUMEL
      WRITE(*,*) SUMC
1001  FORMAT(2I4,2E12.5,F10.2,F7.2,F9.2,F10.3,F10.1,F10.2,F10.3)
1002  FORMAT(5I5,1PE17.10)
1007  FORMAT(2I5,E13.6)
      PG=0.089866
      IDEP0=1
      RETURN
999   CONTINUE
      PRINT *,'NO DATA SET'
      STOP
      END
      FUNCTION ACCUM(X,Y,ELEV1,SLOPE1,SHAPE1,TNSL)
      DIMENSION QI(12),QS(12),TTT(12)
      DATA QI /960.,1036.,1200.,825.,330.,90.,150.,600.,1200.,1020.,
     &         930.,850./
      DATA QS /-0.667,4.6,11.667,9.167,3.667,1.,1.667,6.667,12.,6.333,
     &         0.333,-3.333/
C     DATA AAA /-9.14/, BBB /-.68/, CCC /34.461/
C     DATA WWW /13.05/, XXX /.664/, ZZZ /2.608/
       DATA AAA / -9.62376690     /
       DATA BBB /-0.546917617     /
       DATA CCC /  24.9793854     /
       DATA WWW /  19.1390686     /
       DATA XXX / 0.922791243     /
       DATA ZZZ /-0.738900483     /
C     WRITE(19,*) X,Y,ELEV1,SLOPE1
C CALCULATE LATITUDE
      CALL SETRIG
      CALL RECPOL(X,Y,RLAT,RLONG)
C     PRINT *,'LATITUDE=',RLAT
C ELEVATION (KM)
      ELEV=ELEV1/1000.
C SLOPE (M/KM)
      SLOPE=SLOPE1*1000.
C SHAPE (M/KM/KM)
      SHAPE=SHAPE1*1000.*1000.
      SHAPE=0.
C CALCULATE SURFACE MEAN ANNUAL AIR TEMP
      TS=AAA*ELEV+BBB*RLAT+CCC+TNSL+14.0
C     PRINT *,'SURFACE MEAN ANNUAL AIR TEMP=',TS
C CALCULATE MEAN ANNUAL TEMP OF FREE ATMOSPHERE-ISOTHEMAL LAYER
      TF=0.67*(TS+273.0)+88.9
C     PRINT *,'TEMP FREE AT-ISOTHRMAL LAYER=',TF
C CALCULATE SATURATION VAPOR PRESSURE
      TERM1=-9.09718*(273.16/TF-1.0)
      TERM2=-3.56654*LOG10(273.16/TF)
      TERM3=0.876793*(1.0-TF/273.16)+0.785835
      EXPON=TERM1+TERM2+TERM3
      ES=10.**EXPON
C     PRINT *,'SATURATION VAPOR PRESSURE=',ES
C CALCULATE ACCUMULATION RATE (M/YR)
      TERM1=WWW*ES
      TERM2=XXX*SLOPE
      TERM3=ZZZ
      TERM4=-15.276*SHAPE
      ACC=TERM1+TERM2+TERM3+TERM4
C     WRITE(19,113) TERM1,TERM2,TERM3,ACC
113   FORMAT(4G13.6)
C     PRINT *,'ACCUMULATION=',ACC
C CALCULATE ABLATION
      QY=0.
      DO 10 I=1,12
        QY=QY+QI(I)-QS(I)*RLAT
10    CONTINUE
      QY=QY/12.
      PDD=0.
      DO 20 I=1,12
C       TTT(I)=TS+0.021*((QI(I)+QS(I)*RLAT)-QY)+8.954
        TTT(I)=TS+0.021*((QI(I)-QS(I)*RLAT)-QY)+0.0
C       WRITE(19,*) I,TTT(I),QI(I)-QS(I)*RLAT,QY
        IF(TTT(I).GT.0.0) PDD=PDD+30.*TTT(I)
20    CONTINUE
      ABL=.6*PDD
C     PRINT *,'ABLATION=',ABL
C CALCULATE NET ACCUMULATION
      ACCNET=ACC-ABL
C     ACCNET=ACC
C     PRINT *,'NET ACCUMULATION/ABLATION=',ACCNET
      ACCUM=ACCNET*.01
C     IF(PDD.GT.0.) THEN
C       WRITE(19,111) ELEV,SLOPE,RLAT,TS,TF-273.,ES,PDD,ACC,ABL,ACCNET
C     ENDIF
111   FORMAT(6F7.2,F6.0,3F7.1)
C     WRITE(19,112) (TTT(I),I=1,12)
112   FORMAT(12F6.1)
      END
      SUBROUTINE SETRIG
      IMPLICIT REAL*4(A-H,O-Z)
      COMMON/TRIG/RKMPDEG,RADPDEG,DEGPRAD
      RADIUS=6370.D0*0.53
      PI=4.D0*ATAN(1.D0)
      RADIUS=2.D4/PI
      CIRCUM=2.D0*PI*RADIUS
      RKMPDEG=CIRCUM/360.D0
      RADPDEG=PI/180.D0
      DEGPRAD=180.D0/PI
      END
      SUBROUTINE POLREC(RLAT,RLONG,X,Y)
      IMPLICIT REAL*4(A-H,O-Z)
      COMMON/TRIG/RKMPDEG,RADPDEG,DEGPRAD
      X=1000.D0*(90.D0-RLAT)*RKMPDEG*COS(RLONG*RADPDEG)
      Y=1000.D0*(90.D0-RLAT)*RKMPDEG*SIN(RLONG*RADPDEG)
      END
      SUBROUTINE RECPOL(X,Y,RLAT,RLONG)
      IMPLICIT REAL*4(A-H,O-Z)
      COMMON/TRIG/RKMPDEG,RADPDEG,DEGPRAD
      IF(X.EQ.0.D0) THEN
        IF(Y.GE.0.D0) THEN
          THETA=90.D0
        ELSE
          THETA=270.D0
        ENDIF
      ELSE
        THETA=DEGPRAD*ATAN(Y/X)
        IF(X.LT.0.D0) THEN
          THETA=180.D0+THETA
        ELSEIF(Y.LT.0.D0) THEN
          THETA=360.D0+THETA
        ELSEIF(THETA.GT.360.D0) THEN
          THETA=THETA-360.D0
        ENDIF
      ENDIF
C     R=.001D0*SQRT(X**2+Y**2)
      R=1.D0*SQRT(X**2+Y**2)
      RLAT=90.D0-R/RKMPDEG
      RLONG=THETA
      END
      SUBROUTINE ADJUST(HED, NUMNP, NUMEL, XX, YY, HTICE, ADOT, FRACT,
     &                PSURF, BED, DEPB, FLOWA, SLDGB, THICK, KX, CONST,
     &                  NNODE, KODE, HFIT, NUMCOL, NUMLEV, NUMGBC,
     &                  NDT, INTER, DT, IBFLUX, BFLUX, MXX, IDT,
     &                  SLOPN, AMASS, TIME, NTSTEP, TTIME, VVOL)
      PARAMETER( MMXX=29999 )
      IMPLICIT REAL*4(A-H,O-Z)
      EXTERNAL IFIND
      COMMON /SNO/ SNOLIN,SNO(20)
      DIMENSION TTIME(NTSTEP),VVOL(NTSTEP)
      DIMENSION AMASS(9), KODE(MXX), XX(MXX), YY(MXX), HTICE(MXX),
     &          SLOPN(MXX), ADOT(MXX), FRACT(MXX), PSURF(MXX),
     &          BED(MXX), DEPB(MXX), FLOWA(MXX), SLDGB(MXX),
     &          KX(MXX,4), CONST(MXX), THICK(MXX), NNODE(MXX),
     &          ZZ(MMXX), ICMAP(16), XCHECK(5), YCHECK(5), IFD(MMXX),
     &          IFP(MMXX), HFIT(MXX), KXX(4),
     &          IBFLUX(MXX,2), BFLUX(MXX), NLINE(200), IDT(MXX),
     &          XA(5), YA(5), IDAT(5)
      REAL*4 XA,YA,IGOT
      CHARACTER COORDT*6, COORD1*5, CHAR*1, HED*80, JUNK*80
      DATA ICMAP /0,5,11,4,12,6,13,2,8,7,9,3,10,14,15,1/
      FLOWMIN=1E-6
      SLDGBMIN=1E-9
      IPASS=0
      ICOUNT=NUMNP
      DO 5 I=1,NUMNP
        IFP(I)=I
        ZZ(I)=HTICE(I)-HFIT(I)
5     CONTINUE
      XLEN=100.
      YLEN=100.
      ZLEN=14.
      ILINE=0
      GOTO 2100
C
C BEGIN MAIN PLOT
C
2000  CONTINUE
      IF(IPASS.EQ.0) THEN
        CALL GRSTRT(800,1)
        CALL WINDOW(0.,100.,0.,100.)
        IPASS=1
      ELSE
        CALL NEWPAG
      ENDIF
      CALL OPNSEG(1)
      CALL SCALE3(XX,XLEN,NUMNP,1)
      CALL SCALE3(YY,YLEN,NUMNP,1)
      CALL SCALE2(ZZ,ZLEN,NUMNP,1)
C FOLLOWING TO MAKE X AND Y AXIS SAME SCALE
      IF(XX(NUMNP+2).GT.YY(NUMNP+2)) THEN
        YY(NUMNP+2)=XX(NUMNP+2)
      ELSE
        XX(NUMNP+2)=YY(NUMNP+2)
      ENDIF
C END OF SAME SCALE
      XLL=XX(NUMNP+1)
      YLL=YY(NUMNP+1)
      XD=XX(NUMNP+2)
      YD=YY(NUMNP+2)
      XUR=XLL+XLEN*XD
      YUR=YLL+YLEN*YD
      ZF=ZZ(NUMNP+1)
      ZD=ZZ(NUMNP+2)
      IF(CHAR.EQ.'I') THEN
        ZF=-1300.
        ZD=200.
      ELSEIF(CHAR.EQ.'Z') THEN
        ZF=4.
        ZD=1.
      ELSEIF(CHAR.EQ.'H') THEN
        ZF=-300.
        ZD=300.
      ELSEIF(CHAR.EQ.'U') THEN
        ZF=-300.
        ZD=300.
      ELSEIF(CHAR.EQ.'B') THEN
        ZF=-2000.
        ZD=400.
C       ZD=100
C       ZF=-1400.
      ELSEIF(CHAR.EQ.'D') THEN
        ZF=-.8
        ZD=.20
      ELSEIF(CHAR.EQ.'T') THEN
        ZF=-50.
        ZD=100.
      ENDIF
      CALL TXTCLR(1)
      ILZ=ZLEN+1
      PXX=100.
      PYY=92.5
      RNUM=ZF+ZLEN*ZD
      DO 130 I=1,ILZ
        IZ=ILZ-I+2
        CALL MRKCLR(ICMAP(IZ))
C       CALL MARKER(PXX,PYY,9)
        CALL MARKER(REAL(PXX),REAL(PYY),0)
        CALL MOVE(REAL(PXX+5.),REAL(PYY-2.))
        CALL RNUMBR(REAL(RNUM),3,9)
        RNUM=RNUM-ZD
        PYY=PYY-4.
130   CONTINUE
      RNUM=ZF
      DO 110 I=1,NUMNP
        RZ=(ZZ(I)-ZF)/ZD
        IZ=RZ+2
        IF(IZ.LT.2) IZ=2
        IF(IZ.GT.16) IZ=16
        CALL MRKCLR(ICMAP(IZ))
        XXX=(XX(I)-XLL)/XD
        YYY=(YY(I)-YLL)/YD
        IF(KODE(I).EQ.0) THEN
C         CALL MARKER(XXX,YYY,9)
          CALL MARKER(REAL(XXX),REAL(YYY),0)
        ELSE
          CALL MARKER(REAL(XXX),REAL(YYY),1)
        ENDIF
110   CONTINUE
      IF(ILINE.NE.0) THEN
        IEC=2
        CALL LINCLR(IEC)
C       CALL DASHPT(1)
        DO 120 I=1,NUMEL
          XXX=(XX(KX(I,1))-XLL)/XD
          YYY=(YY(KX(I,1))-YLL)/YD
          CALL MOVE(REAL(XXX),REAL(YYY))
          DO 115 J=2,NNODE(I)
            XXX=(XX(KX(I,J))-XLL)/XD
            YYY=(YY(KX(I,J))-YLL)/YD
            CALL DRAW(REAL(XXX),REAL(YYY))
115       CONTINUE
          XXX=(XX(KX(I,1))-XLL)/XD
          YYY=(YY(KX(I,1))-YLL)/YD
          CALL DRAW(REAL(XXX),REAL(YYY))
          IEC=IEC+1
          IF(IEC.GT.14) IEC=2
          CALL LINCLR(IEC)
120     CONTINUE
        CALL MAKCUR
      ENDIF
C
C READ AND PLOT OUTLINE IN FILE OUTLINE DATA B
      CALL LINCLR(1)
      IREAD=0
123   READ(10,*,END=124) XOUT,YOUT,IREAD
      IF(XOUT.EQ.-99999.) THEN
        READ(10,1000) JUNK
1000  FORMAT(A80)
        GOTO 123
      ENDIF
      IF(XOUT.EQ.-99999.) GOTO 124
      XOUT=(XOUT*1000.-XLL)/XD
      YOUT=(YOUT*1000.-YLL)/YD
      IF(IREAD.EQ.1) THEN
        CALL MOVE(REAL(XOUT),REAL(YOUT))
      ELSE
        CALL DRAW(REAL(XOUT),REAL(YOUT))
      ENDIF
      IREAD=1
      GOTO 123
124   CONTINUE
      REWIND 10
C
      CALL CLOSEG
      CALL MOVE(0.,0.)
      CALL RNUMBR(0.,-1,1)
      CALL MAKCUR
C
C     ******** MAIN MENU ********************************************
C
2100  CONTINUE
      IF(IPASS.EQ.1) CALL MAKCUR
      WRITE(*,*) 'INPUT:'
      WRITE(*,*) '<S>ELECT,<D>ISPLAY,<C>HANGE,<E>LEMENT,<A>GAIN,<M>OVE,'
      WRITE(*,*) 'NEW <N>ODE,NEW E<L>EMENT,<R>EMOVE ELEMENT,<B>ACKUP'
      WRITE(*,*) 'MODIF<Y>,<V>OLUME VS TIME'
      WRITE(*,*) 'L<I>NE,<Q>UIT,<Z>OOM,<U>NZOOM'
      READ(*,1002) CHAR
      IF(CHAR.EQ.'u') THEN
        CHAR='U'
      ELSEIF(CHAR.EQ.'z') THEN
        CHAR='Z'
      ELSEIF(CHAR.EQ.'s') THEN
        CHAR='S'
      ELSEIF(CHAR.EQ.'v') THEN
        CHAR='V'
      ELSEIF(CHAR.EQ.'d') THEN
        CHAR='D'
      ELSEIF(CHAR.EQ.'c') THEN
        CHAR='C'
      ELSEIF(CHAR.EQ.'e') THEN
        CHAR='E'
      ELSEIF(CHAR.EQ.'a') THEN
        CHAR='A'
      ELSEIF(CHAR.EQ.'m') THEN
        CHAR='M'
      ELSEIF(CHAR.EQ.'n') THEN
        CHAR='N'
      ELSEIF(CHAR.EQ.'l') THEN
        CHAR='L'
      ELSEIF(CHAR.EQ.'r') THEN
        CHAR='R'
      ELSEIF(CHAR.EQ.'q') THEN
        CHAR='Q'
      ELSEIF(CHAR.EQ.'b') THEN
        CHAR='B'
      ELSEIF(CHAR.EQ.'i') THEN
        CHAR='I'
      ELSEIF(CHAR.EQ.'y') THEN
        CHAR='Y'
      ENDIF
1002  FORMAT(A1)
      IF(CHAR.EQ.'U') GOTO 2800
      IF(CHAR.EQ.'Z') GOTO 2900
      IF(CHAR.EQ.'S') GOTO 3000
      IF(CHAR.EQ.'V') GOTO 3050
      IF(CHAR.EQ.'D') GOTO 3100
      IF(CHAR.EQ.'C') GOTO 3200
      IF(CHAR.EQ.'E') GOTO 3300
      IF(CHAR.EQ.'A') GOTO 3400
      IF(CHAR.EQ.'M') GOTO 3500
      IF(CHAR.EQ.'N') GOTO 3600
      IF(CHAR.EQ.'L') GOTO 3700
      IF(CHAR.EQ.'R') GOTO 3800
      IF(CHAR.EQ.'Q') GOTO 3900
      IF(CHAR.EQ.'B') GOTO 4000
      IF(CHAR.EQ.'I') GOTO 4100
      IF(CHAR.EQ.'Y') GOTO 4200
      GOTO 2100
C
2800  CONTINUE
C UNZOOM
C     CALL ZOOM(REAL(XMIN),REAL(XMAX),REAL(YMIN),REAL(YMAX))
      GOTO 2100
C
2900  CONTINUE
C SELECT A REGION FOR ZOOM
      WRITE(*,*) 'SELECT REGION FOR ZOOM, LOWER LEFT, THEN UPPER RIGHT'
      CALL SETGIN(1)
      CALL SETINK(1)
      CALL SETRUB(1)
      CALL LOCATE(2,XA,YA,IDAT,IGOT)
      IF(IGOT.NE.2) GOTO 2100
      PXMIN=XA(1)
      PYMIN=YA(1)
      PXMAX=XA(2)
      PYMAX=YA(2)
C     CALL ZOOM(REAL(PXMIN),REAL(PXMAX),REAL(PYMIN),REAL(PYMAX))
C     CALL INUMBR(IERRNM(.TRUE.),10)
      GOTO 2100
C
3000  CONTINUE
C SELECT A REGION FOR MODIFICATION
      WRITE(*,*) 'SELECT A REGION FOR MODIFICATION'
      CALL SETGIN(1)
      CALL SETINK(1)
      CALL SETRUB(1)
      CALL LOCATE(4,XA,YA,IDAT,IGOT)
      ICHAR=IDAT(4)
      DO 20 I=1,4
      XXX=XLL+XA(I)*XD
      YYY=YLL+YA(I)*YD
      XCHECK(I)=XXX
      YCHECK(I)=YYY
20    CONTINUE
      ICOUNT=0
      DO 30 I=1,NUMNP
C FIND RETURNS 1 IF FOUND, 0 IF NOT
        IFOUND=IFIND(XX(I),YY(I),XCHECK,YCHECK)
        IF(IFOUND.EQ.1) THEN
          ICOUNT=ICOUNT+1
          IFP(ICOUNT)=I
          XXX=(XX(I)-XLL)/XD
          YYY=(YY(I)-YLL)/YD
          CALL MRKCLR(1)
          CALL MARKER(REAL(XXX),REAL(YYY),10)
        ENDIF
30    CONTINUE
      GOTO 2100
C
3050  CONTINUE
C PLOT VOLUME VS TIME
      XMIN=1.E30
      YMIN=XMIN
      XMAX=-XMIN
      YMAX=-YMIN
      DO 3051 I=1,NTSTEP
        XMAX=MAX(XMAX,TTIME(I))
        YMAX=MAX(YMAX,1E6*VVOL(I))
        XMIN=MIN(XMIN,TTIME(I))
        YMIN=MIN(YMIN,1E6*VVOL(I))
3051  CONTINUE
      IF(IPASS.EQ.0) CALL GRSTRT(800,1)
      CALL WINDOW(REAL(XMIN),REAL(XMAX),REAL(YMIN),REAL(YMAX))
      CALL MOVE(REAL(TTIME(1)),REAL(1E6*VVOL(I)))
      DO 3052 I=2,NTSTEP
        CALL DRAW(REAL(TTIME(I)),REAL(1E6*VVOL(I)))
3052  CONTINUE
      CALL GRSTOP
      IPASS=0
      GOTO 2100
C
3100  CONTINUE
C DISPLAY ALL FOR CHOSEN REGION
      WRITE(*,*) ' DISPLAY ALL FOR CHOSEN REGION'
      WRITE(*,1050) ICOUNT
      IF(ICOUNT.GT.100) THEN
        WRITE(*,*) 'THIS IS ALOT, DO YOU REALLY WANT THIS?'
        READ(*,1002) CHAR
        IF(CHAR.NE.'Y') GOTO 2100
      ENDIF
1050  FORMAT(' NODAL VALUES FOR ',I5,' NODES',/,
     &'   N  K      X         Y   SURFACE BED  THICK   ACC   FRACT',
     &' PSURF FLOW SLIDE')
      DO 40 I=1,ICOUNT
        J=IFP(I)
        THIK=HTICE(J)-DEPB(J)
        WRITE(*,1100) J, KODE(J), XX(J), YY(J), HTICE(J), DEPB(J),
     &                THIK, ADOT(J), FRACT(J), PSURF(J), FLOWA(J),
     &                SLDGB(J)
1100  FORMAT(1X,I4,I2,2G10.3,3F6.0,2F7.2,F6.0,F5.2,F6.3)
40    CONTINUE
      GOTO 2100
C
3200  CONTINUE
C CHANGE VARIOUS PHYSICAL PARAMETERS
      WRITE(*,*) ' CHANGE VARIOUS PHYSICAL PARAMETERS'
      WRITE(*,*) '<K>ODE,<H>TICE,A<D>OT,PS<U>RF,<P>ERCENT'
      WRITE(*,*) '<S>LIDING,<F>LOW,<B>ED,<C>LEAR CONST,<T>IME STEP'
      WRITE(*,*) 'S<N>OW LINE,SNO<W>LINE GRADIENTS,FLU<X>'
      READ(*,1002) CHAR
      IF(CHAR.EQ.'k') THEN
        CHAR='K'
      ELSEIF(CHAR.EQ.'h') THEN
        CHAR='H'
      ELSEIF(CHAR.EQ.'d') THEN
        CHAR='D'
      ELSEIF(CHAR.EQ.'u') THEN
        CHAR='U'
      ELSEIF(CHAR.EQ.'p') THEN
        CHAR='P'
      ELSEIF(CHAR.EQ.'f') THEN
        CHAR='F'
      ELSEIF(CHAR.EQ.'s') THEN
        CHAR='S'
      ELSEIF(CHAR.EQ.'b') THEN
        CHAR='B'
      ELSEIF(CHAR.EQ.'c') THEN
        CHAR='C'
      ELSEIF(CHAR.EQ.'t') THEN
        CHAR='T'
      ELSEIF(CHAR.EQ.'n') THEN
        CHAR='N'
      ELSEIF(CHAR.EQ.'w') THEN
        CHAR='W'
      ELSEIF(CHAR.EQ.'x') THEN
        CHAR='X'
      ENDIF
      IF(CHAR.EQ.'K') GOTO 3210
      IF(CHAR.EQ.'X') GOTO 3212
      IF(CHAR.EQ.'W') GOTO 3215
      IF(CHAR.EQ.'H') GOTO 3220
      IF(CHAR.EQ.'D') GOTO 3230
      IF(CHAR.EQ.'U') GOTO 3240
      IF(CHAR.EQ.'P') GOTO 3250
      IF(CHAR.EQ.'F') GOTO 3260
      IF(CHAR.EQ.'S') GOTO 3270
      IF(CHAR.EQ.'B') GOTO 3280
      IF(CHAR.EQ.'C') GOTO 3290
      IF(CHAR.EQ.'T') GOTO 3295
      IF(CHAR.EQ.'N') GOTO 3297
      GOTO 3200
C
3210  CONTINUE
      WRITE(*,*) 'INPUT NEW KODE FOR SELECTED NODES'
      READ(*,*) KODEN
      DO 3211 I=1,ICOUNT
        KODE(IFP(I))=KODEN
3211  CONTINUE
      GOTO 2100
C
3212  CONTINUE
      RFLUX=0.
      DO 3213 I=1,NUMGBC
        RFLUX=RFLUX+BFLUX(I)
3213  CONTINUE
      IF(NUMGBC.GT.0) RFLUX=RFLUX/REAL(NUMGBC)
      WRITE(*,*) 'INPUT MULTIPLYING FACTOR FOR FLUXES'
      WRITE(*,*) 'CURRENT AVERAGE FLUX IS',RFLUX
      READ(*,*) RFLUX
      DO 3214 I=1,NUMGBC
        BFLUX(I)=BFLUX(I)*RFLUX
3214  CONTINUE
      GOTO 2100
C
3215  CONTINUE
      WRITE(*,*) 'INPUT POLE SNOWLINE ELEVATION AND LATITUDE GRADIENT,
     &TNSL'
      WRITE(*,*) AMASS(7),AMASS(8),AMASS(9)
      READ(*,*) AMASS(7),AMASS(8),AMASS(9)
      DO 3216 I=1,ICOUNT
        IF(IDT(IFP(I)).GT.0) THEN
            ADOT(IFP(I))=AFUNCT(TIME, IDT(IFP(I)), AMASS, HTICE(IFP(I)),
     &                          BED(IFP(I)), SLOPN(IFP(I)),
     &                          XX(IFP(I)),YY(IFP(I)))
        ENDIF
3216  CONTINUE
      GOTO 2100
C
3220  CONTINUE
      WRITE(*,*) 'INPUT NEW HTICE FOR SELECTED NODES'
      READ(*,*) HTICEN
      DO 3221 I=1,ICOUNT
        HTICE(IFP(I))=MAX(HTICEN,PSURF(IFP(I)))
C       HTICE(IFP(I))=MAX(HTICEN,BED(IFP(I)))
        IF(HTICE(IFP(I)).LT.0.) HTICE(IFP(I))=0.
3221  CONTINUE
      GOTO 2100
C
3230  CONTINUE
      WRITE(*,*) 'INPUT NEW ADOT FOR SELECTED NODES (-999 USES FIT)'
      READ(*,*) ADOTN
      IDOT=0
      IF(ADOTN.EQ.-999) THEN
3231    CONTINUE
        WRITE(*,*) '0 - MULTIPLY, 1 - USER DEFINED'
        WRITE(*,*) 'LINEAR POLAR <= 2-A, 3-B, 4-C'
        WRITE(*,*) 'LINEAR MARITIME <= 5-A, 6-B, 7-C'
        WRITE(*,*) 'SLOPE DEPENDENT 8=>SM,9=>TM,10=>SX,11=>PX,12=>PC'
      WRITE(*,*) 'EXPONENTIAL 13=>SM,14=>TM,15=>SX,16=>PX,17=>PC,18=ANT'
        WRITE(*,*) '19=>GROSSWALD,20=>METEOROLOGICAL'
        READ(*,*) IDOT
        IF(IDOT.LT.0 .OR. IDOT.GT.20) GOTO 3231
        IF(IDOT.EQ.1) THEN
          WRITE(*,*) 'INPUT AMARG,APEAK,ADOME,HE,HP,HM'
          WRITE(*,331) (AMASS(I),I=1,6)
331       FORMAT(3F10.4,3F10.0)
          READ(*,*) (AMASS(I),I=1,6)
        ENDIF
        IF(IDOT.EQ.0) THEN
          WRITE(*,*) 'MULTIPLICATIVE FACTOR'
          READ(*,*) AMULT
        ENDIF
        IF(IDOT.EQ.19) THEN
          WRITE(*,*) 'INPUT POLE SNOWLINE, SNOWLINE GRADIENT, BALANCE'
          WRITE(*,*) 'GRADIENT',(AMASS(II),II=7,9)
          READ(*,*) AMASS(7),AMASS(8),AMASS(9)
        ENDIF
      ENDIF
      DO 3232 I=1,ICOUNT
        IF(ADOTN.NE.-999.) THEN
          IF(BED(IFP(I)).GT.-9999) ADOT(IFP(I))=ADOTN
          IDT(IFP(I))=0
        ELSE
          IF(IDOT.NE.0) THEN
            ADOT(IFP(I))=AFUNCT(TIME, IDOT, AMASS, HTICE(IFP(I)),
     &                          BED(IFP(I)), SLOPN(IFP(I)),
     &                          XX(IFP(I)),YY(IFP(I)))
            IDT(IFP(I))=IDOT
          ELSE
            IF(BED(IFP(I)).GT.-9999.) ADOT(IFP(I))=AMULT*ADOT(IFP(I))
            IDT(IFP(I))=0
          ENDIF
        ENDIF
3232  CONTINUE
      GOTO 2100
C
3240  CONTINUE
      WRITE(*,*) 'INPUT NEW PRESENT SURFACE FOR SELECTED NODES'
      READ(*,*) PSURFN
      DO 3241 I=1,ICOUNT
        PSURF(IFP(I))=PSURFN
3241  CONTINUE
      GOTO 2100
C
3250  CONTINUE
      WRITE(*,*) 'INPUT NEW FRACT (1 FOR ALL SLIDING, 0 FOR NONE)'
      READ(*,*) FRACTN
      DO 3251 I=1,ICOUNT
        FRACT(IFP(I))=FRACTN
3251  CONTINUE
      GOTO 2100
C
3260  CONTINUE
      WRITE(*,*) 'INPUT NEW FLOWA FOR SELECTED NODES'
      READ(*,*) FLOWAN
      DO 3261 I=1,ICOUNT
        FLOWA(IFP(I))=FLOWAN
3261  CONTINUE
      GOTO 2100
C
3270  CONTINUE
      WRITE(*,*) 'INPUT NEW SLDGB FOR SELECTED NODES'
      READ(*,*) SLDGBN
      DO 3271 I=1,ICOUNT
        SLDGB(IFP(I))=SLDGBN
3271  CONTINUE
      GOTO 2100
C
3280  CONTINUE
      WRITE(*,*) 'INPUT NEW BED FOR SELECTED NODES'
      READ(*,*) BEDN
      DO 3281 I=1,ICOUNT
        BED(IFP(I))=BEDN
3281  CONTINUE
      GOTO 2100
C
3290  CONTINUE
      WRITE(*,*) 'SET ALL CONST=0.00'
      DO 3291 I=1,NUMEL
C       CONST(I)=1.92E7
        CONST(I)=0.0
3291  CONTINUE
      GOTO 2100
C
3295  CONTINUE
      WRITE(*,*) 'CURRENT NDT,INTER,DT',NDT,INTER,DT
      READ(*,*) NDT,INTER,DT
      GOTO 2100
C
3297  CONTINUE
      WRITE(*,*)'WHICH SNOW LINE DO YOU WANT TO CHANGE (-1 CHANGES ALL)'
      WRITE(*,*) '13=>SM,14=>TM,15=>SX,16=>PX,17=>PC,18=>??'
      READ(*,*) IDOT
      IF(IDOT.GE.0) THEN
        WRITE(*,*) 'CURRENT SNOW LINE ',SNO(IDOT)
        READ(*,*) SNO(IDOT)
      ELSE
        WRITE(*,*) (I,SNO(I),I=13,18)
        READ(*,*) SNOCON
        DO 3298 I=13,18
          SNO(I)=SNOCON
3298    CONTINUE
      ENDIF
      DO 3299 I=1,NUMNP
        IF(IDT(I).GT.0) ADOT(I)=AFUNCT(TIME, IDT(I), AMASS,
     &                                 HTICE(I), BED(I),
     &                                 SLOPN(I),XX(I),YY(I))
3299  CONTINUE
      GOTO 2100
C
3300  CONTINUE
C TOGGLE ELEMENT DISPLAY ON OR OFF
      WRITE(*,*) ' TOGGLE ELEMENT DISPLAY ON OR OFF'
      IF(ILINE.EQ.0) THEN
        ILINE=1
        WRITE(*,*) 'NOW ON'
      ELSE
        ILINE=0
        WRITE(*,*) 'NOW OFF'
      ENDIF
      GOTO 2100
C
3400  CONTINUE
C REDRAW PICTURE WITH NEW MAIN VARIABLE
      IF(IPASS.EQ.1) CALL DELSEG(1)
      WRITE(*,*) 'INPUT:'
      WRITE(*,*) '<K>ODE, <H>TICE, A<D>OT, PS<U>RF, <P>ERCENT'
      WRITE(*,*) '<S>LIDING, <F>LOW, <B>ED, D<I>FF, <T>HICK'
      WRITE(*,*) 'ADOT<Z>ONE, S<L>OPE'
      READ(*,1002) CHAR
      IF(CHAR.EQ.'k') THEN
        CHAR='K'
      ELSEIF(CHAR.EQ.'h') THEN
        CHAR='H'
      ELSEIF(CHAR.EQ.'t') THEN
        CHAR='T'
      ELSEIF(CHAR.EQ.'d') THEN
        CHAR='D'
      ELSEIF(CHAR.EQ.'z') THEN
        CHAR='Z'
      ELSEIF(CHAR.EQ.'u') THEN
        CHAR='U'
      ELSEIF(CHAR.EQ.'p') THEN
        CHAR='P'
      ELSEIF(CHAR.EQ.'f') THEN
        CHAR='F'
      ELSEIF(CHAR.EQ.'s') THEN
        CHAR='S'
      ELSEIF(CHAR.EQ.'b') THEN
        CHAR='B'
      ELSEIF(CHAR.EQ.'l') THEN
        CHAR='L'
      ELSEIF(CHAR.EQ.'i') THEN
        CHAR='I'
      ELSEIF(CHAR.EQ.'u') THEN
        CHAR='U'
      ENDIF
      IF(CHAR.EQ.'K') GOTO 3410
      IF(CHAR.EQ.'H') GOTO 3420
      IF(CHAR.EQ.'T') GOTO 3425
      IF(CHAR.EQ.'D') GOTO 3430
      IF(CHAR.EQ.'Z') GOTO 3435
      IF(CHAR.EQ.'U') GOTO 3440
      IF(CHAR.EQ.'P') GOTO 3450
      IF(CHAR.EQ.'F') GOTO 3460
      IF(CHAR.EQ.'S') GOTO 3470
      IF(CHAR.EQ.'B') GOTO 3480
      IF(CHAR.EQ.'L') GOTO 3485
      IF(CHAR.EQ.'I') GOTO 3490
      IF(CHAR.EQ.'U') GOTO 3495
      GOTO 3400
C
3410  CONTINUE
      WRITE(*,*) 'KODE '
      DO 3411 I=1,NUMNP
        ZZ(I)=KODE(I)
3411  CONTINUE
      GOTO 2000
C
3420  CONTINUE
      WRITE(*,*) 'HTICE '
      DO 3421 I=1,NUMNP
        ZZ(I)=HTICE(I)
3421  CONTINUE
      GOTO 2000
C
3425  CONTINUE
      WRITE(*,*) 'THICKNESS'
      DO 3426 I=1,NUMNP
        ZZ(I)=HTICE(I)-DEPB(I)
        IF(DEPB(I).LT.0.) THEN
          FLOT=-DEPB(I)*1.03/.917
          IF(HTICE(I).LT.FLOT) THEN
            ZZ(I)=0.
          ENDIF
        ENDIF
3426  CONTINUE
      GOTO 2000
C
3430  CONTINUE
      WRITE(*,*) 'ADOT '
      DO 3431 I=1,NUMNP
        ZZ(I)=ADOT(I)
3431  CONTINUE
      GOTO 2000
C
3435  CONTINUE
      WRITE(*,*) 'ADOT ZONE '
      DO 3436 I=1,NUMNP
        ZZ(I)=IDT(I)
3436  CONTINUE
      GOTO 2000
C
3440  CONTINUE
      WRITE(*,*) 'PSURF '
      DO 3441 I=1,NUMNP
        ZZ(I)=PSURF(I)
3441  CONTINUE
      GOTO 2000
C
3450  CONTINUE
      WRITE(*,*) 'FRACT '
      DO 3451 I=1,NUMNP
        ZZ(I)=FRACT(I)
3451  CONTINUE
      GOTO 2000
C
3460  CONTINUE
      WRITE(*,*) 'FLOWA '
      DO 3461 I=1,NUMNP
        ZZ(I)=FLOWA(I)
3461  CONTINUE
      GOTO 2000
C
3470  CONTINUE
      WRITE(*,*) 'SLDGB '
      DO 3471 I=1,NUMNP
        ZZ(I)=SLDGB(I)
3471  CONTINUE
      GOTO 2000
C
3480  CONTINUE
      WRITE(*,*) 'BED '
      DO 3481 I=1,NUMNP
        ZZ(I)=BED(I)
3481  CONTINUE
      GOTO 2000
C
3485  CONTINUE
      WRITE(*,*) 'SLOPE '
      DO 3486 I=1,NUMNP
        ZZ(I)=SLOPN(I)
3486  CONTINUE
      GOTO 2000
C
3490  CONTINUE
      WRITE(*,*) 'DIFFERENCE '
      DO 3491 I=1,NUMNP
C       ZZ(I)=HTICE(I)-HFIT(I)
        ZZ(I)=HTICE(I)-PSURF(I)
3491  CONTINUE
      GOTO 2000
C
3495  CONTINUE
      WRITE(*,*) 'PRESENT SURFACE '
      DO 3496 I=1,NUMNP
        ZZ(I)=PSURF(I)
3496  CONTINUE
      GOTO 2000
C
3500  CONTINUE
C CHANGE XY COORDINATE
      CALL LOCATE(1,XA,YA,IDAT,IGOT)
      ICHAR=IDAT(1)
      XXX=XLL+XA(1)*XD
      YYY=YLL+YA(1)*YD
      DMIN=1.E30
      DO 3510 I=1,NUMNP
        DIST=(XXX-XX(I))**2+(YYY-YY(I))**2
        IF(DIST.LT.DMIN) THEN
          DMIN=DIST
          IFOUND=I
        ENDIF
3510  CONTINUE
      XXX=(XX(IFOUND)-XLL)/XD
      YYY=(YY(IFOUND)-YLL)/YD
      CALL MRKCLR(1)
C     CALL MARKER(XXX,YYY,9)
      CALL MARKER(REAL(XXX),REAL(YYY),0)
      CALL LOCATE(1,XA,YA,IDAT,IGOT)
      ICHAR=IDAT(1)
      CALL MRKCLR(2)
      CALL MARKER(XA(1),YA(1),10)
      CALL MRKCLR(0)
C     CALL MARKER(REAL(XXX),REAL(YYY),9)
      CALL MARKER(REAL(XXX),REAL(YYY),0)
      CALL MRKCLR(1)
      XXX=XLL+XA(1)*XD
      YYY=YLL+YA(1)*YD
      XX(IFOUND)=XXX
      YY(IFOUND)=YYY
      GOTO 2100
C
3600  CONTINUE
C ADD NEW NODE
      NUMNP=NUMNP+1
      CALL LOCATE(1,XA,YA,IDAT,IGOT)
      ICHAR=IDAT(1)
C     CALL MARKER(XA(1),YA(1),9)
      CALL MARKER(XA(1),YA(1),0)
      XX(NUMNP)=XLL+XA(1)*XD
      YY(NUMNP)=YLL+YA(1)*YD
      KODE(NUMNP)=0
      HTICE(NUMNP)=100.
      ADOT(NUMNP)=1.0
      FRACT(NUMNP)=0.0
      PSURF(NUMNP)=0.0
      BED(NUMNP)=0.0
      FLOWA(NUMNP)=2.0
      SLDGB(NUMNP)=.02
      GOTO 2100
C
3700  CONTINUE
C DEFINE NEW ELEMENT CONNECTIVITY
      NUMEL=NUMEL+1
      WRITE(*,*) 'INPUT NUMBER OF NODES FOR NEW ELEMENT ',NUMEL
      READ(*,*) NN
      NNODE(NUMEL)=NN
      CONST(NUMEL)=1.E6
      KX(NUMEL,4)=0
      WRITE(*,*) 'PICK ENDPOINTS OF NEW ELEMENT BOUNDARY '
      CALL LOCATE(NN+1,XA,YA,IDAT,IGOT)
      ICHAR=IDAT(NN+1)
      DO 3710 I=1,NN
        XXX=XLL+XA(I)*XD
        YYY=YLL+YA(I)*YD
        XCHECK(I)=XXX
        YCHECK(I)=YYY
3710  CONTINUE
      DO 3720 I=1,NN
        DMIN=1.E30
        KX(NUMEL,I)=0
        DO 3720 JJ=1,NUMNP
          DIST=(XCHECK(I)-XX(JJ))**2+(YCHECK(I)-YY(JJ))**2
          IF(DIST.LT.DMIN) THEN
            DMIN=DIST
            KX(NUMEL,I)=JJ
          ENDIF
3720  CONTINUE
      GOTO 2100
C
3800  CONTINUE
      IFOUND=0
      DO 3820 IC=1,ICOUNT
        DO 3810 I=1,NUMEL
          IF(KX(I,4).EQ.0) THEN
            NN=3
          ELSE
            NN=4
          ENDIF
          DO 3809 N=1,NN
            IF(KX(I,N).EQ.IFP(IC)) THEN
C             FOUND ELEMENT CONTAINING NODE, DELETE
              IFOUND=IFOUND+1
              IFD(IFOUND)=I
              GOTO 3810
            ENDIF
3809      CONTINUE
3810    CONTINUE
3820  CONTINUE
      CALL SORTIX(IFOUND,IFD)
      CALL ELMDUP(IFOUND,IFD)
      DO 3830 I=1,IFOUND
        DO 3825 J=IFD(I),NUMEL-1
          DO 3824 K=1,4
            KX(J,K)=KX(J+1,K)
            CONST(J)=CONST(J+1)
3824      CONTINUE
3825    CONTINUE
        NUMEL=NUMEL-1
3830  CONTINUE
      CALL SORTIX(ICOUNT,IFP)
      DO 3833 J=1,ICOUNT
      DO 3833 N=1,NUMEL
      DO 3833 K=1,4
        IF(KX(N,K).GE.IFP(J)) KX(N,K)=KX(N,K)-1
3833  CONTINUE
      DO 3840 I=1,ICOUNT
        DO 3835 J=IFP(I),NUMNP-1
          XX(J)=XX(J+1)
          YY(J)=YY(J+1)
          HTICE(J)=HTICE(J+1)
          ADOT(J)=ADOT(J+1)
          BED(J)=BED(J+1)
          FRACT(J)=FRACT(J+1)
          PSURF(J)=PSURF(J+1)
          FLOWA(J)=FLOWA(J+1)
          SLDGB(J)=SLDGB(J+1)
          THICK(J)=THICK(J+1)
          KODE(J)=KODE(J+1)
          IDT(J)=IDT(J+1)
          ADOT(J)=ADOT(J+1)
3835    CONTINUE
        NUMNP=NUMNP-1
3840  CONTINUE
      GOTO 2100
C
3900  CONTINUE
C QUIT
      GOTO 999
C
4000  CONTINUE
C BACKUP
      REWIND 21
      WRITE(21,100) HED, NUMNP, NUMEL, NUMCOL, NUMLEV, NUMGBC, NDT,
     &              INTER, DT
      DO 4010 NUM=1,NUMNP
        IF(IDT(NUM).EQ.0) THEN
          WRITE(21,200) NUM, KODE(NUM), XX(NUM), YY(NUM), HTICE(NUM),
     &                  ADOT(NUM), FRACT(NUM), PSURF(NUM), BED(NUM),
     &                  FLOWA(NUM), SLDGB(NUM)
        ELSE
          ADDT=-100-IDT(NUM)
          WRITE(21,200) NUM, KODE(NUM), XX(NUM), YY(NUM), HTICE(NUM),
     &                  ADDT, FRACT(NUM), PSURF(NUM), BED(NUM),
     &                  FLOWA(NUM), SLDGB(NUM)
        ENDIF
C       THICK(NUM)=HTICE(NUM)-BED(NUM)
        THICK(NUM)=HTICE(NUM)-DEPB(NUM)
        ZZ(NUM)=HTICE(NUM)
4010  CONTINUE
      DO 4020 I=1,NUMEL
        WRITE(21,300) I, (KX(I,II),II=1,4),CONST(I)
4020  CONTINUE
      DO 4030 N=1,NUMGBC
        WRITE(21,310) IBFLUX(N,1),IBFLUX(N,2),BFLUX(N)
4030  CONTINUE
      GOTO 2100
C
C DEFINE LINE
4100  CONTINUE
      NP=0
4105  CONTINUE
      CALL LOCATE(1,XA,YA,IDAT,IGOT)
      ICHAR=IDAT(1)
      IF(ICHAR.EQ.113 .OR. ICHAR.EQ.81) GOTO 4120
      XXX=XLL+XA(1)*XD
      YYY=YLL+YA(1)*YD
      DMIN=1.E30
      DO 4110 I=1,NUMNP
        DIST=(XXX-XX(I))**2+(YYY-YY(I))**2
        IF(DIST.LT.DMIN) THEN
          DMIN=DIST
          IFOUND=I
        ENDIF
4110  CONTINUE
      CALL MRKCLR(1)
      XXX=(XX(IFOUND)-XLL)/XD
      YYY=(YY(IFOUND)-YLL)/YD
      CALL MARKER(REAL(XXX),REAL(YYY),10)
      NLINE(NP+1)=IFOUND
      NP=NP+1
      GOTO 4105
4120  CONTINUE
      WRITE(20,1010) NP
      WRITE(20,1010) (NLINE(IP),IP=1,NP)
1010  FORMAT(16I5)
      GOTO 2100
C
4200  CONTINUE
C FOLLOWING IS REFERENCE SUURFACE FOR ADJUSTMENTS
      RSURF=-500.
C CHANGE VARIOUS PHYSICAL PARAMETERS
      WRITE(*,*) ' MODIFY VARIOUS PHYSICAL PARAMETERS'
      WRITE(*,*) 'A<D>OT, <P>ERCENT, <S>LIDING, <F>LOW'
      READ(*,1002) CHAR
      IF(CHAR.EQ.'d') THEN
        CHAR='D'
      ELSEIF(CHAR.EQ.'p') THEN
        CHAR='P'
      ELSEIF(CHAR.EQ.'s') THEN
        CHAR='S'
      ELSEIF(CHAR.EQ.'f') THEN
        CHAR='F'
      ENDIF
      WRITE(*,*) 'INPUT RATIO FACTOR'
      READ(*,*) RFACT
      IF(CHAR.EQ.'D') GOTO 4210
      IF(CHAR.EQ.'P') GOTO 4210
      IF(CHAR.EQ.'F') GOTO 4210
      IF(CHAR.EQ.'S') GOTO 4210
      GOTO 4200
C
4210  CONTINUE
      DO 4215 I=1,ICOUNT
        J=IFP(I)
        IF(KODE(J).EQ.1) THEN
          RATIOD=0.
          RATIOF=1.
        ELSE
          RATIOD=1.-(HTICE(J)-RSURF)/(PSURF(J)-RSURF)
          RATIOD=RATIOD*RFACT
          RATIOF=RFACT*(PSURF(J)-RSURF)/(HTICE(J)-RSURF)
        ENDIF
        IF(CHAR.EQ.'D') THEN
          IF(CHAR.EQ.'D') ADOT(J)=ADOT(J)+RATIOD
        ELSEIF(CHAR.EQ.'P') THEN
          FRACT(J)=FRACT(J)-FRACT(J)*RATIOD
          IF(FRACT(J).GT.1.) FRACT(J)=1.
          IF(FRACT(J).LT.0.) FRACT(J)=0.
        ELSEIF(CHAR.EQ.'F') THEN
          FLOWA(J)=FLOWA(J)*RATIOF
          IF(FLOWA(J).GT.10.) FLOWA(J)=10.
          IF(FLOWA(J).LT.FLOWMIN) FLOWA(J)=FLOWMIN
        ELSEIF(CHAR.EQ.'S') THEN
          SLDGB(J)=SLDGB(J)*RATIOF
          IF(SLDGB(J).GT..20) SLDGB(J)=.20
          IF(SLDGB(J).LT.SLDGBMIN) SLDGB(J)=SLDGBMIN
        ENDIF
4215  CONTINUE
      GOTO 2100
C
999   CONTINUE
      IF(IPASS.EQ.1) CALL DELSEG(-1)
      ISKIP=0
      IF(ISKIP.EQ.0) GOTO 9999
      REWIND 21
      WRITE(21,100) HED, NUMNP, NUMEL, NUMCOL, NUMLEV, NUMGBC, NDT,
     &              INTER, DT
100   FORMAT(A80,/,7I5,F8.0)
      DO 101 NUM=1,NUMNP
        IF(IDT(NUM).EQ.0) THEN
          WRITE(21,200) NUM, KODE(NUM), XX(NUM), YY(NUM), HTICE(NUM),
     &                  ADOT(NUM), FRACT(NUM), PSURF(NUM), BED(NUM),
     &                  FLOWA(NUM), SLDGB(NUM)
        ELSE
          ADDT=-100-IDT(NUM)
          WRITE(21,200) NUM, KODE(NUM), XX(NUM), YY(NUM), HTICE(NUM),
     &                  ADDT, FRACT(NUM), PSURF(NUM), BED(NUM),
     &                  FLOWA(NUM), SLDGB(NUM)
        ENDIF
200   FORMAT(2I4,2E12.5,F10.2,F7.2,F9.2,F10.3,F10.1,F10.5,F10.5)
        THICK(NUM)=HTICE(NUM)-DEPB(NUM)
        ZZ(NUM)=HTICE(NUM)
101   CONTINUE
      DO 121 I=1,NUMEL
        WRITE(21,300) I, (KX(I,II),II=1,4),CONST(I)
121   CONTINUE
300   FORMAT(5I5,1PE17.10)
      DO 131 N=1,NUMGBC
        WRITE(21,310) IBFLUX(N,1),IBFLUX(N,2),BFLUX(N)
131   CONTINUE
310   FORMAT(2I5,E13.6)
9999  CALL GRSTOP
      RETURN
      END
      SUBROUTINE PLOTSOL(NUMNP,XX,YY,HTICE,DEPB,KODE)
      PARAMETER(MMXX=29999)
      DIMENSION XX(NUMNP+2),YY(NUMNP+2),ZZ(MMXX)
      DIMENSION HTICE(NUMNP),DEPB(NUMNP),KODE(NUMNP)
      DIMENSION ICMAP(16)
      DATA ICMAP /0,5,11,4,12,6,13,2,8,7,9,3,10,14,15,1/
      RHOI=.917
      RHOW=1.092
      RATIO=RHOW/RHOI
      DO I=1,NUMNP
        ZZ(I)=HTICE(I)-DEPB(I)
        IF(DEPB(I).LT.0.) THEN
          FLOT=-DEPB(I)*RATIO
          IF(HTICE(I).LT.FLOT) THEN
            ZZ(I)=0.
          ENDIF
        ENDIF
      ENDDO
      XLEN=100.
      YLEN=100.
      ZLEN=14.
      CALL SCALE3(XX,XLEN,NUMNP,1)
      CALL SCALE3(YY,YLEN,NUMNP,1)
c     CALL SCALE2(ZZ,ZLEN,NUMNP,1)
C FOLLOWING TO MAKE X AND Y AXIS SAME SCALE
      IF(XX(NUMNP+2).GT.YY(NUMNP+2)) THEN
        YY(NUMNP+2)=XX(NUMNP+2)
      ELSE
        XX(NUMNP+2)=YY(NUMNP+2)
      ENDIF
C END OF SAME SCALE
      XLL=XX(NUMNP+1)
      YLL=YY(NUMNP+1)
      XD=XX(NUMNP+2)
      YD=YY(NUMNP+2)
      XUR=XLL+XLEN*XD
      YUR=YLL+YLEN*YD
c     ZF=ZZ(NUMNP+1)
c     ZD=ZZ(NUMNP+2)
      ZF=-100.
      ZD=200.
      DO I=1,NUMNP
        RZ=(ZZ(I)-ZF)/ZD
        IZ=RZ+2
        IF(IZ.LT.2) IZ=2
        IF(IZ.GT.16) IZ=16
        XXX=(XX(I)-XLL)/XD
        YYY=(YY(I)-YLL)/YD
        IF(ZZ(I).GT.1.) THEN
          CALL MRKCLR(ICMAP(IZ))
          CALL MARKER(REAL(XXX),REAL(YYY),1)
        ELSE
          CALL MRKCLR(0)
          CALL MARKER(REAL(XXX),REAL(YYY),1)
        ENDIF
      ENDDO
      END
      SUBROUTINE DRAWOUT
C
C READ AND PLOT OUTLINE IN FILE OUTLINE DATA B
      REWIND 10
      CALL LINCLR(1)
      IREAD=0
123   CONTINUE
      READ(10,*,END=124) XOUT,YOUT,IREAD
      IF(XOUT.EQ.-99999.) THEN
        READ(10,1000) JUNK
1000  FORMAT(A80)
        GOTO 123
      ENDIF
      IF(XOUT.EQ.-99999.) GOTO 124
      XOUT=(XOUT*1000.-XLL)/XD
      YOUT=(YOUT*1000.-YLL)/YD
      IF(IREAD.EQ.1) THEN
        CALL MOVE(REAL(XOUT),REAL(YOUT))
      ELSE
        CALL DRAW(REAL(XOUT),REAL(YOUT))
      ENDIF
      IREAD=1
      GOTO 123
124   CONTINUE
      REWIND 10
c
      END
