Changeset 2037
- Timestamp:
- Oct 26, 2016 11:15:40 AM (8 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 22 edited
-
advec_ws.f90 (modified) (110 diffs)
-
check_parameters.f90 (modified) (17 diffs)
-
diffusion_e.f90 (modified) (9 diffs)
-
diffusion_s.f90 (modified) (13 diffs)
-
diffusion_u.f90 (modified) (13 diffs)
-
diffusion_v.f90 (modified) (13 diffs)
-
diffusion_w.f90 (modified) (10 diffs)
-
disturb_heatflux.f90 (modified) (4 diffs)
-
flow_statistics.f90 (modified) (58 diffs)
-
header.f90 (modified) (3 diffs)
-
init_3d_model.f90 (modified) (12 diffs)
-
init_grid.f90 (modified) (4 diffs)
-
ls_forcing_mod.f90 (modified) (3 diffs)
-
modules.f90 (modified) (3 diffs)
-
netcdf_interface_mod.f90 (modified) (3 diffs)
-
parin.f90 (modified) (3 diffs)
-
poismg_mod.f90 (modified) (17 diffs)
-
poismg_noopt.f90 (modified) (16 diffs)
-
pres.f90 (modified) (6 diffs)
-
sor.f90 (modified) (7 diffs)
-
surface_layer_fluxes_mod.f90 (modified) (9 diffs)
-
tridia_solver_mod.f90 (modified) (8 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/advec_ws.f90
r2001 r2037 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Anelastic approximation implemented 23 23 ! 24 24 ! Former revisions: … … 902 902 903 903 USE arrays_3d, & 904 ONLY: ddzw, tend, u, v,w904 ONLY: ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw 905 905 906 906 USE constants, & … … 1211 1211 1212 1212 1213 flux_t(k) = w(k,j,i) * (&1213 flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( & 1214 1214 ( 37.0_wp * ibit8 * adv_sca_5 & 1215 1215 + 7.0_wp * ibit7 * adv_sca_3 & … … 1225 1225 ) 1226 1226 1227 diss_t(k) = -ABS( w(k,j,i) ) * (&1227 diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( & 1228 1228 ( 10.0_wp * ibit8 * adv_sca_5 & 1229 1229 + 3.0_wp * ibit7 * adv_sca_3 & … … 1266 1266 fd_1 = ( w(k-1,j,i) * ( sk(k,j,i) + sk(k-1,j,i) ) & 1267 1267 -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) ) & 1268 ) * adv_sca_1 1268 ) * adv_sca_1 * rho_air_zw(k) 1269 1269 1270 1270 ft_1 = ( w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) & 1271 1271 -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) ) & 1272 ) * adv_sca_1 1272 ) * adv_sca_1 * rho_air_zw(k) 1273 1273 ! 1274 1274 !-- Calculate ratio of upwind gradients. Note, Min/Max is just to … … 1386 1386 + IBITS(wall_flags_0(k,j,i-1),2,1) & 1387 1387 ) & 1388 ) * ddx&1388 ) * rho_air(k) * ddx & 1389 1389 + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 ) & 1390 1390 - v(k,j,i) * ( IBITS(wall_flags_0(k,j-1,i),3,1) & … … 1392 1392 + IBITS(wall_flags_0(k,j-1,i),5,1) & 1393 1393 ) & 1394 ) * ddy & 1395 + ( w(k,j,i) * ( ibit6 + ibit7 + ibit8 ) & 1396 - w(k-1,j,i) * ( IBITS(wall_flags_0(k-1,j,i),6,1) & 1394 ) * rho_air(k) * ddy & 1395 + ( w(k,j,i) * rho_air_zw(k) * & 1396 ( ibit6 + ibit7 + ibit8 ) & 1397 - w(k-1,j,i) * rho_air_zw(k-1) * & 1398 ( IBITS(wall_flags_0(k-1,j,i),6,1) & 1397 1399 + IBITS(wall_flags_0(k-1,j,i),7,1) & 1398 1400 + IBITS(wall_flags_0(k-1,j,i),8,1) & … … 1406 1408 + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn) - & 1407 1409 swap_diss_y_local(k,tn) ) * ddy & 1408 + ( flux_t(k) + diss_t(k) - flux_d - diss_d & 1409 ) * ddzw(k) & 1410 + ( ( flux_t(k) + diss_t(k) ) - & 1411 ( flux_d + diss_d ) & 1412 ) * drho_air(k) * ddzw(k) & 1410 1413 ) + sk(k,j,i) * div 1411 1414 … … 1454 1457 k_mm = k - 2 * ibit8 1455 1458 1456 flux_t(k) = w(k,j,i) * (&1459 flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( & 1457 1460 ( 37.0_wp * ibit8 * adv_sca_5 & 1458 1461 + 7.0_wp * ibit7 * adv_sca_3 & … … 1468 1471 ) 1469 1472 1470 diss_t(k) = -ABS( w(k,j,i) ) * (&1473 diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( & 1471 1474 ( 10.0_wp * ibit8 * adv_sca_5 & 1472 1475 + 3.0_wp * ibit7 * adv_sca_3 & … … 1511 1514 fd_1 = ( w(k-1,j,i) * ( sk(k,j,i) + sk(k-1,j,i) ) & 1512 1515 -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) ) & 1513 ) * adv_sca_1 1516 ) * adv_sca_1 * rho_air_zw(k) 1514 1517 1515 1518 ft_1 = ( w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) & 1516 1519 -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) ) & 1517 ) * adv_sca_1 1520 ) * adv_sca_1 * rho_air_zw(k) 1518 1521 ! 1519 1522 !-- Calculate ratio of upwind gradients. Note, Min/Max is just to … … 1626 1629 !-- correction is needed to overcome numerical instabilities introduced 1627 1630 !-- by a not sufficient reduction of divergences near topography. 1628 div = ( u(k,j,i+1) - u(k,j,i) ) * ddx & 1629 + ( v(k,j+1,i) - v(k,j,i) ) * ddy & 1630 + ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 1631 div = ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx & 1632 + ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy & 1633 + ( w(k,j,i) * rho_air_zw(k) - & 1634 w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k) 1631 1635 1632 1636 tend(k,j,i) = tend(k,j,i) - ( & … … 1635 1639 + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn) - & 1636 1640 swap_diss_y_local(k,tn) ) * ddy & 1637 + ( flux_t(k) + diss_t(k) - flux_d - diss_d & 1638 ) * ddzw(k) & 1641 + ( ( flux_t(k) + diss_t(k) ) - & 1642 ( flux_d + diss_d ) & 1643 ) * drho_air(k) * ddzw(k) & 1639 1644 ) + sk(k,j,i) * div 1640 1645 … … 1734 1739 1735 1740 USE arrays_3d, & 1736 ONLY: ddzw, diss_l_u, diss_s_u, flux_l_u, flux_s_u, tend, u, v, w 1741 ONLY: ddzw, diss_l_u, diss_s_u, flux_l_u, flux_s_u, tend, u, v, w,& 1742 drho_air, rho_air, rho_air_zw 1737 1743 1738 1744 USE constants, & … … 1999 2005 2000 2006 w_comp = w(k,j,i) + w(k,j,i-1) 2001 flux_t(k) = w_comp * (&2007 flux_t(k) = w_comp * rho_air_zw(k) * ( & 2002 2008 ( 37.0_wp * ibit17 * adv_mom_5 & 2003 2009 + 7.0_wp * ibit16 * adv_mom_3 & … … 2014 2020 ) 2015 2021 2016 diss_t(k) = - ABS( w_comp ) * (&2022 diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * ( & 2017 2023 ( 10.0_wp * ibit17 * adv_mom_5 & 2018 2024 + 3.0_wp * ibit16 * adv_mom_3 & … … 2037 2043 + IBITS(wall_flags_0(k,j,i-1),10,1) & 2038 2044 + IBITS(wall_flags_0(k,j,i-1),11,1) & 2039 ) & 2040 ) * ddx&2045 ) & 2046 ) * rho_air(k) * ddx & 2041 2047 + ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 ) & 2042 2048 - ( v(k,j,i) + v(k,j,i-1 ) ) & … … 2045 2051 + IBITS(wall_flags_0(k,j-1,i),14,1) & 2046 2052 ) & 2047 ) * ddy&2048 + ( w_comp * ( ibit15 + ibit16 + ibit17 )&2049 - ( w(k-1,j,i) + w(k-1,j,i-1) ) 2053 ) * rho_air(k) * ddy & 2054 + ( w_comp * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 ) & 2055 - ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1) & 2050 2056 * ( IBITS(wall_flags_0(k-1,j,i),15,1) & 2051 2057 + IBITS(wall_flags_0(k-1,j,i),16,1) & … … 2061 2067 + ( flux_n(k) + diss_n(k) & 2062 2068 - flux_s_u(k,tn) - diss_s_u(k,tn) ) * ddy & 2063 + ( flux_t(k) + diss_t(k) & 2064 - flux_d - diss_d ) * ddzw(k) & 2069 + ( ( flux_t(k) + diss_t(k) ) & 2070 - ( flux_d + diss_d ) & 2071 ) * drho_air(k) * ddzw(k) & 2065 2072 ) + div * u(k,j,i) 2066 2073 … … 2127 2134 2128 2135 w_comp = w(k,j,i) + w(k,j,i-1) 2129 flux_t(k) = w_comp * (&2136 flux_t(k) = w_comp * rho_air_zw(k) * ( & 2130 2137 ( 37.0_wp * ibit17 * adv_mom_5 & 2131 2138 + 7.0_wp * ibit16 * adv_mom_3 & … … 2142 2149 ) 2143 2150 2144 diss_t(k) = - ABS( w_comp ) * (&2151 diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * ( & 2145 2152 ( 10.0_wp * ibit17 * adv_mom_5 & 2146 2153 + 3.0_wp * ibit16 * adv_mom_3 & … … 2161 2168 !-- by a not sufficient reduction of divergences near topography. 2162 2169 div = ( ( u_comp(k) - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx & 2170 * rho_air(k)& 2163 2171 + ( v_comp + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy & 2164 + ( w_comp - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) * ddzw(k) & 2172 * rho_air(k)& 2173 + ( w_comp * rho_air_zw(k) - & 2174 ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1) & 2175 ) * ddzw(k) & 2165 2176 ) * 0.5_wp 2166 2177 … … 2170 2181 + ( flux_n(k) + diss_n(k) & 2171 2182 - flux_s_u(k,tn) - diss_s_u(k,tn) ) * ddy & 2172 + ( flux_t(k) + diss_t(k) & 2173 - flux_d - diss_d ) * ddzw(k) & 2183 + ( ( flux_t(k) + diss_t(k) ) & 2184 - ( flux_d + diss_d ) & 2185 ) * drho_air(k) * ddzw(k) & 2174 2186 ) + div * u(k,j,i) 2175 2187 … … 2219 2231 2220 2232 USE arrays_3d, & 2221 ONLY: ddzw, diss_l_v, diss_s_v, flux_l_v, flux_s_v, tend, u, v, w 2233 ONLY: ddzw, diss_l_v, diss_s_v, flux_l_v, flux_s_v, tend, u, v, w, & 2234 drho_air, rho_air, rho_air_zw 2222 2235 2223 2236 USE constants, & … … 2486 2499 2487 2500 w_comp = w(k,j-1,i) + w(k,j,i) 2488 flux_t(k) = w_comp * (&2501 flux_t(k) = w_comp * rho_air_zw(k) * ( & 2489 2502 ( 37.0_wp * ibit26 * adv_mom_5 & 2490 2503 + 7.0_wp * ibit25 * adv_mom_3 & … … 2501 2514 ) 2502 2515 2503 diss_t(k) = - ABS( w_comp ) * (&2516 diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * ( & 2504 2517 ( 10.0_wp * ibit26 * adv_mom_5 & 2505 2518 + 3.0_wp * ibit25 * adv_mom_3 & … … 2525 2538 + IBITS(wall_flags_0(k,j,i-1),19,1) & 2526 2539 + IBITS(wall_flags_0(k,j,i-1),20,1) & 2527 ) & 2528 ) * ddx&2540 ) & 2541 ) * rho_air(k) * ddx & 2529 2542 + ( v_comp(k) & 2530 2543 * ( ibit21 + ibit22 + ibit23 ) & … … 2533 2546 + IBITS(wall_flags_0(k,j-1,i),22,1) & 2534 2547 + IBITS(wall_flags_0(k,j-1,i),23,1) & 2535 ) & 2536 ) * ddy & 2537 + ( w_comp & 2538 * ( ibit24 + ibit25 + ibit26 ) & 2539 - ( w(k-1,j-1,i) + w(k-1,j,i) ) & 2548 ) & 2549 ) * rho_air(k) * ddy & 2550 + ( w_comp * rho_air_zw(k) * ( ibit24 + ibit25 + ibit26 ) & 2551 - ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1) & 2540 2552 * ( IBITS(wall_flags_0(k-1,j,i),24,1) & 2541 2553 + IBITS(wall_flags_0(k-1,j,i),25,1) & … … 2551 2563 + ( flux_n(k) + diss_n(k) & 2552 2564 - flux_s_v(k,tn) - diss_s_v(k,tn) ) * ddy & 2553 + ( flux_t(k) + diss_t(k) & 2554 - flux_d - diss_d ) * ddzw(k) & 2565 + ( ( flux_t(k) + diss_t(k) ) & 2566 - ( flux_d + diss_d ) & 2567 ) * drho_air(k) * ddzw(k) & 2555 2568 ) + v(k,j,i) * div 2556 2569 … … 2622 2635 2623 2636 w_comp = w(k,j-1,i) + w(k,j,i) 2624 flux_t(k) = w_comp * (&2637 flux_t(k) = w_comp * rho_air_zw(k) * ( & 2625 2638 ( 37.0_wp * ibit26 * adv_mom_5 & 2626 2639 + 7.0_wp * ibit25 * adv_mom_3 & … … 2637 2650 ) 2638 2651 2639 diss_t(k) = - ABS( w_comp ) * (&2652 diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * ( & 2640 2653 ( 10.0_wp * ibit26 * adv_mom_5 & 2641 2654 + 3.0_wp * ibit25 * adv_mom_3 & … … 2656 2669 !-- by a not sufficient reduction of divergences near topography. 2657 2670 div = ( ( u_comp + gu - ( u(k,j-1,i) + u(k,j,i) ) ) * ddx & 2671 * rho_air(k)& 2658 2672 + ( v_comp(k) - ( v(k,j,i) + v(k,j-1,i) ) ) * ddy & 2659 + ( w_comp - ( w(k-1,j-1,i) + w(k-1,j,i) ) ) * ddzw(k) & 2673 * rho_air(k)& 2674 + ( w_comp * rho_air_zw(k) - & 2675 ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1) & 2676 ) * ddzw(k) & 2660 2677 ) * 0.5_wp 2661 2678 … … 2665 2682 + ( flux_n(k) + diss_n(k) & 2666 2683 - flux_s_v(k,tn) - diss_s_v(k,tn) ) * ddy & 2667 + ( flux_t(k) + diss_t(k) & 2668 - flux_d - diss_d ) * ddzw(k) & 2684 + ( ( flux_t(k) + diss_t(k) ) & 2685 - ( flux_d + diss_d ) & 2686 ) * drho_air(k) * ddzw(k) & 2669 2687 ) + v(k,j,i) * div 2670 2688 … … 2714 2732 2715 2733 USE arrays_3d, & 2716 ONLY: ddzu, diss_l_w, diss_s_w, flux_l_w, flux_s_w, tend, u, v, w 2734 ONLY: ddzu, diss_l_w, diss_s_w, flux_l_w, flux_s_w, tend, u, v, w,& 2735 drho_air_zw, rho_air, rho_air_zw 2717 2736 2718 2737 USE constants, & … … 2984 3003 2985 3004 w_comp = w(k+1,j,i) + w(k,j,i) 2986 flux_t(k) = w_comp * (&3005 flux_t(k) = w_comp * rho_air(k+1) * ( & 2987 3006 ( 37.0_wp * ibit35 * adv_mom_5 & 2988 3007 + 7.0_wp * ibit34 * adv_mom_3 & … … 2999 3018 ) 3000 3019 3001 diss_t(k) = - ABS( w_comp ) * (&3020 diss_t(k) = - ABS( w_comp ) * rho_air(k+1) * ( & 3002 3021 ( 10.0_wp * ibit35 * adv_mom_5 & 3003 3022 + 3.0_wp * ibit34 * adv_mom_3 & … … 3024 3043 + IBITS(wall_flags_0(k,j,i-1),29,1) & 3025 3044 ) & 3026 ) * ddx&3045 ) * rho_air_zw(k) * ddx & 3027 3046 + ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 ) & 3028 3047 - ( v(k+1,j,i) + v(k,j,i) ) & … … 3031 3050 + IBITS(wall_flags_00(k,j-1,i),0,1) & 3032 3051 ) & 3033 ) * ddy&3034 + ( w_comp * ( ibit33 + ibit34 + ibit35 )&3035 - ( w(k,j,i) + w(k-1,j,i) ) 3052 ) * rho_air_zw(k) * ddy & 3053 + ( w_comp * rho_air(k+1) * ( ibit33 + ibit34 + ibit35 ) & 3054 - ( w(k,j,i) + w(k-1,j,i) ) * rho_air(k) & 3036 3055 * ( IBITS(wall_flags_00(k-1,j,i),1,1) & 3037 3056 + IBITS(wall_flags_00(k-1,j,i),2,1) & … … 3047 3066 + ( flux_n(k) + diss_n(k) & 3048 3067 - flux_s_w(k,tn) - diss_s_w(k,tn) ) * ddy & 3049 + ( flux_t(k) + diss_t(k) & 3050 - flux_d - diss_d ) * ddzu(k+1) & 3068 + ( ( flux_t(k) + diss_t(k) ) & 3069 - ( flux_d + diss_d ) & 3070 ) * drho_air_zw(k) * ddzu(k+1) & 3051 3071 ) + div * w(k,j,i) 3052 3072 … … 3105 3125 3106 3126 w_comp = w(k+1,j,i) + w(k,j,i) 3107 flux_t(k) = w_comp * (&3127 flux_t(k) = w_comp * rho_air(k+1) * ( & 3108 3128 ( 37.0_wp * ibit35 * adv_mom_5 & 3109 3129 + 7.0_wp * ibit34 * adv_mom_3 & … … 3120 3140 ) 3121 3141 3122 diss_t(k) = - ABS( w_comp ) * (&3142 diss_t(k) = - ABS( w_comp ) * rho_air(k+1) * ( & 3123 3143 ( 10.0_wp * ibit35 * adv_mom_5 & 3124 3144 + 3.0_wp * ibit34 * adv_mom_3 & … … 3139 3159 !-- by a not sufficient reduction of divergences near topography. 3140 3160 div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i) ) ) * ddx & 3161 * rho_air_zw(k) & 3141 3162 + ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i) ) ) * ddy & 3142 + ( w_comp - ( w(k,j,i) + w(k-1,j,i) ) ) * ddzu(k+1) & 3163 * rho_air_zw(k) & 3164 + ( w_comp * rho_air(k+1) - & 3165 ( w(k,j,i) + w(k-1,j,i) ) * rho_air(k) & 3166 ) * ddzu(k+1) & 3143 3167 ) * 0.5_wp 3144 3168 … … 3148 3172 + ( flux_n(k) + diss_n(k) & 3149 3173 - flux_s_w(k,tn) - diss_s_w(k,tn) ) * ddy & 3150 + ( flux_t(k) + diss_t(k) & 3151 - flux_d - diss_d ) * ddzu(k+1) & 3174 + ( ( flux_t(k) + diss_t(k) ) & 3175 - ( flux_d + diss_d ) & 3176 ) * drho_air_zw(k) * ddzu(k+1) & 3152 3177 ) + div * w(k,j,i) 3153 3178 … … 3183 3208 3184 3209 USE arrays_3d, & 3185 ONLY: ddzw, tend, u, v,w3210 ONLY: ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw 3186 3211 3187 3212 USE constants, & … … 3484 3509 3485 3510 3486 flux_t(k) = w(k,j,i) * (&3511 flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( & 3487 3512 ( 37.0_wp * ibit8 * adv_sca_5 & 3488 3513 + 7.0_wp * ibit7 * adv_sca_3 & … … 3498 3523 ) 3499 3524 3500 diss_t(k) = -ABS( w(k,j,i) ) * (&3525 diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( & 3501 3526 ( 10.0_wp * ibit8 * adv_sca_5 & 3502 3527 + 3.0_wp * ibit7 * adv_sca_3 & … … 3539 3564 fd_1 = ( w(k-1,j,i) * ( sk(k,j,i) + sk(k-1,j,i) ) & 3540 3565 -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) ) & 3541 ) * adv_sca_1 3566 ) * adv_sca_1 * rho_air_zw(k) 3542 3567 3543 3568 ft_1 = ( w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) & 3544 3569 -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) ) & 3545 ) * adv_sca_1 3570 ) * adv_sca_1 * rho_air_zw(k) 3546 3571 ! 3547 3572 !-- Calculate ratio of upwind gradients. Note, Min/Max is just … … 3659 3684 + IBITS(wall_flags_0(k,j,i-1),2,1) & 3660 3685 ) & 3661 ) * ddx&3686 ) * rho_air(k) * ddx & 3662 3687 + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 ) & 3663 3688 - v(k,j,i) * ( IBITS(wall_flags_0(k,j-1,i),3,1) & … … 3665 3690 + IBITS(wall_flags_0(k,j-1,i),5,1) & 3666 3691 ) & 3667 ) * ddy & 3668 + ( w(k,j,i) * ( ibit6 + ibit7 + ibit8 ) & 3669 - w(k-1,j,i) * ( IBITS(wall_flags_0(k-1,j,i),6,1) & 3692 ) * rho_air(k) * ddy & 3693 + ( w(k,j,i) * rho_air_zw(k) * & 3694 ( ibit6 + ibit7 + ibit8 ) & 3695 - w(k-1,j,i) * rho_air_zw(k-1) * & 3696 ( IBITS(wall_flags_0(k-1,j,i),6,1) & 3670 3697 + IBITS(wall_flags_0(k-1,j,i),7,1) & 3671 3698 + IBITS(wall_flags_0(k-1,j,i),8,1) & … … 3679 3706 + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k) - & 3680 3707 swap_diss_y_local(k) ) * ddy & 3681 + ( flux_t(k) + diss_t(k) - flux_d - diss_d & 3682 ) * ddzw(k) & 3708 + ( ( flux_t(k) + diss_t(k) ) - & 3709 ( flux_d + diss_d ) & 3710 ) * drho_air(k) * ddzw(k) & 3683 3711 ) + sk(k,j,i) * div 3684 3712 … … 3725 3753 3726 3754 3727 flux_t(k) = w(k,j,i) * (&3755 flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( & 3728 3756 ( 37.0_wp * ibit8 * adv_sca_5 & 3729 3757 + 7.0_wp * ibit7 * adv_sca_3 & … … 3739 3767 ) 3740 3768 3741 diss_t(k) = -ABS( w(k,j,i) ) * (&3769 diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( & 3742 3770 ( 10.0_wp * ibit8 * adv_sca_5 & 3743 3771 + 3.0_wp * ibit7 * adv_sca_3 & … … 3780 3808 fd_1 = ( w(k-1,j,i) * ( sk(k,j,i) + sk(k-1,j,i) ) & 3781 3809 -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) ) & 3782 ) * adv_sca_1 3810 ) * adv_sca_1 * rho_air_zw(k) 3783 3811 3784 3812 ft_1 = ( w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) & 3785 3813 -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) ) & 3786 ) * adv_sca_1 3814 ) * adv_sca_1 * rho_air_zw(k) 3787 3815 ! 3788 3816 !-- Calculate ratio of upwind gradients. Note, Min/Max is just … … 3895 3923 !-- correction is needed to overcome numerical instabilities introduced 3896 3924 !-- by a not sufficient reduction of divergences near topography. 3897 div = ( u(k,j,i+1) - u(k,j,i) ) * ddx & 3898 + ( v(k,j+1,i) - v(k,j,i) ) * ddy & 3899 + ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 3925 div = ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx & 3926 + ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy & 3927 + ( w(k,j,i) * rho_air_zw(k) - & 3928 w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k) 3900 3929 3901 3930 tend(k,j,i) = tend(k,j,i) - ( & … … 3904 3933 + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k) - & 3905 3934 swap_diss_y_local(k) ) * ddy & 3906 + ( flux_t(k) + diss_t(k) - flux_d - diss_d & 3907 ) * ddzw(k) & 3935 + ( ( flux_t(k) + diss_t(k) ) - & 3936 ( flux_d + diss_d ) & 3937 ) * drho_air(k) * ddzw(k) & 3908 3938 ) + sk(k,j,i) * div 3909 3939 … … 4004 4034 4005 4035 USE arrays_3d, & 4006 ONLY: ddzw, tend, u, v,w4036 ONLY: ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw 4007 4037 4008 4038 USE constants, & … … 4278 4308 k_mm = k - 2 * ibit8 4279 4309 4280 flux_t = w(k,j,i) * (&4310 flux_t = w(k,j,i) * rho_air_zw(k) * ( & 4281 4311 ( 37.0_wp * ibit8 * adv_sca_5 & 4282 4312 + 7.0_wp * ibit7 * adv_sca_3 & … … 4292 4322 ) 4293 4323 4294 diss_t = -ABS( w(k,j,i) ) * (&4324 diss_t = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( & 4295 4325 ( 10.0_wp * ibit8 * adv_sca_5 & 4296 4326 + 3.0_wp * ibit7 * adv_sca_3 & … … 4333 4363 fd_1 = ( w(k-1,j,i) * ( sk(k,j,i) + sk(k-1,j,i) ) & 4334 4364 -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) ) & 4335 ) * adv_sca_1 4365 ) * adv_sca_1 * rho_air_zw(k) 4336 4366 4337 4367 ft_1 = ( w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) & 4338 4368 -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) ) & 4339 ) * adv_sca_1 4369 ) * adv_sca_1 * rho_air_zw(k) 4340 4370 ! 4341 4371 !-- Calculate ratio of upwind gradients. Note, Min/Max is just … … 4447 4477 + IBITS(wall_flags_0(k,j,i-1),2,1) & 4448 4478 ) & 4449 ) * ddx&4479 ) * rho_air(k) * ddx & 4450 4480 + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 ) & 4451 4481 - v(k,j,i) * ( IBITS(wall_flags_0(k,j-1,i),3,1) & … … 4453 4483 + IBITS(wall_flags_0(k,j-1,i),5,1) & 4454 4484 ) & 4455 ) * ddy & 4456 + ( w(k,j,i) * ( ibit6 + ibit7 + ibit8 ) & 4457 - w(k-1,j,i) * ( IBITS(wall_flags_0(k-1,j,i),6,1) & 4485 ) * rho_air(k) * ddy & 4486 + ( w(k,j,i) * rho_air_zw(k) * & 4487 ( ibit6 + ibit7 + ibit8 ) & 4488 - w(k-1,j,i) * rho_air_zw(k-1) * & 4489 ( IBITS(wall_flags_0(k-1,j,i),6,1) & 4458 4490 + IBITS(wall_flags_0(k-1,j,i),7,1) & 4459 4491 + IBITS(wall_flags_0(k-1,j,i),8,1) & … … 4465 4497 ( flux_r + diss_r - flux_l - diss_l ) * ddx & 4466 4498 + ( flux_n + diss_n - flux_s - diss_s ) * ddy & 4467 + ( flux_t + diss_t - flux_d - diss_d ) * ddzw(k)& 4499 + ( ( flux_t + diss_t ) - & 4500 ( flux_d + diss_d ) & 4501 ) * drho_air(k) * ddzw(k) & 4468 4502 ) + div * sk(k,j,i) 4469 4503 … … 4511 4545 4512 4546 USE arrays_3d, & 4513 ONLY: ddzw, tend, u, v,w4547 ONLY: ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw 4514 4548 4515 4549 USE constants, & … … 4780 4814 4781 4815 w_comp = w(k,j,i) + w(k,j,i-1) 4782 flux_t(k) = w_comp * (&4816 flux_t(k) = w_comp * rho_air_zw(k) * ( & 4783 4817 ( 37.0_wp * ibit17 * adv_mom_5 & 4784 4818 + 7.0_wp * ibit16 * adv_mom_3 & … … 4795 4829 ) 4796 4830 4797 diss_t(k) = - ABS( w_comp ) * (&4831 diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * ( & 4798 4832 ( 10.0_wp * ibit17 * adv_mom_5 & 4799 4833 + 3.0_wp * ibit16 * adv_mom_3 & … … 4818 4852 + IBITS(wall_flags_0(k,j,i-1),10,1) & 4819 4853 + IBITS(wall_flags_0(k,j,i-1),11,1) & 4820 ) & 4821 ) * ddx&4854 ) & 4855 ) * rho_air(k) * ddx & 4822 4856 + ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 ) & 4823 4857 - ( v(k,j,i) + v(k,j,i-1 ) ) & … … 4826 4860 + IBITS(wall_flags_0(k,j-1,i),14,1) & 4827 4861 ) & 4828 ) * ddy&4829 + ( w_comp * ( ibit15 + ibit16 + ibit17 )&4830 - ( w(k-1,j,i) + w(k-1,j,i-1) ) 4862 ) * rho_air(k) * ddy & 4863 + ( w_comp * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 ) & 4864 - ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1) & 4831 4865 * ( IBITS(wall_flags_0(k-1,j,i),15,1) & 4832 4866 + IBITS(wall_flags_0(k-1,j,i),16,1) & … … 4843 4877 + ( flux_n(k) + diss_n(k) & 4844 4878 - swap_flux_y_local_u(k) - swap_diss_y_local_u(k) ) * ddy & 4845 + ( flux_t(k) + diss_t(k)&4846 - flux_d - diss_d&4847 4879 + ( ( flux_t(k) + diss_t(k) ) & 4880 - ( flux_d + diss_d ) & 4881 ) * drho_air(k) * ddzw(k) & 4848 4882 ) + div * u(k,j,i) 4849 4883 … … 4911 4945 4912 4946 w_comp = w(k,j,i) + w(k,j,i-1) 4913 flux_t(k) = w_comp * (&4947 flux_t(k) = w_comp * rho_air_zw(k) * ( & 4914 4948 ( 37.0_wp * ibit17 * adv_mom_5 & 4915 4949 + 7.0_wp * ibit16 * adv_mom_3 & … … 4926 4960 ) 4927 4961 4928 diss_t(k) = - ABS( w_comp ) * (&4962 diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * ( & 4929 4963 ( 10.0_wp * ibit17 * adv_mom_5 & 4930 4964 + 3.0_wp * ibit16 * adv_mom_3 & … … 4945 4979 !-- by a not sufficient reduction of divergences near topography. 4946 4980 div = ( ( u_comp(k) - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx & 4981 * rho_air(k)& 4947 4982 + ( v_comp + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy & 4948 + ( w_comp - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) & 4949 * ddzw(k) & 4983 * rho_air(k)& 4984 + ( w_comp * rho_air_zw(k) - & 4985 ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1) & 4986 ) * ddzw(k) & 4950 4987 ) * 0.5_wp 4951 4988 … … 4955 4992 + ( flux_n(k) + diss_n(k) & 4956 4993 - swap_flux_y_local_u(k) - swap_diss_y_local_u(k) ) * ddy & 4957 + ( flux_t(k) + diss_t(k)&4958 - flux_d - diss_d&4959 4994 + ( ( flux_t(k) + diss_t(k) ) & 4995 - ( flux_d + diss_d ) & 4996 ) * drho_air(k) * ddzw(k) & 4960 4997 ) + div * u(k,j,i) 4961 4998 … … 5004 5041 5005 5042 USE arrays_3d, & 5006 ONLY: ddzw, tend, u, v,w5043 ONLY: ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw 5007 5044 5008 5045 USE constants, & … … 5267 5304 5268 5305 w_comp = w(k,j,i) + w(k,j,i-1) 5269 flux_t = w_comp * (&5306 flux_t = w_comp * rho_air_zw(k) * ( & 5270 5307 ( 37.0_wp * ibit17 * adv_mom_5 & 5271 5308 + 7.0_wp * ibit16 * adv_mom_3 & … … 5282 5319 ) 5283 5320 5284 diss_t = - ABS( w_comp ) * (&5321 diss_t = - ABS( w_comp ) * rho_air_zw(k) * ( & 5285 5322 ( 10.0_wp * ibit17 * adv_mom_5 & 5286 5323 + 3.0_wp * ibit16 * adv_mom_3 & … … 5305 5342 + IBITS(wall_flags_0(k,j,i-1),10,1) & 5306 5343 + IBITS(wall_flags_0(k,j,i-1),11,1) & 5307 ) & 5308 ) * ddx&5344 ) & 5345 ) * rho_air(k) * ddx & 5309 5346 + ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 ) & 5310 5347 - ( v(k,j,i) + v(k,j,i-1 ) ) & … … 5313 5350 + IBITS(wall_flags_0(k,j-1,i),14,1) & 5314 5351 ) & 5315 ) * ddy&5316 + ( w_comp * ( ibit15 + ibit16 + ibit17 )&5317 - ( w(k-1,j,i) + w(k-1,j,i-1) ) 5352 ) * rho_air(k) * ddy & 5353 + ( w_comp * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 ) & 5354 - ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1) & 5318 5355 * ( IBITS(wall_flags_0(k-1,j,i),15,1) & 5319 5356 + IBITS(wall_flags_0(k-1,j,i),16,1) & 5320 5357 + IBITS(wall_flags_0(k-1,j,i),17,1) & 5321 ) & 5358 ) & 5322 5359 ) * ddzw(k) & 5323 5360 ) * 0.5_wp … … 5327 5364 ( flux_r + diss_r - flux_l - diss_l ) * ddx & 5328 5365 + ( flux_n + diss_n - flux_s - diss_s ) * ddy & 5329 + ( flux_t + diss_t - flux_d - diss_d ) * ddzw(k) & 5366 + ( ( flux_t + diss_t ) - & 5367 ( flux_d + diss_d ) & 5368 ) * drho_air(k) * ddzw(k) & 5330 5369 ) + div * u(k,j,i) 5331 5370 … … 5365 5404 5366 5405 USE arrays_3d, & 5367 ONLY: ddzw, tend, u, v,w5406 ONLY: ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw 5368 5407 5369 5408 USE constants, & … … 5633 5672 5634 5673 w_comp = w(k,j-1,i) + w(k,j,i) 5635 flux_t(k) = w_comp * (&5674 flux_t(k) = w_comp * rho_air_zw(k) * ( & 5636 5675 ( 37.0_wp * ibit26 * adv_mom_5 & 5637 5676 + 7.0_wp * ibit25 * adv_mom_3 & … … 5648 5687 ) 5649 5688 5650 diss_t(k) = - ABS( w_comp ) * (&5689 diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * ( & 5651 5690 ( 10.0_wp * ibit26 * adv_mom_5 & 5652 5691 + 3.0_wp * ibit25 * adv_mom_3 & … … 5672 5711 + IBITS(wall_flags_0(k,j,i-1),19,1) & 5673 5712 + IBITS(wall_flags_0(k,j,i-1),20,1) & 5674 ) & 5675 ) * ddx&5713 ) & 5714 ) * rho_air(k) * ddx & 5676 5715 + ( v_comp(k) & 5677 5716 * ( ibit21 + ibit22 + ibit23 ) & … … 5680 5719 + IBITS(wall_flags_0(k,j-1,i),22,1) & 5681 5720 + IBITS(wall_flags_0(k,j-1,i),23,1) & 5682 ) & 5683 ) * ddy&5684 + ( w_comp 5721 ) & 5722 ) * rho_air(k) * ddy & 5723 + ( w_comp * rho_air_zw(k) & 5685 5724 * ( ibit24 + ibit25 + ibit26 ) & 5686 - ( w(k-1,j-1,i) + w(k-1,j,i) ) 5725 - ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1) & 5687 5726 * ( IBITS(wall_flags_0(k-1,j,i),24,1) & 5688 5727 + IBITS(wall_flags_0(k-1,j,i),25,1) & … … 5700 5739 - swap_flux_y_local_v(k) - swap_diss_y_local_v(k) & 5701 5740 ) * ddy & 5702 + ( flux_t(k) + diss_t(k)&5703 - flux_d - diss_d&5704 ) * d dzw(k)&5741 + ( ( flux_t(k) + diss_t(k) ) & 5742 - ( flux_d + diss_d ) & 5743 ) * drho_air(k) * ddzw(k) & 5705 5744 ) + v(k,j,i) * div 5706 5745 … … 5772 5811 5773 5812 w_comp = w(k,j-1,i) + w(k,j,i) 5774 flux_t(k) = w_comp * (&5813 flux_t(k) = w_comp * rho_air_zw(k) * ( & 5775 5814 ( 37.0_wp * ibit26 * adv_mom_5 & 5776 5815 + 7.0_wp * ibit25 * adv_mom_3 & … … 5787 5826 ) 5788 5827 5789 diss_t(k) = - ABS( w_comp ) * (&5828 diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * ( & 5790 5829 ( 10.0_wp * ibit26 * adv_mom_5 & 5791 5830 + 3.0_wp * ibit25 * adv_mom_3 & … … 5806 5845 !-- by a not sufficient reduction of divergences near topography. 5807 5846 div = ( ( u_comp + gu - ( u(k,j-1,i) + u(k,j,i) ) ) * ddx & 5847 * rho_air(k)& 5808 5848 + ( v_comp(k) - ( v(k,j,i) + v(k,j-1,i) ) ) * ddy & 5809 + ( w_comp - ( w(k-1,j-1,i) + w(k-1,j,i) ) ) & 5810 * ddzw(k) & 5849 * rho_air(k)& 5850 + ( w_comp * rho_air_zw(k) - & 5851 ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1) & 5852 ) * ddzw(k) & 5811 5853 ) * 0.5_wp 5812 5854 … … 5818 5860 - swap_flux_y_local_v(k) - swap_diss_y_local_v(k) & 5819 5861 ) * ddy & 5820 + ( flux_t(k) + diss_t(k)&5821 - flux_d - diss_d&5822 ) * d dzw(k)&5862 + ( ( flux_t(k) + diss_t(k) ) & 5863 - ( flux_d + diss_d ) & 5864 ) * drho_air(k) * ddzw(k) & 5823 5865 ) + v(k,j,i) * div 5824 5866 … … 5869 5911 5870 5912 USE arrays_3d, & 5871 ONLY: ddzw, tend, u, v,w5913 ONLY: ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw 5872 5914 5873 5915 USE constants, & … … 6132 6174 6133 6175 w_comp = w(k,j-1,i) + w(k,j,i) 6134 flux_t = w_comp * (&6176 flux_t = w_comp * rho_air_zw(k) * ( & 6135 6177 ( 37.0_wp * ibit26 * adv_mom_5 & 6136 6178 + 7.0_wp * ibit25 * adv_mom_3 & … … 6147 6189 ) 6148 6190 6149 diss_t = - ABS( w_comp ) * (&6191 diss_t = - ABS( w_comp ) * rho_air_zw(k) * ( & 6150 6192 ( 10.0_wp * ibit26 * adv_mom_5 & 6151 6193 + 3.0_wp * ibit25 * adv_mom_3 & … … 6171 6213 + IBITS(wall_flags_0(k,j,i-1),19,1) & 6172 6214 + IBITS(wall_flags_0(k,j,i-1),20,1) & 6173 ) & 6174 ) * ddx&6215 ) & 6216 ) * rho_air(k) * ddx & 6175 6217 + ( v_comp & 6176 6218 * ( ibit21 + ibit22 + ibit23 ) & … … 6179 6221 + IBITS(wall_flags_0(k,j-1,i),22,1) & 6180 6222 + IBITS(wall_flags_0(k,j-1,i),23,1) & 6181 ) & 6182 ) * ddy&6183 + ( w_comp 6223 ) & 6224 ) * rho_air(k) * ddy & 6225 + ( w_comp * rho_air_zw(k) & 6184 6226 * ( ibit24 + ibit25 + ibit26 ) & 6185 - ( w(k-1,j-1,i) + w(k-1,j,i) ) 6227 - ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1) & 6186 6228 * ( IBITS(wall_flags_0(k-1,j,i),24,1) & 6187 6229 + IBITS(wall_flags_0(k-1,j,i),25,1) & … … 6195 6237 ( flux_r + diss_r - flux_l - diss_l ) * ddx & 6196 6238 + ( flux_n + diss_n - flux_s - diss_s ) * ddy & 6197 + ( flux_t + diss_t - flux_d - diss_d ) * ddzw(k) & 6239 + ( ( flux_t + diss_t ) - & 6240 ( flux_d + diss_d ) & 6241 ) * drho_air(k) * ddzw(k) & 6198 6242 ) + div * v(k,j,i) 6199 6243 … … 6235 6279 6236 6280 USE arrays_3d, & 6237 ONLY: ddzu, tend, u, v,w6281 ONLY: ddzu, drho_air_zw, tend, u, v, w, rho_air, rho_air_zw 6238 6282 6239 6283 USE constants, & … … 6509 6553 6510 6554 w_comp = w(k+1,j,i) + w(k,j,i) 6511 flux_t(k) = w_comp * (&6555 flux_t(k) = w_comp * rho_air(k+1) * ( & 6512 6556 ( 37.0_wp * ibit35 * adv_mom_5 & 6513 6557 + 7.0_wp * ibit34 * adv_mom_3 & … … 6524 6568 ) 6525 6569 6526 diss_t(k) = - ABS( w_comp ) * (&6570 diss_t(k) = - ABS( w_comp ) * rho_air(k+1) * ( & 6527 6571 ( 10.0_wp * ibit35 * adv_mom_5 & 6528 6572 + 3.0_wp * ibit34 * adv_mom_3 & … … 6548 6592 + IBITS(wall_flags_0(k,j,i-1),29,1) & 6549 6593 ) & 6550 ) * ddx&6551 + ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 ) & 6594 ) * rho_air_zw(k) * ddx & 6595 + ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 ) & 6552 6596 - ( v(k+1,j,i) + v(k,j,i) ) & 6553 6597 * ( IBITS(wall_flags_0(k,j-1,i),30,1) & … … 6555 6599 + IBITS(wall_flags_00(k,j-1,i),0,1) & 6556 6600 ) & 6557 ) * ddy&6558 + ( w_comp * ( ibit33 + ibit34 + ibit35 )&6559 - ( w(k,j,i) + w(k-1,j,i) ) 6601 ) * rho_air_zw(k) * ddy & 6602 + ( w_comp * rho_air(k+1) * ( ibit33 + ibit34 + ibit35 ) & 6603 - ( w(k,j,i) + w(k-1,j,i) ) * rho_air(k) & 6560 6604 * ( IBITS(wall_flags_00(k-1,j,i),1,1) & 6561 6605 + IBITS(wall_flags_00(k-1,j,i),2,1) & … … 6574 6618 - swap_flux_y_local_w(k) - swap_diss_y_local_w(k) & 6575 6619 ) * ddy & 6576 + ( flux_t(k) + diss_t(k)&6577 - flux_d - diss_d&6578 ) * d dzu(k+1)&6620 + ( ( flux_t(k) + diss_t(k) ) & 6621 - ( flux_d + diss_d ) & 6622 ) * drho_air_zw(k) * ddzu(k+1) & 6579 6623 ) + div * w(k,j,i) 6580 6624 … … 6632 6676 6633 6677 w_comp = w(k+1,j,i) + w(k,j,i) 6634 flux_t(k) = w_comp * (&6678 flux_t(k) = w_comp * rho_air(k+1) * ( & 6635 6679 ( 37.0_wp * ibit35 * adv_mom_5 & 6636 6680 + 7.0_wp * ibit34 * adv_mom_3 & … … 6647 6691 ) 6648 6692 6649 diss_t(k) = - ABS( w_comp ) * (&6693 diss_t(k) = - ABS( w_comp ) * rho_air(k+1) * ( & 6650 6694 ( 10.0_wp * ibit35 * adv_mom_5 & 6651 6695 + 3.0_wp * ibit34 * adv_mom_3 & … … 6666 6710 !-- by a not sufficient reduction of divergences near topography. 6667 6711 div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i) ) ) * ddx & 6712 * rho_air_zw(k) & 6668 6713 + ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i) ) ) * ddy & 6669 + ( w_comp - ( w(k,j,i) + w(k-1,j,i) ) ) & 6670 * ddzu(k+1) & 6714 * rho_air_zw(k) & 6715 + ( w_comp * rho_air(k+1) - & 6716 ( w(k,j,i) + w(k-1,j,i) ) * rho_air(k) & 6717 ) * ddzu(k+1) & 6671 6718 ) * 0.5_wp 6672 6719 … … 6678 6725 - swap_flux_y_local_w(k) - swap_diss_y_local_w(k) & 6679 6726 ) * ddy & 6680 + ( flux_t(k) + diss_t(k)&6681 - flux_d - diss_d&6682 ) * d dzu(k+1)&6727 + ( ( flux_t(k) + diss_t(k) ) & 6728 - ( flux_d + diss_d ) & 6729 ) * drho_air_zw(k) * ddzu(k+1) & 6683 6730 ) + div * w(k,j,i) 6684 6731 … … 6714 6761 6715 6762 USE arrays_3d, & 6716 ONLY: ddzu, tend, u, v,w6763 ONLY: ddzu, drho_air_zw, tend, u, v, w, rho_air, rho_air_zw 6717 6764 6718 6765 USE constants, & … … 6976 7023 6977 7024 w_comp = w(k+1,j,i) + w(k,j,i) 6978 flux_t = w_comp * (&7025 flux_t = w_comp * rho_air(k+1) * ( & 6979 7026 ( 37.0_wp * ibit35 * adv_mom_5 & 6980 7027 + 7.0_wp * ibit34 * adv_mom_3 & … … 6991 7038 ) 6992 7039 6993 diss_t = - ABS( w_comp ) * (&7040 diss_t = - ABS( w_comp ) * rho_air(k+1) * ( & 6994 7041 ( 10.0_wp * ibit35 * adv_mom_5 & 6995 7042 + 3.0_wp * ibit34 * adv_mom_3 & … … 7015 7062 + IBITS(wall_flags_0(k,j,i-1),29,1) & 7016 7063 ) & 7017 ) * ddx&7018 + ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 ) & 7064 ) * rho_air_zw(k) * ddx & 7065 + ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 ) & 7019 7066 - ( v(k+1,j,i) + v(k,j,i) ) & 7020 7067 * ( IBITS(wall_flags_0(k,j-1,i),30,1) & … … 7022 7069 + IBITS(wall_flags_00(k,j-1,i),0,1) & 7023 7070 ) & 7024 ) * ddy&7025 + ( w_comp * ( ibit33 + ibit34 + ibit35 )&7026 - ( w(k,j,i) + w(k-1,j,i) ) 7071 ) * rho_air_zw(k) * ddy & 7072 + ( w_comp * rho_air(k+1) * ( ibit33 + ibit34 + ibit35 ) & 7073 - ( w(k,j,i) + w(k-1,j,i) ) * rho_air(k) & 7027 7074 * ( IBITS(wall_flags_00(k-1,j,i),1,1) & 7028 7075 + IBITS(wall_flags_00(k-1,j,i),2,1) & … … 7036 7083 ( flux_r + diss_r - flux_l - diss_l ) * ddx & 7037 7084 + ( flux_n + diss_n - flux_s - diss_s ) * ddy & 7038 + ( flux_t + diss_t - flux_d - diss_d ) * ddzu(k+1) & 7085 + ( ( flux_t + diss_t ) - & 7086 ( flux_d + diss_d ) * rho_air(k) & 7087 ) * drho_air_zw(k) * ddzu(k+1) & 7039 7088 ) + div * w(k,j,i) 7040 7089 -
palm/trunk/SOURCE/check_parameters.f90
r2032 r2037 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Anelastic approximation implemented 23 23 ! 24 24 ! Former revisions: … … 432 432 USE netcdf_interface, & 433 433 ONLY: dopr_unit, do2d_unit, do3d_unit, netcdf_data_format, & 434 netcdf_data_format_string 434 netcdf_data_format_string, dots_unit, heatflux_output_unit, & 435 waterflux_output_unit, momentumflux_output_unit 435 436 436 437 USE particle_attributes … … 837 838 838 839 ENDIF 840 841 ! 842 !-- Check approximation 843 IF ( TRIM( approximation ) /= 'boussinesq' .AND. & 844 TRIM( approximation ) /= 'anelastic' ) THEN 845 message_string = 'unknown approximation: approximation = "' // & 846 TRIM( approximation ) // '"' 847 CALL message( 'check_parameters', 'PA0000', 1, 2, 0, 6, 0 ) 848 ENDIF 849 850 ! 851 !-- Check approximation requirements 852 IF ( TRIM( approximation ) == 'anelastic' .AND. & 853 TRIM( momentum_advec ) /= 'ws-scheme' ) THEN 854 message_string = 'Anelastic approximation requires: ' // & 855 'momentum_advec = "ws-scheme"' 856 CALL message( 'check_parameters', 'PA0000', 1, 2, 0, 6, 0 ) 857 ENDIF 858 IF ( TRIM( approximation ) == 'anelastic' .AND. & 859 TRIM( psolver ) == 'multigrid' ) THEN 860 message_string = 'Anelastic approximation currently only supports: ' // & 861 'psolver = "poisfft", ' // & 862 'psolver = "sor" and ' // & 863 'psolver = "multigrid_noopt"' 864 CALL message( 'check_parameters', 'PA0000', 1, 2, 0, 6, 0 ) 865 ENDIF 866 IF ( TRIM( approximation ) == 'anelastic' .AND. & 867 conserve_volume_flow ) THEN 868 message_string = 'Anelastic approximation is not allowed with:' // & 869 'conserve_volume_flow = .TRUE.' 870 CALL message( 'check_parameters', 'PA0000', 1, 2, 0, 6, 0 ) 871 ENDIF 872 873 ! 874 !-- Check flux input mode 875 IF ( TRIM( flux_input_mode ) /= 'dynamic' .AND. & 876 TRIM( flux_input_mode ) /= 'kinematic' .AND. & 877 TRIM( flux_input_mode ) /= 'approximation-specific' ) THEN 878 message_string = 'unknown flux input mode: flux_input_mode = "' // & 879 TRIM( flux_input_mode ) // '"' 880 CALL message( 'check_parameters', 'PA0000', 1, 2, 0, 6, 0 ) 881 ENDIF 882 !-- Set flux input mode according to approximation if applicable 883 IF ( TRIM( flux_input_mode ) == 'approximation-specific' ) THEN 884 IF ( TRIM( approximation ) == 'anelastic' ) THEN 885 flux_input_mode = 'dynamic' 886 ELSEIF ( TRIM( approximation ) == 'boussinesq' ) THEN 887 flux_input_mode = 'kinematic' 888 ENDIF 889 ENDIF 890 891 ! 892 !-- Check flux output mode 893 IF ( TRIM( flux_output_mode ) /= 'dynamic' .AND. & 894 TRIM( flux_output_mode ) /= 'kinematic' .AND. & 895 TRIM( flux_output_mode ) /= 'approximation-specific' ) THEN 896 message_string = 'unknown flux output mode: flux_output_mode = "' // & 897 TRIM( flux_output_mode ) // '"' 898 CALL message( 'check_parameters', 'PA0000', 1, 2, 0, 6, 0 ) 899 ENDIF 900 !-- Set flux output mode according to approximation if applicable 901 IF ( TRIM( flux_output_mode ) == 'approximation-specific' ) THEN 902 IF ( TRIM( approximation ) == 'anelastic' ) THEN 903 flux_output_mode = 'dynamic' 904 ELSEIF ( TRIM( approximation ) == 'boussinesq' ) THEN 905 flux_output_mode = 'kinematic' 906 ENDIF 907 ENDIF 908 909 ! 910 !-- set the flux output units according to flux_output_mode 911 IF ( TRIM( flux_output_mode ) == 'kinematic' ) THEN 912 heatflux_output_unit = 'K m/s' 913 waterflux_output_unit = 'kg/kg m/s' 914 momentumflux_output_unit = 'm2/s2' 915 ELSEIF ( TRIM( flux_output_mode ) == 'dynamic' ) THEN 916 heatflux_output_unit = 'W/m2' 917 waterflux_output_unit = 'W/m2' 918 momentumflux_output_unit = 'N/m2' 919 ENDIF 920 921 !-- set time series output units for fluxes 922 dots_unit(14:16) = heatflux_output_unit 923 dots_unit(21) = waterflux_output_unit 924 dots_unit(19:20) = momentumflux_output_unit 839 925 840 926 ! … … 2115 2201 CASE ( 'w"u"' ) 2116 2202 dopr_index(i) = 12 2117 dopr_unit(i) = 'm2/s2'2203 dopr_unit(i) = TRIM ( momentumflux_output_unit ) 2118 2204 hom(:,2,12,:) = SPREAD( zw, 2, statistic_regions+1 ) 2119 2205 IF ( constant_flux_layer ) hom(nzb,2,12,:) = zu(1) … … 2121 2207 CASE ( 'w*u*' ) 2122 2208 dopr_index(i) = 13 2123 dopr_unit(i) = 'm2/s2'2209 dopr_unit(i) = TRIM ( momentumflux_output_unit ) 2124 2210 hom(:,2,13,:) = SPREAD( zw, 2, statistic_regions+1 ) 2125 2211 2126 2212 CASE ( 'w"v"' ) 2127 2213 dopr_index(i) = 14 2128 dopr_unit(i) = 'm2/s2'2214 dopr_unit(i) = TRIM ( momentumflux_output_unit ) 2129 2215 hom(:,2,14,:) = SPREAD( zw, 2, statistic_regions+1 ) 2130 2216 IF ( constant_flux_layer ) hom(nzb,2,14,:) = zu(1) … … 2132 2218 CASE ( 'w*v*' ) 2133 2219 dopr_index(i) = 15 2134 dopr_unit(i) = 'm2/s2'2220 dopr_unit(i) = TRIM ( momentumflux_output_unit ) 2135 2221 hom(:,2,15,:) = SPREAD( zw, 2, statistic_regions+1 ) 2136 2222 2137 2223 CASE ( 'w"pt"' ) 2138 2224 dopr_index(i) = 16 2139 dopr_unit(i) = 'K m/s'2225 dopr_unit(i) = TRIM ( heatflux_output_unit ) 2140 2226 hom(:,2,16,:) = SPREAD( zw, 2, statistic_regions+1 ) 2141 2227 2142 2228 CASE ( 'w*pt*' ) 2143 2229 dopr_index(i) = 17 2144 dopr_unit(i) = 'K m/s'2230 dopr_unit(i) = TRIM ( heatflux_output_unit ) 2145 2231 hom(:,2,17,:) = SPREAD( zw, 2, statistic_regions+1 ) 2146 2232 2147 2233 CASE ( 'wpt' ) 2148 2234 dopr_index(i) = 18 2149 dopr_unit(i) = 'K m/s'2235 dopr_unit(i) = TRIM ( heatflux_output_unit ) 2150 2236 hom(:,2,18,:) = SPREAD( zw, 2, statistic_regions+1 ) 2151 2237 2152 2238 CASE ( 'wu' ) 2153 2239 dopr_index(i) = 19 2154 dopr_unit(i) = 'm2/s2'2240 dopr_unit(i) = TRIM ( momentumflux_output_unit ) 2155 2241 hom(:,2,19,:) = SPREAD( zw, 2, statistic_regions+1 ) 2156 2242 IF ( constant_flux_layer ) hom(nzb,2,19,:) = zu(1) … … 2158 2244 CASE ( 'wv' ) 2159 2245 dopr_index(i) = 20 2160 dopr_unit(i) = 'm2/s2'2246 dopr_unit(i) = TRIM ( momentumflux_output_unit ) 2161 2247 hom(:,2,20,:) = SPREAD( zw, 2, statistic_regions+1 ) 2162 2248 IF ( constant_flux_layer ) hom(nzb,2,20,:) = zu(1) … … 2164 2250 CASE ( 'w*pt*BC' ) 2165 2251 dopr_index(i) = 21 2166 dopr_unit(i) = 'K m/s'2252 dopr_unit(i) = TRIM ( heatflux_output_unit ) 2167 2253 hom(:,2,21,:) = SPREAD( zw, 2, statistic_regions+1 ) 2168 2254 2169 2255 CASE ( 'wptBC' ) 2170 2256 dopr_index(i) = 22 2171 dopr_unit(i) = 'K m/s'2257 dopr_unit(i) = TRIM ( heatflux_output_unit ) 2172 2258 hom(:,2,22,:) = SPREAD( zw, 2, statistic_regions+1 ) 2173 2259 … … 2335 2421 CASE ( 'w"vpt"' ) 2336 2422 dopr_index(i) = 45 2337 dopr_unit(i) = 'K m/s'2423 dopr_unit(i) = TRIM ( heatflux_output_unit ) 2338 2424 hom(:,2,45,:) = SPREAD( zw, 2, statistic_regions+1 ) 2339 2425 2340 2426 CASE ( 'w*vpt*' ) 2341 2427 dopr_index(i) = 46 2342 dopr_unit(i) = 'K m/s'2428 dopr_unit(i) = TRIM ( heatflux_output_unit ) 2343 2429 hom(:,2,46,:) = SPREAD( zw, 2, statistic_regions+1 ) 2344 2430 2345 2431 CASE ( 'wvpt' ) 2346 2432 dopr_index(i) = 47 2347 dopr_unit(i) = 'K m/s'2433 dopr_unit(i) = TRIM ( heatflux_output_unit ) 2348 2434 hom(:,2,47,:) = SPREAD( zw, 2, statistic_regions+1 ) 2349 2435 … … 2356 2442 ELSE 2357 2443 dopr_index(i) = 48 2358 dopr_unit(i) = 'kg/kg m/s'2444 dopr_unit(i) = TRIM ( waterflux_output_unit ) 2359 2445 hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 ) 2360 2446 ENDIF … … 2368 2454 ELSE 2369 2455 dopr_index(i) = 49 2370 dopr_unit(i) = 'kg/kg m/s'2456 dopr_unit(i) = TRIM ( waterflux_output_unit ) 2371 2457 hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 ) 2372 2458 ENDIF … … 2380 2466 ELSE 2381 2467 dopr_index(i) = 50 2382 dopr_unit(i) = 'kg/kg m/s'2468 dopr_unit(i) = TRIM ( waterflux_output_unit ) 2383 2469 hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 ) 2384 2470 ENDIF … … 2423 2509 IF ( humidity .AND. .NOT. cloud_physics ) THEN 2424 2510 dopr_index(i) = 48 2425 dopr_unit(i) = 'kg/kg m/s'2511 dopr_unit(i) = TRIM ( waterflux_output_unit ) 2426 2512 hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 ) 2427 2513 ELSEIF ( humidity .AND. cloud_physics ) THEN 2428 2514 dopr_index(i) = 51 2429 dopr_unit(i) = 'kg/kg m/s'2515 dopr_unit(i) = TRIM ( waterflux_output_unit ) 2430 2516 hom(:,2,51,:) = SPREAD( zw, 2, statistic_regions+1 ) 2431 2517 ELSE … … 2441 2527 THEN 2442 2528 dopr_index(i) = 49 2443 dopr_unit(i) = 'kg/kg m/s'2529 dopr_unit(i) = TRIM ( waterflux_output_unit ) 2444 2530 hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 ) 2445 2531 ELSEIF( humidity .AND. cloud_physics ) THEN 2446 2532 dopr_index(i) = 52 2447 dopr_unit(i) = 'kg/kg m/s'2533 dopr_unit(i) = TRIM ( waterflux_output_unit ) 2448 2534 hom(:,2,52,:) = SPREAD( zw, 2, statistic_regions+1 ) 2449 2535 ELSE … … 2458 2544 IF ( humidity .AND. .NOT. cloud_physics ) THEN 2459 2545 dopr_index(i) = 50 2460 dopr_unit(i) = 'kg/kg m/s'2546 dopr_unit(i) = TRIM ( waterflux_output_unit ) 2461 2547 hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 ) 2462 2548 ELSEIF ( humidity .AND. cloud_physics ) THEN 2463 2549 dopr_index(i) = 53 2464 dopr_unit(i) = 'kg/kg m/s'2550 dopr_unit(i) = TRIM ( waterflux_output_unit ) 2465 2551 hom(:,2,53,:) = SPREAD( zw, 2, statistic_regions+1 ) 2466 2552 ELSE … … 2503 2589 CASE ( 'u"pt"' ) 2504 2590 dopr_index(i) = 58 2505 dopr_unit(i) = 'K m/s'2591 dopr_unit(i) = TRIM ( heatflux_output_unit ) 2506 2592 hom(:,2,58,:) = SPREAD( zu, 2, statistic_regions+1 ) 2507 2593 2508 2594 CASE ( 'u*pt*' ) 2509 2595 dopr_index(i) = 59 2510 dopr_unit(i) = 'K m/s'2596 dopr_unit(i) = TRIM ( heatflux_output_unit ) 2511 2597 hom(:,2,59,:) = SPREAD( zu, 2, statistic_regions+1 ) 2512 2598 2513 2599 CASE ( 'upt_t' ) 2514 2600 dopr_index(i) = 60 2515 dopr_unit(i) = 'K m/s'2601 dopr_unit(i) = TRIM ( heatflux_output_unit ) 2516 2602 hom(:,2,60,:) = SPREAD( zu, 2, statistic_regions+1 ) 2517 2603 2518 2604 CASE ( 'v"pt"' ) 2519 2605 dopr_index(i) = 61 2520 dopr_unit(i) = 'K m/s'2606 dopr_unit(i) = TRIM ( heatflux_output_unit ) 2521 2607 hom(:,2,61,:) = SPREAD( zu, 2, statistic_regions+1 ) 2522 2608 2523 2609 CASE ( 'v*pt*' ) 2524 2610 dopr_index(i) = 62 2525 dopr_unit(i) = 'K m/s'2611 dopr_unit(i) = TRIM ( heatflux_output_unit ) 2526 2612 hom(:,2,62,:) = SPREAD( zu, 2, statistic_regions+1 ) 2527 2613 2528 2614 CASE ( 'vpt_t' ) 2529 2615 dopr_index(i) = 63 2530 dopr_unit(i) = 'K m/s'2616 dopr_unit(i) = TRIM ( heatflux_output_unit ) 2531 2617 hom(:,2,63,:) = SPREAD( zu, 2, statistic_regions+1 ) 2532 2618 … … 2621 2707 CASE ( 'hyp' ) 2622 2708 dopr_index(i) = 72 2623 dopr_unit(i) = ' dbar'2709 dopr_unit(i) = 'hPa' 2624 2710 hom(:,2,72,:) = SPREAD( zu, 2, statistic_regions+1 ) 2711 2712 CASE ( 'rho_air' ) 2713 dopr_index(i) = 121 2714 dopr_unit(i) = 'kg/m3' 2715 hom(:,2,121,:) = SPREAD( zu, 2, statistic_regions+1 ) 2716 2717 CASE ( 'rho_air_zw' ) 2718 dopr_index(i) = 122 2719 dopr_unit(i) = 'kg/m3' 2720 hom(:,2,122,:) = SPREAD( zw, 2, statistic_regions+1 ) 2625 2721 2626 2722 CASE ( 'nr' ) -
palm/trunk/SOURCE/diffusion_e.f90
r2001 r2037 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Anelastic approximation implemented 23 23 ! 24 24 ! Former revisions: … … 128 128 129 129 USE arrays_3d, & 130 ONLY: dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw 130 ONLY: dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw, & 131 drho_air, rho_air_zw 131 132 132 133 USE control_parameters, & … … 222 223 + ( & 223 224 ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) & 225 * rho_air_zw(k) & 224 226 - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) & 225 ) * ddzw(k) & 227 * rho_air_zw(k-1) & 228 ) * ddzw(k) * drho_air(k) & 226 229 - dissipation(k,j) 227 230 … … 294 297 + ( & 295 298 ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) & 299 * rho_air_zw(k) & 296 300 - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) & 297 ) * ddzw(k) & 301 * rho_air_zw(k-1) & 302 ) * ddzw(k) * drho_air(k) & 298 303 - dissipation(k,j) 299 304 … … 339 344 340 345 USE arrays_3d, & 341 ONLY: dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw 346 ONLY: dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw, & 347 drho_air, rho_air_zw 342 348 343 349 USE control_parameters, & … … 430 436 + ( & 431 437 ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) & 438 * rho_air_zw(k) & 432 439 - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) & 433 ) * ddzw(k) & 440 * rho_air_zw(k-1) & 441 ) * ddzw(k) * drho_air(k) & 434 442 - dissipation 435 443 … … 497 505 + ( & 498 506 ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) & 507 * rho_air_zw(k) & 499 508 - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) & 500 ) * ddzw(k) & 509 * rho_air_zw(k-1) & 510 ) * ddzw(k) * drho_air(k) & 501 511 - dissipation 502 512 … … 542 552 543 553 USE arrays_3d, & 544 ONLY: dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw 554 ONLY: dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw, & 555 drho_air, rho_air_zw 545 556 546 557 USE control_parameters, & … … 625 636 + ( & 626 637 ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) & 638 * rho_air_zw(k) & 627 639 - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) & 628 ) * ddzw(k) & 640 * rho_air_zw(k-1) & 641 ) * ddzw(k) * drho_air(k) & 629 642 - dissipation(k) 630 643 -
palm/trunk/SOURCE/diffusion_s.f90
r2001 r2037 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Anelastic approximation implemented 23 23 ! 24 24 ! Former revisions: … … 113 113 114 114 USE arrays_3d, & 115 ONLY: ddzu, ddzw, kh, tend 115 ONLY: ddzu, ddzw, kh, tend, drho_air, rho_air_zw 116 116 117 117 USE control_parameters, & … … 190 190 + 0.5_wp * ( & 191 191 ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) & 192 * rho_air_zw(k) & 192 193 - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k) & 193 ) * ddzw(k) 194 * rho_air_zw(k-1) & 195 ) * ddzw(k) * drho_air(k) 194 196 ENDDO 195 197 … … 205 207 * ( s(k+1,j,i)-s(k,j,i) ) & 206 208 * ddzu(k+1) & 209 * rho_air_zw(k) & 207 210 + s_flux_b(j,i) & 208 ) * ddzw(k) 211 ) * ddzw(k) * drho_air(k) 209 212 210 213 ENDIF … … 222 225 * ( s(k,j,i)-s(k-1,j,i) ) & 223 226 * ddzu(k) & 224 ) * ddzw(k) 227 * rho_air_zw(k-1) & 228 ) * ddzw(k) * drho_air(k) 225 229 226 230 ENDIF … … 240 244 241 245 USE arrays_3d, & 242 ONLY: ddzu, ddzw, kh, tend 246 ONLY: ddzu, ddzw, kh, tend, drho_air, rho_air_zw 243 247 244 248 USE control_parameters, & … … 324 328 + 0.5_wp * ( & 325 329 ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) & 330 * rho_air_zw(k) & 326 331 - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k) & 327 ) * ddzw(k) 332 * rho_air_zw(k-1) & 333 ) * ddzw(k) * drho_air(k) 328 334 ENDIF 329 335 ENDDO … … 338 344 * ( s(k+1,j,i)-s(k,j,i) ) & 339 345 * ddzu(k+1) & 346 * rho_air_zw(k) & 340 347 + s_flux_b(j,i) & 341 ) * ddzw(k) 348 ) * ddzw(k) * drho_air(k) 342 349 ENDIF 343 350 … … 351 358 * ( s(k,j,i)-s(k-1,j,i) ) & 352 359 * ddzu(k) & 353 ) * ddzw(k) 360 * rho_air_zw(k-1) & 361 ) * ddzw(k) * drho_air(k) 354 362 ENDIF 355 363 ENDDO … … 370 378 371 379 USE arrays_3d, & 372 ONLY: ddzu, ddzw, kh, tend 380 ONLY: ddzu, ddzw, kh, tend, drho_air, rho_air_zw 373 381 374 382 USE control_parameters, & … … 446 454 + 0.5_wp * ( & 447 455 ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) & 456 * rho_air_zw(k) & 448 457 - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k) & 449 ) * ddzw(k) 458 * rho_air_zw(k-1) & 459 ) * ddzw(k) * drho_air(k) 450 460 ENDDO 451 461 … … 459 469 * ( s(k+1,j,i)-s(k,j,i) ) & 460 470 * ddzu(k+1) & 471 * rho_air_zw(k) & 461 472 + s_flux_b(j,i) & 462 ) * ddzw(k) 473 ) * ddzw(k) * drho_air(k) 463 474 464 475 ENDIF … … 474 485 * ( s(k,j,i)-s(k-1,j,i) ) & 475 486 * ddzu(k) & 476 ) * ddzw(k) 487 * rho_air_zw(k-1) & 488 ) * ddzw(k) * drho_air(k) 477 489 478 490 ENDIF -
palm/trunk/SOURCE/diffusion_u.f90
r2001 r2037 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Anelastic approximation implemented 23 23 ! 24 24 ! Former revisions: … … 116 116 117 117 USE arrays_3d, & 118 ONLY: ddzu, ddzw, km, tend, u, usws, uswst, v, w 118 ONLY: ddzu, ddzw, km, tend, u, usws, uswst, v, w, & 119 drho_air, rho_air_zw 119 120 120 121 USE control_parameters, & … … 217 218 & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) & 218 219 & + ( w(k,j,i) - w(k,j,i-1) ) * ddx & 219 & ) 220 & ) * rho_air_zw(k) & 220 221 & - kmzm * ( ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k) & 221 222 & + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx & 222 & ) 223 & ) * ddzw(k) 223 & ) * rho_air_zw(k-1) & 224 & ) * ddzw(k) * drho_air(k) 224 225 ENDDO 225 226 … … 242 243 243 244 tend(k,j,i) = tend(k,j,i) & 244 & + ( kmzp * ( w(k,j,i) - w(k,j,i-1) ) * ddx&245 & ) * ddzw(k)&246 & + ( kmzp * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1)&247 & + usws(j,i)&248 & ) * ddzw(k) 245 & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) & 246 & + ( w(k,j,i) - w(k,j,i-1) ) * ddx & 247 & ) * rho_air_zw(k) & 248 & - ( -usws(j,i) ) & 249 & ) * ddzw(k) * drho_air(k) 249 250 ENDIF 250 251 … … 260 261 261 262 tend(k,j,i) = tend(k,j,i) & 262 & - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx&263 & ) * ddzw(k)&264 & + ( -uswst(j,i)&265 & - kmzm * ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k)&266 & ) * ddzw(k) 263 & + ( ( -uswst(j,i) ) & 264 & - kmzm * ( ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k) & 265 & + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx & 266 & ) * rho_air_zw(k-1) & 267 & ) * ddzw(k) * drho_air(k) 267 268 ENDIF 268 269 … … 281 282 282 283 USE arrays_3d, & 283 ONLY: ddzu, ddzw, km, tend, u, usws, uswst, v, w 284 ONLY: ddzu, ddzw, km, tend, u, usws, uswst, v, w, & 285 drho_air, rho_air_zw 284 286 285 287 USE control_parameters, & … … 390 392 & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1)& 391 393 & + ( w(k,j,i) - w(k,j,i-1) ) * ddx & 392 & ) 394 & ) * rho_air_zw(k) & 393 395 & - kmzm * ( ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k)& 394 396 & + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx & 395 & ) 396 & ) * ddzw(k) 397 & ) * rho_air_zw(k-1) & 398 & ) * ddzw(k) * drho_air(k) 397 399 ENDIF 398 400 ENDDO … … 423 425 424 426 tend(k,j,i) = tend(k,j,i) & 425 & + ( kmzp * ( w(k,j,i) - w(k,j,i-1) ) * ddx&426 & ) * ddzw(k)&427 & + ( kmzp * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1)&428 & + usws(j,i)&429 & ) * ddzw(k) 427 & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) & 428 & + ( w(k,j,i) - w(k,j,i-1) ) * ddx & 429 & ) * rho_air_zw(k) & 430 & - ( -usws(j,i) ) & 431 & ) * ddzw(k) * drho_air(k) 430 432 ENDDO 431 433 ENDDO … … 449 451 450 452 tend(k,j,i) = tend(k,j,i) & 451 & - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx&452 & ) * ddzw(k)&453 & + ( -uswst(j,i)&454 & - kmzm * ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k)&455 & ) * ddzw(k) 453 & + ( ( -uswst(j,i) ) & 454 & - kmzm * ( ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k) & 455 & + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx & 456 & ) * rho_air_zw(k-1) & 457 & ) * ddzw(k) * drho_air(k) 456 458 ENDDO 457 459 ENDDO … … 471 473 472 474 USE arrays_3d, & 473 ONLY: ddzu, ddzw, km, tend, u, usws, uswst, v, w 475 ONLY: ddzu, ddzw, km, tend, u, usws, uswst, v, w, & 476 drho_air, rho_air_zw 474 477 475 478 USE control_parameters, & … … 559 562 & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) & 560 563 & + ( w(k,j,i) - w(k,j,i-1) ) * ddx & 561 & ) 564 & ) * rho_air_zw(k) & 562 565 & - kmzm * ( ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k) & 563 566 & + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx & 564 & ) 565 & ) * ddzw(k) 567 & ) * rho_air_zw(k-1) & 568 & ) * ddzw(k) * drho_air(k) 566 569 ENDDO 567 570 … … 583 586 584 587 tend(k,j,i) = tend(k,j,i) & 585 & + ( kmzp * ( w(k,j,i) - w(k,j,i-1) ) * ddx&586 & ) * ddzw(k)&587 & + ( kmzp * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1)&588 & + usws(j,i)&589 & ) * ddzw(k) 588 & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) & 589 & + ( w(k,j,i) - w(k,j,i-1) ) * ddx & 590 & ) * rho_air_zw(k) & 591 & - ( -usws(j,i) ) & 592 & ) * ddzw(k) * drho_air(k) 590 593 ENDIF 591 594 … … 597 600 ! 598 601 !-- Interpolate eddy diffusivities on staggered gridpoints 599 kmzm = 0.25_wp * & 600 ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) ) 602 kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) ) 601 603 602 604 tend(k,j,i) = tend(k,j,i) & 603 & - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx&604 & ) * ddzw(k)&605 & + ( -uswst(j,i)&606 & - kmzm * ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k)&607 & ) * ddzw(k) 605 & + ( ( -uswst(j,i) ) & 606 & - kmzm * ( ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k) & 607 & + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx & 608 & ) * rho_air_zw(k-1) & 609 & ) * ddzw(k) * drho_air(k) 608 610 ENDIF 609 611 -
palm/trunk/SOURCE/diffusion_v.f90
r2001 r2037 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Anelastic approximation implemented 23 23 ! 24 24 ! Former revisions: … … 111 111 112 112 USE arrays_3d, & 113 ONLY: ddzu, ddzw, km, tend, u, v, vsws, vswst, w 113 ONLY: ddzu, ddzw, km, tend, u, v, vsws, vswst, w, & 114 drho_air, rho_air_zw 114 115 115 116 USE control_parameters, & … … 212 213 & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 213 214 & + ( w(k,j,i) - w(k,j-1,i) ) * ddy & 214 & ) 215 & ) * rho_air_zw(k) & 215 216 & - kmzm * ( ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) & 216 217 & + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy & 217 & ) 218 & ) * ddzw(k) 218 & ) * rho_air_zw(k-1) & 219 & ) * ddzw(k) * drho_air(k) 219 220 ENDDO 220 221 … … 237 238 238 239 tend(k,j,i) = tend(k,j,i) & 239 & + ( kmzp * ( w(k,j,i) - w(k,j-1,i) ) * ddy&240 & ) * ddzw(k)&241 & + ( kmzp * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)&242 & + vsws(j,i)&243 & ) * ddzw(k) 240 & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 241 & + ( w(k,j,i) - w(k,j-1,i) ) * ddy & 242 & ) * rho_air_zw(k) & 243 & - ( -vsws(j,i) ) & 244 & ) * ddzw(k) * drho_air(k) 244 245 ENDIF 245 246 … … 255 256 256 257 tend(k,j,i) = tend(k,j,i) & 257 & - ( kmzm * ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy&258 & ) * ddzw(k)&259 & + ( -vswst(j,i)&260 & - kmzm * ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k)&261 & ) * ddzw(k) 258 & + ( ( -vswst(j,i) ) & 259 & - kmzm * ( ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) & 260 & + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy & 261 & ) * rho_air_zw(k-1) & 262 & ) * ddzw(k) * drho_air(k) 262 263 ENDIF 263 264 … … 276 277 277 278 USE arrays_3d, & 278 ONLY: ddzu, ddzw, km, tend, u, v, vsws, vswst, w 279 ONLY: ddzu, ddzw, km, tend, u, v, vsws, vswst, w, & 280 drho_air, rho_air_zw 279 281 280 282 USE control_parameters, & … … 385 387 & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)& 386 388 & + ( w(k,j,i) - w(k,j-1,i) ) * ddy & 387 & ) 389 & ) * rho_air_zw(k) & 388 390 & - kmzm * ( ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k)& 389 391 & + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy & 390 & ) 391 & ) * ddzw(k) 392 & ) * rho_air_zw(k-1) & 393 & ) * ddzw(k) * drho_air(k) 392 394 ENDIF 393 395 ENDDO … … 418 420 419 421 tend(k,j,i) = tend(k,j,i) & 420 & + ( kmzp * ( w(k,j,i) - w(k,j-1,i) ) * ddy&421 & ) * ddzw(k)&422 & + ( kmzp * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)&423 & + vsws(j,i)&424 & ) * ddzw(k) 422 & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 423 & + ( w(k,j,i) - w(k,j-1,i) ) * ddy & 424 & ) * rho_air_zw(k) & 425 & - ( -vsws(j,i) ) & 426 & ) * ddzw(k) * drho_air(k) 425 427 ENDDO 426 428 ENDDO … … 444 446 445 447 tend(k,j,i) = tend(k,j,i) & 446 & - ( kmzm * ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy&447 & ) * ddzw(k)&448 & + ( -vswst(j,i)&449 & - kmzm * ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k)&450 & ) * ddzw(k) 448 & + ( ( -vswst(j,i) ) & 449 & - kmzm * ( ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) & 450 & + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy & 451 & ) * rho_air_zw(k-1) & 452 & ) * ddzw(k) * drho_air(k) 451 453 ENDDO 452 454 ENDDO … … 466 468 467 469 USE arrays_3d, & 468 ONLY: ddzu, ddzw, km, tend, u, v, vsws, vswst, w 470 ONLY: ddzu, ddzw, km, tend, u, v, vsws, vswst, w, & 471 drho_air, rho_air_zw 469 472 470 473 USE control_parameters, & … … 556 559 & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 557 560 & + ( w(k,j,i) - w(k,j-1,i) ) * ddy & 558 & ) 561 & ) * rho_air_zw(k) & 559 562 & - kmzm * ( ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) & 560 563 & + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy & 561 & ) 562 & ) * ddzw(k) 564 & ) * rho_air_zw(k-1) & 565 & ) * ddzw(k) * drho_air(k) 563 566 ENDDO 564 567 … … 580 583 581 584 tend(k,j,i) = tend(k,j,i) & 582 & + ( kmzp * ( w(k,j,i) - w(k,j-1,i) ) * ddy&583 & ) * ddzw(k)&584 & + ( kmzp * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)&585 & + vsws(j,i)&586 & ) * ddzw(k) 585 & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 586 & + ( w(k,j,i) - w(k,j-1,i) ) * ddy & 587 & ) * rho_air_zw(k) & 588 & - ( -vsws(j,i) ) & 589 & ) * ddzw(k) * drho_air(k) 587 590 ENDIF 588 591 … … 594 597 ! 595 598 !-- Interpolate eddy diffusivities on staggered gridpoints 596 kmzm = 0.25_wp * & 597 ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) ) 599 kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) ) 598 600 599 601 tend(k,j,i) = tend(k,j,i) & 600 & - ( kmzm * ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy&601 & ) * ddzw(k)&602 & + ( -vswst(j,i)&603 & - kmzm * ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k)&604 & ) * ddzw(k) 602 & + ( ( -vswst(j,i) ) & 603 & - kmzm * ( ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) & 604 & + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy & 605 & ) * rho_air_zw(k-1) & 606 & ) * ddzw(k) * drho_air(k) 605 607 ENDIF 606 608 -
palm/trunk/SOURCE/diffusion_w.f90
r2001 r2037 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Anelastic approximation implemented 23 23 ! 24 24 ! Former revisions: … … 119 119 120 120 USE arrays_3d, & 121 ONLY : ddzu, ddzw, km, tend, u, v, w 121 ONLY : ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air 122 122 123 123 USE control_parameters, & … … 184 184 & + 2.0_wp * ( & 185 185 & km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) & 186 & * rho_air(k+1) & 186 187 & - km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) & 187 & ) * ddzu(k+1) 188 & * rho_air(k) & 189 & ) * ddzu(k+1) * drho_air_zw(k) 188 190 ENDDO 189 191 … … 227 229 + 2.0_wp * ( & 228 230 km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) & 231 * rho_air(k+1) & 229 232 - km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) & 230 ) * ddzu(k+1) 233 * rho_air(k) & 234 ) * ddzu(k+1) * drho_air_zw(k) 231 235 ENDDO 232 236 ENDIF … … 246 250 247 251 USE arrays_3d, & 248 ONLY : ddzu, ddzw, km, tend, u, v, w 252 ONLY : ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air 249 253 250 254 USE control_parameters, & … … 316 320 & + 2.0_wp * ( & 317 321 & km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) & 322 & * rho_air(k+1) & 318 323 & - km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) & 319 & ) * ddzu(k+1) 324 & * rho_air(k) & 325 & ) * ddzu(k+1) * drho_air_zw(k) 320 326 ENDIF 321 327 ENDDO … … 361 367 + 2.0_wp * ( & 362 368 km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) & 369 * rho_air(k+1) & 363 370 - km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) & 364 ) * ddzu(k+1) 371 * rho_air(k) & 372 ) * ddzu(k+1) * drho_air_zw(k) 365 373 ENDIF 366 374 ENDDO … … 381 389 382 390 USE arrays_3d, & 383 ONLY : ddzu, ddzw, km, tend, u, v, w 391 ONLY : ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air 384 392 385 393 USE control_parameters, & … … 430 438 & + 2.0_wp * ( & 431 439 & km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) & 440 & * rho_air(k+1) & 432 441 & - km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) & 433 & ) * ddzu(k+1) 442 & * rho_air(k) & 443 & ) * ddzu(k+1) * drho_air_zw(k) 434 444 ENDDO 435 445 … … 485 495 + 2.0_wp * ( & 486 496 km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) & 497 * rho_air(k+1) & 487 498 - km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) & 488 ) * ddzu(k+1) 499 * rho_air(k) & 500 ) * ddzu(k+1) * drho_air_zw(k) 489 501 ENDDO 490 502 ENDIF -
palm/trunk/SOURCE/disturb_heatflux.f90
r2001 r2037 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Anelastic approximation implemented 23 23 ! 24 24 ! Former revisions: … … 65 65 66 66 USE arrays_3d, & 67 ONLY: shf 67 ONLY: shf, heatflux_input_conversion 68 68 69 69 USE control_parameters, & … … 76 76 77 77 USE indices, & 78 ONLY: nx, nxl, nxr, ny, nyn, nys, nzb _s_inner78 ONLY: nx, nxl, nxr, ny, nyn, nys, nzb, nzb_s_inner 79 79 80 80 IMPLICIT NONE … … 97 97 THEN 98 98 IF ( nzb_s_inner(j,i) == 0 ) THEN 99 shf(j,i) = randomnumber * surface_heatflux 99 shf(j,i) = randomnumber * surface_heatflux & 100 * heatflux_input_conversion(nzb) 100 101 ELSE 101 102 ! 102 103 !-- Over topography surface_heatflux is replaced by wall_heatflux(0) 103 shf(j,i) = randomnumber * wall_heatflux(0) 104 shf(j,i) = randomnumber * wall_heatflux(0) & 105 * heatflux_input_conversion(nzb_s_inner(j,i)) 104 106 ENDIF 105 107 ENDIF -
palm/trunk/SOURCE/flow_statistics.f90
r2032 r2037 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Anelastic approximation implemented 23 23 ! 24 24 ! Former revisions: … … 214 214 215 215 USE arrays_3d, & 216 ONLY: ddzu, ddzw, e, hyp, km, kh, nr, ol, p, prho, prr, pt, q, qc, ql,& 217 qr, qs, qsws, qswst, rho_ocean, s, sa, ss, ssws, sswst, saswsb, & 218 saswst, shf, td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q, & 219 time_vert, ts, tswst, u, ug, us, usws, uswst, vsws, v, vg, vpt, & 220 vswst, w, w_subs, zw 216 ONLY: ddzu, ddzw, e, heatflux_output_conversion, hyp, km, kh, & 217 momentumflux_output_conversion, nr, ol, p, prho, prr, pt, q, & 218 qc, ql, qr, qs, qsws, qswst, rho_air, rho_air_zw, rho_ocean, s, & 219 sa, ss, ssws, sswst, saswsb, saswst, shf, td_lsa_lpt, td_lsa_q, & 220 td_sub_lpt, td_sub_q, time_vert, ts, tswst, u, ug, us, usws, & 221 uswst, vsws, v, vg, vpt, vswst, w, w_subs, & 222 waterflux_output_conversion, zw 221 223 222 224 USE cloud_parameters, & … … 326 328 sums_l(:,11,:) = sums_l_l(:,sr,:) ! mixing length from diffusivities 327 329 !-- WARNING: next line still has to be adjusted for OpenMP 328 sums_l(:,21,0) = sums_wsts_bc_l(:,sr) ! heat flux from advec_s_bc 330 sums_l(:,21,0) = sums_wsts_bc_l(:,sr) * & 331 heatflux_output_conversion ! heat flux from advec_s_bc 329 332 sums_l(nzb+9,pr_palm,0) = sums_divold_l(sr) ! old divergence from pres 330 333 sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr) ! new divergence from pres … … 358 361 ! 359 362 !-- Swap the turbulent quantities evaluated in advec_ws. 360 sums_l(:,13,i) = sums_wsus_ws_l(:,i) ! w*u* 361 sums_l(:,15,i) = sums_wsvs_ws_l(:,i) ! w*v* 363 sums_l(:,13,i) = sums_wsus_ws_l(:,i) & 364 * momentumflux_output_conversion ! w*u* 365 sums_l(:,15,i) = sums_wsvs_ws_l(:,i) & 366 * momentumflux_output_conversion ! w*v* 362 367 sums_l(:,30,i) = sums_us2_ws_l(:,i) ! u*2 363 368 sums_l(:,31,i) = sums_vs2_ws_l(:,i) ! v*2 … … 373 378 374 379 DO i = 0, threads_per_task-1 375 sums_l(:,17,i) = sums_wspts_ws_l(:,i) ! w*pt* 380 sums_l(:,17,i) = sums_wspts_ws_l(:,i) & 381 * heatflux_output_conversion ! w*pt* 376 382 IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! w*sa* 377 IF ( humidity ) sums_l(:,49,i) = sums_wsqs_ws_l(:,i) ! w*q* 383 IF ( humidity ) sums_l(:,49,i) = sums_wsqs_ws_l(:,i) & 384 * waterflux_output_conversion ! w*q* 378 385 IF ( passive_scalar ) sums_l(:,116,i) = sums_wsss_ws_l(:,i) ! w*s* 379 386 ENDDO … … 722 729 ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) & 723 730 + ( w(k,j,i) - w(k,j,i-1) ) * ddx & 724 ) * rmask(j,i,sr) 731 ) * rmask(j,i,sr) & 732 * rho_air_zw(k) & 733 * momentumflux_output_conversion(k) 725 734 ! 726 735 !-- Momentum flux w"v" … … 730 739 ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 731 740 + ( w(k,j,i) - w(k,j-1,i) ) * ddy & 732 ) * rmask(j,i,sr) 741 ) * rmask(j,i,sr) & 742 * rho_air_zw(k) & 743 * momentumflux_output_conversion(k) 733 744 ! 734 745 !-- Heat flux w"pt" … … 736 747 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )& 737 748 * ( pt(k+1,j,i) - pt(k,j,i) ) & 749 * rho_air_zw(k) & 750 * heatflux_output_conversion(k) & 738 751 * ddzu(k+1) * rmask(j,i,sr) 739 752 … … 754 767 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )& 755 768 * ( vpt(k+1,j,i) - vpt(k,j,i) ) & 769 * rho_air_zw(k) & 770 * heatflux_output_conversion(k) & 756 771 * ddzu(k+1) * rmask(j,i,sr) 757 772 sums_l(k,48,tn) = sums_l(k,48,tn) & 758 773 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )& 759 774 * ( q(k+1,j,i) - q(k,j,i) ) & 775 * rho_air_zw(k) & 776 * waterflux_output_conversion(k)& 760 777 * ddzu(k+1) * rmask(j,i,sr) 761 778 … … 765 782 * ( ( q(k+1,j,i) - ql(k+1,j,i) )& 766 783 - ( q(k,j,i) - ql(k,j,i) ) ) & 784 * rho_air_zw(k) & 785 * waterflux_output_conversion(k)& 767 786 * ddzu(k+1) * rmask(j,i,sr) 768 787 ENDIF … … 784 803 IF ( use_surface_fluxes ) THEN 785 804 sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + & 805 momentumflux_output_conversion(nzb) * & 786 806 usws(j,i) * rmask(j,i,sr) ! w"u" 787 807 sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + & 808 momentumflux_output_conversion(nzb) * & 788 809 vsws(j,i) * rmask(j,i,sr) ! w"v" 789 810 sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + & 811 heatflux_output_conversion(nzb) * & 790 812 shf(j,i) * rmask(j,i,sr) ! w"pt" 791 813 sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + & … … 799 821 IF ( humidity ) THEN 800 822 sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + & 823 waterflux_output_conversion(nzb) * & 801 824 qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv") 802 825 sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + ( & 803 826 ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) * & 804 827 shf(j,i) + 0.61_wp * pt(nzb,j,i) * & 805 qsws(j,i) ) 828 qsws(j,i) ) & 829 * heatflux_output_conversion(nzb) 806 830 IF ( cloud_droplets ) THEN 807 831 sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + ( & 808 832 ( 1.0_wp + 0.61_wp * q(nzb,j,i) - & 809 833 ql(nzb,j,i) ) * shf(j,i) + & 810 0.61_wp * pt(nzb,j,i) * qsws(j,i) ) 834 0.61_wp * pt(nzb,j,i) * qsws(j,i) ) & 835 * heatflux_output_conversion(nzb) 811 836 ENDIF 812 837 IF ( cloud_physics ) THEN 813 838 ! 814 839 !-- Formula does not work if ql(nzb) /= 0.0 815 sums_l(nzb,51,tn) = sums_l(nzb,51,tn) + & 840 sums_l(nzb,51,tn) = sums_l(nzb,51,tn) + & 841 waterflux_output_conversion(nzb) * & 816 842 qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv") 817 843 ENDIF … … 860 886 IF ( use_top_fluxes ) THEN 861 887 sums_l(nzt:nzt+1,12,tn) = sums_l(nzt:nzt+1,12,tn) + & 888 momentumflux_output_conversion(nzt:nzt+1) * & 862 889 uswst(j,i) * rmask(j,i,sr) ! w"u" 863 890 sums_l(nzt:nzt+1,14,tn) = sums_l(nzt:nzt+1,14,tn) + & 891 momentumflux_output_conversion(nzt:nzt+1) * & 864 892 vswst(j,i) * rmask(j,i,sr) ! w"v" 865 893 sums_l(nzt:nzt+1,16,tn) = sums_l(nzt:nzt+1,16,tn) + & 894 heatflux_output_conversion(nzt:nzt+1) * & 866 895 tswst(j,i) * rmask(j,i,sr) ! w"pt" 867 896 sums_l(nzt:nzt+1,58,tn) = sums_l(nzt:nzt+1,58,tn) + & … … 876 905 IF ( humidity ) THEN 877 906 sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + & 907 waterflux_output_conversion(nzt) * & 878 908 qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv") 879 909 sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + ( & 880 910 ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) * & 881 911 tswst(j,i) + 0.61_wp * pt(nzt,j,i) * & 882 qswst(j,i) ) 912 qswst(j,i) ) & 913 * heatflux_output_conversion(nzt) 883 914 IF ( cloud_droplets ) THEN 884 915 sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + ( & 885 916 ( 1.0_wp + 0.61_wp * q(nzt,j,i) - & 886 917 ql(nzt,j,i) ) * tswst(j,i) + & 887 0.61_wp * pt(nzt,j,i) * qswst(j,i) ) 918 0.61_wp * pt(nzt,j,i) * qswst(j,i) )& 919 * heatflux_output_conversion(nzt) 888 920 ENDIF 889 921 IF ( cloud_physics ) THEN … … 891 923 !-- Formula does not work if ql(nzb) /= 0.0 892 924 sums_l(nzt,51,tn) = sums_l(nzt,51,tn) + & ! w"q" (w"qv") 925 waterflux_output_conversion(nzt) * & 893 926 qswst(j,i) * rmask(j,i,sr) 894 927 ENDIF … … 943 976 vpt(k+1,j,i) - hom(k+1,1,44,sr) ) 944 977 sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * & 978 heatflux_output_conversion(k) * & 945 979 rmask(j,i,sr) 946 980 sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * rmask(j,i,sr) … … 953 987 hom(k+1,1,42,sr) ) 954 988 sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) * & 989 waterflux_output_conversion(k) * & 955 990 rmask(j,i,sr) 956 991 sums_l(k,75,tn) = sums_l(k,75,tn) + qc(k,j,i) * & … … 971 1006 vpt(k+1,j,i) - hom(k+1,1,44,sr) ) 972 1007 sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * & 1008 heatflux_output_conversion(k) * & 973 1009 rmask(j,i,sr) 974 1010 ELSE IF ( ws_scheme_sca .AND. sr == 0 ) THEN 975 sums_l(k,46,tn) = ( 1.0_wp + 0.61_wp * & 976 hom(k,1,41,sr) ) * & 977 sums_l(k,17,tn) + & 978 0.61_wp * hom(k,1,4,sr) * & 979 sums_l(k,49,tn) 1011 sums_l(k,46,tn) = ( ( 1.0_wp + 0.61_wp * & 1012 hom(k,1,41,sr) ) * & 1013 sums_l(k,17,tn) + & 1014 0.61_wp * hom(k,1,4,sr) * & 1015 sums_l(k,49,tn) & 1016 ) * heatflux_output_conversion(k) 980 1017 END IF 981 1018 END IF … … 996 1033 sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5_wp * & 997 1034 ( ust**2 + vst**2 + w(k,j,i)**2 ) & 1035 * momentumflux_output_conversion(k) & 998 1036 * rmask(j,i,sr) 999 1037 ENDDO … … 1017 1055 sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5_wp * & 1018 1056 ( w(k,j,i-1) + w(k,j,i) ) & 1057 * momentumflux_output_conversion(k) & 1019 1058 * ust * rmask(j,i,sr) 1020 1059 ! … … 1022 1061 sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5_wp * & 1023 1062 ( w(k,j-1,i) + w(k,j,i) ) & 1063 * momentumflux_output_conversion(k) & 1024 1064 * vst * rmask(j,i,sr) 1025 1065 ENDDO … … 1038 1078 ( pt(k,j,i) - hom(k,1,4,sr) + & 1039 1079 pt(k+1,j,i) - hom(k+1,1,4,sr) ) & 1080 * heatflux_output_conversion(k) & 1040 1081 * w(k,j,i) * rmask(j,i,sr) 1041 1082 IF ( humidity ) THEN … … 1043 1084 q(k+1,j,i) - hom(k+1,1,41,sr) ) 1044 1085 sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * & 1086 waterflux_output_conversion(k) * & 1045 1087 rmask(j,i,sr) 1046 1088 ENDIF … … 1147 1189 ( kh(k,j,i) + kh(k,j,i-1) ) & 1148 1190 * ( pt(k,j,i-1) - pt(k,j,i) ) & 1191 * rho_air_zw(k) & 1192 * heatflux_output_conversion(k) & 1149 1193 * ddx * rmask(j,i,sr) 1150 1194 sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp * & 1151 1195 ( kh(k,j,i) + kh(k,j-1,i) ) & 1152 1196 * ( pt(k,j-1,i) - pt(k,j,i) ) & 1197 * rho_air_zw(k) & 1198 * heatflux_output_conversion(k) & 1153 1199 * ddy * rmask(j,i,sr) 1154 1200 ! … … 1157 1203 ( u(k,j,i) - hom(k,1,1,sr) ) & 1158 1204 * 0.5_wp * ( pt(k,j,i-1) - hom(k,1,4,sr) + & 1159 pt(k,j,i) - hom(k,1,4,sr) ) 1205 pt(k,j,i) - hom(k,1,4,sr) ) & 1206 * heatflux_output_conversion(k) 1160 1207 pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) + & 1161 1208 pt(k,j,i) - hom(k,1,4,sr) ) … … 1163 1210 ( v(k,j,i) - hom(k,1,2,sr) ) & 1164 1211 * 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) + & 1165 pt(k,j,i) - hom(k,1,4,sr) ) 1212 pt(k,j,i) - hom(k,1,4,sr) ) & 1213 * heatflux_output_conversion(k) 1166 1214 ENDDO 1167 1215 ENDDO … … 1489 1537 ENDIF 1490 1538 1539 hom(:,1,121,sr) = rho_air ! rho_air in Kg/m^3 1540 hom(:,1,122,sr) = rho_air_zw ! rho_air_zw in Kg/m^3 1541 1491 1542 hom(:,1,pr_palm,sr) = sums(:,pr_palm) 1492 1543 ! u*, w'u', w'v', t* (in last profile) … … 1597 1648 THEN 1598 1649 hom(nzb+8,1,pr_palm,sr) = & 1599 ( g / hom(k_surface_level+1,1,4,sr) * hom(k_surface_level,1,18,sr)& 1650 ( g / hom(k_surface_level+1,1,4,sr) * & 1651 ( hom(k_surface_level,1,18,sr) / heatflux_output_conversion(nzb) )& 1600 1652 * ABS( z_i(1) - mean_surface_level_height(sr) ) )**0.333333333_wp 1601 1653 ELSE … … 1708 1760 1709 1761 USE arrays_3d, & 1710 ONLY: ddzu, ddzw, e, hyp, km, kh, nr, p, prho, pt, q, qc, ql, qr, qs, & 1711 qsws, qswst, rho_ocean, s, sa, saswsb, saswst, shf, ss, ssws, sswst, & 1712 td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q, time_vert, ts, & 1713 tswst, u, ug, us, usws, uswst, vsws, v, vg, vpt, vswst, w, & 1714 w_subs, zw 1762 ONLY: ddzu, ddzw, e, heatflux_output_conversion, hyp, km, kh, & 1763 momentumflux_output_conversion, nr, p, prho, pt, q, qc, ql, qr, & 1764 qs, qsws, qswst, rho_air, rho_air_zw, rho_ocean, s, sa, saswsb, & 1765 saswst, shf, ss, ssws, sswst, td_lsa_lpt, td_lsa_q, td_sub_lpt, & 1766 td_sub_q, time_vert, ts, tswst, u, ug, us, usws, uswst, vsws, & 1767 v, vg, vpt, vswst, w, w_subs, waterflux_output_conversion, zw 1715 1768 1716 1769 … … 1828 1881 sums_l(:,11,:) = sums_l_l(:,sr,:) ! mixing length from diffusivities 1829 1882 !-- WARNING: next line still has to be adjusted for OpenMP 1830 sums_l(:,21,0) = sums_wsts_bc_l(:,sr) ! heat flux from advec_s_bc 1883 sums_l(:,21,0) = sums_wsts_bc_l(:,sr) * & 1884 heatflux_output_conversion ! heat flux from advec_s_bc 1831 1885 sums_l(nzb+9,pr_palm,0) = sums_divold_l(sr) ! old divergence from pres 1832 1886 sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr) ! new divergence from pres … … 1860 1914 ! 1861 1915 !-- Swap the turbulent quantities evaluated in advec_ws. 1862 sums_l(:,13,i) = sums_wsus_ws_l(:,i) ! w*u* 1863 sums_l(:,15,i) = sums_wsvs_ws_l(:,i) ! w*v* 1916 sums_l(:,13,i) = sums_wsus_ws_l(:,i) & 1917 * momentumflux_output_conversion ! w*u* 1918 sums_l(:,15,i) = sums_wsvs_ws_l(:,i) & 1919 * momentumflux_output_conversion ! w*v* 1864 1920 sums_l(:,30,i) = sums_us2_ws_l(:,i) ! u*2 1865 1921 sums_l(:,31,i) = sums_vs2_ws_l(:,i) ! v*2 … … 1881 1937 1882 1938 DO i = 0, threads_per_task-1 1883 sums_l(:,17,i) = sums_wspts_ws_l(:,i) ! w*pt* from advec_s_ws 1939 sums_l(:,17,i) = sums_wspts_ws_l(:,i) & 1940 * heatflux_output_conversion ! w*pt* from advec_s_ws 1884 1941 IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! w*sa* 1885 IF ( humidity ) sums_l(:,49,i) = sums_wsqs_ws_l(:,i) !w*q* 1942 IF ( humidity ) sums_l(:,49,i) = sums_wsqs_ws_l(:,i) & 1943 * waterflux_output_conversion !w*q* 1886 1944 IF ( passive_scalar ) sums_l(:,116,i) = sums_wsss_ws_l(:,i) !w*s* 1887 1945 ENDDO … … 2378 2436 + ( w(k,j,i) - w(k,j,i-1) ) * ddx & 2379 2437 ) & 2380 * rmask(j,i,sr) * rflags_invers(j,i,k+1) 2438 * rmask(j,i,sr) * rflags_invers(j,i,k+1) & 2439 * rho_air_zw(k) & 2440 * momentumflux_output_conversion(k) 2381 2441 ! 2382 2442 !-- Momentum flux w"v" … … 2387 2447 + ( w(k,j,i) - w(k,j-1,i) ) * ddy & 2388 2448 ) & 2389 * rmask(j,i,sr) * rflags_invers(j,i,k+1) 2449 * rmask(j,i,sr) * rflags_invers(j,i,k+1) & 2450 * rho_air_zw(k) & 2451 * momentumflux_output_conversion(k) 2390 2452 ! 2391 2453 !-- Heat flux w"pt" 2392 2454 s3 = s3 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) ) & 2393 2455 * ( pt(k+1,j,i) - pt(k,j,i) ) & 2456 * rho_air_zw(k) & 2457 * heatflux_output_conversion(k) & 2394 2458 * ddzu(k+1) * rmask(j,i,sr) & 2395 2459 * rflags_invers(j,i,k+1) … … 2435 2499 s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) ) & 2436 2500 * ( vpt(k+1,j,i) - vpt(k,j,i) ) & 2501 * rho_air_zw(k) & 2502 * heatflux_output_conversion(k) & 2437 2503 * ddzu(k+1) * rmask(j,i,sr) & 2438 2504 * rflags_invers(j,i,k+1) 2439 2505 s2 = s2 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) ) & 2440 2506 * ( q(k+1,j,i) - q(k,j,i) ) & 2507 * rho_air_zw(k) & 2508 * waterflux_output_conversion(k) & 2441 2509 * ddzu(k+1) * rmask(j,i,sr) & 2442 2510 * rflags_invers(j,i,k+1) … … 2459 2527 * ( ( q(k+1,j,i) - ql(k+1,j,i) ) & 2460 2528 - ( q(k,j,i) - ql(k,j,i) ) ) & 2529 * rho_air_zw(k) & 2530 * waterflux_output_conversion(k) & 2461 2531 * ddzu(k+1) * rmask(j,i,sr) & 2462 2532 * rflags_invers(j,i,k+1) … … 2506 2576 ! 2507 2577 !-- Subgridscale fluxes in the Prandtl layer 2508 s1 = s1 + usws(j,i) * rmask(j,i,sr) ! w"u" 2509 s2 = s2 + vsws(j,i) * rmask(j,i,sr) ! w"v" 2510 s3 = s3 + shf(j,i) * rmask(j,i,sr) ! w"pt" 2578 s1 = s1 + usws(j,i) * momentumflux_output_conversion(nzb) & 2579 * rmask(j,i,sr) ! w"u" 2580 s2 = s2 + vsws(j,i) * momentumflux_output_conversion(nzb) & 2581 * rmask(j,i,sr) ! w"v" 2582 s3 = s3 + shf(j,i) * heatflux_output_conversion(nzb) & 2583 * rmask(j,i,sr) ! w"pt" 2511 2584 s4 = s4 + 0.0_wp * rmask(j,i,sr) ! u"pt" 2512 2585 s5 = s5 + 0.0_wp * rmask(j,i,sr) ! v"pt" … … 2545 2618 DO i = nxl, nxr 2546 2619 DO j = nys, nyn 2547 s1 = s1 + qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv") 2620 s1 = s1 + qsws(j,i) * waterflux_output_conversion(nzb) & 2621 * rmask(j,i,sr) ! w"q" (w"qv") 2548 2622 s2 = s2 + ( ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) * shf(j,i) & 2549 + 0.61_wp * pt(nzb,j,i) * qsws(j,i) ) 2623 + 0.61_wp * pt(nzb,j,i) * qsws(j,i) ) & 2624 * heatflux_output_conversion(nzb) 2550 2625 ENDDO 2551 2626 ENDDO … … 2564 2639 s1 = s1 + ( ( 1.0_wp + & 2565 2640 0.61_wp * q(nzb,j,i) - ql(nzb,j,i) ) * & 2566 shf(j,i) + 0.61_wp * pt(nzb,j,i) * qsws(j,i) ) 2641 shf(j,i) + 0.61_wp * pt(nzb,j,i) * qsws(j,i) )& 2642 * heatflux_output_conversion(nzb) 2567 2643 ENDDO 2568 2644 ENDDO … … 2582 2658 ! 2583 2659 !-- Formula does not work if ql(nzb) /= 0.0 2584 s1 = s1 + qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv") 2660 s1 = s1 + qsws(j,i) * waterflux_output_conversion(nzb) & 2661 * rmask(j,i,sr) ! w"q" (w"qv") 2585 2662 ENDDO 2586 2663 ENDDO … … 2624 2701 DO i = nxl, nxr 2625 2702 DO j = nys, nyn 2626 s1 = s1 + uswst(j,i) * rmask(j,i,sr) ! w"u" 2627 s2 = s2 + vswst(j,i) * rmask(j,i,sr) ! w"v" 2628 s3 = s3 + tswst(j,i) * rmask(j,i,sr) ! w"pt" 2703 s1 = s1 + uswst(j,i) * momentumflux_output_conversion(nzt:nzt+1) & 2704 * rmask(j,i,sr) ! w"u" 2705 s2 = s2 + vswst(j,i) * momentumflux_output_conversion(nzt:nzt+1) & 2706 * rmask(j,i,sr) ! w"v" 2707 s3 = s3 + tswst(j,i) * heatflux_output_conversion(nzt:nzt+1) & 2708 * rmask(j,i,sr) ! w"pt" 2629 2709 s4 = s4 + 0.0_wp * rmask(j,i,sr) ! u"pt" 2630 2710 s5 = s5 + 0.0_wp * rmask(j,i,sr) ! v"pt" … … 2663 2743 DO i = nxl, nxr 2664 2744 DO j = nys, nyn 2665 s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv") 2745 s1 = s1 + qswst(j,i) * waterflux_output_conversion(nzt) & 2746 * rmask(j,i,sr) ! w"q" (w"qv") 2666 2747 s2 = s2 + ( ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) * tswst(j,i) +& 2667 0.61_wp * pt(nzt,j,i) * qswst(j,i) ) 2748 0.61_wp * pt(nzt,j,i) * qswst(j,i) ) & 2749 * heatflux_output_conversion(nzt) 2668 2750 ENDDO 2669 2751 ENDDO … … 2683 2765 0.61_wp * q(nzt,j,i) - ql(nzt,j,i) ) * & 2684 2766 tswst(j,i) + & 2685 0.61_wp * pt(nzt,j,i) * qswst(j,i) ) 2767 0.61_wp * pt(nzt,j,i) * qswst(j,i) ) & 2768 * heatflux_output_conversion(nzt) 2686 2769 ENDDO 2687 2770 ENDDO … … 2701 2784 ! 2702 2785 !-- Formula does not work if ql(nzb) /= 0.0 2703 s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv") 2786 s1 = s1 + qswst(j,i) * waterflux_output_conversion(nzt) & 2787 * rmask(j,i,sr) ! w"q" (w"qv") 2704 2788 ENDDO 2705 2789 ENDDO … … 2755 2839 !-- Energy flux w*e* (has to be adjusted?) 2756 2840 s3 = s3 + w(k,j,i) * 0.5_wp * ( ust**2 + vst**2 + w(k,j,i)**2 )& 2757 * rmask(j,i,sr) * rflags_invers(j,i,k+1) 2841 * rmask(j,i,sr) * rflags_invers(j,i,k+1) & 2842 * momentumflux_output_conversion(k) 2758 2843 ENDDO 2759 2844 ENDDO … … 2822 2907 s1 = s1 + 0.5_wp * ( vpt(k,j,i) - hom(k,1,44,sr) + & 2823 2908 vpt(k+1,j,i) - hom(k+1,1,44,sr) ) * & 2909 heatflux_output_conversion(k) * & 2824 2910 w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1) 2825 2911 ENDDO … … 2839 2925 s1 = s1 + 0.5_wp * ( ( q(k,j,i) - ql(k,j,i) ) - hom(k,1,42,sr) + & 2840 2926 ( q(k+1,j,i) - ql(k+1,j,i) ) - hom(k+1,1,42,sr) ) & 2927 * waterflux_output_conversion(k) & 2841 2928 * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1) 2842 2929 ENDDO … … 2928 3015 s1 = s1 + 0.5_wp * ( vpt(k,j,i) - hom(k,1,44,sr) + & 2929 3016 vpt(k+1,j,i) - hom(k+1,1,44,sr) ) & 3017 * heatflux_output_conversion(k) & 2930 3018 * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1) 2931 3019 ENDDO … … 2939 3027 !$acc parallel loop present( hom, sums_l ) 2940 3028 DO k = nzb, nzt_diff 2941 sums_l(k,46,tn) = ( 1.0_wp + 0.61_wp * hom(k,1,41,sr) ) * sums_l(k,17,tn) + & 2942 0.61_wp * hom(k,1,4,sr) * sums_l(k,49,tn) 3029 sums_l(k,46,tn) = ( ( 1.0_wp + 0.61_wp * hom(k,1,41,sr) ) * & 3030 sums_l(k,17,tn) + 0.61_wp * & 3031 hom(k,1,4,sr) * sums_l(k,49,tn) & 3032 ) * heatflux_output_conversion(k) 2943 3033 ENDDO 2944 3034 !$acc end parallel loop … … 2993 3083 s1 = s1 + 0.5_wp * ( w(k,j,i-1) + w(k,j,i) ) & 2994 3084 * ust * rmask(j,i,sr) & 3085 * momentumflux_output_conversion(k) & 2995 3086 * rflags_invers(j,i,k+1) 2996 3087 ! … … 2998 3089 s2 = s2 + 0.5_wp * ( w(k,j-1,i) + w(k,j,i) ) & 2999 3090 * vst * rmask(j,i,sr) & 3091 * momentumflux_output_conversion(k) & 3000 3092 * rflags_invers(j,i,k+1) 3001 3093 ENDDO 3002 3094 ENDDO 3003 3095 sums_l(k,13,tn) = s1 3004 sums_l(k,15,tn) = s 13096 sums_l(k,15,tn) = s2 3005 3097 ENDDO 3006 3098 !$acc end parallel loop … … 3021 3113 s1 = s1 + 0.5_wp * ( pt(k,j,i) - hom(k,1,4,sr) + & 3022 3114 pt(k+1,j,i) - hom(k+1,1,4,sr) ) & 3115 * heatflux_output_conversion(k) & 3023 3116 * w(k,j,i) * rmask(j,i,sr) & 3024 3117 * rflags_invers(j,i,k+1) … … 3039 3132 s1 = s1 + 0.5_wp * ( q(k,j,i) - hom(k,1,41,sr) + & 3040 3133 q(k+1,j,i) - hom(k+1,1,41,sr) ) & 3134 * waterflux_output_conversion(k) & 3041 3135 * w(k,j,i) * rmask(j,i,sr) & 3042 3136 * rflags_invers(j,i,k+1) … … 3167 3261 ( kh(k,j,i) + kh(k,j,i-1) ) & 3168 3262 * ( pt(k,j,i-1) - pt(k,j,i) ) & 3263 * rho_air_zw(k) & 3264 * heatflux_output_conversion(k) & 3169 3265 * ddx * rmask(j,i,sr) 3170 3266 sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp * & 3171 3267 ( kh(k,j,i) + kh(k,j-1,i) ) & 3172 3268 * ( pt(k,j-1,i) - pt(k,j,i) ) & 3269 * rho_air_zw(k) & 3270 * heatflux_output_conversion(k) & 3173 3271 * ddy * rmask(j,i,sr) 3174 3272 ! … … 3177 3275 ( u(k,j,i) - hom(k,1,1,sr) ) * 0.5_wp * & 3178 3276 ( pt(k,j,i-1) - hom(k,1,4,sr) + & 3179 pt(k,j,i) - hom(k,1,4,sr) ) 3277 pt(k,j,i) - hom(k,1,4,sr) ) & 3278 * heatflux_output_conversion(k) 3180 3279 pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) + & 3181 3280 pt(k,j,i) - hom(k,1,4,sr) ) … … 3183 3282 ( v(k,j,i) - hom(k,1,2,sr) ) * 0.5_wp * & 3184 3283 ( pt(k,j-1,i) - hom(k,1,4,sr) + & 3185 pt(k,j,i) - hom(k,1,4,sr) ) 3284 pt(k,j,i) - hom(k,1,4,sr) ) & 3285 * heatflux_output_conversion(k) 3186 3286 ENDDO 3187 3287 ENDDO … … 3468 3568 END IF 3469 3569 3570 hom(:,1,121,sr) = rho_air ! rho_air in Kg/m^3 3571 hom(:,1,122,sr) = rho_air_zw ! rho_air_zw in Kg/m^3 3572 3470 3573 hom(:,1,pr_palm,sr) = sums(:,pr_palm) 3471 3574 ! u*, w'u', w'v', t* (in last profile) … … 3576 3679 IF ( hom(nzb,1,18,sr) > 1.0E-8_wp .AND. z_i(1) /= 0.0_wp ) THEN 3577 3680 hom(nzb+8,1,pr_palm,sr) = & 3578 ( g / hom(k_surface_level+1,1,4,sr) * hom(k_surface_level,1,18,sr)& 3681 ( g / hom(k_surface_level+1,1,4,sr) * & 3682 ( hom(k_surface_level,1,18,sr) / heatflux_output_conversion(nzb) )& 3579 3683 * ABS( z_i(1) - mean_surface_level_height(sr) ) )**0.333333333_wp 3580 3684 ELSE -
palm/trunk/SOURCE/header.f90
r2001 r2037 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Anelastic approximation implemented 23 23 ! 24 24 ! Former revisions: … … 552 552 !-- Numerical schemes 553 553 WRITE ( io, 110 ) 554 WRITE ( io, 121 ) TRIM( approximation ) 554 555 IF ( psolver(1:7) == 'poisfft' ) THEN 555 556 WRITE ( io, 111 ) TRIM( fft_method ) … … 1906 1907 110 FORMAT (/' Numerical Schemes:'/ & 1907 1908 ' -----------------'/) 1909 121 FORMAT (' --> Use the ',A,' approximation for the model equations.') 1908 1910 111 FORMAT (' --> Solve perturbation pressure via FFT using ',A,' routines') 1909 1911 112 FORMAT (' --> Solve perturbation pressure via SOR-Red/Black-Schema'/ & -
palm/trunk/SOURCE/init_3d_model.f90
r2032 r2037 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Anelastic approximation implemented 23 23 ! 24 24 ! Former revisions: … … 315 315 USE arrays_3d 316 316 317 USE cloud_parameters, & 318 ONLY: cp, l_v, r_d 319 317 320 USE constants, & 318 321 ONLY: pi … … 324 327 325 328 USE grid_variables, & 326 ONLY: dx, dy 329 ONLY: dx, dy, ddx2_mg, ddy2_mg 327 330 328 331 USE indices … … 393 396 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: ngp_2dh_outer_l !< 394 397 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: ngp_2dh_s_inner_l !< 398 399 REAL(wp) :: t_surface !< air temperature at the surface 400 401 REAL(wp), DIMENSION(:), ALLOCATABLE :: p_hydrostatic !< hydrostatic pressure 402 403 INTEGER(iwp) :: l !< loop variable 404 INTEGER(iwp) :: nzt_l !< index of top PE boundary for multigrid level 405 REAL(wp) :: dx_l !< grid spacing along x on different multigrid level 406 REAL(wp) :: dy_l !< grid spacing along y on different multigrid level 395 407 396 408 REAL(wp), DIMENSION(1:3) :: volume_flow_area_l !< … … 647 659 648 660 ! 661 !-- Allocation of anelastic and Boussinesq approximation specific arrays 662 ALLOCATE( p_hydrostatic(nzb:nzt+1) ) 663 ALLOCATE( rho_air(nzb:nzt+1) ) 664 ALLOCATE( rho_air_zw(nzb:nzt+1) ) 665 ALLOCATE( drho_air(nzb:nzt+1) ) 666 ALLOCATE( drho_air_zw(nzb:nzt+1) ) 667 668 ! 669 !-- Density profile calculation for anelastic approximation 670 IF ( TRIM( approximation ) == 'anelastic' ) THEN 671 t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**( r_d / cp ) 672 DO k = nzb, nzt+1 673 p_hydrostatic(k) = surface_pressure * 100.0_wp * & 674 ( 1 - ( g * zu(k) ) / ( cp * t_surface ) & 675 )**( cp / r_d ) 676 rho_air(k) = ( p_hydrostatic(k) * & 677 ( 100000.0_wp / p_hydrostatic(k) & 678 )**( r_d / cp ) & 679 ) / ( r_d * pt_init(k) ) 680 ENDDO 681 DO k = nzb, nzt 682 rho_air_zw(k) = 0.5_wp * ( rho_air(k) + rho_air(k+1) ) 683 ENDDO 684 rho_air_zw(nzt+1) = rho_air_zw(nzt) & 685 + 2.0_wp * ( rho_air(nzt+1) - rho_air_zw(nzt) ) 686 ELSE 687 rho_air = 1.0_wp 688 rho_air_zw = 1.0_wp 689 ENDIF 690 691 !-- compute the inverse density array in order to avoid expencive divisions 692 drho_air = 1.0_wp / rho_air 693 drho_air_zw = 1.0_wp / rho_air_zw 694 695 ! 696 !-- Allocation of flux conversion arrays 697 ALLOCATE( heatflux_input_conversion(nzb:nzt+1) ) 698 ALLOCATE( waterflux_input_conversion(nzb:nzt+1) ) 699 ALLOCATE( momentumflux_input_conversion(nzb:nzt+1) ) 700 ALLOCATE( heatflux_output_conversion(nzb:nzt+1) ) 701 ALLOCATE( waterflux_output_conversion(nzb:nzt+1) ) 702 ALLOCATE( momentumflux_output_conversion(nzb:nzt+1) ) 703 704 ! 705 !-- calculate flux conversion factors according to approximation and in-/output mode 706 DO k = nzb, nzt+1 707 708 IF ( TRIM( flux_input_mode ) == 'kinematic' ) THEN 709 heatflux_input_conversion(k) = rho_air_zw(k) 710 waterflux_input_conversion(k) = rho_air_zw(k) 711 momentumflux_input_conversion(k) = rho_air_zw(k) 712 ELSEIF ( TRIM( flux_input_mode ) == 'dynamic' ) THEN 713 heatflux_input_conversion(k) = 1.0_wp / cp 714 waterflux_input_conversion(k) = 1.0_wp / l_v 715 momentumflux_input_conversion(k) = 1.0_wp 716 ENDIF 717 718 IF ( TRIM( flux_output_mode ) == 'kinematic' ) THEN 719 heatflux_output_conversion(k) = drho_air_zw(k) 720 waterflux_output_conversion(k) = drho_air_zw(k) 721 momentumflux_output_conversion(k) = drho_air_zw(k) 722 ELSEIF ( TRIM( flux_output_mode ) == 'dynamic' ) THEN 723 heatflux_output_conversion(k) = cp 724 waterflux_output_conversion(k) = l_v 725 momentumflux_output_conversion(k) = 1.0_wp 726 ENDIF 727 728 IF ( .NOT. humidity ) THEN 729 waterflux_input_conversion(k) = 1.0_wp 730 waterflux_output_conversion(k) = 1.0_wp 731 ENDIF 732 733 ENDDO 734 735 ! 736 !-- In case of multigrid method, compute grid lengths and grid factors for the 737 !-- grid levels with respective density on each grid 738 IF ( psolver(1:9) == 'multigrid' ) THEN 739 740 ALLOCATE( ddx2_mg(maximum_grid_level) ) 741 ALLOCATE( ddy2_mg(maximum_grid_level) ) 742 ALLOCATE( dzu_mg(nzb+1:nzt+1,maximum_grid_level) ) 743 ALLOCATE( dzw_mg(nzb+1:nzt+1,maximum_grid_level) ) 744 ALLOCATE( f1_mg(nzb+1:nzt,maximum_grid_level) ) 745 ALLOCATE( f2_mg(nzb+1:nzt,maximum_grid_level) ) 746 ALLOCATE( f3_mg(nzb+1:nzt,maximum_grid_level) ) 747 ALLOCATE( rho_air_mg(nzb:nzt+1,maximum_grid_level) ) 748 ALLOCATE( rho_air_zw_mg(nzb:nzt+1,maximum_grid_level) ) 749 750 dzu_mg(:,maximum_grid_level) = dzu 751 rho_air_mg(:,maximum_grid_level) = rho_air 752 ! 753 !-- Next line to ensure an equally spaced grid. 754 dzu_mg(1,maximum_grid_level) = dzu(2) 755 rho_air_mg(nzb,maximum_grid_level) = rho_air(nzb) + & 756 (rho_air(nzb) - rho_air(nzb+1)) 757 758 dzw_mg(:,maximum_grid_level) = dzw 759 rho_air_zw_mg(:,maximum_grid_level) = rho_air_zw 760 nzt_l = nzt 761 DO l = maximum_grid_level-1, 1, -1 762 dzu_mg(nzb+1,l) = 2.0_wp * dzu_mg(nzb+1,l+1) 763 dzw_mg(nzb+1,l) = 2.0_wp * dzw_mg(nzb+1,l+1) 764 rho_air_mg(nzb,l) = rho_air_mg(nzb,l+1) + (rho_air_mg(nzb,l+1) - rho_air_mg(nzb+1,l+1)) 765 rho_air_zw_mg(nzb,l) = rho_air_zw_mg(nzb,l+1) + (rho_air_zw_mg(nzb,l+1) - rho_air_zw_mg(nzb+1,l+1)) 766 rho_air_mg(nzb+1,l) = rho_air_mg(nzb+1,l+1) 767 rho_air_zw_mg(nzb+1,l) = rho_air_zw_mg(nzb+1,l+1) 768 nzt_l = nzt_l / 2 769 DO k = 2, nzt_l+1 770 dzu_mg(k,l) = dzu_mg(2*k-2,l+1) + dzu_mg(2*k-1,l+1) 771 dzw_mg(k,l) = dzw_mg(2*k-2,l+1) + dzw_mg(2*k-1,l+1) 772 rho_air_mg(k,l) = rho_air_mg(2*k-1,l+1) 773 rho_air_zw_mg(k,l) = rho_air_zw_mg(2*k-1,l+1) 774 ENDDO 775 ENDDO 776 777 nzt_l = nzt 778 dx_l = dx 779 dy_l = dy 780 DO l = maximum_grid_level, 1, -1 781 ddx2_mg(l) = 1.0_wp / dx_l**2 782 ddy2_mg(l) = 1.0_wp / dy_l**2 783 DO k = nzb+1, nzt_l 784 f2_mg(k,l) = rho_air_zw_mg(k,l) / ( dzu_mg(k+1,l) * dzw_mg(k,l) ) 785 f3_mg(k,l) = rho_air_zw_mg(k-1,l) / ( dzu_mg(k,l) * dzw_mg(k,l) ) 786 f1_mg(k,l) = 2.0_wp * ( ddx2_mg(l) + ddy2_mg(l) ) & 787 * rho_air_mg(k,l) + f2_mg(k,l) + f3_mg(k,l) 788 ENDDO 789 nzt_l = nzt_l / 2 790 dx_l = dx_l * 2.0_wp 791 dy_l = dy_l * 2.0_wp 792 ENDDO 793 794 ENDIF 795 796 ! 649 797 !-- 3D-array for storing the dissipation, needed for calculating the sgs 650 798 !-- particle velocities … … 883 1031 vsws = 0.0_wp 884 1032 ENDIF 885 uswst = top_momentumflux_u 886 vswst = top_momentumflux_v 1033 uswst = top_momentumflux_u * momentumflux_input_conversion(nzt+1) 1034 vswst = top_momentumflux_v * momentumflux_input_conversion(nzt+1) 887 1035 888 1036 ! … … 1043 1191 us = 1E-30_wp 1044 1192 usws = 0.0_wp 1045 uswst = top_momentumflux_u 1193 uswst = top_momentumflux_u * momentumflux_input_conversion(nzt+1) 1046 1194 vsws = 0.0_wp 1047 vswst = top_momentumflux_v 1195 vswst = top_momentumflux_v * momentumflux_input_conversion(nzt+1) 1048 1196 IF ( humidity ) qs = 0.0_wp 1049 1197 IF ( passive_scalar ) ss = 0.0_wp … … 1165 1313 CALL disturb_heatflux 1166 1314 ELSE 1167 shf = surface_heatflux 1315 shf = surface_heatflux * heatflux_input_conversion(nzb) 1168 1316 ! 1169 1317 !-- Initialize shf with data from external file LSF_DATA … … 1178 1326 DO j = nysg, nyng 1179 1327 IF ( nzb_s_inner(j,i) /= 0 ) THEN 1180 shf(j,i) = wall_heatflux(0) 1328 shf(j,i) = wall_heatflux(0) & 1329 * heatflux_input_conversion(nzb_s_inner(j,i)) 1181 1330 ENDIF 1182 1331 ENDDO … … 1194 1343 ENDIF 1195 1344 IF ( constant_waterflux ) THEN 1196 qsws = surface_waterflux 1345 qsws = surface_waterflux * waterflux_input_conversion(nzb) 1197 1346 ! 1198 1347 !-- Over topography surface_waterflux is replaced by … … 1203 1352 DO j = nysg, nyng 1204 1353 IF ( nzb_s_inner(j,i) /= 0 ) THEN 1205 qsws(j,i) = wall_qflux(0) 1354 qsws(j,i) = wall_qflux(0) & 1355 * waterflux_input_conversion(nzb_s_inner(j,i)) 1206 1356 ENDIF 1207 1357 ENDDO … … 1241 1391 ! 1242 1392 !-- Prescribe to heat flux 1243 IF ( constant_top_heatflux ) tswst = top_heatflux 1393 IF ( constant_top_heatflux ) tswst = top_heatflux & 1394 * heatflux_input_conversion(nzt+1) 1244 1395 ! 1245 1396 !-- Prescribe zero latent flux at the top -
palm/trunk/SOURCE/init_grid.f90
r2022 r2037 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Anelastic approximation implemented 23 23 ! 24 24 ! Former revisions: … … 222 222 223 223 USE grid_variables, & 224 ONLY: ddx, ddx2, dd x2_mg, ddy, ddy2, ddy2_mg, dx, dx2, dy, dy2, fwxm,&224 ONLY: ddx, ddx2, ddy, ddy2, dx, dx2, dy, dy2, fwxm, & 225 225 fwxp, fwym, fwyp, fxm, fxp, fym, fyp, wall_e_x, wall_e_y, & 226 226 wall_u, wall_v, wall_w_x, wall_w_y, zu_s_inner, zw_w_inner … … 293 293 294 294 REAL(wp) :: dum !< dummy variable to skip columns while reading topography file 295 REAL(wp) :: dx_l !< grid spacing along x on different multigrid level296 REAL(wp) :: dy_l !< grid spacing along y on different multigrid level297 295 REAL(wp) :: dz_stretched !< stretched vertical grid spacing 298 296 … … 431 429 ddzu_pres = ddzu 432 430 ddzu_pres(1) = ddzu_pres(2) ! change for lowest level 433 ENDIF434 435 !436 !-- In case of multigrid method, compute grid lengths and grid factors for the437 !-- grid levels438 IF ( psolver(1:9) == 'multigrid' ) THEN439 440 ALLOCATE( ddx2_mg(maximum_grid_level), ddy2_mg(maximum_grid_level), &441 dzu_mg(nzb+1:nzt+1,maximum_grid_level), &442 dzw_mg(nzb+1:nzt+1,maximum_grid_level), &443 f1_mg(nzb+1:nzt,maximum_grid_level), &444 f2_mg(nzb+1:nzt,maximum_grid_level), &445 f3_mg(nzb+1:nzt,maximum_grid_level) )446 447 dzu_mg(:,maximum_grid_level) = dzu448 !449 !-- Next line to ensure an equally spaced grid.450 dzu_mg(1,maximum_grid_level) = dzu(2)451 452 dzw_mg(:,maximum_grid_level) = dzw453 nzt_l = nzt454 DO l = maximum_grid_level-1, 1, -1455 dzu_mg(nzb+1,l) = 2.0_wp * dzu_mg(nzb+1,l+1)456 dzw_mg(nzb+1,l) = 2.0_wp * dzw_mg(nzb+1,l+1)457 nzt_l = nzt_l / 2458 DO k = 2, nzt_l+1459 dzu_mg(k,l) = dzu_mg(2*k-2,l+1) + dzu_mg(2*k-1,l+1)460 dzw_mg(k,l) = dzw_mg(2*k-2,l+1) + dzw_mg(2*k-1,l+1)461 ENDDO462 ENDDO463 464 nzt_l = nzt465 dx_l = dx466 dy_l = dy467 DO l = maximum_grid_level, 1, -1468 ddx2_mg(l) = 1.0_wp / dx_l**2469 ddy2_mg(l) = 1.0_wp / dy_l**2470 DO k = nzb+1, nzt_l471 f2_mg(k,l) = 1.0_wp / ( dzu_mg(k+1,l) * dzw_mg(k,l) )472 f3_mg(k,l) = 1.0_wp / ( dzu_mg(k,l) * dzw_mg(k,l) )473 f1_mg(k,l) = 2.0_wp * ( ddx2_mg(l) + ddy2_mg(l) ) + &474 f2_mg(k,l) + f3_mg(k,l)475 ENDDO476 nzt_l = nzt_l / 2477 dx_l = dx_l * 2.0_wp478 dy_l = dy_l * 2.0_wp479 ENDDO480 481 431 ENDIF 482 432 -
palm/trunk/SOURCE/ls_forcing_mod.f90
r2001 r2037 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Anelastic approximation implemented 23 23 ! 24 24 ! Former revisions: … … 109 109 ONLY: p_surf, pt_surf, q_surf, qsws_surf, shf_surf, td_lsa_lpt, & 110 110 td_lsa_q, td_sub_lpt, td_sub_q, time_surf, time_vert, & 111 heatflux_input_conversion, waterflux_input_conversion, & 111 112 ug_vert, vg_vert, wsubs_vert, zu 112 113 … … 213 214 ENDDO 214 215 216 shf_surf = shf_surf * heatflux_input_conversion(nzb) 217 qsws_surf = qsws_surf * waterflux_input_conversion(nzb) 215 218 216 219 IF ( time_surf(1) > end_time ) THEN -
palm/trunk/SOURCE/modules.f90
r2032 r2037 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Anelastic approximation implemented 23 23 ! 24 24 ! Former revisions: … … 493 493 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: rif_wall, tri 494 494 495 REAL(wp), DIMENSION(:), ALLOCATABLE :: rho_air !< air density profile on the uv grid 496 REAL(wp), DIMENSION(:), ALLOCATABLE :: rho_air_zw !< air density profile on the w grid 497 REAL(wp), DIMENSION(:), ALLOCATABLE :: drho_air !< inverse air density profile on the uv grid 498 REAL(wp), DIMENSION(:), ALLOCATABLE :: drho_air_zw !< inverse air density profile on the w grid 499 500 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rho_air_mg !< air density profiles on the uv grid for multigrid 501 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rho_air_zw_mg !< air density profiles on the w grid for multigrid 502 503 REAL(wp), DIMENSION(:), ALLOCATABLE :: heatflux_input_conversion !< conversion factor array for heatflux input 504 REAL(wp), DIMENSION(:), ALLOCATABLE :: waterflux_input_conversion !< conversion factor array for waterflux input 505 REAL(wp), DIMENSION(:), ALLOCATABLE :: momentumflux_input_conversion !< conversion factor array for momentumflux input 506 REAL(wp), DIMENSION(:), ALLOCATABLE :: heatflux_output_conversion !< conversion factor array for heatflux output 507 REAL(wp), DIMENSION(:), ALLOCATABLE :: waterflux_output_conversion !< conversion factor array for waterflux output 508 REAL(wp), DIMENSION(:), ALLOCATABLE :: momentumflux_output_conversion !< conversion factor array for momentumflux output 509 495 510 496 511 SAVE … … 616 631 psolver = 'poisfft', & 617 632 scalar_advec = 'ws-scheme' 633 CHARACTER (LEN=20) :: approximation = 'boussinesq' 634 CHARACTER (LEN=40) :: flux_input_mode = 'approximation-specific' 635 CHARACTER (LEN=40) :: flux_output_mode = 'approximation-specific' 618 636 CHARACTER (LEN=20) :: bc_e_b = 'neumann', bc_lr = 'cyclic', & 619 637 bc_ns = 'cyclic', bc_p_b = 'neumann', & -
palm/trunk/SOURCE/netcdf_interface_mod.f90
r2032 r2037 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Anelastic approximation implemented 23 23 ! 24 24 ! Former revisions: … … 252 252 ( 'unknown ', i9 = 1, dots_max-42 ) /) 253 253 254 CHARACTER (LEN=16) :: heatflux_output_unit !< unit for heatflux output 255 CHARACTER (LEN=16) :: waterflux_output_unit !< unit for waterflux output 256 CHARACTER (LEN=16) :: momentumflux_output_unit !< unit for momentumflux output 257 254 258 CHARACTER (LEN=9), DIMENSION(300) :: dopr_unit = 'unknown' 255 259 … … 368 372 id_var_x_fl, id_var_y_fl, id_var_z_fl, nc_stat, & 369 373 netcdf_data_format, netcdf_data_format_string, netcdf_deflate, & 370 netcdf_precision, output_for_t0 371 372 374 netcdf_precision, output_for_t0, heatflux_output_unit, & 375 waterflux_output_unit, momentumflux_output_unit 376 373 377 SAVE 374 378 -
palm/trunk/SOURCE/parin.f90
r2036 r2037 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Anelastic approximation implemented 23 23 ! 24 24 ! Former revisions: … … 313 313 314 314 315 NAMELIST /inipar/ alpha_surface, background_communication, bc_e_b, bc_lr, & 315 NAMELIST /inipar/ alpha_surface, approximation, & 316 background_communication, bc_e_b, bc_lr, & 316 317 bc_ns, bc_p_b, bc_p_t, bc_pt_b, bc_pt_t, bc_q_b, & 317 318 bc_q_t,bc_s_b, bc_s_t, bc_sa_t, bc_uv_b, bc_uv_t, & … … 332 333 dt, dt_pr_1d, dt_run_control_1d, dx, dy, dz, dz_max, & 333 334 dz_stretch_factor, dz_stretch_level, end_time_1d, & 334 ensemble_member_nr, & 335 e_init, e_min, fft_method, galilei_transformation, humidity, & 335 ensemble_member_nr, e_init, e_min, fft_method, & 336 flux_input_mode, flux_output_mode, & 337 galilei_transformation, humidity, & 336 338 inflow_damping_height, inflow_damping_width, & 337 339 inflow_disturbance_begin, inflow_disturbance_end, & -
palm/trunk/SOURCE/poismg_mod.f90
r2022 r2037 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Anelastic approximation implemented (stll error in optimized multigrid) 23 23 ! 24 24 ! Former revisions: … … 262 262 263 263 USE arrays_3d, & 264 ONLY: f1_mg, f2_mg, f3_mg 264 ONLY: f1_mg, f2_mg, f3_mg, rho_air_mg 265 265 266 266 USE control_parameters, & … … 309 309 kp1 = k-ind_even_odd 310 310 r(k,j,i) = f_mg(k,j,i) & 311 - ddx2_mg(l) *&311 - rho_air_mg(k,l) * ddx2_mg(l) * & 312 312 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 313 - ddy2_mg(l) *&313 - rho_air_mg(k,l) * ddy2_mg(l) * & 314 314 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 315 315 - f2_mg_b(k,l) * p_mg(kp1,j,i) & … … 322 322 kp1 = k+ind_even_odd+1 323 323 r(k,j,i) = f_mg(k,j,i) & 324 - ddx2_mg(l) *&324 - rho_air_mg(k,l) * ddx2_mg(l) * & 325 325 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 326 - ddy2_mg(l) *&326 - rho_air_mg(k,l) * ddy2_mg(l) * & 327 327 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 328 328 - f2_mg_b(k,l) * p_mg(kp1,j,i) & … … 693 693 694 694 USE arrays_3d, & 695 ONLY: f1_mg, f2_mg, f3_mg 695 ONLY: f1_mg, f2_mg, f3_mg, rho_air_mg 696 696 697 697 USE control_parameters, & … … 755 755 kp1 = k-ind_even_odd 756 756 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 757 ddx2_mg(l) *&757 rho_air_mg(k,l) * ddx2_mg(l) * & 758 758 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 759 + ddy2_mg(l) *&759 + rho_air_mg(k,l) * ddy2_mg(l) * & 760 760 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 761 761 + f2_mg_b(k,l) * p_mg(kp1,j,i) & … … 774 774 kp1 = k-ind_even_odd 775 775 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 776 ddx2_mg(l) *&776 rho_air_mg(k,l) * ddx2_mg(l) * & 777 777 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 778 + ddy2_mg(l) *&778 + rho_air_mg(k,l) * ddy2_mg(l) * & 779 779 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 780 780 + f2_mg_b(k,l) * p_mg(kp1,j,i) & … … 793 793 kp1 = k+ind_even_odd+1 794 794 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 795 ddx2_mg(l) *&795 rho_air_mg(k,l) * ddx2_mg(l) * & 796 796 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 797 + ddy2_mg(l) *&797 + rho_air_mg(k,l) * ddy2_mg(l) * & 798 798 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 799 799 + f2_mg_b(k,l) * p_mg(kp1,j,i) & … … 812 812 kp1 = k+ind_even_odd+1 813 813 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 814 ddx2_mg(l) *&814 rho_air_mg(k,l) * ddx2_mg(l) * & 815 815 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 816 + ddy2_mg(l) *&816 + rho_air_mg(k,l) * ddy2_mg(l) * & 817 817 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 818 818 + f2_mg_b(k,l) * p_mg(kp1,j,i) & … … 843 843 j = jj 844 844 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 845 ddx2_mg(l) *&845 rho_air_mg(k,l) * ddx2_mg(l) * & 846 846 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 847 + ddy2_mg(l) *&847 + rho_air_mg(k,l) * ddy2_mg(l) * & 848 848 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 849 849 + f2_mg_b(k,l) * p_mg(kp1,j,i) & … … 852 852 j = jj+2 853 853 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 854 ddx2_mg(l) *&854 rho_air_mg(k,l) * ddx2_mg(l) * & 855 855 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 856 + ddy2_mg(l) *&856 + rho_air_mg(k,l) * ddy2_mg(l) * & 857 857 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 858 858 + f2_mg_b(k,l) * p_mg(kp1,j,i) & … … 869 869 j = jj 870 870 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 871 ddx2_mg(l) *&871 rho_air_mg(k,l) * ddx2_mg(l) * & 872 872 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 873 + ddy2_mg(l) *&873 + rho_air_mg(k,l) * ddy2_mg(l) * & 874 874 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 875 875 + f2_mg_b(k,l) * p_mg(kp1,j,i) & … … 878 878 j = jj+2 879 879 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 880 ddx2_mg(l) *&880 rho_air_mg(k,l) * ddx2_mg(l) * & 881 881 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 882 + ddy2_mg(l) *&882 + rho_air_mg(k,l) * ddy2_mg(l) * & 883 883 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 884 884 + f2_mg_b(k,l) * p_mg(kp1,j,i) & … … 895 895 j = jj 896 896 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 897 ddx2_mg(l) *&897 rho_air_mg(k,l) * ddx2_mg(l) * & 898 898 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 899 + ddy2_mg(l) *&899 + rho_air_mg(k,l) * ddy2_mg(l) * & 900 900 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 901 901 + f2_mg_b(k,l) * p_mg(kp1,j,i) & … … 904 904 j = jj+2 905 905 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 906 ddx2_mg(l) *&906 rho_air_mg(k,l) * ddx2_mg(l) * & 907 907 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 908 + ddy2_mg(l) *&908 + rho_air_mg(k,l) * ddy2_mg(l) * & 909 909 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 910 910 + f2_mg_b(k,l) * p_mg(kp1,j,i) & … … 921 921 j = jj 922 922 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 923 ddx2_mg(l) *&923 rho_air_mg(k,l) * ddx2_mg(l) * & 924 924 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 925 + ddy2_mg(l) *&925 + rho_air_mg(k,l) * ddy2_mg(l) * & 926 926 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 927 927 + f2_mg_b(k,l) * p_mg(kp1,j,i) & … … 930 930 j = jj+2 931 931 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 932 ddx2_mg(l) *&932 rho_air_mg(k,l) * ddx2_mg(l) * & 933 933 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 934 + ddy2_mg(l) *&934 + rho_air_mg(k,l) * ddy2_mg(l) * & 935 935 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 936 936 + f2_mg_b(k,l) * p_mg(kp1,j,i) & -
palm/trunk/SOURCE/poismg_noopt.f90
r2022 r2037 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Anelastic approximation implemented 23 23 ! 24 24 ! Former revisions: … … 268 268 269 269 USE arrays_3d, & 270 ONLY: f1_mg, f2_mg, f3_mg 270 ONLY: f1_mg, f2_mg, f3_mg, rho_air_mg 271 271 272 272 USE control_parameters, & … … 339 339 DO k = nzb+1, nzt_mg(l) 340 340 r(k,j,i) = f_mg(k,j,i) & 341 - ddx2_mg(l) *&341 - rho_air_mg(k,l) * ddx2_mg(l) * & 342 342 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 343 343 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 344 344 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 345 345 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 346 - ddy2_mg(l) *&346 - rho_air_mg(k,l) * ddy2_mg(l) * & 347 347 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 348 348 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & … … 745 745 746 746 USE arrays_3d, & 747 ONLY: f1_mg, f2_mg, f3_mg 747 ONLY: f1_mg, f2_mg, f3_mg, rho_air_mg 748 748 749 749 USE control_parameters, & … … 846 846 847 847 p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( & 848 ddx2_mg(l) *&848 rho_air_mg(k,l) * ddx2_mg(l) * & 849 849 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 850 850 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 851 851 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 852 852 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 853 + ddy2_mg(l) *&853 + rho_air_mg(k,l) * ddy2_mg(l) * & 854 854 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 855 855 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & … … 869 869 DO k = nzb+1, nzt_mg(l), 2 870 870 p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( & 871 ddx2_mg(l) *&871 rho_air_mg(k,l) * ddx2_mg(l) * & 872 872 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 873 873 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 874 874 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 875 875 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 876 + ddy2_mg(l) *&876 + rho_air_mg(k,l) * ddy2_mg(l) * & 877 877 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 878 878 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & … … 892 892 DO k = nzb+2, nzt_mg(l), 2 893 893 p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( & 894 ddx2_mg(l) *&894 rho_air_mg(k,l) * ddx2_mg(l) * & 895 895 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 896 896 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 897 897 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 898 898 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 899 + ddy2_mg(l) *&899 + rho_air_mg(k,l) * ddy2_mg(l) * & 900 900 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 901 901 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & … … 915 915 DO k = nzb+2, nzt_mg(l), 2 916 916 p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( & 917 ddx2_mg(l) *&917 rho_air_mg(k,l) * ddx2_mg(l) * & 918 918 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 919 919 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 920 920 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 921 921 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 922 + ddy2_mg(l) *&922 + rho_air_mg(k,l) * ddy2_mg(l) * & 923 923 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 924 924 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & … … 947 947 j = jj 948 948 p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( & 949 ddx2_mg(l) *&949 rho_air_mg(k,l) * ddx2_mg(l) * & 950 950 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 951 951 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 952 952 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 953 953 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 954 + ddy2_mg(l) *&954 + rho_air_mg(k,l) * ddy2_mg(l) * & 955 955 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 956 956 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & … … 964 964 j = jj+2 965 965 p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( & 966 ddx2_mg(l) *&966 rho_air_mg(k,l) * ddx2_mg(l) * & 967 967 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 968 968 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 969 969 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 970 970 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 971 + ddy2_mg(l) *&971 + rho_air_mg(k,l) * ddy2_mg(l) * & 972 972 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 973 973 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & … … 986 986 j =jj 987 987 p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( & 988 ddx2_mg(l) *&988 rho_air_mg(k,l) * ddx2_mg(l) * & 989 989 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 990 990 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 991 991 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 992 992 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 993 + ddy2_mg(l) *&993 + rho_air_mg(k,l) * ddy2_mg(l) * & 994 994 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 995 995 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & … … 1003 1003 j = jj+2 1004 1004 p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( & 1005 ddx2_mg(l) *&1005 rho_air_mg(k,l) * ddx2_mg(l) * & 1006 1006 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 1007 1007 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 1008 1008 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 1009 1009 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 1010 + ddy2_mg(l) *&1010 + rho_air_mg(k,l) * ddy2_mg(l) * & 1011 1011 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 1012 1012 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & … … 1025 1025 j =jj 1026 1026 p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( & 1027 ddx2_mg(l) *&1027 rho_air_mg(k,l) * ddx2_mg(l) * & 1028 1028 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 1029 1029 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 1030 1030 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 1031 1031 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 1032 + ddy2_mg(l) *&1032 + rho_air_mg(k,l) * ddy2_mg(l) * & 1033 1033 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 1034 1034 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & … … 1042 1042 j = jj+2 1043 1043 p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( & 1044 ddx2_mg(l) *&1044 rho_air_mg(k,l) * ddx2_mg(l) * & 1045 1045 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 1046 1046 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 1047 1047 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 1048 1048 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 1049 + ddy2_mg(l) *&1049 + rho_air_mg(k,l) * ddy2_mg(l) * & 1050 1050 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 1051 1051 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & … … 1064 1064 j =jj 1065 1065 p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( & 1066 ddx2_mg(l) *&1066 rho_air_mg(k,l) * ddx2_mg(l) * & 1067 1067 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 1068 1068 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 1069 1069 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 1070 1070 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 1071 + ddy2_mg(l) *&1071 + rho_air_mg(k,l) * ddy2_mg(l) * & 1072 1072 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 1073 1073 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & … … 1081 1081 j = jj+2 1082 1082 p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( & 1083 ddx2_mg(l) *&1083 rho_air_mg(k,l) * ddx2_mg(l) * & 1084 1084 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 1085 1085 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 1086 1086 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 1087 1087 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 1088 + ddy2_mg(l) *&1088 + rho_air_mg(k,l) * ddy2_mg(l) * & 1089 1089 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 1090 1090 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & -
palm/trunk/SOURCE/pres.f90
r2001 r2037 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Anelastic approximation implemented 23 23 ! 24 24 ! Former revisions: … … 128 128 129 129 USE arrays_3d, & 130 ONLY: d, ddzu, ddzu_pres, ddzw, dzw, p, p_loc, tend, u, v, w 130 ONLY: d, ddzu, ddzu_pres, ddzw, dzw, p, p_loc, rho_air, rho_air_zw, & 131 tend, u, v, w 131 132 132 133 USE control_parameters, & … … 399 400 DO j = nys, nyn 400 401 DO k = nzb_s_inner(j,i)+1, nzt 401 d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx + & 402 ( v(k,j+1,i) - v(k,j,i) ) * ddy + & 403 ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d & 404 * d_weight_pres 402 d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + & 403 ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + & 404 ( w(k,j,i) * rho_air_zw(k) - & 405 w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k) & 406 ) * ddt_3d * d_weight_pres 405 407 ENDDO 406 408 ! … … 427 429 DO j = nys, nyn 428 430 DO k = 1, nzt 429 d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx + & 430 ( v(k,j+1,i) - v(k,j,i) ) * ddy + & 431 ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d & 432 * d_weight_pres * rflags_s_inner(k,j,i) 431 d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + & 432 ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + & 433 ( w(k,j,i) * rho_air_zw(k) - & 434 w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k) & 435 ) * ddt_3d * d_weight_pres * rflags_s_inner(k,j,i) 433 436 ENDDO 434 437 ENDDO … … 808 811 DO j = nys, nyn 809 812 DO k = nzb_s_inner(j,i)+1, nzt 810 d(k,j,i) = ( u(k,j,i+1) - u(k,j,i) ) * ddx + & 811 ( v(k,j+1,i) - v(k,j,i) ) * ddy + & 812 ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 813 d(k,j,i) = ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + & 814 ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + & 815 ( w(k,j,i) * rho_air_zw(k) - & 816 w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k) 813 817 ENDDO 814 818 DO k = nzb+1, nzt … … 823 827 DO j = nys, nyn 824 828 DO k = 1, nzt 825 d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx + & 826 ( v(k,j+1,i) - v(k,j,i) ) * ddy + & 827 ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) & 828 ) * rflags_s_inner(k,j,i) 829 d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + & 830 ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + & 831 ( w(k,j,i) * rho_air_zw(k) - & 832 w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k) & 833 ) * rflags_s_inner(k,j,i) 829 834 ENDDO 830 835 ENDDO -
palm/trunk/SOURCE/sor.f90
r2001 r2037 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Anelastic approximation implemented 23 23 ! 24 24 ! Former revisions: … … 59 59 !------------------------------------------------------------------------------! 60 60 SUBROUTINE sor( d, ddzu, ddzw, p ) 61 61 62 USE arrays_3d, & 63 ONLY: rho_air, rho_air_zw 62 64 63 65 USE grid_variables, & … … 101 103 !-- Compute pre-factors. 102 104 DO k = 1, nz 103 f2(k) = ddzu(k+1) * ddzw(k) 104 f3(k) = ddzu(k) * ddzw(k) 105 f1(k) = 2.0_wp * ( ddx2 + ddy2 ) + f2(k) + f3(k)105 f2(k) = ddzu(k+1) * ddzw(k) * rho_air_zw(k) 106 f3(k) = ddzu(k) * ddzw(k) * rho_air_zw(k-1) 107 f1(k) = 2.0_wp * ( ddx2 + ddy2 ) * rho_air(k) + f2(k) + f3(k) 106 108 ENDDO 107 109 … … 131 133 DO k = nzb+1, nzt 132 134 p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( & 133 ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +&134 ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +&135 f2(k) * p(k+1,j,i) +&136 f3(k) * p(k-1,j,i) -&137 d(k,j,i) -&138 135 rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) + & 136 rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) + & 137 f2(k) * p(k+1,j,i) + & 138 f3(k) * p(k-1,j,i) - & 139 d(k,j,i) - & 140 f1(k) * p(k,j,i) ) 139 141 ENDDO 140 142 ENDDO … … 144 146 DO j = nys1, nyn, 2 145 147 DO k = nzb+1, nzt 146 p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( &147 ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +&148 ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +&149 f2(k) * p(k+1,j,i) +&150 f3(k) * p(k-1,j,i) -&151 d(k,j,i) -&152 148 p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( & 149 rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) + & 150 rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) + & 151 f2(k) * p(k+1,j,i) + & 152 f3(k) * p(k-1,j,i) - & 153 d(k,j,i) - & 154 f1(k) * p(k,j,i) ) 153 155 ENDDO 154 156 ENDDO … … 176 178 DO k = nzb+1, nzt 177 179 p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( & 178 ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +&179 ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +&180 f2(k) * p(k+1,j,i) +&181 f3(k) * p(k-1,j,i) -&182 d(k,j,i) -&183 180 rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) + & 181 rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) + & 182 f2(k) * p(k+1,j,i) + & 183 f3(k) * p(k-1,j,i) - & 184 d(k,j,i) - & 185 f1(k) * p(k,j,i) ) 184 186 ENDDO 185 187 ENDDO … … 190 192 DO k = nzb+1, nzt 191 193 p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( & 192 ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +&193 ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +&194 f2(k) * p(k+1,j,i) +&195 f3(k) * p(k-1,j,i) -&196 d(k,j,i) -&197 194 rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) + & 195 rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) + & 196 f2(k) * p(k+1,j,i) + & 197 f3(k) * p(k-1,j,i) - & 198 d(k,j,i) - & 199 f1(k) * p(k,j,i) ) 198 200 ENDDO 199 201 ENDDO -
palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
r2012 r2037 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Anelastic approximation implemented 23 23 ! 24 24 ! Former revisions: … … 158 158 ONLY: e, kh, nr, nrs, nrsws, ol, pt, q, ql, qr, qrs, qrsws, qs, qsws, & 159 159 s, shf, ss, ssws, ts, u, us, usws, v, vpt, vsws, zu, zw, z0, & 160 z0h, z0q 160 z0h, z0q, drho_air_zw, rho_air_zw 161 161 162 162 USE cloud_parameters, & … … 540 540 rib(j,i) = - g * z_mo * ( ( 1.0_wp + 0.61_wp & 541 541 * q(k+1,j,i) ) * shf(j,i) + 0.61_wp & 542 * pt(k+1,j,i) * qsws(j,i) ) &542 * pt(k+1,j,i) * qsws(j,i) ) * drho_air_zw(k) & 543 543 / ( uv_total(j,i)**3 * vpt(k+1,j,i) * kappa**2& 544 544 + 1.0E-20_wp) 545 545 ELSE 546 rib(j,i) = - g * z_mo * shf(j,i) 546 rib(j,i) = - g * z_mo * shf(j,i) * drho_air_zw(k) & 547 547 / ( uv_total(j,i)**3 * pt(k+1,j,i) * kappa**2 & 548 548 + 1.0E-20_wp ) … … 828 828 !-- For a given heat flux in the surface layer: 829 829 !$OMP PARALLEL DO 830 !$acc kernels loop private( j )830 !$acc kernels loop private( j, k ) 831 831 DO i = nxlg, nxrg 832 832 DO j = nysg, nyng 833 ts(j,i) = -shf(j,i) / ( us(j,i) + 1E-30_wp ) 833 k = nzb_s_inner(j,i) 834 ts(j,i) = -shf(j,i) * drho_air_zw(k) / ( us(j,i) + 1E-30_wp ) 834 835 ! 835 836 !-- ts must be limited, because otherwise overflow may occur in case … … 888 889 DO i = nxlg, nxrg 889 890 DO j = nysg, nyng 890 qs(j,i) = -qsws(j,i) / ( us(j,i) + 1E-30_wp ) 891 k = nzb_s_inner(j,i) 892 qs(j,i) = -qsws(j,i) * drho_air_zw(k) / ( us(j,i) + 1E-30_wp ) 891 893 ENDDO 892 894 ENDDO … … 1025 1027 + psi_m( z0(j,i) / ol_mid ) ) 1026 1028 1027 usws(j,i) = -usws(j,i) * 0.5_wp * ( us(j,i-1) + us(j,i) ) 1029 usws(j,i) = -usws(j,i) * 0.5_wp * ( us(j,i-1) + us(j,i) ) & 1030 * rho_air_zw(k) 1028 1031 ENDDO 1029 1032 ENDDO … … 1052 1055 + psi_m( z0(j,i) / ol_mid ) ) 1053 1056 1054 vsws(j,i) = -vsws(j,i) * 0.5_wp * ( us(j,i-1) + us(j,i) ) 1057 vsws(j,i) = -vsws(j,i) * 0.5_wp * ( us(j,i-1) + us(j,i) ) & 1058 * rho_air_zw(k) 1055 1059 1056 1060 ENDDO … … 1074 1078 !$acc loop independent 1075 1079 DO j = nysg, nyng 1076 shf(j,i) = -ts(j,i) * us(j,i) 1080 k = nzb_s_inner(j,i) 1081 shf(j,i) = -ts(j,i) * us(j,i) * rho_air_zw(k) 1077 1082 ENDDO 1078 1083 ENDDO … … 1090 1095 !$acc loop independent 1091 1096 DO j = nysg, nyng 1092 qsws(j,i) = -qs(j,i) * us(j,i) 1097 k = nzb_s_inner(j,i) 1098 qsws(j,i) = -qs(j,i) * us(j,i) * rho_air_zw(k) 1093 1099 ENDDO 1094 1100 ENDDO -
palm/trunk/SOURCE/tridia_solver_mod.f90
r2001 r2037 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Anelastic approximation implemented 23 23 ! 24 24 ! Former revisions: … … 129 129 130 130 USE arrays_3d, & 131 ONLY: ddzu_pres, ddzw 131 ONLY: ddzu_pres, ddzw, rho_air_zw 132 132 133 133 USE kinds … … 140 140 141 141 DO k = 0, nz-1 142 ddzuw(k,1) = ddzu_pres(k+1) * ddzw(k+1) 143 ddzuw(k,2) = ddzu_pres(k+2) * ddzw(k+1) 142 ddzuw(k,1) = ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k) 143 ddzuw(k,2) = ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1) 144 144 ddzuw(k,3) = -1.0_wp * & 145 ( ddzu_pres(k+2) * ddzw(k+1) + ddzu_pres(k+1) * ddzw(k+1) ) 145 ( ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1) + & 146 ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k) ) 146 147 ENDDO 147 148 ! … … 168 169 169 170 USE arrays_3d, & 170 ONLY: tric 171 ONLY: tric, rho_air 171 172 172 173 USE constants, & … … 231 232 DO j = nys_z, nyn_z 232 233 DO i = nxl_z, nxr_z 233 tric(i,j,k) = ddzuw(k,3) - ll(i,j) 234 tric(i,j,k) = ddzuw(k,3) - ll(i,j) * rho_air(k+1) 234 235 ENDDO 235 236 ENDDO … … 491 492 492 493 USE arrays_3d, & 493 ONLY: ddzu_pres, ddzw 494 ONLY: ddzu_pres, ddzw, rho_air, rho_air_zw 494 495 495 496 USE control_parameters, & … … 526 527 DO k = 0, nz-1 527 528 DO i = 0,nx 528 tri_for_1d(2,i,k) = ddzu_pres(k+1) * ddzw(k+1) 529 tri_for_1d(3,i,k) = ddzu_pres(k+2) * ddzw(k+1) 529 tri_for_1d(2,i,k) = ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k) 530 tri_for_1d(3,i,k) = ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1) 530 531 ENDDO 531 532 ENDDO … … 591 592 DO k = 0, nz-1 592 593 DO i = 0, nx 593 a = -1.0_wp * ddzu_pres(k+2) * ddzw(k+1) 594 c = -1.0_wp * ddzu_pres(k+1) * ddzw(k+1) 595 tri_for_1d(1,i,k) = a + c - l(i) 594 a = -1.0_wp * ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1) 595 c = -1.0_wp * ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k) 596 tri_for_1d(1,i,k) = a + c - l(i) * rho_air(k+1) 596 597 ENDDO 597 598 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.
