SUBROUTINE spline_y( vad_in_out, ad_v, var_char ) !------------------------------------------------------------------------------! ! Current revisions: ! ----------------- ! ! ! Former revisions: ! ----------------- ! $Id: spline_y.f90 484 2010-02-05 07:36:54Z maronga $ ! RCS Log replace by Id keyword, revision history cleaned up ! ! Revision 1.9 2004/04/30 12:54:37 raasch ! Names of transpose indices changed, enlarged transposition arrays introduced ! ! Revision 1.1 1999/02/05 09:16:31 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:nya,nxl_y:nxr_ya,nzb_y:nzt_ya), & vad_in_out(0:nya,nxl_y:nxr_ya,nzb_y:nzt_ya) #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,2,:) = 0.0 #if defined( __parallel ) ! !-- Allocate working arrays ALLOCATE( r(-1:ny+1,nxl_y:nxr_y), & vad(-1:ny+1,nxl_y:nxr_y), & wrk_spline(0:ny,nxl_y:nxr_y) ) IF ( long_filter_factor /= 0.0 ) THEN ALLOCATE( tf(0:ny,nxl_y:nxr_y), & wrk_long(0:ny,nxl_y:nxr_y,1:3) ) ENDIF ! !-- Loop over all gridpoints along z DO k = nzb_y, nzt_y ! !-- Store array to be advected on work array and add cyclic boundary along y vad(0:ny,nxl_y:nxr_y) = vad_in_out(0:ny,nxl_y:nxr_y,k) vad(-1,:) = vad(ny,:) vad(ny+1,:) = vad(0,:) ! !-- Calculate right hand side DO i = nxl_y, nxr_y DO j = 0, ny r(j,i) = 3.0 * ( & spl_tri_y(2,j) * ( vad(j,i) - vad(j-1,i) ) * ddy + & spl_tri_y(3,j) * ( vad(j+1,i) - vad(j,i) ) * ddy & ) ENDDO ENDDO ! !-- Forward substitution DO i = nxl_y, nxr_y wrk_spline(0,i) = r(0,i) DO j = 1, ny wrk_spline(j,i) = r(j,i) - spl_tri_y(5,j) * wrk_spline(j-1,i) ENDDO ENDDO ! !-- Backward substitution (sherman-Morrison-formula) DO i = nxl_y, nxr_y r(ny,i) = wrk_spline(ny,i) / spl_tri_y(4,ny) DO j = ny-1, 0, -1 r(j,i) = ( wrk_spline(j,i) - spl_tri_y(3,j) * r(j+1,i) ) / & spl_tri_y(4,j) ENDDO sm_faktor = ( r(0,i) + 0.5 * r(ny,i) / spl_gamma_y ) / & ( 1.0 + spl_z_y(0) + 0.5 * spl_z_y(ny) / spl_gamma_y ) DO j = 0, ny r(j,i) = r(j,i) - sm_faktor * spl_z_y(j) ENDDO ENDDO ! !-- Add cyclic boundary conditions to right hand side r(-1,:) = r(ny,:) r(ny+1,:) = r(0,:) ! !-- Calculate advection along y DO i = nxl_y, nxr_y DO j = 0, ny IF ( ad_v(j,i,k) == 0.0 ) THEN vad_in_out(j,i,k) = vad(j,i) ELSEIF ( ad_v(j,i,k) > 0.0 ) THEN IF ( ABS( vad(j,i) - vad(j-1,i) ) <= ups_limit ) THEN vad_in_out(j,i,k) = vad(j,i) - dt_3d * ad_v(j,i,k) * & ( vad(j,i) - vad(j-1,i) ) * ddy ! !-- Calculate upstream fraction in % (s. flow_statistics) DO sr = 0, statistic_regions sums_up_fraction_l(component,2,sr) = & sums_up_fraction_l(component,2,sr) + 1.0 ENDDO ELSE t1 = ad_v(j,i,k) * dt_3d * ddy t2 = 3.0 * ( vad(j-1,i) - vad(j,i) ) + & ( 2.0 * r(j,i) + r(j-1,i) ) * dy t3 = 2.0 * ( vad(j-1,i) - vad(j,i) ) + & ( r(j,i) + r(j-1,i) ) * dy vad_in_out(j,i,k) = vad(j,i) - r(j,i) * t1 * dy + & t2 * t1**2 - t3 * t1**3 IF ( vad(j-1,i) == vad(j,i) ) THEN vad_in_out(j,i,k) = vad(j,i) ENDIF ENDIF ELSE IF ( ABS( vad(j,i) - vad(j+1,i) ) <= ups_limit ) THEN vad_in_out(j,i,k) = vad(j,i) - dt_3d * ad_v(j,i,k) * & ( vad(j+1,i) - vad(j,i) ) * ddy ! !-- Calculate upstream fraction in % (s. flow_statistics) DO sr = 0, statistic_regions sums_up_fraction_l(component,2,sr) = & sums_up_fraction_l(component,2,sr) + 1.0 ENDDO ELSE t1 = -ad_v(j,i,k) * dt_3d * ddy t2 = 3.0 * ( vad(j,i) - vad(j+1,i) ) + & ( 2.0 * r(j,i) + r(j+1,i) ) * dy t3 = 2.0 * ( vad(j,i) - vad(j+1,i) ) + & ( r(j,i) + r(j+1,i) ) * dy vad_in_out(j,i,k) = vad(j,i) + r(j,i) * t1 * dy - & t2 * t1**2 + t3 * t1**3 IF ( vad(j+1,i) == vad(j,i) ) THEN vad_in_out(j,i,k) = vad(j,i) ENDIF ENDIF ENDIF ENDDO ENDDO ! !-- Limit values in order to prevent overshooting IF ( cut_spline_overshoot ) THEN DO i = nxl_y, nxr_y DO j = 0, ny IF ( ad_v(j,i,k) > 0.0 ) THEN IF ( vad(j,i) > vad(j-1,i) ) THEN vad_in_out(j,i,k) = MIN( vad_in_out(j,i,k), & vad(j,i) + overshoot_limit ) vad_in_out(j,i,k) = MAX( vad_in_out(j,i,k), & vad(j-1,i) - overshoot_limit ) ELSE vad_in_out(j,i,k) = MAX( vad_in_out(j,i,k), & vad(j,i) - overshoot_limit ) vad_in_out(j,i,k) = MIN( vad_in_out(j,i,k), & vad(j-1,i) + overshoot_limit ) ENDIF ELSE IF ( vad(j,i) > vad(j+1,i) ) THEN vad_in_out(j,i,k) = MIN( vad_in_out(j,i,k), & vad(j,i) + overshoot_limit ) vad_in_out(j,i,k) = MAX( vad_in_out(j,i,k), & vad(j+1,i) - overshoot_limit ) ELSE vad_in_out(j,i,k) = MAX( vad_in_out(j,i,k), & vad(j,i) - overshoot_limit ) vad_in_out(j,i,k) = MIN( vad_in_out(j,i,k), & vad(j+1,i) + overshoot_limit ) ENDIF ENDIF ENDDO ENDDO ENDIF ! !-- Long-filter (acting on tendency only) IF ( long_filter_factor /= 0.0 ) THEN ! !-- Compute tendency. Filter only acts on this quantity. DO i = nxl_y, nxr_y DO j = 0, ny tf(j,i) = vad_in_out(j,i,k) - vad(j,i) ENDDO ENDDO ! !-- Apply the filter. DO i = nxl_y, nxr_y wrk_long(0,i,1) = 2.0 * ( 1.0 + long_filter_factor ) wrk_long(0,i,2) = ( 1.0 - long_filter_factor ) / wrk_long(0,i,1) wrk_long(0,i,3) = ( long_filter_factor * tf(ny,i) + & 2.0 * tf(0,i) + tf(1,i) & ) / wrk_long(0,i,1) DO j = 1, ny-1 wrk_long(j,i,1) = 2.0 * ( 1.0 + long_filter_factor ) - & ( 1.0 - long_filter_factor ) * wrk_long(j-1,i,2) wrk_long(j,i,2) = ( 1.0 - long_filter_factor ) / wrk_long(j,i,1) wrk_long(j,i,3) = ( tf(j-1,i) + 2.0 * tf(j,i) + & tf(j+1,i) - ( 1.0 - long_filter_factor ) * & wrk_long(j-1,i,3) ) / wrk_long(j,i,1) ENDDO wrk_long(ny,i,1) = 2.0 * ( 1.0 + long_filter_factor ) - & ( 1.0 - long_filter_factor ) * wrk_long(ny-1,i,2) wrk_long(ny,i,2) = ( 1.0 - long_filter_factor ) / wrk_long(ny,i,1) wrk_long(ny,i,3) = ( tf(ny-1,i) + 2.0 * tf(ny,i) + & long_filter_factor * tf(0,i) - & ( 1.0 - long_filter_factor ) * & wrk_long(ny-1,i,3) & ) / wrk_long(ny,i,1) r(ny,i) = wrk_long(ny,i,3) ENDDO DO j = ny-1, 0, -1 DO i = nxl_y, nxr_y r(j,i) = wrk_long(j,i,3) - wrk_long(j,i,2) * r(j+1,i) ENDDO ENDDO DO i = nxl_y, nxr_y DO j = 0, ny vad_in_out(j,i,k) = vad(j,i) + r(j,i) ENDDO ENDDO ENDIF ! Long filter ENDDO #else ! !-- Allocate working arrays ALLOCATE( r(nzb+1:nzt,nys-1:nyn+1), vad(nzb:nzt+1,nys-1:nyn+1), & wrk_spline(nzb+1:nzt,nys-1:nyn+1) ) IF ( long_filter_factor /= 0.0 ) THEN ALLOCATE( tf(nzb+1:nzt,nys-1:nyn+1), wrk_long(nzb+1:nzt,0:ny,1:3) ) ENDIF ! !-- Loop over all gridpoints along x DO i = nxl, nxr ! !-- Store array to be advected on work array and add cyclic boundary along x vad(:,:) = vad_in_out(:,:,i) vad(:,-1) = vad(:,ny) vad(:,ny+1) = vad(:,0) ! !-- Calculate right hand side DO j = 0, ny DO k = nzb+1, nzt r(k,j) = 3.0 * ( & spl_tri_y(2,j) * ( vad(k,j) - vad(k,j-1) ) * ddy + & spl_tri_y(3,j) * ( vad(k,j+1) - vad(k,j) ) * ddy & ) ENDDO ENDDO ! !-- Forward substitution DO k = nzb+1, nzt wrk_spline(k,0) = r(k,0) ENDDO DO j = 1, ny DO k = nzb+1, nzt wrk_spline(k,j) = r(k,j) - spl_tri_y(5,j) * wrk_spline(k,j-1) ENDDO ENDDO ! !-- Backward substitution (Sherman-Morrison-formula) DO k = nzb+1, nzt r(k,ny) = wrk_spline(k,ny) / spl_tri_y(4,ny) ENDDO DO k = nzb+1, nzt DO j = ny-1, 0, -1 r(k,j) = ( wrk_spline(k,j) - spl_tri_y(3,j) * r(k,j+1) ) / & spl_tri_y(4,j) ENDDO sm_faktor = ( r(k,0) + 0.5 * r(k,ny) / spl_gamma_y ) / & ( 1.0 + spl_z_y(0) + 0.5 * spl_z_y(ny) / spl_gamma_y ) DO j = 0, ny r(k,j) = r(k,j) - sm_faktor * spl_z_y(j) ENDDO ENDDO ! !-- Add cyclic boundary to the right hand side r(:,-1) = r(:,ny) r(:,ny+1) = r(:,0) ! !-- Calculate advection along y DO j = 0, ny DO k = nzb+1, nzt IF ( ad_v(k,j,i) == 0.0 ) THEN vad_in_out(k,j,i) = vad(k,j) ELSEIF ( ad_v(k,j,i) > 0.0 ) THEN IF ( ABS( vad(k,j) - vad(k,j-1) ) <= ups_limit ) THEN vad_in_out(k,j,i) = vad(k,j) - dt_3d * ad_v(k,j,i) * & ( vad(k,j) - vad(k,j-1) ) * ddy ! !-- Calculate upstream fraction in % (s. flow_statistics) DO sr = 0, statistic_regions sums_up_fraction_l(component,2,sr) = & sums_up_fraction_l(component,2,sr) + 1.0 ENDDO ELSE t1 = ad_v(k,j,i) * dt_3d * ddy t2 = 3.0 * ( vad(k,j-1) - vad(k,j) ) + & ( 2.0 * r(k,j) + r(k,j-1) ) * dy t3 = 2.0 * ( vad(k,j-1) - vad(k,j) ) + & ( r(k,j) + r(k,j-1) ) * dy vad_in_out(k,j,i) = vad(k,j) - r(k,j) * t1 * dy + & t2 * t1**2 - t3 * t1**3 IF ( vad(k,j-1) == vad(k,j) ) THEN vad_in_out(k,j,i) = vad(k,j) ENDIF ENDIF ELSE IF ( ABS( vad(k,j) - vad(k,j+1) ) <= ups_limit ) THEN vad_in_out(k,j,i) = vad(k,j) - dt_3d * ad_v(k,j,i) * & ( vad(k,j+1) - vad(k,j) ) * ddy ! !-- Calculate upstream fraction in % (s. flow_statistics) DO sr = 0, statistic_regions sums_up_fraction_l(component,2,sr) = & sums_up_fraction_l(component,2,sr) + 1.0 ENDDO ELSE t1 = -ad_v(k,j,i) * dt_3d * ddy t2 = 3.0 * ( vad(k,j) - vad(k,j+1) ) + & ( 2.0 * r(k,j) + r(k,j+1) ) * dy t3 = 2.0 * ( vad(k,j) - vad(k,j+1) ) + & ( r(k,j) + r(k,j+1) ) * dy vad_in_out(k,j,i) = vad(k,j) + r(k,j) * t1 * dy - & t2 * t1**2 + t3 * t1**3 IF ( vad(k,j+1) == vad(k,j) ) THEN vad_in_out(k,j,i) = vad(k,j) ENDIF ENDIF ENDIF ENDDO ENDDO ! !-- Limit values in order to prevent overshooting IF ( cut_spline_overshoot ) THEN DO j = 0, ny DO k = nzb+1, nzt IF ( ad_v(k,j,i) > 0.0 ) THEN IF ( vad(k,j) > vad(k,j-1) ) THEN vad_in_out(k,j,i) = MIN( vad_in_out(k,j,i), & vad(k,j) + overshoot_limit ) vad_in_out(k,j,i) = MAX( vad_in_out(k,j,i), & vad(k,j-1) - overshoot_limit ) ELSE vad_in_out(k,j,i) = MAX( vad_in_out(k,j,i), & vad(k,j) - overshoot_limit ) vad_in_out(k,j,i) = MIN( vad_in_out(k,j,i), & vad(k,j-1) + overshoot_limit ) ENDIF ELSE IF ( vad(k,j) > vad(k,j+1) ) THEN vad_in_out(k,j,i) = MIN( vad_in_out(k,j,i), & vad(k,j) + overshoot_limit ) vad_in_out(k,j,i) = MAX( vad_in_out(k,j,i), & vad(k,j+1) - overshoot_limit ) ELSE vad_in_out(k,j,i) = MAX( vad_in_out(k,j,i), & vad(k,j) - overshoot_limit ) vad_in_out(k,j,i) = MIN( vad_in_out(k,j,i), & vad(k,j+1) + 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, nyn DO k = nzb+1, nzt tf(k,j) = vad_in_out(k,j,i) - vad(k,j) 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(:,ny) + & 2.0 * tf(:,0) + tf(:,1) ) / wrk_long(:,0,1) DO j = 1, ny-1 DO k = nzb+1, nzt wrk_long(k,j,1) = 2.0 * ( 1.0 + long_filter_factor ) - & ( 1.0 - long_filter_factor ) * & wrk_long(k,j-1,2) wrk_long(k,j,2) = ( 1.0 - long_filter_factor ) / wrk_long(k,j,1) wrk_long(k,j,3) = ( tf(k,j-1) + 2.0 * tf(k,j) + & tf(k,j+1) - ( 1.0 - long_filter_factor ) * & wrk_long(k,j-1,3) ) / wrk_long(k,j,1) ENDDO wrk_long(:,ny,1) = 2.0 * ( 1.0 + long_filter_factor ) - & ( 1.0 - long_filter_factor ) * & wrk_long(:,ny-1,2) wrk_long(:,ny,2) = ( 1.0 - long_filter_factor ) / wrk_long(:,ny,1) wrk_long(:,ny,3) = ( tf(:,ny-1) + 2.0 * tf(:,ny) + & long_filter_factor * tf(:,0) - & ( 1.0 - long_filter_factor ) * & wrk_long(:,ny-1,3) ) / wrk_long(:,ny,1) r(:,ny) = wrk_long(:,ny,3) ENDDO DO j = ny-1, 0, -1 DO k = nzb+1, nzt r(k,j) = wrk_long(k,j,3) - wrk_long(k,j,2) * r(k,j+1) ENDDO ENDDO DO j = 0, ny DO k = nzb+1, nzt vad_in_out(k,j,i) = vad(k,j) + r(k,j) 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_y