tridiag_ground.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
!     #########
       SUBROUTINE tridiag_ground(PA,PB,PC,PY,PX)
!      #########################################
!
!
!!****   *TRIDIAG_GROUND* - routine to solve a time implicit scheme
!!
!!
!!     PURPOSE
!!     -------
!        The purpose of this routine is to resolve the linear system:
!
!       A.X = Y
!
!      where A is a tridiagonal matrix, and X and Y two vertical vectors.
!     However, the computations are performed at the same time for all
!     the verticals where an inversion of the system is necessary.
!     This explain the dimansion of the input variables.
!
!!**   METHOD
!!     ------
!!                      
!!        Then, the classical tridiagonal algorithm is used to invert the 
!!     implicit operator. Its matrix is given by:
!!
!!     (  b(1)      c(1)      0        0        0         0        0        0  )
!!     (  a(2)      b(2)     c(2)      0  ...    0        0        0        0  ) 
!!     (   0        a(3)     b(3)     c(3)       0        0        0        0  ) 
!!      .......................................................................
!!     (   0   ...   0      a(k)      b(k)     c(k)       0   ...  0        0  ) 
!!      .......................................................................
!!     (   0         0        0        0        0 ...  a(n-1)   b(n-1)   c(n-1))
!!     (   0         0        0        0        0 ...     0      a(n)     b(n) )
!!
!!
!!       All these computations are purely vertical and vectorizations are 
!!     easely achieved by processing all the verticals in parallel.
!!
!!     EXTERNAL
!!     --------
!!
!!       NONE
!!
!!     IMPLICIT ARGUMENTS
!!     ------------------
!!
!!     REFERENCE
!!     ---------
!!
!!     AUTHOR
!!     ------
!!       V. Masson
!! 
!!     MODIFICATIONS
!!     -------------
!!       Original        May 13, 1998
!!       Modified :
!!       B. Decharme  08/12 Loop optimization
!! ---------------------------------------------------------------------
!
!*       0. DECLARATIONS
!
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!
!*       0.1 declarations of arguments
!
REAL,    DIMENSION(:,:), INTENT(IN)  :: pa  ! lower diag. elements of A matrix
REAL,    DIMENSION(:,:), INTENT(IN)  :: pb  ! main  diag. elements of A matrix
REAL,    DIMENSION(:,:), INTENT(IN)  :: pc  ! upper diag. elements of A matrix
REAL,    DIMENSION(:,:), INTENT(IN)  :: py  ! r.h.s. term   
!
REAL,    DIMENSION(:,:), INTENT(OUT) :: px  ! solution of A.X = Y 
!
!*       0.2 declarations of local variables
!
INTEGER           :: ji             ! number of point loop control
INTEGER           :: jk             ! vertical loop control
INTEGER           :: ini            ! number of point
INTEGER           :: inl            ! number of vertical levels
!
REAL, DIMENSION(SIZE(PA,1)           ) :: zdet ! work array
REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2)) :: zw   ! work array
REAL(KIND=JPRB) :: zhook_handle
! ---------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('TRIDIAG_GROUND',0,zhook_handle)
ini=SIZE(px,1)
inl=SIZE(px,2)
!
!*       1.  levels going up
!            ---------------
!
!*       1.1 first level
!            -----------
!
zdet(:)   = pb(:,1)

px(:,1) = py(:,1) / zdet(:)
!
!*       1.2 other levels
!            ------------
!
DO jk=2,inl
   DO ji=1,ini
      zw(ji,jk)  = pc(ji,jk-1)/zdet(ji)
      zdet(ji)     = pb(ji,jk  ) - pa(ji,jk)*zw(ji,jk)
      px(ji,jk)  = ( py(ji,jk) - pa(ji,jk)*px(ji,jk-1) ) / zdet(ji)
   END DO
END DO
!
!-------------------------------------------------------------------------------
!
!*       2.  levels going down
!            -----------------
!
DO jk=inl-1,1,-1
   DO ji=1,ini
      px(ji,jk) = px(ji,jk) - zw(ji,jk+1)*px(ji,jk+1)
   END DO
END DO
IF (lhook) CALL dr_hook('TRIDIAG_GROUND',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE tridiag_ground