SUBROUTINE init_advec !------------------------------------------------------------------------------! ! Current revisions: ! ----------------- ! ! ! Former revisions: ! ----------------- ! $Id: init_advec.f90 484 2010-02-05 07:36:54Z maronga $ ! 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