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