Index: palm/trunk/SOURCE/init_advec.f90
===================================================================
--- palm/trunk/SOURCE/init_advec.f90	(revision 4901)
+++ palm/trunk/SOURCE/init_advec.f90	(revision 1)
@@ -1,34 +1,27 @@
-!> @file init_advec.f90
-!--------------------------------------------------------------------------------------------------!
-! This file is part of the PALM model system.
-!
-! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
-! Public License as published by the Free Software Foundation, either version 3 of the License, or
-! (at your option) any later version.
-!
-! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
-! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
-! Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along with PALM. If not, see
-! <http://www.gnu.org/licenses/>.
-!
-! Copyright 1997-2021 Leibniz Universitaet Hannover
-!--------------------------------------------------------------------------------------------------!
-!
-! Current revisions:
+ SUBROUTINE init_advec
+
+!------------------------------------------------------------------------------!
+! Actual revisions:
 ! -----------------
-!
+! 
 !
 ! Former revisions:
 ! -----------------
-! $Id$
-! file re-formatted to follow the PALM coding standard
-!
-! 4360 2020-01-07 11:25:50Z suehring
-! Corrected "Former revisions" section
-!
-! 3655 2019-01-07 16:51:22Z knoop
-! Corrected "Former revisions" section
+! $Log: init_advec.f90,v $
+! Revision 1.6  2004/04/30 11:59:31  raasch
+! impulse_advec renamed momentum_advec
+!
+! Revision 1.5  2001/03/30 07:25:53  raasch
+! Translation of remaining German identifiers (variables, subroutines, etc.)
+!
+! Revision 1.4  2001/01/22 06:57:43  raasch
+! Module test_variables removed
+!
+! Revision 1.3  2000/12/21 16:12:52  letzel
+! All comments translated into English.
+!
+! Revision 1.2  1999/02/17 09:18:45  raasch
+! Fehlerkorrektur, Verfahren muessen unabhaengig voneinander initialisierbar
+! sein
 !
 ! Revision 1.1  1999/02/05 09:07:38  raasch
@@ -38,33 +31,18 @@
 ! Description:
 ! ------------
-!> Initialize constant coefficients and parameters for certain advection schemes.
-!--------------------------------------------------------------------------------------------------!
- SUBROUTINE init_advec
-
-
-    USE advection,                                                                                 &
-        ONLY:  aex, bex, dex, eex
-
-    USE kinds
-
-    USE control_parameters,                                                                        &
-        ONLY:  scalar_advec
+! Initialize constant coefficients and parameters for certain advection schemes.
+!------------------------------------------------------------------------------!
+
+    USE advection
+    USE arrays_3d
+    USE indices
+    USE control_parameters
 
     IMPLICIT NONE
 
-    INTEGER(iwp) ::  i          !<
-    INTEGER(iwp) ::  intervals  !<
-    INTEGER(iwp) ::  j          !<
-
-    REAL(wp) :: delt   !<
-    REAL(wp) :: dn     !<
-    REAL(wp) :: dnneu  !<
-    REAL(wp) :: ex1    !<
-    REAL(wp) :: ex2    !<
-    REAL(wp) :: ex3    !<
-    REAL(wp) :: ex4    !<
-    REAL(wp) :: ex5    !<
-    REAL(wp) :: ex6    !<
-    REAL(wp) :: sterm  !<
+    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
 
 
@@ -76,20 +54,20 @@
        ALLOCATE( aex(intervals), bex(intervals), dex(intervals), eex(intervals) )
 
-       delt  = 1.0_wp / REAL( intervals, KIND=wp )
-       sterm = delt * 0.5_wp
+       delt  = 1.0 / REAL( intervals )
+       sterm = delt * 0.5
 
        DO  i = 1, intervals
 
-          IF ( sterm > 0.5_wp )  THEN
-             dn = -5.0_wp
+          IF ( sterm > 0.5 )  THEN
+             dn = -5.0
           ELSE
-             dn = 5.0_wp
+             dn = 5.0
           ENDIF
 
           DO  j = 1, 15
-             ex1 = dn * EXP( -dn ) - EXP( 0.5_wp * dn ) + EXP( -0.5_wp * dn )
+             ex1 = dn * EXP( -dn ) - EXP( 0.5 * dn ) + EXP( -0.5 * dn )
              ex2 = EXP( dn ) - EXP( -dn )
-             ex3 = EXP( -dn ) * ( 1.0_wp - dn ) - 0.5_wp * EXP(  0.5_wp * dn )                     &
-                                                - 0.5_wp * EXP( -0.5_wp * 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
@@ -99,12 +77,12 @@
           ENDDO
 
-          IF ( sterm < 0.5_wp )  dn = MAX(  2.95E-2_wp, dn )
-          IF ( sterm > 0.5_wp )  dn = MIN( -2.95E-2_wp, dn )
+          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_wp / ex2
+          bex(i) = 1.0 / ex2
           dex(i) = dn
-          eex(i) = EXP( dex(i) * 0.5_wp )
+          eex(i) = EXP( dex(i) * 0.5 )
           sterm = sterm + delt
 
@@ -113,3 +91,157 @@
     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
