cgxstra.for c c

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 '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/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 '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/grdear'
      INCLUDE '/phoenics/d_includ/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