Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

gauss_jordan_method.f90 @ 7346

Last change on this file since 7346 was 3549, checked in by bertrand.guenet, 8 years ago
update between the trunk rev 3340 and SOM it's still not working
File size: 12.4 KB

Line
1	! =================================================================================================================================
2	! MODULE : matrix_resolution
3	!
4	! CONTACT : orchidee-help _at_ ipsl.jussieu.fr
5	!
6	! LICENCE : IPSL (2011)
7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8	!
9	!>\BRIEF This module solves a linear system using the Gauss-Jordan elimination method. It also calculates relative errors
10	!! globally and for the passive pools.
11	!!
12	!!
13	!!\n DESCRIPTION: This module solves a linear system AX = B with the Gauss Jordan elimination method
14	!! (chosen because the system has no particular properties). \n
15	!! The code has originally picked up in Numerical recipes in Fortran 90. \n
16	!! We simplified the code in our case because we solve a basic (7,7) matrix for each point and each pft
17	!! (so we have nptsnvm(7,7) systems to solve each time we call this routine). \n
18	!! We also calculate relative for biomass and passive carbon pools in order to test the threshold error. \n
19	!!
20	!! RECENT CHANGE(S): Didier Solyga - add subroutine for calculating relative error for passive pool.
21	!!
22	!! REFERENCE(S) : None
23	!!
24	!! SVN :
25	!! $HeadURL: $
26	!! $Date: $
27	!! $Revision: $
28	!! \n
29	!_ ================================================================================================================================
30
31	MODULE matrix_resolution
32
33	! modules used
34
35	USE ioipsl ! for precision
36	USE constantes
37
38	IMPLICIT NONE
39
40	CONTAINS
41
42
43	!! ================================================================================================================================
44	!! SUBROUTINE : gauss-jordan_method
45	!!
46	!>\BRIEF This subroutine resolves a linear system by the Gauss-Jordan method.
47	!! (inversion of the system - complexity O(n^3)) .
48	!!
49	!! DESCRIPTION :
50	!!
51	!! RECENT CHANGE(S): None
52	!!
53	!! MAIN OUTPUT VARIABLE(S): vector_b contains the solution of the system.
54	!!
55	!! REFERENCE(S) : None
56	!!
57	!! FLOWCHART : None
58	!! \n
59	!_ ================================================================================================================================
60
61	SUBROUTINE gauss_jordan_method(n,matrix_a,vector_b)
62
63	IMPLICIT NONE
64
65	!! 0. Variables and parameters declaration
66
67	!! 0.1 Input variables
68
69	INTEGER(i_std), INTENT(in) :: n !! size of the system (1-N, unitless)
70
71	!! 0.3 Modified variables
72
73	REAL(r_std), DIMENSION(n,n), INTENT(inout) :: matrix_a !! Matrix A of the linear system A.X = B
74	REAL(r_std), DIMENSION(n), INTENT(inout) :: vector_b !! Vector B in the linear system A.X = B
75
76	!! 0.4 Local Variables
77
78	INTEGER(i_std) :: i,col,row,j,k,ii,jj !! index (unitless)
79	INTEGER(i_std), DIMENSION(n) :: index_pivot !! vector containing the pivot index
80	INTEGER(i_std), DIMENSION(n) :: index_col !! vector containing the columns index
81	INTEGER(i_std), DIMENSION(n) :: index_row !! vector containing the rows index
82	REAL(r_std) :: pivot_max,inv_pivot,temp !! temporary variables
83
84	!_ ================================================================================================================================
85
86	!! Initialization
87	index_pivot(:) = 0
88	col = 0
89	row = 0
90
91	!! Search the pivot (strategy of full pivoting)
92	!! We search the greatest pivot (in order to reduce errors)
93	DO i = 1,n
94	pivot_max = 0.
95	DO j = 1,n
96	IF(index_pivot(j) /= 1) THEN
97	DO k = 1,n
98	IF(index_pivot(k) .EQ. 0) THEN
99	IF(ABS(matrix_a(j,k)) .GE. pivot_max) THEN
100	pivot_max = ABS(matrix_a(j,k))
101	row = j
102	col = k
103	ENDIF
104	ENDIF
105	ENDDO
106	ENDIF
107	ENDDO
108
109	IF (col .EQ. 0) THEN
110	CALL ipslerr_p (3,'gauss_jordan_method','Method failed.','','')
111	ENDIF
112
113	index_pivot(col)=index_pivot(col) + 1
114
115	!! We exchange the rows and the lines if needed
116	IF(row /= col) THEN
117	DO j = 1,n
118	temp = matrix_a(row,j)
119	matrix_a(row,j) = matrix_a(col,j)
120	matrix_a(col,j) = temp
121	ENDDO
122	temp = vector_b(row)
123	vector_b(row) = vector_b(col)
124	vector_b(col) = temp
125	ENDIF
126	index_row(i) = row
127	index_col(i) = col
128	IF(matrix_a(col,col) .EQ. 0.) STOP 'the matrix A is not inversible'
129	inv_pivot = 1./matrix_a(col,col)
130	DO j = 1,n
131	matrix_a(col,j) = matrix_a(col,j) * inv_pivot
132	ENDDO
133	vector_b(col) = vector_b(col) * inv_pivot
134	DO ii = 1,n
135	IF(ii /= col) THEN
136	temp = matrix_a(ii,col)
137	matrix_a(ii,col) = 0.
138	DO jj = 1,n
139	matrix_a(ii,jj) = matrix_a(ii,jj) - matrix_a(col,jj)*temp
140	ENDDO
141	vector_b(ii) = vector_b(ii) - vector_b(col)*temp
142	ENDIF
143	ENDDO
144	ENDDO
145
146	DO j = n,1,-1
147	IF(index_row(j) /= index_col(j)) THEN
148	DO i = 1,n
149	temp = matrix_a(i,index_row(j))
150	matrix_a(i,index_row(j)) = matrix_a(i,index_col(j))
151	matrix_a(i,index_col(j)) = temp
152	ENDDO
153	ENDIF
154	ENDDO
155
156
157	END SUBROUTINE gauss_jordan_method
158
159
160	!! ================================================================================================================================
161	!! SUBROUTINE : error_L1_passive
162	!!
163	!>\BRIEF This subroutine calculates relative errors of a vector by taking the relative error for the passive pool.
164	!!
165	!! DESCRIPTION :
166	!!
167	!! RECENT CHANGE(S): None
168	!!
169	!! MAIN OUTPUT VARIABLE(S): flag is to true if the maximum relative error is less than a threshold chosen by the user.
170	!!
171	!! REFERENCE(S) : None
172	!!
173	!! FLOWCHART : None
174	!! \n
175	!_ ================================================================================================================================
176
177	SUBROUTINE error_L1_passive(npts,nb_veget, nb_pools, current_value, previous_value, veget_max, criterion, flag)
178
179	IMPLICIT NONE
180
181	!! 0. Parameters and variables declaration
182
183	!! 0.1 Input variables
184
185	INTEGER(i_std), INTENT(in) :: npts !! Number of continental grid cells (unitless)
186	INTEGER(i_std), INTENT(in) :: nb_veget !! Number of vegetation types (2-N, unitless)
187	INTEGER(i_std), INTENT(in) :: nb_pools !! Number of carbon pools (unitless)
188	REAL(r_std), DIMENSION(npts,nb_veget,nb_pools), INTENT(in) :: current_value !! Previous values of carbon pools obtained
189	!! by matrix resolution (gC.m^{-2})
190	REAL(r_std), DIMENSION(npts,nb_veget,nb_pools), INTENT(in) :: previous_value !! Current values of carbon pools obtained
191	!! by matrix resolution (gC.m^{-2})
192	REAL(r_std), DIMENSION(npts,nb_veget), INTENT(in) :: veget_max !! Fraction of vegetation (0-1, uniless)
193	REAL(r_std), INTENT(in) :: criterion !! Threshold for the relativ error (0-1, unitless)
194
195	!! 0.2 Output variables
196
197	LOGICAL, DIMENSION(npts), INTENT(out) :: flag !! Logical array used only inside this subroutine (true/false)
198
199	!! 0.4 Local variables
200
201	INTEGER(i_std) :: j !! Index (unitless)
202	REAL(r_std), DIMENSION(npts) :: previous_passive_stock !! Previous value of total passive carbon (gC)
203	REAL(r_std), DIMENSION(npts) :: current_passive_stock !! Current value of total passive carbon (gC)
204	REAL(r_std), DIMENSION(npts) :: error_global !! Temporary arrays containing the relative error for each grid cell
205	!! (unitless)
206	REAL(r_std), DIMENSION(npts) :: temp_diff !! Working array storing difference values between previous_passive_stock
207	!! and current_passive_stock (gC)
208
209	!_ ================================================================================================================================
210
211	!! Initialize flag and error_global
212	flag(:) = .FALSE.
213	error_global(:) = zero
214
215	!! Calculation previous_passive_stock
216	previous_passive_stock(:) = zero
217	DO j = 1, nb_veget
218	previous_passive_stock(:) = previous_passive_stock(:) + previous_value(:,j,ipassive_pool_z1)*veget_max(:,j)
219	ENDDO
220	DO j = 1, nb_veget
221	previous_passive_stock(:) = previous_passive_stock(:) + previous_value(:,j,ipassive_pool_z2)*veget_max(:,j)
222	ENDDO
223	DO j = 1, nb_veget
224	previous_passive_stock(:) = previous_passive_stock(:) + previous_value(:,j,ipassive_pool_z3)*veget_max(:,j)
225	ENDDO
226	DO j = 1, nb_veget
227	previous_passive_stock(:) = previous_passive_stock(:) + previous_value(:,j,ipassive_pool_z4)*veget_max(:,j)
228	ENDDO
229	DO j = 1, nb_veget
230	previous_passive_stock(:) = previous_passive_stock(:) + previous_value(:,j,ipassive_pool_z5)*veget_max(:,j)
231	ENDDO
232	DO j = 1, nb_veget
233	previous_passive_stock(:) = previous_passive_stock(:) + previous_value(:,j,ipassive_pool_z6)*veget_max(:,j)
234	ENDDO
235	DO j = 1, nb_veget
236	previous_passive_stock(:) = previous_passive_stock(:) + previous_value(:,j,ipassive_pool_z7)*veget_max(:,j)
237	ENDDO
238	DO j = 1, nb_veget
239	previous_passive_stock(:) = previous_passive_stock(:) + previous_value(:,j,ipassive_pool_z8)*veget_max(:,j)
240	ENDDO
241	DO j = 1, nb_veget
242	previous_passive_stock(:) = previous_passive_stock(:) + previous_value(:,j,ipassive_pool_z9)*veget_max(:,j)
243	ENDDO
244	DO j = 1, nb_veget
245	previous_passive_stock(:) = previous_passive_stock(:) + previous_value(:,j,ipassive_pool_z10)*veget_max(:,j)
246	ENDDO
247	DO j = 1, nb_veget
248	previous_passive_stock(:) = previous_passive_stock(:) + previous_value(:,j,ipassive_pool_z11)*veget_max(:,j)
249	ENDDO
250
251	!! Calculation current_passive_stock
252	current_passive_stock(:) = zero
253	DO j = 1, nb_veget
254	current_passive_stock(:) = current_passive_stock(:) + current_value(:,j,ipassive_pool_z1)*veget_max(:,j)
255	ENDDO
256	DO j = 1, nb_veget
257	current_passive_stock(:) = current_passive_stock(:) + current_value(:,j,ipassive_pool_z2)*veget_max(:,j)
258	ENDDO
259	DO j = 1, nb_veget
260	current_passive_stock(:) = current_passive_stock(:) + current_value(:,j,ipassive_pool_z3)*veget_max(:,j)
261	ENDDO
262	DO j = 1, nb_veget
263	current_passive_stock(:) = current_passive_stock(:) + current_value(:,j,ipassive_pool_z4)*veget_max(:,j)
264	ENDDO
265	DO j = 1, nb_veget
266	current_passive_stock(:) = current_passive_stock(:) + current_value(:,j,ipassive_pool_z5)*veget_max(:,j)
267	ENDDO
268	DO j = 1, nb_veget
269	current_passive_stock(:) = current_passive_stock(:) + current_value(:,j,ipassive_pool_z6)*veget_max(:,j)
270	ENDDO
271	DO j = 1, nb_veget
272	current_passive_stock(:) = current_passive_stock(:) + current_value(:,j,ipassive_pool_z7)*veget_max(:,j)
273	ENDDO
274	DO j = 1, nb_veget
275	current_passive_stock(:) = current_passive_stock(:) + current_value(:,j,ipassive_pool_z8)*veget_max(:,j)
276	ENDDO
277	DO j = 1, nb_veget
278	current_passive_stock(:) = current_passive_stock(:) + current_value(:,j,ipassive_pool_z9)*veget_max(:,j)
279	ENDDO
280	DO j = 1, nb_veget
281	current_passive_stock(:) = current_passive_stock(:) + current_value(:,j,ipassive_pool_z10)*veget_max(:,j)
282	ENDDO
283	DO j = 1, nb_veget
284	current_passive_stock(:) = current_passive_stock(:) + current_value(:,j,ipassive_pool_z11)*veget_max(:,j)
285	ENDDO
286
287	!! We calculate for the error for the passive pool for each pixel
288	temp_diff(:) = zero
289	temp_diff(:) = current_passive_stock(:) - previous_passive_stock(:)
290	WHERE ( previous_passive_stock(:) > min_stomate )
291	error_global(:) = 100.*ABS(temp_diff(:))/previous_passive_stock(:)
292	ELSEWHERE
293	error_global(:) = ABS(temp_diff(:))
294	ENDWHERE
295
296	!! if the criterion is reached, we can mark the point
297	WHERE (error_global(:) <= criterion)
298	flag = .TRUE.
299	ENDWHERE
300
301	END SUBROUTINE error_L1_passive
302
303
304	END MODULE matrix_resolution

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: branches/publications/ORCHIDEE-MUSLE-r6129/src_global/gauss_jordan_method.f90 @ 7346

Download in other formats: