-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathemprop.f90
151 lines (89 loc) · 3.29 KB
/
emprop.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
MODULE EMPROP
CONTAINS
SUBROUTINE upwind_fund(w_em, omega, l, delta_x)
USE CONSTANTS
USE PARAMETERS
IMPLICIT NONE
REAL(KIND=8), DIMENSION(-n_kv:, :), INTENT(INOUT) :: w_em
REAL(KIND=8), DIMENSION(:):: omega, l, delta_x
REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE:: w_em_new
REAL(KIND=8), DIMENSION(-n_kv:n_kv):: omega_emf_om, v_grp, coeff2
INTEGER(KIND=4)::i, j
ALLOCATE(w_em_new(-n_kv:n_kv, n_xp+3))
omega_emF_om = sqrt(1.+v_c**2/v_t**2 *k_f**2)
v_grp = v_c*v_c*k_f/v_t/omega_emF_om
coeff2 = v_grp*k_f
w_em_new = w_em
DO j=3, n_xp+2
DO i = -n_kv,-1
w_em_new(i,j) = w_em_new(i,j) - dt*v_grp(i)*(w_em(i,j+1) - w_em(i,j))/delta_x(j)
ENDDO
DO i = 1, n_kv
w_em_new(i,j) = w_em_new(i,j) - dt*v_grp(i)*(w_em(i,j) - w_em(i,j-1))/delta_x(j-1)
ENDDO
DO i = -n_kv, -1
w_em_new(i,j) = w_em_new(i,j) + dt *(1./omega_emf_om(i) + coeff2(i))* &
max(-l(j), 0.d0)*( w_em(i,j) - &
w_em(i+1,j))/(k_f(i)-k_f(i+1))
ENDDO
DO i = 1, n_kv
w_em_new(i,j) = w_em_new(i,j) - dt *(1./omega_emf_om(i) + coeff2(i))* &
max(-l(j), 0.d0)*( w_em(i,j) - &
w_em(i-1,j))/(k_f(i)-k_f(i-1))
ENDDO
DO i = -n_kv+1, -1
w_em_new(i,j) = w_em_new(i,j) - dt *(1./omega_emf_om(i) + coeff2(i))* &
max(l(j), 0.d0)*( w_em(i,j) - &
w_em(i-1,j))/(k_f(i)-k_f(i-1))
ENDDO
i= -n_kv
! w_em_new(i,j) = w_em_new(i,j) + dt *(1./omega_emf_om(i) + coeff2(i))* &
! max(l(j), 0.d0)*( w_em(i,j))/(k_f(i+1)-k_f(i))
DO i = 1, n_kv-1
w_em_new(i,j) = w_em_new(i,j) + dt *(1./omega_emf_om(i) + coeff2(i))* &
max(l(j), 0.d0)*( w_em(i+1,j) - &
w_em(i,j))/(k_f(i+1)-k_f(i))
ENDDO
i=n_kv
! w_em_new(i,j) = w_em_new(i,j) + dt *(1./omega_emf_om(i) + coeff2(i))* &
! max(l(j), 0.d0)*(- w_em(i,j))/(k_f(i)-k_f(i-1))
ENDDO
w_em = w_em_new
DEALLOCATE(w_em_new)
END SUBROUTINE upwind_fund
SUBROUTINE upwind_harm(w_em, omega, l, delta_x)
USE CONSTANTS
USE PARAMETERS
IMPLICIT NONE
REAL(KIND=8), DIMENSION(-n_kv:, :), INTENT(INOUT) :: w_em
REAL(KIND=8), DIMENSION(:):: omega, l, delta_x
REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE:: w_em_new
REAL(KIND=8), DIMENSION(-n_kv:n_kv):: omega_emH_om, v_grp, coeff2
INTEGER(KIND=4)::i, j
ALLOCATE(w_em_new(-n_kv:n_kv, n_xp+3))
omega_emH_om = sqrt(1.+v_c**2/v_t**2 *k_h**2)
v_grp = v_c*v_c*k_h/v_t/omega_emH_om
coeff2 = v_grp*k_h
w_em_new = w_em
DO j=3, n_xp+2
w_em_new(:,j) = w_em_new(:,j) - dt*v_grp*(w_em(:,j) - w_em(:,j-1))/delta_x(j-1)
DO i = -n_kv, n_kv-1
w_em_new(i,j) = w_em_new(i,j) + dt * (1./omega_emH_om(i) + coeff2(i))* &
max(l(j), 0.d0)*( w_em(i+1,j) - w_em(i,j))/(k_h(i+1)-k_h(i))
ENDDO
i=n_kv
w_em_new(i,j) = w_em_new(i,j) + dt * (1./omega_emH_om(i) + coeff2(i))* &
max(l(j), 0.d0)*( -1.*w_em(i,j))/(k_h(i)-k_h(i-1))
DO i = -n_kv+1, n_kv
w_em_new(i,j) = w_em_new(i,j) - dt *(1./omega_emH_om(i) + coeff2(i))* &
max(-l(j), 0.d0)*( w_em(i,j) - &
w_em(i-1,j))/(k_h(i)-k_h(i-1))
ENDDO
i=-n_kv
w_em_new(i,j) = w_em_new(i,j) - dt * (1./omega_emH_om(i) + coeff2(i))* &
max(-l(j), 0.d0)*( w_em(i,j))/(k_h(i+1)-k_h(i))
ENDDO
w_em = w_em_new
DEALLOCATE(w_em_new)
END SUBROUTINE upwind_harm
END MODULE EMPROP