Changeset 2019 for palm/trunk/SOURCE/pmc_interface_mod.f90
- Timestamp:
- Sep 30, 2016 1:40:09 PM (8 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/pmc_interface_mod.f90 (modified) (22 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/pmc_interface_mod.f90
r2005 r2019 21 21 ! Current revisions: 22 22 ! ------------------ 23 ! 23 ! Bugfixes mainly in determining the anterpolation index bounds. These errors 24 ! were detected when first time tested using 3:1 grid-spacing. 24 25 ! 25 26 ! Former revisions: … … 126 127 USE arrays_3d, & 127 128 ONLY: dzu, dzw, e, e_p, e_1, e_2, pt, pt_p, pt_1, pt_2, q, q_p, q_1, & 128 q_2, s, s_2, u, u_p, u_1, u_2, v, v_p, v_1, v_2, w, w_p, w_1, w_2,&129 zu, zw, z0129 q_2, s, s_2, u, u_p, u_1, u_2, v, v_p, v_1, v_2, w, w_p, w_1, & 130 w_2, zu, zw, z0 130 131 #endif 131 132 … … 1090 1091 !-- This is needed in the anterpolation phase 1091 1092 nhlr = CEILING( xexr / cg%dx ) 1092 xce = coord_x(nxr) + xexr 1093 xce = coord_x(nxr+1) + xexr 1094 !-- One "extra" layer is taken behind the right boundary 1095 !-- because it may be needed in cases of non-integer grid-spacing ratio 1093 1096 DO i = cg%nx, 0 , -1 1094 1097 IF ( cg%coord_x(i) < xce ) THEN … … 1118 1121 !-- This is needed in the anterpolation phase 1119 1122 nhln = CEILING( yexn / cg%dy ) 1120 yce = coord_y(nyn) + yexn 1123 yce = coord_y(nyn+1) + yexn 1124 !-- One "extra" layer is taken behind the north boundary 1125 !-- because it may be needed in cases of non-integer grid-spacing ratio 1121 1126 DO j = cg%ny, 0, -1 1122 1127 IF ( cg%coord_y(j) < yce ) THEN … … 1358 1363 1359 1364 ! 1360 !-- First horizontal walls 1365 !-- First horizontal walls. Note that also logc_w_? and logc_ratio_w_? need to 1366 !-- be allocated and initialized here. 1361 1367 !-- Left boundary 1362 1368 IF ( nest_bound_l ) THEN … … 1364 1370 ALLOCATE( logc_u_l(1:2,nzb:nzt_topo_nestbc_l,nys:nyn) ) 1365 1371 ALLOCATE( logc_v_l(1:2,nzb:nzt_topo_nestbc_l,nys:nyn) ) 1372 ALLOCATE( logc_w_l(1:2,nzb:nzt_topo_nestbc_l,nys:nyn) ) 1366 1373 ALLOCATE( logc_ratio_u_l(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_l,nys:nyn) ) 1367 1374 ALLOCATE( logc_ratio_v_l(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_l,nys:nyn) ) 1375 ALLOCATE( logc_ratio_w_l(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_l,nys:nyn) ) 1368 1376 logc_u_l = 0 1369 1377 logc_v_l = 0 1378 logc_w_l = 0 1370 1379 logc_ratio_u_l = 1.0_wp 1371 1380 logc_ratio_v_l = 1.0_wp 1381 logc_ratio_w_l = 1.0_wp 1372 1382 direction = 1 1373 1383 inc = 1 … … 1404 1414 !-- Right boundary 1405 1415 IF ( nest_bound_r ) THEN 1406 1416 1407 1417 ALLOCATE( logc_u_r(1:2,nzb:nzt_topo_nestbc_r,nys:nyn) ) 1408 1418 ALLOCATE( logc_v_r(1:2,nzb:nzt_topo_nestbc_r,nys:nyn) ) 1419 ALLOCATE( logc_w_r(1:2,nzb:nzt_topo_nestbc_r,nys:nyn) ) 1409 1420 ALLOCATE( logc_ratio_u_r(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_r,nys:nyn) ) 1410 1421 ALLOCATE( logc_ratio_v_r(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_r,nys:nyn) ) 1422 ALLOCATE( logc_ratio_w_r(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_r,nys:nyn) ) 1411 1423 logc_u_r = 0 1412 1424 logc_v_r = 0 1425 logc_w_r = 0 1413 1426 logc_ratio_u_r = 1.0_wp 1414 1427 logc_ratio_v_r = 1.0_wp 1428 logc_ratio_w_r = 1.0_wp 1415 1429 direction = 1 1416 1430 inc = 1 1431 1417 1432 DO j = nys, nyn 1418 1433 ! … … 1449 1464 ALLOCATE( logc_u_s(1:2,nzb:nzt_topo_nestbc_s,nxl:nxr) ) 1450 1465 ALLOCATE( logc_v_s(1:2,nzb:nzt_topo_nestbc_s,nxl:nxr) ) 1466 ALLOCATE( logc_w_s(1:2,nzb:nzt_topo_nestbc_s,nxl:nxr) ) 1451 1467 ALLOCATE( logc_ratio_u_s(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_s,nxl:nxr) ) 1452 1468 ALLOCATE( logc_ratio_v_s(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_s,nxl:nxr) ) 1469 ALLOCATE( logc_ratio_w_s(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_s,nxl:nxr) ) 1453 1470 logc_u_s = 0 1454 1471 logc_v_s = 0 1472 logc_w_s = 0 1455 1473 logc_ratio_u_s = 1.0_wp 1456 1474 logc_ratio_v_s = 1.0_wp 1475 logc_ratio_w_s = 1.0_wp 1457 1476 direction = 1 1458 1477 inc = 1 … … 1492 1511 ALLOCATE( logc_u_n(1:2,nzb:nzt_topo_nestbc_n,nxl:nxr) ) 1493 1512 ALLOCATE( logc_v_n(1:2,nzb:nzt_topo_nestbc_n,nxl:nxr) ) 1513 ALLOCATE( logc_w_n(1:2,nzb:nzt_topo_nestbc_n,nxl:nxr) ) 1494 1514 ALLOCATE( logc_ratio_u_n(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_n,nxl:nxr) ) 1495 1515 ALLOCATE( logc_ratio_v_n(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_n,nxl:nxr) ) 1516 ALLOCATE( logc_ratio_w_n(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_n,nxl:nxr) ) 1496 1517 logc_u_n = 0 1497 1518 logc_v_n = 0 1519 logc_w_n = 0 1498 1520 logc_ratio_u_n = 1.0_wp 1499 1521 logc_ratio_v_n = 1.0_wp 1522 logc_ratio_w_n = 1.0_wp 1500 1523 direction = 1 1501 1524 inc = 1 … … 1538 1561 IF ( nest_bound_l ) THEN 1539 1562 1540 ALLOCATE( logc_w_l(1:2,nzb:nzt_topo_nestbc_l,nys:nyn) ) 1541 ALLOCATE( logc_ratio_w_l(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_l, & 1542 nys:nyn) ) 1543 1544 logc_w_l = 0 1545 logc_ratio_w_l = 1.0_wp 1546 direction = 2 1563 direction = 2 1564 1547 1565 DO j = nys, nyn 1548 1566 DO k = nzb, nzt_topo_nestbc_l … … 1625 1643 IF ( nest_bound_r ) THEN 1626 1644 1627 ALLOCATE( logc_w_r(1:2,nzb:nzt_topo_nestbc_r,nys:nyn) ) 1628 ALLOCATE( logc_ratio_w_r(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_r, & 1629 nys:nyn) ) 1630 logc_w_r = 0 1631 logc_ratio_w_r = 1.0_wp 1632 direction = 2 1645 direction = 2 1633 1646 i = nxr + 1 1634 1647 … … 1709 1722 IF ( nest_bound_s ) THEN 1710 1723 1711 ALLOCATE( logc_w_s(1:2,nzb:nzt_topo_nestbc_s,nxl:nxr) ) 1712 ALLOCATE( logc_ratio_w_s(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_s, & 1713 nxl:nxr) ) 1714 logc_w_s = 0 1715 logc_ratio_w_s = 1.0_wp 1716 direction = 3 1724 direction = 3 1717 1725 1718 1726 DO i = nxl, nxr … … 1795 1803 IF ( nest_bound_n ) THEN 1796 1804 1797 ALLOCATE( logc_w_n(1:2,nzb:nzt_topo_nestbc_n, nxl:nxr) ) 1798 ALLOCATE( logc_ratio_w_n(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_n, & 1799 nxl:nxr) ) 1800 logc_w_n = 0 1801 logc_ratio_w_n = 1.0_wp 1802 direction = 3 1805 direction = 3 1803 1806 j = nyn + 1 1804 1807 … … 2152 2155 !-- index k 2153 2156 kk = 0 2154 DO WHILE ( cg%zu(kk) < zu(nzt) )2157 DO WHILE ( cg%zu(kk) <= zu(nzt) ) 2155 2158 kk = kk + 1 2156 2159 ENDDO … … 2158 2161 2159 2162 kk = 0 2160 DO WHILE ( cg%zw(kk) < zw(nzt-1) )2163 DO WHILE ( cg%zw(kk) <= zw(nzt-1) ) 2161 2164 kk = kk + 1 2162 2165 ENDDO … … 2182 2185 ! 2183 2186 !-- i-indices of u for each ii-index value 2187 !-- ii=icr is redundant for anterpolation 2184 2188 istart = nxlg 2185 DO ii = icl, icr 2189 DO ii = icl, icr-1 2186 2190 i = istart 2187 2191 DO WHILE ( ( coord_x(i) < cg%coord_x(ii) - 0.5_wp * cg%dx ) .AND. & … … 2190 2194 ENDDO 2191 2195 iflu(ii) = MIN( MAX( i, nxlg ), nxrg ) 2192 DO WHILE ( ( coord_x(i) < cg%coord_x(ii) + 0.5_wp * cg%dx ) .AND.&2193 ( i < nxrg ) )2196 DO WHILE ( ( coord_x(i) <= cg%coord_x(ii) + 0.5_wp * cg%dx ) .AND. & 2197 ( i < nxrg+1 ) ) 2194 2198 i = i + 1 2195 2199 ENDDO 2196 ifuu(ii) = MIN( MAX( i , nxlg), nxrg )2200 ifuu(ii) = MIN( MAX( i-1, iflu(ii) ), nxrg ) 2197 2201 istart = iflu(ii) 2198 2202 ENDDO 2203 iflu(icr) = nxrg 2204 ifuu(icr) = nxrg 2199 2205 2200 2206 ! 2201 2207 !-- i-indices of others for each ii-index value 2208 !-- ii=icr is redundant for anterpolation 2202 2209 istart = nxlg 2203 DO ii = icl, icr 2210 DO ii = icl, icr-1 2204 2211 i = istart 2205 2212 DO WHILE ( ( coord_x(i) + 0.5_wp * dx < cg%coord_x(ii) ) .AND. & … … 2208 2215 ENDDO 2209 2216 iflo(ii) = MIN( MAX( i, nxlg ), nxrg ) 2210 DO WHILE ( ( coord_x(i) + 0.5_wp * dx < cg%coord_x(ii) + cg%dx )&2211 .AND. ( i < nxrg ) )2217 DO WHILE ( ( coord_x(i) + 0.5_wp * dx <= cg%coord_x(ii) + cg%dx ) & 2218 .AND. ( i < nxrg+1 ) ) 2212 2219 i = i + 1 2213 2220 ENDDO 2214 ifuo(ii) = MIN( MAX(i,nxlg),nxrg)2221 ifuo(ii) = MIN( MAX( i-1, iflo(ii) ), nxrg ) 2215 2222 istart = iflo(ii) 2216 2223 ENDDO 2224 iflo(icr) = nxrg 2225 ifuo(icr) = nxrg 2217 2226 2218 2227 ! 2219 2228 !-- j-indices of v for each jj-index value 2229 !-- jj=jcn is redundant for anterpolation 2220 2230 jstart = nysg 2221 DO jj = jcs, jcn 2231 DO jj = jcs, jcn-1 2222 2232 j = jstart 2223 2233 DO WHILE ( ( coord_y(j) < cg%coord_y(jj) - 0.5_wp * cg%dy ) .AND. & … … 2226 2236 ENDDO 2227 2237 jflv(jj) = MIN( MAX( j, nysg ), nyng ) 2228 DO WHILE ( ( coord_y(j) < cg%coord_y(jj) + 0.5_wp * cg%dy ) .AND.&2229 ( j < nyng ) )2238 DO WHILE ( ( coord_y(j) <= cg%coord_y(jj) + 0.5_wp * cg%dy ) .AND. & 2239 ( j < nyng+1 ) ) 2230 2240 j = j + 1 2231 2241 ENDDO 2232 jfuv(jj) = MIN( MAX( j , nysg), nyng )2242 jfuv(jj) = MIN( MAX( j-1, jflv(jj) ), nyng ) 2233 2243 jstart = jflv(jj) 2234 2244 ENDDO 2235 2245 jflv(jcn) = nyng 2246 jfuv(jcn) = nyng 2236 2247 ! 2237 2248 !-- j-indices of others for each jj-index value 2249 !-- jj=jcn is redundant for anterpolation 2238 2250 jstart = nysg 2239 DO jj = jcs, jcn 2251 DO jj = jcs, jcn-1 2240 2252 j = jstart 2241 2253 DO WHILE ( ( coord_y(j) + 0.5_wp * dy < cg%coord_y(jj) ) .AND. & … … 2244 2256 ENDDO 2245 2257 jflo(jj) = MIN( MAX( j, nysg ), nyng ) 2246 DO WHILE ( ( coord_y(j) + 0.5_wp * dy < cg%coord_y(jj) + cg%dy )&2247 .AND. ( j < nyng ) )2258 DO WHILE ( ( coord_y(j) + 0.5_wp * dy <= cg%coord_y(jj) + cg%dy ) & 2259 .AND. ( j < nyng+1 ) ) 2248 2260 j = j + 1 2249 2261 ENDDO 2250 jfuo(jj) = MIN( MAX( j , nysg), nyng )2251 jstart = jfl v(jj)2262 jfuo(jj) = MIN( MAX( j-1, jflo(jj) ), nyng ) 2263 jstart = jflo(jj) 2252 2264 ENDDO 2265 jflo(jcn) = nyng 2266 jfuo(jcn) = nyng 2253 2267 2254 2268 ! … … 2259 2273 DO kk = 1, kctw 2260 2274 k = kstart 2261 DO WHILE ( ( zw(k) < cg%zw(kk) - 0.5_wp * cg%dzw(kk) ) .AND. & 2262 ( k < nzt ) ) 2275 DO WHILE ( ( zw(k) < cg%zu(kk) ) .AND. ( k < nzt ) ) 2263 2276 k = k + 1 2264 2277 ENDDO 2265 2278 kflw(kk) = MIN( MAX( k, 1 ), nzt + 1 ) 2266 DO WHILE ( ( zw(k) < cg%zw(kk) + 0.5_wp * cg%dzw(kk+1) ) .AND. & 2267 ( k < nzt ) ) 2279 DO WHILE ( ( zw(k) <= cg%zu(kk+1) ) .AND. ( k < nzt+1 ) ) 2268 2280 k = k + 1 2269 2281 ENDDO 2270 kfuw(kk) = MIN( MAX( k , 1), nzt + 1 )2282 kfuw(kk) = MIN( MAX( k-1, kflw(kk) ), nzt + 1 ) 2271 2283 kstart = kflw(kk) 2272 2284 ENDDO … … 2279 2291 DO kk = 1, kctu 2280 2292 k = kstart 2281 DO WHILE ( ( zu(k) < cg%zu(kk) - 0.5_wp * cg%dzu(kk) ) .AND. & 2282 ( k < nzt ) ) 2293 DO WHILE ( ( zu(k) < cg%zw(kk-1) ) .AND. ( k < nzt ) ) 2283 2294 k = k + 1 2284 2295 ENDDO 2285 2296 kflo(kk) = MIN( MAX( k, 1 ), nzt + 1 ) 2286 DO WHILE ( ( zu(k) < cg%zu(kk) + 0.5_wp * cg%dzu(kk+1) ) .AND. & 2287 ( k < nzt ) ) 2297 DO WHILE ( ( zu(k) <= cg%zw(kk) ) .AND. ( k < nzt+1 ) ) 2288 2298 k = k + 1 2289 2299 ENDDO 2290 kfuo(kk) = MIN( MAX( k-1, 1), nzt + 1 )2300 kfuo(kk) = MIN( MAX( k-1, kflo(kk) ), nzt + 1 ) 2291 2301 kstart = kflo(kk) 2292 2302 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.
