! =================================================================================================================================
! MODULE       : matrix_resolution
!
! CONTACT      : orchidee-help _at_ listes.ipsl.fr
!
! LICENCE      : IPSL (2011)
! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
!
!>\BRIEF         This module solves a linear system using the Gauss-Jordan elimination method. It also calculates relative errors
!! globally and for the passive pools.        
!! 
!!
!!\n DESCRIPTION:  This module solves a linear system AX = B with the Gauss Jordan elimination method
!!                 (chosen because the system has no particular properties). \n
!!                 The code has originally picked up in Numerical recipes in Fortran 90. \n
!!                 We simplified the code in our case because we solve a basic (7,7) matrix for each point and each pft
!!                 (so we have npts*nvm*(7,7) systems to solve each time we call this routine). \n
!!                 We also calculate relative for biomass and passive carbon pools in order to test the threshold error. \n
!!
!! RECENT CHANGE(S): Didier Solyga - add subroutine for calculating relative error for passive pool.
!!
!! REFERENCE(S)	: None
!!
!! SVN          :
!! $HeadURL: $
!! $Date: $
!! $Revision:  $
!! \n
!_ ================================================================================================================================

MODULE matrix_resolution

  ! modules used 

  USE ioipsl ! for precision
  USE constantes

  IMPLICIT NONE

  CONTAINS
    

!! ================================================================================================================================
!! SUBROUTINE   : gauss-jordan_method 
!!
!>\BRIEF          This subroutine resolves a linear system by the Gauss-Jordan method.
!! (inversion of the system - complexity O(n^3)) .
!!
!! DESCRIPTION  : 
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): vector_b contains the solution of the system.
!!
!! REFERENCE(S) : None
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

    SUBROUTINE gauss_jordan_method(n,matrix_a,vector_b)

      IMPLICIT NONE

      !! 0. Variables and parameters declaration    

      !! 0.1 Input variables 
 
      INTEGER(i_std), INTENT(in) :: n  !! size of the system (1-N, unitless)

      !! 0.3 Modified variables

      REAL(r_std), DIMENSION(n,n), INTENT(inout) :: matrix_a   !! Matrix A of the linear system A.X = B
      REAL(r_std), DIMENSION(n), INTENT(inout)   :: vector_b   !! Vector B in the linear system A.X = B
      
      !! 0.4 Local Variables

      INTEGER(i_std) :: i,col,row,j,k,ii,jj        !! index (unitless)
      INTEGER(i_std), DIMENSION(n) :: index_pivot  !! vector containing the pivot index
      INTEGER(i_std), DIMENSION(n) :: index_col    !! vector containing the columns index
      INTEGER(i_std), DIMENSION(n) :: index_row    !! vector containing the rows index
      REAL(r_std) :: pivot_max,inv_pivot,temp      !! temporary variables
      
!_ ================================================================================================================================
      
      !! Initialization
      index_pivot(:) = 0
      col = 0
      row = 0
      
      !! Search the pivot (strategy of full pivoting)
      !! We search the greatest pivot (in order to reduce errors)
      DO i = 1,n     
         pivot_max = 0.
         DO  j = 1,n
            IF(index_pivot(j) /= 1) THEN
               DO k = 1,n
                  IF(index_pivot(k) .EQ. 0) THEN
                     IF(ABS(matrix_a(j,k)) .GE. pivot_max) THEN  
                        pivot_max = ABS(matrix_a(j,k))           
                        row = j
                        col = k
                     ENDIF
                  ENDIF
               ENDDO
            ENDIF
         ENDDO

         IF (col .EQ. 0) THEN
            CALL ipslerr_p (3,'gauss_jordan_method','Method failed.','','')
         ENDIF

         index_pivot(col)=index_pivot(col) + 1

         !! We exchange the rows and the lines if needed
         IF(row /= col) THEN
            DO j = 1,n
               temp = matrix_a(row,j)
               matrix_a(row,j) = matrix_a(col,j)
               matrix_a(col,j) = temp
            ENDDO
            temp = vector_b(row)
            vector_b(row) = vector_b(col)
            vector_b(col) = temp
         ENDIF
         index_row(i) = row
         index_col(i) = col
         IF(matrix_a(col,col) .EQ. 0.) STOP 'the matrix A is not inversible'
         inv_pivot = 1./matrix_a(col,col)
         DO j = 1,n
            matrix_a(col,j) = matrix_a(col,j) * inv_pivot 
         ENDDO
         vector_b(col) = vector_b(col) * inv_pivot
         DO ii = 1,n
            IF(ii /= col) THEN
               temp = matrix_a(ii,col)
               matrix_a(ii,col) = 0.
               DO jj = 1,n
                  matrix_a(ii,jj) = matrix_a(ii,jj) - matrix_a(col,jj)*temp
               ENDDO
               vector_b(ii) = vector_b(ii) - vector_b(col)*temp
            ENDIF
         ENDDO
      ENDDO
      
      DO j = n,1,-1
         IF(index_row(j) /= index_col(j)) THEN
            DO i = 1,n
               temp = matrix_a(i,index_row(j))
               matrix_a(i,index_row(j)) = matrix_a(i,index_col(j)) 
               matrix_a(i,index_col(j)) = temp
          ENDDO
       ENDIF
    ENDDO
    
    
  END SUBROUTINE gauss_jordan_method


