Changeset 196 for trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2
- Timestamp:
- Jul 5, 2011, 12:50:06 AM (13 years ago)
- Location:
- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/dyn_em
- Files:
-
- 3 edited
-
module_advect_em.F (modified) (132 diffs)
-
module_em.F (modified) (27 diffs)
-
solve_em.F (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_advect_em.F
r11 r196 124 124 SUBROUTINE advect_u ( u, u_old, tendency, & 125 125 ru, rv, rom, & 126 mut, config_flags, & 127 msfu, msfv, msft, & 126 mut, time_step, config_flags, & 127 msfu , msfv , & 128 msft , & 128 129 fzm, fzp, & 129 130 rdx, rdy, rdzw, & … … 161 162 REAL , INTENT(IN ) :: rdx, & 162 163 rdy 164 INTEGER , INTENT(IN ) :: time_step 163 165 164 166 ! Local data … … 192 194 flux3(q_im2, q_im1, q_i, q_ip1, ua) = & 193 195 flux4(q_im2, q_im1, q_i, q_ip1, ua) + & 194 sign(1 .,ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0196 sign(1,time_step)*sign(1.,ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0 195 197 196 198 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = & … … 200 202 flux5(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = & 201 203 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) & 202 -sign(1 .,ua)*(&204 -sign(1,time_step)*sign(1.,ua)*( & 203 205 (q_ip2-q_im3)-5.*(q_ip1-q_im2)+10.*(q_i-q_im1) )/60.0 204 206 … … 227 229 ! of the higher order flux stencils 228 230 229 230 231 232 233 234 235 236 (its > ids+2) ) degrade_xs = .false.237 238 239 240 241 242 (jts > jds+2) ) degrade_ys = .false.243 244 245 (jte < jde-3) ) degrade_ye = .false.231 degrade_xs = .true. 232 degrade_xe = .true. 233 degrade_ys = .true. 234 degrade_ye = .true. 235 236 IF( config_flags%periodic_x .or. & 237 config_flags%symmetric_xs .or. & 238 (its > ids+3) ) degrade_xs = .false. 239 IF( config_flags%periodic_x .or. & 240 config_flags%symmetric_xe .or. & 241 (ite < ide-2) ) degrade_xe = .false. 242 IF( config_flags%periodic_y .or. & 243 config_flags%symmetric_ys .or. & 244 (jts > jds+3) ) degrade_ys = .false. 245 IF( config_flags%periodic_y .or. & 246 config_flags%symmetric_ye .or. & 247 (jte < jde-4) ) degrade_ye = .false. 246 248 247 249 !--------------- y - advection first … … 273 275 ENDIF 274 276 277 275 278 ! compute fluxes, 5th or 6th order 276 279 … … 280 283 j_loop_y_flux_6 : DO j = j_start, j_end+1 281 284 282 283 284 285 286 287 fqy( i, k, jp1 ) = vel*flux6(&288 289 290 291 285 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN ! use full stencil 286 287 DO k=kts,ktf 288 DO i = i_start, i_end 289 vel = 0.5*(rv(i,k,j)+rv(i-1,k,j)) 290 fqy( i, k, jp1 ) = vel*flux6( & 291 u(i,k,j-3), u(i,k,j-2), u(i,k,j-1), & 292 u(i,k,j ), u(i,k,j+1), u(i,k,j+2), vel ) 293 ENDDO 294 ENDDO 292 295 293 296 ! we must be close to some boundary where we need to reduce the order of the stencil 294 297 295 296 297 DO k=kts,ktf298 DO i = i_start, i_end298 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary 299 300 DO k=kts,ktf 301 DO i = i_start, i_end 299 302 fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i-1,k,j)) & 300 *(u(i,k,j)+u(i,k,j-1))301 ENDDO302 ENDDO303 304 305 306 DO k=kts,ktf307 DO i = i_start, i_end303 *(u(i,k,j)+u(i,k,j-1)) 304 ENDDO 305 ENDDO 306 307 ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary 308 309 DO k=kts,ktf 310 DO i = i_start, i_end 308 311 vel = 0.5*(rv(i,k,j)+rv(i-1,k,j)) 309 312 fqy( i, k, jp1 ) = vel*flux4( & 310 313 u(i,k,j-2),u(i,k,j-1), u(i,k,j),u(i,k,j+1),vel ) 311 ENDDO312 ENDDO313 314 315 316 DO k=kts,ktf317 DO i = i_start, i_end314 ENDDO 315 ENDDO 316 317 ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary 318 319 DO k=kts,ktf 320 DO i = i_start, i_end 318 321 fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i-1,k,j)) & 319 322 *(u(i,k,j)+u(i,k,j-1)) 320 ENDDO321 ENDDO322 323 ELSE IF ( j == jde-2 ) THEN ! 3rd or 4thorder flux 2 in from north boundary324 325 DO k=kts,ktf326 DO i = i_start, i_end323 ENDDO 324 ENDDO 325 326 ELSE IF ( j == jde-2 ) THEN ! 3rd order flux 2 in from north boundary 327 328 DO k=kts,ktf 329 DO i = i_start, i_end 327 330 vel = 0.5*(rv(i,k,j)+rv(i-1,k,j)) 328 331 fqy( i, k, jp1 ) = vel*flux4( & 329 332 u(i,k,j-2),u(i,k,j-1), & 330 333 u(i,k,j),u(i,k,j+1),vel ) 331 ENDDO 332 ENDDO 333 334 END IF 335 336 !stopped 334 ENDDO 335 ENDDO 336 337 END IF 337 338 338 339 ! y flux-divergence into tendency … … 348 349 349 350 ENDIF 351 350 352 351 353 … … 475 477 IF( config_flags%periodic_x .or. & 476 478 config_flags%symmetric_xs .or. & 477 (its > ids+ 2) ) degrade_xs = .false.479 (its > ids+3) ) degrade_xs = .false. 478 480 IF( config_flags%periodic_x .or. & 479 481 config_flags%symmetric_xe .or. & … … 481 483 IF( config_flags%periodic_y .or. & 482 484 config_flags%symmetric_ys .or. & 483 (jts > jds+ 2) ) degrade_ys = .false.485 (jts > jds+3) ) degrade_ys = .false. 484 486 IF( config_flags%periodic_y .or. & 485 487 config_flags%symmetric_ye .or. & 486 (jte < jde- 3) ) degrade_ye = .false.488 (jte < jde-4) ) degrade_ye = .false. 487 489 488 490 !--------------- y - advection first … … 514 516 ENDIF 515 517 518 516 519 ! compute fluxes, 5th or 6th order 517 520 … … 562 565 ENDDO 563 566 564 ELSE IF ( j == jde-2 ) THEN ! 3rd or 4th order flux 2 in from north boundary567 ELSE IF ( j == jde-2 ) THEN ! 3rd order flux 2 in from north boundary 565 568 566 569 DO k=kts,ktf … … 587 590 588 591 ENDIF 592 589 593 590 594 … … 709 713 IF( config_flags%periodic_x .or. & 710 714 config_flags%symmetric_xs .or. & 711 (its > ids+ 1) ) degrade_xs = .false.715 (its > ids+2) ) degrade_xs = .false. 712 716 IF( config_flags%periodic_x .or. & 713 717 config_flags%symmetric_xe .or. & … … 715 719 IF( config_flags%periodic_y .or. & 716 720 config_flags%symmetric_ys .or. & 717 (jts > jds+ 1) ) degrade_ys = .false.721 (jts > jds+2) ) degrade_ys = .false. 718 722 IF( config_flags%periodic_y .or. & 719 723 config_flags%symmetric_ye .or. & 720 (jte < jde- 2) ) degrade_ye = .false.724 (jte < jde-3) ) degrade_ye = .false. 721 725 722 726 !--------------- x - advection first … … 817 821 j_end_f = jde-2 818 822 ENDIF 823 819 824 820 825 ! j flux loop for v flux of u momentum … … 835 840 DO k = kts, ktf 836 841 DO i = i_start, i_end 837 fqy(i, k, jp1) = 0.25*(rv(i,k,j_end+1)+rv(i-1,k,j_end+1)) & 838 *(u(i,k,j_end+1)+u(i,k,j_end)) 842 ! Assumes j>j_end_f is ONLY j_end+1 ... 843 ! fqy(i, k, jp1) = 0.25*(rv(i,k,j_end+1)+rv(i-1,k,j_end+1)) & 844 ! *(u(i,k,j_end+1)+u(i,k,j_end)) 845 fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i-1,k,j)) & 846 *(u(i,k,j)+u(i,k,j-1)) 839 847 ENDDO 840 848 ENDDO … … 852 860 END IF 853 861 862 ! y flux-divergence into tendency 863 854 864 IF (j > j_start) THEN 855 856 ! y flux-divergence into tendency857 865 858 866 DO k=kts,ktf … … 865 873 END IF 866 874 875 867 876 jtmp = jp1 868 877 jp1 = jp0 … … 890 899 IF( config_flags%periodic_x .or. & 891 900 config_flags%symmetric_xs .or. & 892 (its > ids+ 1) ) degrade_xs = .false.901 (its > ids+2) ) degrade_xs = .false. 893 902 IF( config_flags%periodic_x .or. & 894 903 config_flags%symmetric_xe .or. & … … 896 905 IF( config_flags%periodic_y .or. & 897 906 config_flags%symmetric_ys .or. & 898 (jts > jds+ 1) ) degrade_ys = .false.907 (jts > jds+2) ) degrade_ys = .false. 899 908 IF( config_flags%periodic_y .or. & 900 909 config_flags%symmetric_ye .or. & 901 (jte < jde- 2) ) degrade_ye = .false.910 (jte < jde-3) ) degrade_ye = .false. 902 911 903 912 !--------------- x - advection first … … 997 1006 j_end_f = jde-2 998 1007 ENDIF 1008 999 1009 1000 1010 ! j flux loop for v flux of u momentum … … 1015 1025 DO k = kts, ktf 1016 1026 DO i = i_start, i_end 1017 fqy(i, k, jp1) = 0.25*(rv(i,k,j_end+1)+rv(i-1,k,j_end+1)) & 1018 *(u(i,k,j_end+1)+u(i,k,j_end)) 1027 ! Assumes j>j_end_f is ONLY j_end+1 ... 1028 ! fqy(i, k, jp1) = 0.25*(rv(i,k,j_end+1)+rv(i-1,k,j_end+1)) & 1029 ! *(u(i,k,j_end+1)+u(i,k,j_end)) 1030 fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i-1,k,j)) & 1031 *(u(i,k,j)+u(i,k,j-1)) 1019 1032 ENDDO 1020 1033 ENDDO … … 1032 1045 END IF 1033 1046 1047 ! y flux-divergence into tendency 1048 1034 1049 IF (j > j_start) THEN 1035 1036 ! y flux-divergence into tendency1037 1050 1038 1051 DO k=kts,ktf … … 1044 1057 1045 1058 END IF 1059 1046 1060 1047 1061 jtmp = jp1 … … 1116 1130 ENDDO 1117 1131 ENDDO 1132 1133 ELSE IF ( horz_order == 0 ) THEN 1134 1135 ! Just in case we want to turn horizontal advection off, we can do it 1118 1136 1119 1137 ELSE … … 1410 1428 SUBROUTINE advect_v ( v, v_old, tendency, & 1411 1429 ru, rv, rom, & 1412 mut, config_flags,&1430 mut, time_step, config_flags, & 1413 1431 msfu, msfv, msft, & 1414 1432 fzm, fzp, & … … 1447 1465 REAL , INTENT(IN ) :: rdx, & 1448 1466 rdy 1467 INTEGER , INTENT(IN ) :: time_step 1468 1449 1469 1450 1470 ! Local data … … 1481 1501 flux3(q_im2, q_im1, q_i, q_ip1, ua) = & 1482 1502 flux4(q_im2, q_im1, q_i, q_ip1, ua) + & 1483 sign(1 .,ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.01503 sign(1,time_step)*sign(1.,ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0 1484 1504 1485 1505 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = & … … 1489 1509 flux5(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = & 1490 1510 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) & 1491 -sign(1 .,ua)*(&1511 -sign(1,time_step)*sign(1.,ua)*( & 1492 1512 (q_ip2-q_im3)-5.*(q_ip1-q_im2)+10.*(q_i-q_im1) )/60.0 1493 1513 … … 1518 1538 ! of the higher order flux stencils 1519 1539 1520 1521 1522 1523 1524 1525 1526 1527 (its > ids+2) ) degrade_xs = .false.1528 1529 1530 1531 1532 1533 (jts > jds+2) ) degrade_ys = .false.1534 1535 1536 (jte < jde-2) ) degrade_ye = .false.1540 degrade_xs = .true. 1541 degrade_xe = .true. 1542 degrade_ys = .true. 1543 degrade_ye = .true. 1544 1545 IF( config_flags%periodic_x .or. & 1546 config_flags%symmetric_xs .or. & 1547 (its > ids+3) ) degrade_xs = .false. 1548 IF( config_flags%periodic_x .or. & 1549 config_flags%symmetric_xe .or. & 1550 (ite < ide-3) ) degrade_xe = .false. 1551 IF( config_flags%periodic_y .or. & 1552 config_flags%symmetric_ys .or. & 1553 (jts > jds+3) ) degrade_ys = .false. 1554 IF( config_flags%periodic_y .or. & 1555 config_flags%symmetric_ye .or. & 1556 (jte < jde-3) ) degrade_ye = .false. 1537 1557 1538 1558 !--------------- y - advection first 1539 1540 ktf=MIN(kte,kde-1)1541 1559 1542 1560 i_start = its … … 1552 1570 1553 1571 IF(degrade_ys) then 1554 1555 1572 j_start = MAX(jts,jds+1) 1573 j_start_f = jds+3 1556 1574 ENDIF 1557 1575 1558 1576 IF(degrade_ye) then 1559 1560 1577 j_end = MIN(jte,jde-1) 1578 j_end_f = jde-2 1561 1579 ENDIF 1562 1580 1563 1581 ! compute fluxes, 5th or 6th order 1564 1582 1565 1566 1567 1568 1569 1570 1571 1572 1573 1574 1575 fqy( i, k, jp1 ) = vel*flux6(&1576 1577 1578 1579 1583 jp1 = 2 1584 jp0 = 1 1585 1586 j_loop_y_flux_6 : DO j = j_start, j_end+1 1587 1588 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN 1589 1590 DO k=kts,ktf 1591 DO i = i_start, i_end 1592 vel = 0.5*(rv(i,k,j)+rv(i,k,j-1)) 1593 fqy( i, k, jp1 ) = vel*flux6( & 1594 v(i,k,j-3), v(i,k,j-2), v(i,k,j-1), & 1595 v(i,k,j ), v(i,k,j+1), v(i,k,j+2), vel ) 1596 ENDDO 1597 ENDDO 1580 1598 1581 1599 ! we must be close to some boundary where we need to reduce the order of the stencil 1582 1600 ! specified uses upstream normal wind at boundaries 1583 1601 1584 1602 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary 1585 1603 1586 1604 DO k=kts,ktf … … 1593 1611 ENDDO 1594 1612 1595 1613 ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary 1596 1614 1597 1615 DO k=kts,ktf 1598 1616 DO i = i_start, i_end 1599 1600 1601 1602 ENDDO 1603 ENDDO 1604 1605 1606 1617 vel = 0.5*(rv(i,k,j)+rv(i,k,j-1)) 1618 fqy( i, k, jp1 ) = vel*flux4( & 1619 v(i,k,j-2),v(i,k,j-1),v(i,k,j),v(i,k,j+1),vel ) 1620 ENDDO 1621 ENDDO 1622 1623 1624 ELSE IF ( j == jde ) THEN ! 2nd order flux next to north boundary 1607 1625 1608 1626 DO k=kts,ktf … … 1615 1633 ENDDO 1616 1634 1617 1635 ELSE IF ( j == jde-1 ) THEN ! 3rd or 4th order flux 2 in from north boundary 1618 1636 1619 1637 DO k=kts,ktf 1620 1638 DO i = i_start, i_end 1621 1622 1623 1624 ENDDO 1625 ENDDO 1626 1627 1639 vel = 0.5*(rv(i,k,j)+rv(i,k,j-1)) 1640 fqy( i, k, jp1 ) = vel*flux4( & 1641 v(i,k,j-2),v(i,k,j-1),v(i,k,j),v(i,k,j+1),vel ) 1642 ENDDO 1643 ENDDO 1644 1645 END IF 1628 1646 1629 1647 ! y flux-divergence into tendency 1630 1648 1631 IF(j > j_start) THEN 1632 1633 DO k=kts,ktf 1634 DO i = i_start, i_end 1635 mrdy=msfv(i,j-1)*rdy 1636 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0)) 1637 ENDDO 1638 ENDDO 1639 1640 ENDIF 1641 1642 jtmp = jp1 1643 jp1 = jp0 1644 jp0 = jtmp 1645 1646 ENDDO j_loop_y_flux_6 1649 1650 IF(j > j_start) THEN 1651 1652 DO k=kts,ktf 1653 DO i = i_start, i_end 1654 mrdy=msfv(i,j-1)*rdy 1655 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0)) 1656 ENDDO 1657 ENDDO 1658 1659 ENDIF 1660 1661 1662 1663 jtmp = jp1 1664 jp1 = jp0 1665 jp0 = jtmp 1666 1667 ENDDO j_loop_y_flux_6 1647 1668 1648 1669 ! next, x - flux divergence … … 1653 1674 j_start = jts 1654 1675 j_end = jte 1676 1655 1677 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts) 1656 1678 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-1,jte) … … 1663 1685 1664 1686 IF(degrade_xs) then 1665 i_start = MAX(ids+1,its) 1666 i_start_f = i_start+2 1687 i_start = MAX(ids+1,its) 1688 ! i_start_f = i_start+2 1689 i_start_f = MIN(i_start+2,ids+3) 1667 1690 ENDIF 1668 1691 1669 1692 IF(degrade_xe) then 1670 1671 1693 i_end = MIN(ide-2,ite) 1694 i_end_f = ide-3 1672 1695 ENDIF 1673 1696 … … 1678 1701 ! 5th or 6th order flux 1679 1702 1680 1681 1682 1683 1684 v(i-1,k,j), v(i ,k,j), &1685 v(i+1,k,j), v(i+2,k,j), &1686 vel )1687 1688 1703 DO k=kts,ktf 1704 DO i = i_start_f, i_end_f 1705 vel = 0.5*(ru(i,k,j)+ru(i,k,j-1)) 1706 fqx( i, k ) = vel*flux6( v(i-3,k,j), v(i-2,k,j), & 1707 v(i-1,k,j), v(i ,k,j), & 1708 v(i+1,k,j), v(i+2,k,j), & 1709 vel ) 1710 ENDDO 1711 ENDDO 1689 1712 1690 1713 ! lower order fluxes close to boundaries (if not periodic or symmetric) 1691 1714 1692 IF( degrade_xs ) THEN 1693 1694 IF( i_start == ids+1 ) THEN ! second order flux next to the boundary 1695 i = ids+1 1696 DO k=kts,ktf 1697 fqx(i,k) = 0.25*(ru(i,k,j)+ru(i,k,j-1)) & 1698 *(v(i,k,j)+v(i-1,k,j)) 1699 ENDDO 1715 IF( degrade_xs ) THEN 1716 1717 DO i=i_start,i_start_f-1 1718 1719 IF(i == ids+1) THEN ! second order 1720 DO k=kts,ktf 1721 fqx(i,k) = 0.25*(ru(i,k,j)+ru(i,k,j-1)) & 1722 *(v(i,k,j)+v(i-1,k,j)) 1723 ENDDO 1700 1724 ENDIF 1701 1725 1702 i = ids+2 1703 DO k=kts,ktf 1704 vel = 0.5*(ru(i,k,j)+ru(i,k,j-1)) 1705 fqx( i,k ) = vel*flux4( v(i-2,k,j), v(i-1,k,j), & 1706 v(i ,k,j), v(i+1,k,j), & 1707 vel ) 1708 ENDDO 1709 1710 ENDIF 1711 1712 IF( degrade_xe ) THEN 1713 1714 IF( i_end == ide-2 ) THEN ! second order flux next to the boundary 1715 i = ide-1 1716 DO k=kts,ktf 1717 fqx(i,k) = 0.25*(ru(i_end+1,k,j)+ru(i_end+1,k,j-1)) & 1718 *(v(i_end+1,k,j)+v(i_end,k,j)) 1719 ENDDO 1726 IF(i == ids+2) THEN ! third order 1727 DO k=kts,ktf 1728 vel = 0.5*(ru(i,k,j)+ru(i,k,j-1)) 1729 fqx( i,k ) = vel*flux4( v(i-2,k,j), v(i-1,k,j), & 1730 v(i ,k,j), v(i+1,k,j), & 1731 vel ) 1732 ENDDO 1720 1733 ENDIF 1721 1734 1722 i = ide-2 1723 DO k=kts,ktf 1724 vel = 0.5*(ru(i,k,j)+ru(i,k,j-1)) 1725 fqx( i,k ) = vel*flux4( v(i-2,k,j), v(i-1,k,j), & 1726 v(i ,k,j), v(i+1,k,j), & 1727 vel ) 1728 ENDDO 1729 1730 ENDIF 1735 ENDDO 1736 1737 ENDIF 1738 1739 IF( degrade_xe ) THEN 1740 1741 DO i = i_end_f+1, i_end+1 1742 1743 IF( i == ide-1 ) THEN ! second order flux next to the boundary 1744 DO k=kts,ktf 1745 fqx(i,k) = 0.25*(ru(i_end+1,k,j)+ru(i_end+1,k,j-1)) & 1746 *(v(i_end+1,k,j)+v(i_end,k,j)) 1747 ENDDO 1748 ENDIF 1749 1750 IF( i == ide-2 ) THEN ! third order flux one in from the boundary 1751 DO k=kts,ktf 1752 vel = 0.5*(ru(i,k,j)+ru(i,k,j-1)) 1753 fqx( i,k ) = vel*flux4( v(i-2,k,j), v(i-1,k,j), & 1754 v(i ,k,j), v(i+1,k,j), & 1755 vel ) 1756 ENDDO 1757 ENDIF 1758 1759 ENDDO 1760 1761 ENDIF 1731 1762 1732 1763 ! x flux-divergence into tendency 1733 1764 1734 1735 1765 DO k=kts,ktf 1766 DO i = i_start, i_end 1736 1767 mrdx=msfv(i,j)*rdx 1737 1768 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k)) 1738 ENDDO1739 ENDDO1740 1741 ENDDO1769 ENDDO 1770 ENDDO 1771 1772 ENDDO 1742 1773 1743 1774 ELSE IF( horz_order == 5 ) THEN … … 1762 1793 IF( config_flags%periodic_x .or. & 1763 1794 config_flags%symmetric_xs .or. & 1764 (its > ids+ 2) ) degrade_xs = .false.1795 (its > ids+3) ) degrade_xs = .false. 1765 1796 IF( config_flags%periodic_x .or. & 1766 1797 config_flags%symmetric_xe .or. & … … 1768 1799 IF( config_flags%periodic_y .or. & 1769 1800 config_flags%symmetric_ys .or. & 1770 (jts > jds+ 2) ) degrade_ys = .false.1801 (jts > jds+3) ) degrade_ys = .false. 1771 1802 IF( config_flags%periodic_y .or. & 1772 1803 config_flags%symmetric_ye .or. & 1773 (jte < jde- 2) ) degrade_ye = .false.1804 (jte < jde-3) ) degrade_ye = .false. 1774 1805 1775 1806 !--------------- y - advection first … … 1875 1906 ENDIF 1876 1907 1908 1877 1909 jtmp = jp1 1878 1910 jp1 = jp0 … … 1899 1931 IF(degrade_xs) then 1900 1932 i_start = MAX(ids+1,its) 1901 i_start_f = i_start+2 1933 ! i_start_f = i_start+2 1934 i_start_f = MIN(i_start+2,ids+3) 1902 1935 ENDIF 1903 1936 … … 1927 1960 IF( degrade_xs ) THEN 1928 1961 1929 IF( i_start == ids+1 ) THEN ! second order flux next to the boundary 1930 i = ids+1 1931 DO k=kts,ktf 1932 fqx(i,k) = 0.25*(ru(i,k,j)+ru(i,k,j-1)) & 1933 *(v(i,k,j)+v(i-1,k,j)) 1934 ENDDO 1935 ENDIF 1936 1937 i = ids+2 1938 DO k=kts,ktf 1939 vel = 0.5*(ru(i,k,j)+ru(i,k,j-1)) 1940 fqx( i,k ) = vel*flux3( v(i-2,k,j), v(i-1,k,j), & 1941 v(i ,k,j), v(i+1,k,j), & 1942 vel ) 1943 ENDDO 1944 1945 ENDIF 1946 1947 IF( degrade_xe ) THEN 1948 1949 IF( i_end == ide-2 ) THEN ! second order flux next to the boundary 1950 i = ide-1 1951 DO k=kts,ktf 1952 fqx(i,k) = 0.25*(ru(i_end+1,k,j)+ru(i_end+1,k,j-1)) & 1953 *(v(i_end+1,k,j)+v(i_end,k,j)) 1954 ENDDO 1955 ENDIF 1956 1957 i = ide-2 1958 DO k=kts,ktf 1959 vel = 0.5*(ru(i,k,j)+ru(i,k,j-1)) 1960 fqx( i,k ) = vel*flux3( v(i-2,k,j), v(i-1,k,j), & 1962 DO i=i_start,i_start_f-1 1963 1964 IF(i == ids+1) THEN ! second order 1965 DO k=kts,ktf 1966 fqx(i,k) = 0.25*(ru(i,k,j)+ru(i,k,j-1)) & 1967 *(v(i,k,j)+v(i-1,k,j)) 1968 ENDDO 1969 ENDIF 1970 1971 IF(i == ids+2) THEN ! third order 1972 DO k=kts,ktf 1973 vel = 0.5*(ru(i,k,j)+ru(i,k,j-1)) 1974 fqx( i,k ) = vel*flux3( v(i-2,k,j), v(i-1,k,j), & 1961 1975 v(i ,k,j), v(i+1,k,j), & 1962 1976 vel ) 1977 ENDDO 1978 ENDIF 1979 1980 ENDDO 1981 1982 ENDIF 1983 1984 IF( degrade_xe ) THEN 1985 1986 DO i = i_end_f+1, i_end+1 1987 1988 IF( i == ide-1 ) THEN ! second order flux next to the boundary 1989 DO k=kts,ktf 1990 fqx(i,k) = 0.25*(ru(i_end+1,k,j)+ru(i_end+1,k,j-1)) & 1991 *(v(i_end+1,k,j)+v(i_end,k,j)) 1992 ENDDO 1993 ENDIF 1994 1995 IF( i == ide-2 ) THEN ! third order flux one in from the boundary 1996 DO k=kts,ktf 1997 vel = 0.5*(ru(i,k,j)+ru(i,k,j-1)) 1998 fqx( i,k ) = vel*flux3( v(i-2,k,j), v(i-1,k,j), & 1999 v(i ,k,j), v(i+1,k,j), & 2000 vel ) 2001 ENDDO 2002 ENDIF 2003 1963 2004 ENDDO 1964 2005 … … 1992 2033 IF( config_flags%periodic_x .or. & 1993 2034 config_flags%symmetric_xs .or. & 1994 (its > ids+ 1) ) degrade_xs = .false.2035 (its > ids+2) ) degrade_xs = .false. 1995 2036 IF( config_flags%periodic_x .or. & 1996 2037 config_flags%symmetric_xe .or. & … … 1998 2039 IF( config_flags%periodic_y .or. & 1999 2040 config_flags%symmetric_ys .or. & 2000 (jts > jds+ 1) ) degrade_ys = .false.2041 (jts > jds+2) ) degrade_ys = .false. 2001 2042 IF( config_flags%periodic_y .or. & 2002 2043 config_flags%symmetric_ye .or. & 2003 (jte < jde- 1) ) degrade_ye = .false.2044 (jte < jde-2) ) degrade_ye = .false. 2004 2045 2005 2046 !--------------- y - advection first … … 2074 2115 ENDDO 2075 2116 ENDDO 2117 2118 2076 2119 END IF 2077 2120 … … 2167 2210 IF( config_flags%periodic_x .or. & 2168 2211 config_flags%symmetric_xs .or. & 2169 (its > ids+ 1) ) degrade_xs = .false.2212 (its > ids+2) ) degrade_xs = .false. 2170 2213 IF( config_flags%periodic_x .or. & 2171 2214 config_flags%symmetric_xe .or. & … … 2173 2216 IF( config_flags%periodic_y .or. & 2174 2217 config_flags%symmetric_ys .or. & 2175 (jts > jds+ 1) ) degrade_ys = .false.2218 (jts > jds+2) ) degrade_ys = .false. 2176 2219 IF( config_flags%periodic_y .or. & 2177 2220 config_flags%symmetric_ye .or. & 2178 (jte < jde- 1) ) degrade_ye = .false.2221 (jte < jde-2) ) degrade_ye = .false. 2179 2222 2180 2223 !--------------- y - advection first … … 2249 2292 ENDDO 2250 2293 ENDDO 2294 2295 2251 2296 END IF 2252 2297 … … 2405 2450 ENDDO 2406 2451 ENDDO 2452 2453 ELSE IF ( horz_order == 0 ) THEN 2454 2455 ! Just in case we want to turn horizontal advection off, we can do it 2407 2456 2408 2457 ELSE … … 2530 2579 ENDDO 2531 2580 2581 ! 2582 ! 2532 2583 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts) 2533 2584 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-1,jte) … … 2570 2621 DO k=kts,ktf 2571 2622 DO i = i_start, i_end 2623 ! 2624 ! 2572 2625 tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k)) 2573 2626 ENDDO … … 2613 2666 DO k=kts,ktf 2614 2667 DO i = i_start, i_end 2668 ! 2669 ! 2615 2670 tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k)) 2616 2671 ENDDO … … 2646 2701 DO k=kts,ktf 2647 2702 DO i = i_start, i_end 2703 ! 2704 ! 2648 2705 tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k)) 2649 2706 ENDDO … … 2679 2736 DO k=kts,ktf 2680 2737 DO i = i_start, i_end 2738 ! 2739 ! 2681 2740 tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k)) 2682 2741 ENDDO … … 2699 2758 DO k=kts,ktf 2700 2759 DO i = i_start, i_end 2760 ! 2761 ! 2701 2762 tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k)) 2702 2703 2763 ENDDO 2704 2764 ENDDO … … 2718 2778 SUBROUTINE advect_scalar ( field, field_old, tendency, & 2719 2779 ru, rv, rom, & 2720 mut, config_flags,&2780 mut, time_step, config_flags, & 2721 2781 msfu, msfv, msft, & 2722 2782 fzm, fzp, & … … 2755 2815 REAL , INTENT(IN ) :: rdx, & 2756 2816 rdy 2817 INTEGER , INTENT(IN ) :: time_step 2818 2757 2819 2758 2820 ! Local data … … 2788 2850 flux3(q_im2, q_im1, q_i, q_ip1, ua) = & 2789 2851 flux4(q_im2, q_im1, q_i, q_ip1, ua) + & 2790 sign(1 .,ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.02852 sign(1,time_step)*sign(1.,ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0 2791 2853 2792 2854 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = & … … 2796 2858 flux5(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = & 2797 2859 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) & 2798 -sign(1 .,ua)*(&2860 -sign(1,time_step)*sign(1.,ua)*( & 2799 2861 (q_ip2-q_im3)-5.*(q_ip1-q_im2)+10.*(q_i-q_im1) )/60.0 2800 2862 … … 2831 2893 IF( config_flags%periodic_x .or. & 2832 2894 config_flags%symmetric_xs .or. & 2833 (its > ids+ 2) ) degrade_xs = .false.2895 (its > ids+3) ) degrade_xs = .false. 2834 2896 IF( config_flags%periodic_x .or. & 2835 2897 config_flags%symmetric_xe .or. & … … 2837 2899 IF( config_flags%periodic_y .or. & 2838 2900 config_flags%symmetric_ys .or. & 2839 (jts > jds+ 2) ) degrade_ys = .false.2901 (jts > jds+3) ) degrade_ys = .false. 2840 2902 IF( config_flags%periodic_y .or. & 2841 2903 config_flags%symmetric_ye .or. & 2842 (jte < jde- 3) ) degrade_ye = .false.2904 (jte < jde-4) ) degrade_ye = .false. 2843 2905 2844 2906 !--------------- y - advection first … … 2865 2927 j_end_f = jde-3 2866 2928 ENDIF 2929 2867 2930 2868 2931 ! compute fluxes, 5th or 6th order … … 2884 2947 ENDDO 2885 2948 2949 2886 2950 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary 2887 2951 … … 2894 2958 ENDDO 2895 2959 2896 ELSE IF ( j == jds+2 ) THEN ! third of4th order flux 2 in from south boundary2960 ELSE IF ( j == jds+2 ) THEN ! 4th order flux 2 in from south boundary 2897 2961 2898 2962 DO k=kts,ktf … … 2939 3003 ENDIF 2940 3004 3005 2941 3006 jtmp = jp1 2942 3007 jp1 = jp0 … … 2961 3026 IF(degrade_xs) then 2962 3027 i_start = MAX(ids+1,its) 2963 i_start_f = i_start+2 3028 ! i_start_f = i_start+2 3029 i_start_f = MIN(i_start+2,ids+3) 2964 3030 ENDIF 2965 3031 … … 2989 3055 IF( degrade_xs ) THEN 2990 3056 2991 IF( i_start == ids+1 ) THEN ! second order flux next to the boundary 2992 i = ids+1 2993 DO k=kts,ktf 2994 fqx(i,k) = 0.5*(ru(i,k,j)) & 2995 *(field(i,k,j)+field(i-1,k,j)) 2996 2997 ENDDO 2998 ENDIF 2999 3000 i = ids+2 3001 DO k=kts,ktf 3002 vel = ru(i,k,j) 3003 fqx( i,k ) = vel*flux4( field(i-2,k,j), field(i-1,k,j), & 3004 field(i ,k,j), field(i+1,k,j), & 3005 vel ) 3057 DO i=i_start,i_start_f-1 3058 3059 IF(i == ids+1) THEN ! second order 3060 DO k=kts,ktf 3061 fqx(i,k) = 0.5*(ru(i,k,j)) & 3062 *(field(i,k,j)+field(i-1,k,j)) 3063 ENDDO 3064 ENDIF 3065 3066 IF(i == ids+2) THEN ! third order 3067 DO k=kts,ktf 3068 vel = ru(i,k,j) 3069 fqx( i,k ) = vel*flux4( field(i-2,k,j), field(i-1,k,j), & 3070 field(i ,k,j), field(i+1,k,j), & 3071 vel ) 3072 ENDDO 3073 END IF 3074 3006 3075 ENDDO 3007 3076 … … 3010 3079 IF( degrade_xe ) THEN 3011 3080 3012 IF( i_end == ide-2 ) THEN ! second order flux next to the boundary 3013 i = ide-1 3014 DO k=kts,ktf 3015 fqx(i,k) = 0.5*(ru(i,k,j)) & 3016 *(field(i,k,j)+field(i-1,k,j)) 3017 ENDDO 3018 ENDIF 3019 3020 i = ide-2 3021 DO k=kts,ktf 3022 vel = ru(i,k,j) 3023 fqx( i,k ) = vel*flux4( field(i-2,k,j), field(i-1,k,j), & 3024 field(i ,k,j), field(i+1,k,j), & 3025 vel ) 3026 ENDDO 3027 3028 ENDIF 3081 DO i = i_end_f+1, i_end+1 3082 3083 IF( i == ide-1 ) THEN ! second order flux next to the boundary 3084 DO k=kts,ktf 3085 fqx(i,k) = 0.5*(ru(i,k,j)) & 3086 *(field(i,k,j)+field(i-1,k,j)) 3087 ENDDO 3088 ENDIF 3089 3090 IF( i == ide-2 ) THEN ! third order flux one in from the boundary 3091 DO k=kts,ktf 3092 vel = ru(i,k,j) 3093 fqx( i,k ) = vel*flux4( field(i-2,k,j), field(i-1,k,j), & 3094 field(i ,k,j), field(i+1,k,j), & 3095 vel ) 3096 ENDDO 3097 ENDIF 3098 3099 ENDDO 3100 3101 ENDIF 3029 3102 3030 3103 ! x flux-divergence into tendency … … 3055 3128 IF( config_flags%periodic_x .or. & 3056 3129 config_flags%symmetric_xs .or. & 3057 (its > ids+ 2) ) degrade_xs = .false.3130 (its > ids+3) ) degrade_xs = .false. 3058 3131 IF( config_flags%periodic_x .or. & 3059 3132 config_flags%symmetric_xe .or. & … … 3061 3134 IF( config_flags%periodic_y .or. & 3062 3135 config_flags%symmetric_ys .or. & 3063 (jts > jds+ 2) ) degrade_ys = .false.3136 (jts > jds+3) ) degrade_ys = .false. 3064 3137 IF( config_flags%periodic_y .or. & 3065 3138 config_flags%symmetric_ye .or. & 3066 (jte < jde- 3) ) degrade_ye = .false.3139 (jte < jde-4) ) degrade_ye = .false. 3067 3140 3068 3141 !--------------- y - advection first … … 3089 3162 j_end_f = jde-3 3090 3163 ENDIF 3164 3091 3165 3092 3166 ! compute fluxes, 5th or 6th order … … 3108 3182 ENDDO 3109 3183 3184 3110 3185 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary 3111 3186 … … 3186 3261 IF(degrade_xs) then 3187 3262 i_start = MAX(ids+1,its) 3188 i_start_f = i_start+2 3263 ! i_start_f = i_start+2 3264 i_start_f = MIN(i_start+2,ids+3) 3189 3265 ENDIF 3190 3266 … … 3214 3290 IF( degrade_xs ) THEN 3215 3291 3216 IF( i_start == ids+1 ) THEN ! second order flux next to the boundary 3217 i = ids+1 3218 DO k=kts,ktf 3219 fqx(i,k) = 0.5*(ru(i,k,j)) & 3220 *(field(i,k,j)+field(i-1,k,j)) 3221 3222 ENDDO 3223 ENDIF 3224 3225 i = ids+2 3226 DO k=kts,ktf 3227 vel = ru(i,k,j) 3228 fqx( i,k ) = vel*flux3( field(i-2,k,j), field(i-1,k,j), & 3229 field(i ,k,j), field(i+1,k,j), & 3230 vel ) 3292 DO i=i_start,i_start_f-1 3293 3294 IF(i == ids+1) THEN ! second order 3295 DO k=kts,ktf 3296 fqx(i,k) = 0.5*(ru(i,k,j)) & 3297 *(field(i,k,j)+field(i-1,k,j)) 3298 ENDDO 3299 ENDIF 3300 3301 IF(i == ids+2) THEN ! third order 3302 DO k=kts,ktf 3303 vel = ru(i,k,j) 3304 fqx( i,k ) = vel*flux3( field(i-2,k,j), field(i-1,k,j), & 3305 field(i ,k,j), field(i+1,k,j), & 3306 vel ) 3307 ENDDO 3308 END IF 3309 3231 3310 ENDDO 3232 3311 … … 3235 3314 IF( degrade_xe ) THEN 3236 3315 3237 IF( i_end == ide-2 ) THEN ! second order flux next to the boundary 3238 i = ide-1 3239 DO k=kts,ktf 3240 fqx(i,k) = 0.5*(ru(i,k,j)) & 3241 *(field(i,k,j)+field(i-1,k,j)) 3242 ENDDO 3243 ENDIF 3244 3245 i = ide-2 3246 DO k=kts,ktf 3247 vel = ru(i,k,j) 3248 fqx( i,k ) = vel*flux3( field(i-2,k,j), field(i-1,k,j), & 3249 field(i ,k,j), field(i+1,k,j), & 3250 vel ) 3251 ENDDO 3252 3253 ENDIF 3316 DO i = i_end_f+1, i_end+1 3317 3318 IF( i == ide-1 ) THEN ! second order flux next to the boundary 3319 DO k=kts,ktf 3320 fqx(i,k) = 0.5*(ru(i,k,j)) & 3321 *(field(i,k,j)+field(i-1,k,j)) 3322 ENDDO 3323 ENDIF 3324 3325 IF( i == ide-2 ) THEN ! third order flux one in from the boundary 3326 DO k=kts,ktf 3327 vel = ru(i,k,j) 3328 fqx( i,k ) = vel*flux3( field(i-2,k,j), field(i-1,k,j), & 3329 field(i ,k,j), field(i+1,k,j), & 3330 vel ) 3331 ENDDO 3332 ENDIF 3333 3334 ENDDO 3335 3336 ENDIF 3254 3337 3255 3338 ! x flux-divergence into tendency … … 3274 3357 IF( config_flags%periodic_x .or. & 3275 3358 config_flags%symmetric_xs .or. & 3276 (its > ids+ 1) ) degrade_xs = .false.3359 (its > ids+2) ) degrade_xs = .false. 3277 3360 IF( config_flags%periodic_x .or. & 3278 3361 config_flags%symmetric_xe .or. & … … 3280 3363 IF( config_flags%periodic_y .or. & 3281 3364 config_flags%symmetric_ys .or. & 3282 (jts > jds+ 1) ) degrade_ys = .false.3365 (jts > jds+2) ) degrade_ys = .false. 3283 3366 IF( config_flags%periodic_y .or. & 3284 3367 config_flags%symmetric_ye .or. & 3285 (jte < jde- 2) ) degrade_ye = .false.3368 (jte < jde-3) ) degrade_ye = .false. 3286 3369 3287 3370 ! begin flux computations … … 3376 3459 j_end_f = jde-2 3377 3460 ENDIF 3461 3378 3462 3379 3463 jp1 = 2 … … 3392 3476 DO k = kts, ktf 3393 3477 DO i = i_start, i_end 3394 fqy(i,k,jp1) = 0.5*rv(i,k,j_end+1) & 3395 *(field(i,k,j_end+1)+field(i,k,j_end)) 3478 ! Assumes j>j_end_f is ONLY j_end+1 ... 3479 ! fqy(i,k,jp1) = 0.5*rv(i,k,j_end+1) & 3480 ! *(field(i,k,j_end+1)+field(i,k,j_end)) 3481 fqy(i,k,jp1) = 0.5*rv(i,k,j) & 3482 *(field(i,k,j)+field(i,k,j-1)) 3396 3483 ENDDO 3397 3484 ENDDO … … 3407 3494 END IF 3408 3495 3496 ! y flux-divergence into tendency 3409 3497 IF ( j > j_start ) THEN 3410 ! y flux-divergence into tendency3411 3498 3412 3499 DO k=kts,ktf … … 3416 3503 ENDDO 3417 3504 ENDDO 3505 3506 3418 3507 END IF 3419 3508 … … 3434 3523 IF( config_flags%periodic_x .or. & 3435 3524 config_flags%symmetric_xs .or. & 3436 (its > ids+ 1) ) degrade_xs = .false.3525 (its > ids+2) ) degrade_xs = .false. 3437 3526 IF( config_flags%periodic_x .or. & 3438 3527 config_flags%symmetric_xe .or. & … … 3440 3529 IF( config_flags%periodic_y .or. & 3441 3530 config_flags%symmetric_ys .or. & 3442 (jts > jds+ 1) ) degrade_ys = .false.3531 (jts > jds+2) ) degrade_ys = .false. 3443 3532 IF( config_flags%periodic_y .or. & 3444 3533 config_flags%symmetric_ye .or. & 3445 (jte < jde- 2) ) degrade_ye = .false.3534 (jte < jde-3) ) degrade_ye = .false. 3446 3535 3447 3536 ! begin flux computations … … 3536 3625 j_end_f = jde-2 3537 3626 ENDIF 3627 3538 3628 3539 3629 jp1 = 2 … … 3552 3642 DO k = kts, ktf 3553 3643 DO i = i_start, i_end 3554 fqy(i,k,jp1) = 0.5*rv(i,k,j_end+1) & 3555 *(field(i,k,j_end+1)+field(i,k,j_end)) 3644 ! Assumes j>j_end_f is ONLY j_end+1 ... 3645 ! fqy(i,k,jp1) = 0.5*rv(i,k,j_end+1) & 3646 ! *(field(i,k,j_end+1)+field(i,k,j_end)) 3647 fqy(i,k,jp1) = 0.5*rv(i,k,j) & 3648 *(field(i,k,j)+field(i,k,j-1)) 3556 3649 ENDDO 3557 3650 ENDDO … … 3567 3660 END IF 3568 3661 3662 ! y flux-divergence into tendency 3569 3663 IF ( j > j_start ) THEN 3570 ! y flux-divergence into tendency3571 3664 3572 3665 DO k=kts,ktf … … 3576 3669 ENDDO 3577 3670 ENDDO 3671 3672 3578 3673 END IF 3579 3674 … … 3610 3705 i_end = MIN(ite,ide-1) 3611 3706 3707 ! 3612 3708 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts) 3613 3709 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte) … … 3624 3720 ENDDO 3625 3721 3722 3723 ELSE IF ( horz_order == 0 ) THEN 3724 3725 ! Just in case we want to turn horizontal advection off, we can do it 3726 3626 3727 ELSE 3627 3728 … … 3859 3960 ENDDO 3860 3961 3861 3862 3962 ELSE IF (vert_order == 2) THEN 3863 3963 … … 3890 3990 SUBROUTINE advect_w ( w, w_old, tendency, & 3891 3991 ru, rv, rom, & 3892 mut, config_flags,&3992 mut, time_step, config_flags, & 3893 3993 msfu, msfv, msft, & 3894 3994 fzm, fzp, & … … 3927 4027 REAL , INTENT(IN ) :: rdx, & 3928 4028 rdy 4029 INTEGER , INTENT(IN ) :: time_step 4030 3929 4031 3930 4032 ! Local data … … 3958 4060 flux3(q_im2, q_im1, q_i, q_ip1, ua) = & 3959 4061 flux4(q_im2, q_im1, q_i, q_ip1, ua) + & 3960 sign(1 .,ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.04062 sign(1,time_step)*sign(1.,ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0 3961 4063 3962 4064 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = & … … 3966 4068 flux5(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = & 3967 4069 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) & 3968 -sign(1 .,ua)*(&4070 -sign(1,time_step)*sign(1.,ua)*( & 3969 4071 (q_ip2-q_im3)-5.*(q_ip1-q_im2)+10.*(q_i-q_im1) )/60.0 3970 4072 … … 4001 4103 IF( config_flags%periodic_x .or. & 4002 4104 config_flags%symmetric_xs .or. & 4105 (its > ids+3) ) degrade_xs = .false. 4106 IF( config_flags%periodic_x .or. & 4107 config_flags%symmetric_xe .or. & 4108 (ite < ide-3) ) degrade_xe = .false. 4109 IF( config_flags%periodic_y .or. & 4110 config_flags%symmetric_ys .or. & 4111 (jts > jds+3) ) degrade_ys = .false. 4112 IF( config_flags%periodic_y .or. & 4113 config_flags%symmetric_ye .or. & 4114 (jte < jde-4) ) degrade_ye = .false. 4115 4116 !--------------- y - advection first 4117 4118 i_start = its 4119 i_end = MIN(ite,ide-1) 4120 j_start = jts 4121 j_end = MIN(jte,jde-1) 4122 4123 ! higher order flux has a 5 or 7 point stencil, so compute 4124 ! bounds so we can switch to second order flux close to the boundary 4125 4126 j_start_f = j_start 4127 j_end_f = j_end+1 4128 4129 IF(degrade_ys) then 4130 j_start = MAX(jts,jds+1) 4131 j_start_f = jds+3 4132 ENDIF 4133 4134 IF(degrade_ye) then 4135 j_end = MIN(jte,jde-2) 4136 j_end_f = jde-3 4137 ENDIF 4138 4139 4140 ! compute fluxes, 5th or 6th order 4141 4142 jp1 = 2 4143 jp0 = 1 4144 4145 j_loop_y_flux_6 : DO j = j_start, j_end+1 4146 4147 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN 4148 4149 DO k=kts+1,ktf 4150 DO i = i_start, i_end 4151 vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j) 4152 fqy( i, k, jp1 ) = vel*flux6( & 4153 w(i,k,j-3), w(i,k,j-2), w(i,k,j-1), & 4154 w(i,k,j ), w(i,k,j+1), w(i,k,j+2), vel ) 4155 ENDDO 4156 ENDDO 4157 4158 k = ktf+1 4159 DO i = i_start, i_end 4160 vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j) 4161 fqy( i, k, jp1 ) = vel*flux6( & 4162 w(i,k,j-3), w(i,k,j-2), w(i,k,j-1), & 4163 w(i,k,j ), w(i,k,j+1), w(i,k,j+2), vel ) 4164 ENDDO 4165 4166 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary 4167 4168 DO k=kts+1,ktf 4169 DO i = i_start, i_end 4170 fqy(i, k, jp1) = 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j))* & 4171 (w(i,k,j)+w(i,k,j-1)) 4172 ENDDO 4173 ENDDO 4174 4175 k = ktf+1 4176 DO i = i_start, i_end 4177 fqy(i, k, jp1) = 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j))* & 4178 (w(i,k,j)+w(i,k,j-1)) 4179 ENDDO 4180 4181 ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary 4182 4183 DO k=kts+1,ktf 4184 DO i = i_start, i_end 4185 vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j) 4186 fqy( i, k, jp1 ) = vel*flux4( & 4187 w(i,k,j-2),w(i,k,j-1),w(i,k,j),w(i,k,j+1),vel ) 4188 ENDDO 4189 ENDDO 4190 4191 k = ktf+1 4192 DO i = i_start, i_end 4193 vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j) 4194 fqy( i, k, jp1 ) = vel*flux4( & 4195 w(i,k,j-2),w(i,k,j-1),w(i,k,j),w(i,k,j+1),vel ) 4196 ENDDO 4197 4198 ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary 4199 4200 DO k=kts+1,ktf 4201 DO i = i_start, i_end 4202 fqy(i, k, jp1) = 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j))* & 4203 (w(i,k,j)+w(i,k,j-1)) 4204 ENDDO 4205 ENDDO 4206 4207 k = ktf+1 4208 DO i = i_start, i_end 4209 fqy(i, k, jp1) = 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j))* & 4210 (w(i,k,j)+w(i,k,j-1)) 4211 ENDDO 4212 4213 ELSE IF ( j == jde-2 ) THEN ! 3rd or 4th order flux 2 in from north boundary 4214 4215 DO k=kts+1,ktf 4216 DO i = i_start, i_end 4217 vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j) 4218 fqy( i, k, jp1 ) = vel*flux4( & 4219 w(i,k,j-2),w(i,k,j-1), & 4220 w(i,k,j),w(i,k,j+1),vel ) 4221 ENDDO 4222 ENDDO 4223 4224 k = ktf+1 4225 DO i = i_start, i_end 4226 vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j) 4227 fqy( i, k, jp1 ) = vel*flux4( & 4228 w(i,k,j-2),w(i,k,j-1), & 4229 w(i,k,j),w(i,k,j+1),vel ) 4230 ENDDO 4231 4232 ENDIF 4233 4234 ! y flux-divergence into tendency 4235 4236 IF(j > j_start) THEN 4237 4238 DO k=kts+1,ktf+1 4239 DO i = i_start, i_end 4240 mrdy=msft(i,j-1)*rdy 4241 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0)) 4242 ENDDO 4243 ENDDO 4244 4245 ENDIF 4246 4247 4248 jtmp = jp1 4249 jp1 = jp0 4250 jp0 = jtmp 4251 4252 ENDDO j_loop_y_flux_6 4253 4254 ! next, x - flux divergence 4255 4256 i_start = its 4257 i_end = MIN(ite,ide-1) 4258 4259 j_start = jts 4260 j_end = MIN(jte,jde-1) 4261 4262 ! higher order flux has a 5 or 7 point stencil, so compute 4263 ! bounds so we can switch to second order flux close to the boundary 4264 4265 i_start_f = i_start 4266 i_end_f = i_end+1 4267 4268 IF(degrade_xs) then 4269 i_start = MAX(ids+1,its) 4270 ! i_start_f = i_start+2 4271 i_start_f = MIN(i_start+2,ids+3) 4272 ENDIF 4273 4274 IF(degrade_xe) then 4275 i_end = MIN(ide-2,ite) 4276 i_end_f = ide-3 4277 ENDIF 4278 4279 ! compute fluxes 4280 4281 DO j = j_start, j_end 4282 4283 ! 5th or 6th order flux 4284 4285 DO k=kts+1,ktf 4286 DO i = i_start_f, i_end_f 4287 vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j) 4288 fqx( i,k ) = vel*flux6( w(i-3,k,j), w(i-2,k,j), & 4289 w(i-1,k,j), w(i ,k,j), & 4290 w(i+1,k,j), w(i+2,k,j), & 4291 vel ) 4292 ENDDO 4293 ENDDO 4294 4295 k = ktf+1 4296 DO i = i_start_f, i_end_f 4297 vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j) 4298 fqx( i,k ) = vel*flux6( w(i-3,k,j), w(i-2,k,j), & 4299 w(i-1,k,j), w(i ,k,j), & 4300 w(i+1,k,j), w(i+2,k,j), & 4301 vel ) 4302 ENDDO 4303 4304 ! lower order fluxes close to boundaries (if not periodic or symmetric) 4305 4306 IF( degrade_xs ) THEN 4307 4308 DO i=i_start,i_start_f-1 4309 4310 IF(i == ids+1) THEN ! second order 4311 DO k=kts+1,ktf 4312 fqx(i,k) = 0.5*(fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)) & 4313 *(w(i,k,j)+w(i-1,k,j)) 4314 ENDDO 4315 k = ktf+1 4316 fqx(i,k) = 0.5*((2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)) & 4317 *(w(i,k,j)+w(i-1,k,j)) 4318 ENDIF 4319 4320 IF(i == ids+2) THEN ! third order 4321 DO k=kts+1,ktf 4322 vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j) 4323 fqx( i,k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j), & 4324 w(i ,k,j), w(i+1,k,j), & 4325 vel ) 4326 ENDDO 4327 k = ktf+1 4328 vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j) 4329 fqx( i,k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j), & 4330 w(i ,k,j), w(i+1,k,j), & 4331 vel ) 4332 END IF 4333 4334 ENDDO 4335 4336 ENDIF 4337 4338 IF( degrade_xe ) THEN 4339 4340 DO i = i_end_f+1, i_end+1 4341 4342 IF( i == ide-1 ) THEN ! second order flux next to the boundary 4343 DO k=kts+1,ktf 4344 fqx(i,k) = 0.5*(fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)) & 4345 *(w(i,k,j)+w(i-1,k,j)) 4346 ENDDO 4347 k = ktf+1 4348 fqx(i,k) = 0.5*((2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)) & 4349 *(w(i,k,j)+w(i-1,k,j)) 4350 ENDIF 4351 4352 IF( i == ide-2 ) THEN ! third order flux one in from the boundary 4353 DO k=kts+1,ktf 4354 vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j) 4355 fqx( i,k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j), & 4356 w(i ,k,j), w(i+1,k,j), & 4357 vel ) 4358 ENDDO 4359 k = ktf+1 4360 vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j) 4361 fqx( i,k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j), & 4362 w(i ,k,j), w(i+1,k,j), & 4363 vel ) 4364 ENDIF 4365 4366 ENDDO 4367 4368 ENDIF 4369 4370 ! x flux-divergence into tendency 4371 4372 DO k=kts+1,ktf+1 4373 DO i = i_start, i_end 4374 mrdx=msft(i,j)*rdx 4375 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k)) 4376 ENDDO 4377 ENDDO 4378 4379 ENDDO 4380 4381 ELSE IF (horz_order == 5 ) THEN 4382 4383 ! determine boundary mods for flux operators 4384 ! We degrade the flux operators from 3rd/4th order 4385 ! to second order one gridpoint in from the boundaries for 4386 ! all boundary conditions except periodic and symmetry - these 4387 ! conditions have boundary zone data fill for correct application 4388 ! of the higher order flux stencils 4389 4390 degrade_xs = .true. 4391 degrade_xe = .true. 4392 degrade_ys = .true. 4393 degrade_ye = .true. 4394 4395 IF( config_flags%periodic_x .or. & 4396 config_flags%symmetric_xs .or. & 4397 (its > ids+3) ) degrade_xs = .false. 4398 IF( config_flags%periodic_x .or. & 4399 config_flags%symmetric_xe .or. & 4400 (ite < ide-3) ) degrade_xe = .false. 4401 IF( config_flags%periodic_y .or. & 4402 config_flags%symmetric_ys .or. & 4403 (jts > jds+3) ) degrade_ys = .false. 4404 IF( config_flags%periodic_y .or. & 4405 config_flags%symmetric_ye .or. & 4406 (jte < jde-4) ) degrade_ye = .false. 4407 4408 !--------------- y - advection first 4409 4410 i_start = its 4411 i_end = MIN(ite,ide-1) 4412 j_start = jts 4413 j_end = MIN(jte,jde-1) 4414 4415 ! higher order flux has a 5 or 7 point stencil, so compute 4416 ! bounds so we can switch to second order flux close to the boundary 4417 4418 j_start_f = j_start 4419 j_end_f = j_end+1 4420 4421 IF(degrade_ys) then 4422 j_start = MAX(jts,jds+1) 4423 j_start_f = jds+3 4424 ENDIF 4425 4426 IF(degrade_ye) then 4427 j_end = MIN(jte,jde-2) 4428 j_end_f = jde-3 4429 ENDIF 4430 4431 4432 ! compute fluxes, 5th or 6th order 4433 4434 jp1 = 2 4435 jp0 = 1 4436 4437 j_loop_y_flux_5 : DO j = j_start, j_end+1 4438 4439 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN 4440 4441 DO k=kts+1,ktf 4442 DO i = i_start, i_end 4443 vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j) 4444 fqy( i, k, jp1 ) = vel*flux5( & 4445 w(i,k,j-3), w(i,k,j-2), w(i,k,j-1), & 4446 w(i,k,j ), w(i,k,j+1), w(i,k,j+2), vel ) 4447 ENDDO 4448 ENDDO 4449 4450 k = ktf+1 4451 DO i = i_start, i_end 4452 vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j) 4453 fqy( i, k, jp1 ) = vel*flux5( & 4454 w(i,k,j-3), w(i,k,j-2), w(i,k,j-1), & 4455 w(i,k,j ), w(i,k,j+1), w(i,k,j+2), vel ) 4456 ENDDO 4457 4458 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary 4459 4460 DO k=kts+1,ktf 4461 DO i = i_start, i_end 4462 fqy(i, k, jp1) = 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j))* & 4463 (w(i,k,j)+w(i,k,j-1)) 4464 ENDDO 4465 ENDDO 4466 4467 k = ktf+1 4468 DO i = i_start, i_end 4469 fqy(i, k, jp1) = 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j))* & 4470 (w(i,k,j)+w(i,k,j-1)) 4471 ENDDO 4472 4473 ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary 4474 4475 DO k=kts+1,ktf 4476 DO i = i_start, i_end 4477 vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j) 4478 fqy( i, k, jp1 ) = vel*flux3( & 4479 w(i,k,j-2),w(i,k,j-1),w(i,k,j),w(i,k,j+1),vel ) 4480 ENDDO 4481 ENDDO 4482 4483 k = ktf+1 4484 DO i = i_start, i_end 4485 vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j) 4486 fqy( i, k, jp1 ) = vel*flux3( & 4487 w(i,k,j-2),w(i,k,j-1),w(i,k,j),w(i,k,j+1),vel ) 4488 ENDDO 4489 4490 ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary 4491 4492 DO k=kts+1,ktf 4493 DO i = i_start, i_end 4494 fqy(i, k, jp1) = 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j))* & 4495 (w(i,k,j)+w(i,k,j-1)) 4496 ENDDO 4497 ENDDO 4498 4499 k = ktf+1 4500 DO i = i_start, i_end 4501 fqy(i, k, jp1) = 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j))* & 4502 (w(i,k,j)+w(i,k,j-1)) 4503 ENDDO 4504 4505 ELSE IF ( j == jde-2 ) THEN ! 3rd or 4th order flux 2 in from north boundary 4506 4507 DO k=kts+1,ktf 4508 DO i = i_start, i_end 4509 vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j) 4510 fqy( i, k, jp1 ) = vel*flux3( & 4511 w(i,k,j-2),w(i,k,j-1), & 4512 w(i,k,j),w(i,k,j+1),vel ) 4513 ENDDO 4514 ENDDO 4515 4516 k = ktf+1 4517 DO i = i_start, i_end 4518 vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j) 4519 fqy( i, k, jp1 ) = vel*flux3( & 4520 w(i,k,j-2),w(i,k,j-1), & 4521 w(i,k,j),w(i,k,j+1),vel ) 4522 ENDDO 4523 4524 ENDIF 4525 4526 ! y flux-divergence into tendency 4527 4528 IF(j > j_start) THEN 4529 4530 DO k=kts+1,ktf+1 4531 DO i = i_start, i_end 4532 mrdy=msft(i,j-1)*rdy 4533 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0)) 4534 ENDDO 4535 ENDDO 4536 4537 ENDIF 4538 4539 4540 jtmp = jp1 4541 jp1 = jp0 4542 jp0 = jtmp 4543 4544 ENDDO j_loop_y_flux_5 4545 4546 ! next, x - flux divergence 4547 4548 i_start = its 4549 i_end = MIN(ite,ide-1) 4550 4551 j_start = jts 4552 j_end = MIN(jte,jde-1) 4553 4554 ! higher order flux has a 5 or 7 point stencil, so compute 4555 ! bounds so we can switch to second order flux close to the boundary 4556 4557 i_start_f = i_start 4558 i_end_f = i_end+1 4559 4560 IF(degrade_xs) then 4561 i_start = MAX(ids+1,its) 4562 ! i_start_f = i_start+2 4563 i_start_f = MIN(i_start+2,ids+3) 4564 ENDIF 4565 4566 IF(degrade_xe) then 4567 i_end = MIN(ide-2,ite) 4568 i_end_f = ide-3 4569 ENDIF 4570 4571 ! compute fluxes 4572 4573 DO j = j_start, j_end 4574 4575 ! 5th or 6th order flux 4576 4577 DO k=kts+1,ktf 4578 DO i = i_start_f, i_end_f 4579 vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j) 4580 fqx( i,k ) = vel*flux5( w(i-3,k,j), w(i-2,k,j), & 4581 w(i-1,k,j), w(i ,k,j), & 4582 w(i+1,k,j), w(i+2,k,j), & 4583 vel ) 4584 ENDDO 4585 ENDDO 4586 4587 k = ktf+1 4588 DO i = i_start_f, i_end_f 4589 vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j) 4590 fqx( i,k ) = vel*flux5( w(i-3,k,j), w(i-2,k,j), & 4591 w(i-1,k,j), w(i ,k,j), & 4592 w(i+1,k,j), w(i+2,k,j), & 4593 vel ) 4594 ENDDO 4595 4596 ! lower order fluxes close to boundaries (if not periodic or symmetric) 4597 4598 IF( degrade_xs ) THEN 4599 4600 DO i=i_start,i_start_f-1 4601 4602 IF(i == ids+1) THEN ! second order 4603 DO k=kts+1,ktf 4604 fqx(i,k) = 0.5*(fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)) & 4605 *(w(i,k,j)+w(i-1,k,j)) 4606 ENDDO 4607 k = ktf+1 4608 fqx(i,k) = 0.5*((2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)) & 4609 *(w(i,k,j)+w(i-1,k,j)) 4610 ENDIF 4611 4612 IF(i == ids+2) THEN ! third order 4613 DO k=kts+1,ktf 4614 vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j) 4615 fqx( i,k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j), & 4616 w(i ,k,j), w(i+1,k,j), & 4617 vel ) 4618 ENDDO 4619 k = ktf+1 4620 vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j) 4621 fqx( i,k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j), & 4622 w(i ,k,j), w(i+1,k,j), & 4623 vel ) 4624 END IF 4625 4626 ENDDO 4627 4628 ENDIF 4629 4630 IF( degrade_xe ) THEN 4631 4632 DO i = i_end_f+1, i_end+1 4633 4634 IF( i == ide-1 ) THEN ! second order flux next to the boundary 4635 DO k=kts+1,ktf 4636 fqx(i,k) = 0.5*(fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)) & 4637 *(w(i,k,j)+w(i-1,k,j)) 4638 ENDDO 4639 k = ktf+1 4640 fqx(i,k) = 0.5*((2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)) & 4641 *(w(i,k,j)+w(i-1,k,j)) 4642 ENDIF 4643 4644 IF( i == ide-2 ) THEN ! third order flux one in from the boundary 4645 DO k=kts+1,ktf 4646 vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j) 4647 fqx( i,k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j), & 4648 w(i ,k,j), w(i+1,k,j), & 4649 vel ) 4650 ENDDO 4651 k = ktf+1 4652 vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j) 4653 fqx( i,k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j), & 4654 w(i ,k,j), w(i+1,k,j), & 4655 vel ) 4656 ENDIF 4657 4658 ENDDO 4659 4660 ENDIF 4661 4662 ! x flux-divergence into tendency 4663 4664 DO k=kts+1,ktf+1 4665 DO i = i_start, i_end 4666 mrdx=msft(i,j)*rdx 4667 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k)) 4668 ENDDO 4669 ENDDO 4670 4671 ENDDO 4672 4673 ELSE IF ( horz_order == 4 ) THEN 4674 4675 degrade_xs = .true. 4676 degrade_xe = .true. 4677 degrade_ys = .true. 4678 degrade_ye = .true. 4679 4680 IF( config_flags%periodic_x .or. & 4681 config_flags%symmetric_xs .or. & 4003 4682 (its > ids+2) ) degrade_xs = .false. 4004 4683 IF( config_flags%periodic_x .or. & 4005 4684 config_flags%symmetric_xe .or. & 4006 (ite < ide- 3) ) degrade_xe = .false.4685 (ite < ide-2) ) degrade_xe = .false. 4007 4686 IF( config_flags%periodic_y .or. & 4008 4687 config_flags%symmetric_ys .or. & … … 4012 4691 (jte < jde-3) ) degrade_ye = .false. 4013 4692 4014 !--------------- y - advection first4015 4016 i_start = its4017 i_end = MIN(ite,ide-1)4018 j_start = jts4019 j_end = MIN(jte,jde-1)4020 4021 ! higher order flux has a 5 or 7 point stencil, so compute4022 ! bounds so we can switch to second order flux close to the boundary4023 4024 j_start_f = j_start4025 j_end_f = j_end+14026 4027 IF(degrade_ys) then4028 j_start = MAX(jts,jds+1)4029 j_start_f = jds+34030 ENDIF4031 4032 IF(degrade_ye) then4033 j_end = MIN(jte,jde-2)4034 j_end_f = jde-34035 ENDIF4036 4037 ! compute fluxes, 5th or 6th order4038 4039 jp1 = 24040 jp0 = 14041 4042 j_loop_y_flux_6 : DO j = j_start, j_end+14043 4044 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN4045 4046 DO k=kts+1,ktf4047 DO i = i_start, i_end4048 vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)4049 fqy( i, k, jp1 ) = vel*flux6( &4050 w(i,k,j-3), w(i,k,j-2), w(i,k,j-1), &4051 w(i,k,j ), w(i,k,j+1), w(i,k,j+2), vel )4052 ENDDO4053 ENDDO4054 4055 k = ktf+14056 DO i = i_start, i_end4057 vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)4058 fqy( i, k, jp1 ) = vel*flux6( &4059 w(i,k,j-3), w(i,k,j-2), w(i,k,j-1), &4060 w(i,k,j ), w(i,k,j+1), w(i,k,j+2), vel )4061 ENDDO4062 4063 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary4064 4065 DO k=kts+1,ktf4066 DO i = i_start, i_end4067 fqy(i, k, jp1) = 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j))* &4068 (w(i,k,j)+w(i,k,j-1))4069 ENDDO4070 ENDDO4071 4072 k = ktf+14073 DO i = i_start, i_end4074 fqy(i, k, jp1) = 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j))* &4075 (w(i,k,j)+w(i,k,j-1))4076 ENDDO4077 4078 ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary4079 4080 DO k=kts+1,ktf4081 DO i = i_start, i_end4082 vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)4083 fqy( i, k, jp1 ) = vel*flux4( &4084 w(i,k,j-2),w(i,k,j-1),w(i,k,j),w(i,k,j+1),vel )4085 ENDDO4086 ENDDO4087 4088 k = ktf+14089 DO i = i_start, i_end4090 vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)4091 fqy( i, k, jp1 ) = vel*flux4( &4092 w(i,k,j-2),w(i,k,j-1),w(i,k,j),w(i,k,j+1),vel )4093 ENDDO4094 4095 ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary4096 4097 DO k=kts+1,ktf4098 DO i = i_start, i_end4099 fqy(i, k, jp1) = 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j))* &4100 (w(i,k,j)+w(i,k,j-1))4101 ENDDO4102 ENDDO4103 4104 k = ktf+14105 DO i = i_start, i_end4106 fqy(i, k, jp1) = 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j))* &4107 (w(i,k,j)+w(i,k,j-1))4108 ENDDO4109 4110 ELSE IF ( j == jde-2 ) THEN ! 3rd or 4th order flux 2 in from north boundary4111 4112 DO k=kts+1,ktf4113 DO i = i_start, i_end4114 vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)4115 fqy( i, k, jp1 ) = vel*flux4( &4116 w(i,k,j-2),w(i,k,j-1), &4117 w(i,k,j),w(i,k,j+1),vel )4118 ENDDO4119 ENDDO4120 4121 k = ktf+14122 DO i = i_start, i_end4123 vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)4124 fqy( i, k, jp1 ) = vel*flux4( &4125 w(i,k,j-2),w(i,k,j-1), &4126 w(i,k,j),w(i,k,j+1),vel )4127 ENDDO4128 4129 ENDIF4130 4131 ! y flux-divergence into tendency4132 4133 IF(j > j_start) THEN4134 4135 DO k=kts+1,ktf+14136 DO i = i_start, i_end4137 mrdy=msft(i,j-1)*rdy4138 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))4139 ENDDO4140 ENDDO4141 4142 ENDIF4143 4144 jtmp = jp14145 jp1 = jp04146 jp0 = jtmp4147 4148 ENDDO j_loop_y_flux_64149 4150 ! next, x - flux divergence4151 4152 i_start = its4153 i_end = MIN(ite,ide-1)4154 4155 j_start = jts4156 j_end = MIN(jte,jde-1)4157 4158 ! higher order flux has a 5 or 7 point stencil, so compute4159 ! bounds so we can switch to second order flux close to the boundary4160 4161 i_start_f = i_start4162 i_end_f = i_end+14163 4164 IF(degrade_xs) then4165 i_start = MAX(ids+1,its)4166 i_start_f = i_start+24167 ENDIF4168 4169 IF(degrade_xe) then4170 i_end = MIN(ide-2,ite)4171 i_end_f = ide-34172 ENDIF4173 4174 ! compute fluxes4175 4176 DO j = j_start, j_end4177 4178 ! 5th or 6th order flux4179 4180 DO k=kts+1,ktf4181 DO i = i_start_f, i_end_f4182 vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)4183 fqx( i,k ) = vel*flux6( w(i-3,k,j), w(i-2,k,j), &4184 w(i-1,k,j), w(i ,k,j), &4185 w(i+1,k,j), w(i+2,k,j), &4186 vel )4187 ENDDO4188 ENDDO4189 4190 k = ktf+14191 DO i = i_start_f, i_end_f4192 vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)4193 fqx( i,k ) = vel*flux6( w(i-3,k,j), w(i-2,k,j), &4194 w(i-1,k,j), w(i ,k,j), &4195 w(i+1,k,j), w(i+2,k,j), &4196 vel )4197 ENDDO4198 4199 ! lower order fluxes close to boundaries (if not periodic or symmetric)4200 4201 IF( degrade_xs ) THEN4202 4203 IF( i_start == ids+1 ) THEN ! second order flux next to the boundary4204 i = ids+14205 DO k=kts+1,ktf4206 fqx(i,k) = 0.5*(fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)) &4207 *(w(i,k,j)+w(i-1,k,j))4208 ENDDO4209 k = ktf+14210 fqx(i,k) = 0.5*((2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)) &4211 *(w(i,k,j)+w(i-1,k,j))4212 ENDIF4213 4214 DO k=kts+1,ktf4215 i = i_start+14216 vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)4217 fqx( i,k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j), &4218 w(i ,k,j), w(i+1,k,j), &4219 vel )4220 ENDDO4221 4222 k = ktf+14223 vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)4224 fqx( i,k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j), &4225 w(i ,k,j), w(i+1,k,j), &4226 vel )4227 ENDIF4228 4229 IF( degrade_xe ) THEN4230 4231 IF( i_end == ide-2 ) THEN ! second order flux next to the boundary4232 i = ide-14233 DO k=kts+1,ktf4234 fqx(i,k) = 0.5*(fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)) &4235 *(w(i,k,j)+w(i-1,k,j))4236 ENDDO4237 k = ktf+14238 fqx(i,k) = 0.5*((2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)) &4239 *(w(i,k,j)+w(i-1,k,j))4240 ENDIF4241 4242 i = ide-24243 DO k=kts+1,ktf4244 vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)4245 fqx( i,k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j), &4246 w(i ,k,j), w(i+1,k,j), &4247 vel )4248 ENDDO4249 4250 k = ktf+14251 vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)4252 fqx( i,k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j), &4253 w(i ,k,j), w(i+1,k,j), &4254 vel )4255 ENDIF4256 4257 ! x flux-divergence into tendency4258 4259 DO k=kts+1,ktf+14260 DO i = i_start, i_end4261 mrdx=msft(i,j)*rdx4262 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))4263 ENDDO4264 ENDDO4265 4266 ENDDO4267 4268 4269 ELSE IF (horz_order == 5 ) THEN4270 4271 ! determine boundary mods for flux operators4272 ! We degrade the flux operators from 3rd/4th order4273 ! to second order one gridpoint in from the boundaries for4274 ! all boundary conditions except periodic and symmetry - these4275 ! conditions have boundary zone data fill for correct application4276 ! of the higher order flux stencils4277 4278 degrade_xs = .true.4279 degrade_xe = .true.4280 degrade_ys = .true.4281 degrade_ye = .true.4282 4283 IF( config_flags%periodic_x .or. &4284 config_flags%symmetric_xs .or. &4285 (its > ids+2) ) degrade_xs = .false.4286 IF( config_flags%periodic_x .or. &4287 config_flags%symmetric_xe .or. &4288 (ite < ide-3) ) degrade_xe = .false.4289 IF( config_flags%periodic_y .or. &4290 config_flags%symmetric_ys .or. &4291 (jts > jds+2) ) degrade_ys = .false.4292 IF( config_flags%periodic_y .or. &4293 config_flags%symmetric_ye .or. &4294 (jte < jde-3) ) degrade_ye = .false.4295 4296 !--------------- y - advection first4297 4298 i_start = its4299 i_end = MIN(ite,ide-1)4300 j_start = jts4301 j_end = MIN(jte,jde-1)4302 4303 ! higher order flux has a 5 or 7 point stencil, so compute4304 ! bounds so we can switch to second order flux close to the boundary4305 4306 j_start_f = j_start4307 j_end_f = j_end+14308 4309 IF(degrade_ys) then4310 j_start = MAX(jts,jds+1)4311 j_start_f = jds+34312 ENDIF4313 4314 IF(degrade_ye) then4315 j_end = MIN(jte,jde-2)4316 j_end_f = jde-34317 ENDIF4318 4319 ! compute fluxes, 5th or 6th order4320 4321 jp1 = 24322 jp0 = 14323 4324 j_loop_y_flux_5 : DO j = j_start, j_end+14325 4326 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN4327 4328 DO k=kts+1,ktf4329 DO i = i_start, i_end4330 vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)4331 fqy( i, k, jp1 ) = vel*flux5( &4332 w(i,k,j-3), w(i,k,j-2), w(i,k,j-1), &4333 w(i,k,j ), w(i,k,j+1), w(i,k,j+2), vel )4334 ENDDO4335 ENDDO4336 4337 k = ktf+14338 DO i = i_start, i_end4339 vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)4340 fqy( i, k, jp1 ) = vel*flux5( &4341 w(i,k,j-3), w(i,k,j-2), w(i,k,j-1), &4342 w(i,k,j ), w(i,k,j+1), w(i,k,j+2), vel )4343 ENDDO4344 4345 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary4346 4347 DO k=kts+1,ktf4348 DO i = i_start, i_end4349 fqy(i, k, jp1) = 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j))* &4350 (w(i,k,j)+w(i,k,j-1))4351 ENDDO4352 ENDDO4353 4354 k = ktf+14355 DO i = i_start, i_end4356 fqy(i, k, jp1) = 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j))* &4357 (w(i,k,j)+w(i,k,j-1))4358 ENDDO4359 4360 ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary4361 4362 DO k=kts+1,ktf4363 DO i = i_start, i_end4364 vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)4365 fqy( i, k, jp1 ) = vel*flux3( &4366 w(i,k,j-2),w(i,k,j-1),w(i,k,j),w(i,k,j+1),vel )4367 ENDDO4368 ENDDO4369 4370 k = ktf+14371 DO i = i_start, i_end4372 vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)4373 fqy( i, k, jp1 ) = vel*flux3( &4374 w(i,k,j-2),w(i,k,j-1),w(i,k,j),w(i,k,j+1),vel )4375 ENDDO4376 4377 ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary4378 4379 DO k=kts+1,ktf4380 DO i = i_start, i_end4381 fqy(i, k, jp1) = 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j))* &4382 (w(i,k,j)+w(i,k,j-1))4383 ENDDO4384 ENDDO4385 4386 k = ktf+14387 DO i = i_start, i_end4388 fqy(i, k, jp1) = 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j))* &4389 (w(i,k,j)+w(i,k,j-1))4390 ENDDO4391 4392 ELSE IF ( j == jde-2 ) THEN ! 3rd or 4th order flux 2 in from north boundary4393 4394 DO k=kts+1,ktf4395 DO i = i_start, i_end4396 vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)4397 fqy( i, k, jp1 ) = vel*flux3( &4398 w(i,k,j-2),w(i,k,j-1), &4399 w(i,k,j),w(i,k,j+1),vel )4400 ENDDO4401 ENDDO4402 4403 k = ktf+14404 DO i = i_start, i_end4405 vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)4406 fqy( i, k, jp1 ) = vel*flux3( &4407 w(i,k,j-2),w(i,k,j-1), &4408 w(i,k,j),w(i,k,j+1),vel )4409 ENDDO4410 4411 ENDIF4412 4413 ! y flux-divergence into tendency4414 4415 IF(j > j_start) THEN4416 4417 DO k=kts+1,ktf+14418 DO i = i_start, i_end4419 mrdy=msft(i,j-1)*rdy4420 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))4421 ENDDO4422 ENDDO4423 4424 ENDIF4425 4426 jtmp = jp14427 jp1 = jp04428 jp0 = jtmp4429 4430 ENDDO j_loop_y_flux_54431 4432 ! next, x - flux divergence4433 4434 i_start = its4435 i_end = MIN(ite,ide-1)4436 4437 j_start = jts4438 j_end = MIN(jte,jde-1)4439 4440 ! higher order flux has a 5 or 7 point stencil, so compute4441 ! bounds so we can switch to second order flux close to the boundary4442 4443 i_start_f = i_start4444 i_end_f = i_end+14445 4446 IF(degrade_xs) then4447 i_start = MAX(ids+1,its)4448 i_start_f = i_start+24449 ENDIF4450 4451 IF(degrade_xe) then4452 i_end = MIN(ide-2,ite)4453 i_end_f = ide-34454 ENDIF4455 4456 ! compute fluxes4457 4458 DO j = j_start, j_end4459 4460 ! 5th or 6th order flux4461 4462 DO k=kts+1,ktf4463 DO i = i_start_f, i_end_f4464 vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)4465 fqx( i,k ) = vel*flux5( w(i-3,k,j), w(i-2,k,j), &4466 w(i-1,k,j), w(i ,k,j), &4467 w(i+1,k,j), w(i+2,k,j), &4468 vel )4469 ENDDO4470 ENDDO4471 4472 k = ktf+14473 DO i = i_start_f, i_end_f4474 vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)4475 fqx( i,k ) = vel*flux5( w(i-3,k,j), w(i-2,k,j), &4476 w(i-1,k,j), w(i ,k,j), &4477 w(i+1,k,j), w(i+2,k,j), &4478 vel )4479 ENDDO4480 4481 ! lower order fluxes close to boundaries (if not periodic or symmetric)4482 4483 IF( degrade_xs ) THEN4484 4485 IF( i_start == ids+1 ) THEN ! second order flux next to the boundary4486 i = ids+14487 DO k=kts+1,ktf4488 fqx(i,k) = 0.5*(fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)) &4489 *(w(i,k,j)+w(i-1,k,j))4490 ENDDO4491 k = ktf+14492 fqx(i,k) = 0.5*((2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)) &4493 *(w(i,k,j)+w(i-1,k,j))4494 ENDIF4495 4496 i = i_start+14497 DO k=kts+1,ktf4498 vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)4499 fqx( i,k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j), &4500 w(i ,k,j), w(i+1,k,j), &4501 vel )4502 ENDDO4503 k = ktf+14504 vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)4505 fqx( i,k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j), &4506 w(i ,k,j), w(i+1,k,j), &4507 vel )4508 4509 ENDIF4510 4511 IF( degrade_xe ) THEN4512 4513 IF( i_end == ide-2 ) THEN ! second order flux next to the boundary4514 i = ide-14515 DO k=kts+1,ktf4516 fqx(i,k) = 0.5*(fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)) &4517 *(w(i,k,j)+w(i-1,k,j))4518 ENDDO4519 k = ktf+14520 fqx(i,k) = 0.5*((2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)) &4521 *(w(i,k,j)+w(i-1,k,j))4522 ENDIF4523 4524 i = ide-24525 DO k=kts+1,ktf4526 vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)4527 fqx( i,k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j), &4528 w(i ,k,j), w(i+1,k,j), &4529 vel )4530 ENDDO4531 k = ktf+14532 vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)4533 fqx( i,k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j), &4534 w(i ,k,j), w(i+1,k,j), &4535 vel )4536 ENDIF4537 4538 ! x flux-divergence into tendency4539 4540 DO k=kts+1,ktf+14541 DO i = i_start, i_end4542 mrdx=msft(i,j)*rdx4543 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))4544 ENDDO4545 ENDDO4546 4547 ENDDO4548 4549 ELSE IF ( horz_order == 4 ) THEN4550 4551 degrade_xs = .true.4552 degrade_xe = .true.4553 degrade_ys = .true.4554 degrade_ye = .true.4555 4556 IF( config_flags%periodic_x .or. &4557 config_flags%symmetric_xs .or. &4558 (its > ids+1) ) degrade_xs = .false.4559 IF( config_flags%periodic_x .or. &4560 config_flags%symmetric_xe .or. &4561 (ite < ide-2) ) degrade_xe = .false.4562 IF( config_flags%periodic_y .or. &4563 config_flags%symmetric_ys .or. &4564 (jts > jds+1) ) degrade_ys = .false.4565 IF( config_flags%periodic_y .or. &4566 config_flags%symmetric_ye .or. &4567 (jte < jde-2) ) degrade_ye = .false.4568 4569 4693 ! begin flux computations 4570 4694 ! start with x flux divergence … … 4676 4800 ENDIF 4677 4801 4802 4678 4803 jp1 = 2 4679 4804 jp0 = 1 … … 4698 4823 DO k = kts+1, ktf 4699 4824 DO i = i_start, i_end 4825 ! Assumes j>j_end_f is ONLY j_end+1 ... 4826 ! fqy(i, k, jp1) = & 4827 ! 0.5*(fzm(k)*rv(i,k,j_end+1)+fzp(k)*rv(i,k-1,j_end+1)) & 4828 ! *(w(i,k,j_end+1)+w(i,k,j_end)) 4700 4829 fqy(i, k, jp1) = & 4701 0.5*(fzm(k)*rv(i,k,j _end+1)+fzp(k)*rv(i,k-1,j_end+1)) &4702 *(w(i,k,j _end+1)+w(i,k,j_end))4830 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)) & 4831 *(w(i,k,j)+w(i,k,j-1)) 4703 4832 ENDDO 4704 4833 ENDDO 4705 4834 k = ktf+1 4706 4835 DO i = i_start, i_end 4836 ! Assumes j>j_end_f is ONLY j_end+1 ... 4837 ! fqy(i, k, jp1) = & 4838 ! 0.5*((2.-fzm(k-1))*rv(i,k-1,j_end+1)-fzp(k-1)*rv(i,k-2,j_end+1)) & 4839 ! *(w(i,k,j_end+1)+w(i,k,j_end)) 4707 4840 fqy(i, k, jp1) = & 4708 0.5*((2.-fzm(k-1))*rv(i,k-1,j _end+1)-fzp(k-1)*rv(i,k-2,j_end+1)) &4709 *(w(i,k,j _end+1)+w(i,k,j_end))4841 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)) & 4842 *(w(i,k,j)+w(i,k,j-1)) 4710 4843 ENDDO 4711 4844 ELSE … … 4728 4861 END IF 4729 4862 4863 ! y flux-divergence into tendency 4730 4864 IF( j > j_start ) THEN 4731 ! y flux-divergence into tendency 4865 4732 4866 DO k = kts+1, ktf+1 4733 4867 DO i = i_start, i_end … … 4736 4870 ENDDO 4737 4871 ENDDO 4872 4873 4738 4874 END IF 4739 4875 … … 4753 4889 IF( config_flags%periodic_x .or. & 4754 4890 config_flags%symmetric_xs .or. & 4755 (its > ids+ 1) ) degrade_xs = .false.4891 (its > ids+2) ) degrade_xs = .false. 4756 4892 IF( config_flags%periodic_x .or. & 4757 4893 config_flags%symmetric_xe .or. & … … 4759 4895 IF( config_flags%periodic_y .or. & 4760 4896 config_flags%symmetric_ys .or. & 4761 (jts > jds+ 1) ) degrade_ys = .false.4897 (jts > jds+2) ) degrade_ys = .false. 4762 4898 IF( config_flags%periodic_y .or. & 4763 4899 config_flags%symmetric_ye .or. & 4764 (jte < jde- 2) ) degrade_ye = .false.4900 (jte < jde-3) ) degrade_ye = .false. 4765 4901 4766 4902 ! begin flux computations … … 4873 5009 ENDIF 4874 5010 5011 4875 5012 jp1 = 2 4876 5013 jp0 = 1 … … 4895 5032 DO k = kts+1, ktf 4896 5033 DO i = i_start, i_end 5034 ! Assumes j>j_end_f is ONLY j_end+1 ... 5035 ! fqy(i, k, jp1) = & 5036 ! 0.5*(fzm(k)*rv(i,k,j_end+1)+fzp(k)*rv(i,k-1,j_end+1)) & 5037 ! *(w(i,k,j_end+1)+w(i,k,j_end)) 4897 5038 fqy(i, k, jp1) = & 4898 0.5*(fzm(k)*rv(i,k,j _end+1)+fzp(k)*rv(i,k-1,j_end+1)) &4899 *(w(i,k,j _end+1)+w(i,k,j_end))5039 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)) & 5040 *(w(i,k,j)+w(i,k,j-1)) 4900 5041 ENDDO 4901 5042 ENDDO 4902 5043 k = ktf+1 4903 5044 DO i = i_start, i_end 5045 ! Assumes j>j_end_f is ONLY j_end+1 ... 5046 ! fqy(i, k, jp1) = & 5047 ! 0.5*((2.-fzm(k-1))*rv(i,k-1,j_end+1)-fzp(k-1)*rv(i,k-2,j_end+1)) & 5048 ! *(w(i,k,j_end+1)+w(i,k,j_end)) 4904 5049 fqy(i, k, jp1) = & 4905 0.5*((2.-fzm(k-1))*rv(i,k-1,j _end+1)-fzp(k-1)*rv(i,k-2,j_end+1)) &4906 *(w(i,k,j _end+1)+w(i,k,j_end))5050 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)) & 5051 *(w(i,k,j)+w(i,k,j-1)) 4907 5052 ENDDO 4908 5053 ELSE … … 4925 5070 END IF 4926 5071 5072 ! y flux-divergence into tendency 4927 5073 IF( j > j_start ) THEN 4928 ! y flux-divergence into tendency 5074 4929 5075 DO k = kts+1, ktf+1 4930 5076 DO i = i_start, i_end … … 4933 5079 ENDDO 4934 5080 ENDDO 5081 5082 4935 5083 END IF 4936 5084 … … 4985 5133 i_start = its 4986 5134 i_end = MIN(ite,ide-1) 5135 ! 4987 5136 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts) 4988 5137 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte) … … 5017 5166 5018 5167 ENDDO 5168 5169 ELSE IF ( horz_order == 0 ) THEN 5170 5171 ! Just in case we want to turn horizontal advection off, we can do it 5019 5172 5020 5173 ELSE … … 5501 5654 flux_upwind(q_im1, q_i, cr ) = 0.5*min( 1.0,(cr+abs(cr)))*q_im1 & 5502 5655 +0.5*max(-1.0,(cr-abs(cr)))*q_i 5656 5657 ! flux_upwind(q_im1, q_i, cr ) = 0.5*(1.+sign(1.,cr))*q_im1 & 5658 ! +0.5*(1.-sign(1.,cr))*q_i 5503 5659 ! flux_upwind(q_im1, q_i, cr ) = 0. 5504 5660 … … 5535 5691 IF( config_flags%periodic_x .or. & 5536 5692 config_flags%symmetric_xs .or. & 5537 (its > ids+ 2) ) degrade_xs = .false.5693 (its > ids+3) ) degrade_xs = .false. 5538 5694 IF( config_flags%periodic_x .or. & 5539 5695 config_flags%symmetric_xe .or. & 5540 (ite < ide- 3) ) degrade_xe = .false.5696 (ite < ide-4) ) degrade_xe = .false. 5541 5697 IF( config_flags%periodic_y .or. & 5542 5698 config_flags%symmetric_ys .or. & 5543 (jts > jds+ 2) ) degrade_ys = .false.5699 (jts > jds+3) ) degrade_ys = .false. 5544 5700 IF( config_flags%periodic_y .or. & 5545 5701 config_flags%symmetric_ye .or. & 5546 (jte < jde- 3) ) degrade_ye = .false.5702 (jte < jde-4) ) degrade_ye = .false. 5547 5703 5548 5704 !--------------- y - advection first … … 5560 5716 !-- modify loop bounds if open or specified 5561 5717 5562 IF(degrade_xs) i_start = its 5563 IF(degrade_xe) i_end = MIN(ite,ide-1) 5718 ! IF(degrade_xs) i_start = MAX(its-1,ids-1) 5719 ! IF(degrade_xe) i_end = MIN(ite+1,ide-2) 5720 IF(degrade_xs) i_start = MAX(its-1,ids) 5721 IF(degrade_xe) i_end = MIN(ite+1,ide-1) 5564 5722 5565 5723 IF(degrade_ys) then 5566 j_start = MAX(jts ,jds+1)5724 j_start = MAX(jts-1,jds+1) 5567 5725 j_start_f = jds+3 5568 5726 ENDIF 5569 5727 5570 5728 IF(degrade_ye) then 5571 j_end = MIN(jte ,jde-2)5729 j_end = MIN(jte+1,jde-2) 5572 5730 j_end_f = jde-3 5573 5731 ENDIF … … 5652 5810 ENDDO 5653 5811 5654 ELSE IF ( j == jde-2 ) THEN ! 4th order flux 2 in from north boundary5812 ELSE IF ( j == jde-2 ) THEN ! 3rd or 4th order flux 2 in from north boundary 5655 5813 5656 5814 DO k=kts,ktf … … 5673 5831 ENDIF 5674 5832 5675 5833 ENDDO j_loop_y_flux_6 5676 5834 5677 5835 ! next, x flux … … 5689 5847 !-- modify loop bounds for open and specified b.c 5690 5848 5691 IF(degrade_ys) j_start = jts 5692 IF(degrade_ye) j_end = MIN(jte,jde-1) 5849 ! IF(degrade_ys) j_start = MAX(jts-1,jds+1) 5850 ! IF(degrade_ye) j_end = MIN(jte+1,jde-2) 5851 IF(degrade_ys) j_start = MAX(jts-1,jds) 5852 IF(degrade_ye) j_end = MIN(jte+1,jde-1) 5693 5853 5694 5854 IF(degrade_xs) then 5695 i_start = MAX(ids+1,its )5696 i_start_f = i _start+25855 i_start = MAX(ids+1,its-1) 5856 i_start_f = ids+3 5697 5857 ENDIF 5698 5858 5699 5859 IF(degrade_xe) then 5700 i_end = MIN(ide-2,ite )5860 i_end = MIN(ide-2,ite+1) 5701 5861 i_end_f = ide-3 5702 5862 ENDIF … … 5706 5866 DO j = j_start, j_end 5707 5867 5708 ! 6th order flux5868 ! 5th order flux 5709 5869 5710 5870 DO k=kts,ktf … … 5730 5890 IF( degrade_xs ) THEN 5731 5891 5732 IF( i_start == ids+1 ) THEN ! second order flux next to the boundary 5733 i = ids+1 5734 DO k=kts,ktf 5735 5736 dx = 2./(msft(i,j)+msft(i-1,j))/rdx 5737 mu = 0.5*(mut(i,j)+mut(i-1,j)) 5738 vel = ru(i,k,j)/mu 5739 cr = vel*dt/dx 5740 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) 5741 5742 fqx(i,k,j) = 0.5*(ru(i,k,j)) & 5743 *(field(i,k,j)+field(i-1,k,j)) 5744 5745 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) 5746 5747 ENDDO 5748 ENDIF 5749 5750 i = ids+2 5751 DO k=kts,ktf 5752 dx = 2./(msft(i,j)+msft(i-1,j))/rdx 5753 mu = 0.5*(mut(i,j)+mut(i-1,j)) 5754 vel = ru(i,k,j) 5755 cr = vel*dt/dx/mu 5756 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) 5757 fqx( i,k,j ) = vel*flux4( field(i-2,k,j), field(i-1,k,j), & 5758 field(i ,k,j), field(i+1,k,j), & 5759 vel ) 5760 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) 5761 5762 ENDDO 5763 5764 ENDIF 5765 5766 IF( degrade_xe ) THEN 5767 5768 IF( i_end == ide-2 ) THEN ! second order flux next to the boundary 5769 i = ide-1 5770 DO k=kts,ktf 5771 dx = 2./(msft(i,j)+msft(i-1,j))/rdx 5772 mu = 0.5*(mut(i,j)+mut(i-1,j)) 5773 vel = ru(i,k,j) 5774 cr = vel*dt/dx/mu 5775 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) 5776 fqx(i,k,j) = 0.5*(ru(i,k,j)) & 5777 *(field(i,k,j)+field(i-1,k,j)) 5778 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) 5779 5780 ENDDO 5781 ENDIF 5782 5783 i = ide-2 5784 DO k=kts,ktf 5785 5786 dx = 2./(msft(i,j)+msft(i-1,j))/rdx 5787 mu = 0.5*(mut(i,j)+mut(i-1,j)) 5788 vel = ru(i,k,j) 5789 cr = vel*dt/dx/mu 5790 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) 5791 fqx( i,k,j ) = vel*flux4( field(i-2,k,j), field(i-1,k,j), & 5892 DO i=i_start,i_start_f-1 5893 5894 IF(i == ids+1) THEN ! second order 5895 DO k=kts,ktf 5896 dx = 2./(msft(i,j)+msft(i-1,j))/rdx 5897 mu = 0.5*(mut(i,j)+mut(i-1,j)) 5898 vel = ru(i,k,j)/mu 5899 cr = vel*dt/dx 5900 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) 5901 fqx(i,k,j) = 0.5*(ru(i,k,j)) & 5902 *(field(i,k,j)+field(i-1,k,j)) 5903 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) 5904 ENDDO 5905 ENDIF 5906 5907 IF(i == ids+2) THEN ! fourth order 5908 DO k=kts,ktf 5909 dx = 2./(msft(i,j)+msft(i-1,j))/rdx 5910 mu = 0.5*(mut(i,j)+mut(i-1,j)) 5911 vel = ru(i,k,j) 5912 cr = vel*dt/dx/mu 5913 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) 5914 fqx( i,k,j ) = vel*flux4( field(i-2,k,j), field(i-1,k,j), & 5792 5915 field(i ,k,j), field(i+1,k,j), & 5793 5916 vel ) 5794 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) 5917 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) 5918 ENDDO 5919 ENDIF 5920 5921 ENDDO 5922 5923 ENDIF 5924 5925 IF( degrade_xe ) THEN 5926 5927 DO i = i_end_f+1, i_end+1 5928 5929 IF( i == ide-1 ) THEN ! second order flux next to the boundary 5930 DO k=kts,ktf 5931 dx = 2./(msft(i,j)+msft(i-1,j))/rdx 5932 mu = 0.5*(mut(i,j)+mut(i-1,j)) 5933 vel = ru(i,k,j) 5934 cr = vel*dt/dx/mu 5935 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) 5936 fqx(i,k,j) = 0.5*(ru(i,k,j)) & 5937 *(field(i,k,j)+field(i-1,k,j)) 5938 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) 5939 ENDDO 5940 ENDIF 5941 5942 5943 IF( i == ide-2 ) THEN ! fourth order flux one in from the boundary 5944 DO k=kts,ktf 5945 dx = 2./(msft(i,j)+msft(i-1,j))/rdx 5946 mu = 0.5*(mut(i,j)+mut(i-1,j)) 5947 vel = ru(i,k,j) 5948 cr = vel*dt/dx/mu 5949 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) 5950 fqx( i,k,j ) = vel*flux4( field(i-2,k,j), field(i-1,k,j), & 5951 field(i ,k,j), field(i+1,k,j), & 5952 vel ) 5953 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) 5954 ENDDO 5955 ENDIF 5795 5956 5796 5957 ENDDO … … 5806 5967 IF( config_flags%periodic_x .or. & 5807 5968 config_flags%symmetric_xs .or. & 5808 (its > ids+ 2) ) degrade_xs = .false.5969 (its > ids+3) ) degrade_xs = .false. 5809 5970 IF( config_flags%periodic_x .or. & 5810 5971 config_flags%symmetric_xe .or. & 5811 (ite < ide- 3) ) degrade_xe = .false.5972 (ite < ide-4) ) degrade_xe = .false. 5812 5973 IF( config_flags%periodic_y .or. & 5813 5974 config_flags%symmetric_ys .or. & 5814 (jts > jds+ 2) ) degrade_ys = .false.5975 (jts > jds+3) ) degrade_ys = .false. 5815 5976 IF( config_flags%periodic_y .or. & 5816 5977 config_flags%symmetric_ye .or. & 5817 (jte < jde- 3) ) degrade_ye = .false.5978 (jte < jde-4) ) degrade_ye = .false. 5818 5979 5819 5980 !--------------- y - advection first … … 5831 5992 !-- modify loop bounds if open or specified 5832 5993 5833 IF(degrade_xs) i_start = its 5834 IF(degrade_xe) i_end = MIN(ite,ide-1) 5994 ! IF(degrade_xs) i_start = MAX(its-1,ids-1) 5995 ! IF(degrade_xe) i_end = MIN(ite+1,ide-2) 5996 IF(degrade_xs) i_start = MAX(its-1,ids) 5997 IF(degrade_xe) i_end = MIN(ite+1,ide-1) 5835 5998 5836 5999 IF(degrade_ys) then 5837 j_start = MAX(jts ,jds+1)6000 j_start = MAX(jts-1,jds+1) 5838 6001 j_start_f = jds+3 5839 6002 ENDIF 5840 6003 5841 6004 IF(degrade_ye) then 5842 j_end = MIN(jte ,jde-2)6005 j_end = MIN(jte+1,jde-2) 5843 6006 j_end_f = jde-3 5844 6007 ENDIF … … 5960 6123 !-- modify loop bounds for open and specified b.c 5961 6124 5962 IF(degrade_ys) j_start = jts 5963 IF(degrade_ye) j_end = MIN(jte,jde-1) 6125 ! IF(degrade_ys) j_start = MAX(jts-1,jds+1) 6126 ! IF(degrade_ye) j_end = MIN(jte+1,jde-2) 6127 IF(degrade_ys) j_start = MAX(jts-1,jds) 6128 IF(degrade_ye) j_end = MIN(jte+1,jde-1) 5964 6129 5965 6130 IF(degrade_xs) then 5966 i_start = MAX(ids+1,its )5967 i_start_f = i _start+26131 i_start = MAX(ids+1,its-1) 6132 i_start_f = ids+3 5968 6133 ENDIF 5969 6134 5970 6135 IF(degrade_xe) then 5971 i_end = MIN(ide-2,ite )6136 i_end = MIN(ide-2,ite+1) 5972 6137 i_end_f = ide-3 5973 6138 ENDIF … … 6001 6166 IF( degrade_xs ) THEN 6002 6167 6003 IF( i_start == ids+1 ) THEN ! second order flux next to the boundary 6004 i = ids+1 6005 DO k=kts,ktf 6006 6007 dx = 2./(msft(i,j)+msft(i-1,j))/rdx 6008 mu = 0.5*(mut(i,j)+mut(i-1,j)) 6009 vel = ru(i,k,j)/mu 6010 cr = vel*dt/dx 6011 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) 6012 6013 fqx(i,k,j) = 0.5*(ru(i,k,j)) & 6014 *(field(i,k,j)+field(i-1,k,j)) 6015 6016 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) 6017 6018 ENDDO 6019 ENDIF 6020 6021 i = ids+2 6022 DO k=kts,ktf 6023 dx = 2./(msft(i,j)+msft(i-1,j))/rdx 6024 mu = 0.5*(mut(i,j)+mut(i-1,j)) 6025 vel = ru(i,k,j) 6026 cr = vel*dt/dx/mu 6027 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) 6028 fqx( i,k,j ) = vel*flux3( field(i-2,k,j), field(i-1,k,j), & 6029 field(i ,k,j), field(i+1,k,j), & 6030 vel ) 6031 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) 6032 6033 ENDDO 6034 6035 ENDIF 6036 6037 IF( degrade_xe ) THEN 6038 6039 IF( i_end == ide-2 ) THEN ! second order flux next to the boundary 6040 i = ide-1 6041 DO k=kts,ktf 6042 dx = 2./(msft(i,j)+msft(i-1,j))/rdx 6043 mu = 0.5*(mut(i,j)+mut(i-1,j)) 6044 vel = ru(i,k,j) 6045 cr = vel*dt/dx/mu 6046 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) 6047 fqx(i,k,j) = 0.5*(ru(i,k,j)) & 6048 *(field(i,k,j)+field(i-1,k,j)) 6049 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) 6050 6051 ENDDO 6052 ENDIF 6053 6054 i = ide-2 6055 DO k=kts,ktf 6056 6057 dx = 2./(msft(i,j)+msft(i-1,j))/rdx 6058 mu = 0.5*(mut(i,j)+mut(i-1,j)) 6059 vel = ru(i,k,j) 6060 cr = vel*dt/dx/mu 6061 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) 6062 fqx( i,k,j ) = vel*flux3( field(i-2,k,j), field(i-1,k,j), & 6168 DO i=i_start,i_start_f-1 6169 6170 IF(i == ids+1) THEN ! second order 6171 DO k=kts,ktf 6172 dx = 2./(msft(i,j)+msft(i-1,j))/rdx 6173 mu = 0.5*(mut(i,j)+mut(i-1,j)) 6174 vel = ru(i,k,j)/mu 6175 cr = vel*dt/dx 6176 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) 6177 fqx(i,k,j) = 0.5*(ru(i,k,j)) & 6178 *(field(i,k,j)+field(i-1,k,j)) 6179 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) 6180 ENDDO 6181 ENDIF 6182 6183 IF(i == ids+2) THEN ! third order 6184 DO k=kts,ktf 6185 dx = 2./(msft(i,j)+msft(i-1,j))/rdx 6186 mu = 0.5*(mut(i,j)+mut(i-1,j)) 6187 vel = ru(i,k,j) 6188 cr = vel*dt/dx/mu 6189 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) 6190 fqx( i,k,j ) = vel*flux3( field(i-2,k,j), field(i-1,k,j), & 6063 6191 field(i ,k,j), field(i+1,k,j), & 6064 6192 vel ) 6065 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) 6193 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) 6194 ENDDO 6195 ENDIF 6196 6197 ENDDO 6198 6199 ENDIF 6200 6201 IF( degrade_xe ) THEN 6202 6203 DO i = i_end_f+1, i_end+1 6204 6205 IF( i == ide-1 ) THEN ! second order flux next to the boundary 6206 DO k=kts,ktf 6207 dx = 2./(msft(i,j)+msft(i-1,j))/rdx 6208 mu = 0.5*(mut(i,j)+mut(i-1,j)) 6209 vel = ru(i,k,j) 6210 cr = vel*dt/dx/mu 6211 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) 6212 fqx(i,k,j) = 0.5*(ru(i,k,j)) & 6213 *(field(i,k,j)+field(i-1,k,j)) 6214 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) 6215 ENDDO 6216 ENDIF 6217 6218 6219 IF( i == ide-2 ) THEN ! third order flux one in from the boundary 6220 DO k=kts,ktf 6221 dx = 2./(msft(i,j)+msft(i-1,j))/rdx 6222 mu = 0.5*(mut(i,j)+mut(i-1,j)) 6223 vel = ru(i,k,j) 6224 cr = vel*dt/dx/mu 6225 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) 6226 fqx( i,k,j ) = vel*flux3( field(i-2,k,j), field(i-1,k,j), & 6227 field(i ,k,j), field(i+1,k,j), & 6228 vel ) 6229 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) 6230 ENDDO 6231 ENDIF 6066 6232 6067 6233 ENDDO … … 6275 6441 IF( config_flags%periodic_x .or. & 6276 6442 config_flags%symmetric_xs .or. & 6277 (its > ids+ 1) ) degrade_xs = .false.6443 (its > ids+2) ) degrade_xs = .false. 6278 6444 IF( config_flags%periodic_x .or. & 6279 6445 config_flags%symmetric_xe .or. & 6280 (ite < ide- 2) ) degrade_xe = .false.6446 (ite < ide-1) ) degrade_xe = .false. 6281 6447 IF( config_flags%periodic_y .or. & 6282 6448 config_flags%symmetric_ys .or. & 6283 (jts > jds+ 1) ) degrade_ys = .false.6449 (jts > jds+2) ) degrade_ys = .false. 6284 6450 IF( config_flags%periodic_y .or. & 6285 6451 config_flags%symmetric_ye .or. & 6286 (jte < jde- 2) ) degrade_ye = .false.6452 (jte < jde-1) ) degrade_ye = .false. 6287 6453 6288 6454 !--------------- y - advection first … … 6475 6641 IF( config_flags%periodic_x .or. & 6476 6642 config_flags%symmetric_xs .or. & 6477 (its > ids ) ) degrade_xs = .false.6643 (its > ids+1) ) degrade_xs = .false. 6478 6644 IF( config_flags%periodic_x .or. & 6479 6645 config_flags%symmetric_xe .or. & 6480 (ite < ide- 1) ) degrade_xe = .false.6646 (ite < ide-2) ) degrade_xe = .false. 6481 6647 IF( config_flags%periodic_y .or. & 6482 6648 config_flags%symmetric_ys .or. & 6483 (jts > jds ) ) degrade_ys = .false.6649 (jts > jds+1) ) degrade_ys = .false. 6484 6650 IF( config_flags%periodic_y .or. & 6485 6651 config_flags%symmetric_ye .or. & 6486 (jte < jde- 1) ) degrade_ye = .false.6652 (jte < jde-2) ) degrade_ye = .false. 6487 6653 6488 6654 !-- y flux compute; these bounds are for periodic and sym b.c. … … 6530 6696 cr = vel*dt/dx/mu 6531 6697 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) 6532 fqx( i,k,j ) = vel*flux3( field(i-2,k,j), field(i-1,k,j),&6533 field(i ,k,j), field(i+1,k,j), &6534 vel ) 6698 fqx( i,k,j ) = 0.5*ru(i,k,j)* & 6699 (field(i,k,j)+field(i-1,k,j)) 6700 6535 6701 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) 6536 6702 ENDDO … … 6538 6704 ENDDO 6539 6705 6540 !--- end of 3nd order horizontal flux calculation6706 !--- end of 2nd order horizontal flux calculation 6541 6707 6542 6708 ELSE … … 6633 6799 !-- loop bounds for open or specified conditions 6634 6800 6635 IF(degrade_xs) i_start = its6636 IF(degrade_xe) i_end = MIN(ite ,ide-1)6637 IF(degrade_ys) j_start = jts6638 IF(degrade_ye) j_end = MIN(jte ,jde-1)6801 IF(degrade_xs) i_start = MAX(its-1,ids) 6802 IF(degrade_xe) i_end = MIN(ite+1,ide-1) 6803 IF(degrade_ys) j_start = MAX(jts-1,jds) 6804 IF(degrade_ye) j_end = MIN(jte+1,jde-1) 6639 6805 6640 6806 vert_order_test : IF (vert_order == 6) THEN … … 6931 7097 !-- loop bounds for open or specified conditions 6932 7098 6933 IF(degrade_xs) i_start = its6934 IF(degrade_xe) i_end = MIN(ite ,ide-1)6935 IF(degrade_ys) j_start = jts6936 IF(degrade_ye) j_end = MIN(jte ,jde-1)7099 IF(degrade_xs) i_start = MAX(its-1,ids) 7100 IF(degrade_xe) i_end = MIN(ite+1,ide-1) 7101 IF(degrade_ys) j_start = MAX(jts-1,jds) 7102 IF(degrade_ye) j_end = MIN(jte+1,jde-1) 6937 7103 6938 7104 IF(config_flags%specified .or. config_flags%nested) THEN 6939 IF (degrade_xs) i_start = MAX(its ,ids+1)6940 IF (degrade_xe) i_end = MIN(ite ,ide-2)6941 IF (degrade_ys) j_start = MAX(jts ,jds+1)6942 IF (degrade_ye) j_end = MIN(jte ,jde-2)7105 IF (degrade_xs) i_start = MAX(its-1,ids+1) 7106 IF (degrade_xe) i_end = MIN(ite+1,ide-2) 7107 IF (degrade_ys) j_start = MAX(jts-1,jds+1) 7108 IF (degrade_ye) j_end = MIN(jte+1,jde-2) 6943 7109 END IF 6944 7110 6945 7111 IF(config_flags%open_xs) THEN 6946 IF (degrade_xs) i_start = MAX(its ,ids+1)7112 IF (degrade_xs) i_start = MAX(its-1,ids+1) 6947 7113 END IF 6948 7114 IF(config_flags%open_xe) THEN 6949 IF (degrade_xe) i_end = MIN(ite ,ide-2)7115 IF (degrade_xe) i_end = MIN(ite+1,ide-2) 6950 7116 END IF 6951 7117 IF(config_flags%open_ys) THEN 6952 IF (degrade_ys) j_start = MAX(jts ,jds+1)7118 IF (degrade_ys) j_start = MAX(jts-1,jds+1) 6953 7119 END IF 6954 7120 IF(config_flags%open_ye) THEN 6955 IF (degrade_ye) j_end = MIN(jte ,jde-2)7121 IF (degrade_ye) j_end = MIN(jte+1,jde-2) 6956 7122 END IF 6957 7123 … … 6962 7128 DO i=i_start, i_end 6963 7129 6964 ph_low = (mub(i,j)+mu_old(i,j))*field_old(i,k,j) & 6965 - dt*( msft(i,j)*( rdx*(fqxl(i+1,k,j)-fqxl(i,k,j)) & 6966 +rdy*(fqyl(i,k,j+1)-fqyl(i,k,j)) ) & 6967 +rdzw(k)*(fqzl(i,k+1,j)-fqzl(i,k,j)) ) 6968 6969 flux_out = dt*(msft(i,j)*( rdx*( max(0.,fqx (i+1,k,j)) & 7130 ph_low = (mub(i,j)+mu_old(i,j))*field_old(i,k,j) & 7131 - dt*( msft(i,j)*( & 7132 rdx*(fqxl(i+1,k,j)-fqxl(i,k,j)) + & 7133 rdy*(fqyl(i,k,j+1)-fqyl(i,k,j)) ) & 7134 + rdzw(k)*(fqzl(i,k+1,j)-fqzl(i,k,j)) ) 7135 7136 flux_out = dt*( msft(i,j) *( & 7137 rdx*( max(0.,fqx (i+1,k,j)) & 6970 7138 -min(0.,fqx (i ,k,j)) ) & 6971 7139 +rdy*( max(0.,fqy (i,k,j+1)) & 6972 7140 -min(0.,fqy (i,k,j )) ) ) & 6973 +rdzw(k)*( min(0.,fqz (i,k+1,j)) &7141 + rdzw(k)*( min(0.,fqz (i,k+1,j)) & 6974 7142 -max(0.,fqz (i,k ,j)) ) ) 6975 7143 … … 7015 7183 ! x flux divergence 7016 7184 ! 7017 IF(degrade_xs) i_start = i_start + 17018 IF(degrade_xe) i_end = i_end - 17185 IF(degrade_xs) i_start = MAX(its,ids+1) 7186 IF(degrade_xe) i_end = MIN(ite,ide-2) 7019 7187 7020 7188 DO j = j_start, j_end … … 7023 7191 7024 7192 tendency (i,k,j) = tendency(i,k,j) & 7025 7193 - msft(i,j)*( rdx*( fqx (i+1,k,j)-fqx (i,k,j) & 7026 7194 +fqxl(i+1,k,j)-fqxl(i,k,j)) ) 7027 7195 … … 7034 7202 i_start = its 7035 7203 i_end = MIN(ite,ide-1) 7036 IF(degrade_ys) j_start = j_start + 17037 IF(degrade_ye) j_end = j_end - 17204 IF(degrade_ys) j_start = MAX(jts,jds+1) 7205 IF(degrade_ye) j_end = MIN(jte,jde-2) 7038 7206 7039 7207 DO j = j_start, j_end … … 7042 7210 7043 7211 tendency (i,k,j) = tendency(i,k,j) & 7044 7212 - msft(i,j)*( rdy*( fqy (i,k,j+1)-fqy (i,k,j) & 7045 7213 +fqyl(i,k,j+1)-fqyl(i,k,j)) ) 7046 7214 … … 7053 7221 !---------------------------------------------------------------- 7054 7222 7223 SUBROUTINE advect_scalar_mono ( field, field_old, tendency, & 7224 ru, rv, rom, & 7225 mut, mub, mu_old, & 7226 config_flags, & 7227 msfu, msfv, & 7228 msft, & 7229 fzm, fzp, & 7230 rdx, rdy, rdzw, dt, & 7231 ids, ide, jds, jde, kds, kde, & 7232 ims, ime, jms, jme, kms, kme, & 7233 its, ite, jts, jte, kts, kte ) 7234 7235 ! monotonic advection option 7236 ! for scalars in WRF RK3 advection. This version is memory intensive -> 7237 ! we save 3d arrays of x, y and z both high and low order fluxes 7238 ! (six in all). Alternatively, we could sweep in a direction 7239 ! and lower the cost considerably. 7240 7241 ! uses the Smolarkiewicz MWR 1989 approach, with addition of first-order 7242 ! fluxes initially 7243 7244 IMPLICIT NONE 7245 7246 ! Input data 7247 7248 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags 7249 7250 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & 7251 ims, ime, jms, jme, kms, kme, & 7252 its, ite, jts, jte, kts, kte 7253 7254 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field, & 7255 field_old, & 7256 ru, & 7257 rv, & 7258 rom 7259 7260 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut, mub, mu_old 7261 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency 7262 7263 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfu, & 7264 msfv, & 7265 msft 7266 7267 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, & 7268 fzp, & 7269 rdzw 7270 7271 REAL , INTENT(IN ) :: rdx, & 7272 rdy, & 7273 dt 7274 7275 ! Local data 7276 7277 INTEGER :: i, j, k, itf, jtf, ktf 7278 INTEGER :: i_start, i_end, j_start, j_end 7279 INTEGER :: i_start_f, i_end_f, j_start_f, j_end_f 7280 INTEGER :: jmin, jmax, jp, jm, imin, imax 7281 7282 REAL :: mrdx, mrdy, ub, vb, uw, vw, mu 7283 REAL , DIMENSION(its:ite, kts:kte) :: vflux 7284 7285 7286 ! storage for high and low order fluxes 7287 7288 REAL, DIMENSION( its-2:ite+2, kts:kte, jts-2:jte+2 ) :: fqx, fqy, fqz 7289 REAL, DIMENSION( its-2:ite+2, kts:kte, jts-2:jte+2 ) :: fqxl, fqyl, fqzl 7290 REAL, DIMENSION( its-2:ite+2, kts:kte, jts-2:jte+2 ) :: qmin, qmax 7291 REAL, DIMENSION( its-2:ite+2, kts:kte, jts-2:jte+2 ) :: scale_in, scale_out 7292 REAL :: ph_upwind 7293 7294 INTEGER :: horz_order, vert_order 7295 7296 LOGICAL :: degrade_xs, degrade_ys 7297 LOGICAL :: degrade_xe, degrade_ye 7298 7299 INTEGER :: jp1, jp0, jtmp 7300 7301 REAL :: flux_out, ph_low, flux_in, ph_hi, scale 7302 REAL, PARAMETER :: eps=1.e-20 7303 7304 7305 ! definition of flux operators, 3rd, 4rth, 5th or 6th order 7306 7307 REAL :: flux3, flux4, flux5, flux6, flux_upwind 7308 REAL :: q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua, vel, cr 7309 7310 flux4(q_im2, q_im1, q_i, q_ip1, ua) = & 7311 (7./12.)*(q_i + q_im1) - (1./12.)*(q_ip1 + q_im2) 7312 7313 flux3(q_im2, q_im1, q_i, q_ip1, ua) = & 7314 flux4(q_im2, q_im1, q_i, q_ip1, ua) + & 7315 sign(1.,ua)*(1./12.)*((q_ip1 - q_im2)-3.*(q_i-q_im1)) 7316 7317 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = & 7318 (37./60.)*(q_i+q_im1) - (2./15.)*(q_ip1+q_im2) & 7319 +(1./60.)*(q_ip2+q_im3) 7320 7321 flux5(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = & 7322 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) & 7323 -sign(1.,ua)*(1./60.)*( & 7324 (q_ip2-q_im3)-5.*(q_ip1-q_im2)+10.*(q_i-q_im1) ) 7325 7326 ! flux_upwind(q_im1, q_i, cr ) = 0. 7327 flux_upwind(q_im1, q_i, cr ) = 0.5*(1.+sign(1.,cr))*q_im1 & 7328 +0.5*(1.-sign(1.,cr))*q_i 7329 7330 LOGICAL, PARAMETER :: mono_limit = .true. 7331 7332 ! set order for the advection schemes 7333 7334 ktf=MIN(kte,kde-1) 7335 horz_order = config_flags%h_sca_adv_order 7336 vert_order = config_flags%v_sca_adv_order 7337 7338 do j=jts-2,jte+2 7339 do k=kts,kte 7340 do i=its-2,ite+2 7341 qmin(i,k,j) = field_old(i,k,j) 7342 qmax(i,k,j) = field_old(i,k,j) 7343 scale_in(i,k,j) = 1. 7344 scale_out(i,k,j) = 1. 7345 fqx(i,k,j) = 0. 7346 fqy(i,k,j) = 0. 7347 fqz(i,k,j) = 0. 7348 fqxl(i,k,j) = 0. 7349 fqyl(i,k,j) = 0. 7350 fqzl(i,k,j) = 0. 7351 enddo 7352 enddo 7353 enddo 7354 7355 ! begin with horizontal flux divergence 7356 ! here is the choice of flux operators 7357 7358 7359 horizontal_order_test : IF( horz_order == 5 ) THEN 7360 7361 ! determine boundary mods for flux operators 7362 ! We degrade the flux operators from 3rd/4rth order 7363 ! to second order one gridpoint in from the boundaries for 7364 ! all boundary conditions except periodic and symmetry - these 7365 ! conditions have boundary zone data fill for correct application 7366 ! of the higher order flux stencils 7367 7368 degrade_xs = .true. 7369 degrade_xe = .true. 7370 degrade_ys = .true. 7371 degrade_ye = .true. 7372 7373 IF( config_flags%periodic_x .or. & 7374 config_flags%symmetric_xs .or. & 7375 (its > ids+3) ) degrade_xs = .false. 7376 IF( config_flags%periodic_x .or. & 7377 config_flags%symmetric_xe .or. & 7378 (ite < ide-4) ) degrade_xe = .false. 7379 IF( config_flags%periodic_y .or. & 7380 config_flags%symmetric_ys .or. & 7381 (jts > jds+3) ) degrade_ys = .false. 7382 IF( config_flags%periodic_y .or. & 7383 config_flags%symmetric_ye .or. & 7384 (jte < jde-4) ) degrade_ye = .false. 7385 7386 !--------------- y - advection first 7387 7388 !-- y flux compute; these bounds are for periodic and sym b.c. 7389 7390 ktf=MIN(kte,kde-1) 7391 i_start = its-1 7392 i_end = MIN(ite,ide-1)+1 7393 j_start = jts-1 7394 j_end = MIN(jte,jde-1)+1 7395 j_start_f = j_start 7396 j_end_f = j_end+1 7397 7398 !-- modify loop bounds if open or specified 7399 7400 ! WCS 20090218 7401 ! IF(degrade_xs) i_start = its 7402 ! IF(degrade_xe) i_end = MIN(ite,ide-1) 7403 IF(degrade_xs) i_start = MAX(its-1,ids) 7404 IF(degrade_xe) i_end = MIN(ite+1,ide-1) 7405 7406 ! WCS 20090218 7407 ! IF(degrade_ys) then 7408 ! j_start = MAX(jts,jds+1) 7409 ! j_start_f = jds+3 7410 ! ENDIF 7411 ! 7412 ! IF(degrade_ye) then 7413 ! j_end = MIN(jte,jde-2) 7414 ! j_end_f = jde-3 7415 ! ENDIF 7416 7417 IF(degrade_ys) then 7418 j_start = MAX(jts-1,jds+1) 7419 j_start_f = jds+3 7420 ENDIF 7421 7422 IF(degrade_ye) then 7423 j_end = MIN(jte+1,jde-2) 7424 j_end_f = jde-3 7425 ENDIF 7426 7427 ! compute fluxes, 5th order 7428 7429 j_loop_y_flux_5 : DO j = j_start, j_end+1 7430 7431 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN ! use full stencil 7432 7433 DO k=kts,ktf 7434 DO i = i_start, i_end 7435 7436 vel = rv(i,k,j) 7437 cr = vel 7438 fqyl(i,k,j) = vel*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), vel) 7439 7440 fqy( i, k, j ) = vel*flux5( & 7441 field(i,k,j-3), field(i,k,j-2), field(i,k,j-1), & 7442 field(i,k,j ), field(i,k,j+1), field(i,k,j+2), vel ) 7443 7444 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j) 7445 7446 if(cr.gt. 0) then 7447 qmax(i,k,j) = amax1(qmax(i,k,j),field_old(i,k,j-1)) 7448 qmin(i,k,j) = amin1(qmin(i,k,j),field_old(i,k,j-1)) 7449 else 7450 qmax(i,k,j-1) = amax1(qmax(i,k,j-1),field_old(i,k,j)) 7451 qmin(i,k,j-1) = amin1(qmin(i,k,j-1),field_old(i,k,j)) 7452 end if 7453 7454 ENDDO 7455 ENDDO 7456 7457 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary 7458 7459 DO k=kts,ktf 7460 DO i = i_start, i_end 7461 7462 vel = rv(i,k,j) 7463 cr = vel 7464 fqyl(i,k,j) = vel*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr) 7465 7466 fqy(i,k, j) = 0.5*rv(i,k,j)* & 7467 (field(i,k,j)+field(i,k,j-1)) 7468 7469 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j) 7470 7471 if(cr.gt. 0) then 7472 qmax(i,k,j) = amax1(qmax(i,k,j),field_old(i,k,j-1)) 7473 qmin(i,k,j) = amin1(qmin(i,k,j),field_old(i,k,j-1)) 7474 else 7475 qmax(i,k,j-1) = amax1(qmax(i,k,j-1),field_old(i,k,j)) 7476 qmin(i,k,j-1) = amin1(qmin(i,k,j-1),field_old(i,k,j)) 7477 end if 7478 7479 ENDDO 7480 ENDDO 7481 7482 ELSE IF ( j == jds+2 ) THEN ! third of 4rth order flux 2 in from south boundary 7483 7484 DO k=kts,ktf 7485 DO i = i_start, i_end 7486 7487 vel = rv(i,k,j) 7488 cr = vel 7489 fqyl(i,k,j) = vel*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr) 7490 7491 fqy( i, k, j ) = vel*flux3( & 7492 field(i,k,j-2),field(i,k,j-1),field(i,k,j),field(i,k,j+1),vel ) 7493 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j) 7494 7495 if(cr.gt. 0) then 7496 qmax(i,k,j) = amax1(qmax(i,k,j),field_old(i,k,j-1)) 7497 qmin(i,k,j) = amin1(qmin(i,k,j),field_old(i,k,j-1)) 7498 else 7499 qmax(i,k,j-1) = amax1(qmax(i,k,j-1),field_old(i,k,j)) 7500 qmin(i,k,j-1) = amin1(qmin(i,k,j-1),field_old(i,k,j)) 7501 end if 7502 7503 ENDDO 7504 ENDDO 7505 7506 ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary 7507 7508 DO k=kts,ktf 7509 DO i = i_start, i_end 7510 7511 vel = rv(i,k,j) 7512 cr = vel 7513 fqyl(i,k,j) = vel*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr) 7514 7515 fqy(i, k, j ) = 0.5*rv(i,k,j)* & 7516 (field(i,k,j)+field(i,k,j-1)) 7517 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j) 7518 7519 if(cr.gt. 0) then 7520 qmax(i,k,j) = amax1(qmax(i,k,j),field_old(i,k,j-1)) 7521 qmin(i,k,j) = amin1(qmin(i,k,j),field_old(i,k,j-1)) 7522 else 7523 qmax(i,k,j-1) = amax1(qmax(i,k,j-1),field_old(i,k,j)) 7524 qmin(i,k,j-1) = amin1(qmin(i,k,j-1),field_old(i,k,j)) 7525 end if 7526 7527 ENDDO 7528 ENDDO 7529 7530 ELSE IF ( j == jde-2 ) THEN ! 3rd or 4rth order flux 2 in from north boundary 7531 7532 DO k=kts,ktf 7533 DO i = i_start, i_end 7534 7535 vel = rv(i,k,j) 7536 cr = vel 7537 fqyl(i,k,j) = vel*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr) 7538 7539 fqy( i, k, j) = vel*flux3( & 7540 field(i,k,j-2),field(i,k,j-1), & 7541 field(i,k,j),field(i,k,j+1),vel ) 7542 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j) 7543 7544 if(cr.gt. 0) then 7545 qmax(i,k,j) = amax1(qmax(i,k,j),field_old(i,k,j-1)) 7546 qmin(i,k,j) = amin1(qmin(i,k,j),field_old(i,k,j-1)) 7547 else 7548 qmax(i,k,j-1) = amax1(qmax(i,k,j-1),field_old(i,k,j)) 7549 qmin(i,k,j-1) = amin1(qmin(i,k,j-1),field_old(i,k,j)) 7550 end if 7551 7552 ENDDO 7553 ENDDO 7554 7555 ENDIF 7556 7557 ENDDO j_loop_y_flux_5 7558 7559 ! next, x flux 7560 7561 !-- these bounds are for periodic and sym conditions 7562 7563 i_start = its-1 7564 i_end = MIN(ite,ide-1)+1 7565 i_start_f = i_start 7566 i_end_f = i_end+1 7567 7568 j_start = jts-1 7569 j_end = MIN(jte,jde-1)+1 7570 7571 !-- modify loop bounds for open and specified b.c 7572 7573 ! WCS 20090218 7574 ! IF(degrade_ys) j_start = jts 7575 ! IF(degrade_ye) j_end = MIN(jte,jde-1) 7576 IF(degrade_ys) j_start = MAX(jts-1,jds) 7577 IF(degrade_ye) j_end = MIN(jte+1,jde-1) 7578 7579 ! WCS 20090218 7580 ! IF(degrade_xs) then 7581 ! i_start = MAX(ids+1,its) 7582 ! i_start_f = i_start+2 7583 ! ENDIF 7584 7585 ! IF(degrade_xe) then 7586 ! i_end = MIN(ide-2,ite) 7587 ! i_end_f = ide-3 7588 ! ENDIF 7589 7590 IF(degrade_xs) then 7591 i_start = MAX(ids+1,its-1) 7592 i_start_f = ids+3 7593 ENDIF 7594 7595 IF(degrade_xe) then 7596 i_end = MIN(ide-2,ite+1) 7597 i_end_f = ide-3 7598 ENDIF 7599 7600 ! compute fluxes 7601 7602 DO j = j_start, j_end 7603 7604 ! 5th or 6th order flux 7605 7606 DO k=kts,ktf 7607 DO i = i_start_f, i_end_f 7608 7609 vel = ru(i,k,j) 7610 cr = vel 7611 fqxl(i,k,j) = vel*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) 7612 7613 fqx( i,k,j ) = vel*flux5( field(i-3,k,j), field(i-2,k,j), & 7614 field(i-1,k,j), field(i ,k,j), & 7615 field(i+1,k,j), field(i+2,k,j), & 7616 vel ) 7617 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) 7618 7619 if(cr.gt. 0) then 7620 qmax(i,k,j) = amax1(qmax(i,k,j),field_old(i-1,k,j)) 7621 qmin(i,k,j) = amin1(qmin(i,k,j),field_old(i-1,k,j)) 7622 else 7623 qmax(i-1,k,j) = amax1(qmax(i-1,k,j),field_old(i,k,j)) 7624 qmin(i-1,k,j) = amin1(qmin(i-1,k,j),field_old(i,k,j)) 7625 end if 7626 7627 ENDDO 7628 ENDDO 7629 7630 ! lower order fluxes close to boundaries (if not periodic or symmetric) 7631 7632 ! WCS 20090218 degrade_xs and xe recoded 7633 7634 IF( degrade_xs ) THEN 7635 7636 DO i=i_start,i_start_f-1 7637 7638 IF(i == ids+1) THEN ! second order 7639 DO k=kts,ktf 7640 vel = ru(i,k,j) 7641 cr = vel 7642 fqxl(i,k,j) = vel*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) 7643 7644 fqx(i,k,j) = 0.5*(ru(i,k,j)) & 7645 *(field(i,k,j)+field(i-1,k,j)) 7646 7647 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) 7648 7649 if(cr.gt. 0) then 7650 qmax(i,k,j) = amax1(qmax(i,k,j),field_old(i-1,k,j)) 7651 qmin(i,k,j) = amin1(qmin(i,k,j),field_old(i-1,k,j)) 7652 else 7653 qmax(i-1,k,j) = amax1(qmax(i-1,k,j),field_old(i,k,j)) 7654 qmin(i-1,k,j) = amin1(qmin(i-1,k,j),field_old(i,k,j)) 7655 end if 7656 ENDDO 7657 ENDIF 7658 7659 IF(i == ids+2) THEN ! third order 7660 DO k=kts,ktf 7661 vel = ru(i,k,j) 7662 cr = vel 7663 fqxl(i,k,j) = vel*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) 7664 fqx( i,k,j ) = vel*flux3( field(i-2,k,j), field(i-1,k,j), & 7665 field(i ,k,j), field(i+1,k,j), & 7666 vel ) 7667 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) 7668 7669 if(cr.gt. 0) then 7670 qmax(i,k,j) = amax1(qmax(i,k,j),field_old(i-1,k,j)) 7671 qmin(i,k,j) = amin1(qmin(i,k,j),field_old(i-1,k,j)) 7672 else 7673 qmax(i-1,k,j) = amax1(qmax(i-1,k,j),field_old(i,k,j)) 7674 qmin(i-1,k,j) = amin1(qmin(i-1,k,j),field_old(i,k,j)) 7675 end if 7676 ENDDO 7677 ENDIF 7678 7679 ENDDO 7680 7681 ENDIF 7682 7683 IF( degrade_xe ) THEN 7684 7685 DO i = i_end_f+1, i_end+1 7686 7687 IF( i == ide-1 ) THEN ! second order flux next to the boundary 7688 DO k=kts,ktf 7689 vel = ru(i,k,j) 7690 cr = vel 7691 fqxl(i,k,j) = vel*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) 7692 fqx(i,k,j) = 0.5*(ru(i,k,j)) & 7693 *(field(i,k,j)+field(i-1,k,j)) 7694 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) 7695 7696 if(cr.gt. 0) then 7697 qmax(i,k,j) = amax1(qmax(i,k,j),field_old(i-1,k,j)) 7698 qmin(i,k,j) = amin1(qmin(i,k,j),field_old(i-1,k,j)) 7699 else 7700 qmax(i-1,k,j) = amax1(qmax(i-1,k,j),field_old(i,k,j)) 7701 qmin(i-1,k,j) = amin1(qmin(i-1,k,j),field_old(i,k,j)) 7702 end if 7703 ENDDO 7704 ENDIF 7705 7706 IF( i == ide-2 ) THEN ! third order flux one in from the boundary 7707 DO k=kts,ktf 7708 vel = ru(i,k,j) 7709 cr = vel 7710 fqxl(i,k,j) = vel*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) 7711 fqx( i,k,j ) = vel*flux3( field(i-2,k,j), field(i-1,k,j), & 7712 field(i ,k,j), field(i+1,k,j), & 7713 vel ) 7714 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) 7715 7716 if(cr.gt. 0) then 7717 qmax(i,k,j) = amax1(qmax(i,k,j),field_old(i-1,k,j)) 7718 qmin(i,k,j) = amin1(qmin(i,k,j),field_old(i-1,k,j)) 7719 else 7720 qmax(i-1,k,j) = amax1(qmax(i-1,k,j),field_old(i,k,j)) 7721 qmin(i-1,k,j) = amin1(qmin(i-1,k,j),field_old(i,k,j)) 7722 end if 7723 ENDDO 7724 ENDIF 7725 ENDDO 7726 ENDIF 7727 7728 ENDDO ! enddo for outer J loop 7729 7730 ELSE 7731 7732 WRITE ( wrf_err_message , * ) 'module_advect: advect_scalar_mono, h_order not known ',horz_order 7733 CALL wrf_error_fatal ( TRIM( wrf_err_message ) ) 7734 7735 ENDIF horizontal_order_test 7736 7737 ! pick up the rest of the horizontal radiation boundary conditions. 7738 ! (these are the computations that don't require 'cb'. 7739 ! first, set to index ranges 7740 7741 i_start = its 7742 i_end = MIN(ite,ide-1) 7743 j_start = jts 7744 j_end = MIN(jte,jde-1) 7745 7746 ! compute x (u) conditions for v, w, or scalar 7747 7748 IF( (config_flags%open_xs) .and. (its == ids) ) THEN 7749 7750 DO j = j_start, j_end 7751 DO k = kts, ktf 7752 ub = MIN( 0.5*(ru(its,k,j)+ru(its+1,k,j)), 0. ) 7753 tendency(its,k,j) = tendency(its,k,j) & 7754 - rdx*( & 7755 ub*( field_old(its+1,k,j) & 7756 - field_old(its ,k,j) ) + & 7757 field(its,k,j)*(ru(its+1,k,j)-ru(its,k,j)) & 7758 ) 7759 ENDDO 7760 ENDDO 7761 7762 ENDIF 7763 7764 IF( (config_flags%open_xe) .and. (ite == ide) ) THEN 7765 7766 DO j = j_start, j_end 7767 DO k = kts, ktf 7768 ub = MAX( 0.5*(ru(ite-1,k,j)+ru(ite,k,j)), 0. ) 7769 tendency(i_end,k,j) = tendency(i_end,k,j) & 7770 - rdx*( & 7771 ub*( field_old(i_end ,k,j) & 7772 - field_old(i_end-1,k,j) ) + & 7773 field(i_end,k,j)*(ru(ite,k,j)-ru(ite-1,k,j)) & 7774 ) 7775 ENDDO 7776 ENDDO 7777 7778 ENDIF 7779 7780 IF( (config_flags%open_ys) .and. (jts == jds) ) THEN 7781 7782 DO i = i_start, i_end 7783 DO k = kts, ktf 7784 vb = MIN( 0.5*(rv(i,k,jts)+rv(i,k,jts+1)), 0. ) 7785 tendency(i,k,jts) = tendency(i,k,jts) & 7786 - rdy*( & 7787 vb*( field_old(i,k,jts+1) & 7788 - field_old(i,k,jts ) ) + & 7789 field(i,k,jts)*(rv(i,k,jts+1)-rv(i,k,jts)) & 7790 ) 7791 ENDDO 7792 ENDDO 7793 7794 ENDIF 7795 7796 IF( (config_flags%open_ye) .and. (jte == jde)) THEN 7797 7798 DO i = i_start, i_end 7799 DO k = kts, ktf 7800 vb = MAX( 0.5*(rv(i,k,jte-1)+rv(i,k,jte)), 0. ) 7801 tendency(i,k,j_end) = tendency(i,k,j_end) & 7802 - rdy*( & 7803 vb*( field_old(i,k,j_end ) & 7804 - field_old(i,k,j_end-1) ) + & 7805 field(i,k,j_end)*(rv(i,k,jte)-rv(i,k,jte-1)) & 7806 ) 7807 ENDDO 7808 ENDDO 7809 7810 ENDIF 7811 7812 !-------------------- vertical advection 7813 7814 !-- loop bounds for periodic or sym conditions 7815 7816 i_start = its-1 7817 i_end = MIN(ite,ide-1)+1 7818 j_start = jts-1 7819 j_end = MIN(jte,jde-1)+1 7820 7821 !-- loop bounds for open or specified conditions 7822 7823 ! WCS 20090218 7824 ! IF(degrade_xs) i_start = its 7825 ! IF(degrade_xe) i_end = MIN(ite,ide-1) 7826 ! IF(degrade_ys) j_start = jts 7827 ! IF(degrade_ye) j_end = MIN(jte,jde-1) 7828 7829 IF(degrade_xs) i_start = MAX(its-1,ids) 7830 IF(degrade_xe) i_end = MIN(ite+1,ide-1) 7831 IF(degrade_ys) j_start = MAX(jts-1,jds) 7832 IF(degrade_ye) j_end = MIN(jte+1,jde-1) 7833 7834 7835 vert_order_test : IF (vert_order == 3) THEN 7836 7837 DO j = j_start, j_end 7838 7839 DO i = i_start, i_end 7840 fqz(i,1,j) = 0. 7841 fqzl(i,1,j) = 0. 7842 fqz(i,kde,j) = 0. 7843 fqzl(i,kde,j) = 0. 7844 ENDDO 7845 7846 DO k=kts+2,ktf-1 7847 DO i = i_start, i_end 7848 7849 vel = rom(i,k,j) 7850 cr = -vel 7851 fqzl(i,k,j) = vel*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr) 7852 7853 fqz(i,k,j) = vel*flux3( & 7854 field(i,k-2,j), field(i,k-1,j), & 7855 field(i,k ,j), field(i,k+1,j), -vel ) 7856 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j) 7857 7858 if(cr.gt. 0) then 7859 qmax(i,k,j) = amax1(qmax(i,k,j),field_old(i,k-1,j)) 7860 qmin(i,k,j) = amin1(qmin(i,k,j),field_old(i,k-1,j)) 7861 else 7862 qmax(i,k-1,j) = amax1(qmax(i,k-1,j),field_old(i,k,j)) 7863 qmin(i,k-1,j) = amin1(qmin(i,k-1,j),field_old(i,k,j)) 7864 end if 7865 7866 ENDDO 7867 ENDDO 7868 7869 DO i = i_start, i_end 7870 7871 k=kts+1 7872 vel = rom(i,k,j) 7873 cr = -vel 7874 fqzl(i,k,j) = vel*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr) 7875 fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j)) 7876 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j) 7877 7878 if(cr.gt. 0) then 7879 qmax(i,k,j) = amax1(qmax(i,k,j),field_old(i,k-1,j)) 7880 qmin(i,k,j) = amin1(qmin(i,k,j),field_old(i,k-1,j)) 7881 else 7882 qmax(i,k-1,j) = amax1(qmax(i,k-1,j),field_old(i,k,j)) 7883 qmin(i,k-1,j) = amin1(qmin(i,k-1,j),field_old(i,k,j)) 7884 end if 7885 7886 k=ktf 7887 vel = rom(i,k,j) 7888 cr = -vel 7889 fqzl(i,k,j) = vel*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr) 7890 fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j)) 7891 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j) 7892 7893 if(cr.gt. 0) then 7894 qmax(i,k,j) = amax1(qmax(i,k,j),field_old(i,k-1,j)) 7895 qmin(i,k,j) = amin1(qmin(i,k,j),field_old(i,k-1,j)) 7896 else 7897 qmax(i,k-1,j) = amax1(qmax(i,k-1,j),field_old(i,k,j)) 7898 qmin(i,k-1,j) = amin1(qmin(i,k-1,j),field_old(i,k,j)) 7899 end if 7900 ENDDO 7901 7902 ENDDO 7903 7904 ELSE 7905 7906 WRITE (wrf_err_message,*) ' advect_scalar_mono, v_order not known ',vert_order 7907 CALL wrf_error_fatal ( wrf_err_message ) 7908 7909 ENDIF vert_order_test 7910 7911 IF (mono_limit) THEN 7912 7913 ! montonic filter 7914 7915 i_start = its-1 7916 i_end = MIN(ite,ide-1)+1 7917 j_start = jts-1 7918 j_end = MIN(jte,jde-1)+1 7919 7920 ! WCS 20090218 7921 7922 !-- loop bounds for open or specified conditions 7923 ! 7924 ! IF(degrade_xs) i_start = its 7925 ! IF(degrade_xe) i_end = MIN(ite,ide-1) 7926 ! IF(degrade_ys) j_start = jts 7927 ! IF(degrade_ye) j_end = MIN(jte,jde-1) 7928 ! 7929 ! IF(config_flags%specified .or. config_flags%nested) THEN 7930 ! IF (degrade_xs) i_start = MAX(its,ids+1) 7931 ! IF (degrade_xe) i_end = MIN(ite,ide-2) 7932 ! IF (degrade_ys) j_start = MAX(jts,jds+1) 7933 ! IF (degrade_ye) j_end = MIN(jte,jde-2) 7934 ! END IF 7935 ! 7936 ! IF(config_flags%open_xs) THEN 7937 ! IF (degrade_xs) i_start = MAX(its,ids+1) 7938 ! END IF 7939 ! IF(config_flags%open_xe) THEN 7940 ! IF (degrade_xe) i_end = MIN(ite,ide-2) 7941 ! END IF 7942 ! IF(config_flags%open_ys) THEN 7943 ! IF (degrade_ys) j_start = MAX(jts,jds+1) 7944 ! END IF 7945 ! IF(config_flags%open_ye) THEN 7946 ! IF (degrade_ye) j_end = MIN(jte,jde-2) 7947 ! END IF 7948 7949 IF(degrade_xs) i_start = MAX(its-1,ids) 7950 IF(degrade_xe) i_end = MIN(ite+1,ide-1) 7951 IF(degrade_ys) j_start = MAX(jts-1,jds) 7952 IF(degrade_ye) j_end = MIN(jte+1,jde-1) 7953 7954 IF(config_flags%specified .or. config_flags%nested) THEN 7955 IF (degrade_xs) i_start = MAX(its-1,ids+1) 7956 IF (degrade_xe) i_end = MIN(ite+1,ide-2) 7957 IF (degrade_ys) j_start = MAX(jts-1,jds+1) 7958 IF (degrade_ye) j_end = MIN(jte+1,jde-2) 7959 END IF 7960 7961 IF(config_flags%open_xs) THEN 7962 IF (degrade_xs) i_start = MAX(its-1,ids+1) 7963 END IF 7964 IF(config_flags%open_xe) THEN 7965 IF (degrade_xe) i_end = MIN(ite+1,ide-2) 7966 END IF 7967 IF(config_flags%open_ys) THEN 7968 IF (degrade_ys) j_start = MAX(jts-1,jds+1) 7969 END IF 7970 IF(config_flags%open_ye) THEN 7971 IF (degrade_ye) j_end = MIN(jte+1,jde-2) 7972 END IF 7973 7974 !-- here is the limiter... 7975 7976 DO j=j_start, j_end 7977 DO k=kts, ktf 7978 DO i=i_start, i_end 7979 7980 ph_upwind = (mub(i,j)+mu_old(i,j))*field_old(i,k,j) & 7981 - dt*( msft(i,j) *( & 7982 rdx*(fqxl(i+1,k,j)-fqxl(i,k,j)) + & 7983 rdy*(fqyl(i,k,j+1)-fqyl(i,k,j)) ) & 7984 + rdzw(k)*(fqzl(i,k+1,j)-fqzl(i,k,j)) ) 7985 7986 flux_in = -dt*( msft(i,j) *( & 7987 rdx*( min(0.,fqx (i+1,k,j)) & 7988 -max(0.,fqx (i ,k,j)) ) & 7989 +rdy*( min(0.,fqy (i,k,j+1)) & 7990 -max(0.,fqy (i,k,j )) ) ) & 7991 + rdzw(k)*( max(0.,fqz (i,k+1,j)) & 7992 -min(0.,fqz (i,k ,j)) ) ) 7993 7994 ph_hi = mut(i,j)*qmax(i,k,j) - ph_upwind 7995 IF( flux_in .gt. ph_hi ) scale_in(i,k,j) = max(0.,ph_hi/(flux_in+eps)) 7996 7997 7998 flux_out = dt*( msft(i,j)*( & 7999 rdx*( max(0.,fqx (i+1,k,j)) & 8000 -min(0.,fqx (i ,k,j)) ) & 8001 +rdy*( max(0.,fqy (i,k,j+1)) & 8002 -min(0.,fqy (i,k,j )) ) ) & 8003 + rdzw(k)*( min(0.,fqz (i,k+1,j)) & 8004 -max(0.,fqz (i,k ,j)) ) ) 8005 8006 ph_low = ph_upwind - mut(i,j)*qmin(i,k,j) 8007 IF( flux_out .gt. ph_low ) scale_out(i,k,j) = max(0.,ph_low/(flux_out+eps)) 8008 8009 ENDDO 8010 ENDDO 8011 ENDDO 8012 8013 DO j=j_start, j_end 8014 DO k=kts, ktf 8015 DO i=i_start, i_end+1 8016 IF( fqx (i,k,j) .gt. 0.) then 8017 fqx(i,k,j) = min(scale_in(i,k,j),scale_out(i-1,k,j))*fqx(i,k,j) 8018 ELSE 8019 fqx(i,k,j) = min(scale_out(i,k,j),scale_in(i-1,k,j))*fqx(i,k,j) 8020 ENDIF 8021 ENDDO 8022 ENDDO 8023 ENDDO 8024 8025 DO j=j_start, j_end+1 8026 DO k=kts, ktf 8027 DO i=i_start, i_end 8028 IF( fqy (i,k,j) .gt. 0.) then 8029 fqy(i,k,j) = min(scale_in(i,k,j),scale_out(i,k,j-1))*fqy(i,k,j) 8030 ELSE 8031 fqy(i,k,j) = min(scale_out(i,k,j),scale_in(i,k,j-1))*fqy(i,k,j) 8032 ENDIF 8033 ENDDO 8034 ENDDO 8035 ENDDO 8036 8037 DO j=j_start, j_end 8038 DO k=kts+1, ktf 8039 DO i=i_start, i_end 8040 IF( fqz (i,k,j) .lt. 0.) then 8041 fqz(i,k,j) = min(scale_in(i,k,j),scale_out(i,k-1,j))*fqz(i,k,j) 8042 ELSE 8043 fqz(i,k,j) = min(scale_out(i,k,j),scale_in(i,k-1,j))*fqz(i,k,j) 8044 ENDIF 8045 ENDDO 8046 ENDDO 8047 ENDDO 8048 8049 END IF 8050 8051 ! add in the mono-limited flux divergence 8052 ! we need to fix this for open b.c set *********** 8053 8054 i_start = its 8055 i_end = MIN(ite,ide-1) 8056 j_start = jts 8057 j_end = MIN(jte,jde-1) 8058 8059 DO j = j_start, j_end 8060 DO k = kts, ktf 8061 DO i = i_start, i_end 8062 8063 tendency (i,k,j) = tendency(i,k,j) & 8064 -rdzw(k)*( fqz (i,k+1,j)-fqz (i,k,j) & 8065 +fqzl(i,k+1,j)-fqzl(i,k,j)) 8066 8067 ENDDO 8068 ENDDO 8069 ENDDO 8070 8071 ! x flux divergence 8072 ! 8073 8074 ! WCS 20090218 8075 ! IF(degrade_xs) i_start = i_start + 1 8076 ! IF(degrade_xe) i_end = i_end - 1 8077 8078 IF(degrade_xs) i_start = MAX(its,ids+1) 8079 IF(degrade_xe) i_end = MIN(ite,ide-2) 8080 8081 DO j = j_start, j_end 8082 DO k = kts, ktf 8083 DO i = i_start, i_end 8084 8085 ! Un-"canceled" map scale factor, ADT Eq. 48 8086 tendency (i,k,j) = tendency(i,k,j) & 8087 - msft(i,j)*( rdx*( fqx (i+1,k,j)-fqx (i,k,j) & 8088 +fqxl(i+1,k,j)-fqxl(i,k,j)) ) 8089 8090 ENDDO 8091 ENDDO 8092 ENDDO 8093 8094 ! y flux divergence 8095 ! 8096 i_start = its 8097 i_end = MIN(ite,ide-1) 8098 8099 ! WCS 20090218 8100 ! IF(degrade_ys) j_start = j_start + 1 8101 ! IF(degrade_ye) j_end = j_end - 1 8102 8103 IF(degrade_ys) j_start = MAX(jts,jds+1) 8104 IF(degrade_ye) j_end = MIN(jte,jde-2) 8105 8106 DO j = j_start, j_end 8107 DO k = kts, ktf 8108 DO i = i_start, i_end 8109 8110 ! Un-"canceled" map scale factor, ADT Eq. 48 8111 tendency (i,k,j) = tendency(i,k,j) & 8112 - msft(i,j)*( rdy*( fqy (i,k,j+1)-fqy (i,k,j) & 8113 +fqyl(i,k,j+1)-fqyl(i,k,j)) ) 8114 8115 ENDDO 8116 ENDDO 8117 ENDDO 8118 8119 END SUBROUTINE advect_scalar_mono 8120 8121 !----------------------------------------------------------- 8122 7055 8123 END MODULE module_advect_em 7056 8124 -
trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_em.F
r11 r196 16 16 u, v, w, t, ph, mu, & 17 17 moist, & 18 ru, rv, rw, ww, php, alt, muu, muv, & 18 ru, rv, rw, ww, php, alt, & 19 muu, muv, & 19 20 mub, mut, phb, pb, p, al, alb, & 20 21 cqu, cqv, cqw, & … … 280 281 REAL :: kdift, khdq, kvdq, cfn, cfn1, cf1, cf2, cf3 281 282 INTEGER :: i,j,k 283 INTEGER :: time_step 282 284 283 285 !<DESCRIPTION> … … 362 364 363 365 ! advection tendencies 366 CALL nl_get_time_step ( 1, time_step ) 364 367 365 368 CALL advect_u ( u, u , ru_tend, ru, rv, ww, & 366 mut, config_flags,&369 mut, time_step, config_flags, & 367 370 msfu, msfv, msft, & 368 371 fnm, fnp, rdx, rdy, rdnw, & … … 372 375 373 376 CALL advect_v ( v, v , rv_tend, ru, rv, ww, & 374 mut, config_flags,&377 mut, time_step, config_flags, & 375 378 msfu, msfv, msft, & 376 379 fnm, fnp, rdx, rdy, rdnw, & … … 381 384 IF (non_hydrostatic) & 382 385 CALL advect_w ( w, w, rw_tend, ru, rv, ww, & 383 mut, config_flags,&386 mut, time_step, config_flags, & 384 387 msfu, msfv, msft, & 385 388 fnm, fnp, rdx, rdy, rdn, & … … 391 394 392 395 CALL advect_scalar ( t, t, t_tend, ru, rv, ww, & 393 mut, config_flags,&396 mut, time_step, config_flags, & 394 397 msfu, msfv, msft, fnm, fnp, & 395 398 rdx, rdy, rdnw, & … … 418 421 ims, ime, jms, jme, kms, kme, & 419 422 its, ite, jts, jte, kts, kte ) 420 421 423 422 424 CALL horizontal_pressure_gradient( ru_tend,rv_tend, & … … 454 456 its, ite, jts, jte, kts, kte ) 455 457 ELSE 456 457 458 CALL coriolis ( ru, rv, rw, & 458 459 ru_tend, rv_tend, rw_tend, & … … 463 464 its, ite, jts, jte, kts, kte ) 464 465 465 END IF 466 466 END IF 467 467 468 468 CALL curvature ( ru, rv, rw, u, v, w, & … … 515 515 its, ite, jts, jte, kts, kte ) 516 516 517 517 518 !!!****MARS: vertical diffusion is done in the physics (TODO: consider the nonhydrostatic case ?) 518 519 pbl_test : IF ( (config_flags%bl_pbl_physics .eq. 0) & … … 548 549 549 550 ENDIF pbl_test 550 551 551 552 552 ! Theta tendency computations. … … 654 654 INTEGER :: i, j, k 655 655 656 657 656 !<DESCRIPTION> 658 657 ! … … 679 678 DO k = kts,kte-1 680 679 DO i = its,ite 680 ! multiply by m to uncouple u 681 681 IF(rk_step == 1)ru_tendf(i,k,j) = ru_tendf(i,k,j) + u_save(i,k,j)*msfu(i,j) 682 ! divide by m to couple u 682 683 ru_tend(i,k,j) = ru_tend(i,k,j) + ru_tendf(i,k,j)/msfu(i,j) 683 684 ENDDO … … 688 689 DO k = kts,kte-1 689 690 DO i = its,MIN(ite,ide-1) 691 ! multiply by m to uncouple v 690 692 IF(rk_step == 1)rv_tendf(i,k,j) = rv_tendf(i,k,j) + v_save(i,k,j)*msfv(i,j) 693 ! divide by m to couple v 691 694 rv_tend(i,k,j) = rv_tend(i,k,j) + rv_tendf(i,k,j)/msfv(i,j) 692 695 ENDDO … … 697 700 DO k = kts,kte 698 701 DO i = its,MIN(ite,ide-1) 702 ! multiply by m to uncouple w 699 703 IF(rk_step == 1)rw_tendf(i,k,j) = rw_tendf(i,k,j) + w_save(i,k,j)*msft(i,j) 704 ! divide by m to couple w 700 705 rw_tend(i,k,j) = rw_tend(i,k,j) + rw_tendf(i,k,j)/msft(i,j) 701 706 IF(rk_step == 1)ph_tendf(i,k,j) = ph_tendf(i,k,j) + ph_save(i,k,j) 707 ! divide by m to couple scalar 702 708 ph_tend(i,k,j) = ph_tend(i,k,j) + ph_tendf(i,k,j)/msft(i,j) 703 709 ENDDO … … 709 715 DO i = its,MIN(ite,ide-1) 710 716 IF(rk_step == 1)t_tendf(i,k,j) = t_tendf(i,k,j) + t_save(i,k,j) 711 t_tend(i,k,j) = t_tend(i,k,j) + t_tendf(i,k,j)/msft(i,j) & 717 ! divide by m to couple theta 718 t_tend(i,k,j) = t_tend(i,k,j) + t_tendf(i,k,j)/msft(i,j) & 712 719 + mut(i,j)*h_diabatic(i,k,j)/msft(i,j) 720 ! divide by m to couple heating 713 721 ENDDO 714 722 ENDDO … … 757 765 758 766 REAL, DIMENSION(ims:ime, kms:kme, jms:jme , scs:sce ), & 759 INTENT(IN OUT) :: scalar, scalar_old767 INTENT(IN ) :: scalar, scalar_old 760 768 761 769 REAL, DIMENSION(ims:ime, kms:kme, jms:jme , scs:sce ), & 762 770 INTENT(INOUT) :: scalar_tends 763 771 764 REAL, DIMENSION(ims:ime, kms:kme, jms:jme ) :: advect_tend772 REAL, DIMENSION(ims:ime, kms:kme, jms:jme ), INTENT(INOUT) :: advect_tend 765 773 766 774 REAL, DIMENSION(ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: RQVFTEN … … 801 809 802 810 INTEGER :: im, i,j,k 811 INTEGER :: time_step 803 812 804 813 REAL :: khdq, kvdq, tendency … … 822 831 its, ite, jts, jte, kts, kte ) 823 832 833 CALL nl_get_time_step ( 1, time_step ) 834 824 835 IF( (rk_step == 3) .and. pd_advection ) THEN 825 836 826 CALL advect_scalar_pd ( scalar(ims,kms,jms,im), & 827 scalar_old(ims,kms,jms,im), & 828 advect_tend(ims,kms,jms), & 829 ru, rv, ww, mut, mub, mu_old, & 830 config_flags, & 831 msfu, msfv, msft, fnm, fnp, & 832 rdx, rdy, rdnw,dt, & 833 ids, ide, jds, jde, kds, kde, & 834 ims, ime, jms, jme, kms, kme, & 835 its, ite, jts, jte, kts, kte ) 836 837 CALL advect_scalar_pd ( scalar(ims,kms,jms,im), & 838 scalar_old(ims,kms,jms,im), & 839 advect_tend(ims,kms,jms), & 840 ru, rv, ww, mut, mub, mu_old, & 841 config_flags, & 842 msfu, msfv, msft, fnm, fnp, & 843 rdx, rdy, rdnw,dt, & 844 ids, ide, jds, jde, kds, kde, & 845 ims, ime, jms, jme, kms, kme, & 846 its, ite, jts, jte, kts, kte ) 847 848 849 850 CALL advect_scalar_mono ( scalar(ims,kms,jms,im), & 851 scalar_old(ims,kms,jms,im), & 852 advect_tend(ims,kms,jms), & 853 ru, rv, ww, mut, mub, mu_old, & 854 config_flags, & 855 msfu, msfv, msft, fnm, fnp, & 856 rdx, rdy, rdnw,dt, & 857 ids, ide, jds, jde, kds, kde, & 858 ims, ime, jms, jme, kms, kme, & 859 its, ite, jts, jte, kts, kte ) 837 860 ELSE 838 861 839 CALL advect_scalar ( scalar(ims,kms,jms,im), & 840 scalar(ims,kms,jms,im), & 841 advect_tend(ims,kms,jms), & 842 ru, rv, ww, mut, config_flags, & 843 msfu, msfv, msft, fnm, fnp, & 844 rdx, rdy, rdnw, & 845 ids, ide, jds, jde, kds, kde, & 846 ims, ime, jms, jme, kms, kme, & 847 its, ite, jts, jte, kts, kte ) 862 CALL advect_scalar ( scalar(ims,kms,jms,im), & 863 scalar(ims,kms,jms,im), & 864 advect_tend(ims,kms,jms), & 865 ru, rv, ww, mut, time_step, & 866 config_flags, & 867 msfu, msfv, msft, fnm, fnp, & 868 rdx, rdy, rdnw, & 869 ids, ide, jds, jde, kds, kde, & 870 ims, ime, jms, jme, kms, kme, & 871 its, ite, jts, jte, kts, kte ) 848 872 END IF 849 873 … … 863 887 scalar_tends(ims,kms,jms,im), mut, & 864 888 config_flags, & 865 msfu, msfv, msft, khdq , xkmhd, rdx, rdy, & 889 msfu, msfv, msft, & 890 khdq , xkmhd, rdx, rdy, & 866 891 ids, ide, jds, jde, kds, kde, & 867 892 ims, ime, jms, jme, kms, kme, & … … 940 965 REAL, DIMENSION(ims:ime, kms:kme, jms:jme , scs:sce), & 941 966 INTENT(INOUT) :: scalar_1, & 942 scalar_2, & 943 sc_tend 967 scalar_2 968 969 REAL, DIMENSION(ims:ime, kms:kme, jms:jme , scs:sce), & 970 INTENT(IN) :: sc_tend 944 971 945 972 REAL, DIMENSION(ims:ime, kms:kme, jms:jme ), & … … 1014 1041 DO k = k_start,k_end 1015 1042 DO i = i_start,i_end 1043 ! scalar was coupled with m 1016 1044 tendency(i,k,j) = advect_tend(i,k,j) * msft(i,j) 1017 1045 ENDDO … … 1064 1092 DO k = k_start,k_end 1065 1093 DO i = i_start,i_end 1094 ! scalar was coupled with m 1066 1095 tendency(i,k,j) = advect_tend(i,k,j) * msft(i,j) 1067 1096 ENDDO … … 1104 1133 SUBROUTINE rk_update_scalar_pd( scs, sce, & 1105 1134 scalar, sc_tend, & 1106 msft, &1107 1135 mu_old, mu_new, mu_base, & 1108 1136 rk_step, dt, spec_zone, & … … 1131 1159 REAL, DIMENSION(ims:ime, jms:jme ), INTENT(IN ) :: mu_old, & 1132 1160 mu_new, & 1133 mu_base, & 1134 msft 1161 mu_base 1135 1162 1136 1163 INTEGER :: i,j,k,im … … 1523 1550 !!****MARS 1524 1551 IF ( (config_flags%bl_pbl_physics .gt. 0) .OR. (config_flags%modif_wrf) ) THEN 1552 1525 1553 DO J=jts,jtf 1526 1554 DO K=kts,ktf … … 1564 1592 1565 1593 ENDIF 1566 1567 1568 1594 1569 1595 ! fdda -
trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/dyn_em/solve_em.F
r156 r196 2046 2046 moist_old(ims,kms,jms,im), & 2047 2047 moist_tend(ims,kms,jms,im), & 2048 grid%msft, &2049 2048 grid%em_mu_1, grid%em_mu_1, grid%em_mub, & 2050 2049 rk_step, dt_rk, grid%spec_zone, & … … 2111 2110 scalar_old(ims,kms,jms,im), & 2112 2111 scalar_tend(ims,kms,jms,im), & 2113 grid%msft, &2114 2112 grid%em_mu_1, grid%em_mu_1, grid%em_mub, & 2115 2113 rk_step, dt_rk, grid%spec_zone, & … … 2178 2176 chem_old(ims,kms,jms,im), & 2179 2177 chem_tend(ims,kms,jms,im), & 2180 grid%msft, &2181 2178 grid%em_mu_1, grid%em_mu_1, grid%em_mub, & 2182 2179 rk_step, dt_rk, grid%spec_zone, & … … 2246 2243 grid%em_tke_1, & 2247 2244 tke_tend(ims,kms,jms), & 2248 grid%msft, &2249 2245 grid%em_mu_1, grid%em_mu_1, grid%em_mub, & 2250 2246 rk_step, dt_rk, grid%spec_zone, &
Note: See TracChangeset
for help on using the changeset viewer.
