source: palm/trunk/SOURCE/init_advec.f90 @ 623

 SUBROUTINE init_advec

!------------------------------------------------------------------------------!
! Current revisions:
! -----------------
! 
!
! Former revisions:
! -----------------
! $Id: init_advec.f90 484 2010-02-05 07:36:54Z raasch $
! RCS Log replace by Id keyword, revision history cleaned up
!
! Revision 1.6  2004/04/30 11:59:31  raasch
! impulse_advec renamed momentum_advec
!
! Revision 1.1  1999/02/05 09:07:38  raasch
! Initial revision
!
!
! Description:
! ------------
! Initialize constant coefficients and parameters for certain advection schemes.
!------------------------------------------------------------------------------!

    USE advection
    USE arrays_3d
    USE indices
    USE control_parameters

    IMPLICIT NONE

    INTEGER :: i, intervals, j, k
    REAL    :: delt, dn, dnneu, ex1, ex2, ex3, ex4, ex5, ex6, spl_alpha, &
               spl_beta, sterm
    REAL, DIMENSION(:), ALLOCATABLE ::  spl_u, temp


    IF ( scalar_advec == 'bc-scheme' )  THEN

!
!--    Compute exponential coefficients for the Bott-Chlond scheme
       intervals = 1000
       ALLOCATE( aex(intervals), bex(intervals), dex(intervals), eex(intervals) )

       delt  = 1.0 / REAL( intervals )
       sterm = delt * 0.5

       DO  i = 1, intervals

          IF ( sterm > 0.5 )  THEN
             dn = -5.0
          ELSE
             dn = 5.0
          ENDIF

          DO  j = 1, 15
             ex1 = dn * EXP( -dn ) - EXP( 0.5 * dn ) + EXP( -0.5 * dn )
             ex2 = EXP( dn ) - EXP( -dn )
             ex3 = EXP( -dn ) * ( 1.0 - dn ) - 0.5 * EXP(  0.5 * dn ) &
                                             - 0.5 * EXP( -0.5 * dn )
             ex4 = EXP( dn ) + EXP( -dn )
             ex5 = dn * sterm + ex1 / ex2
             ex6 = sterm + ( ex3 * ex2 - ex4 * ex1 ) / ( ex2 * ex2 )
             dnneu = dn - ex5 / ex6
             dn  = dnneu
          ENDDO

          IF ( sterm < 0.5 )  dn = MAX(  2.95E-2, dn )
          IF ( sterm > 0.5 )  dn = MIN( -2.95E-2, dn )
          ex1 = EXP( -dn )
          ex2 = EXP( dn ) - ex1
          aex(i) = -ex1 / ex2
          bex(i) = 1.0 / ex2
          dex(i) = dn
          eex(i) = EXP( dex(i) * 0.5 )
          sterm = sterm + delt

       ENDDO

    ENDIF

    IF ( momentum_advec == 'ups-scheme' .OR. scalar_advec  == 'ups-scheme' )  &
    THEN

!
!--    Provide the constant parameters for the Upstream-Spline advection scheme.
!--    In x- und y-direction the Sherman-Morrison formula is applied
!--    (cf. Press et al, 1986 (Numerical Recipes)).
!
!--    Allocate nonlocal arrays
       ALLOCATE( spl_z_x(0:nx), spl_z_y(0:ny), spl_tri_x(5,0:nx), &
                 spl_tri_y(5,0:ny), spl_tri_zu(5,nzb:nzt+1),      &
                 spl_tri_zw(5,nzb:nzt+1) )

!
!--    Provide diagonal elements of the tridiagonal matrices for all 
!--    directions
       spl_tri_x(1,:)  = 2.0
       spl_tri_y(1,:)  = 2.0
       spl_tri_zu(1,:) = 2.0
       spl_tri_zw(1,:) = 2.0

!
!--    Elements of the cyclic tridiagonal matrix
!--    (same for all horizontal directions)
       spl_alpha = 0.5
       spl_beta  = 0.5

!
!--    Sub- and superdiagonal elements, x-direction
       spl_tri_x(2,0:nx) = 0.5
       spl_tri_x(3,0:nx) = 0.5

!
!--    mMdify the diagonal elements (Sherman-Morrison)
       spl_gamma_x    = -spl_tri_x(1,0)
       spl_tri_x(1,0)  = spl_tri_x(1,0) - spl_gamma_x
       spl_tri_x(1,nx) = spl_tri_x(1,nx) - spl_alpha * spl_beta / spl_gamma_x

!
!--    Split the tridiagonal matrix for Thomas algorithm
       spl_tri_x(4,0) = spl_tri_x(1,0)
       DO  i = 1, nx
          spl_tri_x(5,i) = spl_tri_x(2,i) / spl_tri_x(4,i-1)
          spl_tri_x(4,i) = spl_tri_x(1,i) - spl_tri_x(5,i) * spl_tri_x(3,i-1)
       ENDDO

