Changeset 2200 for palm/trunk/SOURCE/advec_ws.f90
- Timestamp:
- Apr 11, 2017 11:37:51 AM (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/advec_ws.f90
r2119 r2200 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! monotonic_adjustment removed 23 23 ! 24 24 ! Former revisions: … … 245 245 USE control_parameters, & 246 246 ONLY: cloud_physics, humidity, loop_optimization, & 247 monotonic_adjustment, passive_scalar, microphysics_seifert,&247 passive_scalar, microphysics_seifert, & 248 248 ocean, ws_scheme_mom, ws_scheme_sca 249 249 … … 258 258 ONLY: sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wsnrs_ws_l,& 259 259 sums_wspts_ws_l, sums_wsqrs_ws_l, sums_wsqs_ws_l, & 260 sums_wsss_ws_l, sums_wssas_ws_l, sums_wsss_ws_l, 260 sums_wsss_ws_l, sums_wssas_ws_l, sums_wsss_ws_l, & 261 261 sums_wsus_ws_l, sums_wsvs_ws_l 262 262 … … 418 418 419 419 420 IF ( scalar_advec == 'ws-scheme' .OR. & 421 scalar_advec == 'ws-scheme-mono' ) THEN 420 IF ( scalar_advec == 'ws-scheme' ) THEN 422 421 ! 423 422 !-- Set flags to steer the degradation of the advection scheme in advec_ws … … 803 802 ! 804 803 !-- Exchange 3D integer wall_flags. 805 IF ( momentum_advec == 'ws-scheme' .OR. scalar_advec == 'ws-scheme' &806 .OR. scalar_advec == 'ws-scheme-mono') THEN804 IF ( momentum_advec == 'ws-scheme' .OR. scalar_advec == 'ws-scheme' & 805 ) THEN 807 806 ! 808 807 !-- Exchange ghost points for advection flags … … 901 900 902 901 USE control_parameters, & 903 ONLY: intermediate_timestep_count, monotonic_adjustment, u_gtrans, & 904 v_gtrans 902 ONLY: intermediate_timestep_count, u_gtrans, v_gtrans 905 903 906 904 USE grid_variables, & … … 946 944 REAL(wp) :: div !< 947 945 REAL(wp) :: flux_d !< 948 REAL(wp) :: fd_1 !<949 REAL(wp) :: fl_1 !<950 REAL(wp) :: fn_1 !<951 REAL(wp) :: fr_1 !<952 REAL(wp) :: fs_1 !<953 REAL(wp) :: ft_1 !<954 REAL(wp) :: phi_d !<955 REAL(wp) :: phi_l !<956 REAL(wp) :: phi_n !<957 REAL(wp) :: phi_r !<958 REAL(wp) :: phi_s !<959 REAL(wp) :: phi_t !<960 REAL(wp) :: rd !<961 REAL(wp) :: rl !<962 REAL(wp) :: rn !<963 REAL(wp) :: rr !<964 REAL(wp) :: rs !<965 REAL(wp) :: rt !<966 946 REAL(wp) :: u_comp !< 967 947 REAL(wp) :: v_comp !< … … 1232 1212 ( sk(k_ppp,j,i) - sk(k_mm,j,i) ) & 1233 1213 ) 1234 !1235 !-- Apply monotonic adjustment.1236 IF ( monotonic_adjustment ) THEN1237 !1238 !-- At first, calculate first order fluxes.1239 u_comp = u(k,j,i) - u_gtrans1240 fl_1 = ( u_comp * ( sk(k,j,i) + sk(k,j,i-1) ) &1241 -ABS( u_comp ) * ( sk(k,j,i) - sk(k,j,i-1) ) &1242 ) * adv_sca_11243 1244 u_comp = u(k,j,i+1) - u_gtrans1245 fr_1 = ( u_comp * ( sk(k,j,i+1) + sk(k,j,i) ) &1246 -ABS( u_comp ) * ( sk(k,j,i+1) - sk(k,j,i) ) &1247 ) * adv_sca_11248 1249 v_comp = v(k,j,i) - v_gtrans1250 fs_1 = ( v_comp * ( sk(k,j,i) + sk(k,j-1,i) ) &1251 -ABS( v_comp ) * ( sk(k,j,i) - sk(k,j-1,i) ) &1252 ) * adv_sca_11253 1254 v_comp = v(k,j+1,i) - v_gtrans1255 fn_1 = ( v_comp * ( sk(k,j+1,i) + sk(k,j,i) ) &1256 -ABS( v_comp ) * ( sk(k,j+1,i) - sk(k,j,i) ) &1257 ) * adv_sca_11258 1259 fd_1 = ( w(k-1,j,i) * ( sk(k,j,i) + sk(k-1,j,i) ) &1260 -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) ) &1261 ) * adv_sca_1 * rho_air_zw(k)1262 1263 ft_1 = ( w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) &1264 -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) ) &1265 ) * adv_sca_1 * rho_air_zw(k)1266 !1267 !-- Calculate ratio of upwind gradients. Note, Min/Max is just to1268 !-- avoid if statements.1269 rl = ( MAX( 0.0_wp, u(k,j,i) - u_gtrans ) * &1270 ABS( ( sk(k,j,i-1) - sk(k,j,i-2) ) / &1271 ( sk(k,j,i) - sk(k,j,i-1) + 1E-20_wp ) &1272 ) + &1273 MIN( 0.0_wp, u(k,j,i) - u_gtrans ) * &1274 ABS( ( sk(k,j,i) - sk(k,j,i+1) ) / &1275 ( sk(k,j,i-1) - sk(k,j,i) + 1E-20_wp ) &1276 ) &1277 ) / ABS( u(k,j,i) - u_gtrans + 1E-20_wp )1278 1279 rr = ( MAX( 0.0_wp, u(k,j,i+1) - u_gtrans ) * &1280 ABS( ( sk(k,j,i) - sk(k,j,i-1) ) / &1281 ( sk(k,j,i+1) - sk(k,j,i) + 1E-20_wp ) &1282 ) + &1283 MIN( 0.0_wp, u(k,j,i+1) - u_gtrans ) * &1284 ABS( ( sk(k,j,i+1) - sk(k,j,i+2) ) / &1285 ( sk(k,j,i) - sk(k,j,i+1) + 1E-20_wp ) &1286 ) &1287 ) / ABS( u(k,j,i+1) - u_gtrans + 1E-20_wp )1288 1289 rs = ( MAX( 0.0_wp, v(k,j,i) - v_gtrans ) * &1290 ABS( ( sk(k,j-1,i) - sk(k,j-2,i) ) / &1291 ( sk(k,j,i) - sk(k,j-1,i) + 1E-20_wp ) &1292 ) + &1293 MIN( 0.0_wp, v(k,j,i) - v_gtrans ) * &1294 ABS( ( sk(k,j,i) - sk(k,j+1,i) ) / &1295 ( sk(k,j-1,i) - sk(k,j,i) + 1E-20_wp ) &1296 ) &1297 ) / ABS( v(k,j,i) - v_gtrans + 1E-20_wp )1298 1299 rn = ( MAX( 0.0_wp, v(k,j+1,i) - v_gtrans ) * &1300 ABS( ( sk(k,j,i) - sk(k,j-1,i) ) / &1301 ( sk(k,j+1,i) - sk(k,j,i) + 1E-20_wp ) &1302 ) + &1303 MIN( 0.0_wp, v(k,j+1,i) - v_gtrans ) * &1304 ABS( ( sk(k,j+1,i) - sk(k,j+2,i) ) / &1305 ( sk(k,j,i) - sk(k,j+1,i) + 1E-20_wp ) &1306 ) &1307 ) / ABS( v(k,j+1,i) - v_gtrans + 1E-20_wp )1308 !1309 !-- Reuse k_mm and compute k_mmm for the vertical gradient ratios.1310 !-- Note, for vertical advection below the third grid point above1311 !-- surface ( or below the model top) rd and rt are set to 0, i.e.1312 !-- use of first order scheme is enforced.1313 k_mmm = k - 3 * ibit81314 1315 rd = ( MAX( 0.0_wp, w(k-1,j,i) ) * &1316 ABS( ( sk(k_mm,j,i) - sk(k_mmm,j,i) ) / &1317 ( sk(k-1,j,i) - sk(k_mm,j,i) + 1E-20_wp ) &1318 ) + &1319 MIN( 0.0_wp, w(k-1,j,i) ) * &1320 ABS( ( sk(k-1,j,i) - sk(k,j,i) ) / &1321 ( sk(k_mm,j,i) - sk(k-1,j,i) + 1E-20_wp ) &1322 ) &1323 ) * ibit8 / ABS( w(k-1,j,i) + 1E-20_wp )1324 1325 rt = ( MAX( 0.0_wp, w(k,j,i) ) * &1326 ABS( ( sk(k,j,i) - sk(k-1,j,i) ) / &1327 ( sk(k+1,j,i) - sk(k,j,i) + 1E-20_wp ) &1328 ) + &1329 MIN( 0.0_wp, w(k,j,i) ) * &1330 ABS( ( sk(k+1,j,i) - sk(k_pp,j,i) ) / &1331 ( sk(k,j,i) - sk(k+1,j,i) + 1E-20_wp ) &1332 ) &1333 ) * ibit8 / ABS( w(k,j,i) + 1E-20_wp )1334 !1335 !-- Calculate empirical limiter function (van Albada2 limiter).1336 phi_l = MIN( 1.0_wp, ( 2.0_wp * ABS( rl ) ) / &1337 ( rl**2 + 1.0_wp ) )1338 phi_r = MIN( 1.0_wp, ( 2.0_wp * ABS( rr ) ) / &1339 ( rr**2 + 1.0_wp ) )1340 phi_s = MIN( 1.0_wp, ( 2.0_wp * ABS( rs ) ) / &1341 ( rs**2 + 1.0_wp ) )1342 phi_n = MIN( 1.0_wp, ( 2.0_wp * ABS( rn ) ) / &1343 ( rn**2 + 1.0_wp ) )1344 phi_d = MIN( 1.0_wp, ( 2.0_wp * ABS( rd ) ) / &1345 ( rd**2 + 1.0_wp ) )1346 phi_t = MIN( 1.0_wp, ( 2.0_wp * ABS( rt ) ) / &1347 ( rt**2 + 1.0_wp ) )1348 !1349 !-- Calculate the resulting monotone flux.1350 swap_flux_x_local(k,j,tn) = fl_1 - phi_l * &1351 ( fl_1 - swap_flux_x_local(k,j,tn) )1352 flux_r(k) = fr_1 - phi_r * &1353 ( fr_1 - flux_r(k) )1354 swap_flux_y_local(k,tn) = fs_1 - phi_s * &1355 ( fs_1 - swap_flux_y_local(k,tn) )1356 flux_n(k) = fn_1 - phi_n * &1357 ( fn_1 - flux_n(k) )1358 flux_d = fd_1 - phi_d * &1359 ( fd_1 - flux_d )1360 flux_t(k) = ft_1 - phi_t * &1361 ( ft_1 - flux_t(k) )1362 !1363 !-- Moreover, modify dissipation flux according to the limiter.1364 swap_diss_x_local(k,j,tn) = swap_diss_x_local(k,j,tn) * phi_l1365 diss_r(k) = diss_r(k) * phi_r1366 swap_diss_y_local(k,tn) = swap_diss_y_local(k,tn) * phi_s1367 diss_n(k) = diss_n(k) * phi_n1368 diss_d = diss_d * phi_d1369 diss_t(k) = diss_t(k) * phi_t1370 1371 ENDIF1372 1214 ! 1373 1215 !-- Calculate the divergence of the velocity field. A respective … … 1478 1320 ( sk(k_ppp,j,i) - sk(k_mm,j,i) ) & 1479 1321 ) 1480 1481 1482 !1483 !-- Apply monotonic adjustment.1484 IF ( monotonic_adjustment ) THEN1485 !1486 !-- At first, calculate first order fluxes.1487 u_comp = u(k,j,i) - u_gtrans1488 fl_1 = ( u_comp * ( sk(k,j,i) + sk(k,j,i-1) ) &1489 -ABS( u_comp ) * ( sk(k,j,i) - sk(k,j,i-1) ) &1490 ) * adv_sca_11491 1492 u_comp = u(k,j,i+1) - u_gtrans1493 fr_1 = ( u_comp * ( sk(k,j,i+1) + sk(k,j,i) ) &1494 -ABS( u_comp ) * ( sk(k,j,i+1) - sk(k,j,i) ) &1495 ) * adv_sca_11496 1497 v_comp = v(k,j,i) - v_gtrans1498 fs_1 = ( v_comp * ( sk(k,j,i) + sk(k,j-1,i) ) &1499 -ABS( v_comp ) * ( sk(k,j,i) - sk(k,j-1,i) ) &1500 ) * adv_sca_11501 1502 v_comp = v(k,j+1,i) - v_gtrans1503 fn_1 = ( v_comp * ( sk(k,j+1,i) + sk(k,j,i) ) &1504 -ABS( v_comp ) * ( sk(k,j+1,i) - sk(k,j,i) ) &1505 ) * adv_sca_11506 1507 fd_1 = ( w(k-1,j,i) * ( sk(k,j,i) + sk(k-1,j,i) ) &1508 -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) ) &1509 ) * adv_sca_1 * rho_air_zw(k)1510 1511 ft_1 = ( w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) &1512 -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) ) &1513 ) * adv_sca_1 * rho_air_zw(k)1514 !1515 !-- Calculate ratio of upwind gradients. Note, Min/Max is just to1516 !-- avoid if statements.1517 rl = ( MAX( 0.0_wp, u(k,j,i) - u_gtrans ) * &1518 ABS( ( sk(k,j,i-1) - sk(k,j,i-2) ) / &1519 ( sk(k,j,i) - sk(k,j,i-1) + 1E-20_wp ) &1520 ) + &1521 MIN( 0.0_wp, u(k,j,i) - u_gtrans ) * &1522 ABS( ( sk(k,j,i) - sk(k,j,i+1) ) / &1523 ( sk(k,j,i-1) - sk(k,j,i) + 1E-20_wp ) &1524 ) &1525 ) / ABS( u(k,j,i) - u_gtrans + 1E-20_wp )1526 1527 rr = ( MAX( 0.0_wp, u(k,j,i+1) - u_gtrans ) * &1528 ABS( ( sk(k,j,i) - sk(k,j,i-1) ) / &1529 ( sk(k,j,i+1) - sk(k,j,i) + 1E-20_wp ) &1530 ) + &1531 MIN( 0.0_wp, u(k,j,i+1) - u_gtrans ) * &1532 ABS( ( sk(k,j,i+1) - sk(k,j,i+2) ) / &1533 ( sk(k,j,i) - sk(k,j,i+1) + 1E-20_wp ) &1534 ) &1535 ) / ABS( u(k,j,i+1) - u_gtrans + 1E-20_wp )1536 1537 rs = ( MAX( 0.0_wp, v(k,j,i) - v_gtrans ) * &1538 ABS( ( sk(k,j-1,i) - sk(k,j-2,i) ) / &1539 ( sk(k,j,i) - sk(k,j-1,i) + 1E-20_wp ) &1540 ) + &1541 MIN( 0.0_wp, v(k,j,i) - v_gtrans ) * &1542 ABS( ( sk(k,j,i) - sk(k,j+1,i) ) / &1543 ( sk(k,j-1,i) - sk(k,j,i) + 1E-20_wp ) &1544 ) &1545 ) / ABS( v(k,j,i) - v_gtrans + 1E-20_wp )1546 1547 rn = ( MAX( 0.0_wp, v(k,j+1,i) - v_gtrans ) * &1548 ABS( ( sk(k,j,i) - sk(k,j-1,i) ) / &1549 ( sk(k,j+1,i) - sk(k,j,i) + 1E-20_wp ) &1550 ) + &1551 MIN( 0.0_wp, v(k,j+1,i) - v_gtrans ) * &1552 ABS( ( sk(k,j+1,i) - sk(k,j+2,i) ) / &1553 ( sk(k,j,i) - sk(k,j+1,i) + 1E-20_wp ) &1554 ) &1555 ) / ABS( v(k,j+1,i) - v_gtrans + 1E-20_wp )1556 !1557 !-- Reuse k_mm and compute k_mmm for the vertical gradient ratios.1558 !-- Note, for vertical advection below the third grid point above1559 !-- surface ( or below the model top) rd and rt are set to 0, i.e.1560 !-- use of first order scheme is enforced.1561 k_mmm = k - 3 * ibit81562 1563 rd = ( MAX( 0.0_wp, w(k-1,j,i) ) * &1564 ABS( ( sk(k_mm,j,i) - sk(k_mmm,j,i) ) / &1565 ( sk(k-1,j,i) - sk(k_mm,j,i) + 1E-20_wp ) &1566 ) + &1567 MIN( 0.0_wp, w(k-1,j,i) ) * &1568 ABS( ( sk(k-1,j,i) - sk(k,j,i) ) / &1569 ( sk(k_mm,j,i) - sk(k-1,j,i) + 1E-20_wp ) &1570 ) &1571 ) * ibit8 / ABS( w(k-1,j,i) + 1E-20_wp )1572 1573 rt = ( MAX( 0.0_wp, w(k,j,i) ) * &1574 ABS( ( sk(k,j,i) - sk(k-1,j,i) ) / &1575 ( sk(k+1,j,i) - sk(k,j,i) + 1E-20_wp ) &1576 ) + &1577 MIN( 0.0_wp, w(k,j,i) ) * &1578 ABS( ( sk(k+1,j,i) - sk(k_pp,j,i) ) / &1579 ( sk(k,j,i) - sk(k+1,j,i) + 1E-20_wp ) &1580 ) &1581 ) * ibit8 / ABS( w(k,j,i) + 1E-20_wp )1582 !1583 !-- Calculate empirical limiter function (van Albada2 limiter).1584 phi_l = MIN( 1.0_wp, ( 2.0_wp * ABS( rl ) ) / &1585 ( rl**2 + 1.0_wp ) )1586 phi_r = MIN( 1.0_wp, ( 2.0_wp * ABS( rr ) ) / &1587 ( rr**2 + 1.0_wp ) )1588 phi_s = MIN( 1.0_wp, ( 2.0_wp * ABS( rs ) ) / &1589 ( rs**2 + 1.0_wp ) )1590 phi_n = MIN( 1.0_wp, ( 2.0_wp * ABS( rn ) ) / &1591 ( rn**2 + 1.0_wp ) )1592 phi_d = MIN( 1.0_wp, ( 2.0_wp * ABS( rd ) ) / &1593 ( rd**2 + 1.0_wp ) )1594 phi_t = MIN( 1.0_wp, ( 2.0_wp * ABS( rt ) ) / &1595 ( rt**2 + 1.0_wp ) )1596 !1597 !-- Calculate the resulting monotone flux.1598 swap_flux_x_local(k,j,tn) = fl_1 - phi_l * &1599 ( fl_1 - swap_flux_x_local(k,j,tn) )1600 flux_r(k) = fr_1 - phi_r * &1601 ( fr_1 - flux_r(k) )1602 swap_flux_y_local(k,tn) = fs_1 - phi_s * &1603 ( fs_1 - swap_flux_y_local(k,tn) )1604 flux_n(k) = fn_1 - phi_n * &1605 ( fn_1 - flux_n(k) )1606 flux_d = fd_1 - phi_d * &1607 ( fd_1 - flux_d )1608 flux_t(k) = ft_1 - phi_t * &1609 ( ft_1 - flux_t(k) )1610 !1611 !-- Moreover, modify dissipation flux according to the limiter.1612 swap_diss_x_local(k,j,tn) = swap_diss_x_local(k,j,tn) * phi_l1613 diss_r(k) = diss_r(k) * phi_r1614 swap_diss_y_local(k,tn) = swap_diss_y_local(k,tn) * phi_s1615 diss_n(k) = diss_n(k) * phi_n1616 diss_d = diss_d * phi_d1617 diss_t(k) = diss_t(k) * phi_t1618 1619 ENDIF1620 1322 ! 1621 1323 !-- Calculate the divergence of the velocity field. A respective … … 3207 2909 3208 2910 USE control_parameters, & 3209 ONLY: intermediate_timestep_count, monotonic_adjustment, u_gtrans,& 3210 v_gtrans 2911 ONLY: intermediate_timestep_count, u_gtrans, v_gtrans 3211 2912 3212 2913 USE grid_variables, & … … 3255 2956 REAL(wp) :: div !< 3256 2957 REAL(wp) :: flux_d !< 3257 REAL(wp) :: fd_1 !<3258 REAL(wp) :: fl_1 !<3259 REAL(wp) :: fn_1 !<3260 REAL(wp) :: fr_1 !<3261 REAL(wp) :: fs_1 !<3262 REAL(wp) :: ft_1 !<3263 REAL(wp) :: phi_d !<3264 REAL(wp) :: phi_l !<3265 REAL(wp) :: phi_n !<3266 REAL(wp) :: phi_r !<3267 REAL(wp) :: phi_s !<3268 REAL(wp) :: phi_t !<3269 REAL(wp) :: rd !<3270 REAL(wp) :: rl !<3271 REAL(wp) :: rn !<3272 REAL(wp) :: rr !<3273 REAL(wp) :: rs !<3274 REAL(wp) :: rt !<3275 2958 REAL(wp) :: u_comp !< 3276 2959 REAL(wp) :: v_comp !< … … 3530 3213 ( sk(k_ppp,j,i) - sk(k_mm,j,i) ) & 3531 3214 ) 3532 !3533 !-- Apply monotonic adjustment.3534 IF ( monotonic_adjustment ) THEN3535 !3536 !-- At first, calculate first order fluxes.3537 u_comp = u(k,j,i) - u_gtrans3538 fl_1 = ( u_comp * ( sk(k,j,i) + sk(k,j,i-1) ) &3539 -ABS( u_comp ) * ( sk(k,j,i) - sk(k,j,i-1) ) &3540 ) * adv_sca_13541 3542 u_comp = u(k,j,i+1) - u_gtrans3543 fr_1 = ( u_comp * ( sk(k,j,i+1) + sk(k,j,i) ) &3544 -ABS( u_comp ) * ( sk(k,j,i+1) - sk(k,j,i) ) &3545 ) * adv_sca_13546 3547 v_comp = v(k,j,i) - v_gtrans3548 fs_1 = ( v_comp * ( sk(k,j,i) + sk(k,j-1,i) ) &3549 -ABS( v_comp ) * ( sk(k,j,i) - sk(k,j-1,i) ) &3550 ) * adv_sca_13551 3552 v_comp = v(k,j+1,i) - v_gtrans3553 fn_1 = ( v_comp * ( sk(k,j+1,i) + sk(k,j,i) ) &3554 -ABS( v_comp ) * ( sk(k,j+1,i) - sk(k,j,i) ) &3555 ) * adv_sca_13556 3557 fd_1 = ( w(k-1,j,i) * ( sk(k,j,i) + sk(k-1,j,i) ) &3558 -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) ) &3559 ) * adv_sca_1 * rho_air_zw(k)3560 3561 ft_1 = ( w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) &3562 -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) ) &3563 ) * adv_sca_1 * rho_air_zw(k)3564 !3565 !-- Calculate ratio of upwind gradients. Note, Min/Max is just3566 !-- to avoid if statements.3567 rl = ( MAX( 0.0_wp, u(k,j,i) - u_gtrans ) * &3568 ABS( ( sk(k,j,i-1) - sk(k,j,i-2) ) /&3569 ( sk(k,j,i) - sk(k,j,i-1) + 1E-20_wp ) &3570 ) + &3571 MIN( 0.0_wp, u(k,j,i) - u_gtrans ) * &3572 ABS( ( sk(k,j,i) - sk(k,j,i+1) ) /&3573 ( sk(k,j,i-1) - sk(k,j,i) + 1E-20_wp ) &3574 ) &3575 ) / ABS( u(k,j,i) - u_gtrans + 1E-20_wp )3576 3577 rr = ( MAX( 0.0_wp, u(k,j,i+1) - u_gtrans ) * &3578 ABS( ( sk(k,j,i) - sk(k,j,i-1) ) /&3579 ( sk(k,j,i+1) - sk(k,j,i) + 1E-20_wp ) &3580 ) + &3581 MIN( 0.0_wp, u(k,j,i+1) - u_gtrans ) * &3582 ABS( ( sk(k,j,i+1) - sk(k,j,i+2) ) /&3583 ( sk(k,j,i) - sk(k,j,i+1) + 1E-20_wp ) &3584 ) &3585 ) / ABS( u(k,j,i+1) - u_gtrans + 1E-20_wp )3586 3587 rs = ( MAX( 0.0_wp, v(k,j,i) - v_gtrans ) * &3588 ABS( ( sk(k,j-1,i) - sk(k,j-2,i) ) /&3589 ( sk(k,j,i) - sk(k,j-1,i) + 1E-20_wp ) &3590 ) + &3591 MIN( 0.0_wp, v(k,j,i) - v_gtrans ) * &3592 ABS( ( sk(k,j,i) - sk(k,j+1,i) ) /&3593 ( sk(k,j-1,i) - sk(k,j,i) + 1E-20_wp ) &3594 ) &3595 ) / ABS( v(k,j,i) - v_gtrans + 1E-20_wp )3596 3597 rn = ( MAX( 0.0_wp, v(k,j+1,i) - v_gtrans ) * &3598 ABS( ( sk(k,j,i) - sk(k,j-1,i) ) /&3599 ( sk(k,j+1,i) - sk(k,j,i) + 1E-20_wp ) &3600 ) + &3601 MIN( 0.0_wp, v(k,j+1,i) - v_gtrans ) * &3602 ABS( ( sk(k,j+1,i) - sk(k,j+2,i) ) /&3603 ( sk(k,j,i) - sk(k,j+1,i) + 1E-20_wp ) &3604 ) &3605 ) / ABS( v(k,j+1,i) - v_gtrans + 1E-20_wp )3606 !3607 !-- Reuse k_mm and compute k_mmm for the vertical gradient ratios.3608 !-- Note, for vertical advection below the third grid point above3609 !-- surface ( or below the model top) rd and rt are set to 0, i.e.3610 !-- use of first order scheme is enforced.3611 k_mmm = k - 3 * ibit83612 3613 rd = ( MAX( 0.0_wp, w(k-1,j,i) ) * &3614 ABS( ( sk(k_mm,j,i) - sk(k_mmm,j,i) ) / &3615 ( sk(k-1,j,i) - sk(k_mm,j,i) + 1E-20_wp ) &3616 ) + &3617 MIN( 0.0_wp, w(k-1,j,i) ) * &3618 ABS( ( sk(k-1,j,i) - sk(k,j,i) ) / &3619 ( sk(k_mm,j,i) - sk(k-1,j,i) + 1E-20_wp ) &3620 ) &3621 ) * ibit8 / ABS( w(k-1,j,i) + 1E-20_wp )3622 3623 rt = ( MAX( 0.0_wp, w(k,j,i) ) * &3624 ABS( ( sk(k,j,i) - sk(k-1,j,i) ) / &3625 ( sk(k+1,j,i) - sk(k,j,i) + 1E-20_wp ) &3626 ) + &3627 MIN( 0.0_wp, w(k,j,i) ) * &3628 ABS( ( sk(k+1,j,i) - sk(k_pp,j,i) ) / &3629 ( sk(k,j,i) - sk(k+1,j,i) + 1E-20_wp ) &3630 ) &3631 ) * ibit8 / ABS( w(k,j,i) + 1E-20_wp )3632 !3633 !-- Calculate empirical limiter function (van Albada2 limiter).3634 phi_l = MIN( 1.0_wp, ( 2.0_wp * ABS( rl ) ) / &3635 ( rl**2 + 1.0_wp ) )3636 phi_r = MIN( 1.0_wp, ( 2.0_wp * ABS( rr ) ) / &3637 ( rr**2 + 1.0_wp ) )3638 phi_s = MIN( 1.0_wp, ( 2.0_wp * ABS( rs ) ) / &3639 ( rs**2 + 1.0_wp ) )3640 phi_n = MIN( 1.0_wp, ( 2.0_wp * ABS( rn ) ) / &3641 ( rn**2 + 1.0_wp ) )3642 phi_d = MIN( 1.0_wp, ( 2.0_wp * ABS( rd ) ) / &3643 ( rd**2 + 1.0_wp ) )3644 phi_t = MIN( 1.0_wp, ( 2.0_wp * ABS( rt ) ) / &3645 ( rt**2 + 1.0_wp ) )3646 !3647 !-- Calculate the resulting monotone flux.3648 swap_flux_x_local(k,j) = fl_1 - phi_l * &3649 ( fl_1 - swap_flux_x_local(k,j) )3650 flux_r(k) = fr_1 - phi_r * &3651 ( fr_1 - flux_r(k) )3652 swap_flux_y_local(k) = fs_1 - phi_s * &3653 ( fs_1 - swap_flux_y_local(k) )3654 flux_n(k) = fn_1 - phi_n * &3655 ( fn_1 - flux_n(k) )3656 flux_d = fd_1 - phi_d * &3657 ( fd_1 - flux_d )3658 flux_t(k) = ft_1 - phi_t * &3659 ( ft_1 - flux_t(k) )3660 !3661 !-- Moreover, modify dissipation flux according to the limiter.3662 swap_diss_x_local(k,j) = swap_diss_x_local(k,j) * phi_l3663 diss_r(k) = diss_r(k) * phi_r3664 swap_diss_y_local(k) = swap_diss_y_local(k) * phi_s3665 diss_n(k) = diss_n(k) * phi_n3666 diss_d = diss_d * phi_d3667 diss_t(k) = diss_t(k) * phi_t3668 3669 ENDIF3670 3215 ! 3671 3216 !-- Calculate the divergence of the velocity field. A respective … … 3774 3319 ( sk(k_ppp,j,i) - sk(k_mm,j,i) ) & 3775 3320 ) 3776 !3777 !-- Apply monotonic adjustment.3778 IF ( monotonic_adjustment ) THEN3779 !3780 !-- At first, calculate first order fluxes.3781 u_comp = u(k,j,i) - u_gtrans3782 fl_1 = ( u_comp * ( sk(k,j,i) + sk(k,j,i-1) ) &3783 -ABS( u_comp ) * ( sk(k,j,i) - sk(k,j,i-1) ) &3784 ) * adv_sca_13785 3786 u_comp = u(k,j,i+1) - u_gtrans3787 fr_1 = ( u_comp * ( sk(k,j,i+1) + sk(k,j,i) ) &3788 -ABS( u_comp ) * ( sk(k,j,i+1) - sk(k,j,i) ) &3789 ) * adv_sca_13790 3791 v_comp = v(k,j,i) - v_gtrans3792 fs_1 = ( v_comp * ( sk(k,j,i) + sk(k,j-1,i) ) &3793 -ABS( v_comp ) * ( sk(k,j,i) - sk(k,j-1,i) ) &3794 ) * adv_sca_13795 3796 v_comp = v(k,j+1,i) - v_gtrans3797 fn_1 = ( v_comp * ( sk(k,j+1,i) + sk(k,j,i) ) &3798 -ABS( v_comp ) * ( sk(k,j+1,i) - sk(k,j,i) ) &3799 ) * adv_sca_13800 3801 fd_1 = ( w(k-1,j,i) * ( sk(k,j,i) + sk(k-1,j,i) ) &3802 -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) ) &3803 ) * adv_sca_1 * rho_air_zw(k)3804 3805 ft_1 = ( w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) &3806 -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) ) &3807 ) * adv_sca_1 * rho_air_zw(k)3808 !3809 !-- Calculate ratio of upwind gradients. Note, Min/Max is just3810 !-- to avoid if statements.3811 rl = ( MAX( 0.0_wp, u(k,j,i) - u_gtrans ) * &3812 ABS( ( sk(k,j,i-1) - sk(k,j,i-2) ) /&3813 ( sk(k,j,i) - sk(k,j,i-1) + 1E-20_wp ) &3814 ) + &3815 MIN( 0.0_wp, u(k,j,i) - u_gtrans ) * &3816 ABS( ( sk(k,j,i) - sk(k,j,i+1) ) /&3817 ( sk(k,j,i-1) - sk(k,j,i) + 1E-20_wp ) &3818 ) &3819 ) / ABS( u(k,j,i) - u_gtrans + 1E-20_wp )3820 3821 rr = ( MAX( 0.0_wp, u(k,j,i+1) - u_gtrans ) * &3822 ABS( ( sk(k,j,i) - sk(k,j,i-1) ) /&3823 ( sk(k,j,i+1) - sk(k,j,i) + 1E-20_wp ) &3824 ) + &3825 MIN( 0.0_wp, u(k,j,i+1) - u_gtrans ) * &3826 ABS( ( sk(k,j,i+1) - sk(k,j,i+2) ) /&3827 ( sk(k,j,i) - sk(k,j,i+1) + 1E-20_wp ) &3828 ) &3829 ) / ABS( u(k,j,i+1) - u_gtrans + 1E-20_wp )3830 3831 rs = ( MAX( 0.0_wp, v(k,j,i) - v_gtrans ) * &3832 ABS( ( sk(k,j-1,i) - sk(k,j-2,i) ) /&3833 ( sk(k,j,i) - sk(k,j-1,i) + 1E-20_wp ) &3834 ) + &3835 MIN( 0.0_wp, v(k,j,i) - v_gtrans ) * &3836 ABS( ( sk(k,j,i) - sk(k,j+1,i) ) /&3837 ( sk(k,j-1,i) - sk(k,j,i) + 1E-20_wp ) &3838 ) &3839 ) / ABS( v(k,j,i) - v_gtrans + 1E-20_wp )3840 3841 rn = ( MAX( 0.0_wp, v(k,j+1,i) - v_gtrans ) * &3842 ABS( ( sk(k,j,i) - sk(k,j-1,i) ) /&3843 ( sk(k,j+1,i) - sk(k,j,i) + 1E-20_wp ) &3844 ) + &3845 MIN( 0.0_wp, v(k,j+1,i) - v_gtrans ) * &3846 ABS( ( sk(k,j+1,i) - sk(k,j+2,i) ) /&3847 ( sk(k,j,i) - sk(k,j+1,i) + 1E-20_wp ) &3848 ) &3849 ) / ABS( v(k,j+1,i) - v_gtrans + 1E-20_wp )3850 !3851 !-- Reuse k_mm and compute k_mmm for the vertical gradient ratios.3852 !-- Note, for vertical advection below the third grid point above3853 !-- surface ( or below the model top) rd and rt are set to 0, i.e.3854 !-- use of first order scheme is enforced.3855 k_mmm = k - 3 * ibit83856 3857 rd = ( MAX( 0.0_wp, w(k-1,j,i) ) * &3858 ABS( ( sk(k_mm,j,i) - sk(k_mmm,j,i) ) / &3859 ( sk(k-1,j,i) - sk(k_mm,j,i) + 1E-20_wp ) &3860 ) + &3861 MIN( 0.0_wp, w(k-1,j,i) ) * &3862 ABS( ( sk(k-1,j,i) - sk(k,j,i) ) / &3863 ( sk(k_mm,j,i) - sk(k-1,j,i) + 1E-20_wp ) &3864 ) &3865 ) * ibit8 / ABS( w(k-1,j,i) + 1E-20_wp )3866 3867 rt = ( MAX( 0.0_wp, w(k,j,i) ) * &3868 ABS( ( sk(k,j,i) - sk(k-1,j,i) ) / &3869 ( sk(k+1,j,i) - sk(k,j,i) + 1E-20_wp ) &3870 ) + &3871 MIN( 0.0_wp, w(k,j,i) ) * &3872 ABS( ( sk(k+1,j,i) - sk(k_pp,j,i) ) / &3873 ( sk(k,j,i) - sk(k+1,j,i) + 1E-20_wp ) &3874 ) &3875 ) * ibit8 / ABS( w(k,j,i) + 1E-20_wp )3876 !3877 !-- Calculate empirical limiter function (van Albada2 limiter).3878 phi_l = MIN( 1.0_wp, ( 2.0_wp * ABS( rl ) ) / &3879 ( rl**2 + 1.0_wp ) )3880 phi_r = MIN( 1.0_wp, ( 2.0_wp * ABS( rr ) ) / &3881 ( rr**2 + 1.0_wp ) )3882 phi_s = MIN( 1.0_wp, ( 2.0_wp * ABS( rs ) ) / &3883 ( rs**2 + 1.0_wp ) )3884 phi_n = MIN( 1.0_wp, ( 2.0_wp * ABS( rn ) ) / &3885 ( rn**2 + 1.0_wp ) )3886 phi_d = MIN( 1.0_wp, ( 2.0_wp * ABS( rd ) ) / &3887 ( rd**2 + 1.0_wp ) )3888 phi_t = MIN( 1.0_wp, ( 2.0_wp * ABS( rt ) ) / &3889 ( rt**2 + 1.0_wp ) )3890 !3891 !-- Calculate the resulting monotone flux.3892 swap_flux_x_local(k,j) = fl_1 - phi_l * &3893 ( fl_1 - swap_flux_x_local(k,j) )3894 flux_r(k) = fr_1 - phi_r * &3895 ( fr_1 - flux_r(k) )3896 swap_flux_y_local(k) = fs_1 - phi_s * &3897 ( fs_1 - swap_flux_y_local(k) )3898 flux_n(k) = fn_1 - phi_n * &3899 ( fn_1 - flux_n(k) )3900 flux_d = fd_1 - phi_d * &3901 ( fd_1 - flux_d )3902 flux_t(k) = ft_1 - phi_t * &3903 ( ft_1 - flux_t(k) )3904 !3905 !-- Moreover, modify dissipation flux according to the limiter.3906 swap_diss_x_local(k,j) = swap_diss_x_local(k,j) * phi_l3907 diss_r(k) = diss_r(k) * phi_r3908 swap_diss_y_local(k) = swap_diss_y_local(k) * phi_s3909 diss_n(k) = diss_n(k) * phi_n3910 diss_d = diss_d * phi_d3911 diss_t(k) = diss_t(k) * phi_t3912 3913 ENDIF3914 3321 ! 3915 3322 !-- Calculate the divergence of the velocity field. A respective
Note: See TracChangeset
for help on using the changeset viewer.