source: branches/publications/ORCHIDEE_gmd-2018-182/src_global/gauss_jordan_method.f90 @ 8550

Last change on this file since 8550 was 1100, checked in by didier.solyga, 12 years ago

Merge Spinup_analytic branch into the trunk.
File size: 9.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 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(r_std) :: i,col,row,j,k,ii,jj        !! index (uniltess)
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     
89      !! Search the pivot (strategy of full pivoting)
90      !! We search the greatest pivot (in order to reduce errors)
91      DO i = 1,n     
92         pivot_max = 0.
93         DO  j = 1,n
94            IF(index_pivot(j) /= 1) THEN
95               DO k = 1,n
96                  IF(index_pivot(k) .EQ. 0) THEN
97                     IF(ABS(matrix_a(j,k)) .GE. pivot_max) THEN 
98                        pivot_max = ABS(matrix_a(j,k))           
99                        row = j
100                        col = k
101                     ENDIF
102                  ENDIF
103               ENDDO
104            ENDIF
105         ENDDO
106         index_pivot(col)=index_pivot(col) + 1
107
108         !! We exchange the rows and the lines if needed
109         IF(row /= col) THEN
110            DO j = 1,n
111               temp = matrix_a(row,j)
112               matrix_a(row,j) = matrix_a(col,j)
113               matrix_a(col,j) = temp
114            ENDDO
115            temp = vector_b(row)
116            vector_b(row) = vector_b(col)
117            vector_b(col) = temp
118         ENDIF
119         index_row(i) = row
120         index_col(i) = col
121         IF(matrix_a(col,col) .EQ. 0.) STOP 'the matrix A is not inversible'
122         inv_pivot = 1./matrix_a(col,col)
123         DO j = 1,n
124            matrix_a(col,j) = matrix_a(col,j) * inv_pivot 
125         ENDDO
126         vector_b(col) = vector_b(col) * inv_pivot
127         DO ii = 1,n
128            IF(ii /= col) THEN
129               temp = matrix_a(ii,col)
130               matrix_a(ii,col) = 0.
131               DO jj = 1,n
132                  matrix_a(ii,jj) = matrix_a(ii,jj) - matrix_a(col,jj)*temp
133               ENDDO
134               vector_b(ii) = vector_b(ii) - vector_b(col)*temp
135            ENDIF
136         ENDDO
137      ENDDO
138     
139      DO j = n,1,-1
140         IF(index_row(j) /= index_col(j)) THEN
141            DO i = 1,n
142               temp = matrix_a(i,index_row(j))
143               matrix_a(i,index_row(j)) = matrix_a(i,index_col(j)) 
144               matrix_a(i,index_col(j)) = temp
145          ENDDO
146       ENDIF
147    ENDDO
148   
149   
150  END SUBROUTINE gauss_jordan_method
151
152
153!! ================================================================================================================================
154!! SUBROUTINE   : error_L1_passive
155!!
156!>\BRIEF          This subroutine calculates relative errors of a vector by taking the relative error for the passive pool.
157!!
158!! DESCRIPTION  :
159!!
160!! RECENT CHANGE(S): None
161!!
162!! MAIN OUTPUT VARIABLE(S): flag is to true if the maximum relative error is less than a threshold chosen by the user.
163!!
164!! REFERENCE(S) : None
165!!
166!! FLOWCHART    : None
167!! \n
168!_ ================================================================================================================================
169
170 SUBROUTINE error_L1_passive(npts,nb_veget, nb_pools, current_value, previous_value, veget_max, criterion, flag)
171 
172    IMPLICIT NONE
173
174    !! 0. Parameters and variables declaration
175   
176    !! 0.1 Input variables   
177   
178    INTEGER(i_std), INTENT(in) :: npts                                             !! Number of continental grid cells (unitless)
179    INTEGER(i_std), INTENT(in) :: nb_veget                                         !! Number of vegetation types (2-N, unitless)
180    INTEGER(i_std), INTENT(in) :: nb_pools                                         !! Number of carbon pools (unitless)
181    REAL(r_std), DIMENSION(npts,nb_veget,nb_pools), INTENT(in) :: current_value    !! Previous values of carbon pools obtained
182                                                                                   !! by matrix resolution (gC.m^{-2})
183    REAL(r_std), DIMENSION(npts,nb_veget,nb_pools), INTENT(in) :: previous_value   !! Current values of carbon pools  obtained
184                                                                                   !! by matrix resolution (gC.m^{-2})   
185    REAL(r_std), DIMENSION(npts,nb_veget), INTENT(in) :: veget_max                 !! Fraction of vegetation (0-1, uniless)
186    REAL(r_std), INTENT(in) :: criterion                                           !! Threshold for the relativ error (0-1, unitless)
187   
188    !! 0.2 Output variables
189   
190    LOGICAL, DIMENSION(npts), INTENT(out) :: flag   !! Logical array used only inside this subroutine (true/false)
191   
192    !! 0.4 Local variables
193   
194    INTEGER(i_std) :: j                                      !! Index (unitless)
195    REAL(r_std), DIMENSION(npts) :: previous_passive_stock   !! Previous value of total passive carbon (gC)
196    REAL(r_std), DIMENSION(npts) :: current_passive_stock    !! Current value of total passive carbon (gC)
197    REAL(r_std), DIMENSION(npts) :: error_global             !! Temporary arrays containing the relative error for each grid cell
198                                                             !! (unitless)         
199    REAL(r_std), DIMENSION(npts) :: temp_diff                !! Working array storing difference values between previous_passive_stock
200                                                             !! and current_passive_stock (gC)
201
202!_ ================================================================================================================================
203
204    !! Initialize flag and error_global
205    flag(:) = .FALSE.
206    error_global(:) = zero
207
208    !! Calculation previous_passive_stock
209    previous_passive_stock(:) = zero
210    DO j = 1, nb_veget
211       previous_passive_stock(:) = previous_passive_stock(:) + previous_value(:,j,ipassive_pool)*veget_max(:,j)
212    ENDDO
213
214    !! Calculation current_passive_stock
215    current_passive_stock(:) = zero
216    DO j = 1, nb_veget
217       current_passive_stock(:) = current_passive_stock(:) + current_value(:,j,ipassive_pool)*veget_max(:,j)
218    ENDDO
219   
220    !! We calculate for the error for the passive pool for each pixel
221    temp_diff(:) = zero
222    temp_diff(:) =  current_passive_stock(:) - previous_passive_stock(:) 
223    WHERE ( previous_passive_stock(:) >  min_stomate )
224       error_global(:) = 100.*ABS(temp_diff(:))/previous_passive_stock(:) 
225    ELSEWHERE
226       error_global(:) = ABS(temp_diff(:))
227    ENDWHERE
228
229    !! if the criterion is reached, we can mark the point
230    WHERE (error_global(:) <= criterion)
231       flag = .TRUE.
232    ENDWHERE
233       
234  END SUBROUTINE error_L1_passive       
235
236
237END MODULE matrix_resolution
Note: See TracBrowser for help on using the repository browser.