SUBROUTINE spline_x( vad_in_out, ad_v, var_char ) !------------------------------------------------------------------------------! ! Actual revisions: ! ----------------- ! ! ! Former revisions: ! ----------------- ! $Id: spline_x.f90 4 2007-02-13 11:33:16Z heinze $ ! RCS Log replace by Id keyword, revision history cleaned up ! ! Revision 1.8 2004/04/30 12:54:20 raasch ! Names of transpose indices changed, enlarged transposition arrays introduced ! ! Revision 1.1 1999/02/05 09:15:59 raasch ! Initial revision ! ! ! Description: ! ------------ ! Upstream-spline advection along x ! ! Input/output parameters: ! ad_v = advecting wind speed component ! vad_in_out = quantity to be advected, excluding ghost- or cyclic boundaries ! result is given to the calling routine in this array ! var_char = string which defines the quantity to be advected ! ! Internal arrays: ! r = 2D-working array (right hand side of linear equation, buffer for ! Long filter) ! tf = tendency field (2D), used for long filter ! vad = quantity to be advected (2D), including ghost- or cyclic ! boundarys along the direction of advection ! wrk_long = working array (long coefficients) ! wrk_spline = working array (spline coefficients) !------------------------------------------------------------------------------! USE advection USE grid_variables USE indices USE statistics USE control_parameters USE transpose_indices IMPLICIT NONE CHARACTER (LEN=*) :: var_char INTEGER :: component, i, j, k, sr REAL :: overshoot_limit, sm_faktor, t1, t2, t3, ups_limit REAL, DIMENSION(:,:), ALLOCATABLE :: r, tf, vad, wrk_spline REAL, DIMENSION(:,:,:), ALLOCATABLE :: wrk_long #if defined( __parallel ) REAL :: ad_v(0:nxa,nys_x:nyn_xa,nzb_x:nzt_xa), & vad_in_out(0:nxa,nys_x:nyn_xa,nzb_x:nzt_xa) #else REAL :: ad_v(nzb+1:nzt,nys:nyn,nxl:nxr), & vad_in_out(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) #endif ! !-- Set criteria for switching between upstream- and upstream-spline-method IF ( var_char == 'u' ) THEN overshoot_limit = overshoot_limit_u ups_limit = ups_limit_u component = 1 ELSEIF ( var_char == 'v' ) THEN overshoot_limit = overshoot_limit_v ups_limit = ups_limit_v component = 2 ELSEIF ( var_char == 'w' ) THEN overshoot_limit = overshoot_limit_w ups_limit = ups_limit_w component = 3 ELSEIF ( var_char == 'pt' ) THEN overshoot_limit = overshoot_limit_pt ups_limit = ups_limit_pt component = 4 ELSEIF ( var_char == 'e' ) THEN overshoot_limit = overshoot_limit_e ups_limit = ups_limit_e component = 5 ENDIF ! !-- Initialize calculation of relative upstream fraction sums_up_fraction_l(component,1,:) = 0.0 #if defined( __parallel ) ! !-- Allocate working arrays ALLOCATE( r(-1:nx+1,nys_x:nyn_x), vad(-1:nx+1,nys_x:nyn_x), & wrk_spline(0:nx,nys_x:nyn_x) ) IF ( long_filter_factor /= 0.0 ) THEN ALLOCATE( tf(0:nx,nys_x:nyn_x), wrk_long(0:nx,nys_x:nyn_x,1:3) ) ENDIF ! !-- Loop over all gridpoints along z DO k = nzb_x, nzt_x ! !-- Store array to be advected on work array and add cyclic boundary along x vad(0:nx,nys_x:nyn_x) = vad_in_out(0:nx,nys_x:nyn_x,k) vad(-1,:) = vad(nx,:) vad(nx+1,:) = vad(0,:) ! !-- Calculate right hand side DO j = nys_x, nyn_x DO i = 0, nx r(i,j) = 3.0 * ( & spl_tri_x(2,i) * ( vad(i,j) - vad(i-1,j) ) * ddx + & spl_tri_x(3,i) * ( vad(i+1,j) - vad(i,j) ) * ddx & ) ENDDO ENDDO ! !-- Forward substitution DO j = nys_x, nyn_x wrk_spline(0,j) = r(0,j) DO i = 1, nx wrk_spline(i,j) = r(i,j) - spl_tri_x(5,i) * wrk_spline(i-1,j) ENDDO ENDDO ! !-- Backward substitution (Sherman-Morrison-formula) DO j = nys_x, nyn_x r(nx,j) = wrk_spline(nx,j) / spl_tri_x(4,nx) DO i = nx-1, 0, -1 r(i,j) = ( wrk_spline(i,j) - spl_tri_x(3,i) * r(i+1,j) ) / & spl_tri_x(4,i) ENDDO sm_faktor = ( r(0,j) + 0.5 * r(nx,j) / spl_gamma_x ) / & ( 1.0 + spl_z_x(0) + 0.5 * spl_z_x(nx) / spl_gamma_x ) DO i = 0, nx r(i,j) = r(i,j) - sm_faktor * spl_z_x(i) ENDDO ENDDO ! !-- Add cyclic boundary to right hand side r(-1,:) = r(nx,:) r(nx+1,:) = r(0,:) ! !-- Calculate advection along x DO j = nys_x, nyn_x DO i = 0, nx IF ( ad_v(i,j,k) == 0.0 ) THEN vad_in_out(i,j,k) = vad(i,j) ELSEIF ( ad_v(i,j,k) > 0.0 ) THEN IF ( ABS( vad(i,j) - vad(i-1,j) ) <= ups_limit ) THEN vad_in_out(i,j,k) = vad(i,j) - dt_3d * ad_v(i,j,k) * & ( vad(i,j) - vad(i-1,j) ) * ddx ! !-- Calculate upstream fraction in % (s. flow_statistics) DO sr = 0, statistic_regions sums_up_fraction_l(component,1,sr) = & sums_up_fraction_l(component,1,sr) + 1.0 ENDDO ELSE t1 = ad_v(i,j,k) * dt_3d * ddx t2 = 3.0 * ( vad(i-1,j) - vad(i,j) ) + & ( 2.0 * r(i,j) + r(i-1,j) ) * dx t3 = 2.0 * ( vad(i-1,j) - vad(i,j) ) + & ( r(i,j) + r(i-1,j) ) * dx vad_in_out(i,j,k) = vad(i,j) - r(i,j) * t1 * dx + & t2 * t1**2 - t3 * t1**3 IF ( vad(i-1,j) == vad(i,j) ) THEN vad_in_out(i,j,k) = vad(i,j) ENDIF ENDIF ELSE IF ( ABS( vad(i,j) - vad(i+1,j) ) <= ups_limit ) THEN vad_in_out(i,j,k) = vad(i,j) - dt_3d * ad_v(i,j,k) * & ( vad(i+1,j) - vad(i,j) ) * ddx ! !-- Calculate upstream fraction in % (s. flow_statistics) DO sr = 0, statistic_regions sums_up_fraction_l(component,1,sr) = & sums_up_fraction_l(component,1,sr) + 1.0 ENDDO ELSE t1 = -ad_v(i,j,k) * dt_3d * ddx t2 = 3.0 * ( vad(i,j) - vad(i+1,j) ) + & ( 2.0 * r(i,j) + r(i+1,j) ) * dx t3 = 2.0 * ( vad(i,j) - vad(i+1,j) ) + & ( r(i,j) + r(i+1,j) ) * dx vad_in_out(i,j,k) = vad(i,j) + r(i,j) * t1 * dx - & t2 * t1**2 + t3 * t1**3 IF ( vad(i+1,j) == vad(i,j) ) THEN vad_in_out(i,j,k) = vad(i,j) ENDIF ENDIF ENDIF ENDDO ENDDO ! !-- Limit values in order to prevent overshooting IF ( cut_spline_overshoot ) THEN DO j = nys_x, nyn_x DO i = 0, nx IF ( ad_v(i,j,k) > 0.0 ) THEN IF ( vad(i,j) > vad(i-1,j) ) THEN vad_in_out(i,j,k) = MIN( vad_in_out(i,j,k), & vad(i,j) + overshoot_limit ) vad_in_out(i,j,k) = MAX( vad_in_out(i,j,k), & vad(i-1,j) - overshoot_limit ) ELSE vad_in_out(i,j,k) = MAX( vad_in_out(i,j,k), & vad(i,j) - overshoot_limit ) vad_in_out(i,j,k) = MIN( vad_in_out(i,j,k), & vad(i-1,j) + overshoot_limit ) ENDIF ELSE IF ( vad(i,j) > vad(i+1,j) ) THEN vad_in_out(i,j,k) = MIN( vad_in_out(i,j,k), & vad(i,j) + overshoot_limit ) vad_in_out(i,j,k) = MAX( vad_in_out(i,j,k), & vad(i+1,j) - overshoot_limit ) ELSE vad_in_out(i,j,k) = MAX( vad_in_out(i,j,k), & vad(i,j) - overshoot_limit ) vad_in_out(i,j,k) = MIN( vad_in_out(i,j,k), & vad(i+1,j) + overshoot_limit ) ENDIF ENDIF ENDDO ENDDO ENDIF ! !-- Long-filter (acting on tendency only) IF ( long_filter_factor /= 0.0 ) THEN ! !-- Compute tendency DO j = nys_x, nyn_x DO i = 0, nx tf(i,j) = vad_in_out(i,j,k) - vad(i,j) ENDDO ENDDO ! !-- Apply the filter. DO j = nys_x, nyn_x wrk_long(0,j,1) = 2.0 * ( 1.0 + long_filter_factor ) wrk_long(0,j,2) = ( 1.0 - long_filter_factor ) / wrk_long(0,j,1) wrk_long(0,j,3) = ( long_filter_factor * tf(nx,j) + & 2.0 * tf(0,j) + tf(1,j) & ) / wrk_long(0,j,1) DO i = 1, nx-1 wrk_long(i,j,1) = 2.0 * ( 1.0 + long_filter_factor ) - & ( 1.0 - long_filter_factor ) * wrk_long(i-1,j,2) wrk_long(i,j,2) = ( 1.0 - long_filter_factor ) / wrk_long(i,j,1) wrk_long(i,j,3) = ( tf(i-1,j) + 2.0 * tf(i,j) + & tf(i+1,j) - ( 1.0 - long_filter_factor ) * & wrk_long(i-1,j,3) ) / wrk_long(i,j,1) ENDDO wrk_long(nx,j,1) = 2.0 * ( 1.0 + long_filter_factor ) - & ( 1.0 - long_filter_factor ) * wrk_long(nx-1,j,2) wrk_long(nx,j,2) = ( 1.0 - long_filter_factor ) / wrk_long(nx,j,1) wrk_long(nx,j,3) = ( tf(nx-1,j) + 2.0 * tf(nx,j) + & long_filter_factor * tf(0,j) - & ( 1.0 - long_filter_factor ) * & wrk_long(nx-1,j,3) & ) / wrk_long(nx,j,1) r(nx,j) = wrk_long(nx,j,3) ENDDO DO i = nx-1, 0, -1 DO j = nys_x, nyn_x r(i,j) = wrk_long(i,j,3) - wrk_long(i,j,2) * r(i+1,j) ENDDO ENDDO DO j = nys_x, nyn_x DO i = 0, nx vad_in_out(i,j,k) = vad(i,j) + r(i,j) ENDDO ENDDO ENDIF ! Long filter ENDDO #else ! !-- Allocate working arrays ALLOCATE( r(nzb+1:nzt,nxl-1:nxr+1), vad(nzb:nzt+1,nxl-1:nxr+1), & wrk_spline(nzb+1:nzt,nxl-1:nxr+1) ) IF ( long_filter_factor /= 0.0 ) THEN ALLOCATE( tf(nzb+1:nzt,nxl-1:nxr+1), wrk_long(nzb+1:nzt,0:nx,1:3) ) ENDIF ! !-- Loop over all gridpoints along y DO j = nys, nyn ! !-- Store array to be advected on work array and add cyclic boundary along x vad(:,:) = vad_in_out(:,j,:) vad(:,-1) = vad(:,nx) vad(:,nx+1) = vad(:,0) ! !-- Calculate right hand side DO i = 0, nx DO k = nzb+1, nzt r(k,i) = 3.0 * ( & spl_tri_x(2,i) * ( vad(k,i) - vad(k,i-1) ) * ddx + & spl_tri_x(3,i) * ( vad(k,i+1) - vad(k,i) ) * ddx & ) ENDDO ENDDO ! !-- Forward substitution DO k = nzb+1, nzt wrk_spline(k,0) = r(k,0) ENDDO DO i = 1, nx DO k = nzb+1, nzt wrk_spline(k,i) = r(k,i) - spl_tri_x(5,i) * wrk_spline(k,i-1) ENDDO ENDDO ! !-- Backward substitution (Sherman-Morrison-formula) DO k = nzb+1, nzt r(k,nx) = wrk_spline(k,nx) / spl_tri_x(4,nx) ENDDO DO k = nzb+1, nzt DO i = nx-1, 0, -1 r(k,i) = ( wrk_spline(k,i) - spl_tri_x(3,i) * r(k,i+1) ) / & spl_tri_x(4,i) ENDDO sm_faktor = ( r(k,0) + 0.5 * r(k,nx) / spl_gamma_x ) / & ( 1.0 + spl_z_x(0) + 0.5 * spl_z_x(nx) / spl_gamma_x ) DO i = 0, nx r(k,i) = r(k,i) - sm_faktor * spl_z_x(i) ENDDO ENDDO ! !-- Add cyclic boundary to the right hand side r(:,-1) = r(:,nx) r(:,nx+1) = r(:,0) ! !-- Calculate advection along x DO i = 0, nx DO k = nzb+1, nzt IF (ad_v(k,j,i) == 0.0) THEN vad_in_out(k,j,i) = vad(k,i) ELSEIF ( ad_v(k,j,i) > 0.0) THEN IF ( ABS( vad(k,i) - vad(k,i-1) ) <= ups_limit ) THEN vad_in_out(k,j,i) = vad(k,i) - dt_3d * ad_v(k,j,i) * & ( vad(k,i) - vad(k,i-1) ) * ddx ! !-- Calculate upstream fraction in % (s. flow_statistics) DO sr = 0, statistic_regions sums_up_fraction_l(component,1,sr) = & sums_up_fraction_l(component,1,sr) + 1.0 ENDDO ELSE t1 = ad_v(k,j,i) * dt_3d * ddx t2 = 3.0 * ( vad(k,i-1) - vad(k,i) ) + & ( 2.0 * r(k,i) + r(k,i-1) ) * dx t3 = 2.0 * ( vad(k,i-1) - vad(k,i) ) + & ( r(k,i) + r(k,i-1) ) * dx vad_in_out(k,j,i) = vad(k,i) - r(k,i) * t1 * dx + & t2 * t1**2 - t3 * t1**3 IF ( vad(k,i-1) == vad(k,i) ) THEN vad_in_out(k,j,i) = vad(k,i) ENDIF ENDIF ELSE IF ( ABS( vad(k,i) - vad(k,i+1) ) <= ups_limit ) THEN vad_in_out(k,j,i) = vad(k,i) - dt_3d * ad_v(k,j,i) * & ( vad(k,i+1) - vad(k,i) ) * ddx ! !-- Calculate upstream fraction in % (s. flow_statistics) DO sr = 0, statistic_regions sums_up_fraction_l(component,1,sr) = & sums_up_fraction_l(component,1,sr) + 1.0 ENDDO ELSE t1 = -ad_v(k,j,i) * dt_3d * ddx t2 = 3.0 * ( vad(k,i) - vad(k,i+1) ) + & ( 2.0 * r(k,i) + r(k,i+1)) * dx t3 = 2.0 * ( vad(k,i) - vad(k,i+1) ) + & ( r(k,i) + r(k,i+1) ) * dx vad_in_out(k,j,i) = vad(k,i) + r(k,i) * t1 * dx - & t2 * t1**2 + t3 * t1**3 IF ( vad(k,i+1) == vad(k,i) ) THEN vad_in_out(k,j,i) = vad(k,i) ENDIF ENDIF ENDIF ENDDO ENDDO ! !-- Limit values in order to prevent overshooting IF ( cut_spline_overshoot ) THEN DO i = 0, nx DO k = nzb+1, nzt IF ( ad_v(k,j,i) > 0.0 ) THEN IF ( vad(k,i) > vad(k,i-1) ) THEN vad_in_out(k,j,i) = MIN( vad_in_out(k,j,i), & vad(k,i) + overshoot_limit ) vad_in_out(k,j,i) = MAX( vad_in_out(k,j,i), & vad(k,i-1) - overshoot_limit ) ELSE vad_in_out(k,j,i) = MAX( vad_in_out(k,j,i), & vad(k,i) - overshoot_limit ) vad_in_out(k,j,i) = MIN( vad_in_out(k,j,i), & vad(k,i-1) + overshoot_limit ) ENDIF ELSE IF ( vad(k,i) > vad(k,i+1) ) THEN vad_in_out(k,j,i) = MIN( vad_in_out(k,j,i), & vad(k,i) + overshoot_limit ) vad_in_out(k,j,i) = MAX( vad_in_out(k,j,i), & vad(k,i+1) - overshoot_limit ) ELSE vad_in_out(k,j,i) = MAX( vad_in_out(k,j,i), & vad(k,i) - overshoot_limit ) vad_in_out(k,j,i) = MIN( vad_in_out(k,j,i), & vad(k,i+1) + overshoot_limit ) ENDIF ENDIF ENDDO ENDDO ENDIF ! !-- Long filter (acting on tendency only) IF ( long_filter_factor /= 0.0 ) THEN ! !-- Compute tendency DO i = nxl, nxr DO k = nzb+1, nzt tf(k,i) = vad_in_out(k,j,i) - vad(k,i) ENDDO ENDDO ! !-- Apply the filter wrk_long(:,0,1) = 2.0 * ( 1.0 + long_filter_factor ) wrk_long(:,0,2) = ( 1.0 - long_filter_factor ) / wrk_long(:,0,1) wrk_long(:,0,3) = ( long_filter_factor * tf(:,nx) + & 2.0 * tf(:,0) + tf(:,1) ) / wrk_long(:,0,1) DO i = 1, nx-1 DO k = nzb+1, nzt wrk_long(k,i,1) = 2.0 * ( 1.0 + long_filter_factor ) - & ( 1.0 - long_filter_factor ) * & wrk_long(k,i-1,2) wrk_long(k,i,2) = ( 1.0 - long_filter_factor ) / wrk_long(k,i,1) wrk_long(k,i,3) = ( tf(k,i-1) + 2.0 * tf(k,i) + & tf(k,i+1) - ( 1.0 - long_filter_factor ) * & wrk_long(k,i-1,3) ) / wrk_long(k,i,1) ENDDO wrk_long(:,nx,1) = 2.0 * ( 1.0 + long_filter_factor ) - & ( 1.0 - long_filter_factor ) * & wrk_long(:,nx-1,2) wrk_long(:,nx,2) = ( 1.0 - long_filter_factor ) / wrk_long(:,nx,1) wrk_long(:,nx,3) = ( tf(:,nx-1) + 2.0 * tf(:,nx) + & long_filter_factor * tf(:,0) - & ( 1.0 - long_filter_factor ) * & wrk_long(:,nx-1,3) ) / wrk_long(:,nx,1) r(:,nx) = wrk_long(:,nx,3) ENDDO DO i = nx-1, 0, -1 DO k = nzb+1, nzt r(k,i) = wrk_long(k,i,3) - wrk_long(k,i,2) * r(k,i+1) ENDDO ENDDO DO i = 0, nx DO k = nzb+1, nzt vad_in_out(k,j,i) = vad(k,i) + r(k,i) ENDDO ENDDO ENDIF ! Long filter ENDDO #endif IF ( long_filter_factor /= 0.0 ) DEALLOCATE( tf, wrk_long ) DEALLOCATE( r, vad, wrk_spline ) END SUBROUTINE spline_x