!
!--    Allocate arrays required locally
       ALLOCATE( temp(0:nx), spl_u(0:nx) )

!
!--    Provide "corrective vector", x-direction
       spl_u(0)      = spl_gamma_x
       spl_u(1:nx-1) = 0.0
       spl_u(nx)     = spl_alpha

!
!--    Solve the system of equations for the corrective vector
!--    (Sherman-Morrison)
       temp(0) = spl_u(0)
       DO  i = 1, nx
          temp(i) = spl_u(i) - spl_tri_x(5,i) * temp(i-1)
       ENDDO
       spl_z_x(nx) = temp(nx) / spl_tri_x(4,nx)
       DO  i = nx-1, 0, -1
          spl_z_x(i) = ( temp(i) - spl_tri_x(3,i) * spl_z_x(i+1) ) / &
                       spl_tri_x(4,i)
       ENDDO

!
!--    Deallocate local arrays, for they are allocated in a different way for
!--    operations in y-direction
       DEALLOCATE( temp, spl_u )    
    
!
!--    Provide sub- and superdiagonal elements, y-direction
       spl_tri_y(2,0:ny) = 0.5
       spl_tri_y(3,0:ny) = 0.5

!
!--    Modify the diagonal elements (Sherman-Morrison)
       spl_gamma_y    = -spl_tri_y(1,0) 
       spl_tri_y(1,0)  = spl_tri_y(1,0) - spl_gamma_y
       spl_tri_y(1,ny) = spl_tri_y(1,ny) - spl_alpha * spl_beta / spl_gamma_y

!
!--    Split the tridiagonal matrix for Thomas algorithm
       spl_tri_y(4,0) = spl_tri_y(1,0)
       DO  j = 1, ny
          spl_tri_y(5,j) = spl_tri_y(2,j) / spl_tri_y(4,j-1)
          spl_tri_y(4,j) = spl_tri_y(1,j) - spl_tri_y(5,j) * spl_tri_y(3,j-1)
       ENDDO

!
!--    Allocate arrays required locally
       ALLOCATE( temp(0:ny), spl_u(0:ny) )

!
!--    Provide "corrective vector", y-direction
       spl_u(0)      = spl_gamma_y
       spl_u(1:ny-1) = 0.0
       spl_u(ny)     = spl_alpha
    
!
!--    Solve the system of equations for the corrective vector
!--    (Sherman-Morrison)
       temp = 0.0
       spl_z_y = 0.0
       temp(0) = spl_u(0)
       DO  j = 1, ny
          temp(j) = spl_u(j) - spl_tri_y(5,j) * temp(j-1)
       ENDDO
       spl_z_y(ny) = temp(ny) / spl_tri_y(4,ny)
       DO  j = ny-1, 0, -1
          spl_z_y(j) = ( temp(j) - spl_tri_y(3,j) * spl_z_y(j+1) ) / &
                       spl_tri_y(4,j)
       ENDDO

!
!--    deallocate local arrays, for they are no longer required
       DEALLOCATE( temp, spl_u )

!
!--    provide sub- and superdiagonal elements, z-direction
       spl_tri_zu(2,nzb)   = 0.0
       spl_tri_zu(2,nzt+1) = 1.0
       spl_tri_zw(2,nzb)   = 0.0
       spl_tri_zw(2,nzt+1) = 1.0

       spl_tri_zu(3,nzb)   = 1.0
       spl_tri_zu(3,nzt+1) = 0.0
       spl_tri_zw(3,nzb)   = 1.0
       spl_tri_zw(3,nzt+1) = 0.0

       DO  k = nzb+1, nzt
          spl_tri_zu(2,k) = dzu(k) / ( dzu(k) + dzu(k+1) )
          spl_tri_zw(2,k) = dzw(k) / ( dzw(k) + dzw(k+1) )
          spl_tri_zu(3,k) = 1.0 - spl_tri_zu(2,k)
          spl_tri_zw(3,k) = 1.0 - spl_tri_zw(2,k)
       ENDDO
    
       spl_tri_zu(4,nzb) = spl_tri_zu(1,nzb)
       spl_tri_zw(4,nzb) = spl_tri_zw(1,nzb)
       DO  k = nzb+1, nzt+1
          spl_tri_zu(5,k) = spl_tri_zu(2,k) / spl_tri_zu(4,k-1)
          spl_tri_zw(5,k) = spl_tri_zw(2,k) / spl_tri_zw(4,k-1)
          spl_tri_zu(4,k) = spl_tri_zu(1,k) - spl_tri_zu(5,k) * &
                            spl_tri_zu(3,k-1)
          spl_tri_zw(4,k) = spl_tri_zw(1,k) - spl_tri_zw(5,k) * &
                            spl_tri_zw(3,k-1)
       ENDDO

    ENDIF

 END SUBROUTINE init_advec