Changeset 59
- Timestamp:
- Mar 9, 2007 2:28:00 PM (18 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/advec_particles.f90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/advec_particles.f90
r57 r59 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! Particle reflection at vertical walls (by Jin Zhang), regarding vertical 7 ! walls in the SGS model, 8 ! + user interface user_advec_particles 6 ! Particle reflection at vertical walls implemented in new subroutine 7 ! particle_boundary_conds, 8 ! vertical walls are regarded in the SGS model, 9 ! + user_advec_particles 9 10 ! TEST: PRINT statements on unit 9 (commented out) 10 11 ! … … 1834 1835 ENDDO ! advection loop 1835 1836 1836 1837 1837 ! 1838 1838 !-- Particle reflection from walls 1839 CALL cpu_log( log_point_s(48), 'advec_part_refle', 'start' ) 1840 1841 DO n = 1, number_of_particles 1842 1843 !------------------------------------------------------- 1844 i2 = (particles(n)%x + 0.5*dx) * ddx 1845 j2 = (particles(n)%y + 0.5*dy) * ddy 1846 k2 = particles(n)%z / dz + 1 1847 1848 particles_x = particles(n)%x 1849 particles_y = particles(n)%y 1850 particles_z = particles(n)%z 1851 1852 1853 IF (k2<=nzb_s_inner(j2,i2).and.nzb_s_inner(j2,i2)/=0) THEN 1854 1855 old_positions_x = particles(n)%x - particles(n)%speed_x * dt_particle 1856 old_positions_y = particles(n)%y - particles(n)%speed_y * dt_particle 1857 old_positions_z = particles(n)%z - particles(n)%speed_z * dt_particle 1858 i1= (old_positions_x +0.5*dx) * ddx 1859 j1 =(old_positions_y +0.5*dy) * ddy 1860 k1 = old_positions_z / dz 1861 1862 1863 !-- CASE 1 1864 1865 IF ( particles(n)%x > old_positions_x .and. & 1866 particles(n)%y > old_positions_y )THEN 1867 1868 1869 1870 1871 t_index = 1 1872 1873 Do i = i1, i2 1874 xline = i*dx+0.5*dx 1875 1876 t(t_index) =(xline-old_positions_x)/(particles(n)%x-old_positions_x) 1877 t_index = t_index + 1 1878 ENDDO 1879 1880 Do j = j1, j2 1881 yline = j*dy+0.5*dy 1882 1883 t(t_index) =(yline-old_positions_y)/(particles(n)%y-old_positions_y) 1884 t_index = t_index + 1 1885 ENDDO 1886 1887 IF ( particles(n)%z < old_positions_z ) THEN 1888 Do k = k1, k2 , -1 1889 zline = k*dz 1890 t(t_index) = (old_positions_z-zline)/(old_positions_z-particles(n)%z) 1891 t_index = t_index + 1 1892 ENDDO 1893 ENDIF 1894 t_index_number = t_index-1 1895 !-- sorting t 1896 inc = 1 1897 jr = 1 1898 DO WHILE ( inc <= t_index_number ) 1899 inc = 3 * inc + 1 1900 ENDDO 1901 1902 DO WHILE ( inc > 1 ) 1903 inc = inc / 3 1904 DO ir = inc+1, t_index_number 1905 tmp_t = t(ir) 1906 jr = ir 1907 1908 DO WHILE ( t(jr-inc) > tmp_t ) 1909 t(jr) = t(jr-inc) 1910 jr = jr - inc 1911 IF ( jr <= inc ) THEN 1912 EXIT 1913 ENDIF 1914 ENDDO 1915 t(jr) = tmp_t 1916 1917 ENDDO 1918 ENDDO 1919 1920 1921 1922 1923 Do t_index = 1, t_index_number 1924 1925 positions_x = old_positions_x + t(t_index)*(particles_x - old_positions_x) 1926 positions_y = old_positions_y + t(t_index)*(particles_y - old_positions_y) 1927 positions_z = old_positions_z + t(t_index)*(particles_z - old_positions_z) 1928 1929 1930 1931 i3 = (positions_x + 0.5*dx) * ddx 1932 j3 = (positions_y + 0.5*dy) * ddy 1933 k3 = positions_z / dz 1934 1935 i5 = positions_x * ddx 1936 j5 = positions_y * ddy 1937 k5 = positions_z / dz 1938 1939 1940 1941 1942 1943 IF (k5 <= nzb_s_inner(j5,i3).and.nzb_s_inner(j5,i3)/=0 ) THEN 1944 1945 IF ( positions_z == nzb_s_inner(j5,i3)*dz.and. & 1946 k3 == nzb_s_inner(j5,i3) ) THEN 1947 1948 particles(n)%z = 2*positions_z - particles_z 1949 particles(n)%speed_z = - particles(n)%speed_z 1950 1951 IF ( use_sgs_for_particles ) THEN 1952 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 1953 ENDIF 1954 GOTO 999 1955 ENDIF 1956 ENDIF 1957 1958 IF (k5 <= nzb_s_inner(j3,i5).and.nzb_s_inner(j3,i5)/=0 ) THEN 1959 1960 IF ( positions_z == nzb_s_inner(j3,i5)*dz.and. & 1961 k3 == nzb_s_inner(j3,i5) ) THEN 1962 1963 particles(n)%z = 2*positions_z - particles_z 1964 particles(n)%speed_z = - particles(n)%speed_z 1965 1966 IF ( use_sgs_for_particles ) THEN 1967 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 1968 ENDIF 1969 1970 GOTO 999 1971 ENDIF 1972 ENDIF 1973 1974 IF (k3 <= nzb_s_inner(j3,i3).and.nzb_s_inner(j3,i3)/=0 ) THEN 1975 1976 IF ( positions_z == nzb_s_inner(j3,i3)*dz.and. & 1977 k3 == nzb_s_inner(j3,i3) ) THEN 1978 1979 particles(n)%z = 2*positions_z - particles_z 1980 particles(n)%speed_z = - particles(n)%speed_z 1981 1982 IF ( use_sgs_for_particles ) THEN 1983 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 1984 ENDIF 1985 1986 GOTO 999 1987 ENDIF 1988 1989 IF ( positions_y == j3*dy-0.5*dy .and. & 1990 positions_z < nzb_s_inner(j3,i3)*dz) THEN 1991 1992 particles(n)%y = 2*positions_y - particles_y 1993 particles(n)%speed_y = - particles(n)%speed_y 1994 1995 IF ( use_sgs_for_particles ) THEN 1996 particles(n)%speed_y_sgs = - particles(n)%speed_y_sgs 1997 ENDIF 1998 1999 GOTO 999 2000 ENDIF 2001 2002 2003 IF ( positions_x == i3*dx-0.5*dx.and. & 2004 positions_z < nzb_s_inner(j3,i3)*dz ) THEN 2005 2006 particles(n)%x = 2*positions_x - particles_x 2007 particles(n)%speed_x = - particles(n)%speed_x 2008 2009 IF ( use_sgs_for_particles ) THEN 2010 particles(n)%speed_x_sgs = - particles(n)%speed_x_sgs 2011 GOTO 999 2012 ENDIF 2013 2014 ENDIF 2015 2016 ENDIF 2017 ENDDO 2018 ENDIF 2019 2020 !-- CASE 2 2021 2022 IF ( particles(n)%x > old_positions_x .and. & 2023 particles(n)%y < old_positions_y )THEN 2024 2025 2026 t_index = 1 2027 2028 Do i = i1, i2 2029 xline = i*dx+0.5*dx 2030 2031 t(t_index) =(xline-old_positions_x)/(particles(n)%x-old_positions_x) 2032 t_index = t_index + 1 2033 ENDDO 2034 2035 Do j = j1, j2,-1 2036 yline = j*dy-0.5*dy 2037 2038 t(t_index) =(old_positions_y-yline)/(old_positions_y-particles(n)%y) 2039 t_index = t_index + 1 2040 ENDDO 2041 2042 IF ( particles(n)%z < old_positions_z ) THEN 2043 Do k = k1, k2 , -1 2044 zline = k*dz 2045 t(t_index) = (old_positions_z-zline)/(old_positions_z-particles(n)%z) 2046 t_index = t_index + 1 2047 ENDDO 2048 ENDIF 2049 t_index_number = t_index-1 2050 2051 !-- sorting t 2052 inc = 1 2053 jr = 1 2054 DO WHILE ( inc <= t_index_number ) 2055 inc = 3 * inc + 1 2056 ENDDO 2057 2058 DO WHILE ( inc > 1 ) 2059 inc = inc / 3 2060 DO ir = inc+1, t_index_number 2061 tmp_t = t(ir) 2062 jr = ir 2063 2064 DO WHILE ( t(jr-inc) > tmp_t ) 2065 t(jr) = t(jr-inc) 2066 jr = jr - inc 2067 IF ( jr <= inc ) THEN 2068 EXIT 2069 ENDIF 2070 ENDDO 2071 t(jr) = tmp_t 2072 2073 ENDDO 2074 ENDDO 2075 2076 2077 2078 2079 2080 2081 Do t_index = 1, t_index_number 2082 2083 positions_x = old_positions_x + t(t_index)*(particles_x - old_positions_x) 2084 positions_y = old_positions_y + t(t_index)*(particles_y - old_positions_y) 2085 positions_z = old_positions_z + t(t_index)*(particles_z - old_positions_z) 2086 2087 2088 2089 i3 = (positions_x + 0.5*dx) * ddx 2090 j3 = (positions_y + 0.5*dy) * ddy 2091 k3 = positions_z / dz 2092 2093 i5 = positions_x * ddx 2094 j5 = positions_y * ddy 2095 k5 = positions_z / dz 2096 2097 2098 2099 IF (k5 <= nzb_s_inner(j3,i5).and.nzb_s_inner(j3,i5)/=0 ) THEN 2100 2101 IF ( positions_z == nzb_s_inner(j3,i5)*dz.and. & 2102 k3 == nzb_s_inner(j3,i5) ) THEN 2103 2104 particles(n)%z = 2*positions_z - particles_z 2105 particles(n)%speed_z = - particles(n)%speed_z 2106 2107 IF ( use_sgs_for_particles ) THEN 2108 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 2109 ENDIF 2110 2111 GOTO 999 2112 ENDIF 2113 ENDIF 2114 2115 IF (k3 <= nzb_s_inner(j3,i3).and.nzb_s_inner(j3,i3)/=0 ) THEN 2116 2117 IF ( positions_z == nzb_s_inner(j3,i3)*dz.and. & 2118 k3 == nzb_s_inner(j3,i3) ) THEN 2119 2120 particles(n)%z = 2*positions_z - particles_z 2121 particles(n)%speed_z = - particles(n)%speed_z 2122 2123 IF ( use_sgs_for_particles ) THEN 2124 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 2125 ENDIF 2126 2127 GOTO 999 2128 ENDIF 2129 2130 2131 IF ( positions_x == i3*dx-0.5*dx.and.& 2132 positions_z < nzb_s_inner(j3,i3)*dz ) THEN 2133 2134 particles(n)%x = 2*positions_x - particles_x 2135 particles(n)%speed_x = - particles(n)%speed_x 2136 2137 IF ( use_sgs_for_particles ) THEN 2138 particles(n)%speed_x_sgs = - particles(n)%speed_x_sgs 2139 ENDIF 2140 2141 GOTO 999 2142 ENDIF 2143 2144 ENDIF 2145 2146 2147 IF (k5 <= nzb_s_inner(j5,i3).and.nzb_s_inner(j5,i3)/=0 ) THEN 2148 2149 IF ( positions_z == nzb_s_inner(j5,i3)*dz.and. & 2150 k3 == nzb_s_inner(j5,i3) ) THEN 2151 2152 particles(n)%z = 2*positions_z - particles_z 2153 particles(n)%speed_z = - particles(n)%speed_z 2154 IF ( use_sgs_for_particles ) THEN 2155 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 2156 ENDIF 2157 2158 GOTO 999 2159 ENDIF 2160 2161 IF ( positions_y == j5 * dy + 0.5*dy.and.& 2162 positions_z < nzb_s_inner(j5,i3)*dz ) THEN 2163 2164 particles(n)%y = 2*positions_y - particles_y 2165 particles(n)%speed_y = - particles(n)%speed_y 2166 2167 IF ( use_sgs_for_particles ) THEN 2168 particles(n)%speed_y_sgs = - particles(n)%speed_y_sgs 2169 ENDIF 2170 GOTO 999 2171 ENDIF 2172 2173 ENDIF 2174 ENDDO 2175 ENDIF 2176 2177 2178 !-- CASE 3 2179 2180 2181 2182 IF ( particles(n)%x < old_positions_x .and. & 2183 particles(n)%y > old_positions_y )THEN 2184 2185 2186 t_index = 1 2187 2188 Do i = i1, i2,-1 2189 xline = i*dx-0.5*dx 2190 2191 t(t_index) =(old_positions_x-xline)/(old_positions_x-particles(n)%x) 2192 t_index = t_index + 1 2193 ENDDO 2194 2195 Do j = j1, j2 2196 yline = j*dy+0.5*dy 2197 2198 t(t_index) =(yline-old_positions_y)/(particles(n)%y-old_positions_y) 2199 t_index = t_index + 1 2200 ENDDO 2201 2202 IF ( particles(n)%z < old_positions_z ) THEN 2203 Do k = k1, k2 , -1 2204 zline = k*dz 2205 t(t_index) = (old_positions_z-zline)/(old_positions_z-particles(n)%z) 2206 t_index = t_index + 1 2207 ENDDO 2208 ENDIF 2209 2210 t_index_number = t_index-1 2211 2212 !-- sorting t 2213 inc = 1 2214 jr = 1 2215 DO WHILE ( inc <= t_index_number ) 2216 inc = 3 * inc + 1 2217 ENDDO 2218 2219 DO WHILE ( inc > 1 ) 2220 inc = inc / 3 2221 DO ir = inc+1, t_index_number 2222 tmp_t = t(ir) 2223 jr = ir 2224 2225 DO WHILE ( t(jr-inc) > tmp_t ) 2226 t(jr) = t(jr-inc) 2227 jr = jr - inc 2228 IF ( jr <= inc ) THEN 2229 EXIT 2230 ENDIF 2231 ENDDO 2232 t(jr) = tmp_t 2233 2234 ENDDO 2235 ENDDO 2236 2237 2238 2239 Do t_index = 1, t_index_number 2240 2241 positions_x = old_positions_x + t(t_index)*(particles_x - old_positions_x) 2242 positions_y = old_positions_y + t(t_index)*(particles_y - old_positions_y) 2243 positions_z = old_positions_z + t(t_index)*(particles_z - old_positions_z) 2244 2245 2246 2247 i3 = (positions_x + 0.5*dx) * ddx 2248 j3 = (positions_y + 0.5*dy) * ddy 2249 k3 = positions_z / dz 2250 2251 i5 = positions_x * ddx 2252 j5 = positions_y * ddy 2253 k5 = positions_z / dz 2254 2255 2256 2257 IF (k5 <= nzb_s_inner(j5,i3).and.nzb_s_inner(j5,i3)/=0 ) THEN 2258 2259 IF ( positions_z == nzb_s_inner(j5,i3)*dz.and. & 2260 k3 == nzb_s_inner(j5,i3) ) THEN 2261 2262 particles(n)%z = 2*positions_z - particles_z 2263 particles(n)%speed_z = - particles(n)%speed_z 2264 2265 IF ( use_sgs_for_particles ) THEN 2266 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 2267 ENDIF 2268 2269 GOTO 999 2270 ENDIF 2271 ENDIF 2272 2273 2274 2275 IF (k3 <= nzb_s_inner(j3,i3).and.nzb_s_inner(j3,i3)/=0 ) THEN 2276 2277 2278 IF ( positions_z == nzb_s_inner(j3,i3)*dz.and. & 2279 k3 == nzb_s_inner(j3,i3) ) THEN 2280 2281 particles(n)%z = 2*positions_z - particles_z 2282 particles(n)%speed_z = - particles(n)%speed_z 2283 2284 IF ( use_sgs_for_particles ) THEN 2285 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 2286 ENDIF 2287 2288 GOTO 999 2289 ENDIF 2290 2291 IF ( positions_y == j3*dy-0.5*dy .and. & 2292 positions_z < nzb_s_inner(j3,i3)*dz) THEN 2293 2294 particles(n)%y = 2*positions_y - particles_y 2295 particles(n)%speed_y = - particles(n)%speed_y 2296 2297 IF ( use_sgs_for_particles ) THEN 2298 particles(n)%speed_y_sgs = - particles(n)%speed_y_sgs 2299 ENDIF 2300 2301 GOTO 999 2302 ENDIF 2303 2304 ENDIF 2305 2306 2307 IF (k5 <= nzb_s_inner(j3,i5).and.nzb_s_inner(j3,i5)/=0 ) THEN 2308 2309 IF ( positions_z == nzb_s_inner(j3,i5)*dz.and. & 2310 k3 == nzb_s_inner(j3,i5) ) THEN 2311 2312 particles(n)%z = 2*positions_z - particles_z 2313 particles(n)%speed_z = - particles(n)%speed_z 2314 2315 IF ( use_sgs_for_particles ) THEN 2316 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 2317 ENDIF 2318 2319 2320 GOTO 999 2321 ENDIF 2322 2323 IF ( positions_x == i5 * dx + 0.5*dx .and. & 2324 positions_z < nzb_s_inner(j3,i5)*dz ) THEN 2325 2326 particles(n)%x= 2*positions_x - particles_x 2327 particles(n)%speed_x = - particles(n)%speed_x 2328 2329 IF ( use_sgs_for_particles ) THEN 2330 particles(n)%speed_x_sgs = - particles(n)%speed_x_sgs 2331 ENDIF 2332 2333 GOTO 999 2334 ENDIF 2335 2336 ENDIF 2337 ENDDO 2338 ENDIF 2339 2340 2341 !-- CASE 4 2342 2343 2344 IF ( particles(n)%x < old_positions_x .and. & 2345 particles(n)%y < old_positions_y )THEN 2346 2347 2348 t_index = 1 2349 2350 Do i = i1, i2,-1 2351 xline = i*dx-0.5*dx 2352 2353 t(t_index) =(old_positions_x-xline)/(old_positions_x-particles(n)%x) 2354 t_index = t_index + 1 2355 2356 2357 2358 ENDDO 2359 2360 Do j = j1, j2,-1 2361 yline = j*dy-0.5*dy 2362 2363 t(t_index) =(old_positions_y-yline)/(old_positions_y-particles(n)%y) 2364 t_index = t_index + 1 2365 2366 2367 ENDDO 2368 2369 IF ( particles(n)%z < old_positions_z ) THEN 2370 Do k = k1, k2 , -1 2371 zline = k*dz 2372 t(t_index) = (old_positions_z-zline)/(old_positions_z-particles(n)%z) 2373 t_index = t_index + 1 2374 2375 2376 ENDDO 2377 ENDIF 2378 2379 t_index_number = t_index-1 2380 2381 !-- sorting t 2382 inc = 1 2383 jr = 1 2384 DO WHILE ( inc <= t_index_number ) 2385 inc = 3 * inc + 1 2386 ENDDO 2387 2388 DO WHILE ( inc > 1 ) 2389 inc = inc / 3 2390 DO ir = inc+1, t_index_number 2391 tmp_t = t(ir) 2392 jr = ir 2393 2394 DO WHILE ( t(jr-inc) > tmp_t ) 2395 t(jr) = t(jr-inc) 2396 jr = jr - inc 2397 IF ( jr <= inc ) THEN 2398 EXIT 2399 ENDIF 2400 ENDDO 2401 t(jr) = tmp_t 2402 2403 ENDDO 2404 ENDDO 2405 2406 2407 Do t_index = 1, t_index_number 2408 2409 positions_x = old_positions_x + t(t_index)*(particles_x - old_positions_x) 2410 positions_y = old_positions_y + t(t_index)*(particles_y - old_positions_y) 2411 positions_z = old_positions_z + t(t_index)*(particles_z - old_positions_z) 2412 2413 2414 2415 i3 = (positions_x + 0.5*dx) * ddx 2416 j3 = (positions_y + 0.5*dy) * ddy 2417 k3 = positions_z / dz 2418 2419 i5 = positions_x * ddx 2420 j5 = positions_y * ddy 2421 k5 = positions_z / dz 2422 2423 2424 2425 2426 IF (k3 <= nzb_s_inner(j3,i3).and.nzb_s_inner(j3,i3)/=0 ) THEN 2427 2428 IF ( positions_z == nzb_s_inner(j3,i3)*dz.and. & 2429 k3 == nzb_s_inner(j3,i3) ) THEN 2430 2431 particles(n)%z = 2*positions_z - particles_z 2432 particles(n)%speed_z = - particles(n)%speed_z 2433 2434 IF ( use_sgs_for_particles ) THEN 2435 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 2436 ENDIF 2437 2438 2439 GOTO 999 2440 ENDIF 2441 ENDIF 2442 2443 2444 IF (k5 <= nzb_s_inner(j3,i5).and.nzb_s_inner(j3,i5)/=0 ) THEN 2445 2446 IF ( positions_z == nzb_s_inner(j3,i5)*dz.and. & 2447 k3 == nzb_s_inner(j3,i5) ) THEN 2448 2449 particles(n)%z = 2*positions_z - particles_z 2450 particles(n)%speed_z = - particles(n)%speed_z 2451 2452 IF ( use_sgs_for_particles ) THEN 2453 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 2454 ENDIF 2455 2456 GOTO 999 2457 ENDIF 2458 2459 IF ( positions_x == i5 * dx + 0.5*dx .and. nzb_s_inner(j3,i5)/=0.and.& 2460 positions_z < nzb_s_inner(j3,i5)*dz) THEN 2461 2462 particles(n)%x= 2*positions_x - particles_x 2463 particles(n)%speed_x = - particles(n)%speed_x 2464 2465 IF ( use_sgs_for_particles ) THEN 2466 particles(n)%speed_x_sgs = - particles(n)%speed_x_sgs 2467 ENDIF 2468 GOTO 999 2469 ENDIF 2470 2471 ENDIF 2472 2473 IF (k5 <= nzb_s_inner(j5,i3).and.nzb_s_inner(j5,i3)/=0) THEN 2474 2475 2476 IF ( positions_z == nzb_s_inner(j5,i3)*dz.and. & 2477 k5 == nzb_s_inner(j5,i3) ) THEN 2478 2479 particles(n)%z = 2*positions_z - particles_z 2480 particles(n)%speed_z = - particles(n)%speed_z 2481 2482 IF ( use_sgs_for_particles ) THEN 2483 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 2484 ENDIF 2485 2486 GOTO 999 2487 ENDIF 2488 2489 IF ( positions_y == j5 * dy + 0.5*dy .and. nzb_s_inner(j5,i3)/=0.and.& 2490 positions_z < nzb_s_inner(j5,i3)*dz ) THEN 2491 2492 particles(n)%y= 2*positions_y - particles_y 2493 particles(n)%speed_y = - particles(n)%speed_y 2494 2495 IF ( use_sgs_for_particles ) THEN 2496 particles(n)%speed_y_sgs = - particles(n)%speed_y_sgs 2497 ENDIF 2498 2499 GOTO 999 2500 ENDIF 2501 2502 2503 ENDIF 2504 ENDDO 2505 ENDIF 2506 2507 999 CONTINUE 2508 2509 ENDIF 2510 2511 !--------------------------------------------------- 2512 2513 ENDDO ! particle reflection from walls 2514 2515 CALL cpu_log( log_point_s(48), 'advec_part_refle', 'stop' ) 1839 CALL particle_boundary_conds 2516 1840 2517 1841 !
Note: See TracChangeset
for help on using the changeset viewer.
