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
 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       PI=4.E0*ATAN(1.E0)
 344
 345 !!!! Remember: XLAMB=C12, XMU=C44 !!!!!
 346
 347       !!! ISO CONSTANTS
 348       BISUM =  2.*EPSA          ! Biaxial strain components
 349       BIZZ  = -2.*EPSA*XLAMB/(XLAMB+2.*XMU)
 350       BIAUX =  EPSA-BIZZ
 351
 352       R=(C11+XLAMB+C13)*EPSA+C13*(EPSC-EPSA)
 353       S=(C33+C13-XLAMB-C11)*EPSA+(C33-C13)*(EPSC-EPSA)
 354
 355       ALPHA=-XMU*C33
 356
 357       R=R/ALPHA; S=S/ALPHA
 358
 359       A1=R*(C33-XMU-C13)-S*(C13+XMU)
 360       A2=R*XMU
 361
 362       B1=S*(C13+2.E0*XMU)-R*(C33-C13-2.E0*XMU)
 363       B2=S*C11-R*(C13+2.E0*XMU-C11)
 364
 365       BETA=( C13*(C13+2.E0*XMU)-C11*C33 )/(-XMU*C33)
 366       GAMA=C11/C33
 367
 368       CN3_1=A1+B1;      CN3_2=A2+B2
 369
 370       P1=2.E0*PZE15; P2=PZE31; P3=PZE33
 371
 372       CI1=2.E0*EPSA*P2+EPSC*P3
 373       CI2=P3
 374       CI3=1.E0
 375
 376       CNE_1=(P1+P2-P3)*A1+(P1/2-P3)*B1
 377       CNE_2=(P1+P2-P3)*A2+(P1/2-P3)*B2
 378
 379       ETA_AUX=SQRT(BETA**2 - 4.E0*GAMA)
 380       ETA1=(-BETA+ETA_AUX)/2.E0
 381       ETA2=(-BETA-ETA_AUX)/2.E0
 382       IF(ETA1.GT.0.E0.OR.ETA2.GT.0.E0) THEN
 383         WRITE(6,*)"ERROR: ETA1 AND ETA2 ARE GREATHER THAN ZERO"
 384         STOP
 385       END IF
 386       ETA1=ABS(ETA1)
 387       ETA2=ABS(ETA2)
 388
 389       ETA_DIF=ETA1-ETA2
 390
 391       CNP31_1=(CN3_2/SQRT(ETA2)-CN3_1*SQRT(ETA2))
 392       CNP31_2=(CN3_2/SQRT(ETA1)-CN3_1*SQRT(ETA1))
 393
 394       CNE2_1=(1.E0-ETA1)*(CNE_2/SQRT(ETA2)-CNE_1*SQRT(ETA2))
 395       CNE2_2=(1.E0-ETA2)*(CNE_2/SQRT(ETA1)-CNE_1*SQRT(ETA1))
 396
 397       CSP=(PSPW-PSPB)
 398       !!! V=q/(4*PI*E0)*PZ=0.9*PZ eV
 399       CPZ=0.9*(4.E0*PI/DIELEC)
 400       CPZ_Q1=(BISUM*PZE31+BIZZ*PZE33)
 401       CPZ_Q2=(BIAUX*(PZE33-(PZE31+2.E0*PZE15)))
 402
 403 !!!   Wetting Layer Parameters
 404
 405       IF(NDots_z.GT.1) THEN
 406         CSPWL=CSP*(C_S-DWL)/C_S
 407         CSPBR=-CSP*(DWL/C_S)
 408       ELSE
 409         CSPWL=CSP
 410         CSPBR=0.E0
 411       END IF
 412
 413         IF(AISO.EQ.0) THEN
 414            CPZWL=(PZE31-PZE33*(C13/C33))*BISUM
 415         ELSE
 416            CPZWL=(BISUM*PZE31+BIZZ*PZE33)
 417         END IF
 418
 419       RETURN
 420
 421     END SUBROUTINE CONSTANTS_PZO
 422
 423     SUBROUTINE CONSTANTS_STR( )
 424
 425       Use Auxiliar_Procedures, ONLY : MTYPE
 426
 427       IMPLICIT NONE
 428
 429       REAL :: ALPHA,BETA,GAMA
 430       REAL :: R,S,A1,A2,B1,B2
 431       REAL :: CN1_1,CN1_2,CN3_1,CN3_2,CN4_1,CN4_2
 432       REAL :: ETA_AUX
 433
 434 !!!! Remember: XLAMB=C12, XMU=C44 !!!!!
 435
 436       !!! ISO CONSTANTS
 437       BISUM =  2.*EPSA          ! Biaxial strain components
 438       BIZZ  = -2.*EPSA*XLAMB/(XLAMB+2.*XMU)
 439       BIAUX =  EPSA-BIZZ
 440       IF(MTYPE.EQ.3) THEN
 441         IF(STR_Action.EQ.0) THEN
 442            BIZZ = -2.E0*EPSA*C13/C33
 443         ELSE
 444            BIZZ=EPSC
 445         END IF
 446       END IF
 447
 448       IF(MTYPE.EQ.1) RETURN
 449
 450       R=(C11+XLAMB+C13)*EPSA+C13*(EPSC-EPSA)
 451       S=(C33+C13-XLAMB-C11)*EPSA+(C33-C13)*(EPSC-EPSA)
 452
 453       ALPHA=-XMU*C33
 454
 455       R=R/ALPHA; S=S/ALPHA
 456
 457       A1=R*(C33-XMU-C13)-S*(C13+XMU)
 458       A2=R*XMU
 459
 460       B1=S*(C13+2.E0*XMU)-R*(C33-C13-2.E0*XMU)
 461       B2=S*C11-R*(C13+2.E0*XMU-C11)
 462
 463       BETA=( C13*(C13+2.E0*XMU)-C11*C33 )/(-XMU*C33)
 464       GAMA=C11/C33
 465
 466       CN1_1=A1;         CN1_2=A2
 467       CN3_1=A1+B1;      CN3_2=A2+B2
 468       CN4_1=A1+B1/2.E0; CN4_2=A2+B2/2.E0
 469
 470       ETA_AUX=SQRT(BETA**2 - 4.E0*GAMA)
 471       ETA1=(-BETA+ETA_AUX)/2.E0
 472       ETA2=(-BETA-ETA_AUX)/2.E0
 473       IF(ETA1.GT.0.E0.OR.ETA2.GT.0.E0) THEN
 474         WRITE(6,*)"ERROR: ETA1 AND ETA2 ARE GREATHER THAN ZERO"
 475         STOP
 476       END IF
 477       ETA1=ABS(ETA1)
 478       ETA2=ABS(ETA2)
 479
 480       ETA_DIF=ETA1-ETA2
 481
 482       CN02_1=(CN1_2/ETA2-CN1_1)
 483       CN02_2=(CN1_2/ETA1-CN1_1)
 484
 485       CN31_1=(CN4_1*SQRT(ETA2)-CN4_2/SQRT(ETA2))
 486       CN31_2=(CN4_1*SQRT(ETA1)-CN4_2/SQRT(ETA1))
 487
 488       CN42_1=(CN3_1*ETA2-CN3_2)
 489       CN42_2=(CN3_1*ETA1-CN3_2)
 490
 491       CN1_2G=CN1_2/GAMA
 492
 493       RETURN
 494
 495     END SUBROUTINE CONSTANTS_STR
 496
 497     SUBROUTINE CONSTANTS_DSP( )
 498
 499       Use Auxiliar_Procedures, ONLY : MTYPE
 500
 501       IMPLICIT NONE
 502
 503       REAL :: ALPHA,BETA,GAMA
 504       REAL :: R,S,A1,A2,B1,B2
 505       REAL :: CN1_1,CN1_2,CN3_1,CN3_2
 506       REAL :: ETA_AUX
 507
 508 !!!! Remember: XLAMB=C12, XMU=C44 !!!!!
 509
 510       !!! ISO CONSTANTS
 511       BISUM =  2.*EPSA          ! Biaxial strain components
 512       BIZZ  = -2.*EPSA*XLAMB/(XLAMB+2.*XMU)
 513       BIAUX =  EPSA-BIZZ
 514       IF(MTYPE.EQ.3) THEN
 515         IF(STR_Action.EQ.0) THEN
 516            BIZZ = -2.E0*EPSA*C13/C33
 517         ELSE
 518            BIZZ=EPSC
 519         END IF
 520       END IF
 521
 522       IF(MTYPE.EQ.1) RETURN
 523
 524       R=(C11+XLAMB+C13)*EPSA+C13*(EPSC-EPSA)
 525       S=(C33+C13-XLAMB-C11)*EPSA+(C33-C13)*(EPSC-EPSA)
 526
 527       ALPHA=-XMU*C33
 528
 529       R=R/ALPHA; S=S/ALPHA
 530
 531       A1=R*(C33-XMU-C13)-S*(C13+XMU)
 532       A2=R*XMU
 533
 534       B1=S*(C13+2.E0*XMU)-R*(C33-C13-2.E0*XMU)
 535       B2=S*C11-R*(C13+2.E0*XMU-C11)
 536
 537       BETA=( C13*(C13+2.E0*XMU)-C11*C33 )/(-XMU*C33)
 538       GAMA=C11/C33
 539
 540       CN1_1=A1;         CN1_2=A2
 541       CN3_1=A1+B1;      CN3_2=A2+B2
 542
 543       ETA_AUX=SQRT(BETA**2 - 4.E0*GAMA)
 544       ETA1=(-BETA+ETA_AUX)/2.E0
 545       ETA2=(-BETA-ETA_AUX)/2.E0
 546       IF(ETA1.GT.0.E0.OR.ETA2.GT.0.E0) THEN
 547         WRITE(6,*)"ERROR: ETA1 AND ETA2 ARE GREATHER THAN ZERO"
 548         STOP
 549       END IF
 550       ETA1=ABS(ETA1)
 551       ETA2=ABS(ETA2)
 552
 553       ETA_DIF=ETA1-ETA2
 554
 555       CN02_1=(CN1_2/ETA2-CN1_1)
 556       CN02_2=(CN1_2/ETA1-CN1_1)
 557
 558       U1=CN3_2/SQRT(ETA2)-CN3_1*SQRT(ETA2)
 559       U2=CN3_2/SQRT(ETA1)-CN3_1*SQRT(ETA1)
 560
 561       CN1_2G=CN1_2/GAMA
 562
 563       RETURN
 564
 565     END SUBROUTINE CONSTANTS_DSP
 566
 567   END MODULE INPUT_DATA
 568