	PROGRAM MAIN
C	MAIN PROGRAM FOR LITHOSPHERIC FLEXURE SOLUTION. LOAD IS ASSUMED TO
C	BE CONVOLVED WITH GREENS FUNCTION (UNIT LOAD) RESPONSE IN ORDER TO
C	CALCULATE THE NET EFFECT OF THE LOAD. IF DESIRED, INCLUSION OF
C       SEDIMENTS AS PART OF THE RESTORING FORCE EFFECT CAN BE CONTROLLED
C       VIA DENSITY ASSIGNED TO DENSITY(2). EACH LOAD ELEMENT CAN HAVE
C       SPECIFIED DENSITY AND SEVERAL LOADINGS EVENTS CAN BE INCORPORATED.
C       SEQUENCE OF LOADING EVENTS IS ASSUMED IN STRATIGRAPHIC SUCCESSION.
C       UP TO 10 SPECIFIED LOADINGS EVENTS OR SEQUENCES CAN BE USED.
C	*****************************************************************
	INTEGER NUMLOAD,POSIT,NUMEVENTS,NUMNODES,NODELOC,HID_POSIT
	REAL K, LAMBDA, LOAD, NODEPOS
	COMMON /PARAM/ NUMEVENTS,NUMLOAD(10),WIDTH_INT,FLEX_RIG,LAMBDA,
     1     K,DENSITY(3)
	COMMON /BASE/ NUMNODES,NODELOC(20),NODEPOS(20),BASELINE(201)
	COMMON /MOD1/ LOAD(10,50,2),POSIT(10,50),GREENF(401),
     1         DEFLECT(201)
	COMMON /MOD2/ NUM_HID, HID_LOAD(50,2), HID_POSIT(50)
C	-----------------------------------------------------------------
      WRITE (*, *) 'DO YOU WANT TO CREATE A NEW INPUT FILE (Y,N)'
      READ (*, 100) RESP
      IF (RESP .EQ. 'Y') THEN
         WRITE (*, *) 'NAME OF INPUT DATA FILE (20 CHARACTER MAX)'
         READ (*, 101) INFIL
         OPEN(8, FILE='INFIL')
         CALL INPUT (INFIL)
      ELSE
         WRITE (*, *) 'NAME OF INPUT DATA FILE'
         READ (*, 101) INFIL
         OPEN(8, FILE='INFIL')
      ENDIF
      WRITE (*, *) 'NAME OF OUTPUT FILE (20 CHARACTERS MAX)'
      READ (*, 101) OUTFIL
      OPEN(7, FILE='OUTFIL')
      REWIND 8
      READ (8, 101) MODEL
      READ (8, 102) NUMNODES
      DO 1 I = 1, NUMNODES
         READ (8, 103) NODELOC(I), NODEPOS(I)
    1 CONTINUE
      READ (8, 102) NUMEVENTS
      READ (8, 104) WIDTH_INT, FLEX_RIG
      DO 3 J = 1, NUMEVENTS
         READ (8, 102) NUMLOAD(J)
         DO 2 I = 1, NUMLOAD(J)
            READ (8, 105) POSIT(J,I), LOAD(J,I,1), LOAD(J,I,2)
    2    CONTINUE
    3 CONTINUE
      READ (8, 102) NUM_HID
      DO 4 I = 1, NUM_HID
         READ (8, 105) HID_POSIT(I), HID_LOAD(I,1), HID_LOAD(I,2)
    4 CONTINUE
      CALL INITIALIZE
      CALL UNITLOAD
      CALL CONVOLVE
      IF (NUM_HID .NE. 0) CALL HIDDEN
      CALL OUTPUT (MODEL)
  100 FORMAT(A1)
  101 FORMAT(A20)
  102   FORMAT(I5)
  103   FORMAT(I5,10X, F10.0)
  104 FORMAT(F10.0,10X,E10.2)
  105 FORMAT(I5,10X,F10.0,10X,F10.0)
      CLOSE(8)
      CLOSE(7)
      STOP 
      END
C	*******************************************************************
C	*******************************************************************
C	PROGRAM INPUT.F
C	PRODUCES INPUT DATA FILE WITH LOAD DISTRIBUTION AND FLEXURAL
C	PARAMETERS
	SUBROUTINE INPUT(FILNM)
C	------------------------------------------------------------------
	CHARACTER FILNM*20, MODEL*20
	INTEGER NUMLOAD,POSIT,NUMEVENTS,NUMNODES,NODELOC,HID_POSIT
	REAL K, LAMBDA, LOAD, NODEPOS
	COMMON /PARAM/ NUMEVENTS,NUMLOAD(10),WIDTH_INT,FLEX_RIG,LAMBDA,
     1     K,DENSITY(3)
	COMMON /BASE/ NUMNODES,NODELOC(20),NODEPOS(20),BASELINE(201)
	COMMON /MOD1/ LOAD(10,50,2), POSIT(10,50), GREENF(401),
     1     DEFLECT(201)
	COMMON /MOD2/ NUM_HID, HID_LOAD(50,2), HID_POSIT(50)