!! ================================================================================================================================
!! SUBROUTINE   : error_L1_passive
!!
!>\BRIEF          This subroutine calculates relative errors of a vector by taking the relative error for the passive pool.
!!
!! DESCRIPTION  : 
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): flag is to true if the maximum relative error is less than a threshold chosen by the user.
!!
!! REFERENCE(S) : None
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

 SUBROUTINE error_L1_passive(npts,nb_veget, nb_pools, current_value, previous_value, veget_max, criterion, flag)
  
    IMPLICIT NONE

    !! 0. Parameters and variables declaration
    
    !! 0.1 Input variables    
    
    INTEGER(i_std), INTENT(in) :: npts                                             !! Number of continental grid cells (unitless)
    INTEGER(i_std), INTENT(in) :: nb_veget                                         !! Number of vegetation types (2-N, unitless)
    INTEGER(i_std), INTENT(in) :: nb_pools                                         !! Number of carbon pools (unitless)
    REAL(r_std), DIMENSION(npts,nb_veget,nb_pools), INTENT(in) :: current_value    !! Previous values of carbon pools obtained
                                                                                   !! by matrix resolution (gC.m^{-2})
    REAL(r_std), DIMENSION(npts,nb_veget,nb_pools), INTENT(in) :: previous_value   !! Current values of carbon pools  obtained
                                                                                   !! by matrix resolution (gC.m^{-2})   
    REAL(r_std), DIMENSION(npts,nb_veget), INTENT(in) :: veget_max                 !! Fraction of vegetation (0-1, uniless)
    REAL(r_std), INTENT(in) :: criterion                                           !! Threshold for the relativ error (0-1, unitless) 
    
    !! 0.2 Output variables
    
    LOGICAL, DIMENSION(npts), INTENT(out) :: flag   !! Logical array used only inside this subroutine (true/false)
    
    !! 0.4 Local variables
    
    INTEGER(i_std) :: j                                      !! Index (unitless)
    REAL(r_std), DIMENSION(npts) :: previous_passive_stock   !! Previous value of total passive carbon (gC) 
    REAL(r_std), DIMENSION(npts) :: current_passive_stock    !! Current value of total passive carbon (gC)
    REAL(r_std), DIMENSION(npts) :: error_global             !! Temporary arrays containing the relative error for each grid cell 
                                                             !! (unitless)         
    REAL(r_std), DIMENSION(npts) :: temp_diff                !! Working array storing difference values between previous_passive_stock
                                                             !! and current_passive_stock (gC)

!_ ================================================================================================================================

    !! Initialize flag and error_global
    flag(:) = .FALSE.
    error_global(:) = zero

    !! Calculation previous_passive_stock
    previous_passive_stock(:) = zero
    DO j = 1, nb_veget
       previous_passive_stock(:) = previous_passive_stock(:) + previous_value(:,j,ipassive_pool)*veget_max(:,j)
    ENDDO

    !! Calculation current_passive_stock
    current_passive_stock(:) = zero
    DO j = 1, nb_veget
       current_passive_stock(:) = current_passive_stock(:) + current_value(:,j,ipassive_pool)*veget_max(:,j)
    ENDDO
    
    !! We calculate for the error for the passive pool for each pixel
    temp_diff(:) = zero
    temp_diff(:) =  current_passive_stock(:) - previous_passive_stock(:) 
    WHERE ( previous_passive_stock(:) >  min_stomate )
       error_global(:) = 100.*ABS(temp_diff(:))/previous_passive_stock(:) 
    ELSEWHERE
       error_global(:) = ABS(temp_diff(:))
    ENDWHERE

    !! if the criterion is reached, we can mark the point 
    WHERE (error_global(:) <= criterion)
       flag = .TRUE.
    ENDWHERE
       
  END SUBROUTINE error_L1_passive       


END MODULE matrix_resolution
