SUBROUTINE spline_z( vad_in_out, ad_v, dz_spline, spline_tri, var_char ) !------------------------------------------------------------------------------! ! Current revisions: ! ----------------- ! ! ! Former revisions: ! ----------------- ! $Id: spline_z.f90 484 2010-02-05 07:36:54Z hoffmann $ ! ! 19 2007-02-23 04:53:48Z raasch ! Boundary condition for pt at top adjusted ! ! RCS Log replace by Id keyword, revision history cleaned up ! ! Revision 1.9 2005/06/29 08:22:56 steinfeld ! Dependency of ug and vg on height considered in the determination of the ! upper boundary condition for vad ! ! Revision 1.1 1999/02/05 09:17:16 raasch ! Initial revision ! ! ! Description: ! ------------ ! Upstream-spline advection along x ! ! Input/output parameters: ! ad_v = advecting wind speed component ! dz_spline = vertical grid spacing (dzu or dzw, depending on quantity to be ! advected) ! spline_tri = grid spacing factors (spl_tri_zu or spl_tri_zw, depending on ! quantity to be advected) ! 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 arrays_3d 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 :: dzwd, dzwu, overshoot_limit, t1, t2, t3, ups_limit REAL :: dz_spline(1:nzt+1) REAL :: spline_tri(5,nzb:nzt+1) REAL :: ad_v(nzb+1:nzta,nys:nyna,nxl:nxra) REAL, DIMENSION(:,:), ALLOCATABLE :: r, tf, vad, wrk_spline REAL, DIMENSION(:,:,:), ALLOCATABLE :: wrk_long #if defined( __parallel ) REAL :: vad_in_out(nzb+1:nzta,nys:nyna,nxl:nxra) #else REAL :: 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 ! !-- Allocate working arrays ALLOCATE( r(nzb:nzt+1,nys:nyn), vad(nzb:nzt+1,nys:nyn), & wrk_spline(nzb:nzt+1,nys:nyn) ) IF ( long_filter_factor /= 0.0 ) THEN ALLOCATE( tf(nzb:nzt+1,nys:nyn), wrk_long(nzb+1:nzt,nys:nyn,1:3) ) ENDIF ! !-- Initialize calculation of relative upstream fraction sums_up_fraction_l(component,3,:) = 0.0 ! !-- Loop over all gridpoints along x DO i = nxl, nxr ! !-- Store array to be advected on work array vad(nzb+1:nzt,:) = vad_in_out(nzb+1:nzt,nys:nyn,i) ! !-- Add boundary conditions along z IF ( var_char == 'u' .OR. var_char == 'v' ) THEN ! !-- Bottom boundary !-- u- and v-component IF ( ibc_uv_b == 0 ) THEN vad(nzb,:) = -vad(nzb+1,:) ELSE vad(nzb,:) = vad(nzb+1,:) ENDIF ! !-- Top boundary !-- Dirichlet condition IF ( ibc_uv_t == 0 .AND. var_char == 'u' ) THEN ! !-- u-component vad(nzt+1,:) = ug(nzt+1) ELSEIF ( ibc_uv_t == 0 .AND. var_char == 'v' ) THEN ! !-- v-component vad(nzt+1,:) = vg(nzt+1) ELSE ! !-- Neumann condition vad(nzt+1,:) = vad(nzt,:) ENDIF ELSEIF ( var_char == 'w' ) THEN ! !-- Bottom and top boundary for w-component vad(nzb,:) = 0.0 vad(nzt+1,:) = 0.0 ELSEIF ( var_char == 'pt' ) THEN ! !-- Bottom boundary for temperature IF ( ibc_pt_b == 1 ) THEN vad(nzb,:) = vad(nzb+1,:) ELSE vad(nzb,:) = pt(nzb,:,i) ENDIF ! !-- Top boundary for temperature IF ( ibc_pt_t == 0 ) THEN vad(nzt+1,:) = pt(nzt+1,nys:nyn,i) ELSEIF ( ibc_pt_t == 1 ) THEN vad(nzt+1,:) = vad(nzt,:) ELSEIF ( ibc_pt_t == 2 ) THEN vad(nzt+1,:) = vad(nzt,:) + bc_pt_t_val * dz_spline(nzt+1) ENDIF ELSEIF ( var_char == 'e' ) THEN ! !-- Boundary conditions for TKE (Neumann in any case) vad(nzb,:) = vad(nzb+1,:) vad(nzt,:) = vad(nzt-1,:) vad(nzt+1,:) = vad(nzt,:) ENDIF ! !-- Calculate right hand side DO j = nys, nyn r(nzb,j) = 3.0 * ( vad(nzb+1,j)-vad(nzb,j) ) / dz_spline(1) r(nzt+1,j) = 3.0 * ( vad(nzt+1,j)-vad(nzt,j) ) / dz_spline(nzt+1) DO k = nzb+1, nzt r(k,j) = 3.0 * ( & spline_tri(2,k) * ( vad(k,j)-vad(k-1,j) ) / dz_spline(k) & + spline_tri(3,k) * ( vad(k+1,j)-vad(k,j) ) / dz_spline(k+1) & ) ENDDO ENDDO ! !-- Forward substitution DO j = nys, nyn wrk_spline(nzb,j) = r(nzb,j) DO k = nzb+1, nzt+1 wrk_spline(k,j) = r(k,j) - spline_tri(5,k) * r(k-1,j) ENDDO ENDDO ! !-- Backward substitution DO j = nys, nyn r(nzt+1,j) = wrk_spline(nzt+1,j) / spline_tri(4,nzt+1) DO k = nzt, nzb, -1 r(k,j) = ( wrk_spline(k,j) - spline_tri(3,k) * r(k+1,j) ) / & spline_tri(4,k) ENDDO ENDDO ! !-- Calculate advection along z DO j = nys, nyn 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-1,j) ) <= ups_limit ) THEN vad_in_out(k,j,i) = vad(k,j) - dt_3d * ad_v(k,j,i) * & ( vad(k,j) - vad(k-1,j) ) * ddzu(k) ! !-- Calculate upstream fraction in % (s. flow_statistics) DO sr = 0, statistic_regions sums_up_fraction_l(component,3,sr) = & sums_up_fraction_l(component,3,sr) + 1.0 ENDDO ELSE t1 = ad_v(k,j,i) * dt_3d / dz_spline(k) t2 = 3.0 * ( vad(k-1,j) - vad(k,j) ) + & ( 2.0 * r(k,j) + r(k-1,j) ) * dz_spline(k) t3 = 2.0 * ( vad(k-1,j) - vad(k,j) ) + & ( r(k,j) + r(k-1,j) ) * dz_spline(k) vad_in_out(k,j,i) = vad(k,j) - r(k,j) * t1* dz_spline(k) + & t2 * t1**2 - t3 * t1**3 IF ( vad(k-1,j) == vad(k,j) ) THEN vad_in_out(k,j,i) = vad(k,j) ENDIF ENDIF ELSE IF( ABS( vad(k,j) - vad(k+1,j) ) <= ups_limit ) THEN vad_in_out(k,j,i) = vad(k,j) - dt_3d * ad_v(k,j,i) * & ( vad(k+1,j) - vad(k,j) ) * ddzu(k+1) ! !-- Calculate upstream fraction in % (s. flow_statistics) DO sr = 0, statistic_regions sums_up_fraction_l(component,3,sr) = & sums_up_fraction_l(component,3,sr) + 1.0 ENDDO ELSE t1 = -ad_v(k,j,i) * dt_3d / dz_spline(k+1) t2 = 3.0 * ( vad(k,j) - vad(k+1,j) ) + & ( 2.0 * r(k,j) + r(k+1,j) ) * dz_spline(k+1) t3 = 2.0 * ( vad(k,j) - vad(k+1,j) ) + & ( r(k,j) + r(k+1,j) ) * dz_spline(k+1) vad_in_out(k,j,i) = vad(k,j) + r(k,j)*t1*dz_spline(k+1) - & t2 * t1**2 + t3 * t1**3 IF ( vad(k+1,j) == 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 = nys, nyn DO k = nzb+1, nzt IF ( ad_v(k,j,i) > 0.0 ) THEN IF ( vad(k,j) > vad(k-1,j) ) 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-1,j) - 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-1,j) + overshoot_limit ) ENDIF ELSE IF ( vad(k,j) > vad(k+1,j) ) 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+1,j) - 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+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, nyn ! !-- Depending on the quantity to be advected, the respective vertical !-- boundary conditions must be applied. IF ( var_char == 'u' .OR. var_char == 'v' ) THEN IF ( ibc_uv_b == 0 ) THEN tf(nzb,j) = - ( vad_in_out(nzb+1,j,i) - vad(nzb+1,j) ) ELSE tf(nzb,j) = vad_in_out(nzb+1,j,i) - vad(nzb+1,j) ENDIF IF ( ibc_uv_t == 0 ) THEN tf(nzt+1,j) = 0.0 ELSE tf(nzt+1,j) = vad_in_out(nzt,j,i) - vad(nzt,j) ENDIF ELSEIF ( var_char == 'w' ) THEN tf(nzb,j) = 0.0 tf(nzt+1,j) = 0.0 ELSEIF ( var_char == 'pt' ) THEN IF ( ibc_pt_b == 1 ) THEN tf(nzb,j) = vad_in_out(nzb+1,j,i) - vad(nzb+1,j) ELSE tf(nzb,j) = 0.0 ENDIF IF ( ibc_pt_t == 1 ) THEN vad_in_out(nzt,j,i) = vad_in_out(nzt-1,j,i) + bc_pt_t_val * & dz_spline(nzt) tf(nzt+1,j) = vad_in_out(nzt,j,i) + bc_pt_t_val * & dz_spline(nzt+1) - vad(nzt+1,j) ELSE vad_in_out(nzt,j,i) = pt(nzt,j,i) tf(nzt+1,j) = 0.0 ENDIF ENDIF DO k = nzb+1, nzt tf(k,j) = vad_in_out(k,j,i) - vad(k,j) ENDDO ENDDO ! !-- Apply the filter. DO j = nys, nyn dzwd = dz_spline(1) / ( dz_spline(1) + dz_spline(2) ) dzwu = dz_spline(2) / ( dz_spline(1) + dz_spline(2) ) wrk_long(nzb+1,j,1) = 2.0 * ( 1.0 + long_filter_factor ) wrk_long(nzb+1,j,2) = ( 1.0 - long_filter_factor ) * dzwd / & wrk_long(nzb+1,j,1) wrk_long(nzb+1,j,3) = ( long_filter_factor * dzwu * tf(nzb,j) + & 2.0 * tf(nzb+1,j) + dzwd * tf(nzb+2,j) & ) / wrk_long(nzb+1,j,1) DO k = nzb+2, nzt-1 dzwd = dz_spline(k) / ( dz_spline(k) + dz_spline(k+1) ) dzwu = dz_spline(k+1) / ( dz_spline(k) + dz_spline(k+1) ) wrk_long(k,j,1) = 2.0 * ( 1.0 + long_filter_factor ) - & ( 1.0 - long_filter_factor ) * dzwu * & wrk_long(k-1,j,2) wrk_long(k,j,2) = ( 1.0 - long_filter_factor ) * dzwd / & wrk_long(k,j,1) wrk_long(k,j,3) = ( dzwu * tf(k-1,j) + 2.0 * tf(k,j) + & dzwd * tf(k+1,j) - & ( 1.0 - long_filter_factor ) * dzwu * & wrk_long(k-1,j,3) & ) / wrk_long(k,j,1) ENDDO dzwd = dz_spline(nzt) / ( dz_spline(nzt) + dz_spline(nzt+1) ) dzwu = dz_spline(nzt+1) / ( dz_spline(nzt) + dz_spline(nzt+1) ) wrk_long(nzt,j,1) = 2.0 * ( 1.0 + long_filter_factor ) - & ( 1.0 - long_filter_factor ) * dzwu * & wrk_long(nzt-1,j,2) wrk_long(nzt,j,2) = ( 1.0 - long_filter_factor ) * dzwd / & wrk_long(nzt,j,1) wrk_long(nzt,j,3) = ( dzwu * tf(nzt-1,j) + 2.0 * tf(nzt,j) + & dzwd * long_filter_factor * tf(nzt+1,j) - & ( 1.0 - long_filter_factor ) * dzwu * & wrk_long(nzt-1,j,3) & ) / wrk_long(nzt,j,1) r(nzt,j) = wrk_long(nzt,j,3) ENDDO DO j = nys, nyn DO k = nzt-1, nzb+1, -1 r(k,j) = wrk_long(k,j,3) - wrk_long(k,j,2) * r(k+1,j) ENDDO ENDDO DO j = nys, nyn DO k = nzb+1, nzt vad_in_out(k,j,i) = vad(k,j) + r(k,j) ENDDO ENDDO ENDIF ! Long filter ENDDO DEALLOCATE( r, vad, wrk_spline ) IF ( long_filter_factor /= 0.0 ) DEALLOCATE( tf, wrk_long ) END SUBROUTINE spline_z