Tesis/ARRESTODINAMICO.F90 at main · gcarretas/Tesis · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
MODULE ARRESTODINAMICO
  USE VARIABLES
  IMPLICIT NONE
  CONTAINS

  SUBROUTINE CRITERIOARRESTO
    USE VARIABLES
    IMPLICIT NONE
    REAL(8)                    :: K, INICIO, TOL, LIMITE
    INTEGER :: I,J,NK
    INTEGER, DIMENSION (NP) :: FLAG
    REAL(8), DIMENSION (NP) :: ERROR, SUMA
    REAL(8), DIMENSION (NP,NP,NKMAX) :: R2RH, CR2R, ML, MLI, CKNI
    REAL(8), DIMENSION (NP,NP) :: NORM, NORMI, R2R, R2RI, FSELF, FSELFI, PSIK, PSIKI
    REAL(8), DIMENSION (NP,NP) :: ALMACENA, GAMMAD, GAMMADSOL
    REAL(8), DIMENSION (NP,10000) :: GAMMALIST

    R2R=0.0D0
    R2RI=0.0D0
    NORM=0.0D0
    DO I=1,NP
      DO J=1,NP
        R2R(I,I)=SQRT(RHO(I))
        R2RI(I,I)=1.0D0/R2R(I,I)
        NORM(I,J)=SQRT(RHO(I)*RHO(J))
      END DO
    END DO

    CALL INVERSION (NORM(:,:),NORMI(:,:),NP)

    OPEN(10, file='CKN.dat')
    OPEN(11, file='HKN.dat')

    DO NK=1,NKMAX
      K=KV1(NK)

      HK(:,:,NK)=MATMUL(NORMI(:,:),(SK(:,:,NK)-KRONEKER))
      WRITE(11,*) SNGL(K), SNGL(HK(1,1,NK)), SNGL(HK(1,2,NK)), SNGL(HK(2,1,NK)), SNGL(HK(2,2,NK))

      CKNI(:,:,NK)=KRONEKER-MATMUL(NORM(:,:),SK(:,:,NK))
      CALL INVERSION(CKNI(:,:,NK),CKN(:,:,NK),NP)
      WRITE(10,*) SNGL(K), SNGL(CKN(1,1,NK)), SNGL(CKN(1,2,NK)), SNGL(CKN(2,1,NK)), SNGL(CKN(2,2,NK))

      R2RH(:,:,NK)=MATMUL(HK(:,:,NK),R2RI)
      CR2R(:,:,NK)=MATMUL(R2RI,CKN(:,:,NK))

      CALL MATRIZLAMBDA(K,NK,ML)
      CALL INVERSION(ML(:,:,NK),MLI(:,:,NK),NP)
    END DO

    CLOSE(10)
    CLOSE(11)

    GAMMALIST=0.0D0
    GAMMAD=0.0D0
    ERROR=1.0D0
    INICIO=1.D-6
    TOL=1.0E-6
    LIMITE=1.0E20

    !VALOR INICIAL DE GAMMA1 Y GAMMA2
    DO I=1,NP
      GAMMAD(I,I)=INICIO
    END DO

    I=0

    LOOP: DO I=1,3000
      SUMA=0.0D0
      DO NK=1,NKMAX
        K=KV1(NK)
        FSELF=KRONEKER+((K**2)*MATMUL(GAMMAD,MLI(:,:,NK)))
        CALL INVERSION(FSELF,FSELFI,NP)
        PSIK=KRONEKER+((K**2)*MATMUL(MATMUL(GAMMAD,MLI(:,:,NK)),SKI(:,:,NK)))
        CALL INVERSION(PSIK,PSIKI,NP)

        ALMACENA=MATMUL(MATMUL(CR2R(:,:,NK),PSIKI),R2RH(:,:,NK))

        DO J=1,NP
          SUMA(J)=SUMA(J)+((K**4)*FSELFI(J,J)*ALMACENA(J,J))
        END DO


      END DO

      SUMA(:)=(SUMA(:)*DK)/(6.0D0*(PI**2))

      DO J=1,NP
        IF (SUMA(J).LT.0.0D0) THEN
          PRINT*, SUMA(J), I
          STOP
        END IF
      END DO

      DO J=1,NP
        GAMMADSOL(J,J)=1.D0/SUMA(J)
        !IF (GAMMADSOL(J,J).GT.1.D20) THEN
        !  GAMMADSOL(J,J)=1.D20
        !END IF
        ERROR(J)=ABS(1.0D0-(GAMMAD(J,J)/GAMMADSOL(J,J)))
      END DO

      GAMMAD=GAMMADSOL

      PRINT *, I, SNGL(GAMMAD(1,1)), SNGL(GAMMAD(2,2))

      IF ((GAMMAD(1,1).GT.1D5).AND.(GAMMAD(2,2).GT.1D5)) EXIT LOOP
      IF ((TOL.GT.ERROR(1).AND.GAMMAD(1,1).LT.LIMITE).AND.(TOL.GT.ERROR(2).AND.GAMMAD(2,2).LT.LIMITE)) EXIT LOOP

    END DO LOOP

      !PRINT *, I, SNGL(GAMMAD(1,1))
    RETURN
  END SUBROUTINE CRITERIOARRESTO

  REAL FUNCTION LAMBDA(A)
    USE VARIABLES
    IMPLICIT NONE

    REAL(8) :: A

    LAMBDA=1.0D0/(1.0D0+(A/KMIN)**2)

  END FUNCTION LAMBDA

  SUBROUTINE MATRIZLAMBDA(K,NK,MLL)
    USE VARIABLES
    IMPLICIT NONE

    REAL(8)                    :: K
    REAL(8), DIMENSION (NP,NP,NKMAX) :: MLL
    INTEGER :: I,NK

    DO I=1,NP
      MLL(I,I,NK)=LAMBDA(K)
    END DO

  END SUBROUTINE MATRIZLAMBDA

END MODULE