	PROGRAM MAIN
C	MAIN PROGRAM FOR LITHOSPHERIC FLEXURE SOLUTION FOR A BROKEN PLATE.
C       LOAD IS ASSUMED TO BE REPRESENTED BY EQUAL WIDTH LOADING ELEMENTS
C       SPECIFIED DISTANCE FROM BROKEN EDGE OF PLATE.
C	INCLUSION OF SEDIMENTS AS PART OF THE RESTORING FORCE EFFECT IS
C       POSSIBLE BY CHOICE OF DENSITY ASSIGNED TO DENSITY (2).
*****************************************************************
	INTEGER NUMLOAD, POSIT, NUMEVENTS
	REAL LOAD, K, LAMBDA
	COMMON /PARAM/ NUMEVENTS,NUMLOAD(10), WIDTH_INT, FLEX_RIG, LAMBDA,
     1     K,DENSITY(3)
	COMMON /MOD1/LOAD(10,201,2),POSIT(10,201),GREENF(201),DEFLECT(201)
	CHARACTER INFIL*20, OUTFIL*20, RESP*1, MODEL*20
      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) NUMEVENTS
      READ (8, 103) WIDTH_INT, FLEX_RIG
      DO 2 J = 1, NUMEVENTS
         READ (8, 102) NUMLOAD(J)
         DO 1 I = 1, NUMLOAD(J)
            READ (8, 104) POSIT(J,I), LOAD(J,I,1), LOAD(J,I,2)
    1    CONTINUE
    2 CONTINUE
      CALL INITIALIZE
      CALL SEMI_INFINITE_BEAM
      CALL OUTPUT (MODEL)
  100 FORMAT(A1)
  101 FORMAT(A20)
  102   FORMAT(I5)
  103 FORMAT(F10.0,10X,E10.2)
  104 FORMAT(I5,10X,F10.0,10X,F10.0)
      CLOSE(8)
      CLOSE(7)
      STOP 
      END
C	*******************************************************************
C	PROGRAM INPUT.F
C	PRODUCES INPUT DATA FILE WITH LOAD DISTRIBUTION AND FLEXURAL
C	PARAMETERS
	SUBROUTINE INPUT(FILNM)
	CHARACTER FILNM*20, MODEL*20
	INTEGER NUMLOAD, POSIT, NUMEVENTS
	REAL LOAD
	COMMON /PARAM/ NUMEVENTS,NUMLOAD(10),WIDTH_INT,FLEX_RIG,LAMBDA,
     1     K,DENSITY(3)
	COMMON /MOD1/LOAD(10,201,2),POSIT(10,201),GREENF(201),DEFLECT(201)
      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 LOADING EVENTS'
      READ (*, 101) NUMEVENTS
      DO 2 J = 1, NUMEVENTS
         WRITE (*, *) 'NUMBER OF LOADING ELEMENTS FOR EVENT', J
         READ (*, 101) NUMLOAD(J)
         DO 1 I = 1, NUMLOAD(J)
            WRITE (*, *) 'FOR LOADING EVENT', J
            WRITE (*, *) 'ENTER POSITION (NODE NUMBER),'
            WRITE (*, *) 'HEIGHT (M) OF LOAD ELEMENT'
            WRITE (*, *) 'AND DENSITY (KG/M**3) OF NODE', I
            READ (*, *) POSIT(J,I), LOAD(J,I,1), LOAD(J,I,2)
    1    CONTINUE
    2 CONTINUE
      WRITE (8, 100) MODEL
      WRITE (8, 101) NUMEVENTS
      WRITE (8, 105) WIDTH_INT, FLEX_RIG
      DO 4 J = 1, NUMEVENTS
         WRITE (8, 101) NUMLOAD(J)
         DO 3 I = 1, NUMLOAD(J)
            WRITE (8, 106) POSIT(J,I), LOAD(J,I,1), LOAD(J,I,2)
    3    CONTINUE
    4 CONTINUE
  100 FORMAT(A20)
  101 FORMAT(I5)
  102 FORMAT(F10.0)
  103 FORMAT(E10.2)
  104 FORMAT(I5, F5.0, F10.0)
  105 FORMAT(F10.0,10X,E10.2)
  106 FORMAT(I5,10X,F10.0,10X,F10.0)
      RETURN 
      END
