Changeset 4079 for palm/trunk/SOURCE/advec_ws.f90
- Timestamp:
- Jul 9, 2019 6:04:41 PM (5 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/advec_ws.f90 (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/advec_ws.f90
r3873 r4079 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Implementation of a flux limiter according to Skamarock (2006) for the 28 ! vertical scalar advection. Please note, this is only implemented for the 29 ! cache-optimized version at the moment. Implementation for the vector- 30 ! optimized version will follow after critical issues concerning 31 ! vectorization are fixed. 32 ! 33 ! 3873 2019-04-08 15:44:30Z knoop 27 34 ! Moved ocean_mode specific code to ocean_mod 28 35 ! … … 273 280 !> the divergence is not sufficiently reduced, resulting in erroneous fluxes 274 281 !> and could lead to numerical instabilities. 282 !> 283 !> @todo Implement monotonic flux limiter also for vector version. 275 284 !------------------------------------------------------------------------------! 276 285 MODULE advec_ws … … 296 305 bc_dirichlet_s, bc_radiation_l, bc_radiation_n, & 297 306 bc_radiation_r, bc_radiation_s, intermediate_timestep_count, & 298 u_gtrans, v_gtrans 307 u_gtrans, v_gtrans, dt_3d 299 308 300 309 USE indices, & … … 1116 1125 SUBROUTINE advec_s_ws_ij( i, j, sk, sk_char, swap_flux_y_local, & 1117 1126 swap_diss_y_local, swap_flux_x_local, & 1118 swap_diss_x_local, i_omp, tn )1127 swap_diss_x_local, i_omp, tn, flux_limitation ) 1119 1128 1120 1129 … … 1126 1135 INTEGER(iwp) :: k !< grid index along z-direction 1127 1136 INTEGER(iwp) :: k_mm !< k-2 index in disretization, can be modified to avoid segmentation faults 1137 INTEGER(iwp) :: k_mmm !< k-3 index in disretization, can be modified to avoid segmentation faults 1128 1138 INTEGER(iwp) :: k_pp !< k+2 index in disretization, can be modified to avoid segmentation faults 1129 1139 INTEGER(iwp) :: k_ppp !< k+3 index in disretization, can be modified to avoid segmentation faults 1130 1140 INTEGER(iwp) :: nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms 1131 1141 INTEGER(iwp) :: tn !< number of OpenMP thread 1132 1133 REAL(wp) :: ibit0 !< flag indicating 1st-order scheme along x-direction 1134 REAL(wp) :: ibit1 !< flag indicating 3rd-order scheme along x-direction 1135 REAL(wp) :: ibit2 !< flag indicating 5th-order scheme along x-direction 1136 REAL(wp) :: ibit3 !< flag indicating 1st-order scheme along y-direction 1137 REAL(wp) :: ibit4 !< flag indicating 3rd-order scheme along y-direction 1138 REAL(wp) :: ibit5 !< flag indicating 5th-order scheme along y-direction 1139 REAL(wp) :: ibit6 !< flag indicating 1st-order scheme along z-direction 1140 REAL(wp) :: ibit7 !< flag indicating 3rd-order scheme along z-direction 1141 REAL(wp) :: ibit8 !< flag indicating 5th-order scheme along z-direction 1142 REAL(wp) :: diss_d !< artificial dissipation term at grid box bottom 1143 REAL(wp) :: div !< diverence on scalar grid 1144 REAL(wp) :: flux_d !< 6th-order flux at grid box bottom 1145 REAL(wp) :: u_comp !< advection velocity along x-direction 1146 REAL(wp) :: v_comp !< advection velocity along y-direction 1142 1143 LOGICAL, OPTIONAL :: flux_limitation !< flag indicating flux limitation of the vertical advection 1144 LOGICAL :: limiter !< control flag indicating the application of flux limitation 1145 1146 REAL(wp) :: diss_d !< artificial dissipation term at grid box bottom 1147 REAL(wp) :: div !< velocity diverence on scalar grid 1148 REAL(wp) :: div_in !< vertical flux divergence of ingoing fluxes 1149 REAL(wp) :: div_out !< vertical flux divergence of outgoing fluxes 1150 REAL(wp) :: f_corr_t !< correction flux at grid-cell top, i.e. the difference between high and low-order flux 1151 REAL(wp) :: f_corr_d !< correction flux at grid-cell bottom, i.e. the difference between high and low-order flux 1152 REAL(wp) :: f_corr_t_in !< correction flux of ingoing flux part at grid-cell top 1153 REAL(wp) :: f_corr_d_in !< correction flux of ingoing flux part at grid-cell bottom 1154 REAL(wp) :: f_corr_t_out !< correction flux of outgoing flux part at grid-cell top 1155 REAL(wp) :: f_corr_d_out !< correction flux of outgoing flux part at grid-cell bottom 1156 REAL(wp) :: fac_correction!< factor to limit the in- and outgoing fluxes 1157 REAL(wp) :: flux_d !< 6th-order flux at grid box bottom 1158 REAL(wp) :: ibit0 !< flag indicating 1st-order scheme along x-direction 1159 REAL(wp) :: ibit1 !< flag indicating 3rd-order scheme along x-direction 1160 REAL(wp) :: ibit2 !< flag indicating 5th-order scheme along x-direction 1161 REAL(wp) :: ibit3 !< flag indicating 1st-order scheme along y-direction 1162 REAL(wp) :: ibit4 !< flag indicating 3rd-order scheme along y-direction 1163 REAL(wp) :: ibit5 !< flag indicating 5th-order scheme along y-direction 1164 REAL(wp) :: ibit6 !< flag indicating 1st-order scheme along z-direction 1165 REAL(wp) :: ibit7 !< flag indicating 3rd-order scheme along z-direction 1166 REAL(wp) :: ibit8 !< flag indicating 5th-order scheme along z-direction 1167 REAL(wp) :: max_val !< maximum value of the quanitity along the numerical stencil (in vertical direction) 1168 REAL(wp) :: min_val !< maximum value of the quanitity along the numerical stencil (in vertical direction) 1169 REAL(wp) :: mon !< monotone solution of the advection equation using 1st-order fluxes 1170 REAL(wp) :: u_comp !< advection velocity along x-direction 1171 REAL(wp) :: v_comp !< advection velocity along y-direction 1147 1172 ! 1148 1173 !-- sk is an array from parameter list. It should not be a pointer, because … … 1153 1178 REAL(wp), INTENT(IN),DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk !< advected scalar 1154 1179 1155 REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !< discretized artificial dissipation at northward-side of the grid box 1156 REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !< discretized artificial dissipation at rightward-side of the grid box 1157 REAL(wp), DIMENSION(nzb:nzt+1) :: diss_t !< discretized artificial dissipation at top of the grid box 1158 REAL(wp), DIMENSION(nzb:nzt+1) :: flux_n !< discretized 6th-order flux at northward-side of the grid box 1159 REAL(wp), DIMENSION(nzb:nzt+1) :: flux_r !< discretized 6th-order flux at rightward-side of the grid box 1160 REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !< discretized 6th-order flux at top of the grid box 1180 REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !< discretized artificial dissipation at northward-side of the grid box 1181 REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !< discretized artificial dissipation at rightward-side of the grid box 1182 REAL(wp), DIMENSION(nzb:nzt+1) :: diss_t !< discretized artificial dissipation at top of the grid box 1183 REAL(wp), DIMENSION(nzb:nzt+1) :: flux_n !< discretized 6th-order flux at northward-side of the grid box 1184 REAL(wp), DIMENSION(nzb:nzt+1) :: flux_r !< discretized 6th-order flux at rightward-side of the grid box 1185 REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !< discretized 6th-order flux at top of the grid box 1186 REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t_1st !< discretized 1st-order flux at top of the grid box 1161 1187 1162 1188 REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: swap_diss_y_local !< discretized artificial dissipation at southward-side of the grid box … … 1178 1204 nzb_max_l = nzb_max 1179 1205 END IF 1206 ! 1207 !-- Set control flag for flux limiter 1208 limiter = .FALSE. 1209 IF ( PRESENT( flux_limitation) ) limiter = flux_limitation 1180 1210 ! 1181 1211 !-- Compute southside fluxes of the respective PE bounds. … … 1231 1261 + ( sk(k,j+2,i) + sk(k,j-3,i) ) & 1232 1262 ) * adv_sca_5 1233 swap_diss_y_local(k,tn) = -ABS( v_comp ) * (&1263 swap_diss_y_local(k,tn) = -ABS( v_comp ) * ( & 1234 1264 10.0_wp * ( sk(k,j,i) - sk(k,j-1,i) ) & 1235 1265 - 5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) ) & 1236 1266 + sk(k,j+2,i) - sk(k,j-3,i) & 1237 1267 ) * adv_sca_5 1238 1268 1239 1269 ENDDO … … 1266 1296 ) 1267 1297 1268 swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * (&1298 swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * ( & 1269 1299 ( 10.0_wp * ibit2 * adv_sca_5 & 1270 1300 + 3.0_wp * ibit1 * adv_sca_3 & … … 1279 1309 ) * & 1280 1310 ( sk(k,j,i+2) - sk(k,j,i-3) ) & 1281 1311 ) 1282 1312 1283 1313 ENDDO … … 1535 1565 ( sk(k_ppp,j,i) - sk(k_mm,j,i) ) & 1536 1566 ) 1537 ENDDO 1567 ENDDO 1568 1569 IF ( limiter ) THEN 1570 ! 1571 !-- Compute monotone first-order fluxes which are required for mononte 1572 !-- flux limitation. 1573 flux_t_1st(nzb) = 0.0_wp 1574 DO k = nzb+1, nzb_max_l 1575 flux_t_1st(k) = ( w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) & 1576 - ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) ) ) & 1577 * rho_air_zw(k) * adv_sca_1 1578 ! 1579 !-- In flux limitation the total flux will be corrected. For the sake 1580 !-- of cleariness the higher-order advective and disspative fluxes 1581 !-- will be merged onto flux_t. 1582 flux_t(k) = flux_t(k) + diss_t(k) 1583 diss_t(k) = 0.0_wp 1584 ENDDO 1585 ! 1586 !-- Flux limitation of vertical fluxes according to Skamarock (2006). 1587 !-- Please note, as flux limitation implies linear dependencies of fluxes, 1588 !-- flux limitation is only made for the vertical advection term. Limitation 1589 !-- of the horizontal terms cannot be parallelized. 1590 !-- Due to the linear dependency, the following loop will not be vectorized. 1591 !-- Further, note that the flux limiter is only applied within the urban 1592 !-- layer, i.e up to the topography top. 1593 DO k = nzb+1, nzb_max_l 1594 ! 1595 !-- Compute one-dimensional divergence along the vertical direction, 1596 !-- which is used to correct the advection discretization. This is 1597 !-- necessary as in one-dimensional space the advection velocity 1598 !-- should actually be constant. 1599 div = ( w(k,j,i) * rho_air_zw(k) & 1600 - w(k-1,j,i) * rho_air_zw(k-1) & 1601 ) * drho_air(k) * ddzw(k) 1602 ! 1603 !-- Compute monotone solution of the advection equation from 1604 !-- 1st-order fluxes. Please note, the advection equation is corrected 1605 !-- by the divergence term (in 1D the advective flow should be divergence 1606 !-- free). Moreover, please note, as time-increment the full timestep 1607 !-- is used, even though a Runge-Kutta scheme will be used. However, 1608 !-- the length of the actual time increment is not important at all 1609 !-- since it cancels out later when the fluxes are limited. 1610 mon = sk(k,j,i) + ( - ( flux_t_1st(k) - flux_t_1st(k-1) ) & 1611 * drho_air(k) * ddzw(k) & 1612 + div * sk(k,j,i) & 1613 ) * dt_3d 1614 ! 1615 !-- Determine minimum and maximum values along the numerical stencil. 1616 k_mmm = MAX( k - 3, nzb + 1 ) 1617 k_ppp = MIN( k + 3, nzt + 1 ) 1618 1619 min_val = MINVAL( sk(k_mmm:k_ppp,j,i) ) 1620 max_val = MAXVAL( sk(k_mmm:k_ppp,j,i) ) 1621 ! 1622 !-- Compute difference between high- and low-order fluxes, which may 1623 !-- act as correction fluxes 1624 f_corr_t = flux_t(k) - flux_t_1st(k) 1625 f_corr_d = flux_t(k-1) - flux_t_1st(k-1) 1626 ! 1627 !-- Determine outgoing fluxes, i.e. the part of the fluxes which can 1628 !-- decrease the value within the grid box 1629 f_corr_t_out = MAX( 0.0_wp, f_corr_t ) 1630 f_corr_d_out = MIN( 0.0_wp, f_corr_d ) 1631 ! 1632 !-- Determine ingoing fluxes, i.e. the part of the fluxes which can 1633 !-- increase the value within the grid box 1634 f_corr_t_in = MIN( 0.0_wp, f_corr_t) 1635 f_corr_d_in = MAX( 0.0_wp, f_corr_d) 1636 ! 1637 !-- Compute divergence of outgoing correction fluxes 1638 div_out = - ( f_corr_t_out - f_corr_d_out ) * drho_air(k) & 1639 * ddzw(k) * dt_3d 1640 ! 1641 !-- Compute divergence of ingoing correction fluxes 1642 div_in = - ( f_corr_t_in - f_corr_d_in ) * drho_air(k) & 1643 * ddzw(k) * dt_3d 1644 ! 1645 !-- Check if outgoing fluxes can lead to undershoots, i.e. values smaller 1646 !-- than the minimum value within the numerical stencil. If so, limit 1647 !-- them. 1648 IF ( mon - min_val < - div_out .AND. ABS( div_out ) > 0.0_wp ) & 1649 THEN 1650 fac_correction = ( mon - min_val ) / ( - div_out ) 1651 f_corr_t_out = f_corr_t_out * fac_correction 1652 f_corr_d_out = f_corr_d_out * fac_correction 1653 ENDIF 1654 ! 1655 !-- Check if ingoing fluxes can lead to overshoots, i.e. values larger 1656 !-- than the maximum value within the numerical stencil. If so, limit 1657 !-- them. 1658 IF ( mon - max_val > - div_in .AND. ABS( div_in ) > 0.0_wp ) & 1659 THEN 1660 fac_correction = ( mon - max_val ) / ( - div_in ) 1661 f_corr_t_in = f_corr_t_in * fac_correction 1662 f_corr_d_in = f_corr_d_in * fac_correction 1663 ENDIF 1664 ! 1665 !-- Finally add the limited fluxes to the original ones. If no 1666 !-- flux limitation was done, the fluxes equal the original ones. 1667 flux_t(k) = flux_t_1st(k) + f_corr_t_out + f_corr_t_in 1668 flux_t(k-1) = flux_t_1st(k-1) + f_corr_d_out + f_corr_d_in 1669 ENDDO 1670 ENDIF 1538 1671 1539 1672 DO k = nzb+1, nzb_max_l
Note: See TracChangeset
for help on using the changeset viewer.
