!----------------------------------------------------------------------
write (*,*) ' '
write (*,*) '***************** COMCOT ****************'
write (*,*) '* *'
write (*,*) '* VERSION=',Version,' *'
write (*,*) '* *'
write (*,*) '********************************************'
!----------------------------------------------------------------------
!///// READ USER-DEFINED OPTIONS (NO. OF LAYERS, COORDINATE, GOVERN EQN.)
CALL SYSTEM_CLOCK(count1, count_rate1, count_max1)
CALL READ_CONFIG (LO,LA,TEND,T_INV,INI_SURF,WAVE_INFO, &
FAULT_INFO,LANDSLIDE_INFO,BCI_INFO, &
START_TYPE,START_TIME,BC_TYPE,OUT_OPT,JOB)
!///// READING BATH/TOPO DATA /////////////////////////////////////////
CALL READ_BATHYMETRY(LO,LA)
!///// READ INITIAL CONDITION /////////////////////////////////////////
IF (BC_TYPE.NE.3) THEN
CALL GET_INI_SURF(LO,LA,INI_SURF,WAVE_INFO,FAULT_INFO, &
LANDSLIDE_INFO)
ENDIF
!///// ADJUST BATHYMETRY, TIME STEP AND CALC DEPENDENT PARAMETERS /////
CALL ADJUST_BATHYMETRY (LO,LA)
CALL CR_CHECK (LO,LA)
CALL ALPHA_CALC (LO,LA)
!CALCULATE PARAMETERS USED FOR SPHERICAL COORD.
CALL SPHERICAL_PARAMETERS(LO,LA)
!///// SET UP BOUNDARY CONDITIONS /////////////////////////////////////
IF (BC_TYPE.EQ.1) CALL SPONGE_COEF (LO)
IF (BC_TYPE.EQ.2) CALL BC_WALL (LO,WAVE_INFO)
IF (BC_TYPE.EQ.3) CALL GET_BC_DATA (BCI_INFO,LO)
!/////DETERMINE STARTING TIME STEP # //////////////////////////////////
! WRITE(*,*) START_STEP
START_STEP = NINT(START_TIME/LO%DT)
ISTART = START_STEP
IF (START_TYPE .EQ. 0) ISTART = 1
!//////// DISPLAY INPUT PARAMETERS ////////////////////////////////////
CALL SCREEN_DISP (LO,LA,FAULT_INFO,LANDSLIDE_INFO,WAVE_INFO, &
INI_SURF,TEND,T_INV,KEND,IPRT,START_TYPE, &
START_TIME,BC_TYPE,BCI_INFO,OUT_OPT,JOB)
!//////// CREATE DATA FILES FOR TIME SERIES RECORDS OUTPUT ////////////
! OPEN DATA FILE TIME.DAT TO WRITE TIME SEQUENCE
OPEN (UNIT=999,FILE='time.dat',STATUS='UNKNOWN')
! TIME HISTORY RECORDS CAN BE CREATED FOR UP TO 9999 LOCATIONS
IF (OUT_OPT.EQ.1 .OR. OUT_OPT.EQ.2) THEN
NUM_REC = 0
TS_LOC = 0.0
TS_DAT = 0.0
OPEN (UNIT=25,FILE='ts_location.dat',STATUS='OLD',IOSTAT=STAT)
IF (STAT /=0) THEN
PRINT *,"ERROR:: CAN'T OPEN TS_LOCATION.DAT; EXITING."
STOP
END IF
WRITE (*,*) ' OUTPUT TIME HISTORY RECORDS AT:'
DO WHILE (.NOT. EOF(25))
NUM_REC = NUM_REC + 1
READ (25,*) TEMP1,TEMP2
ENDDO
REWIND(25)
ALLOCATE(TS_LOC(NUM_REC,2))
ALLOCATE(TS_ID(NUM_REC))
ALLOCATE(FNAME(NUM_REC))
TS_LOC = 0.0
TS_ID = 0
FNAME = ' '
DO K = 1,NUM_REC
READ (25,*) TS_LOC(K,1),TS_LOC(K,2)
IF (LO%LAYCORD.EQ.0 .AND. TS_LOC(K,1).LT.0.0) THEN
TS_LOC(K,1) = 360.0 + TS_LOC(K,1)
ENDIF
WRITE (*,*) ' ',TS_LOC(K,1),TS_LOC(K,2)
CALL GAUGE_LAYERID (ID,TS_LOC(K,1),TS_LOC(K,2),LO,LA)
TS_ID(K) = ID
ENDDO
CLOSE(25)
!OPEN DATA FILES TO WRITE TIME HISTORY DATA
DO NIC = 1,NUM_REC
WRITE (FNAME(NIC),1) NIC
1 FORMAT('ts_record',I4.4,'.dat')
IO_UNIT = 1000+NIC
OPEN (IO_UNIT,FILE=FNAME(NIC),STATUS='UNKNOWN')
ENDDO
END IF
!......................................................................
!////////////////////// SIMULATION BEGINS /////////////////////////////
WRITE (*,*) ' '
WRITE (*,*) '***************** OUTPUT RESULTS ******************'
WRITE (*,*) 'TIMESTEP MINUTE'
TIME = DBLE(ISTART-1)*LO%DT
! T_MINUTES = TIME/60.0
!.....OUTPUT TIME SEQUENCE AND TIME SERIES DATA AT T = 0.0
IF (ISTART.EQ.1) THEN
WRITE (999,*) TIME
CALL ALL_PRT (ISTART-1,LO,LA)
IF (OUT_OPT.EQ.1 .OR. OUT_OPT.EQ.2) THEN
DO NIC = 1,NUM_REC
CALL GET_TS (TS_DAT,TS_LOC(NIC,1),TS_LOC(NIC,2), &
LO,LA,TS_ID(NIC))
!* WRITE(*,*) TS_LOC(NIC,1),TS_LOC(NIC,2),TS_DAT
IO_UNIT = 1000+NIC
WRITE(IO_UNIT,'(F10.5)') TS_DAT
ENDDO
ENDIF
ENDIF
DO K = ISTART,KEND
TIME = TIME+LO%DT
T_MINUTES = TIME/60.0
IF (MOD(K,10) .EQ. 0) THEN
WRITE (*,'(I9,8X,F8.3)') K,T_MINUTES
ENDIF
!.... .. CALL HOT START FUNCTION
IF (K .EQ. START_STEP) THEN
CALL HOT_START (LO,LA,START_TYPE,START_STEP)
ENDIF
!* ! CREATE 5 BACKUPS DURING THE ENTIRE SIMULATION DURATION
!* IF (K.GT.1 .AND. MOD(K,FLOOR(KEND/5.)) .EQ. 0) THEN
!* CALL HOT_START (LO,LA,0,K)
!* ENDIF
!....... CALL MULTIPLE FAULT PLANE MODEL
IF (INI_SURF .EQ. 0 .OR. INI_SURF .EQ. 4) THEN
CALL GET_FLOOR_DEFORM (LO,LA,FAULT_INFO,TIME)
ENDIF
!... ...CALL LANDSLIDE MODEL
IF ( INI_SURF .EQ. 3 .OR. INI_SURF .EQ. 4 ) THEN
IF ( TIME.GE.LANDSLIDE_INFO%T(1) .AND. &
TIME.LE.LANDSLIDE_INFO%T(LANDSLIDE_INFO%NT) ) THEN
CALL LAND_SLIDE(LO,LANDSLIDE_INFO,TIME)
ENDIF
ENDIF
!////////SOLVE MASS CONSERVATION EQN FOR LAYER 1 (THE OUTEST LAYER)
!!....CALL SEDIMENT TRANSPORT MODEL
IF (LO%LAYCORD.NE.0 .AND. LO%SEDI_SWITCH .EQ. 0) THEN
CALL SED_TRANSPORT (LO)
ENDIF
CALL MASS (LO)
!.......SOLVE RADIATION OPEN BOUNDARY
CALL OPEN (LO)
!.......CALL WAVE MAKER TO GENERATE WAVES
IF (INI_SURF .EQ. 2) THEN
CALL WAVE_MAKER (TIME,LO,WAVE_INFO)
ENDIF
!.......APPLY FACTS INPUT CONDITION
IF (BC_TYPE.EQ.3) THEN
IF (TIME.GE.BCI_INFO%T(1) .AND. &
TIME.LT.BCI_INFO%T(BCI_INFO%NT)) THEN
CALL BC_INPUT (BCI_INFO,LO,TIME)
ENDIF
ENDIF
!//////////////////////////////////////////////////////////////////////
! CALL SUBLAYER COUPLED MODEL
!//////////////////////////////////////////////////////////////////////
IF (LO%NUM_CHILD .GE. 1) CALL ALL_GRID (LO,LA)
!//////////////////////////////////////////////////////////////////////
! SOLVE MOMENTUM CONSERVATION EQN FOR LAYER 1
!//////////////////////////////////////////////////////////////////////
CALL MOMENT (LO)
!.......USE SPONGE LAYER .....
IF (BC_TYPE.EQ.1) CALL SPONGE_LAYER(LO)
!*!.......WALL BOUNDARY
IF (BC_TYPE.EQ.2) CALL BC_WALL (LO,WAVE_INFO)
!......................................................................
! OUTPUT RESULTS AT GIVEN TIME INTERVAL (MINUTE)
!......................................................................
IF (MOD(K,IPRT) .EQ. 0) THEN
CALL ALL_PRT (K,LO,LA)
ENDIF
IF (OUT_OPT.EQ.0 .OR. OUT_OPT.EQ.2) THEN
CALL MAX_AMP (LO,LA,TIME,TEND)
ENDIF
!.......UPDATE VARIABLES OF LAYER 01 (LO) FOR NEXT TIME STEP
CALL CHANGE(LO)
!.......OUTPUT TIME SEQUENCE AND TIME HISTORY RECORDS AT T = K*LO%DT
WRITE (999,*) TIME
IF (OUT_OPT.EQ.1 .OR. OUT_OPT.EQ.2) THEN
DO NIC = 1,NUM_REC
CALL GET_TS (TS_DAT,TS_LOC(NIC,1),TS_LOC(NIC,2), &
LO,LA,TS_ID(NIC))
!* WRITE(*,*) TS_LOC(NIC,1),TS_LOC(NIC,2),TS_DAT
IO_UNIT = 1000+NIC
WRITE(IO_UNIT,'(F10.5)') TS_DAT
ENDDO
END IF
END DO
!..... CLOSE ALL OPENNED DATA FILES
CLOSE(999)
IF (OUT_OPT.EQ.1 .OR. OUT_OPT.EQ.2) THEN
DO NIC = 1,NUM_REC
CLOSE(1000+NIC)
ENDDO
END IF
CALL SYSTEM_CLOCK(count2, count_rate2, count_max2)
print *, "total time: ", (count2-count1)*1.0/count_rate2
END
!//////////////// END OF MAIN PROGRAM /////////////////////////////////
!//////////////////////////////////////////////////////////////////////
BLOCK DATA
COMMON /CONS/ ELMAX,GRAV,PI,R_EARTH,GX,EPS,ZERO,ONE,NUM_GRID, &
NUM_FLT,V_LIMIT,RAD_DEG,RAD_MIN
DATA ELMAX / -50.0 / ! UPPER LIMIT OF MAXIMUM RUNUP !
DATA GRAV / 9.807 / ! GRAVITY
DATA PI / 3.14159265358979 /
DATA R_EARTH / 6378000.0 / ! RADIUS OF THE EARTH
DATA GX / 1.0E-5 / ! FOR INUNDATION, DZ=0.0 IF DZ < GX
DATA EPS / 1.0E-10 / ! Z,U,V = 0.0 IF < EPS
DATA ZERO / 0.0 /
DATA ONE / 1.0 /
DATA NUM_GRID / 12 / ! MAX. NUMBER OF SUB-LEVEL GRIDS
DATA NUM_FLT / 999 / ! TOTAL NO. OF FAULT SEGMENTS
DATA V_LIMIT / 20.0 / ! UPPER LIMIT OF MAXIMUM VELOCITY !
DATA RAD_DEG / 0.01745329252 / ! ARC RADIAN OF 1 DEGREE
DATA RAD_MIN / 0.000290888208665721 / ! ARC RADIAN OF 1 MINUTE
END
!//////////////////////////////////////////////////////////////////////
!----------------------------------------------------------------------
SUBROUTINE CHANGE(LO)
!.....TRANSFER INFORMATION (FREE SURFACE ELEVATION, VOLUME FLUX IN
! X, Y DIRECTIONS)FROM LAST STEP TO NEXT STEP (FOR OUTEST LAYER)
!----------------------------------------------------------------------
USE LAYER_PARAMS
TYPE (LAYER) :: LO
DO J=1,LO%NY
DO I=1,LO%NX
LO%Z(I,J,1) = LO%Z(I,J,2)
LO%M(I,J,1) = LO%M(I,J,2)
LO%N(I,J,1) = LO%N(I,J,2)
IF (LO%LAYGOV.GT.1) THEN
LO%M0(I,J) = LO%M(I,J,1)
LO%N0(I,J) = LO%N(I,J,1)
ENDIF
END DO
END DO
!.....UPDATE BATHYMETRY IF TRANSIENT SEAFLOOR MOTION IS ENABLED
IF (LO%INI_SWITCH.EQ.3 .OR. LO%INI_SWITCH.EQ.4) THEN
LO%HT(:,:,1) = LO%HT(:,:,2)
LO%H(:,:) = LO%H(:,:) + LO%HT(:,:,2) - LO%HT(:,:,1)
ENDIF
!
RETURN
END
!----------------------------------------------------------------------
SUBROUTINE ALLOC (L,SWITCH)
! CREATED BY TOM LOGAN, UAF (2005)
! UPDATED ON SEP 17 2006 (XIAOMING WANG,CORNELL UNIVERSITY)
! SWITCH:
! 1 -- TOP-LEVEL GRID (OUTEST LAYER)
! OTHER -- SUB-LEVEL GRID LAYERS
! ...
! LINCHK:
! 0 -- NONLINEAR TERMS WON'T BE CALCULATED IN NSWE
! 1 -- NONLINEAR TERMS WILL BE CALCULATED IN NSWE
! MODSCM:
! 0 -- MODIFIED SCHEME IS IMPLETMENTED (IMAMURA,1987;CHO,1995)
! 1 -- MODIFIED SCHEME IS NOT IMPLETMENTED
!.....UPDATED ON NOV 03 2008 (XIAOMING WANG, GNS)
!----------------------------------------------------------------------
USE LAYER_PARAMS
TYPE(LAYER) :: L
INTEGER :: SWITCH ! 1 - TOP LAYER; ELSE - OTHER LAYER
COMMON /CONS/ ELMAX,GRAV,PI,R_EARTH,GX,EPS,ZERO,ONE,NUM_GRID, &
NUM_FLT,V_LIMIT,RAD_DEG,RAD_MIN
L%RX = L%DT/L%DX
L%RY = L%DT/L%DY
L%GRX = GRAV*L%RX
L%GRY = GRAV*L%RY
IF (SWITCH .EQ. 1) THEN
L%LINCHK = 1
L%MODSCM = 0
L%CORNERS(1) = 1
L%CORNERS(2) = L%NX
L%CORNERS(3) = 1
L%CORNERS(4) = L%NY
ELSE
L%LINCHK = 1
L%MODSCM = 1
END IF
IF ( L%LAYSWITCH.EQ.0 ) THEN
ALLOCATE(L%Z(L%NX,L%NY,2))
ALLOCATE(L%M(L%NX,L%NY,2))
ALLOCATE(L%N(L%NX,L%NY,2))
ALLOCATE(L%H(L%NX,L%NY))
ALLOCATE(L%XFLUX(L%NX,2))
ALLOCATE(L%YFLUX(L%NY,2))
ALLOCATE(L%HP(L%NX,L%NY))
ALLOCATE(L%HQ(L%NX,L%NY))
ALLOCATE(L%DZ(L%NX,L%NY,2))
ALLOCATE(L%DEFORM(L%NX,L%NY))
IF (L%LAYCORD.EQ.0) THEN
ALLOCATE(L%R0(L%NX,L%NY))
ALLOCATE(L%R1(L%NX,L%NY))
ALLOCATE(L%R11(L%NX,L%NY))
ALLOCATE(L%R2(L%NX,L%NY))
ALLOCATE(L%R21(L%NX,L%NY))
ALLOCATE(L%R22(L%NX,L%NY))
ALLOCATE(L%R3(L%NX,L%NY))
ALLOCATE(L%R4(L%NX,L%NY))
ALLOCATE(L%R5(L%NX,L%NY))
ALLOCATE(L%R6(L%NY))
ENDIF
IF (L%INI_SWITCH.EQ.3 .OR. L%INI_SWITCH.EQ.4) THEN
ALLOCATE(L%HT(L%NX,L%NY,2))
L%HT = 0.0
ENDIF
IF (L%LAYGOV.GE.1) THEN
ALLOCATE(L%MASK(L%NX,L%NY))
L%MASK = 0
ENDIF
IF (L%LAYGOV.GT.1) THEN
ALLOCATE(L%M0(L%NX,L%NY))
ALLOCATE(L%N0(L%NX,L%NY))
! ALLOCATE(L%MASK(L%NX,L%NY))
L%M0 = 0.0
L%N0 = 0.0
! L%MASK = 0
ENDIF
IF (L%FRIC_SWITCH.EQ.2) THEN
ALLOCATE(L%FRIC_VCOEF(L%NX,L%NY))
L%FRIC_VCOEF = L%FRIC_COEF
ENDIF
!* IF (L%LAYGOV.EQ.1 .OR. L%LAYGOV.EQ.9) L%MODSCM = 1
IF (MOD(L%LAYGOV,2).NE.0) L%MODSCM = 1
ALLOCATE(L%Z_MAX(L%NX,L%NY))
ALLOCATE(L%Z_MIN(L%NX,L%NY))
!* IF (L%UPZ .EQ. .FALSE.) THEN
!* ALLOCATE(L%POS(L%NX,L%NY,2))
!* ALLOCATE(L%CXY(L%NX,L%NY,4))
!* L%POS = 0
!* L%CXY = 0.0
!* ENDIF
L%Z = 0.0
L%M = 0.0
L%N = 0.0
L%H = 0.0
L%HP = 0.0
L%HQ = 0.0
L%DZ = 0.0
L%DEFORM = 0.0
L%XFLUX = 0.0
L%YFLUX = 0.0
! L%NAME = ''
!.......UPZ - FLAG TO IDENTIFY COORDINATE SYSTEM
! .TRUE. - USE THE SAME COORD. AS ITS PARENT GRID
! .FALSE. - USE DIFFERENT COORD.
! L%UPZ = .TRUE.
L%Z_MAX = 0.0
L%Z_MIN = 0.0
!.......L%X,L%Y,L%DX,L%DY,L%DT ARE ALLOCATED AND
! INITIALIZED IN DX_CALC AND SUBGRID_MATCHING
END IF
IF (L%LAYSWITCH.EQ.0 .AND. L%LAYGOV.GE.2) THEN
ALLOCATE(L%ALPHA(L%NX,L%NY))
ALLOCATE(L%A1X(L%NX,L%NY))
ALLOCATE(L%A2X(L%NX,L%NY))
ALLOCATE(L%A1Y(L%NX,L%NY))
ALLOCATE(L%A2Y(L%NX,L%NY))
ALLOCATE(L%CF(L%NX,L%NY,4))
ALLOCATE(L%CB(L%NX,L%NY,4))
L%ALPHA = 1.0
L%CF = 0.0
L%CB = 0.0
END IF
IF (SWITCH .EQ. 1 .AND. L%BC_TYPE.EQ.1) THEN
ALLOCATE(L%SPONGE_COEFX(L%NX,L%NY))
ALLOCATE(L%SPONGE_COEFY(L%NX,L%NY))
L%SPONGE_COEFX = 0.0
L%SPONGE_COEFY = 0.0
ENDIF
IF (L%FRIC_SWITCH.EQ.3) THEN
L%SEDI_SWITCH = 0
ALLOCATE(L%DH(L%NX,L%NY,2))
L%DH = 0.0
ELSE
L%SEDI_SWITCH = 1
ENDIF
RETURN
END
!----------------------------------------------------------------------
SUBROUTINE SPHERICAL_PARAMETERS(LO,LA)
! CREATED BY TOM LOGAN, UAF (2004)
! ******** PARAMETERS FOR SPHERICAL COORD. *******
! Updated on Dec21 2008 (Xiaoming Wang, GNS)
!----------------------------------------------------------------------
USE LAYER_PARAMS
TYPE (LAYER) :: LO
TYPE (LAYER), DIMENSION(NUM_GRID) :: LA
COMMON /CONS/ ELMAX,GRAV,PI,R_EARTH,GX,EPS,ZERO,ONE,NUM_GRID, &
NUM_FLT,V_LIMIT,RAD_DEG,RAD_MIN
IF (LO%LAYSWITCH .EQ. 0 .AND. LO%LAYCORD .EQ. 0) THEN
WRITE (*,*) ' CALCULATING COEFFICIENTS FOR SPHERICAL COORD.'
CALL PARA(LO)
ENDIF
DO I=1,NUM_GRID
IF (LA(I)%LAYSWITCH .EQ. 0 .AND. LA(I)%LAYCORD .EQ. 0) THEN
CALL PARA(LA(I))
ENDIF
END DO
RETURN
END
!----------------------------------------------------------------------
SUBROUTINE PARA(LO)
!......................................................................
!DESCRIPTION:
! #. CALCULATE PARAMETERS FOR SPHERICAL COORD.
!NOTE:
! #. REVISED SIGNIFICANTLY BY TOM LOGAN, UAF (2005)
! #. UPDATED ON NOV.07, 2008 (XIAOMING WANG, GNS)
! LINES STARTING WITH "!!" ARE REMOVED
! LINES ENDED WITH "!!" ARE ADDED
! #. UPDATED ON FEB06 2009 (XIAOMING WANG, GNS)
!----------------------------------------------------------------------
USE LAYER_PARAMS
TYPE (LAYER) :: LO
COMMON /CONS/ ELMAX,GRAV,PI,R_EARTH,GX,EPS,ZERO,ONE,NUM_GRID, &
NUM_FLT,V_LIMIT,RAD_DEG,RAD_MIN
DATA CORI_W/7.2722E-5/, OSIXTY/0.016666666667/
DATA TWLVTH/0.08333333333333/
!.....PARAMETERS FOR SPHERICAL COORD IN GRID LAYER LO
IF (LO%LAYCORD .EQ. 0) THEN
IF (LO%LAYGOV.LE.1 .OR. LO%LAYGOV.EQ.9) THEN
!* DX = R_EARTH*LO%DX*RAD_MIN
!* RR = LO%DT/DX
RR = LO%DT/R_EARTH
!* RS = 0.5*G*RR
RS = GRAV*RR !!
RT = 0.5*LO%DT*CORI_W
DO J=1,LO%NY
ANG_M = LO%Y(J)*RAD_DEG
ANG_N = (LO%Y(J)+0.5*LO%DEL_Y(J)*OSIXTY)*RAD_DEG
COS_M = COS(ANG_M)
COS_N = COS(ANG_N)
SIN_M = SIN(ANG_M)
SIN_N = SIN(ANG_N)
DO I=1,LO%NX
LO%R0(I,J) = RR/(LO%DEL_Y(J)*RAD_MIN) !!
LO%R1(I,J) = RR/COS_M/(LO%DX*RAD_MIN)
LO%R11(I,J) = RR/COS_M/(LO%DEL_Y(J)*RAD_MIN)
IF (LO%LAYGOV.EQ.0) THEN
LO%R2(I,J) = RS/COS_M/(LO%DX*RAD_MIN) &
*AMAX1(LO%HP(I,J),ZERO)
ENDIF
IF (LO%LAYGOV.EQ.1) THEN
LO%R2(I,J) = RS/COS_M/(LO%DX*RAD_MIN)
ENDIF
LO%R21(I,J) = RR/COS_M/(LO%DX*RAD_MIN)
LO%R22(I,J) = RR/COS_N/(LO%DX*RAD_MIN)
LO%R3(I,J) = RT*SIN_M
IF (LO%LAYGOV.EQ.0) THEN
LO%R4(I,J) = RS/(LO%DEL_Y(J)*RAD_MIN) &
*AMAX1(LO%HQ(I,J),ZERO)
ENDIF
IF (LO%LAYGOV.EQ.1) THEN
LO%R4(I,J) = RS/(LO%DEL_Y(J)*RAD_MIN)
ENDIF
LO%R5(I,J) = RT*SIN_N
END DO
LO%R6(J) = COS_N
END DO
ELSEIF (LO%LAYGOV.EQ.2 .OR. LO%LAYGOV.EQ.3) THEN
!* DX = R_EARTH*LO%DX*RAD_MIN
!* RR = LO%DT/DX
RR = LO%DT/R_EARTH
!* RS = 0.5*G*RR
RS = GRAV*RR !!
RT = 0.5*LO%DT*CORI_W
DO J=1,LO%NY
ANG_M = LO%Y(J)*RAD_DEG
ANG_N = (LO%Y(J)+0.5*LO%DEL_Y(J)*OSIXTY)*RAD_DEG
COS_M = COS(ANG_M)
COS_N = COS(ANG_N)
SIN_M = SIN(ANG_M)
SIN_N = SIN(ANG_N)
DO I=1,LO%NX
IF (LO%MASK(I,J).EQ.0) THEN
DX = LO%DX
DY = LO%DEL_Y(J)
ENDIF
IF (LO%MASK(I,J).EQ.1) THEN
DX = LO%DX*LO%ALPHA(I,J)
DY = LO%DEL_Y(J)*LO%ALPHA(I,J)
ENDIF
LO%R0(I,J) = RR/(DY*RAD_MIN) !!
LO%R1(I,J) = RR/COS_M/(DX*RAD_MIN)
LO%R11(I,J) = RR/COS_M/(DY*RAD_MIN)
IF (LO%LAYGOV.EQ.2) THEN
LO%R2(I,J) = RS/COS_M/(DX*RAD_MIN) &
*AMAX1(LO%HP(I,J),ZERO)
ENDIF
IF (LO%LAYGOV.EQ.3) THEN
LO%R2(I,J) = RS/COS_M/(DX*RAD_MIN)
ENDIF
LO%R21(I,J) = RR/COS_M/(DX*RAD_MIN)
LO%R22(I,J) = RR/COS_N/(DX*RAD_MIN)
LO%R3(I,J) = RT*SIN_M
IF (LO%LAYGOV.EQ.2) THEN
LO%R4(I,J) = RS/(DY*RAD_MIN) &
*AMAX1(LO%HQ(I,J),ZERO)
ENDIF
IF (LO%LAYGOV.EQ.3) THEN
LO%R4(I,J) = RS/(DY*RAD_MIN)
ENDIF
LO%R5(I,J) = RT*SIN_N
END DO
LO%R6(J) = COS_N
END DO
ENDIF
ENDIF
RETURN
END
!----------------------------------------------------------------------
SUBROUTINE SED_TRANSPORT (LO)
!.....TENTATIVE TRY TO INCLUDING SEDIMENT TRANSPORT MODEL (APR. 2007)
! SEDIMENT FLUX IS COMPUTED ACCORDING TO RIBBERINK (1998)
!----------------------------------------------------------------------
USE LAYER_PARAMS
TYPE (LAYER) :: LO
REAL D, DS, THETA_C, THETA_P,THETA_Q,G,RHO,RHO_S
REAL H(LO%NX,LO%NY),DH(LO%NX,LO%NY)
REAL QX(LO%NX,LO%NY),QY(LO%NX,LO%NY)
CHARACTER(LEN=30) FNAME
G = 9.807
RHO = 1000.0
RHO_S = 2650.0
DS = 0.0002
POROSITY = 0.3
THETA_C = 0.04
S = RHO_S/RHO
! Q(:,:) = 0.0
QX(:,:) = 0.0
QY(:,:) = 0.0
! H(:,:) = -LO%H(:,:)
DH(:,:) = LO%H(:,:) + LO%Z(:,:,1)
DO I = 1,LO%NX
IM1 = I-1
IP1 = I+1
IF (IM1.LE.0) IM1=1
IF (IP1.GE.LO%NX) IP1=LO%NX
DO J = 1,LO%NY
JM1 = J-1
JP1 = J+1
IF (JM1.LE.1) JM1=1
IF (JP1.GE.LO%NY) JP1=LO%NY
!TOTAL WATER DEPTH AT DISCHARGE LOCATION P
HP = 0.5*(DH(I,J)+DH(IP1,J))
!TOTAL WATER DEPTH AT DISCHARGE LOCATION Q
HQ = 0.5*(DH(I,J)+DH(I,JP1))
!SEDIMENT FLUX AT DISCHARGE LOCATION P
IF (HP .GT. 0.05) THEN
COEF_P = RHO*G*(LO%FRIC_COEF**2)/(HP**2.33333)
!Y COMPONENT OF VOLUME FLUX AT DISCHARGE LOCATION P
XQQ = 0.25*(LO%N(I,J,1)+LO%N(IP1,J,1) &
+ LO%N(I,JM1,1)+LO%N(IP1,JM1,1))
!TOTAL VOLUME FLUX AT DISCHARGE LOCATION P
FLUX_P = SQRT(LO%M(I,J,1)**2+XQQ**2)
!X COMPONENT OF BOTTOM SHEAR STRESS AT DISCHARGE LOCATION P
TPX = COEF_P*LO%M(I,J,1)*FLUX_P
!Y COMPONENT OF BOTTOM SHEAR STRESS AT DISCHARGE LOCATION P
TPY = COEF_P*XQQ*FLUX_P
! TOTAL SHEAR STRESS AT DISCHARGE LOCATION P
TPB = SQRT(TPX**2+TPY**2)
!SHIELDS PARAMETER AT DISCHARGE LOCATION P
THETA_P = TPB/((RHO_S-RHO)*G*DS)
IF (THETA_P .GT. THETA_C) THEN
!RIBBERINK (1998)
QS_P = 11.0*SQRT((S-1.0)*G*DS**3) &
* ((THETA_P-THETA_C)**1.65)
IF (QS_P.GT.0.02) QS_P = 0.02
!X COMPONENT OF SAND FLUX AT DISCHARGE LOCATION P
IF (FLUX_P .GT. 1.0E-10) &
QX(I,J) = QS_P * LO%M(I,J,1)/FLUX_P
ENDIF
ENDIF
!SEDIMENT FLUX AT DISCHARGE LOCATION Q
IF (HQ .GT. 0.05) THEN
COEF_Q = RHO*G*(LO%FRIC_COEF**2)/(HQ**2.33333)
XPP = 0.25*(LO%M(I,J,1)+LO%M(I,JP1,1)+LO%M(IM1,J,1) &
+ LO%M(IM1,JP1,1))
FLUX_Q = SQRT(XPP**2+LO%N(I,J,1)**2)
TQX = COEF_Q*XPP*FLUX_Q
TQY = COEF_Q*LO%N(I,J,1)*FLUX_Q
TQB = SQRT(TQX**2+TQY**2)
THETA_Q = TQB/((RHO_S-RHO)*G*DS)
IF (THETA_Q .GT. THETA_C) THEN
!RIBBERINK (1998)
QS_Q = 11.0*SQRT((S-1.0)*G*DS**3) &
*((THETA_Q-THETA_C)**1.65)
IF (QS_Q.GT.0.02) QS_Q = 0.02
IF (FLUX_Q .GT. 1.0E-10) &
QY(I,J) = QS_Q*LO%N(I,J,1)/FLUX_Q
ENDIF
ENDIF
ENDDO
ENDDO
DO I = 2,LO%NX-1
DO J = 2,LO%NY-1
LO%DH(I,J,2) = LO%DH(I,J,1)-LO%RX*(1/(1-POROSITY)) &
*(QX(I,J)-QX(I-1,J)+QY(I,J)-QY(I,J-1))
ENDDO
ENDDO
LO%H(:,:) = LO%H(:,:)-LO%DH(:,:,2)
RETURN
END
write (*,*) ' '
write (*,*) '***************** COMCOT ****************'
write (*,*) '* *'
write (*,*) '* VERSION=',Version,' *'
write (*,*) '* *'
write (*,*) '********************************************'
!----------------------------------------------------------------------
!///// READ USER-DEFINED OPTIONS (NO. OF LAYERS, COORDINATE, GOVERN EQN.)
CALL SYSTEM_CLOCK(count1, count_rate1, count_max1)
CALL READ_CONFIG (LO,LA,TEND,T_INV,INI_SURF,WAVE_INFO, &
FAULT_INFO,LANDSLIDE_INFO,BCI_INFO, &
START_TYPE,START_TIME,BC_TYPE,OUT_OPT,JOB)
!///// READING BATH/TOPO DATA /////////////////////////////////////////
CALL READ_BATHYMETRY(LO,LA)
!///// READ INITIAL CONDITION /////////////////////////////////////////
IF (BC_TYPE.NE.3) THEN
CALL GET_INI_SURF(LO,LA,INI_SURF,WAVE_INFO,FAULT_INFO, &
LANDSLIDE_INFO)
ENDIF
!///// ADJUST BATHYMETRY, TIME STEP AND CALC DEPENDENT PARAMETERS /////
CALL ADJUST_BATHYMETRY (LO,LA)
CALL CR_CHECK (LO,LA)
CALL ALPHA_CALC (LO,LA)
!CALCULATE PARAMETERS USED FOR SPHERICAL COORD.
CALL SPHERICAL_PARAMETERS(LO,LA)
!///// SET UP BOUNDARY CONDITIONS /////////////////////////////////////
IF (BC_TYPE.EQ.1) CALL SPONGE_COEF (LO)
IF (BC_TYPE.EQ.2) CALL BC_WALL (LO,WAVE_INFO)
IF (BC_TYPE.EQ.3) CALL GET_BC_DATA (BCI_INFO,LO)
!/////DETERMINE STARTING TIME STEP # //////////////////////////////////
! WRITE(*,*) START_STEP
START_STEP = NINT(START_TIME/LO%DT)
ISTART = START_STEP
IF (START_TYPE .EQ. 0) ISTART = 1
!//////// DISPLAY INPUT PARAMETERS ////////////////////////////////////
CALL SCREEN_DISP (LO,LA,FAULT_INFO,LANDSLIDE_INFO,WAVE_INFO, &
INI_SURF,TEND,T_INV,KEND,IPRT,START_TYPE, &
START_TIME,BC_TYPE,BCI_INFO,OUT_OPT,JOB)
!//////// CREATE DATA FILES FOR TIME SERIES RECORDS OUTPUT ////////////
! OPEN DATA FILE TIME.DAT TO WRITE TIME SEQUENCE
OPEN (UNIT=999,FILE='time.dat',STATUS='UNKNOWN')
! TIME HISTORY RECORDS CAN BE CREATED FOR UP TO 9999 LOCATIONS
IF (OUT_OPT.EQ.1 .OR. OUT_OPT.EQ.2) THEN
NUM_REC = 0
TS_LOC = 0.0
TS_DAT = 0.0
OPEN (UNIT=25,FILE='ts_location.dat',STATUS='OLD',IOSTAT=STAT)
IF (STAT /=0) THEN
PRINT *,"ERROR:: CAN'T OPEN TS_LOCATION.DAT; EXITING."
STOP
END IF
WRITE (*,*) ' OUTPUT TIME HISTORY RECORDS AT:'
DO WHILE (.NOT. EOF(25))
NUM_REC = NUM_REC + 1
READ (25,*) TEMP1,TEMP2
ENDDO
REWIND(25)
ALLOCATE(TS_LOC(NUM_REC,2))
ALLOCATE(TS_ID(NUM_REC))
ALLOCATE(FNAME(NUM_REC))
TS_LOC = 0.0
TS_ID = 0
FNAME = ' '
DO K = 1,NUM_REC
READ (25,*) TS_LOC(K,1),TS_LOC(K,2)
IF (LO%LAYCORD.EQ.0 .AND. TS_LOC(K,1).LT.0.0) THEN
TS_LOC(K,1) = 360.0 + TS_LOC(K,1)
ENDIF
WRITE (*,*) ' ',TS_LOC(K,1),TS_LOC(K,2)
CALL GAUGE_LAYERID (ID,TS_LOC(K,1),TS_LOC(K,2),LO,LA)
TS_ID(K) = ID
ENDDO
CLOSE(25)
!OPEN DATA FILES TO WRITE TIME HISTORY DATA
DO NIC = 1,NUM_REC
WRITE (FNAME(NIC),1) NIC
1 FORMAT('ts_record',I4.4,'.dat')
IO_UNIT = 1000+NIC
OPEN (IO_UNIT,FILE=FNAME(NIC),STATUS='UNKNOWN')
ENDDO
END IF
!......................................................................
!////////////////////// SIMULATION BEGINS /////////////////////////////
WRITE (*,*) ' '
WRITE (*,*) '***************** OUTPUT RESULTS ******************'
WRITE (*,*) 'TIMESTEP MINUTE'
TIME = DBLE(ISTART-1)*LO%DT
! T_MINUTES = TIME/60.0
!.....OUTPUT TIME SEQUENCE AND TIME SERIES DATA AT T = 0.0
IF (ISTART.EQ.1) THEN
WRITE (999,*) TIME
CALL ALL_PRT (ISTART-1,LO,LA)
IF (OUT_OPT.EQ.1 .OR. OUT_OPT.EQ.2) THEN
DO NIC = 1,NUM_REC
CALL GET_TS (TS_DAT,TS_LOC(NIC,1),TS_LOC(NIC,2), &
LO,LA,TS_ID(NIC))
!* WRITE(*,*) TS_LOC(NIC,1),TS_LOC(NIC,2),TS_DAT
IO_UNIT = 1000+NIC
WRITE(IO_UNIT,'(F10.5)') TS_DAT
ENDDO
ENDIF
ENDIF
DO K = ISTART,KEND
TIME = TIME+LO%DT
T_MINUTES = TIME/60.0
IF (MOD(K,10) .EQ. 0) THEN
WRITE (*,'(I9,8X,F8.3)') K,T_MINUTES
ENDIF
!.... .. CALL HOT START FUNCTION
IF (K .EQ. START_STEP) THEN
CALL HOT_START (LO,LA,START_TYPE,START_STEP)
ENDIF
!* ! CREATE 5 BACKUPS DURING THE ENTIRE SIMULATION DURATION
!* IF (K.GT.1 .AND. MOD(K,FLOOR(KEND/5.)) .EQ. 0) THEN
!* CALL HOT_START (LO,LA,0,K)
!* ENDIF
!....... CALL MULTIPLE FAULT PLANE MODEL
IF (INI_SURF .EQ. 0 .OR. INI_SURF .EQ. 4) THEN
CALL GET_FLOOR_DEFORM (LO,LA,FAULT_INFO,TIME)
ENDIF
!... ...CALL LANDSLIDE MODEL
IF ( INI_SURF .EQ. 3 .OR. INI_SURF .EQ. 4 ) THEN
IF ( TIME.GE.LANDSLIDE_INFO%T(1) .AND. &
TIME.LE.LANDSLIDE_INFO%T(LANDSLIDE_INFO%NT) ) THEN
CALL LAND_SLIDE(LO,LANDSLIDE_INFO,TIME)
ENDIF
ENDIF
!////////SOLVE MASS CONSERVATION EQN FOR LAYER 1 (THE OUTEST LAYER)
!!....CALL SEDIMENT TRANSPORT MODEL
IF (LO%LAYCORD.NE.0 .AND. LO%SEDI_SWITCH .EQ. 0) THEN
CALL SED_TRANSPORT (LO)
ENDIF
CALL MASS (LO)
!.......SOLVE RADIATION OPEN BOUNDARY
CALL OPEN (LO)
!.......CALL WAVE MAKER TO GENERATE WAVES
IF (INI_SURF .EQ. 2) THEN
CALL WAVE_MAKER (TIME,LO,WAVE_INFO)
ENDIF
!.......APPLY FACTS INPUT CONDITION
IF (BC_TYPE.EQ.3) THEN
IF (TIME.GE.BCI_INFO%T(1) .AND. &
TIME.LT.BCI_INFO%T(BCI_INFO%NT)) THEN
CALL BC_INPUT (BCI_INFO,LO,TIME)
ENDIF
ENDIF
!//////////////////////////////////////////////////////////////////////
! CALL SUBLAYER COUPLED MODEL
!//////////////////////////////////////////////////////////////////////
IF (LO%NUM_CHILD .GE. 1) CALL ALL_GRID (LO,LA)
!//////////////////////////////////////////////////////////////////////
! SOLVE MOMENTUM CONSERVATION EQN FOR LAYER 1
!//////////////////////////////////////////////////////////////////////
CALL MOMENT (LO)
!.......USE SPONGE LAYER .....
IF (BC_TYPE.EQ.1) CALL SPONGE_LAYER(LO)
!*!.......WALL BOUNDARY
IF (BC_TYPE.EQ.2) CALL BC_WALL (LO,WAVE_INFO)
!......................................................................
! OUTPUT RESULTS AT GIVEN TIME INTERVAL (MINUTE)
!......................................................................
IF (MOD(K,IPRT) .EQ. 0) THEN
CALL ALL_PRT (K,LO,LA)
ENDIF
IF (OUT_OPT.EQ.0 .OR. OUT_OPT.EQ.2) THEN
CALL MAX_AMP (LO,LA,TIME,TEND)
ENDIF
!.......UPDATE VARIABLES OF LAYER 01 (LO) FOR NEXT TIME STEP
CALL CHANGE(LO)
!.......OUTPUT TIME SEQUENCE AND TIME HISTORY RECORDS AT T = K*LO%DT
WRITE (999,*) TIME
IF (OUT_OPT.EQ.1 .OR. OUT_OPT.EQ.2) THEN
DO NIC = 1,NUM_REC
CALL GET_TS (TS_DAT,TS_LOC(NIC,1),TS_LOC(NIC,2), &
LO,LA,TS_ID(NIC))
!* WRITE(*,*) TS_LOC(NIC,1),TS_LOC(NIC,2),TS_DAT
IO_UNIT = 1000+NIC
WRITE(IO_UNIT,'(F10.5)') TS_DAT
ENDDO
END IF
END DO
!..... CLOSE ALL OPENNED DATA FILES
CLOSE(999)
IF (OUT_OPT.EQ.1 .OR. OUT_OPT.EQ.2) THEN
DO NIC = 1,NUM_REC
CLOSE(1000+NIC)
ENDDO
END IF
CALL SYSTEM_CLOCK(count2, count_rate2, count_max2)
print *, "total time: ", (count2-count1)*1.0/count_rate2
END
!//////////////// END OF MAIN PROGRAM /////////////////////////////////
!//////////////////////////////////////////////////////////////////////
BLOCK DATA
COMMON /CONS/ ELMAX,GRAV,PI,R_EARTH,GX,EPS,ZERO,ONE,NUM_GRID, &
NUM_FLT,V_LIMIT,RAD_DEG,RAD_MIN
DATA ELMAX / -50.0 / ! UPPER LIMIT OF MAXIMUM RUNUP !
DATA GRAV / 9.807 / ! GRAVITY
DATA PI / 3.14159265358979 /
DATA R_EARTH / 6378000.0 / ! RADIUS OF THE EARTH
DATA GX / 1.0E-5 / ! FOR INUNDATION, DZ=0.0 IF DZ < GX
DATA EPS / 1.0E-10 / ! Z,U,V = 0.0 IF < EPS
DATA ZERO / 0.0 /
DATA ONE / 1.0 /
DATA NUM_GRID / 12 / ! MAX. NUMBER OF SUB-LEVEL GRIDS
DATA NUM_FLT / 999 / ! TOTAL NO. OF FAULT SEGMENTS
DATA V_LIMIT / 20.0 / ! UPPER LIMIT OF MAXIMUM VELOCITY !
DATA RAD_DEG / 0.01745329252 / ! ARC RADIAN OF 1 DEGREE
DATA RAD_MIN / 0.000290888208665721 / ! ARC RADIAN OF 1 MINUTE
END
!//////////////////////////////////////////////////////////////////////
!----------------------------------------------------------------------
SUBROUTINE CHANGE(LO)
!.....TRANSFER INFORMATION (FREE SURFACE ELEVATION, VOLUME FLUX IN
! X, Y DIRECTIONS)FROM LAST STEP TO NEXT STEP (FOR OUTEST LAYER)
!----------------------------------------------------------------------
USE LAYER_PARAMS
TYPE (LAYER) :: LO
DO J=1,LO%NY
DO I=1,LO%NX
LO%Z(I,J,1) = LO%Z(I,J,2)
LO%M(I,J,1) = LO%M(I,J,2)
LO%N(I,J,1) = LO%N(I,J,2)
IF (LO%LAYGOV.GT.1) THEN
LO%M0(I,J) = LO%M(I,J,1)
LO%N0(I,J) = LO%N(I,J,1)
ENDIF
END DO
END DO
!.....UPDATE BATHYMETRY IF TRANSIENT SEAFLOOR MOTION IS ENABLED
IF (LO%INI_SWITCH.EQ.3 .OR. LO%INI_SWITCH.EQ.4) THEN
LO%HT(:,:,1) = LO%HT(:,:,2)
LO%H(:,:) = LO%H(:,:) + LO%HT(:,:,2) - LO%HT(:,:,1)
ENDIF
!
RETURN
END
!----------------------------------------------------------------------
SUBROUTINE ALLOC (L,SWITCH)
! CREATED BY TOM LOGAN, UAF (2005)
! UPDATED ON SEP 17 2006 (XIAOMING WANG,CORNELL UNIVERSITY)
! SWITCH:
! 1 -- TOP-LEVEL GRID (OUTEST LAYER)
! OTHER -- SUB-LEVEL GRID LAYERS
! ...
! LINCHK:
! 0 -- NONLINEAR TERMS WON'T BE CALCULATED IN NSWE
! 1 -- NONLINEAR TERMS WILL BE CALCULATED IN NSWE
! MODSCM:
! 0 -- MODIFIED SCHEME IS IMPLETMENTED (IMAMURA,1987;CHO,1995)
! 1 -- MODIFIED SCHEME IS NOT IMPLETMENTED
!.....UPDATED ON NOV 03 2008 (XIAOMING WANG, GNS)
!----------------------------------------------------------------------
USE LAYER_PARAMS
TYPE(LAYER) :: L
INTEGER :: SWITCH ! 1 - TOP LAYER; ELSE - OTHER LAYER
COMMON /CONS/ ELMAX,GRAV,PI,R_EARTH,GX,EPS,ZERO,ONE,NUM_GRID, &
NUM_FLT,V_LIMIT,RAD_DEG,RAD_MIN
L%RX = L%DT/L%DX
L%RY = L%DT/L%DY
L%GRX = GRAV*L%RX
L%GRY = GRAV*L%RY
IF (SWITCH .EQ. 1) THEN
L%LINCHK = 1
L%MODSCM = 0
L%CORNERS(1) = 1
L%CORNERS(2) = L%NX
L%CORNERS(3) = 1
L%CORNERS(4) = L%NY
ELSE
L%LINCHK = 1
L%MODSCM = 1
END IF
IF ( L%LAYSWITCH.EQ.0 ) THEN
ALLOCATE(L%Z(L%NX,L%NY,2))
ALLOCATE(L%M(L%NX,L%NY,2))
ALLOCATE(L%N(L%NX,L%NY,2))
ALLOCATE(L%H(L%NX,L%NY))
ALLOCATE(L%XFLUX(L%NX,2))
ALLOCATE(L%YFLUX(L%NY,2))
ALLOCATE(L%HP(L%NX,L%NY))
ALLOCATE(L%HQ(L%NX,L%NY))
ALLOCATE(L%DZ(L%NX,L%NY,2))
ALLOCATE(L%DEFORM(L%NX,L%NY))
IF (L%LAYCORD.EQ.0) THEN
ALLOCATE(L%R0(L%NX,L%NY))
ALLOCATE(L%R1(L%NX,L%NY))
ALLOCATE(L%R11(L%NX,L%NY))
ALLOCATE(L%R2(L%NX,L%NY))
ALLOCATE(L%R21(L%NX,L%NY))
ALLOCATE(L%R22(L%NX,L%NY))
ALLOCATE(L%R3(L%NX,L%NY))
ALLOCATE(L%R4(L%NX,L%NY))
ALLOCATE(L%R5(L%NX,L%NY))
ALLOCATE(L%R6(L%NY))
ENDIF
IF (L%INI_SWITCH.EQ.3 .OR. L%INI_SWITCH.EQ.4) THEN
ALLOCATE(L%HT(L%NX,L%NY,2))
L%HT = 0.0
ENDIF
IF (L%LAYGOV.GE.1) THEN
ALLOCATE(L%MASK(L%NX,L%NY))
L%MASK = 0
ENDIF
IF (L%LAYGOV.GT.1) THEN
ALLOCATE(L%M0(L%NX,L%NY))
ALLOCATE(L%N0(L%NX,L%NY))
! ALLOCATE(L%MASK(L%NX,L%NY))
L%M0 = 0.0
L%N0 = 0.0
! L%MASK = 0
ENDIF
IF (L%FRIC_SWITCH.EQ.2) THEN
ALLOCATE(L%FRIC_VCOEF(L%NX,L%NY))
L%FRIC_VCOEF = L%FRIC_COEF
ENDIF
!* IF (L%LAYGOV.EQ.1 .OR. L%LAYGOV.EQ.9) L%MODSCM = 1
IF (MOD(L%LAYGOV,2).NE.0) L%MODSCM = 1
ALLOCATE(L%Z_MAX(L%NX,L%NY))
ALLOCATE(L%Z_MIN(L%NX,L%NY))
!* IF (L%UPZ .EQ. .FALSE.) THEN
!* ALLOCATE(L%POS(L%NX,L%NY,2))
!* ALLOCATE(L%CXY(L%NX,L%NY,4))
!* L%POS = 0
!* L%CXY = 0.0
!* ENDIF
L%Z = 0.0
L%M = 0.0
L%N = 0.0
L%H = 0.0
L%HP = 0.0
L%HQ = 0.0
L%DZ = 0.0
L%DEFORM = 0.0
L%XFLUX = 0.0
L%YFLUX = 0.0
! L%NAME = ''
!.......UPZ - FLAG TO IDENTIFY COORDINATE SYSTEM
! .TRUE. - USE THE SAME COORD. AS ITS PARENT GRID
! .FALSE. - USE DIFFERENT COORD.
! L%UPZ = .TRUE.
L%Z_MAX = 0.0
L%Z_MIN = 0.0
!.......L%X,L%Y,L%DX,L%DY,L%DT ARE ALLOCATED AND
! INITIALIZED IN DX_CALC AND SUBGRID_MATCHING
END IF
IF (L%LAYSWITCH.EQ.0 .AND. L%LAYGOV.GE.2) THEN
ALLOCATE(L%ALPHA(L%NX,L%NY))
ALLOCATE(L%A1X(L%NX,L%NY))
ALLOCATE(L%A2X(L%NX,L%NY))
ALLOCATE(L%A1Y(L%NX,L%NY))
ALLOCATE(L%A2Y(L%NX,L%NY))
ALLOCATE(L%CF(L%NX,L%NY,4))
ALLOCATE(L%CB(L%NX,L%NY,4))
L%ALPHA = 1.0
L%CF = 0.0
L%CB = 0.0
END IF
IF (SWITCH .EQ. 1 .AND. L%BC_TYPE.EQ.1) THEN
ALLOCATE(L%SPONGE_COEFX(L%NX,L%NY))
ALLOCATE(L%SPONGE_COEFY(L%NX,L%NY))
L%SPONGE_COEFX = 0.0
L%SPONGE_COEFY = 0.0
ENDIF
IF (L%FRIC_SWITCH.EQ.3) THEN
L%SEDI_SWITCH = 0
ALLOCATE(L%DH(L%NX,L%NY,2))
L%DH = 0.0
ELSE
L%SEDI_SWITCH = 1
ENDIF
RETURN
END
!----------------------------------------------------------------------
SUBROUTINE SPHERICAL_PARAMETERS(LO,LA)
! CREATED BY TOM LOGAN, UAF (2004)
! ******** PARAMETERS FOR SPHERICAL COORD. *******
! Updated on Dec21 2008 (Xiaoming Wang, GNS)
!----------------------------------------------------------------------
USE LAYER_PARAMS
TYPE (LAYER) :: LO
TYPE (LAYER), DIMENSION(NUM_GRID) :: LA
COMMON /CONS/ ELMAX,GRAV,PI,R_EARTH,GX,EPS,ZERO,ONE,NUM_GRID, &
NUM_FLT,V_LIMIT,RAD_DEG,RAD_MIN
IF (LO%LAYSWITCH .EQ. 0 .AND. LO%LAYCORD .EQ. 0) THEN
WRITE (*,*) ' CALCULATING COEFFICIENTS FOR SPHERICAL COORD.'
CALL PARA(LO)
ENDIF
DO I=1,NUM_GRID
IF (LA(I)%LAYSWITCH .EQ. 0 .AND. LA(I)%LAYCORD .EQ. 0) THEN
CALL PARA(LA(I))
ENDIF
END DO
RETURN
END
!----------------------------------------------------------------------
SUBROUTINE PARA(LO)
!......................................................................
!DESCRIPTION:
! #. CALCULATE PARAMETERS FOR SPHERICAL COORD.
!NOTE:
! #. REVISED SIGNIFICANTLY BY TOM LOGAN, UAF (2005)
! #. UPDATED ON NOV.07, 2008 (XIAOMING WANG, GNS)
! LINES STARTING WITH "!!" ARE REMOVED
! LINES ENDED WITH "!!" ARE ADDED
! #. UPDATED ON FEB06 2009 (XIAOMING WANG, GNS)
!----------------------------------------------------------------------
USE LAYER_PARAMS
TYPE (LAYER) :: LO
COMMON /CONS/ ELMAX,GRAV,PI,R_EARTH,GX,EPS,ZERO,ONE,NUM_GRID, &
NUM_FLT,V_LIMIT,RAD_DEG,RAD_MIN
DATA CORI_W/7.2722E-5/, OSIXTY/0.016666666667/
DATA TWLVTH/0.08333333333333/
!.....PARAMETERS FOR SPHERICAL COORD IN GRID LAYER LO
IF (LO%LAYCORD .EQ. 0) THEN
IF (LO%LAYGOV.LE.1 .OR. LO%LAYGOV.EQ.9) THEN
!* DX = R_EARTH*LO%DX*RAD_MIN
!* RR = LO%DT/DX
RR = LO%DT/R_EARTH
!* RS = 0.5*G*RR
RS = GRAV*RR !!
RT = 0.5*LO%DT*CORI_W
DO J=1,LO%NY
ANG_M = LO%Y(J)*RAD_DEG
ANG_N = (LO%Y(J)+0.5*LO%DEL_Y(J)*OSIXTY)*RAD_DEG
COS_M = COS(ANG_M)
COS_N = COS(ANG_N)
SIN_M = SIN(ANG_M)
SIN_N = SIN(ANG_N)
DO I=1,LO%NX
LO%R0(I,J) = RR/(LO%DEL_Y(J)*RAD_MIN) !!
LO%R1(I,J) = RR/COS_M/(LO%DX*RAD_MIN)
LO%R11(I,J) = RR/COS_M/(LO%DEL_Y(J)*RAD_MIN)
IF (LO%LAYGOV.EQ.0) THEN
LO%R2(I,J) = RS/COS_M/(LO%DX*RAD_MIN) &
*AMAX1(LO%HP(I,J),ZERO)
ENDIF
IF (LO%LAYGOV.EQ.1) THEN
LO%R2(I,J) = RS/COS_M/(LO%DX*RAD_MIN)
ENDIF
LO%R21(I,J) = RR/COS_M/(LO%DX*RAD_MIN)
LO%R22(I,J) = RR/COS_N/(LO%DX*RAD_MIN)
LO%R3(I,J) = RT*SIN_M
IF (LO%LAYGOV.EQ.0) THEN
LO%R4(I,J) = RS/(LO%DEL_Y(J)*RAD_MIN) &
*AMAX1(LO%HQ(I,J),ZERO)
ENDIF
IF (LO%LAYGOV.EQ.1) THEN
LO%R4(I,J) = RS/(LO%DEL_Y(J)*RAD_MIN)
ENDIF
LO%R5(I,J) = RT*SIN_N
END DO
LO%R6(J) = COS_N
END DO
ELSEIF (LO%LAYGOV.EQ.2 .OR. LO%LAYGOV.EQ.3) THEN
!* DX = R_EARTH*LO%DX*RAD_MIN
!* RR = LO%DT/DX
RR = LO%DT/R_EARTH
!* RS = 0.5*G*RR
RS = GRAV*RR !!
RT = 0.5*LO%DT*CORI_W
DO J=1,LO%NY
ANG_M = LO%Y(J)*RAD_DEG
ANG_N = (LO%Y(J)+0.5*LO%DEL_Y(J)*OSIXTY)*RAD_DEG
COS_M = COS(ANG_M)
COS_N = COS(ANG_N)
SIN_M = SIN(ANG_M)
SIN_N = SIN(ANG_N)
DO I=1,LO%NX
IF (LO%MASK(I,J).EQ.0) THEN
DX = LO%DX
DY = LO%DEL_Y(J)
ENDIF
IF (LO%MASK(I,J).EQ.1) THEN
DX = LO%DX*LO%ALPHA(I,J)
DY = LO%DEL_Y(J)*LO%ALPHA(I,J)
ENDIF
LO%R0(I,J) = RR/(DY*RAD_MIN) !!
LO%R1(I,J) = RR/COS_M/(DX*RAD_MIN)
LO%R11(I,J) = RR/COS_M/(DY*RAD_MIN)
IF (LO%LAYGOV.EQ.2) THEN
LO%R2(I,J) = RS/COS_M/(DX*RAD_MIN) &
*AMAX1(LO%HP(I,J),ZERO)
ENDIF
IF (LO%LAYGOV.EQ.3) THEN
LO%R2(I,J) = RS/COS_M/(DX*RAD_MIN)
ENDIF
LO%R21(I,J) = RR/COS_M/(DX*RAD_MIN)
LO%R22(I,J) = RR/COS_N/(DX*RAD_MIN)
LO%R3(I,J) = RT*SIN_M
IF (LO%LAYGOV.EQ.2) THEN
LO%R4(I,J) = RS/(DY*RAD_MIN) &
*AMAX1(LO%HQ(I,J),ZERO)
ENDIF
IF (LO%LAYGOV.EQ.3) THEN
LO%R4(I,J) = RS/(DY*RAD_MIN)
ENDIF
LO%R5(I,J) = RT*SIN_N
END DO
LO%R6(J) = COS_N
END DO
ENDIF
ENDIF
RETURN
END
!----------------------------------------------------------------------
SUBROUTINE SED_TRANSPORT (LO)
!.....TENTATIVE TRY TO INCLUDING SEDIMENT TRANSPORT MODEL (APR. 2007)
! SEDIMENT FLUX IS COMPUTED ACCORDING TO RIBBERINK (1998)
!----------------------------------------------------------------------
USE LAYER_PARAMS
TYPE (LAYER) :: LO
REAL D, DS, THETA_C, THETA_P,THETA_Q,G,RHO,RHO_S
REAL H(LO%NX,LO%NY),DH(LO%NX,LO%NY)
REAL QX(LO%NX,LO%NY),QY(LO%NX,LO%NY)
CHARACTER(LEN=30) FNAME
G = 9.807
RHO = 1000.0
RHO_S = 2650.0
DS = 0.0002
POROSITY = 0.3
THETA_C = 0.04
S = RHO_S/RHO
! Q(:,:) = 0.0
QX(:,:) = 0.0
QY(:,:) = 0.0
! H(:,:) = -LO%H(:,:)
DH(:,:) = LO%H(:,:) + LO%Z(:,:,1)
DO I = 1,LO%NX
IM1 = I-1
IP1 = I+1
IF (IM1.LE.0) IM1=1
IF (IP1.GE.LO%NX) IP1=LO%NX
DO J = 1,LO%NY
JM1 = J-1
JP1 = J+1
IF (JM1.LE.1) JM1=1
IF (JP1.GE.LO%NY) JP1=LO%NY
!TOTAL WATER DEPTH AT DISCHARGE LOCATION P
HP = 0.5*(DH(I,J)+DH(IP1,J))
!TOTAL WATER DEPTH AT DISCHARGE LOCATION Q
HQ = 0.5*(DH(I,J)+DH(I,JP1))
!SEDIMENT FLUX AT DISCHARGE LOCATION P
IF (HP .GT. 0.05) THEN
COEF_P = RHO*G*(LO%FRIC_COEF**2)/(HP**2.33333)
!Y COMPONENT OF VOLUME FLUX AT DISCHARGE LOCATION P
XQQ = 0.25*(LO%N(I,J,1)+LO%N(IP1,J,1) &
+ LO%N(I,JM1,1)+LO%N(IP1,JM1,1))
!TOTAL VOLUME FLUX AT DISCHARGE LOCATION P
FLUX_P = SQRT(LO%M(I,J,1)**2+XQQ**2)
!X COMPONENT OF BOTTOM SHEAR STRESS AT DISCHARGE LOCATION P
TPX = COEF_P*LO%M(I,J,1)*FLUX_P
!Y COMPONENT OF BOTTOM SHEAR STRESS AT DISCHARGE LOCATION P
TPY = COEF_P*XQQ*FLUX_P
! TOTAL SHEAR STRESS AT DISCHARGE LOCATION P
TPB = SQRT(TPX**2+TPY**2)
!SHIELDS PARAMETER AT DISCHARGE LOCATION P
THETA_P = TPB/((RHO_S-RHO)*G*DS)
IF (THETA_P .GT. THETA_C) THEN
!RIBBERINK (1998)
QS_P = 11.0*SQRT((S-1.0)*G*DS**3) &
* ((THETA_P-THETA_C)**1.65)
IF (QS_P.GT.0.02) QS_P = 0.02
!X COMPONENT OF SAND FLUX AT DISCHARGE LOCATION P
IF (FLUX_P .GT. 1.0E-10) &
QX(I,J) = QS_P * LO%M(I,J,1)/FLUX_P
ENDIF
ENDIF
!SEDIMENT FLUX AT DISCHARGE LOCATION Q
IF (HQ .GT. 0.05) THEN
COEF_Q = RHO*G*(LO%FRIC_COEF**2)/(HQ**2.33333)
XPP = 0.25*(LO%M(I,J,1)+LO%M(I,JP1,1)+LO%M(IM1,J,1) &
+ LO%M(IM1,JP1,1))
FLUX_Q = SQRT(XPP**2+LO%N(I,J,1)**2)
TQX = COEF_Q*XPP*FLUX_Q
TQY = COEF_Q*LO%N(I,J,1)*FLUX_Q
TQB = SQRT(TQX**2+TQY**2)
THETA_Q = TQB/((RHO_S-RHO)*G*DS)
IF (THETA_Q .GT. THETA_C) THEN
!RIBBERINK (1998)
QS_Q = 11.0*SQRT((S-1.0)*G*DS**3) &
*((THETA_Q-THETA_C)**1.65)
IF (QS_Q.GT.0.02) QS_Q = 0.02
IF (FLUX_Q .GT. 1.0E-10) &
QY(I,J) = QS_Q*LO%N(I,J,1)/FLUX_Q
ENDIF
ENDIF
ENDDO
ENDDO
DO I = 2,LO%NX-1
DO J = 2,LO%NY-1
LO%DH(I,J,2) = LO%DH(I,J,1)-LO%RX*(1/(1-POROSITY)) &
*(QX(I,J)-QX(I-1,J)+QY(I,J)-QY(I,J-1))
ENDDO
ENDDO
LO%H(:,:) = LO%H(:,:)-LO%DH(:,:,2)
RETURN
END