- Timestamp:
- Jan 9, 2018 5:44:02 PM (7 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/advec_ws.f90 (modified) (79 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/advec_ws.f90
r2718 r2731 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Enable loop vectorization by splitting the k-loop 28 ! 29 ! 2718 2018-01-02 08:49:38Z maronga 27 30 ! Corrected "Former revisions" section 28 31 ! … … 1353 1356 flux_t(0) = 0.0_wp 1354 1357 diss_t(0) = 0.0_wp 1355 flux_d = 0.0_wp1356 diss_d = 0.0_wp1357 1358 ! 1358 1359 !-- Now compute the fluxes and tendency terms for the horizontal and … … 1473 1474 ( sk(k_ppp,j,i) - sk(k_mm,j,i) ) & 1474 1475 ) 1476 ENDDO 1477 1478 DO k = nzb+1, nzb_max 1479 1480 flux_d = flux_t(k-1) 1481 diss_d = diss_t(k-1) 1475 1482 ! 1476 1483 !-- Calculate the divergence of the velocity field. A respective … … 1513 1520 swap_flux_x_local(k,j,tn) = flux_r(k) 1514 1521 swap_diss_x_local(k,j,tn) = diss_r(k) 1515 flux_d = flux_t(k)1516 diss_d = diss_t(k)1517 1522 1518 1523 ENDDO … … 1581 1586 ( sk(k_ppp,j,i) - sk(k_mm,j,i) ) & 1582 1587 ) 1588 1589 ENDDO 1590 1591 DO k = nzb_max+1, nzt 1592 1593 flux_d = flux_t(k-1) 1594 diss_d = diss_t(k-1) 1595 1583 1596 ! 1584 1597 !-- Calculate the divergence of the velocity field. A respective … … 1606 1619 swap_flux_x_local(k,j,tn) = flux_r(k) 1607 1620 swap_diss_x_local(k,j,tn) = diss_r(k) 1608 flux_d = flux_t(k)1609 diss_d = diss_t(k)1610 1621 1611 1622 ENDDO … … 1765 1776 REAL(wp) :: gv !< 1766 1777 REAL(wp) :: u_comp_l !< 1767 REAL(wp) :: v_comp !<1768 REAL(wp) :: w_comp !<1769 1778 1770 1779 REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !< … … 1775 1784 REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !< 1776 1785 REAL(wp), DIMENSION(nzb:nzt+1) :: u_comp !< 1786 REAL(wp), DIMENSION(nzb:nzt+1) :: v_comp !< 1787 REAL(wp), DIMENSION(nzb:nzt+1) :: w_comp !< 1777 1788 1778 1789 gu = 2.0_wp * u_gtrans … … 1788 1799 ibit12 = IBITS(advc_flags_1(k,j-1,i),12,1) 1789 1800 1790 v_comp = v(k,j,i) + v(k,j,i-1) - gv1791 flux_s_u(k,tn) = v_comp * (&1801 v_comp(k) = v(k,j,i) + v(k,j,i-1) - gv 1802 flux_s_u(k,tn) = v_comp(k) * ( & 1792 1803 ( 37.0_wp * ibit14 * adv_mom_5 & 1793 1804 + 7.0_wp * ibit13 * adv_mom_3 & … … 1802 1813 ) * & 1803 1814 ( u(k,j+2,i) + u(k,j-3,i) ) & 1804 )1805 1806 diss_s_u(k,tn) = - ABS ( v_comp ) * (&1815 ) 1816 1817 diss_s_u(k,tn) = - ABS ( v_comp(k) ) * ( & 1807 1818 ( 10.0_wp * ibit14 * adv_mom_5 & 1808 1819 + 3.0_wp * ibit13 * adv_mom_3 & … … 1817 1828 ) * & 1818 1829 ( u(k,j+2,i) - u(k,j-3,i) ) & 1819 )1830 ) 1820 1831 1821 1832 ENDDO … … 1823 1834 DO k = nzb_max+1, nzt 1824 1835 1825 v_comp 1826 flux_s_u(k,tn) = v_comp * (&1836 v_comp(k) = v(k,j,i) + v(k,j,i-1) - gv 1837 flux_s_u(k,tn) = v_comp(k) * ( & 1827 1838 37.0_wp * ( u(k,j,i) + u(k,j-1,i) ) & 1828 1839 - 8.0_wp * ( u(k,j+1,i) + u(k,j-2,i) ) & 1829 1840 + ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5 1830 diss_s_u(k,tn) = - ABS(v_comp ) * (&1841 diss_s_u(k,tn) = - ABS(v_comp(k)) * ( & 1831 1842 10.0_wp * ( u(k,j,i) - u(k,j-1,i) ) & 1832 1843 - 5.0_wp * ( u(k,j+1,i) - u(k,j-2,i) ) & … … 1897 1908 flux_t(0) = 0.0_wp 1898 1909 diss_t(0) = 0.0_wp 1899 flux_d = 0.0_wp 1900 diss_d = 0.0_wp 1910 w_comp(0) = 0.0_wp 1901 1911 ! 1902 1912 !-- Now compute the fluxes tendency terms for the horizontal and … … 1943 1953 ibit12 = IBITS(advc_flags_1(k,j,i),12,1) 1944 1954 1945 v_comp 1946 flux_n(k) = v_comp * (&1955 v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv 1956 flux_n(k) = v_comp(k) * ( & 1947 1957 ( 37.0_wp * ibit14 * adv_mom_5 & 1948 1958 + 7.0_wp * ibit13 * adv_mom_3 & … … 1957 1967 ) * & 1958 1968 ( u(k,j+3,i) + u(k,j-2,i) ) & 1959 )1960 1961 diss_n(k) = - ABS ( v_comp ) * (&1969 ) 1970 1971 diss_n(k) = - ABS ( v_comp(k) ) * ( & 1962 1972 ( 10.0_wp * ibit14 * adv_mom_5 & 1963 1973 + 3.0_wp * ibit13 * adv_mom_3 & … … 1972 1982 ) * & 1973 1983 ( u(k,j+3,i) - u(k,j-2,i) ) & 1974 )1984 ) 1975 1985 ! 1976 1986 !-- k index has to be modified near bottom and top, else array … … 1984 1994 k_mm = k - 2 * ibit17 1985 1995 1986 w_comp 1987 flux_t(k) = w_comp * rho_air_zw(k) * (&1996 w_comp(k) = w(k,j,i) + w(k,j,i-1) 1997 flux_t(k) = w_comp(k) * rho_air_zw(k) * ( & 1988 1998 ( 37.0_wp * ibit17 * adv_mom_5 & 1989 1999 + 7.0_wp * ibit16 * adv_mom_3 & … … 1998 2008 ) * & 1999 2009 ( u(k_ppp,j,i) + u(k_mm,j,i) ) & 2000 )2001 2002 diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (&2010 ) 2011 2012 diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( & 2003 2013 ( 10.0_wp * ibit17 * adv_mom_5 & 2004 2014 + 3.0_wp * ibit16 * adv_mom_3 & … … 2013 2023 ) * & 2014 2024 ( u(k_ppp,j,i) - u(k_mm,j,i) ) & 2015 ) 2025 ) 2026 ENDDO 2027 2028 DO k = nzb+1, nzb_max 2029 2030 flux_d = flux_t(k-1) 2031 diss_d = diss_t(k-1) 2016 2032 ! 2017 2033 !-- Calculate the divergence of the velocity field. A respective … … 2025 2041 ) & 2026 2042 ) * ddx & 2027 + ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 )&2043 + ( ( v_comp(k) + gv ) * ( ibit12 + ibit13 + ibit14 ) & 2028 2044 - ( v(k,j,i) + v(k,j,i-1 ) ) & 2029 2045 * ( IBITS(advc_flags_1(k,j-1,i),12,1) & … … 2032 2048 ) & 2033 2049 ) * ddy & 2034 + ( w_comp * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 )&2035 - ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1)&2050 + ( w_comp(k) * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 )& 2051 - w_comp(k-1) * rho_air_zw(k-1) & 2036 2052 * ( IBITS(advc_flags_1(k-1,j,i),15,1) & 2037 2053 + IBITS(advc_flags_1(k-1,j,i),16,1) & … … 2056 2072 flux_s_u(k,tn) = flux_n(k) 2057 2073 diss_s_u(k,tn) = diss_n(k) 2058 flux_d = flux_t(k)2059 diss_d = diss_t(k)2060 2074 ! 2061 2075 !-- Statistical Evaluation of u'u'. The factor has to be applied for … … 2073 2087 sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) & 2074 2088 + ( flux_t(k) & 2075 * ( w_comp - 2.0_wp * hom(k,1,3,0)) &2076 / ( w_comp + SIGN( 1.0E-20_wp, w_comp )) &2089 * ( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) & 2090 / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) ) ) & 2077 2091 + diss_t(k) & 2078 * ABS( w_comp - 2.0_wp * hom(k,1,3,0)) &2079 / ( ABS( w_comp ) + 1.0E-20_wp) &2092 * ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) & 2093 / ( ABS( w_comp(k) ) + 1.0E-20_wp ) & 2080 2094 ) * weight_substep(intermediate_timestep_count) 2081 2095 ENDDO … … 2093 2107 + ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5 2094 2108 2095 v_comp 2096 flux_n(k) = v_comp * (&2109 v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv 2110 flux_n(k) = v_comp(k) * ( & 2097 2111 37.0_wp * ( u(k,j+1,i) + u(k,j,i) ) & 2098 2112 - 8.0_wp * ( u(k,j+2,i) + u(k,j-1,i) ) & 2099 2113 + ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5 2100 diss_n(k) = - ABS( v_comp ) * (&2114 diss_n(k) = - ABS( v_comp(k) ) * ( & 2101 2115 10.0_wp * ( u(k,j+1,i) - u(k,j,i) ) & 2102 2116 - 5.0_wp * ( u(k,j+2,i) - u(k,j-1,i) ) & … … 2113 2127 k_mm = k - 2 * ibit17 2114 2128 2115 w_comp 2116 flux_t(k) = w_comp * rho_air_zw(k) * (&2129 w_comp(k) = w(k,j,i) + w(k,j,i-1) 2130 flux_t(k) = w_comp(k) * rho_air_zw(k) * ( & 2117 2131 ( 37.0_wp * ibit17 * adv_mom_5 & 2118 2132 + 7.0_wp * ibit16 * adv_mom_3 & … … 2127 2141 ) * & 2128 2142 ( u(k_ppp,j,i) + u(k_mm,j,i) ) & 2129 )2130 2131 diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (&2143 ) 2144 2145 diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( & 2132 2146 ( 10.0_wp * ibit17 * adv_mom_5 & 2133 2147 + 3.0_wp * ibit16 * adv_mom_3 & … … 2142 2156 ) * & 2143 2157 ( u(k_ppp,j,i) - u(k_mm,j,i) ) & 2144 ) 2158 ) 2159 2160 ENDDO 2161 2162 DO k = nzb_max+1, nzt 2163 2164 flux_d = flux_t(k-1) 2165 diss_d = diss_t(k-1) 2145 2166 ! 2146 2167 !-- Calculate the divergence of the velocity field. A respective 2147 2168 !-- correction is needed to overcome numerical instabilities introduced 2148 2169 !-- by a not sufficient reduction of divergences near topography. 2149 div = ( ( u_comp(k) - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx&2150 + ( v_comp + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy&2151 + ( w_comp * rho_air_zw(k) -&2152 ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1)&2170 div = ( ( u_comp(k) - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx & 2171 + ( v_comp(k) + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy & 2172 + ( w_comp(k) * rho_air_zw(k) - & 2173 w_comp(k-1) * rho_air_zw(k-1) & 2153 2174 ) * drho_air(k) * ddzw(k) & 2154 2175 ) * 0.5_wp … … 2168 2189 flux_s_u(k,tn) = flux_n(k) 2169 2190 diss_s_u(k,tn) = diss_n(k) 2170 flux_d = flux_t(k)2171 diss_d = diss_t(k)2172 2191 ! 2173 2192 !-- Statistical Evaluation of u'u'. The factor has to be applied for … … 2185 2204 sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) & 2186 2205 + ( flux_t(k) & 2187 * ( w_comp - 2.0_wp * hom(k,1,3,0) )&2188 / ( w_comp + SIGN( 1.0E-20_wp, w_comp ) )&2206 * ( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) & 2207 / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) ) ) & 2189 2208 + diss_t(k) & 2190 * ABS( w_comp - 2.0_wp * hom(k,1,3,0) )&2191 / ( ABS( w_comp ) + 1.0E-20_wp )&2209 * ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) & 2210 / ( ABS( w_comp(k) ) + 1.0E-20_wp ) & 2192 2211 ) * weight_substep(intermediate_timestep_count) 2193 2212 ENDDO … … 2254 2273 REAL(wp) :: gu !< 2255 2274 REAL(wp) :: gv !< 2256 REAL(wp) :: u_comp !<2257 2275 REAL(wp) :: v_comp_l !< 2258 REAL(wp) :: w_comp !<2259 2276 2260 2277 REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !< … … 2264 2281 REAL(wp), DIMENSION(nzb:nzt+1) :: flux_r !< 2265 2282 REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !< 2283 REAL(wp), DIMENSION(nzb:nzt+1) :: u_comp !< 2266 2284 REAL(wp), DIMENSION(nzb:nzt+1) :: v_comp !< 2285 REAL(wp), DIMENSION(nzb:nzt+1) :: w_comp !< 2267 2286 2268 2287 gu = 2.0_wp * u_gtrans … … 2279 2298 ibit18 = IBITS(advc_flags_1(k,j,i-1),18,1) 2280 2299 2281 u_comp 2282 flux_l_v(k,j,tn) = u_comp * (&2300 u_comp(k) = u(k,j-1,i) + u(k,j,i) - gu 2301 flux_l_v(k,j,tn) = u_comp(k) * ( & 2283 2302 ( 37.0_wp * ibit20 * adv_mom_5 & 2284 2303 + 7.0_wp * ibit19 * adv_mom_3 & … … 2293 2312 ) * & 2294 2313 ( v(k,j,i+2) + v(k,j,i-3) ) & 2295 )2296 2297 diss_l_v(k,j,tn) = - ABS( u_comp ) * (&2314 ) 2315 2316 diss_l_v(k,j,tn) = - ABS( u_comp(k) ) * ( & 2298 2317 ( 10.0_wp * ibit20 * adv_mom_5 & 2299 2318 + 3.0_wp * ibit19 * adv_mom_3 & … … 2308 2327 ) * & 2309 2328 ( v(k,j,i+2) - v(k,j,i-3) ) & 2310 )2329 ) 2311 2330 2312 2331 ENDDO … … 2314 2333 DO k = nzb_max+1, nzt 2315 2334 2316 u_comp = u(k,j-1,i) + u(k,j,i) - gu2317 flux_l_v(k,j,tn) = u_comp * (&2335 u_comp(k) = u(k,j-1,i) + u(k,j,i) - gu 2336 flux_l_v(k,j,tn) = u_comp(k) * ( & 2318 2337 37.0_wp * ( v(k,j,i) + v(k,j,i-1) ) & 2319 2338 - 8.0_wp * ( v(k,j,i+1) + v(k,j,i-2) ) & 2320 2339 + ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5 2321 diss_l_v(k,j,tn) = - ABS( u_comp ) * (&2340 diss_l_v(k,j,tn) = - ABS( u_comp(k) ) * ( & 2322 2341 10.0_wp * ( v(k,j,i) - v(k,j,i-1) ) & 2323 2342 - 5.0_wp * ( v(k,j,i+1) - v(k,j,i-2) ) & … … 2388 2407 flux_t(0) = 0.0_wp 2389 2408 diss_t(0) = 0.0_wp 2390 flux_d = 0.0_wp 2391 diss_d = 0.0_wp 2409 w_comp(0) = 0.0_wp 2392 2410 ! 2393 2411 !-- Now compute the fluxes and tendency terms for the horizontal and … … 2399 2417 ibit18 = IBITS(advc_flags_1(k,j,i),18,1) 2400 2418 2401 u_comp 2402 flux_r(k) = u_comp * (&2419 u_comp(k) = u(k,j-1,i+1) + u(k,j,i+1) - gu 2420 flux_r(k) = u_comp(k) * ( & 2403 2421 ( 37.0_wp * ibit20 * adv_mom_5 & 2404 2422 + 7.0_wp * ibit19 * adv_mom_3 & … … 2413 2431 ) * & 2414 2432 ( v(k,j,i+3) + v(k,j,i-2) ) & 2415 )2416 2417 diss_r(k) = - ABS( u_comp ) * (&2433 ) 2434 2435 diss_r(k) = - ABS( u_comp(k) ) * ( & 2418 2436 ( 10.0_wp * ibit20 * adv_mom_5 & 2419 2437 + 3.0_wp * ibit19 * adv_mom_3 & … … 2428 2446 ) * & 2429 2447 ( v(k,j,i+3) - v(k,j,i-2) ) & 2430 )2448 ) 2431 2449 2432 2450 ibit23 = IBITS(advc_flags_1(k,j,i),23,1) … … 2476 2494 k_mm = k - 2 * ibit26 2477 2495 2478 w_comp 2479 flux_t(k) = w_comp * rho_air_zw(k) * (&2496 w_comp(k) = w(k,j-1,i) + w(k,j,i) 2497 flux_t(k) = w_comp(k) * rho_air_zw(k) * ( & 2480 2498 ( 37.0_wp * ibit26 * adv_mom_5 & 2481 2499 + 7.0_wp * ibit25 * adv_mom_3 & … … 2490 2508 ) * & 2491 2509 ( v(k_ppp,j,i) + v(k_mm,j,i) ) & 2492 )2493 2494 diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (&2510 ) 2511 2512 diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( & 2495 2513 ( 10.0_wp * ibit26 * adv_mom_5 & 2496 2514 + 3.0_wp * ibit25 * adv_mom_3 & … … 2505 2523 ) * & 2506 2524 ( v(k_ppp,j,i) - v(k_mm,j,i) ) & 2507 ) 2525 ) 2526 ENDDO 2527 2528 DO k = nzb+1, nzb_max 2529 2530 flux_d = flux_t(k-1) 2531 diss_d = diss_t(k-1) 2508 2532 ! 2509 2533 !-- Calculate the divergence of the velocity field. A respective 2510 2534 !-- correction is needed to overcome numerical instabilities introduced 2511 2535 !-- by a not sufficient reduction of divergences near topography. 2512 div = ( ( ( u_comp + gu )&2536 div = ( ( ( u_comp(k) + gu ) & 2513 2537 * ( ibit18 + ibit19 + ibit20 ) & 2514 2538 - ( u(k,j-1,i) + u(k,j,i) ) & … … 2526 2550 ) & 2527 2551 ) * ddy & 2528 + ( w_comp * rho_air_zw(k) * ( ibit24 + ibit25 + ibit26 )&2529 - ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1)&2552 + ( w_comp(k) * rho_air_zw(k) * ( ibit24 + ibit25 + ibit26 )& 2553 - w_comp(k-1) * rho_air_zw(k-1) & 2530 2554 * ( IBITS(advc_flags_1(k-1,j,i),24,1) & 2531 2555 + IBITS(advc_flags_1(k-1,j,i),25,1) & … … 2550 2574 flux_s_v(k,tn) = flux_n(k) 2551 2575 diss_s_v(k,tn) = diss_n(k) 2552 flux_d = flux_t(k)2553 diss_d = diss_t(k)2554 2576 2555 2577 ! … … 2568 2590 sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) & 2569 2591 + ( flux_t(k) & 2570 * ( w_comp - 2.0_wp * hom(k,1,3,0) )&2571 / ( w_comp + SIGN( 1.0E-20_wp, w_comp ) )&2592 * ( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) & 2593 / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) ) ) & 2572 2594 + diss_t(k) & 2573 * ABS( w_comp - 2.0_wp * hom(k,1,3,0) )&2574 / ( ABS( w_comp ) + 1.0E-20_wp )&2595 * ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) & 2596 / ( ABS( w_comp(k) ) + 1.0E-20_wp ) & 2575 2597 ) * weight_substep(intermediate_timestep_count) 2576 2598 … … 2579 2601 DO k = nzb_max+1, nzt 2580 2602 2581 u_comp 2582 flux_r(k) = u_comp * (&2603 u_comp(k) = u(k,j-1,i+1) + u(k,j,i+1) - gu 2604 flux_r(k) = u_comp(k) * ( & 2583 2605 37.0_wp * ( v(k,j,i+1) + v(k,j,i) ) & 2584 2606 - 8.0_wp * ( v(k,j,i+2) + v(k,j,i-1) ) & 2585 2607 + ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5 2586 2608 2587 diss_r(k) = - ABS( u_comp ) * (&2609 diss_r(k) = - ABS( u_comp(k) ) * ( & 2588 2610 10.0_wp * ( v(k,j,i+1) - v(k,j,i) ) & 2589 2611 - 5.0_wp * ( v(k,j,i+2) - v(k,j,i-1) ) & … … 2612 2634 k_mm = k - 2 * ibit26 2613 2635 2614 w_comp 2615 flux_t(k) = w_comp * rho_air_zw(k) * (&2636 w_comp(k) = w(k,j-1,i) + w(k,j,i) 2637 flux_t(k) = w_comp(k) * rho_air_zw(k) * ( & 2616 2638 ( 37.0_wp * ibit26 * adv_mom_5 & 2617 2639 + 7.0_wp * ibit25 * adv_mom_3 & … … 2626 2648 ) * & 2627 2649 ( v(k_ppp,j,i) + v(k_mm,j,i) ) & 2628 )2629 2630 diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (&2650 ) 2651 2652 diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( & 2631 2653 ( 10.0_wp * ibit26 * adv_mom_5 & 2632 2654 + 3.0_wp * ibit25 * adv_mom_3 & … … 2641 2663 ) * & 2642 2664 ( v(k_ppp,j,i) - v(k_mm,j,i) ) & 2643 ) 2665 ) 2666 ENDDO 2667 2668 DO k = nzb_max+1, nzt 2669 2670 flux_d = flux_t(k-1) 2671 diss_d = diss_t(k-1) 2644 2672 ! 2645 2673 !-- Calculate the divergence of the velocity field. A respective 2646 2674 !-- correction is needed to overcome numerical instabilities introduced 2647 2675 !-- by a not sufficient reduction of divergences near topography. 2648 div = ( ( u_comp + gu - ( u(k,j-1,i) + u(k,j,i) ) ) * ddx&2649 + ( v_comp(k) - ( v(k,j,i) + v(k,j-1,i) ) ) * ddy&2650 + ( w_comp * rho_air_zw(k) -&2651 ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1)&2676 div = ( ( u_comp(k) + gu - ( u(k,j-1,i) + u(k,j,i) ) ) * ddx & 2677 + ( v_comp(k) - ( v(k,j,i) + v(k,j-1,i) ) ) * ddy & 2678 + ( w_comp(k) * rho_air_zw(k) - & 2679 w_comp(k-1) * rho_air_zw(k-1) & 2652 2680 ) * drho_air(k) * ddzw(k) & 2653 2681 ) * 0.5_wp … … 2667 2695 flux_s_v(k,tn) = flux_n(k) 2668 2696 diss_s_v(k,tn) = diss_n(k) 2669 flux_d = flux_t(k)2670 diss_d = diss_t(k)2671 2672 2697 ! 2673 2698 !-- Statistical Evaluation of v'v'. The factor has to be applied for … … 2685 2710 sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) & 2686 2711 + ( flux_t(k) & 2687 * ( w_comp - 2.0_wp * hom(k,1,3,0) )&2688 / ( w_comp + SIGN( 1.0E-20_wp, w_comp ) )&2712 * ( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) & 2713 / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) ) ) & 2689 2714 + diss_t(k) & 2690 * ABS( w_comp - 2.0_wp * hom(k,1,3,0) )&2691 / ( ABS( w_comp ) + 1.0E-20_wp )&2715 * ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) & 2716 / ( ABS( w_comp(k) ) + 1.0E-20_wp ) & 2692 2717 ) * weight_substep(intermediate_timestep_count) 2693 2718 … … 2754 2779 REAL(wp) :: gu !< 2755 2780 REAL(wp) :: gv !< 2756 REAL(wp) :: u_comp !<2757 REAL(wp) :: v_comp !<2758 REAL(wp) :: w_comp !<2759 2781 2760 2782 REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !< … … 2764 2786 REAL(wp), DIMENSION(nzb:nzt+1) :: flux_r !< 2765 2787 REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !< 2788 REAL(wp), DIMENSION(nzb:nzt+1) :: u_comp !< 2789 REAL(wp), DIMENSION(nzb:nzt+1) :: v_comp !< 2790 REAL(wp), DIMENSION(nzb:nzt+1) :: w_comp !< 2766 2791 2767 2792 gu = 2.0_wp * u_gtrans … … 2777 2802 ibit30 = IBITS(advc_flags_1(k,j-1,i),30,1) 2778 2803 2779 v_comp 2780 flux_s_w(k,tn) = v_comp * (&2804 v_comp(k) = v(k+1,j,i) + v(k,j,i) - gv 2805 flux_s_w(k,tn) = v_comp(k) * ( & 2781 2806 ( 37.0_wp * ibit32 * adv_mom_5 & 2782 2807 + 7.0_wp * ibit31 * adv_mom_3 & … … 2791 2816 ) * & 2792 2817 ( w(k,j+2,i) + w(k,j-3,i) ) & 2793 )2794 2795 diss_s_w(k,tn) = - ABS( v_comp ) * (&2818 ) 2819 2820 diss_s_w(k,tn) = - ABS( v_comp(k) ) * ( & 2796 2821 ( 10.0_wp * ibit32 * adv_mom_5 & 2797 2822 + 3.0_wp * ibit31 * adv_mom_3 & … … 2806 2831 ) * & 2807 2832 ( w(k,j+2,i) - w(k,j-3,i) ) & 2808 )2833 ) 2809 2834 2810 2835 ENDDO … … 2812 2837 DO k = nzb_max+1, nzt 2813 2838 2814 v_comp 2815 flux_s_w(k,tn) = v_comp * (&2839 v_comp(k) = v(k+1,j,i) + v(k,j,i) - gv 2840 flux_s_w(k,tn) = v_comp(k) * ( & 2816 2841 37.0_wp * ( w(k,j,i) + w(k,j-1,i) ) & 2817 2842 - 8.0_wp * ( w(k,j+1,i) +w(k,j-2,i) ) & 2818 2843 + ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5 2819 diss_s_w(k,tn) = - ABS( v_comp ) * (&2844 diss_s_w(k,tn) = - ABS( v_comp(k) ) * ( & 2820 2845 10.0_wp * ( w(k,j,i) - w(k,j-1,i) ) & 2821 2846 - 5.0_wp * ( w(k,j+1,i) - w(k,j-2,i) ) & … … 2835 2860 ibit27 = IBITS(advc_flags_1(k,j,i-1),27,1) 2836 2861 2837 u_comp 2838 flux_l_w(k,j,tn) = u_comp * (&2862 u_comp(k) = u(k+1,j,i) + u(k,j,i) - gu 2863 flux_l_w(k,j,tn) = u_comp(k) * ( & 2839 2864 ( 37.0_wp * ibit29 * adv_mom_5 & 2840 2865 + 7.0_wp * ibit28 * adv_mom_3 & … … 2849 2874 ) * & 2850 2875 ( w(k,j,i+2) + w(k,j,i-3) ) & 2851 )2852 2853 diss_l_w(k,j,tn) = - ABS( u_comp ) * (&2876 ) 2877 2878 diss_l_w(k,j,tn) = - ABS( u_comp(k) ) * ( & 2854 2879 ( 10.0_wp * ibit29 * adv_mom_5 & 2855 2880 + 3.0_wp * ibit28 * adv_mom_3 & … … 2864 2889 ) * & 2865 2890 ( w(k,j,i+2) - w(k,j,i-3) ) & 2866 )2891 ) 2867 2892 2868 2893 ENDDO … … 2870 2895 DO k = nzb_max+1, nzt 2871 2896 2872 u_comp 2873 flux_l_w(k,j,tn) = u_comp * (&2897 u_comp(k) = u(k+1,j,i) + u(k,j,i) - gu 2898 flux_l_w(k,j,tn) = u_comp(k) * ( & 2874 2899 37.0_wp * ( w(k,j,i) + w(k,j,i-1) ) & 2875 2900 - 8.0_wp * ( w(k,j,i+1) + w(k,j,i-2) ) & 2876 2901 + ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5 2877 diss_l_w(k,j,tn) = - ABS( u_comp ) * (&2902 diss_l_w(k,j,tn) = - ABS( u_comp(k) ) * ( & 2878 2903 10.0_wp * ( w(k,j,i) - w(k,j,i-1) ) & 2879 2904 - 5.0_wp * ( w(k,j,i+1) - w(k,j,i-2) ) & … … 2888 2913 !-- advc_flags_1. 2889 2914 k = nzb + 1 2890 w_comp = w(k,j,i) + w(k-1,j,i) 2891 flux_t(0) = w_comp * ( w(k,j,i) + w(k-1,j,i) ) * adv_mom_1 2892 diss_t(0) = -ABS(w_comp) * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1 2893 flux_d = flux_t(0) 2894 diss_d = diss_t(0) 2915 w_comp(k) = w(k,j,i) + w(k-1,j,i) 2916 flux_t(0) = w_comp(k) * ( w(k,j,i) + w(k-1,j,i) ) * adv_mom_1 2917 diss_t(0) = -ABS(w_comp(k)) * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1 2895 2918 ! 2896 2919 !-- Now compute the fluxes and tendency terms for the horizontal … … 2902 2925 ibit27 = IBITS(advc_flags_1(k,j,i),27,1) 2903 2926 2904 u_comp 2905 flux_r(k) = u_comp * (&2927 u_comp(k) = u(k+1,j,i+1) + u(k,j,i+1) - gu 2928 flux_r(k) = u_comp(k) * ( & 2906 2929 ( 37.0_wp * ibit29 * adv_mom_5 & 2907 2930 + 7.0_wp * ibit28 * adv_mom_3 & … … 2916 2939 ) * & 2917 2940 ( w(k,j,i+3) + w(k,j,i-2) ) & 2918 )2919 2920 diss_r(k) = - ABS( u_comp ) * (&2941 ) 2942 2943 diss_r(k) = - ABS( u_comp(k) ) * ( & 2921 2944 ( 10.0_wp * ibit29 * adv_mom_5 & 2922 2945 + 3.0_wp * ibit28 * adv_mom_3 & … … 2931 2954 ) * & 2932 2955 ( w(k,j,i+3) - w(k,j,i-2) ) & 2933 )2956 ) 2934 2957 2935 2958 ibit32 = IBITS(advc_flags_2(k,j,i),0,1) … … 2937 2960 ibit30 = IBITS(advc_flags_1(k,j,i),30,1) 2938 2961 2939 v_comp 2940 flux_n(k) = v_comp * (&2962 v_comp(k) = v(k+1,j+1,i) + v(k,j+1,i) - gv 2963 flux_n(k) = v_comp(k) * ( & 2941 2964 ( 37.0_wp * ibit32 * adv_mom_5 & 2942 2965 + 7.0_wp * ibit31 * adv_mom_3 & … … 2951 2974 ) * & 2952 2975 ( w(k,j+3,i) + w(k,j-2,i) ) & 2953 )2954 2955 diss_n(k) = - ABS( v_comp ) * (&2976 ) 2977 2978 diss_n(k) = - ABS( v_comp(k) ) * ( & 2956 2979 ( 10.0_wp * ibit32 * adv_mom_5 & 2957 2980 + 3.0_wp * ibit31 * adv_mom_3 & … … 2966 2989 ) * & 2967 2990 ( w(k,j+3,i) - w(k,j-2,i) ) & 2968 )2991 ) 2969 2992 ! 2970 2993 !-- k index has to be modified near bottom and top, else array … … 2978 3001 k_mm = k - 2 * ibit35 2979 3002 2980 w_comp 2981 flux_t(k) = w_comp * rho_air(k+1) * (&3003 w_comp(k) = w(k+1,j,i) + w(k,j,i) 3004 flux_t(k) = w_comp(k) * rho_air(k+1) * ( & 2982 3005 ( 37.0_wp * ibit35 * adv_mom_5 & 2983 3006 + 7.0_wp * ibit34 * adv_mom_3 & … … 2992 3015 ) * & 2993 3016 ( w(k_ppp,j,i) + w(k_mm,j,i) ) & 2994 )2995 2996 diss_t(k) = - ABS( w_comp ) * rho_air(k+1) * (&3017 ) 3018 3019 diss_t(k) = - ABS( w_comp(k) ) * rho_air(k+1) * ( & 2997 3020 ( 10.0_wp * ibit35 * adv_mom_5 & 2998 3021 + 3.0_wp * ibit34 * adv_mom_3 & … … 3007 3030 ) * & 3008 3031 ( w(k_ppp,j,i) - w(k_mm,j,i) ) & 3009 ) 3010 3032 ) 3033 ENDDO 3034 3035 DO k = nzb+1, nzb_max 3036 3037 flux_d = flux_t(k-1) 3038 diss_d = diss_t(k-1) 3011 3039 ! 3012 3040 !-- Calculate the divergence of the velocity field. A respective 3013 3041 !-- correction is needed to overcome numerical instabilities introduced 3014 3042 !-- by a not sufficient reduction of divergences near topography. 3015 div = ( ( ( u_comp + gu ) * ( ibit27 + ibit28 + ibit29 )&3043 div = ( ( ( u_comp(k) + gu ) * ( ibit27 + ibit28 + ibit29 ) & 3016 3044 - ( u(k+1,j,i) + u(k,j,i) ) & 3017 3045 * ( IBITS(advc_flags_1(k,j,i-1),27,1) & … … 3020 3048 ) & 3021 3049 ) * ddx & 3022 + ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 )&3050 + ( ( v_comp(k) + gv ) * ( ibit30 + ibit31 + ibit32 ) & 3023 3051 - ( v(k+1,j,i) + v(k,j,i) ) & 3024 3052 * ( IBITS(advc_flags_1(k,j-1,i),30,1) & … … 3027 3055 ) & 3028 3056 ) * ddy & 3029 + ( w_comp * rho_air(k+1) * ( ibit33 + ibit34 + ibit35 ) & 3030 - ( w(k,j,i) + w(k-1,j,i) ) * rho_air(k) & 3057 + ( w_comp(k) * rho_air(k+1) & 3058 * ( ibit33 + ibit34 + ibit35 ) & 3059 - ( w(k,j,i) + w(k-1,j,i) ) * rho_air(k) & 3031 3060 * ( IBITS(advc_flags_2(k-1,j,i),1,1) & 3032 3061 + IBITS(advc_flags_2(k-1,j,i),2,1) & … … 3051 3080 flux_s_w(k,tn) = flux_n(k) 3052 3081 diss_s_w(k,tn) = diss_n(k) 3053 flux_d = flux_t(k)3054 diss_d = diss_t(k)3055 3082 ! 3056 3083 !-- Statistical Evaluation of w'w'. 3057 3084 sums_ws2_ws_l(k,tn) = sums_ws2_ws_l(k,tn) & 3058 3085 + ( flux_t(k) & 3059 * ( w_comp - 2.0_wp * hom(k,1,3,0)) &3060 / ( w_comp + SIGN( 1.0E-20_wp, w_comp )) &3086 * ( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) & 3087 / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) ) ) & 3061 3088 + diss_t(k) & 3062 * ABS( w_comp - 2.0_wp * hom(k,1,3,0)) &3063 / ( ABS( w_comp ) + 1.0E-20_wp) &3089 * ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) & 3090 / ( ABS( w_comp(k) ) + 1.0E-20_wp ) & 3064 3091 ) * weight_substep(intermediate_timestep_count) 3065 3092 … … 3068 3095 DO k = nzb_max+1, nzt 3069 3096 3070 u_comp 3071 flux_r(k) = u_comp * (&3097 u_comp(k) = u(k+1,j,i+1) + u(k,j,i+1) - gu 3098 flux_r(k) = u_comp(k) * ( & 3072 3099 37.0_wp * ( w(k,j,i+1) + w(k,j,i) ) & 3073 3100 - 8.0_wp * ( w(k,j,i+2) + w(k,j,i-1) ) & 3074 3101 + ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5 3075 3102 3076 diss_r(k) = - ABS( u_comp ) * (&3103 diss_r(k) = - ABS( u_comp(k) ) * ( & 3077 3104 10.0_wp * ( w(k,j,i+1) - w(k,j,i) ) & 3078 3105 - 5.0_wp * ( w(k,j,i+2) - w(k,j,i-1) ) & 3079 3106 + ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5 3080 3107 3081 v_comp 3082 flux_n(k) = v_comp * (&3108 v_comp(k) = v(k+1,j+1,i) + v(k,j+1,i) - gv 3109 flux_n(k) = v_comp(k) * ( & 3083 3110 37.0_wp * ( w(k,j+1,i) + w(k,j,i) ) & 3084 3111 - 8.0_wp * ( w(k,j+2,i) + w(k,j-1,i) ) & 3085 3112 + ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5 3086 3113 3087 diss_n(k) = - ABS( v_comp ) * (&3114 diss_n(k) = - ABS( v_comp(k) ) * ( & 3088 3115 10.0_wp * ( w(k,j+1,i) - w(k,j,i) ) & 3089 3116 - 5.0_wp * ( w(k,j+2,i) - w(k,j-1,i) ) & … … 3100 3127 k_mm = k - 2 * ibit35 3101 3128 3102 w_comp 3103 flux_t(k) = w_comp * rho_air(k+1) * (&3129 w_comp(k) = w(k+1,j,i) + w(k,j,i) 3130 flux_t(k) = w_comp(k) * rho_air(k+1) * ( & 3104 3131 ( 37.0_wp * ibit35 * adv_mom_5 & 3105 3132 + 7.0_wp * ibit34 * adv_mom_3 & … … 3114 3141 ) * & 3115 3142 ( w(k_ppp,j,i) + w(k_mm,j,i) ) & 3116 )3117 3118 diss_t(k) = - ABS( w_comp ) * rho_air(k+1) * (&3143 ) 3144 3145 diss_t(k) = - ABS( w_comp(k) ) * rho_air(k+1) * ( & 3119 3146 ( 10.0_wp * ibit35 * adv_mom_5 & 3120 3147 + 3.0_wp * ibit34 * adv_mom_3 & … … 3129 3156 ) * & 3130 3157 ( w(k_ppp,j,i) - w(k_mm,j,i) ) & 3131 ) 3158 ) 3159 ENDDO 3160 3161 DO k = nzb_max+1, nzt 3162 3163 flux_d = flux_t(k-1) 3164 diss_d = diss_t(k-1) 3132 3165 ! 3133 3166 !-- Calculate the divergence of the velocity field. A respective 3134 3167 !-- correction is needed to overcome numerical instabilities introduced 3135 3168 !-- by a not sufficient reduction of divergences near topography. 3136 div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i) ) ) * ddx&3137 + ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i) ) ) * ddy&3138 + ( w_comp * rho_air(k+1) -&3139 ( w(k,j,i) + w(k-1,j,i) ) * rho_air(k)&3169 div = ( ( u_comp(k) + gu - ( u(k+1,j,i) + u(k,j,i) ) ) * ddx & 3170 + ( v_comp(k) + gv - ( v(k+1,j,i) + v(k,j,i) ) ) * ddy & 3171 + ( w_comp(k) * rho_air(k+1) - & 3172 ( w(k,j,i) + w(k-1,j,i) ) * rho_air(k) & 3140 3173 ) * drho_air_zw(k) * ddzu(k+1) & 3141 3174 ) * 0.5_wp … … 3155 3188 flux_s_w(k,tn) = flux_n(k) 3156 3189 diss_s_w(k,tn) = diss_n(k) 3157 flux_d = flux_t(k)3158 diss_d = diss_t(k)3159 3190 ! 3160 3191 !-- Statistical Evaluation of w'w'. 3161 3192 sums_ws2_ws_l(k,tn) = sums_ws2_ws_l(k,tn) & 3162 3193 + ( flux_t(k) & 3163 * ( w_comp - 2.0_wp * hom(k,1,3,0)) &3164 / ( w_comp + SIGN( 1.0E-20_wp, w_comp )) &3194 * ( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) & 3195 / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) ) ) & 3165 3196 + diss_t(k) & 3166 * ABS( w_comp - 2.0_wp * hom(k,1,3,0)) &3167 / ( ABS( w_comp ) + 1.0E-20_wp) &3197 * ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) & 3198 / ( ABS( w_comp(k) ) + 1.0E-20_wp ) & 3168 3199 ) * weight_substep(intermediate_timestep_count) 3169 3200
Note: See TracChangeset
for help on using the changeset viewer.
