c
c
C.... File Name .............. GXSTRA.FOR ....................... 271107
SUBROUTINE GXSTRA(IDUMMY)
C-----------------------------------------------------------------------
C SUBROUTINE GXSTRA is called from EARTH once only, in order to set
c indices and print warnings if necessary.
C
C-----------------------------------------------------------------------
INCLUDE 'satear'
INCLUDE 'grdear'
COMMON/GENI/NXNY,IGFIL1(48),ITEM1,IGFIL2(5),IPRPS,IGFIL3(4)
c COMMON/ISTRA2/LBEX,LBEY,LBEZ,LBET,LBSX,LBSY,LBSZ,LBSXY,LBSYZ,LBSZX
c VIA : 29.06.08
COMMON/ISTRA2/LBEX,LBEY,LBEZ,LBET,LBSX,LBSY,LBSZ,LBSXY,LBSYZ,LBSXZ
COMMON/XYZCON/XCONST,YCONST,ZCONST
COMMON/NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
COMMON/DBE/dbear
LOGICAL dbear
COMMON/IGRND/IG(100)
C
c----VIA 06.05.07 :
COMMON /STRA05/ LCDM, ISLDIM, BCONST1, BCONST2, LCDMTIME,
1 cdmDTime1, cdmDTime2
LOGICAL LCDM, LCDMTIME
COMMON /STRA07/ KSU_LOWX, KSU_LOWY, KSU_LOWZ
COMMON/USPGLOB/USP, USPDBG
LOGICAL USP,USPDBG
C*****************************************************************
NAMSUB = 'GXSTRA'
C
C--- GROUP 1. Run title and other preliminaries
C
IF(IGR.EQ.1.AND.ISC.EQ.1) THEN
IF(IPRPS.EQ.0) THEN
CALL WRIT40('PRPS must be stored for solid stress')
CALL SET_ERR(512,
1 'Error. PRPS must be stored for solid stress.',1)
RETURN
ENDIF
XCONST=0.0 ! constants indicating restraining force
YCONST=0.0 ! exerted by container when NX, NY or NZ = 1
ZCONST=0.0 ! Default zero signifies no restraint
IF(NX.EQ.1) CALL GETSDR('BOUNDARY','XCONST',XCONST)
IF(NY.EQ.1) CALL GETSDR('BOUNDARY','YCONST',YCONST)
IF(NZ.EQ.1) CALL GETSDR('BOUNDARY','ZCONST',ZCONST)
IF(.NOT.NULLPR) THEN
CALL WRIT40('Stress-in-solids-related data')
CALL WRIT40('boundary-stress constants')
CALL WRIT3R('xconst ',xconst,'yconst ',yconst,
1 'zconst ',zconst)
CALL WRITST
CALL WRITBL
ENDIF
c
LBEX=LBNAME('EPSX') ! determine which of the post-processed
LBEY=LBNAME('EPSY') ! variables are required for printing or
LBEZ=LBNAME('EPSZ') ! plotting
LBET=LBNAME('EPST')
IF(LBET.GT.0) THEN
IF(ITEM1.LE.0) THEN
WRITE(14,*) 'store tem1 if thermal stress needed'
WRITE(14,*) 'if not, do not store epst'
WRITE(6,*) 'store tem1 if thermal stress needed'
WRITE(6,*) 'if not, do not store epst'
CALL SET_ERR(513, 'Error. See result file.',1)
RETURN
ENDIF
ENDIF
c----VIA 16.04.07 ---------------------------------------------------
LCDM = LBNAME('DISX').NE.0.OR.LBNAME('DISY').NE.0
1 .OR.LBNAME('DISZ').NE.0
if(.NOT.LCDM.AND.USP) then
WRITE(14,*) 'USP: For STRA=T CDM model is needed'
CALL SET_ERR(1050, 'USP: For STRA=T CDM model is needed',1)
RETURN
endif
if(LCDM) then
if(NZ.GT.1) then
IF(LBEX.LE.0.OR.LBEY.LE.0.OR.LBEZ.LE.0) THEN
WRITE(14,*) 'store EPSX/EPSY/EPSZ is needed'
CALL SET_ERR(513, 'Error. See result file.',1)
RETURN
ENDIF
endif
else
IF(LBEX.LE.0) THEN
WRITE(14,*) 'store EPSX is needed'
CALL SET_ERR(513, 'Error. See result file.',1)
RETURN
ENDIF
IF(LBEY.LE.0) THEN
WRITE(14,*) 'store EPSY is needed'
CALL SET_ERR(513, 'Error. See result file.',1)
RETURN
ENDIF
IF(LBEZ.LE.0) THEN
WRITE(14,*) 'store EPSZ is needed'
CALL SET_ERR(513, 'Error. See result file.',1)
RETURN
ENDIF
endif
LBSX=LBNAME('STRX')
LBSY=LBNAME('STRY')
LBSZ=LBNAME('STRZ')
LBSXY=LBNAME('STXY')
LBSYZ=LBNAME('STYZ')
LBSXZ=LBNAME('STXZ')
LCDMTIME = (LBNAME('VDSX').NE.0.OR.LBNAME('VDSY').NE.0
1 .OR.LBNAME('VDSZ').NE.0) .AND.
2 (.NOT.STEADY).AND. LCDM
IF(.NOT.NULLPR) THEN
CALL WRIT40('Stress-in-solids model:')
IF(LCDM) THEN
CALL WRIT40('Collocated-displacements method')
IF(LCDMTIME)
1 CALL WRIT40('and Time-dependent term (T-scheme)')
ELSE
CALL WRIT40('Staggered-displacements method')
ENDIF
CALL WRITST
ENDIF
c----- calculate variable for function CorrZ
IF(LCDM) THEN
ISLDIM = 0
BCONST1 = 0.0
BCONST2 = 0.0
if(NX.EQ.1.AND. NY.NE.1.AND.NZ.NE.1) then
if(XCONST.LE.0.9E20) then
ISLDIM = 1
BCONST1 = XCONST*1.e-11
endif
elseif(NY.EQ.1.AND. NX.NE.1.AND.NZ.NE.1) then
if(YCONST.LE.0.9E20) then
ISLDIM = 1
BCONST1 = YCONST*1.e-11
endif
elseif(NZ.EQ.1.AND. NX.NE.1.AND.NY.NE.1) then
if(ZCONST.LE.0.9E20) then
ISLDIM = 1
BCONST1 = ZCONST*1.e-11
endif
elseif(NX.EQ.1.AND.NY.EQ.1 .AND.NZ.NE.1) then
if(XCONST.LE.0.9E20.and.YCONST.LE.0.9E20) then
ISLDIM = 2
BCONST1 = XCONST*1.e-11
BCONST2 = YCONST*1.e-11
elseif(XCONST.GE.0.9E20.and.YCONST.LE.0.9E20) then
ISLDIM = 1
BCONST1 = YCONST*1.e-11
elseif(XCONST.LE.0.9E20.and.YCONST.GE.0.9E20) then
ISLDIM = 1
BCONST1 = XCONST*1.e-11
endif
elseif(NX.EQ.1.AND.NY.NE.1 .AND.NZ.EQ.1) then
if(XCONST.LE.0.9E20.and.ZCONST.LE.0.9E20) then
ISLDIM = 2
BCONST1 = XCONST*1.e-11
BCONST2 = ZCONST*1.e-11
elseif(XCONST.GE.0.9E20.and.ZCONST.LE.0.9E20) then
ISLDIM = 1
BCONST1 = ZCONST*1.e-11
elseif(XCONST.LE.0.9E20.and.ZCONST.GE.0.9E20) then
ISLDIM = 1
BCONST1 = XCONST*1.e-11
endif
elseif(NX.NE.1.AND.NY.EQ.1 .AND.NZ.EQ.1) then
if(YCONST.LE.0.9E20.and.ZCONST.LE.0.9E20) then
ISLDIM = 2
BCONST1 = YCONST*1.e-11
BCONST2 = ZCONST*1.e-11
elseif(YCONST.GE.0.9E20.and.ZCONST.LE.0.9E20) then
ISLDIM = 1
BCONST1 = ZCONST*1.e-11
elseif(YCONST.LE.0.9E20.and.ZCONST.GE.0.9E20) then
ISLDIM = 1
BCONST1 = ZCONST*1.e-11
endif
endif
if(USP) CALL USPSTRINIT(1)
ENDIF
ELSEIF(IGR.EQ.1.AND.ISC.EQ.3) THEN ! slab-wise storage
if(LCDM.AND.NZ.GT.1) then
if(USP) RETURN
N = NX*NY
if(NX.GT.1) CALL GXMAKE(KSU_LOWX,N,'SL_X')
if(NY.GT.1) CALL GXMAKE(KSU_LOWY,N,'SL_Y')
CALL GXMAKE(KSU_LOWZ,N,'SL_Z')
endif
ENDIF
C********************************************************************
C
C--- GROUP 19. Special calls to GROUND from EARTH
C
NAMSUB = 'gxstra'
END
***********************************************************************
c
c STRCAL is called from EARTH at print-out time in order to compute
c strains and stresses from the within-EARTH computed displacements.
C
C The displacements use the storage locations which, for the fluid
C region, are occupied by the velocities.
SUBROUTINE STRCAL
C----------------------------------------------------------------------
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'grdear'
INCLUDE 'prpcmn'
COMMON/F01/M0F(600)
cc COMMON/F0/IFIL1(33),KXYDZG,IFIL2(5),KXYDXG,KXYDYG,IFIL3(263)
COMMON/F0/IFIL0(7),KYR,IFIL1(25),KXYDZG,IFIL2(5),KXYDXG,KXYDYG,
1 IFIL3(263)
COMMON/GENI/NXNY,IGFIL1(5),NYST,IGFIL2(2),NFM,IGFIL3(50)
COMMON /NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
COMMON/ISTRA2/LBEX,LBEY,LBEZ,LBET,LBSX,LBSY,LBSZ,LBSXY,LBSYZ,LBSXZ
LOGICAL SLD
REAL LAMDA
c----VIA 16.04.07 :
COMMON /STRA05/ LCDM, ISLDIM, BCONST1, BCONST2, LCDMTIME,
1 cdmDTime1, cdmDTime2
LOGICAL LCDM, LCDMTIME
c----VIA 16.04.07 :
NAMSUB = 'STRCAL'
PRINTNUL=-1.234E-11 ! this prints as 'none'
C----------------------------------------------------------------------
C.... Compute zero locations, and other preliminaries.
C----------------------------------------------------------------------
if(LCDM) then
CALL DilateCDM
return
endif
CALL SUB4(L0EX,0,L0EY,0,L0EZ,0,L0ET,0)
CALL SUB4(L0SX,0,L0SY,0,L0SZ,0,L0TEM,0)
CALL SUB3(L0SXY,0,L0SYZ,0,L0SXZ,0)
IF(LBEX.NE.0) L0EX = L0F(LBEX)
IF(LBEY.NE.0) L0EY = L0F(LBEY)
IF(LBEZ.NE.0) L0EZ = L0F(LBEZ)
IF(LBET.NE.0) L0ET = L0F(LBET)
IF(LBSX.NE.0) L0SX = L0F(LBSX)
IF(LBSY.NE.0) L0SY = L0F(LBSY)
IF(LBSZ.NE.0) L0SZ = L0F(LBSZ)
IF(LBSXY.NE.0) L0SXY = L0F(LBSXY)
IF(LBSYZ.NE.0) L0SYZ = L0F(LBSYZ)
IF(LBSXZ.NE.0) L0SXZ = L0F(LBSXZ)
IF(IDRH1.NE.0) L0DRH1= L0F(IDRH1)
IF((L0EX+L0EY+L0EZ+L0ET+L0SX+L0SY+L0SZ).EQ.0) RETURN
C----------------------------------------------------------------------
C.... Calculate normal stresses.
C----------------------------------------------------------------------
DO IDIR =1,3 ! loop over all three directions
! direct stresses
IF(IDIR.EQ.1) THEN ! x-direction.
IF(L0SX.NE.0.AND.L0EX.NE.0) THEN
L0S=L0SX
L0E=L0EX
GO TO 10 ! compute stresses
ELSE
GO TO 11 ! skip
ENDIF
ELSEIF(IDIR.EQ.2) THEN ! y-direction.
IF(L0SY.NE.0.AND.L0EY.NE.0) THEN
L0S=L0SY
L0E=L0EY
GO TO 10 ! compute stresses
ELSE
GO TO 11 ! skip
ENDIF
ELSEIF(IDIR.EQ.3) THEN ! Z-direction.
IF(L0SZ.NE.0.AND.L0EZ.NE.0) THEN
L0S=L0SZ
L0E=L0EZ
GO TO 10 ! compute stresses
ELSE
GO TO 11 ! skip
ENDIF
ENDIF
10 DO I=1,NXNY ! compute stresses
IF(SLD(I)) THEN
DILSTR = LAMDA(I)*F(M0F(1)+I)
EPSSTR = TWOG(I)*F(L0E+I)
STRESS = DILSTR + EPSSTR
TESTSTR= ABS(DILSTR)+ABS(EPSSTR)
IF(LBET.GT.0) THEN
THESTR = - AITCH(I)*T1VALI(I)
STRESS = STRESS + THESTR
TESTSTR= TESTSTR + ABS(THESTR)
ENDIF
IF(ABS(STRESS).LT.1.0E-4*TESTSTR) ! set to
1 STRESS=0.0 ! zero resultant stresses when affected by round-off
F(L0S+I) = STRESS * 1.E11 ! restore correct magnitude
ELSE
F(L0S+I) = PRINTNUL
F(L0E+I) = PRINTNUL
ENDIF
ENDDO
11 CONTINUE
! shear stresses
IF(IDIR.EQ.1) THEN ! x-direction.
IF(L0SYZ.NE.0) THEN
c L0NUM1=M0F(5) ! NUM refers to numerator, always a displacement
c L0NUM2=M0F(7)
c L0DEN1=KXYDZG ! NUM refers to denominator, always a distance
c L0DEN2=KXYDYG ! interval
c Str = dV/DZ+dW/dY
c IADD1=NYST
c IADD2=NFM
L0NUM2=M0F(5) ! NUM refers to numerator, always a displacement
L0NUM1=M0F(7)
L0DEN2=KXYDZG ! NUM refers to denominator, always a distance
L0DEN1=KXYDYG ! interval
IADD2=NFM
IADD1=1
CALL CONNXY(KXYDZG,DZG)
L0S=L0SYZ
GO TO 20 ! compute stresses
ELSE
GO TO 21 ! skip
ENDIF
ELSEIF(IDIR.EQ.2) THEN ! y-direction.
IF(L0SXZ.NE.0) THEN
L0NUM1=M0F(7)
L0NUM2=M0F(3)
L0DEN1=KXYDXG
L0DEN2=KXYDZG
c Str = dU/DZ+dW/dX
c IADD1=NFM
c IADD2=1
IADD1=NY
IADD2=NFM
CALL CONNXY(KXYDZG,DZG)
L0S=L0SXZ
GO TO 20 ! compute stresses
ELSE
GO TO 21 ! skip
ENDIF
ELSEIF(IDIR.EQ.3) THEN ! Z-direction.
IF(L0SXY.NE.0) THEN
c Str = dU/DY+dV/dX
L0NUM1=M0F(3)
L0NUM2=M0F(5)
L0DEN1=KXYDYG
L0DEN2=KXYDXG
IADD1=1
IADD2=NY
L0S=L0SXY
GO TO 20 ! compute stresses
ELSE
GO TO 21 ! skip
ENDIF
ENDIF
c 20 DO I=1,NXNY ! compute stresses
c IF(SLD(I).AND.SLD(I+IADD1)) THEN
c F(L0S+I) = ( F(L0NUM1+I+IADD1) - F(L0NUM1+I) )/F(L0DEN1+I)
c 1 + ( F(L0NUM2+I+IADD2) - F(L0NUM2+I) )/F(L0DEN2+I)
c F(L0S+I) = F(L0S+I) * ONEG(I) * 1.E11
c ELSE
20 CONTINUE
IF(L0NUM1.EQ.0.OR.L0NUM2.EQ.0) GO TO 21
DO IX=1,NX
IADD = (IX-1)*NY
DO IY=1,NY
I = IY+IADD
F(L0S+I) = PRINTNUL
IF(IDIR.EQ.3.AND.(IX.EQ.NX.OR.IY.EQ.NY)) GO TO 30
IF(IDIR.EQ.1.AND.(IZ.EQ.NZ.OR.IY.EQ.NY)) GO TO 30
IF(IDIR.EQ.2.AND.(IZ.EQ.NZ.OR.IX.EQ.NX)) GO TO 30
Stress = 0.0
NoDiriv = 0
IF(SLD(I).AND.SLD(I+IADD1)) THEN
NoDiriv = 1
Stress = ( F(L0NUM1+I+IADD1) - F(L0NUM1+I) )/F(L0DEN1+I)
if(.NOT.CARTES.AND.IDIR.EQ.3) then
Stress = Stress - 0.5*(F(L0NUM1+I+IADD1)/F(KYR+IY+1)
1 +F(L0NUM1+I)/F(KYR+IY))
endif
ENDIF
IF(SLD(I).AND.SLD(I+IADD2)) THEN
NoDiriv = NoDiriv+1
Stress = Stress+
1 (F(L0NUM2+I+IADD2)-F(L0NUM2+I))/F(L0DEN2+I)
ENDIF
IF(NoDiriv.NE.0) F(L0S+I) = Stress * ONEG(I) * 1.E11
30 CONTINUE
ENDDO
ENDDO
21 CONTINUE
ENDDO
NAMSUB = 'STRCAL'
END
C
C-----------------------------------------------------------------------
c