C	------------------------------------------------------------------
      WRITE (*, *) 'NAME OF MODEL RUN (20 CHARACTERS MAXIMUM)'
      READ (*, 100) MODEL
      WRITE (*, *) 'WIDTH OF LOADING ELEMENT (METERS)'
      READ (*, 102) WIDTH_INT
      WRITE (*, *) 'VALUE FOR FLEXURAL RIGIDITY (Nm)'
      READ (*, 103) FLEX_RIG
      WRITE (*, *) 'NUMBER OF NODES DESCRIBING BASELINE POSITION'
      READ (*, 101) NUMNODES
      DO 1 I = 1, NUMNODES
         WRITE (*, *) 'FOR BASELINE NODE', I
         WRITE (*, *) 'ENTER POSITION (NODE NUMBER), '
         WRITE (*, *) 'AND ELEVATION (METERS + UP)'
         READ (*, 104) NODELOC(I), NODEPOS(I)
    1 CONTINUE
      WRITE (*, *) 'NUMBER OF LOADING EVENTS'
      READ (*, 101) NUMEVENTS
      DO 3 J = 1, NUMEVENTS
         WRITE (*, *) 'NUMBER OF LOADING ELEMENTS FOR EVENT', J
         READ (*, 101) NUMLOAD(J)
         DO 2 I = 1, NUMLOAD(J)
            WRITE (*, *) 'FOR LOADING EVENT', J
            WRITE (*, *) 'ENTER POSITION (NODE NUMBER), '
            WRITE (*, *) 'HEIGHT (METERS) OF LOAD ELEMENT'
            WRITE (*, *) 'AND DENSITY (KG/M**3) OF NODE', I
            READ (*, 105) POSIT(J,I), LOAD(J,I,1), LOAD(J,I,2)
    2    CONTINUE
    3 CONTINUE
      WRITE (*, *) ' NUMBER OF HIDDEN LOAD ELEMENTS'
      READ (*, 101) NUM_HID
      DO 4 I = 1, NUM_HID
         WRITE (*, *) 'ENTER POSITION (NODE NUMBER), '
         WRITE (*, *) 'THICKNESS (METERS) OF LOAD ELEMENT'
         WRITE (*, *) 'AND DENSITY (KG/M**3) OF HIDDEN LOAD', I
         READ (*, 105) HID_POSIT(I), HID_LOAD(I,1), HID_LOAD(I,2)
    4 CONTINUE
      WRITE (8, 100) MODEL
      WRITE (8, 101) NUMNODES
      DO 5 I = 1, NUMNODES
         WRITE (8, 106) NODELOC(I), NODEPOS(I)
    5 CONTINUE
      WRITE (8, 101) NUMEVENTS
      WRITE (8, 107) WIDTH_INT, FLEX_RIG
      DO 7 J = 1, NUMEVENTS
         WRITE (8, 101) NUMLOAD(J)
         DO 6 I = 1, NUMLOAD(J)
            WRITE (8, 108) POSIT(J,I), LOAD(J,I,1), LOAD(J,I,2)
    6    CONTINUE
    7 CONTINUE
      WRITE (8, 101) NUM_HID
      DO 8 I = 1, NUM_HID
         WRITE (8, 108) HID_POSIT(I), HID_LOAD(I,1), HID_LOAD(I,2)
    8 CONTINUE
  100 FORMAT(A20)
  101 FORMAT(I5)
  102 FORMAT(F10.0)
  103 FORMAT(E10.2)
  104   FORMAT(I5, F10.0)
  105 FORMAT(I5, F10.0, F10.0)
  106   FORMAT(I5, 10X, F10.0)
  107 FORMAT(F10.0,10X,E10.2)
  108 FORMAT(I5,10X,F10.0,10X,F10.0)
      RETURN 
      END
