source: branches/publications/ORCHIDEE_IPSLCM5A2.1.r5307/ORCHIDEE/src_global/gauss_jordan_method.f90 @ 7683

Last change on this file since 7683 was 2901, checked in by fabienne.maignan, 9 years ago

Correction of a spinup problem in gauss_jordan_method: initialization of variables and error handling (#153)
File size: 9.5 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 npts*nvm*(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
31MODULE 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)*veget_max(:,j)
219    ENDDO
220
221    !! Calculation current_passive_stock
222    current_passive_stock(:) = zero
223    DO j = 1, nb_veget
224       current_passive_stock(:) = current_passive_stock(:) + current_value(:,j,ipassive_pool)*veget_max(:,j)
225    ENDDO
226   
227    !! We calculate for the error for the passive pool for each pixel
228    temp_diff(:) = zero
229    temp_diff(:) =  current_passive_stock(:) - previous_passive_stock(:) 
230    WHERE ( previous_passive_stock(:) >  min_stomate )
231       error_global(:) = 100.*ABS(temp_diff(:))/previous_passive_stock(:) 
232    ELSEWHERE
233       error_global(:) = ABS(temp_diff(:))
234    ENDWHERE
235
236    !! if the criterion is reached, we can mark the point
237    WHERE (error_global(:) <= criterion)
238       flag = .TRUE.
239    ENDWHERE
240       
241  END SUBROUTINE error_L1_passive       
242
243
244END MODULE matrix_resolution
Note: See TracBrowser for help on using the repository browser.