str_mod.f90

   1 MODULE STRAIN_CALC
   2
   3   IMPLICIT NONE
   4
   5 CONTAINS
   6
   7    SUBROUTINE STRISO(EESUM,EEDIF,EEZZ,EERZ,CHI)
   8         Use Input_Data
   9      IMPLICIT NONE
  10      REAL, INTENT(IN   ) :: CHI
  11      REAL, INTENT(  OUT) :: EESUM,EEDIF,EEZZ,EERZ
  12
  13      DOUBLE PRECISION :: AASUM,AADIF,AARZ
  14
  15      !!! COMMON
  16      REAL :: RHO,ZETA,ETA
  17      COMMON /QAGON/RHO,ZETA,ETA
  18
  19      ETA=1.E0
  20      CALL PREST(AASUM,AADIF,AARZ)
  21
  22      EESUM = - BIAUX*(AASUM+CHI)+BISUM*CHI
  23      EEZZ  = + BIAUX*(AASUM+CHI)+BIZZ*CHI
  24      EEDIF = BIAUX*(AADIF-CHI)
  25      EERZ  = BIAUX*AARZ
  26
  27    END SUBROUTINE STRISO
  28
  29    SUBROUTINE STRANISO(EESUM,EEDIF,EEZZ,EERZ,CHI)
  30         Use Input_Data
  31      IMPLICIT NONE
  32      REAL, INTENT(IN   ) :: CHI
  33      REAL, INTENT(  OUT) :: EESUM,EEDIF,EEZZ,EERZ
  34
  35      DOUBLE PRECISION :: AASUM1,AADIF1,AARZ1,&
  36                          AASUM2,AADIF2,AARZ2
  37
  38      !!! COMMON
  39      REAL :: RHO,ZETA,ETA
  40      COMMON /QAGON/RHO,ZETA,ETA
  41
  42      ETA=ETA1
  43      CALL PREST(AASUM1,AADIF1,AARZ1)
  44      ETA=ETA2
  45      CALL PREST(AASUM2,AADIF2,AARZ2)
  46
  47      EESUM = (2.E0*EPSA + CN1_2G)*CHI + &
  48              (CN02_1*AASUM2 - CN02_2*AASUM1)/ETA_DIF
  49      EEDIF = +CN1_2G*CHI - &
  50              (CN02_1*AADIF2 - CN02_2*AADIF1)/ETA_DIF
  51      EEZZ  = EPSC*CHI + &
  52              (CN42_1*AASUM2 - CN42_2*AASUM1)/ETA_DIF
  53      EERZ  = (CN31_1*AARZ2 - CN31_2*AARZ1)/ETA_DIF
  54
  55 !    WRITE(12,'(10(E15.8,1X))')ZETA,RHO,AASUM1,AASUM2,AADIF1,AADIF2
  56
  57      RETURN
  58
  59    END SUBROUTINE STRANISO
  60
  61    SUBROUTINE PREST(AASUM,AADIF,AARZ)
  62
  63       Use Input_Data, ONLY: RD,RC,ZC
  64
  65    IMPLICIT NONE
  66
  67 !! 'dummy' and local variables !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  68
  69    DOUBLE PRECISION AASUM,AADIF,AARZ,RO,ZE
  70
  71 !! GAG Variables  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  72
  73    INTEGER, PARAMETER :: LIMITE=100
  74    DOUBLE PRECISION AL,AU,ABSERR,EPSABS,EPSREL,WORK(4*LIMITE)
  75    INTEGER IER,IWORK(LIMITE),KEY,LAST,LENW,LIMIT,NEVAL
  76
  77    DOUBLE PRECISION, EXTERNAL :: I_FA00,I_FA10,I_FA20
  78
  79    !!! COMMON
  80    REAL :: RHO,ZETA,ETA
  81    COMMON /QAGON/RHO,ZETA,ETA
  82
  83 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  84
  85     RO=RHO
  86     ZE=ZETA
  87     AL=0.D0
  88     AU=DBLE(RD)
  89
  90     EPSABS=0.D0
  91     EPSREL=1.D-4
  92     KEY=6
  93     LIMIT=LIMITE
  94     LENW=LIMIT*4
  95
  96     CALL DQAG(I_FA00,AL,AU,EPSABS,EPSREL,KEY,AASUM,ABSERR,NEVAL, &
  97              IER,LIMIT,LENW,LAST,IWORK,WORK)
  98     IF (IER.NE.0) &
  99        CALL DQAG_CHECK("FA00","AASUM",AASUM,ABSERR,RO*RC,ZE*ZC,NEVAL,IER,LAST)
 100
 101
 102     CALL DQAG(I_FA10,AL,AU,EPSABS,EPSREL,KEY,AARZ,ABSERR,NEVAL, &
 103              IER,LIMIT,LENW,LAST,IWORK,WORK)
 104
 105     IF (IER.NE.0) &
 106        CALL DQAG_CHECK("FA10","AARZ",AARZ,ABSERR,RO*RC,ZE*ZC,NEVAL,IER,LAST)
 107
 108     CALL DQAG(I_FA20,AL,AU,EPSABS,EPSREL,KEY,AADIF,ABSERR,NEVAL, &
 109              IER,LIMIT,LENW,LAST,IWORK,WORK)
 110
 111     IF (IER.NE.0) &
 112        CALL DQAG_CHECK("FA20","AADIF",AADIF,ABSERR,RO*RC,ZE*ZC,NEVAL,IER,LAST)
 113
 114    RETURN
 115    END SUBROUTINE PREST
 116
 117    SUBROUTINE DQAG_CHECK(FUNC,RESNAME,RES,ABSERR,RO,ZE,NEVAL,IER,LAST)
 118    IMPLICIT NONE
 119    CHARACTER(LEN=*)    :: FUNC,RESNAME
 120    DOUBLE PRECISION    :: RES,ABSERR,RO,ZE
 121    INTEGER             :: NEVAL,IER,LAST
 122      WRITE(6,*)"ERROR AT: RHO: ",RO," ZETA: ",ZE
 123      WRITE(6,*)"QAG ",TRIM(FUNC)," ",TRIM(RESNAME)," RESULT=",RES
 124      WRITE(6,*)"QAG ABSERR=",ABSERR
 125      WRITE(6,*)"QAG NEVAL=",NEVAL
 126      WRITE(6,*)"QAG IER=",IER
 127      WRITE(6,*)"QAG LAST=",LAST
 128
 129      RETURN
 130    END SUBROUTINE DQAG_CHECK
 131
 132   FUNCTION FA00(RHO,ZETA,RHOP,ZM,ETA) RESULT(A00)
 133         Use Input_Data, ONLY: RC,ZC
 134         Use ELIPT_MOD
 135
 136     IMPLICIT NONE
 137     DOUBLE PRECISION, INTENT(IN) :: RHO,ZETA,RHOP,ZM,ETA
 138     DOUBLE PRECISION :: A00
 139
 140     DOUBLE PRECISION ::  PI,DPI
 141     DOUBLE PRECISION :: SRRHOP,RRHOP15,SM,S0,DELTAM,DELTA0,&
 142                         KM,K0,CKM,CK0,KKM,KK0,EKM,EK0
 143
 144       PI  = DATAN(1.D0)*4.D0
 145       DPI = 2.D0*PI
 146
 147       SRRHOP=SQRT(RHO*RHOP)
 148       RRHOP15=SQRT(RHO*RHOP)**3
 149
 150       SM = DSIGN(1.D0,ZETA-ZM)
 151       S0 = DSIGN(1.D0,ZETA)
 152
 153       DELTAM = DSQRT(ETA)*DABS(ZETA-ZM)*ZC/RC
 154       DELTA0 = DSQRT(ETA)*DABS(ZETA)*ZC/RC
 155
 156       KM = 2.D0*DSQRT( RHO*RHOP/((RHO+RHOP)**2+DELTAM**2) )
 157       K0 = 2.D0*DSQRT( RHO*RHOP/((RHO+RHOP)**2+DELTA0**2) )
 158
 159       CKM = (1.D0-KM**2)
 160       CK0 = (1.D0-K0**2)
 161
 162       CALL ELIPT0(KM,KKM,EKM)
 163       CALL ELIPT0(K0,KK0,EK0)
 164
 165        IF (RHO .EQ. 0.D0) THEN
 166           A00 = ( SM*DELTAM/SQRT(RHOP**2+DELTAM**2)**3) &
 167                -(S0*DELTA0/SQRT(RHOP**2+DELTA0**2)**3 )
 168           A00 = A00/2.D0
 169           RETURN
 170        END IF
 171
 172       A00 = SM*DELTAM*KM**3/(2.D0*DPI*CKM*RRHOP15)*EKM &
 173            -S0*DELTA0*K0**3/(2.D0*DPI*CK0*RRHOP15)*EK0
 174
 175       A00 = A00/2.D0
 176
 177       RETURN
 178
 179
 180   END FUNCTION FA00
 181
 182   FUNCTION FA10(RHO,ZETA,RHOP,ZM,ETA) RESULT(A10)
 183         Use Input_Data, ONLY: RC,ZC
 184         Use ELIPT_MOD
 185
 186     IMPLICIT NONE
 187     DOUBLE PRECISION, INTENT(IN) :: RHO,ZETA,RHOP,ZM,ETA
 188     DOUBLE PRECISION :: A10
 189
 190     DOUBLE PRECISION ::  PI,DPI
 191     DOUBLE PRECISION :: SRRHOP,RRHOP15,SM,S0,DELTAM,DELTA0,&
 192                         KM,K0,CKM,CK0,KKM,KK0,EKM,EK0
 193
 194       PI  = DATAN(1.D0)*4.D0
 195       DPI = 2.D0*PI
 196
 197       SRRHOP=DSQRT(RHO*RHOP)
 198       RRHOP15=DSQRT(RHO*RHOP)**3
 199
 200       SM = DSIGN(1.D0,ZETA-ZM)
 201       S0 = DSIGN(1.D0,ZETA)
 202
 203       DELTAM = DSQRT(ETA)*DABS(ZETA-ZM)*ZC/RC
 204       DELTA0 = DSQRT(ETA)*DABS(ZETA)*ZC/RC
 205
 206       KM = 2.D0*DSQRT( RHO*RHOP/((RHO+RHOP)**2+DELTAM**2) )
 207       K0 = 2.D0*DSQRT( RHO*RHOP/((RHO+RHOP)**2+DELTA0**2) )
 208
 209       CKM = (1.D0-KM**2)
 210       CK0 = (1.D0-K0**2)
 211
 212       CALL ELIPT0(KM,KKM,EKM)
 213       CALL ELIPT0(K0,KK0,EK0)
 214
 215        IF (RHO .EQ. 0.D0) THEN
 216           A10 = 0.D0
 217           RETURN
 218        END IF
 219
 220        A10 = KM/(4.D0*DPI*CKM*RHO*RRHOP15)  &
 221                * ( KM**2*(RHO**2-RHOP**2-DELTAM**2)*EKM   &
 222                  + 4.D0*RHO*RHOP*CKM*KKM )   &
 223           - K0/(4.D0*DPI*CK0*RHO*RRHOP15)  &
 224                * ( K0**2*(RHO**2-RHOP**2-DELTA0**2)*EK0 &
 225                  + 4.D0*RHO*RHOP*CK0*KK0 )
 226
 227        A10 = A10/2.D0
 228
 229       RETURN
 230
 231   END FUNCTION FA10
 232
 233   FUNCTION FA20(RHO,ZETA,RHOP,ZM,ETA) RESULT(A20)
 234         Use Input_Data, ONLY: RC,ZC
 235         Use ELIPT_MOD
 236
 237     IMPLICIT NONE
 238     DOUBLE PRECISION, INTENT(IN) :: RHO,ZETA,RHOP,ZM,ETA
 239     DOUBLE PRECISION :: A20
 240
 241     DOUBLE PRECISION ::  PI,DPI
 242     DOUBLE PRECISION :: SRRHOP,RRHOP15,SM,S0,DELTAM,DELTA0,&
 243                         KM,K0,CKM,CK0,KKM,KK0,EKM,EK0
 244
 245     DOUBLE PRECISION ::  SIRHO, PSIM, PSI0
 246
 247       PI  = DATAN(1.D0)*4.D0
 248       DPI = 2.D0*PI
 249
 250       SRRHOP=DSQRT(RHO*RHOP)
 251       RRHOP15=DSQRT(RHO*RHOP)**3
 252
 253       SM = DSIGN(1.D0,ZETA-ZM)
 254       S0 = DSIGN(1.D0,ZETA)
 255
 256       DELTAM = DSQRT(ETA)*DABS(ZETA-ZM)*ZC/RC
 257       DELTA0 = DSQRT(ETA)*DABS(ZETA)*ZC/RC
 258
 259       SIRHO = DSIGN(1.D0,RHO-RHOP)
 260
 261       KM = 2.D0*DSQRT( RHO*RHOP/((RHO+RHOP)**2+DELTAM**2) )
 262       K0 = 2.D0*DSQRT( RHO*RHOP/((RHO+RHOP)**2+DELTA0**2) )
 263
 264       CKM = (1.D0-KM**2)
 265       CK0 = (1.D0-K0**2)
 266
 267       CALL ELIPT0(KM,KKM,EKM)
 268       CALL ELIPT0(K0,KK0,EK0)
 269
 270       PSIM = DASIN( DELTAM/DSQRT((RHO-RHOP)**2+DELTAM**2) )
 271       PSI0 = DASIN( DELTA0/DSQRT((RHO-RHOP)**2+DELTA0**2) )
 272
 273
 274        IF (RHO .EQ. 0.D0) THEN
 275           A20 = 0.D0
 276           RETURN
 277        END IF
 278
 279        A20 = SM*( -DELTAM*KM**3/(4.D0*PI*CKM*RRHOP15) * EKM &
 280                  -( DELTAM*KM/(PI*SRRHOP)*KKM  &
 281                   + SIRHO*LAMBDA0(PSIM,KM) )/RHO**2 )  &
 282            - S0*( -DELTA0*K0**3/(4.D0*PI*CK0*RRHOP15) * EK0 &
 283                  -( DELTA0*K0/(PI*SRRHOP)*KK0  &
 284                   + SIRHO*LAMBDA0(PSI0,K0) )/RHO**2 )
 285
 286        A20 = A20/2.D0
 287
 288       RETURN
 289
 290   END FUNCTION FA20
 291
 292 END MODULE STRAIN_CALC