Changeset 4079 for palm/trunk/SOURCE/advec_ws.f90
 Timestamp:
 Jul 9, 2019 6:04:41 PM (4 years ago)
 File:

 1 edited
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 ! cacheoptimized version at the moment. Implementation for the vector 30 ! optimized version will follow after critical issues concerning 31 ! vectorization are fixed. 32 ! 33 ! 3873 20190408 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 zdirection 1127 1136 INTEGER(iwp) :: k_mm !< k2 index in disretization, can be modified to avoid segmentation faults 1137 INTEGER(iwp) :: k_mmm !< k3 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 1storder scheme along xdirection 1134 REAL(wp) :: ibit1 !< flag indicating 3rdorder scheme along xdirection 1135 REAL(wp) :: ibit2 !< flag indicating 5thorder scheme along xdirection 1136 REAL(wp) :: ibit3 !< flag indicating 1storder scheme along ydirection 1137 REAL(wp) :: ibit4 !< flag indicating 3rdorder scheme along ydirection 1138 REAL(wp) :: ibit5 !< flag indicating 5thorder scheme along ydirection 1139 REAL(wp) :: ibit6 !< flag indicating 1storder scheme along zdirection 1140 REAL(wp) :: ibit7 !< flag indicating 3rdorder scheme along zdirection 1141 REAL(wp) :: ibit8 !< flag indicating 5thorder scheme along zdirection 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 !< 6thorder flux at grid box bottom 1145 REAL(wp) :: u_comp !< advection velocity along xdirection 1146 REAL(wp) :: v_comp !< advection velocity along ydirection 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 gridcell top, i.e. the difference between high and loworder flux 1151 REAL(wp) :: f_corr_d !< correction flux at gridcell bottom, i.e. the difference between high and loworder flux 1152 REAL(wp) :: f_corr_t_in !< correction flux of ingoing flux part at gridcell top 1153 REAL(wp) :: f_corr_d_in !< correction flux of ingoing flux part at gridcell bottom 1154 REAL(wp) :: f_corr_t_out !< correction flux of outgoing flux part at gridcell top 1155 REAL(wp) :: f_corr_d_out !< correction flux of outgoing flux part at gridcell bottom 1156 REAL(wp) :: fac_correction!< factor to limit the in and outgoing fluxes 1157 REAL(wp) :: flux_d !< 6thorder flux at grid box bottom 1158 REAL(wp) :: ibit0 !< flag indicating 1storder scheme along xdirection 1159 REAL(wp) :: ibit1 !< flag indicating 3rdorder scheme along xdirection 1160 REAL(wp) :: ibit2 !< flag indicating 5thorder scheme along xdirection 1161 REAL(wp) :: ibit3 !< flag indicating 1storder scheme along ydirection 1162 REAL(wp) :: ibit4 !< flag indicating 3rdorder scheme along ydirection 1163 REAL(wp) :: ibit5 !< flag indicating 5thorder scheme along ydirection 1164 REAL(wp) :: ibit6 !< flag indicating 1storder scheme along zdirection 1165 REAL(wp) :: ibit7 !< flag indicating 3rdorder scheme along zdirection 1166 REAL(wp) :: ibit8 !< flag indicating 5thorder scheme along zdirection 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 1storder fluxes 1170 REAL(wp) :: u_comp !< advection velocity along xdirection 1171 REAL(wp) :: v_comp !< advection velocity along ydirection 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 northwardside of the grid box 1156 REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !< discretized artificial dissipation at rightwardside 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 6thorder flux at northwardside of the grid box 1159 REAL(wp), DIMENSION(nzb:nzt+1) :: flux_r !< discretized 6thorder flux at rightwardside of the grid box 1160 REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !< discretized 6thorder flux at top of the grid box 1180 REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !< discretized artificial dissipation at northwardside of the grid box 1181 REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !< discretized artificial dissipation at rightwardside 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 6thorder flux at northwardside of the grid box 1184 REAL(wp), DIMENSION(nzb:nzt+1) :: flux_r !< discretized 6thorder flux at rightwardside of the grid box 1185 REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !< discretized 6thorder flux at top of the grid box 1186 REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t_1st !< discretized 1storder flux at top of the grid box 1161 1187 1162 1188 REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task1) :: swap_diss_y_local !< discretized artificial dissipation at southwardside 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,j3,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,j1,i) ) & 1235 1265  5.0_wp * ( sk(k,j+1,i)  sk(k,j2,i) ) & 1236 1266 + sk(k,j+2,i)  sk(k,j3,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,i3) ) & 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 firstorder 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 higherorder 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 onedimensional divergence along the vertical direction, 1596 ! which is used to correct the advection discretization. This is 1597 ! necessary as in onedimensional space the advection velocity 1598 ! should actually be constant. 1599 div = ( w(k,j,i) * rho_air_zw(k) & 1600  w(k1,j,i) * rho_air_zw(k1) & 1601 ) * drho_air(k) * ddzw(k) 1602 ! 1603 ! Compute monotone solution of the advection equation from 1604 ! 1storder 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 timeincrement the full timestep 1607 ! is used, even though a RungeKutta 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(k1) ) & 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 loworder 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(k1)  flux_t_1st(k1) 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(k1) = flux_t_1st(k1) + 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.