- Timestamp:
- Aug 26, 2015 4:54:34 PM (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/advec_ws.f90
r1586 r1628 493 493 DO k = nzb+1, nzb_max 494 494 495 ibit5 = IBITS(wall_flags_0(k,j ,i),5,1)496 ibit4 = IBITS(wall_flags_0(k,j ,i),4,1)497 ibit3 = IBITS(wall_flags_0(k,j ,i),3,1)495 ibit5 = IBITS(wall_flags_0(k,j-1,i),5,1) 496 ibit4 = IBITS(wall_flags_0(k,j-1,i),4,1) 497 ibit3 = IBITS(wall_flags_0(k,j-1,i),3,1) 498 498 499 499 v_comp = v(k,j,i) - v_gtrans … … 554 554 DO k = nzb+1, nzb_max 555 555 556 ibit2 = IBITS(wall_flags_0(k,j,i ),2,1)557 ibit1 = IBITS(wall_flags_0(k,j,i ),1,1)558 ibit0 = IBITS(wall_flags_0(k,j,i ),0,1)556 ibit2 = IBITS(wall_flags_0(k,j,i-1),2,1) 557 ibit1 = IBITS(wall_flags_0(k,j,i-1),1,1) 558 ibit0 = IBITS(wall_flags_0(k,j,i-1),0,1) 559 559 560 560 u_comp = u(k,j,i) - u_gtrans … … 874 874 !-- correction is needed to overcome numerical instabilities caused 875 875 !-- by a not sufficient reduction of divergences near topography. 876 div = ( u(k,j,i+1) - u(k,j,i) ) * ddx & 877 + ( v(k,j+1,i) - v(k,j,i) ) * ddy & 878 + ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 876 div = ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 ) & 877 - u(k,j,i) * ( IBITS(wall_flags_0(k,j,i-1),0,1) & 878 + IBITS(wall_flags_0(k,j,i-1),1,1) & 879 + IBITS(wall_flags_0(k,j,i-1),2,1) & 880 ) & 881 ) * ddx & 882 + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 ) & 883 - v(k,j,i) * ( IBITS(wall_flags_0(k,j-1,i),3,1) & 884 + IBITS(wall_flags_0(k,j-1,i),4,1) & 885 + IBITS(wall_flags_0(k,j-1,i),5,1) & 886 ) & 887 ) * ddy & 888 + ( w(k,j,i) * ( ibit6 + ibit7 + ibit8 ) & 889 - w(k-1,j,i) * ( IBITS(wall_flags_0(k-1,j,i),6,1) & 890 + IBITS(wall_flags_0(k-1,j,i),7,1) & 891 + IBITS(wall_flags_0(k-1,j,i),8,1) & 892 ) & 893 ) * ddzw(k) 894 879 895 880 896 tend(k,j,i) = tend(k,j,i) - ( & … … 1247 1263 DO k = nzb+1, nzb_max 1248 1264 1249 ibit14 = IBITS(wall_flags_0(k,j ,i),14,1)1250 ibit13 = IBITS(wall_flags_0(k,j ,i),13,1)1251 ibit12 = IBITS(wall_flags_0(k,j ,i),12,1)1265 ibit14 = IBITS(wall_flags_0(k,j-1,i),14,1) 1266 ibit13 = IBITS(wall_flags_0(k,j-1,i),13,1) 1267 ibit12 = IBITS(wall_flags_0(k,j-1,i),12,1) 1252 1268 1253 1269 v_comp = v(k,j,i) + v(k,j,i-1) - gv … … 1305 1321 DO k = nzb+1, nzb_max 1306 1322 1307 ibit11 = IBITS(wall_flags_0(k,j,i ),11,1)1308 ibit10 = IBITS(wall_flags_0(k,j,i ),10,1)1309 ibit9 = IBITS(wall_flags_0(k,j,i ),9,1)1323 ibit11 = IBITS(wall_flags_0(k,j,i-1),11,1) 1324 ibit10 = IBITS(wall_flags_0(k,j,i-1),10,1) 1325 ibit9 = IBITS(wall_flags_0(k,j,i-1),9,1) 1310 1326 1311 1327 u_comp_l = u(k,j,i) + u(k,j,i-1) - gu … … 1481 1497 !-- correction is needed to overcome numerical instabilities introduced 1482 1498 !-- by a not sufficient reduction of divergences near topography. 1483 div = ( ( u_comp(k) - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx & 1484 + ( v_comp + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy & 1485 + ( w_comp - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) * ddzw(k) & 1499 div = ( ( u_comp(k) * ( ibit9 + ibit10 + ibit11 ) & 1500 - ( u(k,j,i) + u(k,j,i-1) ) & 1501 * ( IBITS(wall_flags_0(k,j,i-1),9,1) & 1502 + IBITS(wall_flags_0(k,j,i-1),10,1) & 1503 + IBITS(wall_flags_0(k,j,i-1),11,1) & 1504 ) & 1505 ) * ddx & 1506 + ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 ) & 1507 - ( v(k,j,i) + v(k,j,i-1 ) ) & 1508 * ( IBITS(wall_flags_0(k,j-1,i),12,1) & 1509 + IBITS(wall_flags_0(k,j-1,i),13,1) & 1510 + IBITS(wall_flags_0(k,j-1,i),14,1) & 1511 ) & 1512 ) * ddy & 1513 + ( w_comp * ( ibit15 + ibit16 + ibit17 ) & 1514 - ( w(k-1,j,i) + w(k-1,j,i-1) ) & 1515 * ( IBITS(wall_flags_0(k-1,j,i),15,1) & 1516 + IBITS(wall_flags_0(k-1,j,i),16,1) & 1517 + IBITS(wall_flags_0(k-1,j,i),17,1) & 1518 ) & 1519 ) * ddzw(k) & 1486 1520 ) * 0.5_wp 1521 1487 1522 1488 1523 tend(k,j,i) = tend(k,j,i) - ( & … … 1702 1737 DO k = nzb+1, nzb_max 1703 1738 1704 ibit20 = IBITS(wall_flags_0(k,j,i ),20,1)1705 ibit19 = IBITS(wall_flags_0(k,j,i ),19,1)1706 ibit18 = IBITS(wall_flags_0(k,j,i ),18,1)1739 ibit20 = IBITS(wall_flags_0(k,j,i-1),20,1) 1740 ibit19 = IBITS(wall_flags_0(k,j,i-1),19,1) 1741 ibit18 = IBITS(wall_flags_0(k,j,i-1),18,1) 1707 1742 1708 1743 u_comp = u(k,j-1,i) + u(k,j,i) - gu … … 1760 1795 DO k = nzb+1, nzb_max 1761 1796 1762 ibit23 = IBITS(wall_flags_0(k,j ,i),23,1)1763 ibit22 = IBITS(wall_flags_0(k,j ,i),22,1)1764 ibit21 = IBITS(wall_flags_0(k,j ,i),21,1)1797 ibit23 = IBITS(wall_flags_0(k,j-1,i),23,1) 1798 ibit22 = IBITS(wall_flags_0(k,j-1,i),22,1) 1799 ibit21 = IBITS(wall_flags_0(k,j-1,i),21,1) 1765 1800 1766 1801 v_comp_l = v(k,j,i) + v(k,j-1,i) - gv … … 1937 1972 !-- correction is needed to overcome numerical instabilities introduced 1938 1973 !-- by a not sufficient reduction of divergences near topography. 1939 div = ( ( u_comp + gu - ( u(k,j-1,i) + u(k,j,i) ) ) * ddx & 1940 + ( v_comp(k) - ( v(k,j,i) + v(k,j-1,i) ) ) * ddy & 1941 + ( w_comp - ( w(k-1,j-1,i) + w(k-1,j,i) ) ) * ddzw(k) & 1974 div = ( ( ( u_comp + gu ) & 1975 * ( ibit18 + ibit19 + ibit20 ) & 1976 - ( u(k,j-1,i) + u(k,j,i) ) & 1977 * ( IBITS(wall_flags_0(k,j,i-1),18,1) & 1978 + IBITS(wall_flags_0(k,j,i-1),19,1) & 1979 + IBITS(wall_flags_0(k,j,i-1),20,1) & 1980 ) & 1981 ) * ddx & 1982 + ( v_comp(k) & 1983 * ( ibit21 + ibit22 + ibit23 ) & 1984 - ( v(k,j,i) + v(k,j-1,i) ) & 1985 * ( IBITS(wall_flags_0(k,j-1,i),21,1) & 1986 + IBITS(wall_flags_0(k,j-1,i),22,1) & 1987 + IBITS(wall_flags_0(k,j-1,i),23,1) & 1988 ) & 1989 ) * ddy & 1990 + ( w_comp & 1991 * ( ibit24 + ibit25 + ibit26 ) & 1992 - ( w(k-1,j-1,i) + w(k-1,j,i) ) & 1993 * ( IBITS(wall_flags_0(k-1,j,i),24,1) & 1994 + IBITS(wall_flags_0(k-1,j,i),25,1) & 1995 + IBITS(wall_flags_0(k-1,j,i),26,1) & 1996 ) & 1997 ) * ddzw(k) & 1942 1998 ) * 0.5_wp 1999 1943 2000 1944 2001 tend(k,j,i) = tend(k,j,i) - ( & … … 2162 2219 2163 2220 DO k = nzb+1, nzb_max 2164 ibit32 = IBITS(wall_flags_00(k,j ,i),0,1)2165 ibit31 = IBITS(wall_flags_0(k,j ,i),31,1)2166 ibit30 = IBITS(wall_flags_0(k,j ,i),30,1)2221 ibit32 = IBITS(wall_flags_00(k,j-1,i),0,1) 2222 ibit31 = IBITS(wall_flags_0(k,j-1,i),31,1) 2223 ibit30 = IBITS(wall_flags_0(k,j-1,i),30,1) 2167 2224 2168 2225 v_comp = v(k+1,j,i) + v(k,j,i) - gv … … 2220 2277 DO k = nzb+1, nzb_max 2221 2278 2222 ibit29 = IBITS(wall_flags_0(k,j,i ),29,1)2223 ibit28 = IBITS(wall_flags_0(k,j,i ),28,1)2224 ibit27 = IBITS(wall_flags_0(k,j,i ),27,1)2279 ibit29 = IBITS(wall_flags_0(k,j,i-1),29,1) 2280 ibit28 = IBITS(wall_flags_0(k,j,i-1),28,1) 2281 ibit27 = IBITS(wall_flags_0(k,j,i-1),27,1) 2225 2282 2226 2283 u_comp = u(k+1,j,i) + u(k,j,i) - gu … … 2402 2459 !-- correction is needed to overcome numerical instabilities introduced 2403 2460 !-- by a not sufficient reduction of divergences near topography. 2404 div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i) ) ) * ddx & 2405 + ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i) ) ) * ddy & 2406 + ( w_comp - ( w(k,j,i) + w(k-1,j,i) ) ) * ddzu(k+1) & 2461 div = ( ( ( u_comp + gu ) * ( ibit27 + ibit28 + ibit29 ) & 2462 - ( u(k+1,j,i) + u(k,j,i) ) & 2463 * ( IBITS(wall_flags_0(k,j,i-1),27,1) & 2464 + IBITS(wall_flags_0(k,j,i-1),28,1) & 2465 + IBITS(wall_flags_0(k,j,i-1),29,1) & 2466 ) & 2467 ) * ddx & 2468 + ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 ) & 2469 - ( v(k+1,j,i) + v(k,j,i) ) & 2470 * ( IBITS(wall_flags_0(k,j-1,i),30,1) & 2471 + IBITS(wall_flags_0(k,j-1,i),31,1) & 2472 + IBITS(wall_flags_00(k,j-1,i),0,1) & 2473 ) & 2474 ) * ddy & 2475 + ( w_comp * ( ibit33 + ibit34 + ibit35 ) & 2476 - ( w(k,j,i) + w(k-1,j,i) ) & 2477 * ( IBITS(wall_flags_00(k-1,j,i),1,1) & 2478 + IBITS(wall_flags_00(k-1,j,i),2,1) & 2479 + IBITS(wall_flags_00(k-1,j,i),3,1) & 2480 ) & 2481 ) * ddzu(k+1) & 2407 2482 ) * 0.5_wp 2483 2408 2484 2409 2485 tend(k,j,i) = tend(k,j,i) - ( & … … 2631 2707 DO k = nzb+1, nzb_max 2632 2708 2633 ibit2 = IBITS(wall_flags_0(k,j,i ),2,1)2634 ibit1 = IBITS(wall_flags_0(k,j,i ),1,1)2635 ibit0 = IBITS(wall_flags_0(k,j,i ),0,1)2709 ibit2 = IBITS(wall_flags_0(k,j,i-1),2,1) 2710 ibit1 = IBITS(wall_flags_0(k,j,i-1),1,1) 2711 ibit0 = IBITS(wall_flags_0(k,j,i-1),0,1) 2636 2712 2637 2713 u_comp = u(k,j,i) - u_gtrans … … 2692 2768 DO k = nzb+1, nzb_max 2693 2769 2694 ibit5 = IBITS(wall_flags_0(k,j ,i),5,1)2695 ibit4 = IBITS(wall_flags_0(k,j ,i),4,1)2696 ibit3 = IBITS(wall_flags_0(k,j ,i),3,1)2770 ibit5 = IBITS(wall_flags_0(k,j-1,i),5,1) 2771 ibit4 = IBITS(wall_flags_0(k,j-1,i),4,1) 2772 ibit3 = IBITS(wall_flags_0(k,j-1,i),3,1) 2697 2773 2698 2774 v_comp = v(k,j,i) - v_gtrans … … 3006 3082 !-- correction is needed to overcome numerical instabilities caused 3007 3083 !-- by a not sufficient reduction of divergences near topography. 3008 div = ( u(k,j,i+1) - u(k,j,i) ) * ddx & 3009 + ( v(k,j+1,i) - v(k,j,i) ) * ddy & 3010 + ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 3084 div = ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 ) & 3085 - u(k,j,i) * ( IBITS(wall_flags_0(k,j,i-1),0,1) & 3086 + IBITS(wall_flags_0(k,j,i-1),1,1) & 3087 + IBITS(wall_flags_0(k,j,i-1),2,1) & 3088 ) & 3089 ) * ddx & 3090 + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 ) & 3091 - v(k,j,i) * ( IBITS(wall_flags_0(k,j-1,i),3,1) & 3092 + IBITS(wall_flags_0(k,j-1,i),4,1) & 3093 + IBITS(wall_flags_0(k,j-1,i),5,1) & 3094 ) & 3095 ) * ddy & 3096 + ( w(k,j,i) * ( ibit6 + ibit7 + ibit8 ) & 3097 - w(k-1,j,i) * ( IBITS(wall_flags_0(k-1,j,i),6,1) & 3098 + IBITS(wall_flags_0(k-1,j,i),7,1) & 3099 + IBITS(wall_flags_0(k-1,j,i),8,1) & 3100 ) & 3101 ) * ddzw(k) 3102 3011 3103 3012 3104 tend(k,j,i) = tend(k,j,i) - ( & … … 3388 3480 DO k = nzb+1, nzt 3389 3481 3390 ibit2 = IBITS(wall_flags_0(k,j,i ),2,1)3391 ibit1 = IBITS(wall_flags_0(k,j,i ),1,1)3392 ibit0 = IBITS(wall_flags_0(k,j,i ),0,1)3482 ibit2 = IBITS(wall_flags_0(k,j,i-1),2,1) 3483 ibit1 = IBITS(wall_flags_0(k,j,i-1),1,1) 3484 ibit0 = IBITS(wall_flags_0(k,j,i-1),0,1) 3393 3485 3394 3486 u_comp = u(k,j,i) - u_gtrans … … 3423 3515 ) 3424 3516 3517 ibit2 = IBITS(wall_flags_0(k,j,i),2,1) 3518 ibit1 = IBITS(wall_flags_0(k,j,i),1,1) 3519 ibit0 = IBITS(wall_flags_0(k,j,i),0,1) 3520 3425 3521 u_comp = u(k,j,i+1) - u_gtrans 3426 3522 flux_r = u_comp * ( & … … 3454 3550 ) 3455 3551 3456 ibit5 = IBITS(wall_flags_0(k,j ,i),5,1)3457 ibit4 = IBITS(wall_flags_0(k,j ,i),4,1)3458 ibit3 = IBITS(wall_flags_0(k,j ,i),3,1)3552 ibit5 = IBITS(wall_flags_0(k,j-1,i),5,1) 3553 ibit4 = IBITS(wall_flags_0(k,j-1,i),4,1) 3554 ibit3 = IBITS(wall_flags_0(k,j-1,i),3,1) 3459 3555 3460 3556 v_comp = v(k,j,i) - v_gtrans … … 3489 3585 ) 3490 3586 3587 ibit5 = IBITS(wall_flags_0(k,j,i),5,1) 3588 ibit4 = IBITS(wall_flags_0(k,j,i),4,1) 3589 ibit3 = IBITS(wall_flags_0(k,j,i),3,1) 3491 3590 3492 3591 v_comp = v(k,j+1,i) - v_gtrans … … 3732 3831 !-- correction is needed to overcome numerical instabilities caused 3733 3832 !-- by a not sufficient reduction of divergences near topography. 3734 div = ( u(k,j,i+1) - u(k,j,i) ) * ddx & 3735 + ( v(k,j+1,i) - v(k,j,i) ) * ddy & 3736 + ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 3833 div = ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 ) & 3834 - u(k,j,i) * ( IBITS(wall_flags_0(k,j,i-1),0,1) & 3835 + IBITS(wall_flags_0(k,j,i-1),1,1) & 3836 + IBITS(wall_flags_0(k,j,i-1),2,1) & 3837 ) & 3838 ) * ddx & 3839 + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 ) & 3840 - v(k,j,i) * ( IBITS(wall_flags_0(k,j-1,i),3,1) & 3841 + IBITS(wall_flags_0(k,j-1,i),4,1) & 3842 + IBITS(wall_flags_0(k,j-1,i),5,1) & 3843 ) & 3844 ) * ddy & 3845 + ( w(k,j,i) * ( ibit6 + ibit7 + ibit8 ) & 3846 - w(k-1,j,i) * ( IBITS(wall_flags_0(k-1,j,i),6,1) & 3847 + IBITS(wall_flags_0(k-1,j,i),7,1) & 3848 + IBITS(wall_flags_0(k-1,j,i),8,1) & 3849 ) & 3850 ) * ddzw(k) 3851 3737 3852 3738 3853 tend(k,j,i) = - ( & … … 3852 3967 DO k = nzb+1, nzb_max 3853 3968 3854 ibit11 = IBITS(wall_flags_0(k,j,i ),11,1)3855 ibit10 = IBITS(wall_flags_0(k,j,i ),10,1)3856 ibit9 = IBITS(wall_flags_0(k,j,i ),9,1)3969 ibit11 = IBITS(wall_flags_0(k,j,i-1),11,1) 3970 ibit10 = IBITS(wall_flags_0(k,j,i-1),10,1) 3971 ibit9 = IBITS(wall_flags_0(k,j,i-1),9,1) 3857 3972 3858 3973 u_comp(k) = u(k,j,i) + u(k,j,i-1) - gu … … 3910 4025 DO k = nzb+1, nzb_max 3911 4026 3912 ibit14 = IBITS(wall_flags_0(k,j ,i),14,1)3913 ibit13 = IBITS(wall_flags_0(k,j ,i),13,1)3914 ibit12 = IBITS(wall_flags_0(k,j ,i),12,1)4027 ibit14 = IBITS(wall_flags_0(k,j-1,i),14,1) 4028 ibit13 = IBITS(wall_flags_0(k,j-1,i),13,1) 4029 ibit12 = IBITS(wall_flags_0(k,j-1,i),12,1) 3915 4030 3916 4031 v_comp = v(k,j,i) + v(k,j,i-1) - gv … … 4085 4200 !-- correction is needed to overcome numerical instabilities caused 4086 4201 !-- by a not sufficient reduction of divergences near topography. 4087 div = ( ( u_comp(k) - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx & 4088 + ( v_comp + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy & 4089 + ( w_comp - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) & 4090 * ddzw(k) & 4091 ) * 0.5_wp 4202 div = ( ( u_comp(k) * ( ibit9 + ibit10 + ibit11 ) & 4203 - ( u(k,j,i) + u(k,j,i-1) ) & 4204 * ( IBITS(wall_flags_0(k,j,i-1),9,1) & 4205 + IBITS(wall_flags_0(k,j,i-1),10,1) & 4206 + IBITS(wall_flags_0(k,j,i-1),11,1) & 4207 ) & 4208 ) * ddx & 4209 + ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 ) & 4210 - ( v(k,j,i) + v(k,j,i-1 ) ) & 4211 * ( IBITS(wall_flags_0(k,j-1,i),12,1) & 4212 + IBITS(wall_flags_0(k,j-1,i),13,1) & 4213 + IBITS(wall_flags_0(k,j-1,i),14,1) & 4214 ) & 4215 ) * ddy & 4216 + ( w_comp * ( ibit15 + ibit16 + ibit17 ) & 4217 - ( w(k-1,j,i) + w(k-1,j,i-1) ) & 4218 * ( IBITS(wall_flags_0(k-1,j,i),15,1) & 4219 + IBITS(wall_flags_0(k-1,j,i),16,1) & 4220 + IBITS(wall_flags_0(k-1,j,i),17,1) & 4221 ) & 4222 ) * ddzw(k) & 4223 ) * 0.5_wp 4224 4225 4092 4226 4093 4227 tend(k,j,i) = tend(k,j,i) - ( & … … 4316 4450 DO k = nzb+1, nzt 4317 4451 4318 ibit11 = IBITS(wall_flags_0(k,j,i ),11,1)4319 ibit10 = IBITS(wall_flags_0(k,j,i ),10,1)4320 ibit9 = IBITS(wall_flags_0(k,j,i ),9,1)4452 ibit11 = IBITS(wall_flags_0(k,j,i-1),11,1) 4453 ibit10 = IBITS(wall_flags_0(k,j,i-1),10,1) 4454 ibit9 = IBITS(wall_flags_0(k,j,i-1),9,1) 4321 4455 4322 4456 u_comp_l = u(k,j,i) + u(k,j,i-1) - gu … … 4351 4485 ) 4352 4486 4487 ibit11 = IBITS(wall_flags_0(k,j,i),11,1) 4488 ibit10 = IBITS(wall_flags_0(k,j,i),10,1) 4489 ibit9 = IBITS(wall_flags_0(k,j,i),9,1) 4490 4353 4491 u_comp = u(k,j,i+1) + u(k,j,i) 4354 4492 flux_r = ( u_comp - gu ) * ( & … … 4382 4520 ) 4383 4521 4384 ibit14 = IBITS(wall_flags_0(k,j ,i),14,1)4385 ibit13 = IBITS(wall_flags_0(k,j ,i),13,1)4386 ibit12 = IBITS(wall_flags_0(k,j ,i),12,1)4522 ibit14 = IBITS(wall_flags_0(k,j-1,i),14,1) 4523 ibit13 = IBITS(wall_flags_0(k,j-1,i),13,1) 4524 ibit12 = IBITS(wall_flags_0(k,j-1,i),12,1) 4387 4525 4388 4526 v_comp_s = v(k,j,i) + v(k,j,i-1) - gv … … 4418 4556 4419 4557 4558 ibit14 = IBITS(wall_flags_0(k,j,i),14,1) 4559 ibit13 = IBITS(wall_flags_0(k,j,i),13,1) 4560 ibit12 = IBITS(wall_flags_0(k,j,i),12,1) 4561 4420 4562 v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv 4421 4563 flux_n = v_comp * ( & … … 4532 4674 !-- correction is needed to overcome numerical instabilities caused 4533 4675 !-- by a not sufficient reduction of divergences near topography. 4534 div = ( ( u_comp - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx & 4535 + ( v_comp + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy & 4536 + ( w_comp - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) & 4537 * ddzw(k) & 4538 ) * 0.5_wp 4676 div = ( ( u_comp * ( ibit9 + ibit10 + ibit11 ) & 4677 - ( u(k,j,i) + u(k,j,i-1) ) & 4678 * ( IBITS(wall_flags_0(k,j,i-1),9,1) & 4679 + IBITS(wall_flags_0(k,j,i-1),10,1) & 4680 + IBITS(wall_flags_0(k,j,i-1),11,1) & 4681 ) & 4682 ) * ddx & 4683 + ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 ) & 4684 - ( v(k,j,i) + v(k,j,i-1 ) ) & 4685 * ( IBITS(wall_flags_0(k,j-1,i),12,1) & 4686 + IBITS(wall_flags_0(k,j-1,i),13,1) & 4687 + IBITS(wall_flags_0(k,j-1,i),14,1) & 4688 ) & 4689 ) * ddy & 4690 + ( w_comp * ( ibit15 + ibit16 + ibit17 ) & 4691 - ( w(k-1,j,i) + w(k-1,j,i-1) ) & 4692 * ( IBITS(wall_flags_0(k-1,j,i),15,1) & 4693 + IBITS(wall_flags_0(k-1,j,i),16,1) & 4694 + IBITS(wall_flags_0(k-1,j,i),17,1) & 4695 ) & 4696 ) * ddzw(k) & 4697 ) * 0.5_wp 4698 4539 4699 4540 4700 tend(k,j,i) = - ( & … … 4646 4806 DO k = nzb+1, nzb_max 4647 4807 4648 ibit20 = IBITS(wall_flags_0(k,j,i ),20,1)4649 ibit19 = IBITS(wall_flags_0(k,j,i ),19,1)4650 ibit18 = IBITS(wall_flags_0(k,j,i ),18,1)4808 ibit20 = IBITS(wall_flags_0(k,j,i-1),20,1) 4809 ibit19 = IBITS(wall_flags_0(k,j,i-1),19,1) 4810 ibit18 = IBITS(wall_flags_0(k,j,i-1),18,1) 4651 4811 4652 4812 u_comp = u(k,j-1,i) + u(k,j,i) - gu … … 4704 4864 DO k = nzb+1, nzb_max 4705 4865 4706 ibit23 = IBITS(wall_flags_0(k,j ,i),23,1)4707 ibit22 = IBITS(wall_flags_0(k,j ,i),22,1)4708 ibit21 = IBITS(wall_flags_0(k,j ,i),21,1)4866 ibit23 = IBITS(wall_flags_0(k,j-1,i),23,1) 4867 ibit22 = IBITS(wall_flags_0(k,j-1,i),22,1) 4868 ibit21 = IBITS(wall_flags_0(k,j-1,i),21,1) 4709 4869 4710 4870 v_comp(k) = v(k,j,i) + v(k,j-1,i) - gv … … 4878 5038 !-- correction is needed to overcome numerical instabilities caused 4879 5039 !-- by a not sufficient reduction of divergences near topography. 4880 div = ( ( u_comp + gu - ( u(k,j-1,i) + u(k,j,i) ) ) * ddx & 4881 + ( v_comp(k) - ( v(k,j,i) + v(k,j-1,i) ) ) * ddy & 4882 + ( w_comp - ( w(k-1,j-1,i) + w(k-1,j,i) ) & 4883 ) * ddzw(k) & 4884 ) * 0.5_wp 5040 div = ( ( ( u_comp + gu ) & 5041 * ( ibit18 + ibit19 + ibit20 ) & 5042 - ( u(k,j-1,i) + u(k,j,i) ) & 5043 * ( IBITS(wall_flags_0(k,j,i-1),18,1) & 5044 + IBITS(wall_flags_0(k,j,i-1),19,1) & 5045 + IBITS(wall_flags_0(k,j,i-1),20,1) & 5046 ) & 5047 ) * ddx & 5048 + ( v_comp(k) & 5049 * ( ibit21 + ibit22 + ibit23 ) & 5050 - ( v(k,j,i) + v(k,j-1,i) ) & 5051 * ( IBITS(wall_flags_0(k,j-1,i),21,1) & 5052 + IBITS(wall_flags_0(k,j-1,i),22,1) & 5053 + IBITS(wall_flags_0(k,j-1,i),23,1) & 5054 ) & 5055 ) * ddy & 5056 + ( w_comp & 5057 * ( ibit24 + ibit25 + ibit26 ) & 5058 - ( w(k-1,j-1,i) + w(k-1,j,i) ) & 5059 * ( IBITS(wall_flags_0(k-1,j,i),24,1) & 5060 + IBITS(wall_flags_0(k-1,j,i),25,1) & 5061 + IBITS(wall_flags_0(k-1,j,i),26,1) & 5062 ) & 5063 ) * ddzw(k) & 5064 ) * 0.5_wp 5065 4885 5066 4886 5067 tend(k,j,i) = tend(k,j,i) - ( & … … 5120 5301 DO k = nzb+1, nzt 5121 5302 5122 ibit20 = IBITS(wall_flags_0(k,j,i ),20,1)5123 ibit19 = IBITS(wall_flags_0(k,j,i ),19,1)5124 ibit18 = IBITS(wall_flags_0(k,j,i ),18,1)5303 ibit20 = IBITS(wall_flags_0(k,j,i-1),20,1) 5304 ibit19 = IBITS(wall_flags_0(k,j,i-1),19,1) 5305 ibit18 = IBITS(wall_flags_0(k,j,i-1),18,1) 5125 5306 5126 5307 u_comp_l = u(k,j-1,i) + u(k,j,i) - gu … … 5155 5336 ) 5156 5337 5338 ibit20 = IBITS(wall_flags_0(k,j,i),20,1) 5339 ibit19 = IBITS(wall_flags_0(k,j,i),19,1) 5340 ibit18 = IBITS(wall_flags_0(k,j,i),18,1) 5341 5157 5342 u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu 5158 5343 flux_r = u_comp * ( & … … 5186 5371 ) 5187 5372 5188 ibit23 = IBITS(wall_flags_0(k,j ,i),23,1)5189 ibit22 = IBITS(wall_flags_0(k,j ,i),22,1)5190 ibit21 = IBITS(wall_flags_0(k,j ,i),21,1)5373 ibit23 = IBITS(wall_flags_0(k,j-1,i),23,1) 5374 ibit22 = IBITS(wall_flags_0(k,j-1,i),22,1) 5375 ibit21 = IBITS(wall_flags_0(k,j-1,i),21,1) 5191 5376 5192 5377 … … 5222 5407 ) 5223 5408 5409 ibit23 = IBITS(wall_flags_0(k,j,i),23,1) 5410 ibit22 = IBITS(wall_flags_0(k,j,i),22,1) 5411 ibit21 = IBITS(wall_flags_0(k,j,i),21,1) 5412 5224 5413 v_comp = v(k,j+1,i) + v(k,j,i) 5225 5414 flux_n = ( v_comp - gv ) * ( & … … 5336 5525 !-- correction is needed to overcome numerical instabilities caused 5337 5526 !-- by a not sufficient reduction of divergences near topography. 5338 div = ( ( u_comp + gu - ( u(k,j-1,i) + u(k,j,i) ) ) * ddx & 5339 + ( v_comp - ( v(k,j,i) + v(k,j-1,i) ) ) * ddy & 5340 + ( w_comp - ( w(k-1,j-1,i) + w(k-1,j,i) ) & 5341 ) * ddzw(k) & 5342 ) * 0.5_wp 5527 div = ( ( ( u_comp + gu ) & 5528 * ( ibit18 + ibit19 + ibit20 ) & 5529 - ( u(k,j-1,i) + u(k,j,i) ) & 5530 * ( IBITS(wall_flags_0(k,j,i-1),18,1) & 5531 + IBITS(wall_flags_0(k,j,i-1),19,1) & 5532 + IBITS(wall_flags_0(k,j,i-1),20,1) & 5533 ) & 5534 ) * ddx & 5535 + ( v_comp & 5536 * ( ibit21 + ibit22 + ibit23 ) & 5537 - ( v(k,j,i) + v(k,j-1,i) ) & 5538 * ( IBITS(wall_flags_0(k,j-1,i),21,1) & 5539 + IBITS(wall_flags_0(k,j-1,i),22,1) & 5540 + IBITS(wall_flags_0(k,j-1,i),23,1) & 5541 ) & 5542 ) * ddy & 5543 + ( w_comp & 5544 * ( ibit24 + ibit25 + ibit26 ) & 5545 - ( w(k-1,j-1,i) + w(k-1,j,i) ) & 5546 * ( IBITS(wall_flags_0(k-1,j,i),24,1) & 5547 + IBITS(wall_flags_0(k-1,j,i),25,1) & 5548 + IBITS(wall_flags_0(k-1,j,i),26,1) & 5549 ) & 5550 ) * ddzw(k) & 5551 ) * 0.5_wp 5552 5343 5553 5344 5554 tend(k,j,i) = - ( & … … 5452 5662 DO k = nzb+1, nzb_max 5453 5663 5454 ibit29 = IBITS(wall_flags_0(k,j,i ),29,1)5455 ibit28 = IBITS(wall_flags_0(k,j,i ),28,1)5456 ibit27 = IBITS(wall_flags_0(k,j,i ),27,1)5664 ibit29 = IBITS(wall_flags_0(k,j,i-1),29,1) 5665 ibit28 = IBITS(wall_flags_0(k,j,i-1),28,1) 5666 ibit27 = IBITS(wall_flags_0(k,j,i-1),27,1) 5457 5667 5458 5668 u_comp = u(k+1,j,i) + u(k,j,i) - gu … … 5510 5720 DO k = nzb+1, nzb_max 5511 5721 5512 ibit32 = IBITS(wall_flags_00(k,j ,i),0,1)5513 ibit31 = IBITS(wall_flags_0(k,j ,i),31,1)5514 ibit30 = IBITS(wall_flags_0(k,j ,i),30,1)5722 ibit32 = IBITS(wall_flags_00(k,j-1,i),0,1) 5723 ibit31 = IBITS(wall_flags_0(k,j-1,i),31,1) 5724 ibit30 = IBITS(wall_flags_0(k,j-1,i),30,1) 5515 5725 5516 5726 v_comp = v(k+1,j,i) + v(k,j,i) - gv … … 5690 5900 !-- correction is needed to overcome numerical instabilities caused 5691 5901 !-- by a not sufficient reduction of divergences near topography. 5692 div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i) ) ) * ddx & 5693 + ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i) ) ) * ddy & 5694 + ( w_comp - ( w(k,j,i) + w(k-1,j,i) ) ) & 5695 * ddzu(k+1) & 5696 ) * 0.5_wp 5902 div = ( ( ( u_comp + gu ) * ( ibit27 + ibit28 + ibit29 ) & 5903 - ( u(k+1,j,i) + u(k,j,i) ) & 5904 * ( IBITS(wall_flags_0(k,j,i-1),27,1) & 5905 + IBITS(wall_flags_0(k,j,i-1),28,1) & 5906 + IBITS(wall_flags_0(k,j,i-1),29,1) & 5907 ) & 5908 ) * ddx & 5909 + ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 ) & 5910 - ( v(k+1,j,i) + v(k,j,i) ) & 5911 * ( IBITS(wall_flags_0(k,j-1,i),30,1) & 5912 + IBITS(wall_flags_0(k,j-1,i),31,1) & 5913 + IBITS(wall_flags_00(k,j-1,i),0,1) & 5914 ) & 5915 ) * ddy & 5916 + ( w_comp * ( ibit33 + ibit34 + ibit35 ) & 5917 - ( w(k,j,i) + w(k-1,j,i) ) & 5918 * ( IBITS(wall_flags_00(k-1,j,i),1,1) & 5919 + IBITS(wall_flags_00(k-1,j,i),2,1) & 5920 + IBITS(wall_flags_00(k-1,j,i),3,1) & 5921 ) & 5922 ) * ddzu(k+1) & 5923 ) * 0.5_wp 5924 5925 5697 5926 5698 5927 tend(k,j,i) = tend(k,j,i) - ( & … … 5903 6132 DO k = nzb+1, nzt 5904 6133 5905 ibit27 = IBITS(wall_flags_0(k,j,i ),27,1)5906 ibit28 = IBITS(wall_flags_0(k,j,i ),28,1)5907 ibit29 = IBITS(wall_flags_0(k,j,i ),29,1)6134 ibit27 = IBITS(wall_flags_0(k,j,i-1),27,1) 6135 ibit28 = IBITS(wall_flags_0(k,j,i-1),28,1) 6136 ibit29 = IBITS(wall_flags_0(k,j,i-1),29,1) 5908 6137 5909 6138 u_comp_l = u(k+1,j,i) + u(k,j,i) - gu … … 5938 6167 ) 5939 6168 6169 ibit27 = IBITS(wall_flags_0(k,j,i),27,1) 6170 ibit28 = IBITS(wall_flags_0(k,j,i),28,1) 6171 ibit29 = IBITS(wall_flags_0(k,j,i),29,1) 6172 5940 6173 u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu 5941 6174 flux_r = u_comp * ( & … … 5968 6201 ( w(k,j,i+3) - w(k,j,i-2) ) & 5969 6202 ) 5970 ibit32 = IBITS(wall_flags_00(k,j ,i),0,1)5971 ibit31 = IBITS(wall_flags_0(k,j ,i),31,1)5972 ibit30 = IBITS(wall_flags_0(k,j ,i),30,1)6203 ibit32 = IBITS(wall_flags_00(k,j-1,i),0,1) 6204 ibit31 = IBITS(wall_flags_0(k,j-1,i),31,1) 6205 ibit30 = IBITS(wall_flags_0(k,j-1,i),30,1) 5973 6206 5974 6207 v_comp_s = v(k+1,j,i) + v(k,j,i) - gv … … 6003 6236 ) 6004 6237 6238 ibit32 = IBITS(wall_flags_00(k,j,i),0,1) 6239 ibit31 = IBITS(wall_flags_0(k,j,i),31,1) 6240 ibit30 = IBITS(wall_flags_0(k,j,i),30,1) 6241 6005 6242 v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv 6006 6243 flux_n = v_comp * ( & … … 6118 6355 !-- correction is needed to overcome numerical instabilities caused 6119 6356 !-- by a not sufficient reduction of divergences near topography. 6120 div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i) ) ) * ddx & 6121 + ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i) ) ) * ddy & 6122 + ( w_comp - ( w(k,j,i) + w(k-1,j,i) ) ) & 6123 * ddzu(k+1) & 6124 ) * 0.5_wp 6357 div = ( ( ( u_comp + gu ) * ( ibit27 + ibit28 + ibit29 ) & 6358 - ( u(k+1,j,i) + u(k,j,i) ) & 6359 * ( IBITS(wall_flags_0(k,j,i-1),27,1) & 6360 + IBITS(wall_flags_0(k,j,i-1),28,1) & 6361 + IBITS(wall_flags_0(k,j,i-1),29,1) & 6362 ) & 6363 ) * ddx & 6364 + ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 ) & 6365 - ( v(k+1,j,i) + v(k,j,i) ) & 6366 * ( IBITS(wall_flags_0(k,j-1,i),30,1) & 6367 + IBITS(wall_flags_0(k,j-1,i),31,1) & 6368 + IBITS(wall_flags_00(k,j-1,i),0,1) & 6369 ) & 6370 ) * ddy & 6371 + ( w_comp * ( ibit33 + ibit34 + ibit35 ) & 6372 - ( w(k,j,i) + w(k-1,j,i) ) & 6373 * ( IBITS(wall_flags_00(k-1,j,i),1,1) & 6374 + IBITS(wall_flags_00(k-1,j,i),2,1) & 6375 + IBITS(wall_flags_00(k-1,j,i),3,1) & 6376 ) & 6377 ) * ddzu(k+1) & 6378 ) * 0.5_wp 6379 6125 6380 6126 6381 tend(k,j,i) = - ( &
Note: See TracChangeset
for help on using the changeset viewer.