ecearth3/sources/nemo-3.6/NEMO/TOP_SRC/PISCES/SED/sedmat.F90

MODULE sedmat
#if defined key_sed
   !!======================================================================
   !!              ***  MODULE  sedmat  ***
   !!    Sediment : linear system of equations
   !!=====================================================================
   !! * Modules used
   !!----------------------------------------------------------------------

   USE sed     ! sediment global variable

   IMPLICIT NONE
   PRIVATE

   PUBLIC sed_mat 

   INTERFACE sed_mat
      MODULE PROCEDURE sed_mat_dsr, sed_mat_btb
   END INTERFACE

   INTEGER, PARAMETER :: nmax = 30 


   !! $Id: sedmat.F90 2355 2015-05-20 07:11:50Z ufla $
 CONTAINS

    SUBROUTINE sed_mat_dsr( nvar, ndim, nlev, preac, psol )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE sed_mat_dsr  ***
       !!
       !! ** Purpose :  solves tridiagonal system of linear equations 
       !!
       !! ** Method  : 
       !!        1 - computes left hand side of linear system of equations
       !!            for dissolution reaction
       !!         For mass balance in kbot+sediment :
       !!              dz3d  (:,1) = dz(1) = 0.5 cm
       !!              volw3d(:,1) = dzkbot ( see sedini.F90 ) 
       !!              dz(2)       = 0.3 cm 
       !!              dz3d(:,2)   = 0.3 + dzdep   ( see seddsr.F90 )     
       !!              volw3d(:,2) and vols3d(l,2) are thickened ( see seddsr.F90 ) 
       !!
       !!         2 - forward/backward substitution. 
       !!
       !!   History :
       !!        !  04-10 (N. Emprin, M. Gehlen ) original
       !!        !  06-04 (C. Ethe)  Module Re-organization
       !!----------------------------------------------------------------------
       !! * Arguments
       INTEGER , INTENT(in) ::  nvar  ! number of variables
       INTEGER , INTENT(in) ::  ndim  ! number of points
       INTEGER , INTENT(in) ::  nlev  ! number of sediment levels

       REAL(wp), DIMENSION(ndim,nlev,nvar), INTENT(in   ) :: preac  ! reaction rates
       REAL(wp), DIMENSION(ndim,nlev,nvar), INTENT(inout) :: psol   ! solution ( undersaturation values )
 
       !---Local declarations
       INTEGER  ::  ji, jk, jn
       REAL(wp), DIMENSION(ndim,nlev) :: za, zb, zc, zr
       REAL(wp), DIMENSION(ndim)      :: zbet
       REAL(wp), DIMENSION(ndim,nmax) :: zgamm

       REAL(wp) ::  aplus,aminus  
       REAL(wp) ::  rplus,rminus   
       REAL(wp) ::  dxplus,dxminus

       !----------------------------------------------------------------------


       ! Computation left hand side of linear system of 
       ! equations for dissolution reaction
       !---------------------------------------------


       ! first sediment level          
       DO ji = 1, ndim
          aplus  = ( ( volw3d(ji,1) / dz3d(ji,1) ) + &
                     ( volw3d(ji,2) / dz3d(ji,2) ) ) / 2.
          dxplus = ( dz3d(ji,1) + dz3d(ji,2) ) / 2.
          rplus  = ( dtsed / volw3d(ji,1) ) * diff(1) * aplus / dxplus 

          za(ji,1) = 0.
          zb(ji,1) = 1. + rplus
          zc(ji,1) = -rplus
       ENDDO
 
       DO jk = 2, nlev - 1
          DO ji = 1, ndim
             aminus  = ( ( volw3d(ji,jk-1) / dz3d(ji,jk-1) ) + &
                &        ( volw3d(ji,jk  ) / dz3d(ji,jk  ) ) ) / 2.
             dxminus = ( dz3d(ji,jk-1) + dz3d(ji,jk) ) / 2.

             aplus   = ( ( volw3d(ji,jk  ) / dz3d(ji,jk  ) ) + &
                &        ( volw3d(ji,jk+1) / dz3d(ji,jk+1) ) ) / 2.
             dxplus  = ( dz3d(ji,jk) + dz3d(ji,jk+1) ) / 2
             !
             rminus  = ( dtsed / volw3d(ji,jk) ) * diff(jk-1) * aminus / dxminus
             rplus   = ( dtsed / volw3d(ji,jk) ) * diff(jk)   * aplus / dxplus
             !     
             za(ji,jk) = -rminus
             zb(ji,jk) = 1. + rminus + rplus 
             zc(ji,jk) = -rplus
          END DO
       END DO

       DO ji = 1, ndim
          aminus  = ( ( volw3d(ji,nlev-1) / dz3d(ji,nlev-1) ) + &
             &        ( volw3d(ji,nlev)  / dz3d(ji,nlev) ) ) / 2.
          dxminus = ( dz3d(ji,nlev-1) + dz3d(ji,nlev) ) / 2.
          rminus  = ( dtsed / volw3d(ji,nlev) ) * diff(nlev-1) * aminus / dxminus
          !
          za(ji,nlev) = -rminus
          zb(ji,nlev) = 1. + rminus
          zc(ji,nlev) = 0.
       END DO


       ! solves tridiagonal system of linear equations 
       ! -----------------------------------------------

       DO jn = 1, nvar

          zr  (:,:)    = psol(:,:,jn)
          zbet(:  )    = zb(:,1) + preac(:,1,jn)
          psol(:,1,jn) = zr(:,1) / zbet(:)
          ! 
          DO jk = 2, nlev
             DO ji = 1, ndim
                zgamm(ji,jk)   =  zc(ji,jk-1) / zbet(ji)
                zbet(ji)       = ( zb(ji,jk) + preac(ji,jk,jn) ) - za(ji,jk) * zgamm(ji,jk)
                psol(ji,jk,jn) = ( zr(ji,jk) - za(ji,jk) * psol(ji,jk-1,jn) ) / zbet(ji)
             END DO
          ENDDO
          ! 
          DO jk = nlev - 1, 1, -1
             DO ji = 1,ndim
                psol(ji,jk,jn) = psol(ji,jk,jn) - zgamm(ji,jk+1) * psol(ji,jk+1,jn)
             END DO
          ENDDO

       ENDDO


    END SUBROUTINE sed_mat_dsr

    
    SUBROUTINE sed_mat_btb( nvar, ndim, nlev, psol )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE sed_mat_btb  ***
       !!
       !! ** Purpose :  solves tridiagonal system of linear equations 
       !!
       !! ** Method  : 
       !!        1 - computes left hand side of linear system of equations
       !!            for dissolution reaction
       !!
       !!         2 - forward/backward substitution. 
       !!
       !!   History :
       !!        !  04-10 (N. Emprin, M. Gehlen ) original
       !!        !  06-04 (C. Ethe)  Module Re-organization
       !!----------------------------------------------------------------------
       !! * Arguments
       INTEGER , INTENT(in) :: &
          nvar , &  ! number of variables
          ndim , &  ! number of points
          nlev      ! number of sediment levels

      REAL(wp), DIMENSION(ndim,nlev,nvar), INTENT(inout) :: &
          psol      ! solution

       !---Local declarations
       INTEGER  ::  &
          ji, jk, jn

       REAL(wp) ::  &
          aplus,aminus   ,  & 
          rplus,rminus   ,  &
          dxplus,dxminus 

       REAL(wp), DIMENSION(nlev)      :: za, zb, zc
       REAL(wp), DIMENSION(ndim,nlev) :: zr
       REAL(wp), DIMENSION(nmax)      :: zgamm 
       REAL(wp) ::  zbet

 
       !----------------------------------------------------------------------

       ! Computation left hand side of linear system of 
       ! equations for dissolution reaction
       !---------------------------------------------


       ! first sediment level          
       aplus  = ( ( vols(2) / dz(2) ) + ( vols(3) / dz(3) ) ) / 2.
       dxplus = ( dz(2) + dz(3) ) / 2.
       rplus  = ( dtsed / vols(2) ) * db * aplus / dxplus 

       za(1) = 0.
       zb(1) = 1. + rplus
       zc(1) = -rplus

             
       DO jk = 2, nlev - 1
          aminus  = ( ( vols(jk) / dz(jk) ) + ( vols(jk+1) / dz(jk+1) ) ) / 2.
          dxminus = ( dz(jk) + dz(jk+1) ) / 2.
          rminus  = ( dtsed / vols(jk+1) ) * db * aminus / dxminus
          !
          aplus   = ( ( vols(jk+1) / dz(jk+1  ) ) + ( vols(jk+2) / dz(jk+2) ) ) / 2.
          dxplus  = ( dz(jk+1) + dz(jk+2) ) / 2.
          rplus   = ( dtsed / vols(jk+1) ) * db * aplus / dxplus
          !     
          za(jk) = -rminus
          zb(jk) = 1. + rminus + rplus 
          zc(jk) = -rplus
       ENDDO
 
       aminus  = ( ( vols(nlev) / dz(nlev) ) + ( vols(nlev+1) / dz(nlev+1) ) ) / 2.
       dxminus = ( dz(nlev) + dz(nlev+1) ) / 2.
       rminus  = ( dtsed / vols(nlev+1) ) * db * aminus / dxminus
       !
 
       za(nlev) = -rminus
       zb(nlev) = 1. + rminus
       zc(nlev) = 0.


       ! solves tridiagonal system of linear equations 
       ! -----------------------------------------------    
       DO jn = 1, nvar
          
          zr  (:,:)    = psol(:,:,jn)
          zbet         = zb(1)
          psol(:,1,jn) = zr(:,1) / zbet
          ! 
          DO jk = 2, nlev
             zgamm(jk) =  zc(jk-1) / zbet
             zbet      =  zb(jk) - za(jk) * zgamm(jk)
             DO ji = 1, ndim
                psol(ji,jk,jn) = ( zr(ji,jk) - za(jk) * psol(ji,jk-1,jn) ) / zbet
             END DO
          ENDDO
          ! 
          DO jk = nlev - 1, 1, -1
             DO ji = 1,ndim
                psol(ji,jk,jn) = psol(ji,jk,jn) - zgamm(jk+1) * psol(ji,jk+1,jn)
             END DO
          ENDDO

       ENDDO

       
    END SUBROUTINE sed_mat_btb


#else
   !!======================================================================
   !! MODULE sedmat  :   Dummy module 
   !!======================================================================
   !! $Id: sedmat.F90 2355 2015-05-20 07:11:50Z ufla $
CONTAINS
   SUBROUTINE sed_mat         ! Empty routine
   END SUBROUTINE sed_mat
   !!======================================================================
#endif

 END MODULE sedmat