C	*******************************************************************
C	PROGRAM INITIALIZE.F
C	SETS UP LOAD DISTRIBUTION, VARIOUS FLEXURAL PARAMETERS, AND
C	PREPARES FOR SEMI_INFINITE_BEAM
	SUBROUTINE INITIALIZE
	INTEGER NUMLOAD, POSIT, NUMEVENTS
	REAL K, LAMBDA, LOAD, GRAV
	COMMON /PARAM/ NUMEVENTS,NUMLOAD(10), WIDTH_INT, FLEX_RIG, LAMBDA,
     1     K, DENSITY(3)
	COMMON /MOD1/LOAD(10,201,2),POSIT(10,201),GREENF(201),DEFLECT(201)
      DENSITY(1) = 3300.0
      DENSITY(2) = 0.0
      DENSITY(3) = 2800.0
      K = (DENSITY(1)-DENSITY(2))/(DENSITY(3)-DENSITY(2))
      GRAV = 9.8
      LAMBDA = ABS((DENSITY(1)-DENSITY(2))*GRAV/(4.0*FLEX_RIG))**0.25
      RETURN 
      END
C	*******************************************************************
C       PROGRAM SEMI_INFINITE_BEAM.F
C       CALCULATES DEFLECTION FOR A SEMI-INFINITE BEAM WITH A FREE END
C       HAVING A FIX REFERENCE FRAME AT THE END OF THE BEAM AND
C       MOVING THE POINT LOADING ALONG IT
	SUBROUTINE SEMI_INFINITE_BEAM
	INTEGER NUMLOAD, POSIT, NUMEVENTS
	REAL K, LAMBDA, LOAD, LAM, LAM1, LAM2, C, C1, LAM_WID, LAM0
	COMMON /PARAM/ NUMEVENTS,NUMLOAD(10), WIDTH_INT, FLEX_RIG, LAMBDA,
     1     K, DENSITY(3)
	COMMON /MOD1/LOAD(10,201,2),POSIT(10,201),GREENF(201),DEFLECT(201)
      DO 1 I = 1, 201
         DEFLECT(I) = 0.0
    1 CONTINUE
      DO 4 I = 1, 201
         DO 3 JJ = 1, NUMEVENTS
            DO 2 J = 1, NUMLOAD(JJ)
               C1 = POSIT(JJ,J)*WIDTH_INT - WIDTH_INT
               C = (C1*0.01745)/111300
               LAM_WID = LAMBDA*WIDTH_INT
               LAM0 = LAMBDA*C1
               LAM = LAMBDA*C
               LAM1 = LAMBDA*((WIDTH_INT*I-WIDTH_INT)-C1)
               LAM2 = LAMBDA*(I*WIDTH_INT-WIDTH_INT)
               GREENF(I) = -(((LAM_WID/(2.*K))*EXP((-ABS(LAM1))))*(COS(
     #            ABS(LAM1))+SIN(ABS(LAM1)))+((LAM_WID/K)*EXP((-LAM2))
     #            *EXP((-LAM0)))*(COS(LAM)*COS(LAM2)-0.5*COS(LAM2)*(SIN(
     #            LAM)-COS(LAM))+0.5*SIN(LAM2)*(SIN(LAM)-COS(LAM))))
               DEFLECT(I) = LOAD(JJ,J,1)*GREENF(I)*(LOAD(JJ,J,2)/2800)
     #             + DEFLECT(I)
    2       CONTINUE
    3    CONTINUE
    4 CONTINUE
      RETURN 
      END
C	*******************************************************************
C	PROGRAM OUTPUT.F
C	WRITES RESULTS OF MODEL TO FILE (OUTFIL) FOR LATER USE
	SUBROUTINE OUTPUT(MODEL)
	CHARACTER MODEL*20
	INTEGER NUMLOAD, POSIT, NUMEVENTS
	REAL K, LAMBDA, LOAD, DIF(201)
	COMMON /PARAM/ NUMEVENTS, NUMLOAD(10), WIDTH_INT,FLEX_RIG, LAMBDA,
     1     K, DENSITY(3)
	COMMON /MOD1/LOAD(10,201,2),POSIT(10,201),GREENF(201),DEFLECT(201)
      WRITE (7, 100) MODEL
      WRITE (7, 101) NUMEVENTS
      WRITE (7, 102) WIDTH_INT, FLEX_RIG
      DO 1 I = 1, 201
         DIF(I) = DEFLECT(I)
    1 CONTINUE
      DO 2 I = 1, 201, 3
         WRITE (7, 105) (II, DEFLECT(II), II = I, I + 2)
    2 CONTINUE
      DO 5 J = 1, NUMEVENTS
         DO 3 I = 1, NUMLOAD(J)
            DIF(POSIT(J,I)) = DIF(POSIT(J,I)) + LOAD(J,I,1)
    3    CONTINUE
         DO 4 JJ = 1, 201, 3
            WRITE (7, 105) (II, DIF(II), II = JJ, JJ + 2)
    4    CONTINUE
    5 CONTINUE
  100   FORMAT(A20)
  101   FORMAT(I5)
  102   FORMAT(F10.0,10X,E10.2)
  103   FORMAT(I5,10X,F10.0,10X,E10.2)
  104 FORMAT(I5,5X, I5,5X,F10.0,5X,F10.0)
  105 FORMAT(3(I5, F12.0,5x))
      RETURN 
      END
 
