Changeset 3661 for palm/trunk/SOURCE/advec_ws.f90
- Timestamp:
- Jan 8, 2019 6:22:50 PM (6 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/advec_ws.f90 (modified) (91 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/advec_ws.f90
r3655 r3661 25 25 ! ----------------- 26 26 ! $Id$ 27 ! - Minor bugfix in divergence correction (only has implications at 28 ! downward-facing wall surfaces) 29 ! - Remove setting of Neumann condition for horizontal velocity variances 30 ! - Split loops for tendency calculation and divergence correction in order to 31 ! reduce bit queries 32 ! - Introduce new parameter nzb_max_l to better control order degradation at 33 ! non-cyclic boundaries 34 ! 35 ! 3655 2019-01-07 16:51:22Z knoop 27 36 ! OpenACC port for SPEC 28 37 ! … … 247 256 !> degraded. 248 257 !> A divergence correction is applied. It is necessary for topography, since 249 !> the divergence is not sufficiently reduced, resulting in erroneous fluxes and250 !> partlynumerical instabilities.258 !> the divergence is not sufficiently reduced, resulting in erroneous fluxes 259 !> and could lead to numerical instabilities. 251 260 !-----------------------------------------------------------------------------! 252 261 MODULE advec_ws … … 1192 1201 1193 1202 USE control_parameters, & 1194 ONLY: intermediate_timestep_count, u_gtrans, v_gtrans 1203 ONLY: bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, & 1204 bc_dirichlet_s, bc_radiation_l, bc_radiation_n, & 1205 bc_radiation_r, bc_radiation_s, intermediate_timestep_count, & 1206 u_gtrans, v_gtrans 1195 1207 1196 1208 USE grid_variables, & … … 1218 1230 CHARACTER (LEN = *), INTENT(IN) :: sk_char !< string identifier, used for assign fluxes to the correct dimension in the analysis array 1219 1231 1220 INTEGER(iwp) :: i !< grid index along x-direction 1221 INTEGER(iwp) :: i_omp !< leftmost index on subdomain, or in case of OpenMP, on thread 1222 INTEGER(iwp) :: j !< grid index along y-direction 1223 INTEGER(iwp) :: k !< grid index along z-direction 1224 INTEGER(iwp) :: k_mm !< k-2 index in disretization, can be modified to avoid segmentation faults 1225 INTEGER(iwp) :: k_pp !< k+2 index in disretization, can be modified to avoid segmentation faults 1226 INTEGER(iwp) :: k_ppp !< k+3 index in disretization, can be modified to avoid segmentation faults 1227 INTEGER(iwp) :: tn !< number of OpenMP thread 1232 INTEGER(iwp) :: i !< grid index along x-direction 1233 INTEGER(iwp) :: i_omp !< leftmost index on subdomain, or in case of OpenMP, on thread 1234 INTEGER(iwp) :: j !< grid index along y-direction 1235 INTEGER(iwp) :: k !< grid index along z-direction 1236 INTEGER(iwp) :: k_mm !< k-2 index in disretization, can be modified to avoid segmentation faults 1237 INTEGER(iwp) :: k_pp !< k+2 index in disretization, can be modified to avoid segmentation faults 1238 INTEGER(iwp) :: k_ppp !< k+3 index in disretization, can be modified to avoid segmentation faults 1239 INTEGER(iwp) :: nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms 1240 INTEGER(iwp) :: tn !< number of OpenMP thread 1228 1241 1229 1242 REAL(wp) :: ibit0 !< flag indicating 1st-order scheme along x-direction … … 1261 1274 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: swap_diss_x_local !< discretized artificial dissipation at leftward-side of the grid box 1262 1275 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: swap_flux_x_local !< discretized 6th-order flux at leftward-side of the grid box 1263 1264 1276 ! 1277 !-- Used local modified copy of nzb_max (used to degrade order of 1278 !-- discretization) at non-cyclic boundaries. Modify only at relevant points 1279 !-- instead of the entire subdomain. This should lead to better 1280 !-- load balance between boundary and non-boundary PEs. 1281 IF( ( bc_dirichlet_l .OR. bc_radiation_l ) .AND. i <= nxl + 2 .OR. & 1282 ( bc_dirichlet_r .OR. bc_radiation_r ) .AND. i >= nxr - 2 .OR. & 1283 ( bc_dirichlet_s .OR. bc_radiation_s ) .AND. j <= nys + 2 .OR. & 1284 ( bc_dirichlet_n .OR. bc_radiation_n ) .AND. j >= nyn - 2 ) THEN 1285 nzb_max_l = nzt 1286 ELSE 1287 nzb_max_l = nzb_max 1288 END IF 1265 1289 ! 1266 1290 !-- Compute southside fluxes of the respective PE bounds. … … 1268 1292 ! 1269 1293 !-- Up to the top of the highest topography. 1270 DO k = nzb+1, nzb_max 1294 DO k = nzb+1, nzb_max_l 1271 1295 1272 1296 ibit5 = REAL( IBITS(advc_flags_1(k,j-1,i),5,1), KIND = wp ) … … 1329 1353 IF ( i == i_omp ) THEN 1330 1354 1331 DO k = nzb+1, nzb_max 1355 DO k = nzb+1, nzb_max_l 1332 1356 1333 1357 ibit2 = REAL( IBITS(advc_flags_1(k,j,i-1),2,1), KIND = wp ) … … 1368 1392 ENDDO 1369 1393 1370 DO k = nzb_max +1, nzt1394 DO k = nzb_max_l+1, nzt 1371 1395 1372 1396 u_comp = u(k,j,i) - u_gtrans + u_stokes_zu(k) … … 1389 1413 !-- Now compute the fluxes and tendency terms for the horizontal and 1390 1414 !-- vertical parts up to the top of the highest topography. 1391 DO k = nzb+1, nzb_max 1415 DO k = nzb+1, nzb_max_l 1392 1416 ! 1393 1417 !-- Note: It is faster to conduct all multiplications explicitly, e.g. … … 1469 1493 !-- vertical parts above the top of the highest topography. No degradation 1470 1494 !-- for the horizontal parts, but for the vertical it is stell needed. 1471 DO k = nzb_max +1, nzt1495 DO k = nzb_max_l+1, nzt 1472 1496 1473 1497 u_comp = u(k,j,i+1) - u_gtrans + u_stokes_zu(k) … … 1622 1646 ENDDO 1623 1647 1624 DO k = nzb+1, nz t1648 DO k = nzb+1, nzb_max_l 1625 1649 1626 1650 flux_d = flux_t(k-1) 1627 1651 diss_d = diss_t(k-1) 1628 1652 1653 ibit2 = REAL( IBITS(advc_flags_1(k,j,i),2,1), KIND = wp ) 1654 ibit1 = REAL( IBITS(advc_flags_1(k,j,i),1,1), KIND = wp ) 1655 ibit0 = REAL( IBITS(advc_flags_1(k,j,i),0,1), KIND = wp ) 1656 1657 ibit5 = REAL( IBITS(advc_flags_1(k,j,i),5,1), KIND = wp ) 1658 ibit4 = REAL( IBITS(advc_flags_1(k,j,i),4,1), KIND = wp ) 1659 ibit3 = REAL( IBITS(advc_flags_1(k,j,i),3,1), KIND = wp ) 1660 1661 ibit8 = REAL( IBITS(advc_flags_1(k,j,i),8,1), KIND = wp ) 1662 ibit7 = REAL( IBITS(advc_flags_1(k,j,i),7,1), KIND = wp ) 1663 ibit6 = REAL( IBITS(advc_flags_1(k,j,i),6,1), KIND = wp ) 1629 1664 ! 1630 1665 !-- Calculate the divergence of the velocity field. A respective … … 1654 1689 ) & 1655 1690 ) * drho_air(k) * ddzw(k) 1691 1692 tend(k,j,i) = tend(k,j,i) - ( & 1693 ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - & 1694 swap_diss_x_local(k,j,tn) ) * ddx & 1695 + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn) - & 1696 swap_diss_y_local(k,tn) ) * ddy & 1697 + ( ( flux_t(k) + diss_t(k) ) - & 1698 ( flux_d + diss_d ) & 1699 ) * drho_air(k) * ddzw(k) & 1700 ) + sk(k,j,i) * div 1701 1702 1703 swap_flux_y_local(k,tn) = flux_n(k) 1704 swap_diss_y_local(k,tn) = diss_n(k) 1705 swap_flux_x_local(k,j,tn) = flux_r(k) 1706 swap_diss_x_local(k,j,tn) = diss_r(k) 1707 1708 ENDDO 1709 1710 DO k = nzb_max_l+1, nzt 1711 1712 flux_d = flux_t(k-1) 1713 diss_d = diss_t(k-1) 1714 ! 1715 !-- Calculate the divergence of the velocity field. A respective 1716 !-- correction is needed to overcome numerical instabilities introduced 1717 !-- by a not sufficient reduction of divergences near topography. 1718 div = ( u(k,j,i+1) - u(k,j,i) ) * ddx & 1719 + ( v(k,j+1,i) - v(k,j,i) ) * ddy & 1720 + ( w(k,j,i) * rho_air_zw(k) & 1721 - w(k-1,j,i) * rho_air_zw(k-1) & 1722 ) * drho_air(k) * ddzw(k) 1656 1723 1657 1724 tend(k,j,i) = tend(k,j,i) - ( & … … 1795 1862 SUBROUTINE advec_u_ws_ij( i, j, i_omp, tn ) 1796 1863 1797 USE arrays_3d, &1798 ONLY: ddzw, diss_l_u, diss_s_u, flux_l_u, flux_s_u, tend, u, v, w, &1864 USE arrays_3d, & 1865 ONLY: ddzw, diss_l_u, diss_s_u, flux_l_u, flux_s_u, tend, u, v, w, & 1799 1866 drho_air, rho_air_zw 1800 1867 1801 USE control_parameters, & 1802 ONLY: intermediate_timestep_count, u_gtrans, v_gtrans 1803 1804 USE grid_variables, & 1868 USE control_parameters, & 1869 ONLY: bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, & 1870 bc_dirichlet_s, bc_radiation_l, bc_radiation_n, & 1871 bc_radiation_r, bc_radiation_s, intermediate_timestep_count, & 1872 u_gtrans, v_gtrans 1873 1874 USE grid_variables, & 1805 1875 ONLY: ddx, ddy 1806 1876 1807 USE indices, &1808 ONLY: n xlu, nys, nzb, nzb_max, nzt, advc_flags_11877 USE indices, & 1878 ONLY: nyn, nys, nxl, nxlu, nxr, nzb, nzb_max, nzt, advc_flags_1 1809 1879 1810 1880 USE kinds 1811 1881 1812 USE statistics, &1882 USE statistics, & 1813 1883 ONLY: hom, sums_us2_ws_l, sums_wsus_ws_l, weight_substep 1814 1884 1815 1885 IMPLICIT NONE 1816 1886 1817 INTEGER(iwp) :: i !< grid index along x-direction 1818 INTEGER(iwp) :: i_omp !< leftmost index on subdomain, or in case of OpenMP, on thread 1819 INTEGER(iwp) :: j !< grid index along y-direction 1820 INTEGER(iwp) :: k !< grid index along z-direction 1821 INTEGER(iwp) :: k_mm !< k-2 index in disretization, can be modified to avoid segmentation faults 1822 INTEGER(iwp) :: k_pp !< k+2 index in disretization, can be modified to avoid segmentation faults 1823 INTEGER(iwp) :: k_ppp !< k+3 index in disretization, can be modified to avoid segmentation faults 1824 INTEGER(iwp) :: tn !< number of OpenMP thread 1887 INTEGER(iwp) :: i !< grid index along x-direction 1888 INTEGER(iwp) :: i_omp !< leftmost index on subdomain, or in case of OpenMP, on thread 1889 INTEGER(iwp) :: j !< grid index along y-direction 1890 INTEGER(iwp) :: k !< grid index along z-direction 1891 INTEGER(iwp) :: k_mm !< k-2 index in disretization, can be modified to avoid segmentation faults 1892 INTEGER(iwp) :: k_pp !< k+2 index in disretization, can be modified to avoid segmentation faults 1893 INTEGER(iwp) :: k_ppp !< k+3 index in disretization, can be modified to avoid segmentation faults 1894 INTEGER(iwp) :: nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms 1895 INTEGER(iwp) :: tn !< number of OpenMP thread 1825 1896 1826 1897 REAL(wp) :: ibit9 !< flag indicating 1st-order scheme along x-direction … … 1849 1920 REAL(wp), DIMENSION(nzb:nzt+1) :: v_comp !< advection velocity along y 1850 1921 REAL(wp), DIMENSION(nzb:nzt+1) :: w_comp !< advection velocity along z 1851 1922 ! 1923 !-- Used local modified copy of nzb_max (used to degrade order of 1924 !-- discretization) at non-cyclic boundaries. Modify only at relevant points 1925 !-- instead of the entire subdomain. This should lead to better 1926 !-- load balance between boundary and non-boundary PEs. 1927 IF( ( bc_dirichlet_l .OR. bc_radiation_l ) .AND. i <= nxl + 2 .OR. & 1928 ( bc_dirichlet_r .OR. bc_radiation_r ) .AND. i >= nxr - 2 .OR. & 1929 ( bc_dirichlet_s .OR. bc_radiation_s ) .AND. j <= nys + 2 .OR. & 1930 ( bc_dirichlet_n .OR. bc_radiation_n ) .AND. j >= nyn - 2 ) THEN 1931 nzb_max_l = nzt 1932 ELSE 1933 nzb_max_l = nzb_max 1934 END IF 1935 1852 1936 gu = 2.0_wp * u_gtrans 1853 1937 gv = 2.0_wp * v_gtrans … … 1856 1940 IF ( j == nys ) THEN 1857 1941 1858 DO k = nzb+1, nzb_max 1942 DO k = nzb+1, nzb_max_l 1859 1943 1860 1944 ibit14 = REAL( IBITS(advc_flags_1(k,j-1,i),14,1), KIND = wp ) … … 1895 1979 ENDDO 1896 1980 1897 DO k = nzb_max +1, nzt1981 DO k = nzb_max_l+1, nzt 1898 1982 1899 1983 v_comp(k) = v(k,j,i) + v(k,j,i-1) - gv … … 1914 1998 IF ( i == i_omp .OR. i == nxlu ) THEN 1915 1999 1916 DO k = nzb+1, nzb_max 2000 DO k = nzb+1, nzb_max_l 1917 2001 1918 2002 ibit11 = REAL( IBITS(advc_flags_1(k,j,i-1),11,1), KIND = wp ) … … 1953 2037 ENDDO 1954 2038 1955 DO k = nzb_max +1, nzt2039 DO k = nzb_max_l+1, nzt 1956 2040 1957 2041 u_comp_l = u(k,j,i) + u(k,j,i-1) - gu … … 1971 2055 !-- Now compute the fluxes tendency terms for the horizontal and 1972 2056 !-- vertical parts. 1973 DO k = nzb+1, nzb_max 2057 DO k = nzb+1, nzb_max_l 1974 2058 1975 2059 ibit11 = REAL( IBITS(advc_flags_1(k,j,i),11,1), KIND = wp ) … … 2044 2128 ENDDO 2045 2129 2046 DO k = nzb_max +1, nzt2130 DO k = nzb_max_l+1, nzt 2047 2131 2048 2132 u_comp(k) = u(k,j,i+1) + u(k,j,i) … … 2205 2289 ENDDO 2206 2290 2207 DO k = nzb+1, nz t2291 DO k = nzb+1, nzb_max_l 2208 2292 2209 2293 flux_d = flux_t(k-1) 2210 2294 diss_d = diss_t(k-1) 2295 2296 ibit11 = REAL( IBITS(advc_flags_1(k,j,i),11,1), KIND = wp ) 2297 ibit10 = REAL( IBITS(advc_flags_1(k,j,i),10,1), KIND = wp ) 2298 ibit9 = REAL( IBITS(advc_flags_1(k,j,i),9,1), KIND = wp ) 2299 2300 ibit14 = REAL( IBITS(advc_flags_1(k,j,i),14,1), KIND = wp ) 2301 ibit13 = REAL( IBITS(advc_flags_1(k,j,i),13,1), KIND = wp ) 2302 ibit12 = REAL( IBITS(advc_flags_1(k,j,i),12,1), KIND = wp ) 2303 2304 ibit17 = REAL( IBITS(advc_flags_1(k,j,i),17,1), KIND = wp ) 2305 ibit16 = REAL( IBITS(advc_flags_1(k,j,i),16,1), KIND = wp ) 2306 ibit15 = REAL( IBITS(advc_flags_1(k,j,i),15,1), KIND = wp ) 2211 2307 ! 2212 2308 !-- Calculate the divergence of the velocity field. A respective … … 2239 2335 ) * 0.5_wp 2240 2336 2241 tend(k,j,i) = tend(k,j,i) - (&2337 tend(k,j,i) = tend(k,j,i) - ( & 2242 2338 ( flux_r(k) + diss_r(k) & 2243 2339 - flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx & … … 2249 2345 ) + div * u(k,j,i) 2250 2346 2251 2252 2253 2254 2255 ! 2256 !-- 2257 !-- 2258 sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)&2347 flux_l_u(k,j,tn) = flux_r(k) 2348 diss_l_u(k,j,tn) = diss_r(k) 2349 flux_s_u(k,tn) = flux_n(k) 2350 diss_s_u(k,tn) = diss_n(k) 2351 ! 2352 !-- Statistical Evaluation of u'u'. The factor has to be applied for 2353 !-- right evaluation when gallilei_trans = .T. . 2354 sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) & 2259 2355 + ( flux_r(k) & 2260 2356 * ( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) & … … 2265 2361 ) * weight_substep(intermediate_timestep_count) 2266 2362 ! 2267 !-- 2268 sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)&2363 !-- Statistical Evaluation of w'u'. 2364 sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) & 2269 2365 + ( flux_t(k) & 2270 2366 * ( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) & … … 2275 2371 ) * weight_substep(intermediate_timestep_count) 2276 2372 ENDDO 2277 2278 sums_us2_ws_l(nzb,tn) = sums_us2_ws_l(nzb+1,tn) 2373 2374 DO k = nzb_max_l+1, nzt 2375 2376 flux_d = flux_t(k-1) 2377 diss_d = diss_t(k-1) 2378 ! 2379 !-- Calculate the divergence of the velocity field. A respective 2380 !-- correction is needed to overcome numerical instabilities introduced 2381 !-- by a not sufficient reduction of divergences near topography. 2382 div = ( ( u_comp(k) - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx & 2383 + ( v_comp(k) + gv - ( v(k,j,i) + v(k,j,i-1) ) ) * ddy & 2384 + ( w_comp(k) * rho_air_zw(k) & 2385 - w_comp(k-1) * rho_air_zw(k-1) & 2386 ) * drho_air(k) * ddzw(k) & 2387 ) * 0.5_wp 2388 2389 tend(k,j,i) = tend(k,j,i) - ( & 2390 ( flux_r(k) + diss_r(k) & 2391 - flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx & 2392 + ( flux_n(k) + diss_n(k) & 2393 - flux_s_u(k,tn) - diss_s_u(k,tn) ) * ddy & 2394 + ( ( flux_t(k) + diss_t(k) ) & 2395 - ( flux_d + diss_d ) & 2396 ) * drho_air(k) * ddzw(k) & 2397 ) + div * u(k,j,i) 2398 2399 flux_l_u(k,j,tn) = flux_r(k) 2400 diss_l_u(k,j,tn) = diss_r(k) 2401 flux_s_u(k,tn) = flux_n(k) 2402 diss_s_u(k,tn) = diss_n(k) 2403 ! 2404 !-- Statistical Evaluation of u'u'. The factor has to be applied for 2405 !-- right evaluation when gallilei_trans = .T. . 2406 sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) & 2407 + ( flux_r(k) & 2408 * ( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) & 2409 / ( u_comp(k) - gu + SIGN( 1.0E-20_wp, u_comp(k) - gu ) ) & 2410 + diss_r(k) & 2411 * ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) & 2412 / ( ABS( u_comp(k) - gu ) + 1.0E-20_wp ) & 2413 ) * weight_substep(intermediate_timestep_count) 2414 ! 2415 !-- Statistical Evaluation of w'u'. 2416 sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) & 2417 + ( flux_t(k) & 2418 * ( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) & 2419 / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) ) ) & 2420 + diss_t(k) & 2421 * ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) & 2422 / ( ABS( w_comp(k) ) + 1.0E-20_wp ) & 2423 ) * weight_substep(intermediate_timestep_count) 2424 ENDDO 2279 2425 2280 2426 … … 2296 2442 2297 2443 USE control_parameters, & 2298 ONLY: intermediate_timestep_count, u_gtrans, v_gtrans 2444 ONLY: bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, & 2445 bc_dirichlet_s, bc_radiation_l, bc_radiation_n, & 2446 bc_radiation_r, bc_radiation_s, intermediate_timestep_count, & 2447 u_gtrans, v_gtrans 2299 2448 2300 2449 USE grid_variables, & … … 2302 2451 2303 2452 USE indices, & 2304 ONLY: ny sv, nzb, nzb_max, nzt, advc_flags_12453 ONLY: nyn, nys, nysv, nxl, nxr, nzb, nzb_max, nzt, advc_flags_1 2305 2454 2306 2455 USE kinds … … 2311 2460 IMPLICIT NONE 2312 2461 2313 INTEGER(iwp) :: i !< grid index along x-direction 2314 INTEGER(iwp) :: i_omp !< leftmost index on subdomain, or in case of OpenMP, on thread 2315 INTEGER(iwp) :: j !< grid index along y-direction 2316 INTEGER(iwp) :: k !< grid index along z-direction 2317 INTEGER(iwp) :: k_mm !< k-2 index in disretization, can be modified to avoid segmentation faults 2318 INTEGER(iwp) :: k_pp !< k+2 index in disretization, can be modified to avoid segmentation faults 2319 INTEGER(iwp) :: k_ppp !< k+3 index in disretization, can be modified to avoid segmentation faults 2320 INTEGER(iwp) :: tn !< number of OpenMP thread 2462 INTEGER(iwp) :: i !< grid index along x-direction 2463 INTEGER(iwp) :: i_omp !< leftmost index on subdomain, or in case of OpenMP, on thread 2464 INTEGER(iwp) :: j !< grid index along y-direction 2465 INTEGER(iwp) :: k !< grid index along z-direction 2466 INTEGER(iwp) :: k_mm !< k-2 index in disretization, can be modified to avoid segmentation faults 2467 INTEGER(iwp) :: k_pp !< k+2 index in disretization, can be modified to avoid segmentation faults 2468 INTEGER(iwp) :: k_ppp !< k+3 index in disretization, can be modified to avoid segmentation faults 2469 INTEGER(iwp) :: nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms 2470 INTEGER(iwp) :: tn !< number of OpenMP thread 2321 2471 2322 2472 REAL(wp) :: ibit18 !< flag indicating 1st-order scheme along x-direction … … 2345 2495 REAL(wp), DIMENSION(nzb:nzt+1) :: v_comp !< advection velocity along y 2346 2496 REAL(wp), DIMENSION(nzb:nzt+1) :: w_comp !< advection velocity along z 2347 2497 ! 2498 !-- Used local modified copy of nzb_max (used to degrade order of 2499 !-- discretization) at non-cyclic boundaries. Modify only at relevant points 2500 !-- instead of the entire subdomain. This should lead to better 2501 !-- load balance between boundary and non-boundary PEs. 2502 IF( ( bc_dirichlet_l .OR. bc_radiation_l ) .AND. i <= nxl + 2 .OR. & 2503 ( bc_dirichlet_r .OR. bc_radiation_r ) .AND. i >= nxr - 2 .OR. & 2504 ( bc_dirichlet_s .OR. bc_radiation_s ) .AND. j <= nys + 2 .OR. & 2505 ( bc_dirichlet_n .OR. bc_radiation_n ) .AND. j >= nyn - 2 ) THEN 2506 nzb_max_l = nzt 2507 ELSE 2508 nzb_max_l = nzb_max 2509 END IF 2510 2348 2511 gu = 2.0_wp * u_gtrans 2349 2512 gv = 2.0_wp * v_gtrans … … 2353 2516 IF ( i == i_omp ) THEN 2354 2517 2355 DO k = nzb+1, nzb_max 2518 DO k = nzb+1, nzb_max_l 2356 2519 2357 2520 ibit20 = REAL( IBITS(advc_flags_1(k,j,i-1),20,1), KIND = wp ) … … 2392 2555 ENDDO 2393 2556 2394 DO k = nzb_max +1, nzt2557 DO k = nzb_max_l+1, nzt 2395 2558 2396 2559 u_comp(k) = u(k,j-1,i) + u(k,j,i) - gu … … 2411 2574 IF ( j == nysv ) THEN 2412 2575 2413 DO k = nzb+1, nzb_max 2576 DO k = nzb+1, nzb_max_l 2414 2577 2415 2578 ibit23 = REAL( IBITS(advc_flags_1(k,j-1,i),23,1), KIND = wp ) … … 2450 2613 ENDDO 2451 2614 2452 DO k = nzb_max +1, nzt2615 DO k = nzb_max_l+1, nzt 2453 2616 2454 2617 v_comp_l = v(k,j,i) + v(k,j-1,i) - gv … … 2468 2631 !-- Now compute the fluxes and tendency terms for the horizontal and 2469 2632 !-- verical parts. 2470 DO k = nzb+1, nzb_max 2633 DO k = nzb+1, nzb_max_l 2471 2634 2472 2635 ibit20 = REAL( IBITS(advc_flags_1(k,j,i),20,1), KIND = wp ) … … 2542 2705 ENDDO 2543 2706 2544 DO k = nzb_max +1, nzt2707 DO k = nzb_max_l+1, nzt 2545 2708 2546 2709 u_comp(k) = u(k,j-1,i+1) + u(k,j,i+1) - gu … … 2706 2869 ENDDO 2707 2870 2708 DO k = nzb+1, nz t2871 DO k = nzb+1, nzb_max_l 2709 2872 2710 2873 flux_d = flux_t(k-1) 2711 2874 diss_d = diss_t(k-1) 2875 2876 ibit20 = REAL( IBITS(advc_flags_1(k,j,i),20,1), KIND = wp ) 2877 ibit19 = REAL( IBITS(advc_flags_1(k,j,i),19,1), KIND = wp ) 2878 ibit18 = REAL( IBITS(advc_flags_1(k,j,i),18,1), KIND = wp ) 2879 2880 ibit23 = REAL( IBITS(advc_flags_1(k,j,i),23,1), KIND = wp ) 2881 ibit22 = REAL( IBITS(advc_flags_1(k,j,i),22,1), KIND = wp ) 2882 ibit21 = REAL( IBITS(advc_flags_1(k,j,i),21,1), KIND = wp ) 2883 2884 ibit26 = REAL( IBITS(advc_flags_1(k,j,i),26,1), KIND = wp ) 2885 ibit25 = REAL( IBITS(advc_flags_1(k,j,i),25,1), KIND = wp ) 2886 ibit24 = REAL( IBITS(advc_flags_1(k,j,i),24,1), KIND = wp ) 2712 2887 ! 2713 2888 !-- Calculate the divergence of the velocity field. A respective … … 2752 2927 ) + v(k,j,i) * div 2753 2928 2754 2755 2756 2757 2758 ! 2759 !-- 2760 !-- 2761 sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn)&2929 flux_l_v(k,j,tn) = flux_r(k) 2930 diss_l_v(k,j,tn) = diss_r(k) 2931 flux_s_v(k,tn) = flux_n(k) 2932 diss_s_v(k,tn) = diss_n(k) 2933 ! 2934 !-- Statistical Evaluation of v'v'. The factor has to be applied for 2935 !-- right evaluation when gallilei_trans = .T. . 2936 sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn) & 2762 2937 + ( flux_n(k) & 2763 2938 * ( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) & … … 2768 2943 ) * weight_substep(intermediate_timestep_count) 2769 2944 ! 2770 !-- 2771 sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)&2945 !-- Statistical Evaluation of w'u'. 2946 sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) & 2772 2947 + ( flux_t(k) & 2773 2948 * ( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) & … … 2779 2954 2780 2955 ENDDO 2781 sums_vs2_ws_l(nzb,tn) = sums_vs2_ws_l(nzb+1,tn) 2956 2957 DO k = nzb_max_l+1, nzt 2958 2959 flux_d = flux_t(k-1) 2960 diss_d = diss_t(k-1) 2961 ! 2962 !-- Calculate the divergence of the velocity field. A respective 2963 !-- correction is needed to overcome numerical instabilities introduced 2964 !-- by a not sufficient reduction of divergences near topography. 2965 div = ( ( u_comp(k) + gu - ( u(k,j-1,i) + u(k,j,i) ) ) * ddx & 2966 + ( v_comp(k) - ( v(k,j,i) + v(k,j-1,i) ) ) * ddy & 2967 + ( w_comp(k) * rho_air_zw(k) & 2968 - w_comp(k-1) * rho_air_zw(k-1) & 2969 ) * drho_air(k) * ddzw(k) & 2970 ) * 0.5_wp 2971 2972 tend(k,j,i) = tend(k,j,i) - ( & 2973 ( flux_r(k) + diss_r(k) & 2974 - flux_l_v(k,j,tn) - diss_l_v(k,j,tn) ) * ddx & 2975 + ( flux_n(k) + diss_n(k) & 2976 - flux_s_v(k,tn) - diss_s_v(k,tn) ) * ddy & 2977 + ( ( flux_t(k) + diss_t(k) ) & 2978 - ( flux_d + diss_d ) & 2979 ) * drho_air(k) * ddzw(k) & 2980 ) + v(k,j,i) * div 2981 2982 flux_l_v(k,j,tn) = flux_r(k) 2983 diss_l_v(k,j,tn) = diss_r(k) 2984 flux_s_v(k,tn) = flux_n(k) 2985 diss_s_v(k,tn) = diss_n(k) 2986 ! 2987 !-- Statistical Evaluation of v'v'. The factor has to be applied for 2988 !-- right evaluation when gallilei_trans = .T. . 2989 sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn) & 2990 + ( flux_n(k) & 2991 * ( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) & 2992 / ( v_comp(k) - gv + SIGN( 1.0E-20_wp, v_comp(k) - gv ) ) & 2993 + diss_n(k) & 2994 * ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) & 2995 / ( ABS( v_comp(k) - gv ) + 1.0E-20_wp ) & 2996 ) * weight_substep(intermediate_timestep_count) 2997 ! 2998 !-- Statistical Evaluation of w'u'. 2999 sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) & 3000 + ( flux_t(k) & 3001 * ( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) & 3002 / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) ) ) & 3003 + diss_t(k) & 3004 * ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) & 3005 / ( ABS( w_comp(k) ) + 1.0E-20_wp ) & 3006 ) * weight_substep(intermediate_timestep_count) 3007 3008 ENDDO 2782 3009 2783 3010 … … 2793 3020 SUBROUTINE advec_w_ws_ij( i, j, i_omp, tn ) 2794 3021 2795 USE arrays_3d, &2796 ONLY: ddzu, diss_l_w, diss_s_w, flux_l_w, flux_s_w, tend, u, v, w, &3022 USE arrays_3d, & 3023 ONLY: ddzu, diss_l_w, diss_s_w, flux_l_w, flux_s_w, tend, u, v, w, & 2797 3024 drho_air_zw, rho_air 2798 3025 2799 USE control_parameters, & 2800 ONLY: intermediate_timestep_count, u_gtrans, v_gtrans 2801 2802 USE grid_variables, & 3026 USE control_parameters, & 3027 ONLY: bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, & 3028 bc_dirichlet_s, bc_radiation_l, bc_radiation_n, & 3029 bc_radiation_r, bc_radiation_s,intermediate_timestep_count, & 3030 u_gtrans, v_gtrans 3031 3032 USE grid_variables, & 2803 3033 ONLY: ddx, ddy 2804 3034 2805 USE indices, & 2806 ONLY: nys, nzb, nzb_max, nzt, advc_flags_1, advc_flags_2 3035 USE indices, & 3036 ONLY: nyn, nys, nxl, nxr, nzb, nzb_max, nzt, advc_flags_1, & 3037 advc_flags_2 2807 3038 2808 3039 USE kinds 2809 3040 2810 USE statistics, &3041 USE statistics, & 2811 3042 ONLY: hom, sums_ws2_ws_l, weight_substep 2812 3043 2813 3044 IMPLICIT NONE 2814 3045 2815 INTEGER(iwp) :: i !< grid index along x-direction 2816 INTEGER(iwp) :: i_omp !< leftmost index on subdomain, or in case of OpenMP, on thread 2817 INTEGER(iwp) :: j !< grid index along y-direction 2818 INTEGER(iwp) :: k !< grid index along z-direction 2819 INTEGER(iwp) :: k_mm !< k-2 index in disretization, can be modified to avoid segmentation faults 2820 INTEGER(iwp) :: k_pp !< k+2 index in disretization, can be modified to avoid segmentation faults 2821 INTEGER(iwp) :: k_ppp !< k+3 index in disretization, can be modified to avoid segmentation faults 2822 INTEGER(iwp) :: tn !< number of OpenMP thread 3046 INTEGER(iwp) :: i !< grid index along x-direction 3047 INTEGER(iwp) :: i_omp !< leftmost index on subdomain, or in case of OpenMP, on thread 3048 INTEGER(iwp) :: j !< grid index along y-direction 3049 INTEGER(iwp) :: k !< grid index along z-direction 3050 INTEGER(iwp) :: k_mm !< k-2 index in disretization, can be modified to avoid segmentation faults 3051 INTEGER(iwp) :: k_pp !< k+2 index in disretization, can be modified to avoid segmentation faults 3052 INTEGER(iwp) :: k_ppp !< k+3 index in disretization, can be modified to avoid segmentation faults 3053 INTEGER(iwp) :: nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms 3054 INTEGER(iwp) :: tn !< number of OpenMP thread 2823 3055 2824 3056 REAL(wp) :: ibit27 !< flag indicating 1st-order scheme along x-direction … … 2846 3078 REAL(wp), DIMENSION(nzb:nzt+1) :: v_comp !< advection velocity along y 2847 3079 REAL(wp), DIMENSION(nzb:nzt+1) :: w_comp !< advection velocity along z 2848 3080 ! 3081 !-- Used local modified copy of nzb_max (used to degrade order of 3082 !-- discretization) at non-cyclic boundaries. Modify only at relevant points 3083 !-- instead of the entire subdomain. This should lead to better 3084 !-- load balance between boundary and non-boundary PEs. 3085 IF( ( bc_dirichlet_l .OR. bc_radiation_l ) .AND. i <= nxl + 2 .OR. & 3086 ( bc_dirichlet_r .OR. bc_radiation_r ) .AND. i >= nxr - 2 .OR. & 3087 ( bc_dirichlet_s .OR. bc_radiation_s ) .AND. j <= nys + 2 .OR. & 3088 ( bc_dirichlet_n .OR. bc_radiation_n ) .AND. j >= nyn - 2 ) THEN 3089 nzb_max_l = nzt 3090 ELSE 3091 nzb_max_l = nzb_max 3092 END IF 3093 2849 3094 gu = 2.0_wp * u_gtrans 2850 3095 gv = 2.0_wp * v_gtrans 2851 2852 3096 ! 2853 3097 !-- Compute southside fluxes for the respective boundary. 2854 3098 IF ( j == nys ) THEN 2855 3099 2856 DO k = nzb+1, nzb_max 3100 DO k = nzb+1, nzb_max_l 2857 3101 ibit32 = REAL( IBITS(advc_flags_2(k,j-1,i),0,1), KIND = wp ) 2858 3102 ibit31 = REAL( IBITS(advc_flags_1(k,j-1,i),31,1), KIND = wp ) … … 2892 3136 ENDDO 2893 3137 2894 DO k = nzb_max +1, nzt3138 DO k = nzb_max_l+1, nzt 2895 3139 2896 3140 v_comp(k) = v(k+1,j,i) + v(k,j,i) - gv … … 2911 3155 IF ( i == i_omp ) THEN 2912 3156 2913 DO k = nzb+1, nzb_max 3157 DO k = nzb+1, nzb_max_l 2914 3158 2915 3159 ibit29 = REAL( IBITS(advc_flags_1(k,j,i-1),29,1), KIND = wp ) … … 2950 3194 ENDDO 2951 3195 2952 DO k = nzb_max +1, nzt3196 DO k = nzb_max_l+1, nzt 2953 3197 2954 3198 u_comp(k) = u(k+1,j,i) + u(k,j,i) - gu … … 2968 3212 !-- Now compute the fluxes and tendency terms for the horizontal 2969 3213 !-- and vertical parts. 2970 DO k = nzb+1, nzb_max 3214 DO k = nzb+1, nzb_max_l 2971 3215 2972 3216 ibit29 = REAL( IBITS(advc_flags_1(k,j,i),29,1), KIND = wp ) … … 3041 3285 ENDDO 3042 3286 3043 DO k = nzb_max +1, nzt3287 DO k = nzb_max_l+1, nzt 3044 3288 3045 3289 u_comp(k) = u(k+1,j,i+1) + u(k,j,i+1) - gu … … 3206 3450 ENDDO 3207 3451 3208 DO k = nzb+1, nz t3452 DO k = nzb+1, nzb_max_l 3209 3453 3210 3454 flux_d = flux_t(k-1) 3211 3455 diss_d = diss_t(k-1) 3456 3457 ibit29 = REAL( IBITS(advc_flags_1(k,j,i),29,1), KIND = wp ) 3458 ibit28 = REAL( IBITS(advc_flags_1(k,j,i),28,1), KIND = wp ) 3459 ibit27 = REAL( IBITS(advc_flags_1(k,j,i),27,1), KIND = wp ) 3460 3461 ibit32 = REAL( IBITS(advc_flags_2(k,j,i),0,1), KIND = wp ) 3462 ibit31 = REAL( IBITS(advc_flags_1(k,j,i),31,1), KIND = wp ) 3463 ibit30 = REAL( IBITS(advc_flags_1(k,j,i),30,1), KIND = wp ) 3464 3465 ibit35 = REAL( IBITS(advc_flags_2(k,j,i),3,1), KIND = wp ) 3466 ibit34 = REAL( IBITS(advc_flags_2(k,j,i),2,1), KIND = wp ) 3467 ibit33 = REAL( IBITS(advc_flags_2(k,j,i),1,1), KIND = wp ) 3212 3468 ! 3213 3469 !-- Calculate the divergence of the velocity field. A respective … … 3267 3523 3268 3524 ENDDO 3269 3525 3526 DO k = nzb_max_l+1, nzt 3527 3528 flux_d = flux_t(k-1) 3529 diss_d = diss_t(k-1) 3530 ! 3531 !-- Calculate the divergence of the velocity field. A respective 3532 !-- correction is needed to overcome numerical instabilities introduced 3533 !-- by a not sufficient reduction of divergences near topography. 3534 div = ( ( u_comp(k) + gu - ( u(k+1,j,i) + u(k,j,i) ) ) * ddx & 3535 + ( v_comp(k) + gv - ( v(k+1,j,i) + v(k,j,i) ) ) * ddy & 3536 + ( w_comp(k) * rho_air(k+1) & 3537 - ( w(k,j,i) + w(k-1,j,i) ) * rho_air(k) & 3538 ) * drho_air_zw(k) * ddzu(k+1) & 3539 ) * 0.5_wp 3540 3541 tend(k,j,i) = tend(k,j,i) - ( & 3542 ( flux_r(k) + diss_r(k) & 3543 - flux_l_w(k,j,tn) - diss_l_w(k,j,tn) ) * ddx & 3544 + ( flux_n(k) + diss_n(k) & 3545 - flux_s_w(k,tn) - diss_s_w(k,tn) ) * ddy & 3546 + ( ( flux_t(k) + diss_t(k) ) & 3547 - ( flux_d + diss_d ) & 3548 ) * drho_air_zw(k) * ddzu(k+1) & 3549 ) + div * w(k,j,i) 3550 3551 flux_l_w(k,j,tn) = flux_r(k) 3552 diss_l_w(k,j,tn) = diss_r(k) 3553 flux_s_w(k,tn) = flux_n(k) 3554 diss_s_w(k,tn) = diss_n(k) 3555 ! 3556 !-- Statistical Evaluation of w'w'. 3557 sums_ws2_ws_l(k,tn) = sums_ws2_ws_l(k,tn) & 3558 + ( flux_t(k) & 3559 * ( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) & 3560 / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) ) ) & 3561 + diss_t(k) & 3562 * ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) & 3563 / ( ABS( w_comp(k) ) + 1.0E-20_wp ) & 3564 ) * weight_substep(intermediate_timestep_count) 3565 3566 ENDDO 3270 3567 3271 3568 END SUBROUTINE advec_w_ws_ij … … 3284 3581 3285 3582 USE control_parameters, & 3286 ONLY: intermediate_timestep_count, u_gtrans, v_gtrans 3583 ONLY: bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, & 3584 bc_dirichlet_s, bc_radiation_l, bc_radiation_n, & 3585 bc_radiation_r, bc_radiation_s, intermediate_timestep_count, & 3586 u_gtrans, v_gtrans 3287 3587 3288 3588 USE grid_variables, & … … 3310 3610 INTEGER(iwp) :: sk_num !< integer identifier, used for assign fluxes to the correct dimension in the analysis array 3311 3611 3312 INTEGER(iwp) :: i !< grid index along x-direction 3313 INTEGER(iwp) :: j !< grid index along y-direction 3314 INTEGER(iwp) :: k !< grid index along z-direction 3315 INTEGER(iwp) :: k_mm !< k-2 index in disretization, can be modified to avoid segmentation faults 3316 INTEGER(iwp) :: k_pp !< k+2 index in disretization, can be modified to avoid segmentation faults 3317 INTEGER(iwp) :: k_ppp !< k+3 index in disretization, can be modified to avoid segmentation faults 3318 INTEGER(iwp) :: tn = 0 !< number of OpenMP thread 3612 INTEGER(iwp) :: i !< grid index along x-direction 3613 INTEGER(iwp) :: j !< grid index along y-direction 3614 INTEGER(iwp) :: k !< grid index along z-direction 3615 INTEGER(iwp) :: k_mm !< k-2 index in disretization, can be modified to avoid segmentation faults 3616 INTEGER(iwp) :: k_pp !< k+2 index in disretization, can be modified to avoid segmentation faults 3617 INTEGER(iwp) :: k_ppp !< k+3 index in disretization, can be modified to avoid segmentation faults 3618 INTEGER(iwp) :: nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms 3619 INTEGER(iwp) :: tn = 0 !< number of OpenMP thread 3319 3620 3320 3621 ! … … 3378 3679 #endif 3379 3680 3681 ! 3682 !-- Set local version of nzb_max. At non-cyclic boundaries the order of the 3683 !-- advection need to be degraded near the boundary. Please note, in contrast 3684 !-- to the cache-optimized routines, nzb_max_l is set constantly for the 3685 !-- entire subdomain, in order to avoid unsymmetric loops which might be 3686 !-- an issue for GPUs. 3687 IF( bc_dirichlet_l .OR. bc_radiation_l .OR. & 3688 bc_dirichlet_r .OR. bc_radiation_r .OR. & 3689 bc_dirichlet_s .OR. bc_radiation_s .OR. & 3690 bc_dirichlet_n .OR. bc_radiation_n ) THEN 3691 nzb_max_l = nzt 3692 ELSE 3693 nzb_max_l = nzb_max 3694 END IF 3695 3380 3696 SELECT CASE ( sk_char ) 3381 3697 … … 3407 3723 DO j = nys, nyn 3408 3724 3409 DO k = nzb+1, nzb_max 3725 DO k = nzb+1, nzb_max_l 3410 3726 3411 3727 ibit2 = REAL( IBITS(advc_flags_1(k,j,i-1),2,1), KIND = wp ) … … 3446 3762 ENDDO 3447 3763 3448 DO k = nzb_max +1, nzt3764 DO k = nzb_max_l+1, nzt 3449 3765 3450 3766 u_comp = u(k,j,i) - u_gtrans + u_stokes_zu(k) … … 3480 3796 !$ACC PRESENT(drho_air, rho_air_zw, ddzw) & 3481 3797 !$ACC PRESENT(tend) & 3482 !$ACC PRESENT(hom(nzb+1:nzb_max ,1,1:3,0)) &3798 !$ACC PRESENT(hom(nzb+1:nzb_max_l,1,1:3,0)) & 3483 3799 !$ACC PRESENT(weight_substep(intermediate_timestep_count)) & 3484 3800 !$ACC PRESENT(sums_wspts_ws_l, sums_wssas_ws_l) & … … 3491 3807 #ifndef _OPENACC 3492 3808 j = nys 3493 DO k = nzb+1, nzb_max 3809 DO k = nzb+1, nzb_max_l 3494 3810 3495 3811 ibit5 = REAL( IBITS(advc_flags_1(k,j-1,i),5,1), KIND = wp ) … … 3531 3847 ! 3532 3848 !-- Above to the top of the highest topography. No degradation necessary. 3533 DO k = nzb_max +1, nzt3849 DO k = nzb_max_l+1, nzt 3534 3850 3535 3851 v_comp = v(k,j,i) - v_gtrans + v_stokes_zu(k) … … 3553 3869 diss_d = 0.0_wp 3554 3870 3555 DO k = nzb+1, nzb_max 3871 DO k = nzb+1, nzb_max_l 3556 3872 3557 3873 ibit2 = REAL( IBITS(advc_flags_1(k,j,i),2,1), KIND = wp ) … … 3895 4211 ENDDO 3896 4212 3897 DO k = nzb_max +1, nzt4213 DO k = nzb_max_l+1, nzt 3898 4214 3899 4215 u_comp = u(k,j,i+1) - u_gtrans + u_stokes_zu(k) … … 3968 4284 3969 4285 3970 flux_t = w(k,j,i) * rho_air_zw(k) * ( &4286 flux_t = w(k,j,i) * rho_air_zw(k) * ( & 3971 4287 ( 37.0_wp * ibit8 * adv_sca_5 & 3972 4288 + 7.0_wp * ibit7 * adv_sca_3 & … … 3982 4298 ) 3983 4299 3984 diss_t = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( &4300 diss_t = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( & 3985 4301 ( 10.0_wp * ibit8 * adv_sca_5 & 3986 4302 + 3.0_wp * ibit7 * adv_sca_3 & … … 4140 4456 4141 4457 USE control_parameters, & 4142 ONLY: intermediate_timestep_count, u_gtrans, v_gtrans 4458 ONLY: bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, & 4459 bc_dirichlet_s, bc_radiation_l, bc_radiation_n, & 4460 bc_radiation_r, bc_radiation_s, intermediate_timestep_count, & 4461 u_gtrans, v_gtrans 4143 4462 4144 4463 USE grid_variables, & … … 4155 4474 IMPLICIT NONE 4156 4475 4157 INTEGER(iwp) :: i !< grid index along x-direction 4158 INTEGER(iwp) :: j !< grid index along y-direction 4159 INTEGER(iwp) :: k !< grid index along z-direction 4160 INTEGER(iwp) :: k_mm !< k-2 index in disretization, can be modified to avoid segmentation faults 4161 INTEGER(iwp) :: k_pp !< k+2 index in disretization, can be modified to avoid segmentation faults 4162 INTEGER(iwp) :: k_ppp !< k+3 index in disretization, can be modified to avoid segmentation faults 4163 INTEGER(iwp) :: tn = 0 !< number of OpenMP thread 4476 INTEGER(iwp) :: i !< grid index along x-direction 4477 INTEGER(iwp) :: j !< grid index along y-direction 4478 INTEGER(iwp) :: k !< grid index along z-direction 4479 INTEGER(iwp) :: k_mm !< k-2 index in disretization, can be modified to avoid segmentation faults 4480 INTEGER(iwp) :: k_pp !< k+2 index in disretization, can be modified to avoid segmentation faults 4481 INTEGER(iwp) :: k_ppp !< k+3 index in disretization, can be modified to avoid segmentation faults 4482 INTEGER(iwp) :: nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms 4483 INTEGER(iwp) :: tn = 0 !< number of OpenMP thread 4164 4484 4165 4485 REAL(wp) :: ibit9 !< flag indicating 1st-order scheme along x-direction … … 4217 4537 REAL(wp) :: u_comp_l !< 4218 4538 #endif 4539 ! 4540 !-- Set local version of nzb_max. At non-cyclic boundaries the order of the 4541 !-- advection need to be degraded near the boundary. Please note, in contrast 4542 !-- to the cache-optimized routines, nzb_max_l is set constantly for the 4543 !-- entire subdomain, in order to avoid unsymmetric loops which might be 4544 !-- an issue for GPUs. 4545 IF( bc_dirichlet_l .OR. bc_radiation_l .OR. & 4546 bc_dirichlet_r .OR. bc_radiation_r .OR. & 4547 bc_dirichlet_s .OR. bc_radiation_s .OR. & 4548 bc_dirichlet_n .OR. bc_radiation_n ) THEN 4549 nzb_max_l = nzt 4550 ELSE 4551 nzb_max_l = nzb_max 4552 END IF 4219 4553 4220 4554 gu = 2.0_wp * u_gtrans … … 4226 4560 i = nxlu 4227 4561 DO j = nys, nyn 4228 DO k = nzb+1, nzb_max 4562 DO k = nzb+1, nzb_max_l 4229 4563 4230 4564 ibit11 = REAL( IBITS(advc_flags_1(k,j,i-1),11,1), KIND = wp ) … … 4265 4599 ENDDO 4266 4600 4267 DO k = nzb_max +1, nzt4601 DO k = nzb_max_l+1, nzt 4268 4602 4269 4603 u_comp = u(k,j,i) + u(k,j,i-1) - gu … … 4295 4629 !$ACC PRESENT(drho_air, rho_air_zw, ddzw) & 4296 4630 !$ACC PRESENT(tend) & 4297 !$ACC PRESENT(hom(nzb+1:nzb_max ,1,1:3,0)) &4631 !$ACC PRESENT(hom(nzb+1:nzb_max_l,1,1:3,0)) & 4298 4632 !$ACC PRESENT(weight_substep(intermediate_timestep_count)) & 4299 4633 !$ACC PRESENT(sums_us2_ws_l, sums_wsus_ws_l) … … 4303 4637 !-- The following loop computes the fluxes for the south boundary points 4304 4638 j = nys 4305 DO k = nzb+1, nzb_max 4639 DO k = nzb+1, nzb_max_l 4306 4640 4307 4641 ibit14 = REAL( IBITS(advc_flags_1(k,j-1,i),14,1), KIND = wp ) … … 4342 4676 ENDDO 4343 4677 4344 DO k = nzb_max +1, nzt4678 DO k = nzb_max_l+1, nzt 4345 4679 4346 4680 v_comp = v(k,j,i) + v(k,j,i-1) - gv … … 4364 4698 diss_d = 0.0_wp 4365 4699 4366 DO k = nzb+1, nzb_max 4700 DO k = nzb+1, nzb_max_l 4367 4701 4368 4702 ibit11 = REAL( IBITS(advc_flags_1(k,j,i),11,1), KIND = wp ) … … 4636 4970 ENDDO 4637 4971 4638 DO k = nzb_max +1, nzt4972 DO k = nzb_max_l+1, nzt 4639 4973 4640 4974 u_comp = u(k,j,i+1) + u(k,j,i) … … 4787 5121 ENDDO 4788 5122 ENDDO 4789 sums_us2_ws_l(nzb,tn) = sums_us2_ws_l(nzb+1,tn)4790 4791 5123 4792 5124 END SUBROUTINE advec_u_ws … … 4804 5136 4805 5137 USE control_parameters, & 4806 ONLY: intermediate_timestep_count, u_gtrans, v_gtrans 5138 ONLY: bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, & 5139 bc_dirichlet_s, bc_radiation_l, bc_radiation_n, & 5140 bc_radiation_r, bc_radiation_s, intermediate_timestep_count, & 5141 u_gtrans, v_gtrans 4807 5142 4808 5143 USE grid_variables, & … … 4820 5155 4821 5156 4822 INTEGER(iwp) :: i !< grid index along x-direction 4823 INTEGER(iwp) :: j !< grid index along y-direction 4824 INTEGER(iwp) :: k !< grid index along z-direction 4825 INTEGER(iwp) :: k_mm !< k-2 index in disretization, can be modified to avoid segmentation faults 4826 INTEGER(iwp) :: k_pp !< k+2 index in disretization, can be modified to avoid segmentation faults 4827 INTEGER(iwp) :: k_ppp !< k+3 index in disretization, can be modified to avoid segmentation faults 4828 INTEGER(iwp) :: tn = 0 !< number of OpenMP thread 5157 INTEGER(iwp) :: i !< grid index along x-direction 5158 INTEGER(iwp) :: j !< grid index along y-direction 5159 INTEGER(iwp) :: k !< grid index along z-direction 5160 INTEGER(iwp) :: k_mm !< k-2 index in disretization, can be modified to avoid segmentation faults 5161 INTEGER(iwp) :: k_pp !< k+2 index in disretization, can be modified to avoid segmentation faults 5162 INTEGER(iwp) :: k_ppp !< k+3 index in disretization, can be modified to avoid segmentation faults 5163 INTEGER(iwp) :: nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms 5164 INTEGER(iwp) :: tn = 0 !< number of OpenMP thread 4829 5165 4830 5166 REAL(wp) :: ibit18 !< flag indicating 1st-order scheme along x-direction … … 4882 5218 REAL(wp) :: v_comp_s !< 4883 5219 #endif 4884 5220 ! 5221 !-- Set local version of nzb_max. At non-cyclic boundaries the order of the 5222 !-- advection need to be degraded near the boundary. Please note, in contrast 5223 !-- to the cache-optimized routines, nzb_max_l is set constantly for the 5224 !-- entire subdomain, in order to avoid unsymmetric loops which might be 5225 !-- an issue for GPUs. 5226 IF( bc_dirichlet_l .OR. bc_radiation_l .OR. & 5227 bc_dirichlet_r .OR. bc_radiation_r .OR. & 5228 bc_dirichlet_s .OR. bc_radiation_s .OR. & 5229 bc_dirichlet_n .OR. bc_radiation_n ) THEN 5230 nzb_max_l = nzt 5231 ELSE 5232 nzb_max_l = nzb_max 5233 END IF 5234 4885 5235 gu = 2.0_wp * u_gtrans 4886 5236 gv = 2.0_wp * v_gtrans … … 4891 5241 i = nxl 4892 5242 DO j = nysv, nyn 4893 DO k = nzb+1, nzb_max 5243 DO k = nzb+1, nzb_max_l 4894 5244 4895 5245 ibit20 = REAL( IBITS(advc_flags_1(k,j,i-1),20,1), KIND = wp ) … … 4930 5280 ENDDO 4931 5281 4932 DO k = nzb_max +1, nzt5282 DO k = nzb_max_l+1, nzt 4933 5283 4934 5284 u_comp = u(k,j-1,i) + u(k,j,i) - gu … … 4961 5311 !$ACC PRESENT(drho_air, rho_air_zw, ddzw) & 4962 5312 !$ACC PRESENT(tend) & 4963 !$ACC PRESENT(hom(nzb+1:nzb_max ,1,2:3,0)) &5313 !$ACC PRESENT(hom(nzb+1:nzb_max_l,1,2:3,0)) & 4964 5314 !$ACC PRESENT(weight_substep(intermediate_timestep_count)) & 4965 5315 !$ACC PRESENT(sums_vs2_ws_l, sums_wsvs_ws_l) … … 4968 5318 #ifndef _OPENACC 4969 5319 j = nysv 4970 DO k = nzb+1, nzb_max 5320 DO k = nzb+1, nzb_max_l 4971 5321 4972 5322 ibit23 = REAL( IBITS(advc_flags_1(k,j-1,i),23,1), KIND = wp ) … … 5007 5357 ENDDO 5008 5358 5009 DO k = nzb_max +1, nzt5359 DO k = nzb_max_l+1, nzt 5010 5360 5011 5361 v_comp = v(k,j,i) + v(k,j-1,i) - gv … … 5027 5377 diss_d = 0.0_wp 5028 5378 5029 DO k = nzb+1, nzb_max 5379 DO k = nzb+1, nzb_max_l 5030 5380 5031 5381 ibit20 = REAL( IBITS(advc_flags_1(k,j,i),20,1), KIND = wp ) … … 5302 5652 ENDDO 5303 5653 5304 DO k = nzb_max +1, nzt5654 DO k = nzb_max_l+1, nzt 5305 5655 5306 5656 u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu … … 5457 5807 ENDDO 5458 5808 ENDDO 5459 !$ACC UPDATE HOST(sums_vs2_ws_l(nzb+1,tn))5460 sums_vs2_ws_l(nzb,tn) = sums_vs2_ws_l(nzb+1,tn)5461 !$ACC UPDATE DEVICE(sums_vs2_ws_l(nzb,tn))5462 5463 5809 5464 5810 END SUBROUTINE advec_v_ws … … 5476 5822 5477 5823 USE control_parameters, & 5478 ONLY: intermediate_timestep_count, u_gtrans, v_gtrans 5824 ONLY: bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, & 5825 bc_dirichlet_s, bc_radiation_l, bc_radiation_n, & 5826 bc_radiation_r, bc_radiation_s, intermediate_timestep_count, & 5827 u_gtrans, v_gtrans 5479 5828 5480 5829 USE grid_variables, & … … 5492 5841 IMPLICIT NONE 5493 5842 5494 INTEGER(iwp) :: i !< grid index along x-direction 5495 INTEGER(iwp) :: j !< grid index along y-direction 5496 INTEGER(iwp) :: k !< grid index along z-direction 5497 INTEGER(iwp) :: k_mm !< k-2 index in disretization, can be modified to avoid segmentation faults 5498 INTEGER(iwp) :: k_pp !< k+2 index in disretization, can be modified to avoid segmentation faults 5499 INTEGER(iwp) :: k_ppp !< k+3 index in disretization, can be modified to avoid segmentation faults 5500 INTEGER(iwp) :: tn = 0 !< number of OpenMP thread 5843 INTEGER(iwp) :: i !< grid index along x-direction 5844 INTEGER(iwp) :: j !< grid index along y-direction 5845 INTEGER(iwp) :: k !< grid index along z-direction 5846 INTEGER(iwp) :: k_mm !< k-2 index in disretization, can be modified to avoid segmentation faults 5847 INTEGER(iwp) :: k_pp !< k+2 index in disretization, can be modified to avoid segmentation faults 5848 INTEGER(iwp) :: k_ppp !< k+3 index in disretization, can be modified to avoid segmentation faults 5849 INTEGER(iwp) :: nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms 5850 INTEGER(iwp) :: tn = 0 !< number of OpenMP thread 5501 5851 5502 5852 REAL(wp) :: ibit27 !< flag indicating 1st-order scheme along x-direction … … 5555 5905 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local_w !< discretized 6th-order flux at leftward-side of the grid box 5556 5906 #endif 5557 5907 ! 5908 !-- Set local version of nzb_max. At non-cyclic boundaries the order of the 5909 !-- advection need to be degraded near the boundary. Please note, in contrast 5910 !-- to the cache-optimized routines, nzb_max_l is set constantly for the 5911 !-- entire subdomain, in order to avoid unsymmetric loops which might be 5912 !-- an issue for GPUs. 5913 IF( bc_dirichlet_l .OR. bc_radiation_l .OR. & 5914 bc_dirichlet_r .OR. bc_radiation_r .OR. & 5915 bc_dirichlet_s .OR. bc_radiation_s .OR. & 5916 bc_dirichlet_n .OR. bc_radiation_n ) THEN 5917 nzb_max_l = nzt 5918 ELSE 5919 nzb_max_l = nzb_max 5920 END IF 5921 5558 5922 gu = 2.0_wp * u_gtrans 5559 5923 gv = 2.0_wp * v_gtrans … … 5564 5928 i = nxl 5565 5929 DO j = nys, nyn 5566 DO k = nzb+1, nzb_max 5930 DO k = nzb+1, nzb_max_l 5567 5931 5568 5932 ibit29 = REAL( IBITS(advc_flags_1(k,j,i-1),29,1), KIND = wp ) … … 5603 5967 ENDDO 5604 5968 5605 DO k = nzb_max +1, nzt5969 DO k = nzb_max_l+1, nzt 5606 5970 5607 5971 u_comp = u(k+1,j,i) + u(k,j,i) - gu … … 5634 5998 !$ACC PRESENT(rho_air, drho_air_zw, ddzu) & 5635 5999 !$ACC PRESENT(tend) & 5636 !$ACC PRESENT(hom(nzb+1:nzb_max ,1,3,0)) &6000 !$ACC PRESENT(hom(nzb+1:nzb_max_l,1,3,0)) & 5637 6001 !$ACC PRESENT(weight_substep(intermediate_timestep_count)) & 5638 !$ACC PRESENT(sums_ws2_ws_l(nzb+1:nzb_max ,0))6002 !$ACC PRESENT(sums_ws2_ws_l(nzb+1:nzb_max_l,0)) 5639 6003 DO i = nxl, nxr 5640 6004 5641 6005 #ifndef _OPENACC 5642 6006 j = nys 5643 DO k = nzb+1, nzb_max 6007 DO k = nzb+1, nzb_max_l 5644 6008 5645 6009 ibit32 = REAL( IBITS(advc_flags_2(k,j-1,i),0,1), KIND = wp ) … … 5680 6044 ENDDO 5681 6045 5682 DO k = nzb_max +1, nzt6046 DO k = nzb_max_l+1, nzt 5683 6047 5684 6048 v_comp = v(k+1,j,i) + v(k,j,i) - gv … … 5706 6070 diss_d = -ABS(w_comp) * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1 5707 6071 5708 DO k = nzb+1, nzb_max 6072 DO k = nzb+1, nzb_max_l 5709 6073 5710 6074 ibit29 = REAL( IBITS(advc_flags_1(k,j,i),29,1), KIND = wp ) … … 5966 6330 ENDDO 5967 6331 5968 DO k = nzb_max +1, nzt6332 DO k = nzb_max_l+1, nzt 5969 6333 5970 6334 u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu
Note: See TracChangeset
for help on using the changeset viewer.
