Changeset 1015 for palm/trunk/SOURCE/prognostic_equations.f90
- Timestamp:
- Sep 27, 2012 9:23:24 AM (12 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/prognostic_equations.f90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/prognostic_equations.f90
r1002 r1015 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! new branch prognostic_equations_acc 7 ! OpenACC statements added + code changes required for GPU optimization 7 8 ! 8 9 ! Former revisions: … … 136 137 PRIVATE 137 138 PUBLIC prognostic_equations_noopt, prognostic_equations_cache, & 138 prognostic_equations_vector 139 prognostic_equations_vector, prognostic_equations_acc 139 140 140 141 INTERFACE prognostic_equations_noopt … … 149 150 MODULE PROCEDURE prognostic_equations_vector 150 151 END INTERFACE prognostic_equations_vector 152 153 INTERFACE prognostic_equations_acc 154 MODULE PROCEDURE prognostic_equations_acc 155 END INTERFACE prognostic_equations_acc 151 156 152 157 … … 1900 1905 1901 1906 1907 SUBROUTINE prognostic_equations_acc 1908 1909 !------------------------------------------------------------------------------! 1910 ! Version for accelerator boards 1911 !------------------------------------------------------------------------------! 1912 1913 IMPLICIT NONE 1914 1915 CHARACTER (LEN=9) :: time_to_string 1916 INTEGER :: i, j, k, runge_step 1917 REAL :: sbt 1918 1919 ! 1920 !-- Set switch for intermediate Runge-Kutta step 1921 runge_step = 0 1922 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1923 IF ( intermediate_timestep_count == 1 ) THEN 1924 runge_step = 1 1925 ELSEIF ( intermediate_timestep_count < & 1926 intermediate_timestep_count_max ) THEN 1927 runge_step = 2 1928 ENDIF 1929 ENDIF 1930 1931 ! 1932 !-- Calculate those variables needed in the tendency terms which need 1933 !-- global communication 1934 IF ( .NOT. neutral ) CALL calc_mean_profile( pt, 4 ) 1935 IF ( ocean ) CALL calc_mean_profile( rho, 64 ) 1936 IF ( humidity ) CALL calc_mean_profile( vpt, 44 ) 1937 IF ( ( ws_scheme_mom .OR. ws_scheme_sca ) .AND. & 1938 intermediate_timestep_count == 1 ) CALL ws_statistics 1939 1940 ! 1941 !-- u-velocity component 1942 !++ Statistics still not ported to accelerators 1943 !$acc update device( hom ) 1944 CALL cpu_log( log_point(5), 'u-equation', 'start' ) 1945 1946 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1947 IF ( ws_scheme_mom ) THEN 1948 CALL advec_u_ws_acc 1949 ELSE 1950 tend = 0.0 ! to be removed later?? 1951 CALL advec_u_pw 1952 ENDIF 1953 ELSE 1954 CALL advec_u_up 1955 ENDIF 1956 CALL diffusion_u_acc 1957 CALL coriolis_acc( 1 ) 1958 IF ( sloping_surface .AND. .NOT. neutral ) THEN 1959 CALL buoyancy( pt, pt_reference, 1, 4 ) 1960 ENDIF 1961 1962 ! 1963 !-- Drag by plant canopy 1964 IF ( plant_canopy ) CALL plant_canopy_model( 1 ) 1965 1966 ! 1967 !-- External pressure gradient 1968 IF ( dp_external ) THEN 1969 DO i = nxlu, nxr 1970 DO j = nys, nyn 1971 DO k = dp_level_ind_b+1, nzt 1972 tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k) 1973 ENDDO 1974 ENDDO 1975 ENDDO 1976 ENDIF 1977 1978 CALL user_actions( 'u-tendency' ) 1979 1980 ! 1981 !-- Prognostic equation for u-velocity component 1982 !$acc kernels present( nzb_u_inner, rdf, tend, tu_m, u, ug, u_p ) 1983 !$acc loop 1984 DO i = nxlu, nxr 1985 DO j = nys, nyn 1986 !$acc loop vector( 32 ) 1987 DO k = 1, nzt 1988 IF ( k > nzb_u_inner(j,i) ) THEN 1989 u_p(k,j,i) = u(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + & 1990 tsc(3) * tu_m(k,j,i) ) & 1991 - tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) ) 1992 ! 1993 !-- Tendencies for the next Runge-Kutta step 1994 IF ( runge_step == 1 ) THEN 1995 tu_m(k,j,i) = tend(k,j,i) 1996 ELSEIF ( runge_step == 2 ) THEN 1997 tu_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tu_m(k,j,i) 1998 ENDIF 1999 ENDIF 2000 ENDDO 2001 ENDDO 2002 ENDDO 2003 !$acc end kernels 2004 2005 CALL cpu_log( log_point(5), 'u-equation', 'stop' ) 2006 !$acc update host( u_p ) 2007 2008 ! 2009 !-- v-velocity component 2010 CALL cpu_log( log_point(6), 'v-equation', 'start' ) 2011 2012 IF ( timestep_scheme(1:5) == 'runge' ) THEN 2013 IF ( ws_scheme_mom ) THEN 2014 CALL advec_v_ws_acc 2015 ELSE 2016 tend = 0.0 ! to be removed later?? 2017 CALL advec_v_pw 2018 END IF 2019 ELSE 2020 CALL advec_v_up 2021 ENDIF 2022 CALL diffusion_v_acc 2023 CALL coriolis_acc( 2 ) 2024 2025 ! 2026 !-- Drag by plant canopy 2027 IF ( plant_canopy ) CALL plant_canopy_model( 2 ) 2028 2029 ! 2030 !-- External pressure gradient 2031 IF ( dp_external ) THEN 2032 DO i = nxl, nxr 2033 DO j = nysv, nyn 2034 DO k = dp_level_ind_b+1, nzt 2035 tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k) 2036 ENDDO 2037 ENDDO 2038 ENDDO 2039 ENDIF 2040 2041 CALL user_actions( 'v-tendency' ) 2042 2043 ! 2044 !-- Prognostic equation for v-velocity component 2045 !$acc kernels present( nzb_v_inner, rdf, tend, tv_m, v, vg, v_p ) 2046 !$acc loop 2047 DO i = nxl, nxr 2048 DO j = nysv, nyn 2049 !$acc loop vector( 32 ) 2050 DO k = 1, nzt 2051 IF ( k > nzb_v_inner(j,i) ) THEN 2052 v_p(k,j,i) = v(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + & 2053 tsc(3) * tv_m(k,j,i) ) & 2054 - tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) ) 2055 ! 2056 !-- Tendencies for the next Runge-Kutta step 2057 IF ( runge_step == 1 ) THEN 2058 tv_m(k,j,i) = tend(k,j,i) 2059 ELSEIF ( runge_step == 2 ) THEN 2060 tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i) 2061 ENDIF 2062 ENDIF 2063 ENDDO 2064 ENDDO 2065 ENDDO 2066 !$acc end kernels 2067 2068 CALL cpu_log( log_point(6), 'v-equation', 'stop' ) 2069 !$acc update host( v_p ) 2070 2071 ! 2072 !-- w-velocity component 2073 CALL cpu_log( log_point(7), 'w-equation', 'start' ) 2074 2075 IF ( timestep_scheme(1:5) == 'runge' ) THEN 2076 IF ( ws_scheme_mom ) THEN 2077 CALL advec_w_ws_acc 2078 ELSE 2079 tend = 0.0 ! to be removed later?? 2080 CALL advec_w_pw 2081 ENDIF 2082 ELSE 2083 CALL advec_w_up 2084 ENDIF 2085 CALL diffusion_w_acc 2086 CALL coriolis_acc( 3 ) 2087 2088 IF ( .NOT. neutral ) THEN 2089 IF ( ocean ) THEN 2090 CALL buoyancy( rho, rho_reference, 3, 64 ) 2091 ELSE 2092 IF ( .NOT. humidity ) THEN 2093 CALL buoyancy_acc( pt, pt_reference, 3, 4 ) 2094 ELSE 2095 CALL buoyancy( vpt, pt_reference, 3, 44 ) 2096 ENDIF 2097 ENDIF 2098 ENDIF 2099 2100 ! 2101 !-- Drag by plant canopy 2102 IF ( plant_canopy ) CALL plant_canopy_model( 3 ) 2103 2104 CALL user_actions( 'w-tendency' ) 2105 2106 ! 2107 !-- Prognostic equation for w-velocity component 2108 !$acc kernels present( nzb_w_inner, rdf, tend, tw_m, w, w_p ) 2109 !$acc loop 2110 DO i = nxl, nxr 2111 DO j = nys, nyn 2112 !$acc loop vector( 32 ) 2113 DO k = 1, nzt-1 2114 IF ( k > nzb_w_inner(j,i) ) THEN 2115 w_p(k,j,i) = w(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + & 2116 tsc(3) * tw_m(k,j,i) ) & 2117 - tsc(5) * rdf(k) * w(k,j,i) 2118 ! 2119 !-- Tendencies for the next Runge-Kutta step 2120 IF ( runge_step == 1 ) THEN 2121 tw_m(k,j,i) = tend(k,j,i) 2122 ELSEIF ( runge_step == 2 ) THEN 2123 tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i) 2124 ENDIF 2125 ENDIF 2126 ENDDO 2127 ENDDO 2128 ENDDO 2129 !$acc end kernels 2130 2131 CALL cpu_log( log_point(7), 'w-equation', 'stop' ) 2132 !$acc update host( w_p ) 2133 2134 2135 ! 2136 !-- If required, compute prognostic equation for potential temperature 2137 IF ( .NOT. neutral ) THEN 2138 2139 CALL cpu_log( log_point(13), 'pt-equation', 'start' ) 2140 2141 ! 2142 !-- pt-tendency terms with communication 2143 sbt = tsc(2) 2144 IF ( scalar_advec == 'bc-scheme' ) THEN 2145 2146 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 2147 ! 2148 !-- Bott-Chlond scheme always uses Euler time step. Thus: 2149 sbt = 1.0 2150 ENDIF 2151 tend = 0.0 2152 CALL advec_s_bc( pt, 'pt' ) 2153 2154 ENDIF 2155 2156 ! 2157 !-- pt-tendency terms with no communication 2158 IF ( scalar_advec /= 'bc-scheme' ) THEN 2159 tend = 0.0 2160 IF ( timestep_scheme(1:5) == 'runge' ) THEN 2161 IF ( ws_scheme_sca ) THEN 2162 CALL advec_s_ws_acc( pt, 'pt' ) 2163 ELSE 2164 tend = 0.0 ! to be removed later?? 2165 CALL advec_s_pw( pt ) 2166 ENDIF 2167 ELSE 2168 CALL advec_s_up( pt ) 2169 ENDIF 2170 ENDIF 2171 2172 CALL diffusion_s_acc( pt, shf, tswst, wall_heatflux ) 2173 2174 ! 2175 !-- If required compute heating/cooling due to long wave radiation processes 2176 IF ( radiation ) THEN 2177 CALL calc_radiation 2178 ENDIF 2179 2180 ! 2181 !-- If required compute impact of latent heat due to precipitation 2182 IF ( precipitation ) THEN 2183 CALL impact_of_latent_heat 2184 ENDIF 2185 2186 ! 2187 !-- Consideration of heat sources within the plant canopy 2188 IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN 2189 CALL plant_canopy_model( 4 ) 2190 ENDIF 2191 2192 ! 2193 !-- If required compute influence of large-scale subsidence/ascent 2194 IF ( large_scale_subsidence ) THEN 2195 CALL subsidence( tend, pt, pt_init ) 2196 ENDIF 2197 2198 CALL user_actions( 'pt-tendency' ) 2199 2200 ! 2201 !-- Prognostic equation for potential temperature 2202 !$acc kernels present( nzb_s_inner, rdf_sc, ptdf_x, ptdf_y, pt_init ) & 2203 !$acc present( tend, tpt_m, pt, pt_p ) 2204 !$acc loop 2205 DO i = nxl, nxr 2206 DO j = nys, nyn 2207 !$acc loop vector( 32 ) 2208 DO k = 1, nzt 2209 IF ( k > nzb_s_inner(j,i) ) THEN 2210 pt_p(k,j,i) = pt(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + & 2211 tsc(3) * tpt_m(k,j,i) ) & 2212 - tsc(5) * ( pt(k,j,i) - pt_init(k) ) *& 2213 ( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) ) 2214 ! 2215 !-- Tendencies for the next Runge-Kutta step 2216 IF ( runge_step == 1 ) THEN 2217 tpt_m(k,j,i) = tend(k,j,i) 2218 ELSEIF ( runge_step == 2 ) THEN 2219 tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tpt_m(k,j,i) 2220 ENDIF 2221 ENDIF 2222 ENDDO 2223 ENDDO 2224 ENDDO 2225 !$acc end kernels 2226 2227 CALL cpu_log( log_point(13), 'pt-equation', 'stop' ) 2228 !$acc update host( pt_p ) 2229 2230 ENDIF 2231 2232 ! 2233 !-- If required, compute prognostic equation for salinity 2234 IF ( ocean ) THEN 2235 2236 CALL cpu_log( log_point(37), 'sa-equation', 'start' ) 2237 2238 ! 2239 !-- sa-tendency terms with communication 2240 sbt = tsc(2) 2241 IF ( scalar_advec == 'bc-scheme' ) THEN 2242 2243 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 2244 ! 2245 !-- Bott-Chlond scheme always uses Euler time step. Thus: 2246 sbt = 1.0 2247 ENDIF 2248 tend = 0.0 2249 CALL advec_s_bc( sa, 'sa' ) 2250 2251 ENDIF 2252 2253 ! 2254 !-- sa-tendency terms with no communication 2255 IF ( scalar_advec /= 'bc-scheme' ) THEN 2256 tend = 0.0 2257 IF ( timestep_scheme(1:5) == 'runge' ) THEN 2258 IF ( ws_scheme_sca ) THEN 2259 CALL advec_s_ws( sa, 'sa' ) 2260 ELSE 2261 CALL advec_s_pw( sa ) 2262 ENDIF 2263 ELSE 2264 CALL advec_s_up( sa ) 2265 ENDIF 2266 ENDIF 2267 2268 CALL diffusion_s( sa, saswsb, saswst, wall_salinityflux ) 2269 2270 CALL user_actions( 'sa-tendency' ) 2271 2272 ! 2273 !-- Prognostic equation for salinity 2274 DO i = nxl, nxr 2275 DO j = nys, nyn 2276 DO k = nzb_s_inner(j,i)+1, nzt 2277 sa_p(k,j,i) = sa(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + & 2278 tsc(3) * tsa_m(k,j,i) ) & 2279 - tsc(5) * rdf_sc(k) * & 2280 ( sa(k,j,i) - sa_init(k) ) 2281 IF ( sa_p(k,j,i) < 0.0 ) sa_p(k,j,i) = 0.1 * sa(k,j,i) 2282 ! 2283 !-- Tendencies for the next Runge-Kutta step 2284 IF ( runge_step == 1 ) THEN 2285 tsa_m(k,j,i) = tend(k,j,i) 2286 ELSEIF ( runge_step == 2 ) THEN 2287 tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tsa_m(k,j,i) 2288 ENDIF 2289 ENDDO 2290 ENDDO 2291 ENDDO 2292 2293 CALL cpu_log( log_point(37), 'sa-equation', 'stop' ) 2294 2295 ! 2296 !-- Calculate density by the equation of state for seawater 2297 CALL cpu_log( log_point(38), 'eqns-seawater', 'start' ) 2298 CALL eqn_state_seawater 2299 CALL cpu_log( log_point(38), 'eqns-seawater', 'stop' ) 2300 2301 ENDIF 2302 2303 ! 2304 !-- If required, compute prognostic equation for total water content / scalar 2305 IF ( humidity .OR. passive_scalar ) THEN 2306 2307 CALL cpu_log( log_point(29), 'q/s-equation', 'start' ) 2308 2309 ! 2310 !-- Scalar/q-tendency terms with communication 2311 sbt = tsc(2) 2312 IF ( scalar_advec == 'bc-scheme' ) THEN 2313 2314 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 2315 ! 2316 !-- Bott-Chlond scheme always uses Euler time step. Thus: 2317 sbt = 1.0 2318 ENDIF 2319 tend = 0.0 2320 CALL advec_s_bc( q, 'q' ) 2321 2322 ENDIF 2323 2324 ! 2325 !-- Scalar/q-tendency terms with no communication 2326 IF ( scalar_advec /= 'bc-scheme' ) THEN 2327 tend = 0.0 2328 IF ( timestep_scheme(1:5) == 'runge' ) THEN 2329 IF ( ws_scheme_sca ) THEN 2330 CALL advec_s_ws( q, 'q' ) 2331 ELSE 2332 CALL advec_s_pw( q ) 2333 ENDIF 2334 ELSE 2335 CALL advec_s_up( q ) 2336 ENDIF 2337 ENDIF 2338 2339 CALL diffusion_s( q, qsws, qswst, wall_qflux ) 2340 2341 ! 2342 !-- If required compute decrease of total water content due to 2343 !-- precipitation 2344 IF ( precipitation ) THEN 2345 CALL calc_precipitation 2346 ENDIF 2347 2348 ! 2349 !-- Sink or source of scalar concentration due to canopy elements 2350 IF ( plant_canopy ) CALL plant_canopy_model( 5 ) 2351 2352 ! 2353 !-- If required compute influence of large-scale subsidence/ascent 2354 IF ( large_scale_subsidence ) THEN 2355 CALL subsidence( tend, q, q_init ) 2356 ENDIF 2357 2358 CALL user_actions( 'q-tendency' ) 2359 2360 ! 2361 !-- Prognostic equation for total water content / scalar 2362 DO i = nxl, nxr 2363 DO j = nys, nyn 2364 DO k = nzb_s_inner(j,i)+1, nzt 2365 q_p(k,j,i) = q(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + & 2366 tsc(3) * tq_m(k,j,i) ) & 2367 - tsc(5) * rdf_sc(k) * & 2368 ( q(k,j,i) - q_init(k) ) 2369 IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i) 2370 ! 2371 !-- Tendencies for the next Runge-Kutta step 2372 IF ( runge_step == 1 ) THEN 2373 tq_m(k,j,i) = tend(k,j,i) 2374 ELSEIF ( runge_step == 2 ) THEN 2375 tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i) 2376 ENDIF 2377 ENDDO 2378 ENDDO 2379 ENDDO 2380 2381 CALL cpu_log( log_point(29), 'q/s-equation', 'stop' ) 2382 2383 ENDIF 2384 2385 ! 2386 !-- If required, compute prognostic equation for turbulent kinetic 2387 !-- energy (TKE) 2388 IF ( .NOT. constant_diffusion ) THEN 2389 2390 CALL cpu_log( log_point(16), 'tke-equation', 'start' ) 2391 2392 ! 2393 !-- TKE-tendency terms with communication 2394 CALL production_e_init 2395 2396 sbt = tsc(2) 2397 IF ( .NOT. use_upstream_for_tke ) THEN 2398 IF ( scalar_advec == 'bc-scheme' ) THEN 2399 2400 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 2401 ! 2402 !-- Bott-Chlond scheme always uses Euler time step. Thus: 2403 sbt = 1.0 2404 ENDIF 2405 tend = 0.0 2406 CALL advec_s_bc( e, 'e' ) 2407 2408 ENDIF 2409 ENDIF 2410 2411 ! 2412 !-- TKE-tendency terms with no communication 2413 IF ( scalar_advec /= 'bc-scheme' .OR. use_upstream_for_tke ) THEN 2414 IF ( use_upstream_for_tke ) THEN 2415 tend = 0.0 2416 CALL advec_s_up( e ) 2417 ELSE 2418 IF ( timestep_scheme(1:5) == 'runge' ) THEN 2419 IF ( ws_scheme_sca ) THEN 2420 CALL advec_s_ws_acc( e, 'e' ) 2421 ELSE 2422 tend = 0.0 ! to be removed later?? 2423 CALL advec_s_pw( e ) 2424 ENDIF 2425 ELSE 2426 tend = 0.0 ! to be removed later?? 2427 CALL advec_s_up( e ) 2428 ENDIF 2429 ENDIF 2430 ENDIF 2431 2432 IF ( .NOT. humidity ) THEN 2433 IF ( ocean ) THEN 2434 CALL diffusion_e( prho, prho_reference ) 2435 ELSE 2436 CALL diffusion_e_acc( pt, pt_reference ) 2437 ENDIF 2438 ELSE 2439 CALL diffusion_e( vpt, pt_reference ) 2440 ENDIF 2441 2442 CALL production_e_acc 2443 2444 ! 2445 !-- Additional sink term for flows through plant canopies 2446 IF ( plant_canopy ) CALL plant_canopy_model( 6 ) 2447 CALL user_actions( 'e-tendency' ) 2448 2449 ! 2450 !-- Prognostic equation for TKE. 2451 !-- Eliminate negative TKE values, which can occur due to numerical 2452 !-- reasons in the course of the integration. In such cases the old TKE 2453 !-- value is reduced by 90%. 2454 !$acc kernels present( e, e_p, nzb_s_inner, tend, te_m ) 2455 !$acc loop 2456 DO i = nxl, nxr 2457 DO j = nys, nyn 2458 !$acc loop vector( 32 ) 2459 DO k = 1, nzt 2460 IF ( k > nzb_s_inner(j,i) ) THEN 2461 e_p(k,j,i) = e(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + & 2462 tsc(3) * te_m(k,j,i) ) 2463 IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i) 2464 ! 2465 !-- Tendencies for the next Runge-Kutta step 2466 IF ( runge_step == 1 ) THEN 2467 te_m(k,j,i) = tend(k,j,i) 2468 ELSEIF ( runge_step == 2 ) THEN 2469 te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i) 2470 ENDIF 2471 ENDIF 2472 ENDDO 2473 ENDDO 2474 ENDDO 2475 !$acc end kernels 2476 2477 CALL cpu_log( log_point(16), 'tke-equation', 'stop' ) 2478 !$acc update host( e_p ) 2479 2480 ENDIF 2481 2482 2483 END SUBROUTINE prognostic_equations_acc 2484 2485 1902 2486 END MODULE prognostic_equations_mod
Note: See TracChangeset
for help on using the changeset viewer.