C	*******************************************************************
C	*******************************************************************
C	PROGRAM INITIALIZE.F
C	SETS UP LOAD DISTRIBUTION, VARIOUS FLEXURAL PARAMETERS, AND
C	PREPARES FOR UNITLOAD AND CONVOLVE
	SUBROUTINE INITIALIZE
	INTEGER NUMLOAD,POSIT,NUMEVENTS,NUMNODES,NODELOC,HID_POSIT
	REAL K,LAMBDA,LOAD,NODEPOS,GRAV,NODEDIF
	COMMON /PARAM/ NUMEVENTS,NUMLOAD(10),WIDTH_INT,FLEX_RIG,LAMBDA,
     1     K,DENSITY(3)
	COMMON /BASE/ NUMNODES,NODELOC(20),NODEPOS(20),BASELINE(201)
	COMMON /MOD1/ LOAD(10,50,2), POSIT(10,50), GREENF(401),
     1     DEFLECT(201)
	COMMON /MOD2/ NUM_HID, HID_LOAD(50,2), HID_POSIT(50)
      DENSITY(1) = 3300.0
      DENSITY(2) = 0.0
      DENSITY(3) = 2800.0
      K = (DENSITY(3)-DENSITY(2))/(DENSITY(1)-DENSITY(2))
      GRAV = 9.8
      LAMBDA = ABS((DENSITY(1)-DENSITY(2))*GRAV/(4.0*FLEX_RIG))**0.25
      DO 1 I = 1, 201
         BASELINE(I) = 0.0
    1 CONTINUE
      IF (NUMNODES .GT. 1) THEN
         DO 3 I = 1, NUMNODES - 1
            II = I + 1
            DO 2 J = NODELOC(I), NODELOC(II)
               NODEDIF = (NODEPOS(II)-NODEPOS(I))/(NODELOC(II)-NODELOC(I
     #            ))
               NODEDIST = J - NODELOC(I)
               BASELINE(J) = NODEPOS(I) + NODEDIF*NODEDIST
    2       CONTINUE
    3    CONTINUE
      ENDIF
      RETURN 
      END
C	*******************************************************************
C	*******************************************************************
C	PROGRAM UNITLOAD.F
C	CALCULATES THE UNIT RESPONSE TO LOADING FOR FLEXURAL PARAMETERS
C	ASSUMED. RESULT IS USED IN CONVOLUTION PERFORMED IN CONVOLVE.F
	SUBROUTINE UNITLOAD
	INTEGER NUMLOAD,POSIT,NUMEVENTS,NUMNODES,NODELOC,HID_POSIT
	REAL K, LAMBDA, LOAD, NODEPOS, LAM_WID
	COMMON /PARAM/ NUMEVENTS,NUMLOAD(10),WIDTH_INT,FLEX_RIG,LAMBDA,
     1     K,DENSITY(3)
	COMMON /BASE/ NUMNODES,NODELOC(20),NODEPOS(20),BASELINE(201)
	COMMON /MOD1/ LOAD(10,50,2),POSIT(10,50),GREENF(401),
     1     DEFLECT(201)
	COMMON /MOD2/ NUM_HID, HID_LOAD(50,2), HID_POSIT(50)
      LAM_WID = LAMBDA*WIDTH_INT
      DO 1 I = 1, 200
         P = 201 - I
         GREENF(I) = K/2.*(EXP((-LAM_WID*(P+0.5)))*COS(LAM_WID*(P+0.5))-
     #      EXP((-LAM_WID*(P-0.5)))*COS(LAM_WID*(P-0.5)))
    1 CONTINUE
      GREENF(201) = (-K/2.)*(2-2*EXP((-LAM_WID/2.))*COS(LAM_WID/2.))
      DO 2 I = 202, 401
         P = I - 201
         GREENF(I) = K/2.*(EXP((-LAM_WID*(P+0.5)))*COS(LAM_WID*(P+0.5))-
     #      EXP((-LAM_WID*(P-0.5)))*COS(LAM_WID*(P-0.5)))
    2 CONTINUE
      RETURN 
      END
C	*******************************************************************
C	*******************************************************************
C	PROGRAM CONVOLVE.F
C	CALCULATES THE RESPONSE TO A SPECIFIED LOAD DISTRIBUTION BY THE
C	CONVOLUTION OF THE LOAD DISTRIBUTION WITH THE RESPONSE TO A UNIT
C	LOAD AS CALCULATED IN UNITLOAD.F
	SUBROUTINE CONVOLVE
	INTEGER NUMLOAD,POSIT,NUMEVENTS,NUMNODES,NODELOC,SA,HID_POSIT
	REAL K, LAMBDA, LOAD, NODEPOS
	COMMON /PARAM/ NUMEVENTS,NUMLOAD(10),WIDTH_INT,FLEX_RIG,LAMBDA,
     1     K,DENSITY(3)
	COMMON /BASE/ NUMNODES,NODELOC(20),NODEPOS(20),BASELINE(201)
	COMMON /MOD1/ LOAD(10,50,2),POSIT(10,50),GREENF(401),
     1     DEFLECT(201)
	COMMON /MOD2/ NUM_HID, HID_LOAD(50,2), HID_POSIT(50)
      DO 1 I = 1, 201
         DEFLECT(I) = 0.0
    1 CONTINUE
      DO 4 I = 1, 401
         DO 3 JJ = 1, NUMEVENTS
            DO 2 J = 1, NUMLOAD(JJ)
               SA = (-I) + POSIT(JJ,J) + 201
               TEMPFLEX = LOAD(JJ,J,1)*(LOAD(JJ,J,2)/2800.0)*GREENF(I)
               IF (SA .GE. 1) THEN
                  IF (SA .LE. 201) DEFLECT(SA) = DEFLECT(SA) + TEMPFLEX
               ENDIF
    2       CONTINUE
    3    CONTINUE
    4 CONTINUE
      RETURN 
      END
C	*******************************************************************
C	*******************************************************************
C	PROGRAM HIDDEN.F
C	CALCULATES THE RESPONSE TO A HIDDEN LOAD DISTRIBUTION BY THE
C	CONVOLUTION OF THE LOAD DISTRIBUTION WITH THE RESPONSE TO A UNIT
C	LOAD AS CALCULATED IN UNITLOAD.F
	SUBROUTINE HIDDEN
	INTEGER NUMLOAD,POSIT,NUMEVENTS,NUMNODES,NODELOC,SA,HID_POSIT
	REAL K, LAMBDA, LOAD, NODEPOS
	COMMON /PARAM/ NUMEVENTS,NUMLOAD(10),WIDTH_INT,FLEX_RIG,LAMBDA,
     1     K,DENSITY(3)
	COMMON /BASE/ NUMNODES,NODELOC(20),NODEPOS(20),BASELINE(201)
	COMMON /MOD1/ LOAD(10,50,2),POSIT(10,50),GREENF(401),
     1     DEFLECT(201)
	COMMON /MOD2/ NUM_HID, HID_LOAD(50,2), HID_POSIT(50)
      DO 2 I = 1, 401
         DO 1 J = 1, NUM_HID
            SA = (-I) + HID_POSIT(J) + 201
            HID_NET = HID_LOAD(J,2) - DENSITY(3)
            TEMPFLEX = HID_LOAD(J,1)*(HID_NET/DENSITY(3))*GREENF(I)
            IF (SA .GE. 1) THEN
               IF (SA .LE. 201) DEFLECT(SA) = DEFLECT(SA) + TEMPFLEX
            ENDIF
    1    CONTINUE
    2 CONTINUE
      RETURN 
      END
C	*******************************************************************
C	*******************************************************************
C	PROGRAM OUTPUT.F
C	WRITES RESULTS OF MODEL TO FILE (OUTFIL) FOR LATER USE IN
C	PLOTTING ROUTINES
	SUBROUTINE OUTPUT(MODEL)
	CHARACTER MODEL*20
	INTEGER NUMLOAD,POSIT,NUMEVENTS,NUMNODES,NODELOC,HID_POSIT
	REAL K, LAMBDA, LOAD, NODEPOS, DIF(201)
	COMMON /PARAM/ NUMEVENTS,NUMLOAD(10),WIDTH_INT,FLEX_RIG,LAMBDA,
     1     K,DENSITY(3)
	COMMON /BASE/ NUMNODES,NODELOC(20),NODEPOS(20),BASELINE(201)
	COMMON /MOD1/ LOAD(10,50,2),POSIT(10,50),GREENF(401),
     1    DEFLECT(201)
	COMMON /MOD2/ NUM_HID, HID_LOAD(50,2), HID_POSIT(50)
      WRITE (7, 100) MODEL
      WRITE (7, 101) NUMEVENTS
      WRITE (7, 102) WIDTH_INT, FLEX_RIG
      DO 1 I = 1, 201
         DIF(I) = DEFLECT(I) + BASELINE(I)
    1 CONTINUE
      DO 2 I = 1, 201, 3
         WRITE (7, 104) (II, BASELINE(II), II = I, I + 2)
    2 CONTINUE
      DO 3 I = 1, 201, 3
         WRITE (7, 104) (II, DEFLECT(II), II = I, I + 2)
    3 CONTINUE
      DO 4 I = 1, 201, 3
         WRITE (7, 104) (II, DIF(II), II = I, I + 2)
    4 CONTINUE
      DO 7 J = 1, NUMEVENTS
         WRITE (7, 101) NUMLOAD(J)
         DO 5 I = 1, NUMLOAD(J)
            WRITE (7, 103) J, I, POSIT(J,I), LOAD(J,I,1), LOAD(J,I,2)
            DIF(POSIT(J,I)) = DIF(POSIT(J,I)) + LOAD(J,I,1)
    5    CONTINUE
         DO 6 JJ = 1, 201, 3
            WRITE (7, 104) (II, DIF(II), II = JJ, JJ + 2)
    6    CONTINUE
    7 CONTINUE
  100 FORMAT(A20)
  101   FORMAT(I5)
  102 FORMAT(F10.0,10X,E10.2)
  103 FORMAT(I5,5X,I5,5X, I5,5X,F10.0,5X,F10.0)
  104 FORMAT(3(I5, F12.3,5x))
      RETURN 
      END
 
