tridiag_ground_snowcro.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.
MODULE modi_tridiag_ground_snowcro 
!
INTERFACE tridiag_ground_snowcro
!
SUBROUTINE tridiag_ground_snowcro_1d(PA,PB,PC,PY,PX,KNLVLS_USE,KDIFLOOP)
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 
!
INTEGER,    DIMENSION(:), INTENT(IN) :: knlvls_use ! number of effective layers 
!
INTEGER, INTENT(IN) :: kdifloop       ! shift in control loops: 0 or 1
END SUBROUTINE tridiag_ground_snowcro_1d
!
SUBROUTINE tridiag_ground_snowcro_2d(PA,PB,PC,PY,PX,KNLVLS_USE,KDIFLOOP)
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 
!
INTEGER,    DIMENSION(:,:), INTENT(IN) :: knlvls_use ! number of effective layers 
!
INTEGER, INTENT(IN) :: kdifloop       ! shift in control loops: 0 or 1
END SUBROUTINE tridiag_ground_snowcro_2d
!
END INTERFACE tridiag_ground_snowcro
!
END MODULE modi_tridiag_ground_snowcro
!
!###################################################################
       SUBROUTINE tridiag_ground_snowcro_1d(PA,PB,PC,PY,PX,KNLVLS_USE,KDIFLOOP)
!      #########################################
!
!
!!****   *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
!!       05/2011: Brun  Special treatment to tackle the variable number
!!                      of snow layers 
!!                      In case of second call, a shift of 1 snow layer
!!                      is applied in the control loops.
!! ---------------------------------------------------------------------
!
!*       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 
!
INTEGER,    DIMENSION(:), INTENT(IN) :: knlvls_use ! number of effective layers 
!
INTEGER, INTENT(IN) :: kdifloop       ! shift in control loops: 0 or 1
!
!*       0.2 declarations of local variables
!
REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2)) :: zw   ! work array
REAL, DIMENSION(SIZE(PA,1)           ) :: zdet ! work array
!
!
INTEGER           :: jl, jj            ! vertical loop control
INTEGER           :: inl               ! number of vertical levels
!
REAL(KIND=JPRB) :: zhook_handle
! ---------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('TRIDIAG_GROUND_SNOWCRO_1D',0,zhook_handle)
!
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 jl = 2,inl
  !
  DO jj = 1,SIZE(px,1)
    !
    IF ( jl<=knlvls_use(jj)-kdifloop ) THEN
      !
      zw(jj,jl) = pc(jj,jl-1)/zdet(jj)
      zdet(jj)  = pb(jj,jl  ) - pa(jj,jl)*zw(jj,jl)
      !
      px(jj,jl) = ( py(jj,jl) - pa(jj,jl)*px(jj,jl-1) ) / zdet(jj)
      !
    ENDIF
    !
  ENDDO
  !
END DO
!
!-------------------------------------------------------------------------------
!
!*       2.  levels going down
!            -----------------
!
DO jl = inl-1,1,-1
  !
  DO jj = 1,SIZE(px,1)
    !
    IF ( jl<=knlvls_use(jj)-1-kdifloop ) THEN
      !
      px(jj,jl) = px(jj,jl) - zw(jj,jl+1)*px(jj,jl+1)
      !
    ENDIF
    !
  ENDDO
  !
END DO
!
IF (lhook) CALL dr_hook('TRIDIAG_GROUND_SNOWCRO_1D',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE tridiag_ground_snowcro_1d
!
!###################################################################
       SUBROUTINE tridiag_ground_snowcro_2d(PA,PB,PC,PY,PX,KNLVLS_USE,KDIFLOOP)
!      #########################################
!
!
!!****   *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
!!       05/2011: Brun  Special treatment to tackle the variable number
!!                      of snow layers 
!!                      In case of second call, a shift of 1 snow layer
!!                      is applied in the control loops.
!! ---------------------------------------------------------------------
!
!*       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 
!
INTEGER,    DIMENSION(:,:), INTENT(IN) :: knlvls_use ! number of effective layers 
!
INTEGER, INTENT(IN) :: kdifloop       ! shift in control loops: 0 or 1
!
!*       0.2 declarations of local variables
!
REAL, DIMENSION(SIZE(PA,2)) :: zw   ! work array
REAL :: zdet ! work array
!
!
INTEGER           :: jl, jj, jb     ! vertical loop control
INTEGER           :: inl, inb            ! number of vertical levels
!
REAL(KIND=JPRB) :: zhook_handle
! ---------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('TRIDIAG_GROUND_SNOWCRO_2D',0,zhook_handle)
!
inl = SIZE(px,2)
inb = SIZE(px,3)
!
!*       1.  levels going up
!            ---------------
!
!*       1.1 first level
!            -----------
!
DO jb = 1,inb
  !
  DO jj = 1,SIZE(px,1)
    !
    zdet = pb(jj,1,jb)
    !
    px(jj,1,jb) = py(jj,1,jb) / zdet
    !
    !*       1.2 other levels
    !            ------------
    !
    DO jl = 2,inl
      !
      IF ( jl<=knlvls_use(jj,jb)-kdifloop ) THEN
        !
        zw(jl) = pc(jj,jl-1,jb)  /zdet
        zdet   = pb(jj,jl  ,jb) - pa(jj,jl,jb) * zw(jl)
        !
        px(jj,jl,jb) = ( py(jj,jl,jb) - pa(jj,jl,jb) * px(jj,jl-1,jb) ) / zdet
        !
      ENDIF
      !
    ENDDO
    !
    !-------------------------------------------------------------------------------
    !
    !*       2.  levels going down
    !            -----------------
    !
    DO jl = inl-1,1,-1
      !
      IF ( jl<=knlvls_use(jj,jb)-1-kdifloop ) THEN
        !
        px(jj,jl,jb) = px(jj,jl,jb) - zw(jl+1) * px(jj,jl+1,jb)
        !
      ENDIF
      !
    ENDDO
    !
  END DO
  !
ENDDO
!
IF (lhook) CALL dr_hook('TRIDIAG_GROUND_SNOWCRO_2D',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE tridiag_ground_snowcro_2d