input.f90

   1   MODULE Input_Data
   2   IMPLICIT NONE
   3
   4   REAL, SAVE  :: A_S,B_S,C_S ! -> Superlattice constants
   5   REAL, DIMENSION(3), SAVE :: A1_S,A2_S,A3_S ! -> Superlattice vectors
   6   INTEGER, SAVE :: NMin_X,NMax_X,NMin_Y,NMax_Y,NMin_Z,NMax_Z
   7   REAL, SAVE    :: X_Min,X_Max,Y_Min,Y_Max,Z_Min,Z_Max
   8   REAL, SAVE    :: X_Inc,Y_Inc,Z_Inc
   9   INTEGER, SAVE :: XDim,YDim,ZDim
  10   INTEGER, SAVE :: KCOOR
  11   INTEGER,PRIVATE :: NDots_X,NDots_Y,NDots_Z
  12
  13   !!! Material Parameters
  14   REAL, SAVE    :: XLAMB,XMU,C13,C33,C11,EPSC,EPSA
  15   REAL, PRIVATE :: PZE15,PZE31,PZE33,PSPB,PSPW,DIELEC
  16   REAL, SAVE    :: ACW,AVW,BW,DW
  17   REAL, SAVE    :: VWE,VWH,VWSO,DW1,DW2,DW3
  18   REAL, SAVE    :: VBE,VBH,VBSO,DB1,DB2,DB3
  19   REAL, SAVE    :: C1,C2,D1,D2
  20   REAL, SAVE    :: D3,D4,D5,D6
  21
  22   !!! Calculation Constants
  23   REAL, SAVE   :: ETA1,ETA2,ETA_DIF,CN02_1,CN02_2,CN42_1, &
  24                   CN42_2,CN31_1,CN31_2,CNP31_1,CNP31_2, &
  25                   CN1_2G,BISUM,BIZZ,BIAUX,CARAC_MAX,&
  26                   U1,U2
  27   REAL, SAVE    :: CI1,CI2,CI3
  28   REAL, SAVE    :: CNE2_1,CNE2_2,CNE_1,CNE_2,CSP,CPZ
  29   REAL, SAVE    :: CSPWL,CSPBR,CPZWL,CPZ_Q1,CPZ_Q2
  30
  31   !!! Normalized Dot Parameters
  32   REAL, SAVE :: RC,ZC,D,RD,HD
  33
  34   INTEGER,SAVE              :: STR_Action,PZO_Action,POT_Action, &
  35                                DIAG_Action,DPL_Action
  36   CHARACTER (LEN=120), SAVE :: STR_Filename,PZO_Filename,POT_Filename,&
  37                                DPL_Filename
  38
  39   CONTAINS
  40
  41   SUBROUTINE READ_INPUT()
  42
  43       Use Dot_Geometry
  44       Use Auxiliar_Procedures, ONLY : AISO, MTYPE
  45
  46   IMPLICIT NONE
  47
  48       INTEGER :: G_Method
  49       REAL    :: D_Z
  50       REAL,DIMENSION(3) :: V_AUX
  51
  52
  53
  54 !!!!! read data from data file !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  55
  56       READ (*,*)
  57       READ (*,*)
  58       READ (*,*)
  59       READ (*,*)
  60 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  61       READ (*,*)
  62       READ (*,*) DWL            ! Wetting layer thickness (A)
  63       READ (*,*) ISHAPE         ! Shape of the dot
  64       READ (*,*) HQD            ! Height of the quantum dot (A)
  65       READ (*,*) Rqd_Base       ! Base Radius (A)
  66       READ (*,*) Rqd_Top        ! Top Radius (A)
  67 !!! Superlattice Begins !!!!
  68       READ (*,*)
  69       READ (*,*) A_S,B_S,C_S
  70       READ (*,*) A1_S(1),A1_S(2),A1_S(3)
  71       READ (*,*) A2_S(1),A2_S(2),A2_S(3)
  72       READ (*,*) A3_S(1),A3_S(2),A3_S(3)
  73       READ (*,*) NDots_X,NDots_Y,NDots_Z
  74       V_AUX=(/A_S,B_S,C_S/)
  75       A1_S=A1_S*V_AUX
  76       A2_S=A2_S*V_AUX
  77       A3_S=A3_S*V_AUX
  78       IF(MOD(NDots_X,2).EQ.0) THEN
  79               NDots_X=NDots_X+1
  80               WRITE(6,'(A,1X,I3)')"Warning: NDots_X Even -> NDots_X=",NDots_X
  81       END IF
  82       IF(MOD(NDots_Y,2).EQ.0) THEN
  83               NDots_Y=NDots_Y+1
  84               WRITE(6,'(A,1X,I3)')"Warning: NDots_Y Even -> NDots_Y=",NDots_Y
  85       END IF
  86       IF(MOD(NDots_Z,2).EQ.0) THEN
  87               NDots_Z=NDots_Z+1
  88               WRITE(6,'(A,1X,I3)')"Warning: NDots_Z Even -> NDots_Z=",NDots_Z
  89       END IF
  90       IF(NDots_X.EQ.1) THEN
  91               NMin_X=0; NMax_X=0
  92       ELSE
  93               NMin_X=-(NDots_X-1)/2; NMax_X=-NMin_X
  94       END IF
  95       IF(NDots_Y.EQ.1) THEN
  96               NMin_Y=0; NMax_Y=0
  97       ELSE
  98               NMin_Y=-(NDots_Y-1)/2; NMax_Y=-NMin_Y
  99       END IF
 100       IF(NDots_Z.EQ.1) THEN
 101               NMin_Z=0; NMax_Z=0
 102       ELSE
 103               NMin_Z=-(NDots_Z-1)/2; NMax_Z=-NMin_Z
 104       END IF
 105 !!! Grid Begins !!!!
 106       READ (*,*)
 107       READ (*,*) G_Method
 108       READ (*,*) X_Min,X_Max
 109       READ (*,*) Y_Min,Y_Max
 110       READ (*,*) Z_Min,Z_Max
 111       READ (*,*) XDim,YDim,ZDim
 112
 113       IF(G_Method.EQ.1) THEN
 114            IF(XDim.EQ.1) THEN
 115               X_Max=X_Min
 116            ELSE
 117               X_Max=A1_S(1)/2.E0
 118               X_Min=0.E0
 119            END IF
 120            IF(YDim.EQ.1) THEN
 121               Y_Max=Y_Min
 122            ELSE
 123               Y_Max=A2_S(2)/2.E0
 124               Y_Min=0.E0
 125            END IF
 126            IF(ZDim.EQ.1) THEN
 127               Z_Max=Z_Min
 128            ELSE
 129               Z_Max=A3_S(3)/2.E0
 130               D_Z=Z_Max-(Hqd+DWL)/2.E0
 131               Z_Max=Hqd+D_Z
 132               Z_Min=-DWL-D_Z
 133            END IF
 134       END IF
 135
 136       IF (XDim.EQ.1) THEN
 137               X_Inc=0.E0
 138       ELSE
 139               X_Inc=(X_Max-X_Min)/REAL(XDim-1)
 140       END IF
 141       IF (YDim.EQ.1) THEN
 142               Y_Inc=0.E0
 143       ELSE
 144               Y_Inc=(Y_Max-Y_Min)/REAL(YDim-1)
 145       END IF
 146       IF (ZDim.EQ.1) THEN
 147               Z_Inc=0.E0
 148       ELSE
 149               Z_Inc=(Z_Max-Z_Min)/REAL(ZDim-1)
 150       END IF
 151
 152 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 153
 154       READ (*,*)
 155       READ (*,*) MTYPE          ! MTYPE: 1-> ZBI, 2-> WZI, 3-> WZA
 156       SELECT CASE (MTYPE)
 157        CASE(1)
 158          READ(*,*)VWE,VWH,VWSO
 159          READ(*,*)VBE,VBH,VBSO
 160          READ(*,*)ACW,AVW,BW,DW
 161          READ(*,*)XLAMB,XMU
 162          READ(*,*)EPSA
 163          READ(*,*)
 164         CASE(2:)
 165          READ(*,*)VWE,VWH,DW1,DW2,DW3
 166          READ(*,*)VBE,VBH,DB1,DB2,DB3
 167          READ(*,*)C1,C2,D1,D2
 168          READ(*,*)D3,D4,D5,D6
 169          READ(*,*)C11,C13,C33    ! Elastic modulus
 170          READ(*,*)XLAMB,XMU      ! Lame constants
 171          READ(*,*)EPSA           ! Misfit strain: EPS0
 172          READ(*,*)EPSC           ! Misfit strain: EPSC
 173          READ(*,*)PZE15,PZE31,PZE33 !Piezoelectric Constants
 174          READ(*,*)PSPB,PSPW      ! Spontaneous Polarization
 175          READ(*,*)DIELEC         ! Barrier Dielectric Constant
 176          READ(*,*)
 177       END SELECT
 178
 179       READ(*,*)CARAC_MAX
 180       READ(*,*)STR_Action     ! 0-> Nothing 1-> Caculate 2-> Read
 181       READ(*,*)DPL_Action     ! 0-> Nothing 1-> Caculate 2-> Read
 182       READ(*,*)PZO_Action     ! 0-> Nothing 1-> Caculate 2-> Read
 183       READ(*,*)POT_Action     ! 0-> Nothing 1-> Caculate
 184       READ(*,*)DIAG_Action    ! 0-> Only Diagonal, 1->Full Hamiltonian
 185       READ(*,*)
 186       IF(PZO_Action.GT.0.AND.MTYPE.EQ.1) THEN
 187         WRITE(6,*)"Error: Electrostatic Potential for ZB not implemented"
 188         STOP
 189       END IF
 190       IF (MTYPE.LT.3) THEN
 191        AISO=1 ! Isotropic Material
 192       ELSE
 193        AISO=0 ! Anisotropic Material
 194       END IF
 195       READ(*,*) STR_Filename
 196       READ(*,*) DPL_Filename
 197       READ(*,*) PZO_Filename
 198       READ(*,*) POT_Filename
 199       READ(*,*) KCOOR          ! KCOOR=0 Cartesian Coordinates
 200                                ! KCOOR=1 Cylindrical Coordinates
 201       IF(KCOOR.EQ.1.AND.YDIM.GT.1) THEN
 202          WRITE(6,*)"ERROR: Selected coordinates are cylindrical and YDIM.GT.1"
 203          WRITE(6,*)"       Exiting Program"
 204          STOP
 205       END IF
 206       READ (*,*)
 207 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 208
 209       RETURN
 210
 211     END SUBROUTINE READ_INPUT
 212
 213     SUBROUTINE CONSTANTS( )
 214         Use Dot_Geometry, ONLY: DWL
 215         Use Auxiliar_Procedures, ONLY : AISO, MTYPE
 216       IMPLICIT NONE
 217
 218       REAL :: ALPHA,BETA,GAMA
 219       REAL :: R,S,A1,A2,B1,B2
 220       REAL :: CN1_1,CN1_2,CN3_1,CN3_2,CN4_1,CN4_2
 221       REAL :: P1,P2,P3
 222       REAL :: ETA_AUX, PI
 223
 224       PI=4.E0*ATAN(1.E0)
 225
 226 !!!! Remember: XLAMB=C12, XMU=C44 !!!!!
 227
 228       !!! ISO CONSTANTS
 229       BISUM =  2.*EPSA          ! Biaxial strain components
 230       BIAUX =  EPSA-BIZZ
 231       BIZZ  = -2.*EPSA*XLAMB/(XLAMB+2.*XMU)
 232       IF(MTYPE.EQ.3) THEN
 233         IF(STR_Action.EQ.0) THEN
 234            BIZZ = -2.E0*EPSA*C13/C33
 235         ELSE
 236            BIZZ=EPSC
 237         END IF
 238       END IF
 239
 240       IF(MTYPE.EQ.1) RETURN
 241
 242       R=(C11+XLAMB+C13)*EPSA+C13*(EPSC-EPSA)
 243       S=(C33+C13-XLAMB-C11)*EPSA+(C33-C13)*(EPSC-EPSA)
 244
 245       ALPHA=-XMU*C33
 246
 247       R=R/ALPHA; S=S/ALPHA
 248
 249       A1=R*(C33-XMU-C13)-S*(C13+XMU)
 250       A2=R*XMU
 251
 252       B1=S*(C13+2.E0*XMU)-R*(C33-C13-2.E0*XMU)
 253       B2=S*C11-R*(C13+2.E0*XMU-C11)
 254
 255       BETA=( C13*(C13+2.E0*XMU)-C11*C33 )/(-XMU*C33)
 256       GAMA=C11/C33
 257
 258       CN1_1=A1;         CN1_2=A2
 259       CN3_1=A1+B1;      CN3_2=A2+B2
 260       CN4_1=A1+B1/2.E0; CN4_2=A2+B2/2.E0
 261
 262       P1=2.E0*PZE15; P2=PZE31; P3=PZE33
 263
 264       CI1=2.E0*EPSA*P2+EPSC*P3
 265       CI2=P3
 266       CI3=1.E0
 267
 268       CNE_1=(P1+P2-P3)*A1+(P1/2-P3)*B1
 269       CNE_2=(P1+P2-P3)*A2+(P1/2-P3)*B2
 270
 271       ETA_AUX=SQRT(BETA**2 - 4.E0*GAMA)
 272       ETA1=(-BETA+ETA_AUX)/2.E0
 273       ETA2=(-BETA-ETA_AUX)/2.E0
 274       IF(ETA1.GT.0.E0.OR.ETA2.GT.0.E0) THEN
 275         WRITE(6,*)"ERROR: ETA1 AND ETA2 ARE GREATHER THAN ZERO"
 276         STOP
 277       END IF
 278       ETA1=ABS(ETA1)
 279       ETA2=ABS(ETA2)
 280
 281       ETA_DIF=ETA1-ETA2
 282
 283       CN02_1=(CN1_2/ETA2-CN1_1)
 284       CN02_2=(CN1_2/ETA1-CN1_1)
 285
 286       CN31_1=(CN4_1*SQRT(ETA2)-CN4_2/SQRT(ETA2))
 287       CN31_2=(CN4_1*SQRT(ETA1)-CN4_2/SQRT(ETA1))
 288
 289       CNP31_1=(CN3_2/SQRT(ETA2)-CN3_1*SQRT(ETA2))
 290       CNP31_2=(CN3_2/SQRT(ETA1)-CN3_1*SQRT(ETA1))
 291
 292       CN42_1=(CN3_1*ETA2-CN3_2)
 293       CN42_2=(CN3_1*ETA1-CN3_2)
 294
 295       U1=CN3_2/SQRT(ETA2)-CN3_1*SQRT(ETA2)
 296       U2=CN3_2/SQRT(ETA1)-CN3_1*SQRT(ETA1)
 297
 298       CN1_2G=CN1_2/GAMA
 299
 300       CNE2_1=(1.E0-ETA1)*(CNE_2/SQRT(ETA2)-CNE_1*SQRT(ETA2))
 301       CNE2_2=(1.E0-ETA2)*(CNE_2/SQRT(ETA1)-CNE_1*SQRT(ETA1))
 302
 303       CSP=(PSPW-PSPB)
 304       !!! V=q/(4*PI*E0)*PZ=0.9*PZ eV
 305       CPZ=0.9*(4.E0*PI/DIELEC)
 306       CPZ_Q1=(BISUM*PZE31+BIZZ*PZE33)
 307       CPZ_Q2=(BIAUX*(PZE33-(PZE31+2.E0*PZE15)))
 308
 309 !!!   Wetting Layer Parameters
 310
 311       IF(NDots_z.GT.1) THEN
 312         CSPWL=CSP*(C_S-DWL)/C_S
 313         CSPBR=-CSP*(DWL/C_S)
 314       ELSE
 315         CSPWL=CSP
 316         CSPBR=0.E0
 317       END IF
 318
 319         IF(AISO.EQ.0) THEN
 320            CPZWL=(PZE31-PZE33*(C13/C33))*BISUM
 321         ELSE
 322            CPZWL=(BISUM*PZE31+BIZZ*PZE33)
 323         END IF
 324
 325       RETURN
 326
 327     END SUBROUTINE CONSTANTS
 328
 329     SUBROUTINE CONSTANTS_PZO( )
 330         Use Dot_Geometry, ONLY: DWL
 331         Use Auxiliar_Procedures, ONLY : AISO
 332       IMPLICIT NONE
 333
 334       REAL :: ALPHA,BETA,GAMA
 335       REAL :: R,S,A1,A2,B1,B2
 336       REAL :: CN3_1,CN3_2,P1,P2,P3
 337       REAL :: ETA_AUX, PI
 338
 339       PI=4.E0*ATAN(1.E0)
 340
 341 !!!! Remember: XLAMB=C12, XMU=C44 !!!!!
 342
 343       !!! ISO CONSTANTS
 344       BISUM =  2.*EPSA          ! Biaxial strain components
 345       BIZZ  = -2.*EPSA*XLAMB/(XLAMB+2.*XMU)
 346       BIAUX =  EPSA-BIZZ
 347
 348       R=(C11+XLAMB+C13)*EPSA+C13*(EPSC-EPSA)
 349       S=(C33+C13-XLAMB-C11)*EPSA+(C33-C13)*(EPSC-EPSA)
 350
 351       ALPHA=-XMU*C33
 352
 353       R=R/ALPHA; S=S/ALPHA
 354
 355       A1=R*(C33-XMU-C13)-S*(C13+XMU)
 356       A2=R*XMU
 357
 358       B1=S*(C13+2.E0*XMU)-R*(C33-C13-2.E0*XMU)
 359       B2=S*C11-R*(C13+2.E0*XMU-C11)
 360
 361       BETA=( C13*(C13+2.E0*XMU)-C11*C33 )/(-XMU*C33)
 362       GAMA=C11/C33
 363
 364       CN3_1=A1+B1;      CN3_2=A2+B2
 365
 366       P1=2.E0*PZE15; P2=PZE31; P3=PZE33
 367
 368       CI1=2.E0*EPSA*P2+EPSC*P3
 369       CI2=P3
 370       CI3=1.E0
 371
 372       CNE_1=(P1+P2-P3)*A1+(P1/2-P3)*B1
 373       CNE_2=(P1+P2-P3)*A2+(P1/2-P3)*B2
 374
 375       ETA_AUX=SQRT(BETA**2 - 4.E0*GAMA)
 376       ETA1=(-BETA+ETA_AUX)/2.E0
 377       ETA2=(-BETA-ETA_AUX)/2.E0
 378       IF(ETA1.GT.0.E0.OR.ETA2.GT.0.E0) THEN
 379         WRITE(6,*)"ERROR: ETA1 AND ETA2 ARE GREATHER THAN ZERO"
 380         STOP
 381       END IF
 382       ETA1=ABS(ETA1)
 383       ETA2=ABS(ETA2)
 384
 385       ETA_DIF=ETA1-ETA2
 386
 387       CNP31_1=(CN3_2/SQRT(ETA2)-CN3_1*SQRT(ETA2))
 388       CNP31_2=(CN3_2/SQRT(ETA1)-CN3_1*SQRT(ETA1))
 389
 390       CNE2_1=(1.E0-ETA1)*(CNE_2/SQRT(ETA2)-CNE_1*SQRT(ETA2))
 391       CNE2_2=(1.E0-ETA2)*(CNE_2/SQRT(ETA1)-CNE_1*SQRT(ETA1))
 392
 393       CSP=(PSPW-PSPB)
 394       !!! V=q/(4*PI*E0)*PZ=0.9*PZ eV
 395       CPZ=0.9*(4.E0*PI/DIELEC)
 396       CPZ_Q1=(BISUM*PZE31+BIZZ*PZE33)
 397       CPZ_Q2=(BIAUX*(PZE33-(PZE31+2.E0*PZE15)))
 398
 399 !!!   Wetting Layer Parameters
 400
 401       IF(NDots_z.GT.1) THEN
 402         CSPWL=CSP*(C_S-DWL)/C_S
 403         CSPBR=-CSP*(DWL/C_S)
 404       ELSE
 405         CSPWL=CSP
 406         CSPBR=0.E0
 407       END IF
 408
 409         IF(AISO.EQ.0) THEN
 410            CPZWL=(PZE31-PZE33*(C13/C33))*BISUM
 411         ELSE
 412            CPZWL=(BISUM*PZE31+BIZZ*PZE33)
 413         END IF
 414
 415       RETURN
 416
 417     END SUBROUTINE CONSTANTS_PZO
 418
 419     SUBROUTINE CONSTANTS_STR( )
 420
 421       Use Auxiliar_Procedures, ONLY : MTYPE
 422
 423       IMPLICIT NONE
 424
 425       REAL :: ALPHA,BETA,GAMA
 426       REAL :: R,S,A1,A2,B1,B2
 427       REAL :: CN1_1,CN1_2,CN3_1,CN3_2,CN4_1,CN4_2
 428       REAL :: ETA_AUX
 429
 430 !!!! Remember: XLAMB=C12, XMU=C44 !!!!!
 431
 432       !!! ISO CONSTANTS
 433       BISUM =  2.*EPSA          ! Biaxial strain components
 434       BIZZ  = -2.*EPSA*XLAMB/(XLAMB+2.*XMU)
 435       BIAUX =  EPSA-BIZZ
 436       IF(MTYPE.EQ.3) THEN
 437         IF(STR_Action.EQ.0) THEN
 438            BIZZ = -2.E0*EPSA*C13/C33
 439         ELSE
 440            BIZZ=EPSC
 441         END IF
 442       END IF
 443
 444       IF(MTYPE.EQ.1) RETURN
 445
 446       R=(C11+XLAMB+C13)*EPSA+C13*(EPSC-EPSA)
 447       S=(C33+C13-XLAMB-C11)*EPSA+(C33-C13)*(EPSC-EPSA)
 448
 449       ALPHA=-XMU*C33
 450
 451       R=R/ALPHA; S=S/ALPHA
 452
 453       A1=R*(C33-XMU-C13)-S*(C13+XMU)
 454       A2=R*XMU
 455
 456       B1=S*(C13+2.E0*XMU)-R*(C33-C13-2.E0*XMU)
 457       B2=S*C11-R*(C13+2.E0*XMU-C11)
 458
 459       BETA=( C13*(C13+2.E0*XMU)-C11*C33 )/(-XMU*C33)
 460       GAMA=C11/C33
 461
 462       CN1_1=A1;         CN1_2=A2
 463       CN3_1=A1+B1;      CN3_2=A2+B2
 464       CN4_1=A1+B1/2.E0; CN4_2=A2+B2/2.E0
 465
 466       ETA_AUX=SQRT(BETA**2 - 4.E0*GAMA)
 467       ETA1=(-BETA+ETA_AUX)/2.E0
 468       ETA2=(-BETA-ETA_AUX)/2.E0
 469       IF(ETA1.GT.0.E0.OR.ETA2.GT.0.E0) THEN
 470         WRITE(6,*)"ERROR: ETA1 AND ETA2 ARE GREATHER THAN ZERO"
 471         STOP
 472       END IF
 473       ETA1=ABS(ETA1)
 474       ETA2=ABS(ETA2)
 475
 476       ETA_DIF=ETA1-ETA2
 477
 478       CN02_1=(CN1_2/ETA2-CN1_1)
 479       CN02_2=(CN1_2/ETA1-CN1_1)
 480
 481       CN31_1=(CN4_1*SQRT(ETA2)-CN4_2/SQRT(ETA2))
 482       CN31_2=(CN4_1*SQRT(ETA1)-CN4_2/SQRT(ETA1))
 483
 484       CN42_1=(CN3_1*ETA2-CN3_2)
 485       CN42_2=(CN3_1*ETA1-CN3_2)
 486
 487       CN1_2G=CN1_2/GAMA
 488
 489       RETURN
 490
 491     END SUBROUTINE CONSTANTS_STR
 492
 493     SUBROUTINE CONSTANTS_DSP( )
 494
 495       Use Auxiliar_Procedures, ONLY : MTYPE
 496
 497       IMPLICIT NONE
 498
 499       REAL :: ALPHA,BETA,GAMA
 500       REAL :: R,S,A1,A2,B1,B2
 501       REAL :: CN1_1,CN1_2,CN3_1,CN3_2
 502       REAL :: ETA_AUX
 503
 504 !!!! Remember: XLAMB=C12, XMU=C44 !!!!!
 505
 506       !!! ISO CONSTANTS
 507       BISUM =  2.*EPSA          ! Biaxial strain components
 508       BIZZ  = -2.*EPSA*XLAMB/(XLAMB+2.*XMU)
 509       BIAUX =  EPSA-BIZZ
 510       IF(MTYPE.EQ.3) THEN
 511         IF(STR_Action.EQ.0) THEN
 512            BIZZ = -2.E0*EPSA*C13/C33
 513         ELSE
 514            BIZZ=EPSC
 515         END IF
 516       END IF
 517
 518       IF(MTYPE.EQ.1) RETURN
 519
 520       R=(C11+XLAMB+C13)*EPSA+C13*(EPSC-EPSA)
 521       S=(C33+C13-XLAMB-C11)*EPSA+(C33-C13)*(EPSC-EPSA)
 522
 523       ALPHA=-XMU*C33
 524
 525       R=R/ALPHA; S=S/ALPHA
 526
 527       A1=R*(C33-XMU-C13)-S*(C13+XMU)
 528       A2=R*XMU
 529
 530       B1=S*(C13+2.E0*XMU)-R*(C33-C13-2.E0*XMU)
 531       B2=S*C11-R*(C13+2.E0*XMU-C11)
 532
 533       BETA=( C13*(C13+2.E0*XMU)-C11*C33 )/(-XMU*C33)
 534       GAMA=C11/C33
 535
 536       CN1_1=A1;         CN1_2=A2
 537       CN3_1=A1+B1;      CN3_2=A2+B2
 538
 539       ETA_AUX=SQRT(BETA**2 - 4.E0*GAMA)
 540       ETA1=(-BETA+ETA_AUX)/2.E0
 541       ETA2=(-BETA-ETA_AUX)/2.E0
 542       IF(ETA1.GT.0.E0.OR.ETA2.GT.0.E0) THEN
 543         WRITE(6,*)"ERROR: ETA1 AND ETA2 ARE GREATHER THAN ZERO"
 544         STOP
 545       END IF
 546       ETA1=ABS(ETA1)
 547       ETA2=ABS(ETA2)
 548
 549       ETA_DIF=ETA1-ETA2
 550
 551       CN02_1=(CN1_2/ETA2-CN1_1)
 552       CN02_2=(CN1_2/ETA1-CN1_1)
 553
 554       U1=CN3_2/SQRT(ETA2)-CN3_1*SQRT(ETA2)
 555       U2=CN3_2/SQRT(ETA1)-CN3_1*SQRT(ETA1)
 556
 557       CN1_2G=CN1_2/GAMA
 558
 559       RETURN
 560
 561     END SUBROUTINE CONSTANTS_DSP
 562
 563   END MODULE INPUT_DATA
 564