-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathModule_Conservation.f90
175 lines (141 loc) · 5.15 KB
/
Module_Conservation.f90
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
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
MODULE MODULE_CONSERVATION
CONTAINS
SUBROUTINE MASS_CONSERVATION(H,Dt,dist,ray,k,BV_a,BV_b,V_t1,V_t2,delta0)
!*****************************************************************
! Control the mass conservation at each time step pour le code balmforth
!*****************************************************************
IMPLICIT NONE
! Tableaux
DOUBLE PRECISION, DIMENSION(:,:) ,INTENT(IN) :: H
DOUBLE PRECISION, DIMENSION(:), INTENT(IN) :: dist,ray
! Parametre du model
DOUBLE PRECISION ,INTENT(INOUT) :: BV_a,BV_b
DOUBLE PRECISION ,INTENT(INOUT) :: V_t1,V_t2
INTEGER ,INTENT(IN) :: k
! Nombre sans dimension
DOUBLE PRECISION ,INTENT(IN) :: Dt,delta0
!Parametre du sous programme
DOUBLE PRECISION :: A_t1,Int_t1,A_t2,Int_t2,Dr
INTEGER :: i,N
N = COUNT(H(:,1)>delta0)
Dr = ray(1)
! temps t
A_t1 = H(1,1)
Int_t1 = A_t1*Dr**2
DO i=2,N-1
A_t1 = H(i,1)
Int_t1 = Int_t1+A_t1*(ray(i)**2-ray(i-1)**2)
ENDDO
A_t1 = H(N,1)
Int_t1 = Int_t1+A_t1*(dist(N)**2-ray(N-1)**2)
! temps t+dt
A_t2 = H(1,3)
Int_t2 = A_t2*Dr**2
DO i=2,N-1
A_t2 = H(i,3)
Int_t2 = Int_t2+A_t2*(ray(i)**2-ray(i-1)**2)
ENDDO
A_t2 = H(N,3)
Int_t2 = Int_t2+A_t2*(dist(N)**2-ray(N-1)**2)
! Volume de l'energie au temps t+Dt
BV_a = Int_t2
! Volume de l'intrusion au temp t + ce qu'on a rajoute
BV_b = Int_t1+ Dt
! Volume temps t
V_t1 = Int_t1
! Volume temps t+dt
V_t2 = Int_t2
END SUBROUTINE MASS_CONSERVATION
SUBROUTINE ENERGY_CONSERVATION(H,BL,T,Ts,Pe,Dt,dist,ray,k,psi,BE_a,BE_b,En_t1,En_t2,Phi_s,Phi_l,delta0,sigma&
&,Rheology)
!*****************************************************************
! Control the energy conservation at each time step pour le code balmforth
!*****************************************************************
IMPLICIT NONE
! Tableaux
DOUBLE PRECISION, DIMENSION(:,:) ,INTENT(IN) :: Ts,H,T,BL
DOUBLE PRECISION, DIMENSION(:), INTENT(IN) :: dist,ray
! Parametre du model
DOUBLE PRECISION ,INTENT(INOUT) :: BE_a,BE_b
DOUBLE PRECISION ,INTENT(INOUT) :: En_t1,En_t2,Phi_s,Phi_l
INTEGER ,INTENT(IN) :: k
! Nombre sans dimension
DOUBLE PRECISION ,INTENT(IN) :: Pe,Dt,psi,delta0,sigma
INTEGER ,INTENT(IN) :: Rheology
!Parametre du sous programme
DOUBLE PRECISION :: tbar_t1,tbar_t2,A_1_t1,A_1_t2,A_2,A_3,A_4
DOUBLE PRECISION :: Int_1_t1,Int_1_t2,Int_2,Int_3,Int_4
DOUBLE PRECISION :: E_ta,E_tb,Dr
INTEGER :: i,N
N=0
DO i=1,COUNT(H(:,3)>0D0),1
IF (H(i,3)-delta0>0.d0 .OR. N /= 0) THEN
CYCLE
ELSE
N = i
END IF
ENDDO
Dr = ray(1)
! Calcule au temps t
tbar_t1 = -(2*(T(1,1)-Ts(1,1))*BL(1,1))/(3.d0*H(1,1))+T(1,1)
! Calcule au temps t
tbar_t2 = -(2*(T(1,3)-Ts(1,3))*BL(1,3))/(3.d0*H(1,3))+T(1,3)
A_1_t1 = tbar_t1*H(1,1)
A_1_t2 = tbar_t2*H(1,3)
A_2 = (tbar_t2*H(1,3)-tbar_t1*H(1,1))/Dt-Ts(1,1)*(H(1,3)-H(1,1))/Dt
A_3 = ((T(1,1)-Ts(1,1))/BL(1,1))
A_4 = 2D0/(sigma**4)*(sigma**2-dist(1)**2)
Int_1_t1 = A_1_t1*Dr**2
Int_1_t2 = A_1_t2*Dr**2
Int_2 = A_2*Dr**2
Int_3 = A_3*Dr**2
Int_4 = A_4*Dr**2
DO i=2,N-1
tbar_t1 = -(2*(T(i,1)-Ts(i,1))*BL(i,1))/(3.d0*H(i,1))+T(i,1)
tbar_t2 = -(2*(T(i,3)-Ts(i,3))*BL(i,3))/(3.d0*H(i,3))+T(i,3)
A_1_t1 = tbar_t1*H(i,1)
A_1_t2 = tbar_t2*H(i,3)
A_2 = (tbar_t2*H(i,3)-tbar_t1*H(i,1))/Dt-Ts(i,1)*(H(i,3)-H(i,1))/Dt
A_3 = ((T(i,1)-Ts(i,1))/BL(i,1))
IF (dist(i)<sigma) THEN
A_4 = 2D0/(sigma**4)*(sigma**2-dist(i)**2)
ELSE
A_4 =0D0
ENDIF
Int_1_t1 = Int_1_t1 + A_1_t1*(ray(i)**2-ray(i-1)**2)
Int_1_t2 = Int_1_t2 + A_1_t2*(ray(i)**2-ray(i-1)**2)
Int_2 = Int_2 + A_2*(ray(i)**2-ray(i-1)**2)
Int_3 = Int_3 + A_3*(ray(i)**2-ray(i-1)**2)
Int_4 = Int_4 + A_4*(ray(i)**2-ray(i-1)**2)
ENDDO
tbar_t1 = -(2*(T(N,1)-Ts(N,1))*BL(N,1))/(3.d0*H(N,1))+T(N,1)
tbar_t2 = -(2*(T(N,3)-Ts(N,3))*BL(N,3))/(3.d0*H(N,3))+T(N,3)
A_1_t1 = tbar_t1*H(N,1)
A_1_t2 = tbar_t2*H(N,3)
A_2 = (tbar_t2*H(N,3)-tbar_t1*H(N,1))/Dt-Ts(N,1)*(H(N,3)-H(N,1))/Dt
A_3 = ((T(N,1)-Ts(N,1))/BL(N,1))
IF (dist(N)<sigma) THEN
A_4 = 2D0/(sigma**4)*(sigma**2-dist(N)**2)
ELSE
A_4 =0D0
ENDIF
Int_1_t1 = Int_1_t1 + A_1_t1*(dist(N)**2-ray(N-1)**2)
Int_1_t2 = Int_1_t2 + A_1_t2*(dist(N)**2-ray(N-1)**2)
Int_2 = Int_2 + A_2*(dist(N)**2-ray(N-1)**2)
Int_3 = Int_3 + A_3*(dist(N)**2-ray(N-1)**2)
Int_4 = Int_4 + A_4*(dist(N)**2-ray(N-1)**2)
! Bilan d'energie
!E_t2a: Energie Billan calculer a partir de l'intrusin au temps t+dt
BE_a = Int_1_t2
!E_t2b: Energie Bilan calculer a partir de l'intrusion aut temps t + gane perdu
BE_b = Int_1_t1+(Int_4-psi*Int_2-4D0*Pe*Int_3)*Dt
! Energie dans l'intrusion au temps t
En_t1 = Int_1_t1
! Energie dans l'intrusion au temps t+dt
En_t2 = Int_1_t2
! Energie sources en J s
Phi_s = Int_4-psi*Int_2
! Energie lost en J s
Phi_l = 4D0*Pe*Int_3
END SUBROUTINE ENERGY_CONSERVATION
END MODULE MODULE_CONSERVATION