Changeset 713
- Timestamp:
- Mar 30, 2011 2:21:21 PM (14 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 2 edited
-
advec_ws.f90 (modified) (100 diffs)
-
init_3d_model.f90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/advec_ws.f90
r711 r713 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! File reformatted. 7 ! Bugfix in vertical advection of w concerning the optimized version for 8 ! vector architecture. 9 ! Constants adv_mom_3, adv_mom_5, adv_sca_5, adv_sca_3 reformulated as 10 ! broken numbers. 7 11 ! 8 12 ! Former revisions: … … 13 17 ! formatting adjustments 14 18 ! 15 ! 705 2011-03-25 11:21:43 Z suehring 19 ! 705 2011-03-25 11:21:43 Z suehring $ 16 20 ! Bugfix in declaration of logicals concerning outflow boundaries. 17 21 ! … … 111 115 ! 112 116 !-- Set the appropriate factors for scalar and momentum advection. 113 adv_sca_5 = 0.016666666666666114 adv_sca_3 = 0.083333333333333115 adv_mom_5 = 0.0083333333333333116 adv_mom_3 = 0.041666666666666117 adv_sca_5 = 1./60. 118 adv_sca_3 = 1./12. 119 adv_mom_5 = 1./120. 120 adv_mom_3 = 1./24. 117 121 118 122 ! … … 317 321 CASE ( 3 ) 318 322 319 DO k = nzb_s_inner(j,i) +1, nzt323 DO k = nzb_s_inner(j,i)+1, nzt 320 324 v_comp = v(k,j+1,i) - v_gtrans 321 325 flux_n(k) = v_comp * ( & … … 473 477 IF ( j == nys .AND. ( .NOT. degraded_s ) ) THEN 474 478 475 DO k = nzb_s_inner(j,i) +1, nzt479 DO k = nzb_s_inner(j,i)+1, nzt 476 480 v_comp = v(k,j,i) - v_gtrans 477 481 swap_flux_y_local(k) = v_comp * ( & … … 534 538 ! 535 539 !-- 2nd-order scheme (bottom) 536 k = nzb_s_inner(j,i) +1540 k = nzb_s_inner(j,i)+1 537 541 flux_d = flux_t(k-1) 538 542 diss_d = diss_t(k-1) … … 684 688 685 689 CASE ( 1 ) 686 DO k = nzb_u_inner(j,i) +1, nzt690 DO k = nzb_u_inner(j,i)+1, nzt 687 691 u_comp(k) = u(k,j,i+1) + u(k,j,i) 688 692 flux_r(k) = ( u_comp(k) - gu ) * ( & 689 7. * (u(k,j,i+1) + u(k,j,i) )&690 -( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3691 diss_r(k) = - abs(u_comp(k) - gu) * (&692 3. * (u(k,j,i+1) - u(k,j,i) )&693 -( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3693 7.0 * ( u(k,j,i+1) + u(k,j,i) ) & 694 - ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3 695 diss_r(k) = - ABS( u_comp(k) - gu ) * ( & 696 3.0 * ( u(k,j,i+1) - u(k,j,i) ) & 697 - ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3 694 698 ENDDO 695 699 696 700 CASE ( 2 ) 697 DO k = nzb_u_inner(j,i) +1, nzt701 DO k = nzb_u_inner(j,i)+1, nzt 698 702 u_comp(k) = u(k,j,i+1) + u(k,j,i) 699 flux_r(k) = (u_comp(k) - gu) * ( u(k,j,i+1) + u(k,j,i) )*0.25 703 flux_r(k) = ( u_comp(k) - gu ) * ( & 704 u(k,j,i+1) + u(k,j,i) ) * 0.25 700 705 diss_r(k) = diss_2nd( u(k,j,i+1) ,u(k,j,i+1), u(k,j,i), & 701 706 u(k,j,i-1), u(k,j,i-2), u_comp(k), & … … 704 709 705 710 CASE ( 3 ) 706 DO k = nzb_u_inner(j,i) +1, nzt711 DO k = nzb_u_inner(j,i)+1, nzt 707 712 v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv 708 713 flux_n(k) = v_comp * ( & 709 7. * ( u(k,j+1,i) + u(k,j,i) )&710 -( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3711 diss_n(k) = - abs(v_comp) * (&712 3. * ( u(k,j+1,i) - u(k,j,i) )&713 -( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3714 7.0 * ( u(k,j+1,i) + u(k,j,i) ) & 715 - ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3 716 diss_n(k) = - ABS( v_comp ) * ( & 717 3.0 * ( u(k,j+1,i) - u(k,j,i) ) & 718 - ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3 714 719 ENDDO 715 720 716 721 CASE ( 4 ) 717 DO k = nzb_u_inner(j,i) +1, nzt722 DO k = nzb_u_inner(j,i)+1, nzt 718 723 v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv 719 724 flux_n(k) = v_comp * ( u(k,j+1,i) + u(k,j,i) ) * 0.25 … … 724 729 725 730 CASE ( 5 ) 726 DO k = nzb_u_inner(j,i) +1, nzt731 DO k = nzb_u_inner(j,i)+1, nzt 727 732 ! 728 733 !-- Compute leftside fluxes for the left boundary of PE domain 729 734 u_comp(k) = u(k,j,i) + u(k,j,i-1) - gu 730 735 flux_l_u(k,j) = u_comp(k) * ( u(k,j,i) + u(k,j,i-1) ) * 0.25 731 diss_l_u(k,j) = diss_2nd( u(k,j,i+2), u(k,j,i+1), u(k,j,i),&732 u(k,j,i-1), u(k,j,i-1), u_comp(k),&733 0.25, ddx )736 diss_l_u(k,j) = diss_2nd( u(k,j,i+2), u(k,j,i+1), u(k,j,i), & 737 u(k,j,i-1), u(k,j,i-1), u_comp(k),& 738 0.25, ddx ) 734 739 735 740 u_comp(k) = u(k,j,i+1) + u(k,j,i) 736 741 flux_r(k) = ( u_comp(k) - gu ) * ( & 737 7. * (u(k,j,i+1) + u(k,j,i) )&738 -( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3739 diss_r(k) = - abs( u_comp(k) -gu ) * ( &740 3. * ( u(k,j,i+1) - u(k,j,i) )&741 -( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3742 7.0 * ( u(k,j,i+1) + u(k,j,i) ) & 743 - ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3 744 diss_r(k) = - ABS( u_comp(k) -gu ) * ( & 745 3.0 * ( u(k,j,i+1) - u(k,j,i) ) & 746 - ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3 742 747 ENDDO 743 748 degraded_l = .TRUE. 744 749 745 750 CASE ( 7 ) 746 DO k = nzb_u_inner(j,i) +1, nzt751 DO k = nzb_u_inner(j,i)+1, nzt 747 752 v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv 748 753 flux_n(k) = v_comp * ( & 749 7. * ( u(k,j+1,i) + u(k,j,i) )&750 -( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3751 diss_n(k) = - abs(v_comp) * (&752 3. * ( u(k,j+1,i) - u(k,j,i) )&753 -( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3754 7.0 * ( u(k,j+1,i) + u(k,j,i) ) & 755 - ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3 756 diss_n(k) = - ABS( v_comp ) * ( & 757 3.0 * ( u(k,j+1,i) - u(k,j,i) ) & 758 - ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3 754 759 ENDDO 755 760 756 761 CASE ( 8 ) 757 DO k = nzb_u_inner(j,i) +1, nzt762 DO k = nzb_u_inner(j,i)+1, nzt 758 763 ! 759 764 !-- Compute southside fluxes for the south boundary of PE domain … … 766 771 v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv 767 772 flux_n(k) = v_comp * ( & 768 7. * ( u(k,j+1,i) + u(k,j,i) )&769 -( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3770 diss_n(k) = - abs(v_comp) * (&771 3. * ( u(k,j+1,i) - u(k,j,i) )&772 -( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3773 7.0 * ( u(k,j+1,i) + u(k,j,i) ) & 774 - ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3 775 diss_n(k) = - ABS( v_comp ) * ( & 776 3.0 * ( u(k,j+1,i) - u(k,j,i) ) & 777 - ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3 773 778 ENDDO 774 779 degraded_s = .TRUE. … … 779 784 ! 780 785 !-- Compute the crosswise 5th order fluxes at the outflow 781 IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 782 .OR.boundary_flags(j,i) == 5 ) THEN786 IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 .OR. & 787 boundary_flags(j,i) == 5 ) THEN 783 788 784 789 DO k = nzb_u_inner(j,i)+1, nzt 785 790 v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv 786 791 flux_n(k) = v_comp * ( & 787 37. 788 - 8.* ( u(k,j+2,i) + u(k,j-1,i) ) &789 + ( u(k,j+3,i) + u(k,j-2,i) ) ) *adv_mom_5790 diss_n(k) = - abs(v_comp) * (&791 10. 792 - 5.* ( u(k,j+2,i) - u(k,j-1,i) ) &793 +( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5792 37.0 * ( u(k,j+1,i) + u(k,j,i) ) & 793 - 8.0 * ( u(k,j+2,i) + u(k,j-1,i) ) & 794 + ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5 795 diss_n(k) = - ABS ( v_comp ) * ( & 796 10.0 * ( u(k,j+1,i) - u(k,j,i) ) & 797 - 5.0 * ( u(k,j+2,i) - u(k,j-1,i) ) & 798 + ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5 794 799 ENDDO 795 800 796 801 ELSE 797 802 798 DO k = nzb_u_inner(j,i) +1, nzt803 DO k = nzb_u_inner(j,i)+1, nzt 799 804 u_comp(k) = u(k,j,i+1) + u(k,j,i) 800 805 flux_r(k) = ( u_comp(k) - gu ) * ( & 801 37. 802 - 8.* ( u(k,j,i+2) + u(k,j,i-1) ) &803 +( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5804 diss_r(k) = - abs(u_comp(k) - gu) * (&805 10. 806 - 5.* ( u(k,j,i+2) - u(k,j,i-1) ) &807 + ( u(k,j,i+3) - u(k,j,i-2) ) ) *adv_mom_5806 37.0 * ( u(k,j,i+1) + u(k,j,i) ) & 807 - 8.0 * ( u(k,j,i+2) + u(k,j,i-1) ) & 808 + ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5 809 diss_r(k) = - ABS( u_comp(k) - gu ) * ( & 810 10.0 * ( u(k,j,i+1) - u(k,j,i) ) & 811 - 5.0 * ( u(k,j,i+2) - u(k,j,i-1) ) & 812 + ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5 808 813 ENDDO 809 814 … … 813 818 ! 814 819 !-- Compute the fifth order fluxes for the interior of PE domain. 815 DO k = nzb_u_inner(j,i) +1, nzt820 DO k = nzb_u_inner(j,i)+1, nzt 816 821 u_comp(k) = u(k,j,i+1) + u(k,j,i) 817 822 flux_r(k) = ( u_comp(k) - gu ) * ( & 818 37. 819 - 8.* ( u(k,j,i+2) + u(k,j,i-1) ) &820 +( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5821 diss_r(k) = - abs(u_comp(k) - gu) * (&822 10. 823 - 5.* ( u(k,j,i+2) - u(k,j,i-1) ) &824 +( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5823 37.0 * ( u(k,j,i+1) + u(k,j,i) ) & 824 - 8.0 * ( u(k,j,i+2) + u(k,j,i-1) ) & 825 + ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5 826 diss_r(k) = - ABS( u_comp(k) - gu ) * ( & 827 10.0 * ( u(k,j,i+1) - u(k,j,i) ) & 828 - 5.0 * ( u(k,j,i+2) - u(k,j,i-1) ) & 829 + ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5 825 830 826 831 v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv 827 832 flux_n(k) = v_comp * ( & 828 37. 829 - 8.* ( u(k,j+2,i) + u(k,j-1,i) ) &830 +( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5831 diss_n(k) = - abs(v_comp) * (&832 10. 833 - 5.* ( u(k,j+2,i) - u(k,j-1,i) ) &834 +( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5833 37.0 * ( u(k,j+1,i) + u(k,j,i) ) & 834 - 8.0 * ( u(k,j+2,i) + u(k,j-1,i) ) & 835 + ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5 836 diss_n(k) = - ABS( v_comp ) * ( & 837 10.0 * ( u(k,j+1,i) - u(k,j,i) ) & 838 - 5.0 * ( u(k,j+2,i) - u(k,j-1,i) ) & 839 + ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5 835 840 ENDDO 836 841 … … 840 845 IF ( j == nys .AND. ( .NOT. degraded_s ) ) THEN 841 846 842 DO k = nzb_u_inner(j,i) +1, nzt847 DO k = nzb_u_inner(j,i)+1, nzt 843 848 v_comp = v(k,j,i) + v(k,j,i-1) - gv 844 849 flux_s_u(k) = v_comp * ( & 845 37. 846 - 8.* ( u(k,j+1,i) + u(k,j-2,i) ) &847 +( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5848 diss_s_u(k) = - abs(v_comp) * ( &849 10. * ( u(k,j,i) - u(k,j-1,i) )&850 - 5. * ( u(k,j+1,i) - u(k,j-2,i) )&851 + (u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5850 37.0 * ( u(k,j,i) + u(k,j-1,i) ) & 851 - 8.0 * ( u(k,j+1,i) + u(k,j-2,i) ) & 852 + ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5 853 diss_s_u(k) = - ABS(v_comp) * ( & 854 10.0 * ( u(k,j,i) - u(k,j-1,i) ) & 855 - 5.0 * ( u(k,j+1,i) - u(k,j-2,i) ) & 856 + ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5 852 857 ENDDO 853 858 … … 859 864 u_comp_l = u(k,j,i)+u(k,j,i-1)-gu 860 865 flux_l_u(k,j) = u_comp_l * ( & 861 37. 862 - 8.* ( u(k,j,i+1) + u(k,j,i-2) ) &863 +( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5864 diss_l_u(k,j) = - abs(u_comp_l) * ( &865 10. 866 - 5.* ( u(k,j,i+1) - u(k,j,i-2) ) &867 +( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5866 37.0 * ( u(k,j,i) + u(k,j,i-1) ) & 867 - 8.0 * ( u(k,j,i+1) + u(k,j,i-2) ) & 868 + ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5 869 diss_l_u(k,j) = - ABS(u_comp_l) * ( & 870 10.0 * ( u(k,j,i) - u(k,j,i-1) ) & 871 - 5.0 * ( u(k,j,i+1) - u(k,j,i-2) ) & 872 + ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5 868 873 ENDDO 869 874 … … 871 876 ! 872 877 !-- Now compute the tendency terms for the horizontal parts. 873 DO k = nzb_u_inner(j,i) +1, nzt878 DO k = nzb_u_inner(j,i)+1, nzt 874 879 tend(k,j,i) = tend(k,j,i) - ( & 875 880 ( flux_r(k) + diss_r(k) & 876 - ( flux_l_u(k,j) + diss_l_u(k,j)) ) * ddx &881 - flux_l_u(k,j) - diss_l_u(k,j) ) * ddx & 877 882 + ( flux_n(k) + diss_n(k) & 878 - ( flux_s_u(k) + diss_s_u(k)) ) * ddy )883 - flux_s_u(k) - diss_s_u(k) ) * ddy ) 879 884 880 885 flux_l_u(k,j) = flux_r(k) … … 887 892 sums_us2_ws_l(k,:) = sums_us2_ws_l(k,:) & 888 893 + ( flux_r(k) * & 889 ( u_comp(k) - 2. * hom(k,1,1,:) ) / ( u_comp(k) - gu + 1.0E-20 ) & 890 + diss_r(k) & 891 * abs(u_comp(k) - 2. * hom(k,1,1,:) ) & 892 / (abs(u_comp(k) - gu) + 1.0E-20) ) & 893 * weight_substep(intermediate_timestep_count) * rmask(j,i,:) 894 ( u_comp(k) - 2.0 * hom(k,1,1,:) ) & 895 / ( u_comp(k) - gu + 1.0E-20 ) & 896 + diss_r(k) * & 897 ABS( u_comp(k) - 2.0 * hom(k,1,1,:) ) & 898 / ( ABS( u_comp(k) - gu ) + 1.0E-20 ) ) & 899 * weight_substep(intermediate_timestep_count) * rmask(j,i,:) 894 900 ENDDO 895 sums_us2_ws_l(nzb_u_inner(j,i),:) = &896 sums_us2_ws_l(nzb_u_inner(j,i)+1,:)901 sums_us2_ws_l(nzb_u_inner(j,i),:) = sums_us2_ws_l(nzb_u_inner(j,i)+1,:) 902 897 903 898 904 ! … … 903 909 diss_t(nzb_u_inner(j,i)) = 0. 904 910 ! 905 !-- 2nd order scheme 906 k = nzb_u_inner(j,i) +1911 !-- 2nd order scheme (bottom) 912 k = nzb_u_inner(j,i)+1 907 913 flux_d = flux_t(k-1) 908 914 diss_d = diss_t(k-1) … … 913 919 914 920 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 915 - ( flux_d + diss_d )) * ddzw(k)916 ! 917 !-- WS3 as an intermediate step 918 k = nzb_u_inner(j,i) +2921 - flux_d - diss_d ) * ddzw(k) 922 ! 923 !-- WS3 as an intermediate step (bottom) 924 k = nzb_u_inner(j,i)+2 919 925 flux_d = flux_t(k-1) 920 926 diss_d = diss_t(k-1) 921 927 w_comp = w(k,j,i) + w(k,j,i-1) 922 flux_t(k) = w_comp *(&923 7. * ( u(k+1,j,i) + u(k,j,i) )&924 -( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3925 diss_t(k) = - abs(w_comp)*(&926 3. * ( u(k+1,j,i) - u(k,j,i) )&927 -( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3928 flux_t(k) = w_comp * ( & 929 7.0 * ( u(k+1,j,i) + u(k,j,i) ) & 930 - ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3 931 diss_t(k) = -ABS( w_comp) * ( & 932 3.0 * ( u(k+1,j,i) - u(k,j,i) ) & 933 - ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3 928 934 929 935 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 930 - ( flux_d + diss_d ) ) * ddzw(k) 931 932 DO k = nzb_u_inner(j,i) + 3, nzt - 2 936 - flux_d - diss_d ) * ddzw(k) 937 ! 938 !-- WS5 939 DO k = nzb_u_inner(j,i)+3, nzt-2 933 940 flux_d = flux_t(k-1) 934 941 diss_d = diss_t(k-1) 935 942 w_comp = w(k,j,i) + w(k,j,i-1) 936 flux_t(k) = w_comp *(&937 37. *( u(k+1,j,i) + u(k,j,i) ) &938 - 8.* ( u(k+2,j,i) + u(k-1,j,i) ) &939 +( u(k+3,j,i) + u(k-2,j,i) ) ) * adv_mom_5940 diss_t(k) = - abs(w_comp) * (&941 10. * ( u(k+1,j,i) - u(k,j,i) )&942 - 5.* ( u(k+2,j,i) - u(k-1,j,i) ) &943 +( u(k+3,j,i) - u(k-2,j,i) ) ) * adv_mom_5943 flux_t(k) = w_comp * ( & 944 37.0 * ( u(k+1,j,i) + u(k,j,i) ) & 945 - 8.0 * ( u(k+2,j,i) + u(k-1,j,i) ) & 946 + ( u(k+3,j,i) + u(k-2,j,i) ) ) * adv_mom_5 947 diss_t(k) = - ABS( w_comp ) * ( & 948 10.0 * ( u(k+1,j,i) - u(k,j,i) ) & 949 - 5.0 * ( u(k+2,j,i) - u(k-1,j,i) ) & 950 + ( u(k+3,j,i) - u(k-2,j,i) ) ) * adv_mom_5 944 951 945 952 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 946 - ( flux_d + diss_d )) * ddzw(k)953 - flux_d - diss_d ) * ddzw(k) 947 954 ENDDO 948 955 ! 949 !-- WS3 as an intermediate step 956 !-- WS3 as an intermediate step (top) 950 957 k = nzt - 1 951 958 flux_d = flux_t(k-1) … … 953 960 w_comp = w(k,j,i) + w(k,j,i-1) 954 961 flux_t(k) = w_comp * ( & 955 7. * ( u(k+1,j,i) + u(k,j,i) )&956 -( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3957 diss_t(k) = - abs(w_comp) * (&958 3. * ( u(k+1,j,i) - u(k,j,i) )&959 -( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3962 7.0 * ( u(k+1,j,i) + u(k,j,i) ) & 963 - ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3 964 diss_t(k) = - ABS( w_comp ) * ( & 965 3.0 * ( u(k+1,j,i) - u(k,j,i) ) & 966 - ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3 960 967 961 968 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 962 - ( flux_d + diss_d )) * ddzw(k)969 - flux_d - diss_d ) * ddzw(k) 963 970 964 971 ! 965 !-- 2nd order scheme 972 !-- 2nd order scheme (top) 966 973 k = nzt 967 974 flux_d = flux_t(k-1) … … 973 980 974 981 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 975 - ( flux_d + diss_d )) * ddzw(k)982 - flux_d - diss_d ) * ddzw(k) 976 983 977 984 ! … … 980 987 sums_wsus_ws_l(k,:) = sums_wsus_ws_l(k,:) & 981 988 + ( flux_t(k) + diss_t(k) ) & 982 * weight_substep(intermediate_timestep_count) * rmask(j,i,:)989 * weight_substep(intermediate_timestep_count) * rmask(j,i,:) 983 990 ENDDO 984 991 … … 1021 1028 1022 1029 CASE ( 1 ) 1023 DO k = nzb_v_inner(j,i) +1, nzt1030 DO k = nzb_v_inner(j,i)+1, nzt 1024 1031 u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu 1025 1032 flux_r(k) = u_comp * ( & 1026 7. * (v(k,j,i+1) + v(k,j,i) )&1027 -( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_31028 diss_r(k) = - abs(u_comp) * (&1029 3. * ( v(k,j,i+1) - v(k,j,i) )&1030 -( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_31033 7.0 * (v(k,j,i+1) + v(k,j,i) ) & 1034 - ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3 1035 diss_r(k) = - ABS( u_comp ) * ( & 1036 3.0 * ( v(k,j,i+1) - v(k,j,i) ) & 1037 - ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3 1031 1038 ENDDO 1032 1039 1033 1040 CASE ( 2 ) 1034 DO k = nzb_v_inner(j,i) +1, nzt1041 DO k = nzb_v_inner(j,i)+1, nzt 1035 1042 u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu 1036 1043 flux_r(k) = u_comp * ( v(k,j,i+1) + v(k,j,i) ) * 0.25 … … 1041 1048 1042 1049 CASE ( 3 ) 1043 DO k = nzb_v_inner(j,i) +1, nzt1050 DO k = nzb_v_inner(j,i)+1, nzt 1044 1051 v_comp(k) = v(k,j+1,i) + v(k,j,i) 1045 1052 flux_n(k) = ( v_comp(k)- gv ) * ( & 1046 7. * ( v(k,j+1,i) + v(k,j,i) )&1047 -( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_31048 diss_n(k) = - abs(v_comp(k) - gv) * (&1049 3. * ( v(k,j+1,i) - v(k,j,i) )&1050 -( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_31053 7.0 * ( v(k,j+1,i) + v(k,j,i) ) & 1054 - ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3 1055 diss_n(k) = - ABS( v_comp(k) - gv ) * ( & 1056 3.0 * ( v(k,j+1,i) - v(k,j,i) ) & 1057 - ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3 1051 1058 ENDDO 1052 1059 1053 1060 CASE ( 4 ) 1054 DO k = nzb_v_inner(j,i) +1, nzt1061 DO k = nzb_v_inner(j,i)+1, nzt 1055 1062 v_comp(k) = v(k,j+1,i) + v(k,j,i) 1056 1063 flux_n(k) = ( v_comp(k) - gv ) * & … … 1062 1069 1063 1070 CASE ( 5 ) 1064 DO k = nzb_w_inner(j,i) +1, nzt1071 DO k = nzb_w_inner(j,i)+1, nzt 1065 1072 u_comp = u(k,j-1,i) + u(k,j,i) - gu 1066 1073 flux_r(k) = u_comp * ( & 1067 7. * ( v(k,j,i+1) + v(k,j,i) )&1068 -( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_31069 diss_r(k) = - abs(u_comp) * (&1070 3. * ( w(k,j,i+1) - w(k,j,i) )&1071 -( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_31074 7.0 * ( v(k,j,i+1) + v(k,j,i) ) & 1075 - ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3 1076 diss_r(k) = - ABS( u_comp ) * ( & 1077 3.0 * ( w(k,j,i+1) - w(k,j,i) ) & 1078 - ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3 1072 1079 ENDDO 1073 1080 1074 1081 CASE ( 6 ) 1075 DO k = nzb_v_inner(j,i) +1, nzt1082 DO k = nzb_v_inner(j,i)+1, nzt 1076 1083 u_comp = u(k,j-1,i) + u(k,j,i) - gu 1077 1084 flux_l_v(k,j) = u_comp * ( v(k,j,i) + v(k,j,i-1) ) * 0.25 … … 1082 1089 u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu 1083 1090 flux_r(k) = u_comp * ( & 1084 7. * ( v(k,j,i+1) + v(k,j,i) )&1085 -( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_31086 diss_r(k) = - abs(u_comp) * (&1087 3. * ( v(k,j,i+1) - v(k,j,i) )&1088 -( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_31091 7.0 * ( v(k,j,i+1) + v(k,j,i) ) & 1092 - ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3 1093 diss_r(k) = - ABS( u_comp ) * ( & 1094 3.0 * ( v(k,j,i+1) - v(k,j,i) ) & 1095 - ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3 1089 1096 ENDDO 1090 1097 degraded_l = .TRUE. 1091 1098 1092 1099 CASE ( 7 ) 1093 DO k = nzb_v_inner(j,i) +1, nzt1100 DO k = nzb_v_inner(j,i)+1, nzt 1094 1101 v_comp(k) = v(k,j,i) + v(k,j-1,i) - gv 1095 1102 flux_s_v(k) = v_comp(k) * ( v(k,j,i) + v(k,j-1,i) ) * 0.25 … … 1100 1107 v_comp(k) = v(k,j+1,i) + v(k,j,i) 1101 1108 flux_n(k) = ( v_comp(k) - gv ) * ( & 1102 7. * ( v(k,j+1,i) + v(k,j,i) )&1103 -( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_31104 diss_n(k) = - abs(v_comp(k) - gv) * (&1105 3. * ( v(k,j+1,i) - v(k,j,i) )&1106 -( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_31109 7.0 * ( v(k,j+1,i) + v(k,j,i) ) & 1110 - ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3 1111 diss_n(k) = - ABS( v_comp(k) - gv ) * ( & 1112 3.0 * ( v(k,j+1,i) - v(k,j,i) ) & 1113 - ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3 1107 1114 ENDDO 1108 1115 degraded_s = .TRUE. … … 1113 1120 ! 1114 1121 !-- Compute the crosswise 5th order fluxes at the outflow 1115 IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 1116 .OR.boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 ) THEN1122 IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 .OR. & 1123 boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 ) THEN 1117 1124 1118 DO k = nzb_v_inner(j,i) +1, nzt1125 DO k = nzb_v_inner(j,i)+1, nzt 1119 1126 v_comp(k) = v(k,j+1,i) + v(k,j,i) 1120 1127 flux_n(k) = ( v_comp(k) - gv ) * ( & 1121 37. 1122 - 8.* ( v(k,j+2,i) + v(k,j-1,i) ) &1123 +( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_51124 diss_n(k) = - abs(v_comp(k) - gv ) * (&1125 10. 1126 - 5.* ( v(k,j+2,i) - v(k,j-1,i) ) &1127 + ( v(k,j+3,i) - v(k,j-2,i) ) ) *adv_mom_51128 37.0 * ( v(k,j+1,i) + v(k,j,i) ) & 1129 - 8.0 * ( v(k,j+2,i) + v(k,j-1,i) ) & 1130 + ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5 1131 diss_n(k) = - ABS( v_comp(k) - gv ) * ( & 1132 10.0 * ( v(k,j+1,i) - v(k,j,i) ) & 1133 - 5.0 * ( v(k,j+2,i) - v(k,j-1,i) ) & 1134 + ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5 1128 1135 ENDDO 1129 1136 1130 1137 ELSE 1131 1138 1132 DO k = nzb_v_inner(j,i) +1, nzt1139 DO k = nzb_v_inner(j,i)+1, nzt 1133 1140 u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu 1134 1141 flux_r(k) = u_comp * ( & 1135 37. 1136 - 8.* ( v(k,j,i+2) + v(k,j,i-1) ) &1137 +( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_51138 diss_r(k) = - abs(u_comp) * (&1139 10. * ( v(k,j,i+1) - v(k,j,i) )&1140 -5. * ( v(k,j,i+2) - v(k,j,i-1) )&1141 +( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_51142 37.0 * ( v(k,j,i+1) + v(k,j,i) ) & 1143 - 8.0 * ( v(k,j,i+2) + v(k,j,i-1) ) & 1144 + ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5 1145 diss_r(k) = - ABS( u_comp ) * ( & 1146 10.0 * ( v(k,j,i+1) - v(k,j,i) ) & 1147 - 5.0 * ( v(k,j,i+2) - v(k,j,i-1) ) & 1148 + ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5 1142 1149 ENDDO 1143 1150 … … 1147 1154 ! 1148 1155 !-- Compute the fifth order fluxes for the interior of PE domain. 1149 DO k = nzb_v_inner(j,i) +1, nzt1156 DO k = nzb_v_inner(j,i)+1, nzt 1150 1157 u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu 1151 1158 flux_r(k) = u_comp * ( & 1152 37. 1153 - 8.* ( v(k,j,i+2) + v(k,j,i-1) ) &1154 +( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_51155 diss_r(k) = - abs(u_comp) * (&1156 10. 1157 - 5.* ( v(k,j,i+2) - v(k,j,i-1) ) &1158 +( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_51159 37.0 * ( v(k,j,i+1) + v(k,j,i) ) & 1160 - 8.0 * ( v(k,j,i+2) + v(k,j,i-1) ) & 1161 + ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5 1162 diss_r(k) = - ABS( u_comp ) * ( & 1163 10.0 * ( v(k,j,i+1) - v(k,j,i) ) & 1164 - 5.0 * ( v(k,j,i+2) - v(k,j,i-1) ) & 1165 + ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5 1159 1166 1160 1167 v_comp(k) = v(k,j+1,i) + v(k,j,i) 1161 1168 flux_n(k) = ( v_comp(k) - gv ) * ( & 1162 37. 1163 - 8.* ( v(k,j+2,i) + v(k,j-1,i) ) &1169 37.0 * ( v(k,j+1,i) + v(k,j,i) ) & 1170 - 8.0 * ( v(k,j+2,i) + v(k,j-1,i) ) & 1164 1171 + ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5 1165 diss_n(k) = - abs(v_comp(k) - gv) * (&1166 10. 1167 - 5.* ( v(k,j+2,i) - v(k,j-1,i) ) &1168 +( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_51172 diss_n(k) = - ABS( v_comp(k) - gv ) * ( & 1173 10.0 * ( v(k,j+1,i) - v(k,j,i) ) & 1174 - 5.0 * ( v(k,j+2,i) - v(k,j-1,i) ) & 1175 + ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5 1169 1176 ENDDO 1170 1177 … … 1173 1180 !-- Compute left- and southside fluxes for the respective boundary 1174 1181 IF ( i == nxl .AND. ( .NOT. degraded_l ) ) THEN 1175 DO k = nzb_v_inner(j,i) +1, nzt1182 DO k = nzb_v_inner(j,i)+1, nzt 1176 1183 u_comp = u(k,j-1,i) + u(k,j,i) - gu 1177 1184 flux_l_v(k,j) = u_comp * ( & 1178 37. 1179 - 8.* ( v(k,j,i+1) + v(k,j,i-2) ) &1180 +( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_51181 diss_l_v(k,j) = - abs(u_comp) * (&1182 10. 1183 - 5.* ( v(k,j,i+1) - v(k,j,i-2) ) &1184 +( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_51185 37.0 * ( v(k,j,i) + v(k,j,i-1) ) & 1186 - 8.0 * ( v(k,j,i+1) + v(k,j,i-2) ) & 1187 + ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5 1188 diss_l_v(k,j) = - ABS( u_comp ) * ( & 1189 10.0 * ( v(k,j,i) - v(k,j,i-1) ) & 1190 - 5.0 * ( v(k,j,i+1) - v(k,j,i-2) ) & 1191 + ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5 1185 1192 ENDDO 1186 1193 … … 1189 1196 IF ( j == nysv .AND. ( .NOT. degraded_s ) ) THEN 1190 1197 1191 DO k = nzb_v_inner(j,i) +1, nzt1198 DO k = nzb_v_inner(j,i)+1, nzt 1192 1199 v_comp_l = v(k,j,i) + v(k,j-1,i) - gv 1193 1200 flux_s_v(k) = v_comp_l * ( & 1194 37. 1195 - 8.* ( v(k,j+1,i) + v(k,j-2,i) ) &1196 +( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_51197 diss_s_v(k) = - abs(v_comp_l) * (&1198 10. 1199 - 5.* ( v(k,j+1,i) - v(k,j-2,i) ) &1200 +( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_51201 37.0 * ( v(k,j,i) + v(k,j-1,i) ) & 1202 - 8.0 * ( v(k,j+1,i) + v(k,j-2,i) ) & 1203 + ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5 1204 diss_s_v(k) = - ABS( v_comp_l ) * ( & 1205 10.0 * ( v(k,j,i) - v(k,j-1,i) ) & 1206 - 5.0 * ( v(k,j+1,i) - v(k,j-2,i) ) & 1207 + ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5 1201 1208 ENDDO 1202 1209 … … 1204 1211 ! 1205 1212 !-- Now compute the tendency terms for the horizontal parts. 1206 DO k = nzb_v_inner(j,i) +1, nzt1213 DO k = nzb_v_inner(j,i)+1, nzt 1207 1214 tend(k,j,i) = tend(k,j,i) - ( & 1208 1215 ( flux_r(k) + diss_r(k) & 1209 - ( flux_l_v(k,j) + diss_l_v(k,j) )) * ddx &1216 - flux_l_v(k,j) - diss_l_v(k,j) ) * ddx & 1210 1217 + ( flux_n(k) + diss_n(k) & 1211 - ( flux_s_v(k) + diss_s_v(k) )) * ddy )1218 - flux_s_v(k) - diss_s_v(k) ) * ddy ) 1212 1219 1213 1220 flux_l_v(k,j) = flux_r(k) … … 1222 1229 sums_vs2_ws_l(k,:) = sums_vs2_ws_l(k,:) & 1223 1230 + ( flux_n(k) & 1224 * ( v_comp(k) - 2.* hom(k,1,2,:) ) &1231 * ( v_comp(k) - 2.0 * hom(k,1,2,:) ) & 1225 1232 / ( v_comp(k) - gv + 1.0E-20 ) & 1226 + diss_n(k)&1227 * abs( v_comp(k) - 2. * hom(k,1,2,:) )&1228 / ( abs( v_comp(k) - gv ) +1.0E-20 ) ) &1229 * weight_substep(intermediate_timestep_count) * rmask(j,i,:)1233 + diss_n(k) & 1234 * ABS( v_comp(k) - 2.0 * hom(k,1,2,:) ) & 1235 / ( ABS( v_comp(k) - gv ) +1.0E-20 ) ) & 1236 * weight_substep(intermediate_timestep_count) * rmask(j,i,:) 1230 1237 1231 1238 ENDDO 1232 sums_vs2_ws_l(nzb_v_inner(j,i),:) = &1233 sums_vs2_ws_l(nzb_v_inner(j,i)+1,:)1239 sums_vs2_ws_l(nzb_v_inner(j,i),:) = sums_vs2_ws_l(nzb_v_inner(j,i)+1,:) 1240 1234 1241 ! 1235 1242 !-- Vertical advection, degradation of order near surface and top. 1236 1243 !-- The fluxes flux_d and diss_d at the surface are 0. Due to reasons of 1237 1244 !-- statistical evaluation the top flux at the surface should be 0 1238 flux_t(nzb_v_inner(j,i)) = 0. !statistical reasons1239 diss_t(nzb_v_inner(j,i)) = 0. 1240 ! 1241 !-- 2nd order scheme 1242 k = nzb_v_inner(j,i) +11245 flux_t(nzb_v_inner(j,i)) = 0.0 !statistical reasons 1246 diss_t(nzb_v_inner(j,i)) = 0.0 1247 ! 1248 !-- 2nd order scheme (bottom) 1249 k = nzb_v_inner(j,i)+1 1243 1250 flux_d = flux_t(k-1) 1244 1251 diss_d = diss_t(k-1) 1245 1252 w_comp = w(k,j-1,i) + w(k,j,i) 1246 1253 flux_t(k) = w_comp * ( v(k+1,j,i) + v(k,j,i) ) * 0.25 1247 diss_t(k) = diss_2nd( v(k+2,j,i), v(k+1,j,i), v(k,j,i), 0. , 0., w_comp,&1248 0.25, ddzw(k) )1254 diss_t(k) = diss_2nd( v(k+2,j,i), v(k+1,j,i), v(k,j,i), 0.0, 0.0, & 1255 w_comp, 0.25, ddzw(k) ) 1249 1256 1250 1257 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 1251 - ( flux_d + diss_d )) * ddzw(k)1258 - flux_d - diss_d ) * ddzw(k) 1252 1259 1253 1260 ! 1254 !-- WS3 as an intermediate step 1255 k = nzb_v_inner(j,i) +21261 !-- WS3 as an intermediate step (bottom) 1262 k = nzb_v_inner(j,i)+2 1256 1263 flux_d = flux_t(k-1) 1257 1264 diss_d = diss_t(k-1) 1258 1265 w_comp = w(k,j-1,i) + w(k,j,i) 1259 1266 flux_t(k) = w_comp * ( & 1260 7. * ( v(k+1,j,i) + v(k,j,i) )&1261 -( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_31262 diss_t(k) = - abs(w_comp) * (&1263 3. * ( v(k+1,j,i) - v(k,j,i) )&1264 -( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_31267 7.0 * ( v(k+1,j,i) + v(k,j,i) ) & 1268 - ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3 1269 diss_t(k) = - ABS( w_comp ) * ( & 1270 3.0 * ( v(k+1,j,i) - v(k,j,i) ) & 1271 - ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3 1265 1272 1266 1273 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 1267 - ( flux_d + diss_d ) ) * ddzw(k) 1274 - flux_d - diss_d ) * ddzw(k) 1275 ! 1268 1276 !-- WS5 1269 DO k = nzb_v_inner(j,i) + 3, nzt -21277 DO k = nzb_v_inner(j,i)+3, nzt-2 1270 1278 flux_d = flux_t(k-1) 1271 1279 diss_d = diss_t(k-1) 1272 1280 w_comp = w(k,j-1,i) + w(k,j,i) 1273 1281 flux_t(k) = w_comp * ( & 1274 37. 1275 - 8.* ( v(k+2,j,i) + v(k-1,j,i) ) &1276 1277 diss_t(k) = - abs(w_comp) * (&1278 10. 1279 - 5.* ( v(k+2,j,i) - v(k-1,j,i) ) &1280 +( v(k+3,j,i) - v(k-2,j,i) ) ) * adv_mom_51282 37.0 * ( v(k+1,j,i) + v(k,j,i) ) & 1283 - 8.0 * ( v(k+2,j,i) + v(k-1,j,i) ) & 1284 + ( v(k+3,j,i) + v(k-2,j,i) ) ) * adv_mom_5 1285 diss_t(k) = - ABS( w_comp ) * ( & 1286 10.0 * ( v(k+1,j,i) - v(k,j,i) ) & 1287 - 5.0 * ( v(k+2,j,i) - v(k-1,j,i) ) & 1288 + ( v(k+3,j,i) - v(k-2,j,i) ) ) * adv_mom_5 1281 1289 1282 1290 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 1283 - ( flux_d + diss_d )) * ddzw(k)1291 - flux_d - diss_d ) * ddzw(k) 1284 1292 ENDDO 1285 1293 ! 1286 !-- WS3 as an intermediate step 1294 !-- WS3 as an intermediate step (top) 1287 1295 k = nzt - 1 1288 1296 flux_d = flux_t(k-1) … … 1290 1298 w_comp = w(k,j-1,i) + w(k,j,i) 1291 1299 flux_t(k) = w_comp * ( & 1292 7. * ( v(k+1,j,i) + v(k,j,i) )&1293 -( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_31294 diss_t(k) = - abs(w_comp) * (&1295 3. * ( v(k+1,j,i) - v(k,j,i) )&1296 -( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_31300 7.0 * ( v(k+1,j,i) + v(k,j,i) ) & 1301 - ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3 1302 diss_t(k) = - ABS( w_comp ) * ( & 1303 3.0 * ( v(k+1,j,i) - v(k,j,i) ) & 1304 - ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3 1297 1305 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 1298 - ( flux_d + diss_d )) * ddzw(k)1299 ! 1300 !-- 2nd order scheme 1306 - flux_d - diss_d ) * ddzw(k) 1307 ! 1308 !-- 2nd order scheme (top) 1301 1309 k = nzt 1302 1310 flux_d = flux_t(k-1) … … 1308 1316 1309 1317 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 1310 - ( flux_d + diss_d )) * ddzw(k)1318 - flux_d - diss_d ) * ddzw(k) 1311 1319 1312 1320 DO k = nzb_v_inner(j,i), nzt 1313 1321 sums_wsvs_ws_l(k,:) = sums_wsvs_ws_l(k,:) & 1314 + ( flux_t(k) + diss_t(k) )&1315 *weight_substep(intermediate_timestep_count) * rmask(j,i,:)1322 + ( flux_t(k) + diss_t(k) ) & 1323 * weight_substep(intermediate_timestep_count) * rmask(j,i,:) 1316 1324 ENDDO 1317 1325 … … 1352 1360 1353 1361 CASE ( 1 ) 1354 DO k = nzb_w_inner(j,i) +1, nzt1362 DO k = nzb_w_inner(j,i)+1, nzt 1355 1363 u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu 1356 1364 flux_r(k) = u_comp * ( & 1357 7. * ( w(k,j,i+1) + w(k,j,i) )&1358 -( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_31359 diss_r(k) = - abs(u_comp) * (&1360 3. * ( w(k,j,i+1) - w(k,j,i) )&1361 -( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_31365 7.0 * ( w(k,j,i+1) + w(k,j,i) ) & 1366 - ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3 1367 diss_r(k) = -ABS( u_comp ) * ( & 1368 3.0 * ( w(k,j,i+1) - w(k,j,i) ) & 1369 - ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3 1362 1370 ENDDO 1363 1371 1364 1372 CASE ( 2 ) 1365 DO k = nzb_w_inner(j,i) +1, nzt1373 DO k = nzb_w_inner(j,i)+1, nzt 1366 1374 u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu 1367 1375 flux_r(k) = u_comp * ( w(k,j,i+1) + w(k,j,i) ) * 0.25 … … 1372 1380 1373 1381 CASE ( 3 ) 1374 DO k = nzb_w_inner(j,i) +1, nzt1382 DO k = nzb_w_inner(j,i)+1, nzt 1375 1383 v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv 1376 1384 flux_n(k) = v_comp * ( & 1377 7. * ( w(k,j+1,i) + w(k,j,i) )&1378 -( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_31379 diss_n(k) = - abs(v_comp) * (&1380 3. * ( w(k,j+1,i) - w(k,j,i) )&1381 -( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_31385 7.0 * ( w(k,j+1,i) + w(k,j,i) ) & 1386 - ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3 1387 diss_n(k) = -ABS( v_comp ) * ( & 1388 3.0 * ( w(k,j+1,i) - w(k,j,i) ) & 1389 - ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3 1382 1390 ENDDO 1383 1391 1384 1392 CASE ( 4 ) 1385 DO k = nzb_w_inner(j,i) +1, nzt1393 DO k = nzb_w_inner(j,i)+1, nzt 1386 1394 v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv 1387 1395 flux_n(k) = v_comp * ( w(k,j+1,i) + w(k,j,i) ) * 0.25 … … 1392 1400 1393 1401 CASE ( 5 ) 1394 DO k = nzb_w_inner(j,i) +1, nzt1402 DO k = nzb_w_inner(j,i)+1, nzt 1395 1403 u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu 1396 1404 flux_r(k) = u_comp * ( & 1397 7. * ( w(k,j,i+1) + w(k,j,i) )&1398 -( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_31399 diss_r(k) = - abs(u_comp) * (&1400 3. * ( w(k,j,i+1) - w(k,j,i) )&1401 -( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_31405 7.0 * ( w(k,j,i+1) + w(k,j,i) ) & 1406 - ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3 1407 diss_r(k) = - ABS( u_comp ) * ( & 1408 3.0 * ( w(k,j,i+1) - w(k,j,i) ) & 1409 - ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3 1402 1410 ENDDO 1403 1411 1404 1412 CASE ( 6 ) 1405 DO k = nzb_w_inner(j,i) +1, nzt1413 DO k = nzb_w_inner(j,i)+1, nzt 1406 1414 ! 1407 1415 !-- Compute leftside fluxes for the left boundary of PE domain 1408 1416 u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu 1409 1417 flux_r(k) = u_comp *( & 1410 7. * ( w(k,j,i+1) + w(k,j,i) )&1411 -( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_31412 diss_r(k) = - abs(u_comp) * (&1413 3. * ( w(k,j,i+1) - w(k,j,i) )&1414 -( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_31418 7.0 * ( w(k,j,i+1) + w(k,j,i) ) & 1419 - ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3 1420 diss_r(k) = - ABS( u_comp ) * ( & 1421 3.0 * ( w(k,j,i+1) - w(k,j,i) ) & 1422 - ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3 1415 1423 1416 1424 u_comp = u(k+1,j,i) + u(k,j,i) - gu … … 1418 1426 diss_l_w(k,j) = diss_2nd( w(k,j,i+2), w(k,j,i+1), w(k,j,i), & 1419 1427 w(k,j,i-1), w(k,j,i-1), u_comp, & 1420 0.25, ddx)1428 0.25, ddx ) 1421 1429 ENDDO 1422 1430 degraded_l = .TRUE. 1423 1431 1424 1432 CASE ( 7 ) 1425 DO k = nzb_w_inner(j,i) +1, nzt1433 DO k = nzb_w_inner(j,i)+1, nzt 1426 1434 v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv 1427 1435 flux_n(k) = v_comp *( & 1428 7. * ( w(k,j+1,i) + w(k,j,i) )&1429 -( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_31430 diss_n(k) = - abs(v_comp) * (&1431 3. * ( w(k,j+1,i) - w(k,j,i) )&1432 -( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_31436 7.0 * ( w(k,j+1,i) + w(k,j,i) ) & 1437 - ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3 1438 diss_n(k) = - ABS( v_comp ) * ( & 1439 3.0 * ( w(k,j+1,i) - w(k,j,i) ) & 1440 - ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3 1433 1441 ENDDO 1434 1442 1435 1443 CASE ( 8 ) 1436 DO k = nzb_w_inner(j,i) +1, nzt1444 DO k = nzb_w_inner(j,i)+1, nzt 1437 1445 v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv 1438 1446 flux_n(k) = v_comp * ( & 1439 7. * ( w(k,j+1,i) + w(k,j,i) )&1440 -( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_31441 diss_n(k) = - abs(v_comp) * (&1442 3. * ( w(k,j+1,i) - w(k,j,i) )&1443 -( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_31447 7.0 * ( w(k,j+1,i) + w(k,j,i) ) & 1448 - ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3 1449 diss_n(k) = - ABS( v_comp ) * ( & 1450 3.0 * ( w(k,j+1,i) - w(k,j,i) ) & 1451 - ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3 1444 1452 ! 1445 1453 !-- Compute southside fluxes for the south boundary of PE domain … … 1457 1465 ! 1458 1466 !-- Compute the crosswise 5th order fluxes at the outflow 1459 IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 1460 .OR.boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 ) THEN1467 IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 .OR. & 1468 boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 ) THEN 1461 1469 1462 DO k = nzb_w_inner(j,i) +1, nzt1470 DO k = nzb_w_inner(j,i)+1, nzt 1463 1471 v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv 1464 1472 flux_n(k) = v_comp * ( & 1465 37. 1466 - 8.* ( w(k,j+2,i) + w(k,j-1,i) ) &1467 +( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_51468 diss_n(k) = - abs(v_comp) * (&1469 10. 1470 - 5.* ( w(k,j+2,i) - w(k,j-1,i) ) &1471 +( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_51473 37.0 * ( w(k,j+1,i) + w(k,j,i) ) & 1474 - 8.0 * ( w(k,j+2,i) + w(k,j-1,i) ) & 1475 + ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5 1476 diss_n(k) = - ABS( v_comp ) * ( & 1477 10.0 * ( w(k,j+1,i) - w(k,j,i) ) & 1478 - 5.0 * ( w(k,j+2,i) - w(k,j-1,i) ) & 1479 + ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5 1472 1480 ENDDO 1473 1481 1474 1482 ELSE 1475 1483 1476 DO k = nzb_w_inner(j,i) +1, nzt1484 DO k = nzb_w_inner(j,i)+1, nzt 1477 1485 u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu 1478 1486 flux_r(k) = u_comp * ( & 1479 37. 1480 - 8.* ( w(k,j,i+2) + w(k,j,i-1) ) &1481 +( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_51482 diss_r(k) = - abs(u_comp) * (&1483 10. 1484 - 5.* ( w(k,j,i+2) - w(k,j,i-1) ) &1485 +( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_51487 37.0 * ( w(k,j,i+1) + w(k,j,i) ) & 1488 - 8.0 * ( w(k,j,i+2) + w(k,j,i-1) ) & 1489 + ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5 1490 diss_r(k) = - ABS( u_comp ) * ( & 1491 10.0 * ( w(k,j,i+1) - w(k,j,i) ) & 1492 - 5.0 * ( w(k,j,i+2) - w(k,j,i-1) ) & 1493 + ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5 1486 1494 ENDDO 1487 1495 … … 1491 1499 ! 1492 1500 !-- Compute the fifth order fluxes for the interior of PE domain. 1493 DO k = nzb_w_inner(j,i) +1, nzt1501 DO k = nzb_w_inner(j,i)+1, nzt 1494 1502 u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu 1495 1503 flux_r(k) = u_comp * ( & 1496 37. 1497 - 8.* ( w(k,j,i+2) + w(k,j,i-1) ) &1498 +( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_51499 diss_r(k) = - abs(u_comp) * (&1500 10. 1501 - 5.* ( w(k,j,i+2) - w(k,j,i-1) ) &1502 +( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_51504 37.0 * ( w(k,j,i+1) + w(k,j,i) ) & 1505 - 8.0 * ( w(k,j,i+2) + w(k,j,i-1) ) & 1506 + ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5 1507 diss_r(k) = - ABS( u_comp ) * ( & 1508 10.0 * ( w(k,j,i+1) - w(k,j,i) ) & 1509 - 5.0 * ( w(k,j,i+2) - w(k,j,i-1) ) & 1510 + ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5 1503 1511 1504 1512 v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv 1505 1513 flux_n(k) = v_comp * ( & 1506 37. 1507 - 8.* ( w(k,j+2,i) + w(k,j-1,i) ) &1508 +( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_51509 diss_n(k) = - abs(v_comp) * (&1510 10. 1511 - 5.* ( w(k,j+2,i) - w(k,j-1,i) ) &1512 +( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_51514 37.0 * ( w(k,j+1,i) + w(k,j,i) ) & 1515 - 8.0 * ( w(k,j+2,i) + w(k,j-1,i) ) & 1516 + ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5 1517 diss_n(k) = - ABS( v_comp ) * ( & 1518 10.0 * ( w(k,j+1,i) - w(k,j,i) ) & 1519 - 5.0 * ( w(k,j+2,i) - w(k,j-1,i) ) & 1520 + ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5 1513 1521 ENDDO 1514 1522 … … 1518 1526 IF ( j == nys .AND. ( .NOT. degraded_s ) ) THEN 1519 1527 1520 DO k = nzb_w_inner(j,i) +1, nzt1528 DO k = nzb_w_inner(j,i)+1, nzt 1521 1529 v_comp = v(k+1,j,i) + v(k,j,i) - gv 1522 1530 flux_s_w(k) = v_comp * ( & 1523 37. 1524 - 8.* ( w(k,j+1,i) +w(k,j-2,i) ) &1525 +( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_51526 diss_s_w(k) = - abs(v_comp) * (&1527 10. 1528 - 5.* ( w(k,j+1,i) - w(k,j-2,i) ) &1529 +( w(k,j+2,i) - w(k,j-3,i) ) ) * adv_mom_51531 37.0 * ( w(k,j,i) + w(k,j-1,i) ) & 1532 - 8.0 * ( w(k,j+1,i) +w(k,j-2,i) ) & 1533 + ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5 1534 diss_s_w(k) = - ABS( v_comp ) * ( & 1535 10.0 * ( w(k,j,i) - w(k,j-1,i) ) & 1536 - 5.0 * ( w(k,j+1,i) - w(k,j-2,i) ) & 1537 + ( w(k,j+2,i) - w(k,j-3,i) ) ) * adv_mom_5 1530 1538 ENDDO 1531 1539 … … 1534 1542 IF ( i == nxl .AND. ( .NOT. degraded_l ) ) THEN 1535 1543 1536 DO k = nzb_w_inner(j,i) +1, nzt1544 DO k = nzb_w_inner(j,i)+1, nzt 1537 1545 u_comp = u(k+1,j,i) + u(k,j,i) - gu 1538 1546 flux_l_w(k,j) = u_comp * ( & 1539 37. 1540 - 8.* ( w(k,j,i+1) + w(k,j,i-2) ) &1541 +( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_51542 diss_l_w(k,j) = - abs(u_comp) * (&1543 10. 1544 - 5.* ( w(k,j,i+1) - w(k,j,i-2) ) &1545 +( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_51547 37.0 * ( w(k,j,i) + w(k,j,i-1) ) & 1548 - 8.0 * ( w(k,j,i+1) + w(k,j,i-2) ) & 1549 + ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5 1550 diss_l_w(k,j) = - ABS( u_comp ) * ( & 1551 10.0 * ( w(k,j,i) - w(k,j,i-1) ) & 1552 - 5.0 * ( w(k,j,i+1) - w(k,j,i-2) ) & 1553 + ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5 1546 1554 ENDDO 1547 1555 … … 1549 1557 ! 1550 1558 !-- Now compute the tendency terms for the horizontal parts. 1551 DO k = nzb_w_inner(j,i) +1, nzt1559 DO k = nzb_w_inner(j,i)+1, nzt 1552 1560 tend(k,j,i) = tend(k,j,i) - ( & 1553 1561 ( flux_r(k) + diss_r(k) & 1554 - ( flux_l_w(k,j) + diss_l_w(k,j) )) * ddx &1562 - flux_l_w(k,j) - diss_l_w(k,j) ) * ddx & 1555 1563 + ( flux_n(k) + diss_n(k) & 1556 - ( flux_s_w(k) + diss_s_w(k) )) * ddy )1564 - flux_s_w(k) - diss_s_w(k) ) * ddy ) 1557 1565 1558 1566 flux_l_w(k,j) = flux_r(k) … … 1566 1574 !-- The fluxes flux_d and diss_d at the surface are 0. Due to reasons of 1567 1575 !-- statistical evaluation the top flux at the surface should be 0 1568 flux_t(nzb_w_inner(j,i)) = 0. !statistical reasons1569 diss_t(nzb_w_inner(j,i)) = 0. 1570 ! 1571 !-- 2nd order scheme 1572 k = nzb_w_inner(j,i) +11576 flux_t(nzb_w_inner(j,i)) = 0.0 !statistical reasons 1577 diss_t(nzb_w_inner(j,i)) = 0.0 1578 ! 1579 !-- 2nd order scheme (bottom) 1580 k = nzb_w_inner(j,i)+1 1573 1581 flux_d = flux_t(k-1) 1574 1582 diss_d = diss_t(k-1) 1575 1583 w_comp = w(k+1,j,i) + w(k,j,i) 1576 1584 flux_t(k) = w_comp * ( w(k+1,j,i) + w(k,j,i) ) * 0.25 1577 diss_t(k) = diss_2nd( w(k+2,j,i), w(k+1,j,i), w(k,j,i), 0. , 0., &1585 diss_t(k) = diss_2nd( w(k+2,j,i), w(k+1,j,i), w(k,j,i), 0.0, 0.0, & 1578 1586 w_comp, 0.25, ddzu(k+1) ) 1579 1587 1580 1588 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 1581 - ( flux_d + diss_d )) * ddzu(k+1)1582 ! 1583 !-- WS3 as an intermediate step 1584 k = nzb_w_inner(j,i) +21589 - flux_d - diss_d ) * ddzu(k+1) 1590 ! 1591 !-- WS3 as an intermediate step (bottom) 1592 k = nzb_w_inner(j,i)+2 1585 1593 flux_d = flux_t(k-1) 1586 1594 diss_d = diss_t(k-1) 1587 1595 w_comp = w(k+1,j,i) + w(k,j,i) 1588 1596 flux_t(k) = w_comp * ( & 1589 7. * ( w(k+1,j,i) + w(k,j,i) )&1590 -( w(k+2,j,i) + w(k-1,j,i) ) ) * adv_mom_31591 diss_t(k) = - abs(w_comp) * (&1592 3. * ( w(k+1,j,i) - w(k,j,i) )&1593 -( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_31597 7.0 * ( w(k+1,j,i) + w(k,j,i) ) & 1598 - ( w(k+2,j,i) + w(k-1,j,i) ) ) * adv_mom_3 1599 diss_t(k) = - ABS( w_comp ) * ( & 1600 3.0 * ( w(k+1,j,i) - w(k,j,i) ) & 1601 - ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3 1594 1602 1595 1603 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 1596 - ( flux_d + diss_d )) * ddzu(k+1)1604 - flux_d - diss_d ) * ddzu(k+1) 1597 1605 ! 1598 1606 !-- WS5 1599 DO k = nzb_w_inner(j,i) + 3, nzt-21607 DO k = nzb_w_inner(j,i)+3, nzt-2 1600 1608 flux_d = flux_t(k-1) 1601 1609 diss_d = diss_t(k-1) 1602 1610 w_comp = w(k+1,j,i) + w(k,j,i) 1603 1611 flux_t(k) = w_comp * ( & 1604 37. 1605 - 8.* ( w(k+2,j,i) + w(k-1,j,i) ) &1606 +( w(k+3,j,i) + w(k-2,j,i) ) ) * adv_mom_51607 diss_t(k) = - abs(w_comp) * (&1608 10. 1609 - 5.* ( w(k+2,j,i) - w(k-1,j,i) ) &1610 +( w(k+3,j,i) - w(k-2,j,i) ) ) * adv_mom_51612 37.0 * ( w(k+1,j,i) + w(k,j,i) ) & 1613 - 8.0 * ( w(k+2,j,i) + w(k-1,j,i) ) & 1614 + ( w(k+3,j,i) + w(k-2,j,i) ) ) * adv_mom_5 1615 diss_t(k) = - ABS( w_comp ) * ( & 1616 10.0 * ( w(k+1,j,i) - w(k,j,i) ) & 1617 - 5.0 * ( w(k+2,j,i) - w(k-1,j,i) ) & 1618 + ( w(k+3,j,i) - w(k-2,j,i) ) ) * adv_mom_5 1611 1619 1612 1620 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 1613 - ( flux_d + diss_d )) * ddzu(k+1)1621 - flux_d - diss_d ) * ddzu(k+1) 1614 1622 ENDDO 1615 !-- WS3 as an intermediate step 1623 !-- WS3 as an intermediate step (top) 1616 1624 k = nzt - 1 1617 1625 flux_d = flux_t(k-1) … … 1619 1627 w_comp = w(k+1,j,i) + w(k,j,i) 1620 1628 flux_t(k) = w_comp * ( & 1621 7. * ( w(k+1,j,i) + w(k,j,i) )&1622 -( w(k+2,j,i) + w(k-1,j,i) ) ) *adv_mom_31623 diss_t(k) = - abs(w_comp) * (&1624 3. * ( w(k+1,j,i) - w(k,j,i) )&1625 -( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_31629 7.0 * ( w(k+1,j,i) + w(k,j,i) ) & 1630 - ( w(k+2,j,i) + w(k-1,j,i) ) ) *adv_mom_3 1631 diss_t(k) = - ABS( w_comp ) * ( & 1632 3.0 * ( w(k+1,j,i) - w(k,j,i) ) & 1633 - ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3 1626 1634 1627 1635 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 1628 - ( flux_d + diss_d )) * ddzu(k+1)1629 ! 1630 !-- 2nd order scheme 1636 - flux_d - diss_d ) * ddzu(k+1) 1637 ! 1638 !-- 2nd order scheme (top) 1631 1639 k = nzt 1632 1640 flux_d = flux_t(k-1) … … 1638 1646 1639 1647 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 1640 - ( flux_d + diss_d )) * ddzu(k+1)1648 - flux_d - diss_d ) * ddzu(k+1) 1641 1649 1642 1650 DO k = nzb_w_inner(j,i), nzt 1643 1651 sums_ws2_ws_l(k,:) = sums_ws2_ws_l(k,:) & 1644 + ( flux_t(k) + diss_t(k) )&1645 *weight_substep(intermediate_timestep_count) * rmask(j,i,:)1652 + ( flux_t(k) + diss_t(k) ) & 1653 * weight_substep(intermediate_timestep_count) * rmask(j,i,:) 1646 1654 ENDDO 1647 1655 … … 1649 1657 1650 1658 1651 ! 1659 !------------------------------------------------------------------------------! 1652 1660 ! Scalar advection - Call for all grid points 1653 1661 !------------------------------------------------------------------------------! … … 1680 1688 IF ( boundary_flags(j,i) == 6 ) THEN 1681 1689 1682 DO k = nzb_s_inner(j,i) +1, nzt1690 DO k = nzb_s_inner(j,i)+1, nzt 1683 1691 u_comp = u(k,j,i) - u_gtrans 1684 1692 swap_flux_x_local(k,j) = u_comp * ( & … … 1692 1700 ELSE 1693 1701 1694 DO k = nzb_s_inner(j,i) +1, nzt1702 DO k = nzb_s_inner(j,i)+1, nzt 1695 1703 u_comp = u(k,j,i) - u_gtrans 1696 swap_flux_x_local(k,j) = u_comp*( &1697 37. * (sk(k,j,i)+sk(k,j,i-1) )&1698 -8.* ( sk(k,j,i+1) + sk(k,j,i-2) ) &1699 +( sk(k,j,i+2) + sk(k,j,i-3) ) )&1700 1701 swap_diss_x_local(k,j) = - abs(u_comp) * (&1702 10. * (sk(k,j,i) - sk(k,j,i-1) ) &1703 -5.* ( sk(k,j,i+1) - sk(k,j,i-2) ) &1704 +( sk(k,j,i+2) - sk(k,j,i-3) ) )&1705 1704 swap_flux_x_local(k,j) = u_comp*( & 1705 37.0 * ( sk(k,j,i)+sk(k,j,i-1) ) & 1706 - 8.0 * ( sk(k,j,i+1) + sk(k,j,i-2) ) & 1707 + ( sk(k,j,i+2) + sk(k,j,i-3) ) )& 1708 * adv_sca_5 1709 swap_diss_x_local(k,j) = - ABS( u_comp ) * ( & 1710 10.0 * (sk(k,j,i) - sk(k,j,i-1) ) & 1711 - 5.0 * ( sk(k,j,i+1) - sk(k,j,i-2) ) & 1712 + ( sk(k,j,i+2) - sk(k,j,i-3) ) )& 1713 * adv_sca_5 1706 1714 ENDDO 1707 1715 ENDIF … … 1715 1723 IF ( boundary_flags(j,i) == 8 ) THEN 1716 1724 1717 DO k = nzb_s_inner(j,i) +1, nzt1725 DO k = nzb_s_inner(j,i)+1, nzt 1718 1726 v_comp = v(k,j,i) - v_gtrans 1719 swap_flux_y_local(k) = v_comp * &1727 swap_flux_y_local(k) = v_comp * & 1720 1728 ( sk(k,j,i) + sk(k,j-1,i) ) * 0.5 1721 swap_diss_y_local(k) = diss_2nd( sk(k,j+2,i), sk(k,j+1,i), &1722 sk(k,j,i), sk(k,j-1,i), &1723 sk(k,j-1,i), v_comp, 0.5, ddy )1729 swap_diss_y_local(k) = diss_2nd( sk(k,j+2,i), sk(k,j+1,i), & 1730 sk(k,j,i), sk(k,j-1,i), & 1731 sk(k,j-1,i), v_comp, 0.5, ddy ) 1724 1732 ENDDO 1725 1733 1726 1734 ELSE 1727 1735 1728 DO k = nzb_s_inner(j,i) +1, nzt1736 DO k = nzb_s_inner(j,i)+1, nzt 1729 1737 v_comp = v(k,j,i) - v_gtrans 1730 swap_flux_y_local(k) = v_comp * ( &1731 37. * ( sk(k,j,i) + sk(k,j-1,i) )&1732 - 8. * ( sk(k,j+1,i) + sk(k,j-2,i) )&1733 + ( sk(k,j+2,i) + sk(k,j-3,i) ) )&1734 1735 swap_diss_y_local(k)= - abs(v_comp) * (&1736 10. * ( sk(k,j,i) - sk(k,j-1,i) )&1737 - 5. * ( sk(k,j+1,i) - sk(k,j-2,i) )&1738 + ( sk(k,j+2,i)-sk(k,j-3,i) ) )&1739 1738 swap_flux_y_local(k) = v_comp * ( & 1739 37.0 * ( sk(k,j,i) + sk(k,j-1,i) ) & 1740 - 8.0 * ( sk(k,j+1,i) + sk(k,j-2,i) ) & 1741 + ( sk(k,j+2,i) + sk(k,j-3,i) ) ) & 1742 * adv_sca_5 1743 swap_diss_y_local(k)= - ABS( v_comp ) * ( & 1744 10.0 * ( sk(k,j,i) - sk(k,j-1,i) ) & 1745 - 5.0 * ( sk(k,j+1,i) - sk(k,j-2,i) ) & 1746 + ( sk(k,j+2,i)-sk(k,j-3,i) ) ) & 1747 * adv_sca_5 1740 1748 ENDDO 1741 1749 … … 1749 1757 1750 1758 CASE ( 1 ) 1751 DO k = nzb_s_inner(j,i) +1, nzt1759 DO k = nzb_s_inner(j,i)+1, nzt 1752 1760 u_comp = u(k,j,i+1) - u_gtrans 1753 flux_r(k) = u_comp * ( & 1754 7. * ( sk(k,j,i+1) + sk(k,j,i) ) & 1755 - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3 1756 diss_r(k) = - abs(u_comp) * ( & 1757 3. * ( sk(k,j,i+1) - sk(k,j,i) ) & 1758 - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3 1761 flux_r(k) = u_comp * ( & 1762 7.0 * ( sk(k,j,i+1) + sk(k,j,i) ) & 1763 - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) & 1764 * adv_sca_3 1765 diss_r(k) = - ABS( u_comp ) * ( & 1766 3.0 * ( sk(k,j,i+1) - sk(k,j,i) ) & 1767 - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) & 1768 * adv_sca_3 1759 1769 ENDDO 1760 1770 1761 1771 CASE ( 2 ) 1762 DO k = nzb_s_inner(j,i) +1, nzt1772 DO k = nzb_s_inner(j,i)+1, nzt 1763 1773 u_comp = u(k,j,i+1) - u_gtrans 1764 1774 flux_r(k) = u_comp * ( sk(k,j,i+1) + sk(k,j,i) ) * 0.5 1765 diss_r(k) = diss_2nd( sk(k,j,i+1), sk(k,j,i+1), &1766 sk(k,j,i), sk(k,j,i-1), &1775 diss_r(k) = diss_2nd( sk(k,j,i+1), sk(k,j,i+1), & 1776 sk(k,j,i), sk(k,j,i-1), & 1767 1777 sk(k,j,i-2), u_comp, 0.5, ddx ) 1768 1778 ENDDO 1769 1779 1770 1780 CASE ( 3 ) 1771 DO k = nzb_s_inner(j,i) +1, nzt1781 DO k = nzb_s_inner(j,i)+1, nzt 1772 1782 v_comp = v(k,j+1,i) - v_gtrans 1773 flux_n(k) = v_comp * ( & 1774 7. * ( sk(k,j+1,i) + sk(k,j,i) ) & 1775 - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3 1776 diss_n(k) = - abs(v_comp) * ( & 1777 3. * ( sk(k,j+1,i) - sk(k,j,i) ) & 1778 - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3 1783 flux_n(k) = v_comp * ( & 1784 7.0 * ( sk(k,j+1,i) + sk(k,j,i) ) & 1785 - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) & 1786 * adv_sca_3 1787 diss_n(k) = - ABS( v_comp ) * ( & 1788 3.0 * ( sk(k,j+1,i) - sk(k,j,i) ) & 1789 - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) & 1790 * adv_sca_3 1779 1791 ENDDO 1780 1792 1781 1793 CASE ( 4 ) 1782 DO k = nzb_s_inner(j,i) +1, nzt1794 DO k = nzb_s_inner(j,i)+1, nzt 1783 1795 v_comp = v(k,j+1,i) - v_gtrans 1784 1796 flux_n(k) = v_comp * ( sk(k,j+1,i) + sk(k,j,i) ) * 0.5 1785 diss_n(k) = diss_2nd( sk(k,j+1,i), sk(k,j+1,i), &1786 sk(k,j,i), sk(k,j-1,i), &1797 diss_n(k) = diss_2nd( sk(k,j+1,i), sk(k,j+1,i), & 1798 sk(k,j,i), sk(k,j-1,i), & 1787 1799 sk(k,j-2,i), v_comp, 0.5, ddy ) 1788 1800 ENDDO 1789 1801 1790 1802 CASE ( 5 ) 1791 DO k = nzb_w_inner(j,i) +1, nzt1803 DO k = nzb_w_inner(j,i)+1, nzt 1792 1804 u_comp = u(k,j,i+1) - u_gtrans 1793 flux_r(k) = u_comp * ( & 1794 7. * ( sk(k,j,i+1) + sk(k,j,i) ) & 1795 - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3 1796 diss_r(k) = - abs(u_comp) * ( & 1797 3. * ( sk(k,j,i+1) - sk(k,j,i) ) & 1798 - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3 1805 flux_r(k) = u_comp * ( & 1806 7.0 * ( sk(k,j,i+1) + sk(k,j,i) ) & 1807 - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) & 1808 * adv_sca_3 1809 diss_r(k) = - ABS( u_comp ) * ( & 1810 3.0 * ( sk(k,j,i+1) - sk(k,j,i) ) & 1811 - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) & 1812 * adv_sca_3 1799 1813 ENDDO 1800 1814 1801 1815 CASE ( 6 ) 1802 DO k = nzb_s_inner(j,i) +1, nzt1816 DO k = nzb_s_inner(j,i)+1, nzt 1803 1817 u_comp = u(k,j,i+1) - u_gtrans 1804 flux_r(k) = u_comp * ( & 1805 7. * ( sk(k,j,i+1) + sk(k,j,i) ) & 1806 - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3 1807 diss_r(k) = - abs(u_comp) * ( & 1808 3. * ( sk(k,j,i+1) - sk(k,j,i) ) & 1809 - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3 1818 flux_r(k) = u_comp * ( & 1819 7.0 * ( sk(k,j,i+1) + sk(k,j,i) ) & 1820 - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) & 1821 * adv_sca_3 1822 diss_r(k) = - ABS( u_comp ) * ( & 1823 3.0 * ( sk(k,j,i+1) - sk(k,j,i) ) & 1824 - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) & 1825 * adv_sca_3 1810 1826 ENDDO 1811 1827 1812 1828 CASE ( 7 ) 1813 DO k = nzb_s_inner(j,i) +1, nzt1829 DO k = nzb_s_inner(j,i)+1, nzt 1814 1830 v_comp = v(k,j+1,i) - v_gtrans 1815 flux_n(k) = v_comp * ( & 1816 7. * ( sk(k,j+1,i) + sk(k,j,i) ) & 1817 - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3 1818 diss_n(k) = - abs(v_comp) * ( & 1819 3. * ( sk(k,j+1,i) - sk(k,j,i) ) & 1820 - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3 1831 flux_n(k) = v_comp * ( & 1832 7.0 * ( sk(k,j+1,i) + sk(k,j,i) ) & 1833 - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) & 1834 * adv_sca_3 1835 diss_n(k) = - ABS( v_comp ) * ( & 1836 3.0 * ( sk(k,j+1,i) - sk(k,j,i) ) & 1837 - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) & 1838 * adv_sca_3 1821 1839 ENDDO 1822 1840 1823 1841 CASE ( 8 ) 1824 DO k = nzb_s_inner(j,i) +1, nzt1842 DO k = nzb_s_inner(j,i)+1, nzt 1825 1843 v_comp = v(k,j+1,i) - v_gtrans 1826 flux_n(k) = v_comp * ( & 1827 7. * ( sk(k,j+1,i) + sk(k,j,i) ) & 1828 - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3 1829 diss_n(k) = - abs(v_comp) * ( & 1830 3. * ( sk(k,j+1,i) - sk(k,j,i) ) & 1831 - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3 1844 flux_n(k) = v_comp * ( & 1845 7.0 * ( sk(k,j+1,i) + sk(k,j,i) ) & 1846 - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) & 1847 * adv_sca_3 1848 diss_n(k) = - ABS( v_comp ) * ( & 1849 3.0 * ( sk(k,j+1,i) - sk(k,j,i) ) & 1850 - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) & 1851 * adv_sca_3 1832 1852 ENDDO 1833 1853 … … 1836 1856 END SELECT 1837 1857 1838 IF (boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 & 1839 .OR. boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6) THEN 1858 IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 .OR.& 1859 boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 ) & 1860 THEN 1840 1861 1841 DO k = nzb_s_inner(j,i) +1, nzt1862 DO k = nzb_s_inner(j,i)+1, nzt 1842 1863 v_comp = v(k,j+1,i) - v_gtrans 1843 flux_n(k) = v_comp * ( &1844 37. * ( sk(k,j+1,i) + sk(k,j,i)) &1845 - 8. * (sk(k,j+2,i) + sk(k,j-1,i) ) &1846 + ( sk(k,j+3,i) + sk(k,j-2,i) ) )&1847 1848 diss_n(k) = - abs(v_comp) * (&1849 10. * ( sk(k,j+1,i) - sk(k,j,i) )&1850 - 5. * ( sk(k,j+2,i) - sk(k,j-1,i) )&1851 + ( sk(k,j+3,i) - sk(k,j-2,i) ) )&1852 1864 flux_n(k) = v_comp * ( & 1865 37.0 * ( sk(k,j+1,i) + sk(k,j,i) ) & 1866 - 8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) ) & 1867 + ( sk(k,j+3,i) + sk(k,j-2,i) ) ) & 1868 * adv_sca_5 1869 diss_n(k) = - ABS( v_comp ) * ( & 1870 10.0 * ( sk(k,j+1,i) - sk(k,j,i) ) & 1871 - 5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) ) & 1872 + ( sk(k,j+3,i) - sk(k,j-2,i) ) ) & 1873 * adv_sca_5 1853 1874 ENDDO 1854 1875 1855 1876 ELSE 1856 1877 1857 DO k = nzb_s_inner(j,i) +1, nzt1878 DO k = nzb_s_inner(j,i)+1, nzt 1858 1879 u_comp = u(k,j,i+1) - u_gtrans 1859 flux_r(k) = u_comp * ( &1860 37. * ( sk(k,j,i+1) + sk(k,j,i) )&1861 - 8. * ( sk(k,j,i+2) + sk(k,j,i-1) )&1862 + ( sk(k,j,i+3) + sk(k,j,i-2) ) )&1863 1864 diss_r(k) = - abs(u_comp) * (&1865 10. * ( sk(k,j,i+1) - sk(k,j,i) )&1866 - 5. * ( sk(k,j,i+2) - sk(k,j,i-1) )&1867 + ( sk(k,j,i+3) - sk(k,j,i-2) ) )&1868 1880 flux_r(k) = u_comp * ( & 1881 37.0 * ( sk(k,j,i+1) + sk(k,j,i) ) & 1882 - 8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) ) & 1883 + ( sk(k,j,i+3) + sk(k,j,i-2) ) ) & 1884 * adv_sca_5 1885 diss_r(k) = - ABS( u_comp ) * ( & 1886 10.0 * ( sk(k,j,i+1) - sk(k,j,i) ) & 1887 - 5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) ) & 1888 + ( sk(k,j,i+3) - sk(k,j,i-2) ) ) & 1889 * adv_sca_5 1869 1890 ENDDO 1870 1891 … … 1873 1894 ELSE 1874 1895 1875 DO k = nzb_s_inner(j,i) +1, nzt1896 DO k = nzb_s_inner(j,i)+1, nzt 1876 1897 u_comp = u(k,j,i+1) - u_gtrans 1877 flux_r(k) = u_comp * ( &1878 37. * ( sk(k,j,i+1) + sk(k,j,i) )&1879 - 8. * ( sk(k,j,i+2) + sk(k,j,i-1) )&1880 + ( sk(k,j,i+3) + sk(k,j,i-2) ) )&1881 1882 diss_r(k) = - abs(u_comp) * (&1883 10. * ( sk(k,j,i+1) - sk(k,j,i) )&1884 - 5. * ( sk(k,j,i+2) - sk(k,j,i-1) )&1885 + ( sk(k,j,i+3) - sk(k,j,i-2) ) )&1886 1898 flux_r(k) = u_comp * ( & 1899 37.0 * ( sk(k,j,i+1) + sk(k,j,i) ) & 1900 - 8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) ) & 1901 + ( sk(k,j,i+3) + sk(k,j,i-2) ) ) & 1902 * adv_sca_5 1903 diss_r(k) = - ABS( u_comp ) * ( & 1904 10.0 * ( sk(k,j,i+1) - sk(k,j,i) ) & 1905 - 5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) ) & 1906 + ( sk(k,j,i+3) - sk(k,j,i-2) ) ) & 1907 * adv_sca_5 1887 1908 1888 1909 v_comp = v(k,j+1,i) - v_gtrans 1889 flux_n(k) = v_comp * ( &1890 37. * ( sk(k,j+1,i) + sk(k,j,i) )&1891 - 8. * ( sk(k,j+2,i) + sk(k,j-1,i) )&1892 + ( sk(k,j+3,i) + sk(k,j-2,i) ) )&1893 1894 diss_n(k) = - abs(v_comp) * (&1895 10. * ( sk(k,j+1,i) - sk(k,j,i) )&1896 - 5. * ( sk(k,j+2,i) - sk(k,j-1,i) )&1897 + ( sk(k,j+3,i) - sk(k,j-2,i) ) )&1898 1910 flux_n(k) = v_comp * ( & 1911 37.0 * ( sk(k,j+1,i) + sk(k,j,i) ) & 1912 - 8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) ) & 1913 + ( sk(k,j+3,i) + sk(k,j-2,i) ) ) & 1914 * adv_sca_5 1915 diss_n(k) = - ABS( v_comp ) * ( & 1916 10.0 * ( sk(k,j+1,i) - sk(k,j,i) ) & 1917 - 5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) ) & 1918 + ( sk(k,j+3,i) - sk(k,j-2,i) ) ) & 1919 * adv_sca_5 1899 1920 ENDDO 1900 1921 1901 1922 ENDIF 1902 1923 1903 DO k = nzb_s_inner(j,i) +1, nzt1924 DO k = nzb_s_inner(j,i)+1, nzt 1904 1925 tend(k,j,i) = tend(k,j,i) - ( & 1905 1926 ( flux_r(k) + diss_r(k) & 1906 - ( swap_flux_x_local(k,j) + swap_diss_x_local(k,j) )) * ddx &1927 - swap_flux_x_local(k,j) - swap_diss_x_local(k,j) ) * ddx & 1907 1928 + ( flux_n(k) + diss_n(k) & 1908 - ( swap_flux_y_local(k) + swap_diss_y_local(k) )) * ddy)1929 - swap_flux_y_local(k) - swap_diss_y_local(k) ) * ddy) 1909 1930 1910 1931 swap_flux_x_local(k,j) = flux_r(k) … … 1923 1944 DO j = nys, nyn 1924 1945 ! 1925 !-- 2nd order scheme 1926 k=nzb_s_inner(j,i) +11946 !-- 2nd order scheme (bottom) 1947 k=nzb_s_inner(j,i)+1 1927 1948 ! 1928 1949 !-- The fluxes flux_d and diss_d at the surface are 0. Due to static 1929 1950 !-- reasons the top flux at the surface should be 0. 1930 flux_t(nzb_s_inner(j,i)) = 0. 1931 diss_t(nzb_s_inner(j,i)) = 0. 1951 flux_t(nzb_s_inner(j,i)) = 0.0 1952 diss_t(nzb_s_inner(j,i)) = 0.0 1932 1953 flux_d = flux_t(k-1) 1933 1954 diss_d = diss_t(k-1) 1934 flux_t(k) = w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) *0.51955 flux_t(k) = w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) * 0.5 1935 1956 diss_t(k) = diss_2nd( sk(k+2,j,i), sk(k+1,j,i), sk(k,j,i), & 1936 1957 sk(k,j,i), sk(k,j,i), w(k,j,i), & 1937 1958 0.5, ddzw(k) ) 1938 1959 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 1939 - ( flux_d + diss_d )) * ddzw(k)1940 ! 1941 !-- WS3 as an intermediate step 1942 k = nzb_s_inner(j,i) +21960 - flux_d - diss_d ) * ddzw(k) 1961 ! 1962 !-- WS3 as an intermediate step (bottom) 1963 k = nzb_s_inner(j,i)+2 1943 1964 flux_d = flux_t(k-1) 1944 1965 diss_d = diss_t(k-1) 1945 1966 flux_t(k) = w(k,j,i) * ( & 1946 7. * ( sk(k+1,j,i) + sk(k,j,i) )&1967 7.0 * ( sk(k+1,j,i) + sk(k,j,i) ) & 1947 1968 - ( sk(k+2,j,i) + sk(k-1,j,i) ) ) * adv_sca_3 1948 diss_t(k) = - abs(w(k,j,i)) * (&1949 3. * ( sk(k+1,j,i) - sk(k,j,i) )&1950 -( sk(k+2,j,i) - sk(k-1,j,i) ) ) * adv_sca_31969 diss_t(k) = - ABS( w(k,j,i) ) * ( & 1970 3.0 * ( sk(k+1,j,i) - sk(k,j,i) ) & 1971 - ( sk(k+2,j,i) - sk(k-1,j,i) ) ) * adv_sca_3 1951 1972 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 1952 - ( flux_d + diss_d )) * ddzw(k)1973 - flux_d - diss_d ) * ddzw(k) 1953 1974 ! 1954 1975 !-- WS5 1955 DO k = nzb_s_inner(j,i) + 3, nzt -21976 DO k = nzb_s_inner(j,i)+3, nzt-2 1956 1977 flux_d = flux_t(k-1) 1957 1978 diss_d = diss_t(k-1) 1958 1979 flux_t(k) = w(k,j,i) * ( & 1959 37. 1960 - 8.* ( sk(k+2,j,i) + sk(k-1,j,i) ) &1961 +( sk(k+3,j,i) + sk(k-2,j,i) ) ) * adv_sca_51962 diss_t(k) = - abs(w(k,j,i)) * ( &1963 10. * ( sk(k+1,j,i) -sk(k,j,i) )&1964 - 5.* ( sk(k+2,j,i) - sk(k-1,j,i) ) &1965 +( sk(k+3,j,i) - sk(k-2,j,i) ) ) * adv_sca_51980 37.0 * ( sk(k+1,j,i) + sk(k,j,i) ) & 1981 - 8.0 * ( sk(k+2,j,i) + sk(k-1,j,i) ) & 1982 + ( sk(k+3,j,i) + sk(k-2,j,i) ) ) * adv_sca_5 1983 diss_t(k) = - ABS(w(k,j,i)) * ( & 1984 10.0 * ( sk(k+1,j,i) -sk(k,j,i) ) & 1985 - 5.0 * ( sk(k+2,j,i) - sk(k-1,j,i) ) & 1986 + ( sk(k+3,j,i) - sk(k-2,j,i) ) ) * adv_sca_5 1966 1987 1967 1988 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 1968 - ( flux_d + diss_d )) * ddzw(k)1989 - flux_d - diss_d ) * ddzw(k) 1969 1990 ENDDO 1970 1991 ! 1971 !-- WS3 as an intermediate step 1992 !-- WS3 as an intermediate step (top) 1972 1993 k = nzt - 1 1973 1994 flux_d = flux_t(k-1) 1974 1995 diss_d = diss_t(k-1) 1975 1996 flux_t(k) = w(k,j,i) * ( & 1976 7. * ( sk(k+1,j,i) + sk(k,j,i) )&1977 -( sk(k+2,j,i) + sk(k-1,j,i) ) ) * adv_sca_31978 diss_t(k) = - abs(w(k,j,i)) * ( &1979 3. * ( sk(k+1,j,i) - sk(k,j,i) )&1980 -( sk(k+2,j,i) - sk(k-1,j,i) ) ) * adv_sca_31997 7.0 * ( sk(k+1,j,i) + sk(k,j,i) ) & 1998 - ( sk(k+2,j,i) + sk(k-1,j,i) ) ) * adv_sca_3 1999 diss_t(k) = - ABS(w(k,j,i)) * ( & 2000 3.0 * ( sk(k+1,j,i) - sk(k,j,i) ) & 2001 - ( sk(k+2,j,i) - sk(k-1,j,i) ) ) * adv_sca_3 1981 2002 1982 2003 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 1983 - ( flux_d + diss_d )) * ddzw(k)1984 ! 1985 !-- 2nd order scheme 2004 - flux_d - diss_d ) * ddzw(k) 2005 ! 2006 !-- 2nd order scheme (top) 1986 2007 k = nzt 1987 2008 flux_d = flux_t(k-1) … … 1993 2014 1994 2015 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 1995 - ( flux_d + diss_d )) * ddzw(k)2016 - flux_d - diss_d ) * ddzw(k) 1996 2017 ! 1997 2018 !-- evaluation of statistics … … 2000 2021 CASE ( 'pt' ) 2001 2022 DO k = nzb_s_inner(j,i), nzt 2002 sums_wspts_ws_l(k,:) = sums_wspts_ws_l(k,:) 2003 + ( flux_t(k) + diss_t(k) ) 2004 * weight_substep(intermediate_timestep_count)&2005 * rmask(j,i,:)2023 sums_wspts_ws_l(k,:) = sums_wspts_ws_l(k,:) & 2024 + ( flux_t(k) + diss_t(k) ) & 2025 * weight_substep(intermediate_timestep_count) & 2026 * rmask(j,i,:) 2006 2027 ENDDO 2007 2028 CASE ( 'sa' ) 2008 2029 DO k = nzb_s_inner(j,i), nzt 2009 sums_wssas_ws_l(k,:) = sums_wssas_ws_l(k,:) 2010 + ( flux_t(k) + diss_t(k) ) 2011 * weight_substep(intermediate_timestep_count)&2012 * rmask(j,i,:)2030 sums_wssas_ws_l(k,:) = sums_wssas_ws_l(k,:) & 2031 + ( flux_t(k) + diss_t(k) ) & 2032 * weight_substep(intermediate_timestep_count) & 2033 * rmask(j,i,:) 2013 2034 ENDDO 2014 2035 CASE ( 'q' ) 2015 2036 DO k = nzb_s_inner(j,i), nzt 2016 sums_wsqs_ws_l(k,:) = sums_wsqs_ws_l(k,:) 2017 + ( flux_t(k) + diss_t(k) ) 2018 * weight_substep(intermediate_timestep_count)&2019 * rmask(j,i,:)2037 sums_wsqs_ws_l(k,:) = sums_wsqs_ws_l(k,:) & 2038 + ( flux_t(k) + diss_t(k) ) & 2039 * weight_substep(intermediate_timestep_count) & 2040 * rmask(j,i,:) 2020 2041 ENDDO 2021 2042 … … 2028 2049 2029 2050 2030 ! 2051 !------------------------------------------------------------------------------! 2031 2052 ! Advection of u - Call for all grid points 2032 2053 !------------------------------------------------------------------------------! … … 2059 2080 IF( boundary_flags(j,i) == 5 ) THEN 2060 2081 2061 DO k = nzb_u_inner(j,i) +1, nzt2082 DO k = nzb_u_inner(j,i)+1, nzt 2062 2083 u_comp(k) = u(k,j,i) + u(k,j,i-1) - gu 2063 swap_flux_x_local_u(k,j) = u_comp(k) * &2084 swap_flux_x_local_u(k,j) = u_comp(k) * & 2064 2085 ( u(k,j,i) + u(k,j,i-1) ) * 0.25 2065 swap_diss_x_local_u(k,j) = diss_2nd( u(k,j,i+2), u(k,j,i+1), &2066 u(k,j,i), u(k,j,i-1), &2067 u(k,j,i-1), u_comp(k), &2086 swap_diss_x_local_u(k,j) = diss_2nd( u(k,j,i+2), u(k,j,i+1), & 2087 u(k,j,i), u(k,j,i-1), & 2088 u(k,j,i-1), u_comp(k), & 2068 2089 0.25, ddx ) 2069 2090 ENDDO … … 2071 2092 ELSE 2072 2093 2073 DO k = nzb_u_inner(j,i) +1, nzt2094 DO k = nzb_u_inner(j,i)+1, nzt 2074 2095 u_comp(k) = u(k,j,i) + u(k,j,i-1) - gu 2075 swap_flux_x_local_u(k,j) = u_comp(k) * ( &2076 37. * ( u(k,j,i) + u(k,j,i-1) )&2077 - 8. * ( u(k,j,i+1) + u(k,j,i-2) )&2078 + (u(k,j,i+2)+u(k,j,i-3) ) )&2096 swap_flux_x_local_u(k,j) = u_comp(k) * ( & 2097 37.0 * ( u(k,j,i) + u(k,j,i-1) ) & 2098 - 8.0 * ( u(k,j,i+1) + u(k,j,i-2) ) & 2099 + (u(k,j,i+2)+u(k,j,i-3) ) ) & 2079 2100 * adv_mom_5 2080 swap_diss_x_local_u(k,j) = - abs(u_comp(k)) * (&2081 10. * ( u(k,j,i) - u(k,j,i-1) )&2082 - 5. * ( u(k,j,i+1) - u(k,j,i-2) )&2083 + ( u(k,j,i+2) - u(k,j,i-3) ) )&2101 swap_diss_x_local_u(k,j) = - ABS(u_comp(k)) * ( & 2102 10.0 * ( u(k,j,i) - u(k,j,i-1) ) & 2103 - 5.0 * ( u(k,j,i+1) - u(k,j,i-2) ) & 2104 + ( u(k,j,i+2) - u(k,j,i-3) ) )& 2084 2105 * adv_mom_5 2085 2106 ENDDO … … 2096 2117 ! 2097 2118 !-- Compute southside fluxes for the south boundary of PE domain 2098 DO k = nzb_u_inner(j,i) +1, nzt2119 DO k = nzb_u_inner(j,i)+1, nzt 2099 2120 v_comp = v(k,j,i) + v(k,j,i-1) - gv 2100 2121 swap_flux_y_local_u(k) = v_comp * & … … 2108 2129 ELSE 2109 2130 2110 DO k = nzb_u_inner(j,i) +1, nzt2131 DO k = nzb_u_inner(j,i)+1, nzt 2111 2132 v_comp = v(k,j,i) + v(k,j,i-1) - gv 2112 swap_flux_y_local_u(k) = v_comp * ( &2113 37. * ( u(k,j,i) + u(k,j-1,i) )&2114 - 8. * ( u(k,j+1,i) + u(k,j-2,i) )&2115 + ( u(k,j+2,i) + u(k,j-3,i) ) )&2133 swap_flux_y_local_u(k) = v_comp * ( & 2134 37.0 * ( u(k,j,i) + u(k,j-1,i) ) & 2135 - 8.0 * ( u(k,j+1,i) + u(k,j-2,i) ) & 2136 + ( u(k,j+2,i) + u(k,j-3,i) ) ) & 2116 2137 * adv_mom_5 2117 swap_diss_y_local_u(k) = - abs(v_comp) * (&2118 10. * ( u(k,j,i) - u(k,j-1,i) )&2119 - 5. * ( u(k,j+1,i) - u(k,j-2,i) )&2120 + ( u(k,j+2,i) - u(k,j-3,i) ) )&2138 swap_diss_y_local_u(k) = - ABS( v_comp ) * ( & 2139 10.0 * ( u(k,j,i) - u(k,j-1,i) ) & 2140 - 5.0 * ( u(k,j+1,i) - u(k,j-2,i) ) & 2141 + ( u(k,j+2,i) - u(k,j-3,i) ) ) & 2121 2142 * adv_mom_5 2122 2143 ENDDO … … 2132 2153 2133 2154 CASE ( 1 ) 2134 DO k = nzb_u_inner(j,i) +1, nzt2155 DO k = nzb_u_inner(j,i)+1, nzt 2135 2156 u_comp(k) = u(k,j,i+1) + u(k,j,i) 2136 flux_r(k) = ( u_comp(k) - gu ) * ( & 2137 7. * ( u(k,j,i+1) + u(k,j,i) ) & 2138 - ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3 2139 diss_r(k) = - abs(u_comp(k) - gu) * ( & 2140 3. * ( u(k,j,i+1) - u(k,j,i) ) & 2141 - ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3 2157 flux_r(k) = ( u_comp(k) - gu ) * ( & 2158 7.0 * ( u(k,j,i+1) + u(k,j,i) ) & 2159 - ( u(k,j,i+2) + u(k,j,i-1) ) ) & 2160 * adv_mom_3 2161 diss_r(k) = - ABS( u_comp(k) - gu ) * ( & 2162 3.0 * ( u(k,j,i+1) - u(k,j,i) ) & 2163 - ( u(k,j,i+2) - u(k,j,i-1) ) ) & 2164 * adv_mom_3 2142 2165 ENDDO 2143 2166 2144 2167 CASE ( 2 ) 2145 DO k = nzb_u_inner(j,i) +1, nzt2168 DO k = nzb_u_inner(j,i)+1, nzt 2146 2169 u_comp(k) = u(k,j,i+1) + u(k,j,i) 2147 2170 flux_r(k) = ( u_comp(k) - gu ) * & … … 2153 2176 2154 2177 CASE ( 3 ) 2155 DO k = nzb_u_inner(j,i) +1, nzt2178 DO k = nzb_u_inner(j,i)+1, nzt 2156 2179 v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv 2157 flux_n(k) = v_comp * ( & 2158 7. * ( u(k,j+1,i) + u(k,j,i) ) & 2159 - ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3 2160 diss_n(k) = - abs(v_comp) * ( & 2161 3. * ( u(k,j+1,i) - u(k,j,i) ) & 2162 - ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3 2180 flux_n(k) = v_comp * ( & 2181 7.0 * ( u(k,j+1,i) + u(k,j,i) ) & 2182 - ( u(k,j+2,i) + u(k,j-1,i) ) ) & 2183 * adv_mom_3 2184 diss_n(k) = - ABS( v_comp ) * ( & 2185 3.0 * ( u(k,j+1,i) - u(k,j,i) ) & 2186 - ( u(k,j+2,i) - u(k,j-1,i) ) ) & 2187 * adv_mom_3 2163 2188 ENDDO 2164 2189 2165 2190 CASE ( 4 ) 2166 DO k = nzb_u_inner(j,i) +1, nzt2191 DO k = nzb_u_inner(j,i)+1, nzt 2167 2192 v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv 2168 2193 flux_n(k) = v_comp * ( u(k,j+1,i) + u(k,j,i) ) * 0.25 … … 2173 2198 2174 2199 CASE ( 5 ) 2175 DO k = nzb_u_inner(j,i) +1, nzt2200 DO k = nzb_u_inner(j,i)+1, nzt 2176 2201 u_comp(k) = u(k,j,i+1) + u(k,j,i) 2177 flux_r(k) = ( u_comp(k) - gu ) * ( & 2178 7. * ( u(k,j,i+1) + u(k,j,i) ) & 2179 - ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3 2180 diss_r(k) = - abs(u_comp(k) - gu) * ( & 2181 3. * ( u(k,j,i+1) - u(k,j,i) ) & 2182 - ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3 2202 flux_r(k) = ( u_comp(k) - gu ) * ( & 2203 7.0 * ( u(k,j,i+1) + u(k,j,i) ) & 2204 - ( u(k,j,i+2) + u(k,j,i-1) ) ) & 2205 * adv_mom_3 2206 diss_r(k) = - ABS( u_comp(k) - gu ) * ( & 2207 3.0 * ( u(k,j,i+1) - u(k,j,i) ) & 2208 - ( u(k,j,i+2) - u(k,j,i-1) ) ) & 2209 * adv_mom_3 2183 2210 ENDDO 2184 2211 2185 2212 CASE ( 7 ) 2186 DO k = nzb_u_inner(j,i) +1, nzt2213 DO k = nzb_u_inner(j,i)+1, nzt 2187 2214 v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv 2188 flux_n(k) = v_comp * ( & 2189 7. * ( u(k,j+1,i) + u(k,j,i) ) & 2190 - ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3 2191 diss_n(k) = - abs(v_comp) * ( & 2192 3. * ( u(k,j+1,i) - u(k,j,i) ) & 2193 - ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3 2215 flux_n(k) = v_comp * ( & 2216 7.0 * ( u(k,j+1,i) + u(k,j,i) ) & 2217 - ( u(k,j+2,i) + u(k,j-1,i) ) ) & 2218 * adv_mom_3 2219 diss_n(k) = - ABS( v_comp ) * ( & 2220 3.0 * ( u(k,j+1,i) - u(k,j,i) ) & 2221 - ( u(k,j+2,i) - u(k,j-1,i) ) ) & 2222 * adv_mom_3 2194 2223 ENDDO 2195 2224 2196 2225 CASE ( 8 ) 2197 DO k = nzb_u_inner(j,i) +1, nzt2226 DO k = nzb_u_inner(j,i)+1, nzt 2198 2227 v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv 2199 flux_n(k) = v_comp * ( & 2200 7. * ( u(k,j+1,i) + u(k,j,i) ) & 2201 - ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3 2202 diss_n(k) = - abs(v_comp) * ( & 2203 3. * ( u(k,j+1,i) - u(k,j,i) ) & 2204 - ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3 2228 flux_n(k) = v_comp * ( & 2229 7.0 * ( u(k,j+1,i) + u(k,j,i) ) & 2230 - ( u(k,j+2,i) + u(k,j-1,i) ) ) & 2231 * adv_mom_3 2232 diss_n(k) = - ABS( v_comp ) * ( & 2233 3.0 * ( u(k,j+1,i) - u(k,j,i) ) & 2234 - ( u(k,j+2,i) - u(k,j-1,i) ) ) & 2235 * adv_mom_3 2205 2236 ENDDO 2206 2237 … … 2210 2241 ! 2211 2242 !-- Compute the crosswise 5th order fluxes at the outflow 2212 IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2&2213 .OR. boundary_flags(j,i) == 5) THEN2243 IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 .OR.& 2244 boundary_flags(j,i) == 5 ) THEN 2214 2245 2215 DO k = nzb_u_inner(j,i) +1, nzt2246 DO k = nzb_u_inner(j,i)+1, nzt 2216 2247 v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv 2217 flux_n(k) = v_comp * ( & 2218 37. * ( u(k,j+1,i) + u(k,j,i) ) & 2219 - 8. * ( u(k,j+2,i) +u(k,j-1,i) ) & 2220 + ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5 2221 diss_n(k) = - abs(v_comp) * ( & 2222 10. * ( u(k,j+1,i) - u(k,j,i) ) & 2223 - 5. * ( u(k,j+2,i) - u(k,j-1,i) ) & 2224 + ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5 2248 flux_n(k) = v_comp * ( & 2249 37.0 * ( u(k,j+1,i) + u(k,j,i) ) & 2250 - 8.0 * ( u(k,j+2,i) +u(k,j-1,i) ) & 2251 + ( u(k,j+3,i) + u(k,j-2,i) ) ) & 2252 * adv_mom_5 2253 diss_n(k) = - ABS( v_comp ) * ( & 2254 10.0 * ( u(k,j+1,i) - u(k,j,i) ) & 2255 - 5.0 * ( u(k,j+2,i) - u(k,j-1,i) ) & 2256 + ( u(k,j+3,i) - u(k,j-2,i) ) ) & 2257 * adv_mom_5 2225 2258 ENDDO 2226 2259 2227 2260 ELSE 2228 2261 2229 DO k = nzb_u_inner(j,i) +1, nzt2262 DO k = nzb_u_inner(j,i)+1, nzt 2230 2263 u_comp(k) = u(k,j,i+1) + u(k,j,i) 2231 flux_r(k) = ( u_comp(k) - gu ) * ( & 2232 37. * ( u(k,j,i+1) + u(k,j,i) ) & 2233 - 8. * ( u(k,j,i+2) + u(k,j,i-1) ) & 2234 + ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5 2235 diss_r(k) = - abs(u_comp(k) - gu) * ( & 2236 10. * ( u(k,j,i+1) - u(k,j,i) ) & 2237 - 5. * ( u(k,j,i+2) - u(k,j,i-1) ) & 2238 + ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5 2264 flux_r(k) = ( u_comp(k) - gu ) * ( & 2265 37.0 * ( u(k,j,i+1) + u(k,j,i) ) & 2266 - 8.0 * ( u(k,j,i+2) + u(k,j,i-1) ) & 2267 + ( u(k,j,i+3) + u(k,j,i-2) ) ) & 2268 * adv_mom_5 2269 diss_r(k) = - ABS( u_comp(k) - gu ) * ( & 2270 10.0 * ( u(k,j,i+1) - u(k,j,i) ) & 2271 - 5.0 * ( u(k,j,i+2) - u(k,j,i-1) ) & 2272 + ( u(k,j,i+3) - u(k,j,i-2) ) ) & 2273 * adv_mom_5 2239 2274 ENDDO 2240 2275 … … 2243 2278 ELSE 2244 2279 2245 DO k = nzb_u_inner(j,i) +1, nzt2280 DO k = nzb_u_inner(j,i)+1, nzt 2246 2281 u_comp(k) = u(k,j,i+1) + u(k,j,i) 2247 flux_r(k) = ( u_comp(k) - gu ) * ( & 2248 37. * ( u(k,j,i+1) + u(k,j,i) ) & 2249 - 8. * ( u(k,j,i+2) + u(k,j,i-1) ) & 2250 + ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5 2251 diss_r(k) = - abs(u_comp(k) - gu) * ( & 2252 10. * ( u(k,j,i+1) - u(k,j,i) ) & 2253 - 5. * ( u(k,j,i+2) - u(k,j,i-1) ) & 2254 + ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5 2282 flux_r(k) = ( u_comp(k) - gu ) * ( & 2283 37.0 * ( u(k,j,i+1) + u(k,j,i) ) & 2284 - 8.0 * ( u(k,j,i+2) + u(k,j,i-1) ) & 2285 + ( u(k,j,i+3) + u(k,j,i-2) ) ) & 2286 * adv_mom_5 2287 diss_r(k) = - ABS( u_comp(k) - gu ) * ( & 2288 10.0 * ( u(k,j,i+1) - u(k,j,i) ) & 2289 - 5.0 * ( u(k,j,i+2) - u(k,j,i-1) ) & 2290 + ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5 2255 2291 2256 2292 v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv 2257 flux_n(k) = v_comp * ( &2258 37. * ( u(k,j+1,i) + u(k,j,i) )&2259 - 8. * ( u(k,j+2,i) + u(k,j-1,i) )&2260 + ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_52261 diss_n(k) = - abs(v_comp) * (&2262 10. * ( u(k,j+1,i) - u(k,j,i) )&2263 - 5. * ( u(k,j+2,i) - u(k,j-1,i) )&2264 + ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_52293 flux_n(k) = v_comp * ( & 2294 37.0 * ( u(k,j+1,i) + u(k,j,i) ) & 2295 - 8.0 * ( u(k,j+2,i) + u(k,j-1,i) ) & 2296 + ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5 2297 diss_n(k) = - ABS( v_comp ) * ( & 2298 10.0 * ( u(k,j+1,i) - u(k,j,i) ) & 2299 - 5.0 * ( u(k,j+2,i) - u(k,j-1,i) ) & 2300 + ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5 2265 2301 2266 2302 ENDDO … … 2268 2304 ENDIF 2269 2305 2270 DO k = nzb_u_inner(j,i) +1, nzt2306 DO k = nzb_u_inner(j,i)+1, nzt 2271 2307 2272 tend(k,j,i) = tend(k,j,i) - ( &2273 ( flux_r(k) + diss_r(k)&2274 - ( swap_flux_x_local_u(k,j) + swap_diss_x_local_u(k,j) ) ) * ddx&2275 + ( flux_n(k) + diss_n(k) &2276 - ( swap_flux_y_local_u(k) + swap_diss_y_local_u(k) )) * ddy )2308 tend(k,j,i) = tend(k,j,i) - ( & 2309 ( flux_r(k) + diss_r(k) & 2310 - swap_flux_x_local_u(k,j) - swap_diss_x_local_u(k,j) ) * ddx & 2311 + ( flux_n(k) + diss_n(k) & 2312 - swap_flux_y_local_u(k) - swap_diss_y_local_u(k) ) * ddy ) 2277 2313 2278 2314 swap_flux_x_local_u(k,j) = flux_r(k) … … 2281 2317 swap_diss_y_local_u(k) = diss_n(k) 2282 2318 2283 sums_us2_ws_l(k,:) = sums_us2_ws_l(k,:) &2284 + ( flux_r(k) &2285 * ( u_comp(k) - 2. * hom(k,1,1,:) ) &2286 / ( u_comp(k) - gu + 1.0E-20 ) &2287 + diss_r(k)&2288 * abs(u_comp(k) - 2. * hom(k,1,1,:) )&2289 / ( abs(u_comp(k) - gu) + 1.0E-20) )&2290 * weight_substep(intermediate_timestep_count) * rmask(j,i,:)2319 sums_us2_ws_l(k,:) = sums_us2_ws_l(k,:) & 2320 + ( flux_r(k) & 2321 * ( u_comp(k) - 2.0 * hom(k,1,1,:) ) & 2322 / ( u_comp(k) - gu + 1.0E-20 ) & 2323 + diss_r(k) & 2324 * ABS( u_comp(k) - 2.0 * hom(k,1,1,:) ) & 2325 / ( ABS( u_comp(k) - gu) + 1.0E-20) ) & 2326 * weight_substep(intermediate_timestep_count) * rmask(j,i,:) 2291 2327 ENDDO 2292 sums_us2_ws_l(nzb_u_inner(j,i),:) = &2328 sums_us2_ws_l(nzb_u_inner(j,i),:) = & 2293 2329 sums_us2_ws_l(nzb_u_inner(j,i)+1,:) 2294 2330 ENDDO … … 2301 2337 DO i = nxlu, nxr 2302 2338 DO j = nys, nyn 2303 k = nzb_u_inner(j,i) +12339 k = nzb_u_inner(j,i)+1 2304 2340 ! 2305 2341 !-- The fluxes flux_d and diss_d at the surface are 0. Due to static 2306 2342 !-- reasons the top flux at the surface should be 0. 2307 flux_t(nzb_u_inner(j,i)) = 0. 2308 diss_t(nzb_u_inner(j,i)) = 0. 2343 flux_t(nzb_u_inner(j,i)) = 0.0 2344 diss_t(nzb_u_inner(j,i)) = 0.0 2309 2345 flux_d = flux_t(k-1) 2310 2346 diss_d = diss_t(k-1) 2311 2347 ! 2312 !-- 2nd order scheme 2348 !-- 2nd order scheme (bottom) 2313 2349 w_comp = w(k,j,i) + w(k,j,i-1) 2314 2350 flux_t(k) = w_comp * ( u(k+1,j,i) + u(k,j,i) ) * 0.25 2315 diss_t(k) = diss_2nd( u(k+2,j,i), u(k+1,j,i), u(k,j,i), &2316 0. , 0., w_comp, 0.25, ddzw(k) )2351 diss_t(k) = diss_2nd( u(k+2,j,i), u(k+1,j,i), u(k,j,i), & 2352 0.0, 0.0, w_comp, 0.25, ddzw(k) ) 2317 2353 2318 2354 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 2319 - ( flux_d + diss_d )) * ddzw(k)2320 ! 2321 !-- WS3 as an intermediate step 2322 k = nzb_u_inner(j,i) +22355 - flux_d - diss_d ) * ddzw(k) 2356 ! 2357 !-- WS3 as an intermediate step (bottom) 2358 k = nzb_u_inner(j,i)+2 2323 2359 flux_d = flux_t(k-1) 2324 2360 diss_d = diss_t(k-1) 2325 2361 w_comp = w(k,j,i) + w(k,j,i-1) 2326 2362 flux_t(k) = w_comp * ( & 2327 7. * (u(k+1,j,i) + u(k,j,i) )&2328 -( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_32329 diss_t(k) = - abs(w_comp) * (&2330 3. * ( u(k+1,j,i) - u(k,j,i) )&2331 -( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_32363 7.0 * ( u(k+1,j,i) + u(k,j,i) ) & 2364 - ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3 2365 diss_t(k) = - ABS( w_comp ) * ( & 2366 3.0 * ( u(k+1,j,i) - u(k,j,i) ) & 2367 - ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3 2332 2368 2333 2369 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 2334 - ( flux_d + diss_d )) * ddzw(k)2370 - flux_d - diss_d ) * ddzw(k) 2335 2371 ! 2336 2372 !WS5 2337 DO k = nzb_u_inner(j,i) + 3, nzt - 32373 DO k = nzb_u_inner(j,i)+3, nzt-2 2338 2374 2339 2375 flux_d = flux_t(k-1) … … 2341 2377 w_comp = w(k,j,i) + w(k,j,i-1) 2342 2378 flux_t(k) = w_comp * ( & 2343 37. 2344 - 8.* ( u(k+2,j,i) + u(k-1,j,i) ) &2345 +( u(k+3,j,i) + u(k-2,j,i) ) ) * adv_mom_52346 diss_t(k) = - abs(w_comp) * (&2347 10. 2348 - 5.* ( u(k+2,j,i) - u(k-1,j,i) ) &2349 +( u(k+3,j,i) - u(k-2,j,i) ) ) * adv_mom_52379 37.0 * ( u(k+1,j,i) + u(k,j,i) ) & 2380 - 8.0 * ( u(k+2,j,i) + u(k-1,j,i) ) & 2381 + ( u(k+3,j,i) + u(k-2,j,i) ) ) * adv_mom_5 2382 diss_t(k) = - ABS( w_comp ) * ( & 2383 10.0 * ( u(k+1,j,i) - u(k,j,i) ) & 2384 - 5.0 * ( u(k+2,j,i) - u(k-1,j,i) ) & 2385 + ( u(k+3,j,i) - u(k-2,j,i) ) ) * adv_mom_5 2350 2386 2351 2387 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 2352 - ( flux_d + diss_d )) * ddzw(k)2388 - flux_d - diss_d ) * ddzw(k) 2353 2389 2354 2390 ENDDO 2355 2391 ! 2356 !-- WS3 as an intermediate step2357 DO k = nzt - 2, nzt -12358 flux_d= flux_t(k-1)2359 diss_d= diss_t(k-1)2360 w_comp= w(k,j,i) + w(k,j,i-1)2361 flux_t(k) = w_comp * (&2362 7. * ( u(k+1,j,i) + u(k,j,i) )&2363 -( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_32364 diss_t(k) = - abs(w_comp) * (&2365 3. * ( u(k+1,j,i) - u(k,j,i) )&2366 -( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_32392 !-- WS3 as an intermediate step (top) 2393 k = nzt-1 2394 flux_d = flux_t(k-1) 2395 diss_d = diss_t(k-1) 2396 w_comp = w(k,j,i) + w(k,j,i-1) 2397 flux_t(k) = w_comp * ( & 2398 7.0 * ( u(k+1,j,i) + u(k,j,i) ) & 2399 - ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3 2400 diss_t(k) = - ABS( w_comp ) * ( & 2401 3.0 * ( u(k+1,j,i) - u(k,j,i) ) & 2402 - ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3 2367 2403 2368 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 2369 - ( flux_d + diss_d ) ) * ddzw(k) 2370 ENDDO 2371 ! 2372 !-- 2nd order scheme 2404 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 2405 - flux_d - diss_d ) * ddzw(k) 2406 ! 2407 !-- 2nd order scheme (top) 2373 2408 k = nzt 2374 2409 flux_d = flux_t(k-1) … … 2381 2416 2382 2417 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 2383 - ( flux_d + diss_d )) * ddzw(k)2418 - flux_d - diss_d ) * ddzw(k) 2384 2419 ! 2385 2420 !-- at last vertical momentum flux is accumulated 2386 2421 DO k = nzb_u_inner(j,i), nzt 2387 2422 sums_wsus_ws_l(k,:) = sums_wsus_ws_l(k,:) & 2388 + ( flux_t(k) + diss_t(k) )&2389 * weight_substep(intermediate_timestep_count)&2390 *rmask(j,i,:)2423 + ( flux_t(k) + diss_t(k) ) & 2424 * weight_substep(intermediate_timestep_count) & 2425 * rmask(j,i,:) 2391 2426 ENDDO 2392 2427 ENDDO … … 2397 2432 2398 2433 2399 ! 2434 !------------------------------------------------------------------------------! 2400 2435 ! Advection of v - Call for all grid points 2401 2436 !------------------------------------------------------------------------------! … … 2429 2464 2430 2465 IF ( boundary_flags(j,i) == 6 ) THEN 2431 DO k = nzb_v_inner(j,i) +1, nzt2466 DO k = nzb_v_inner(j,i)+1, nzt 2432 2467 u_comp = u(k,j-1,i) + u(k,j,i) - gu 2433 2468 swap_flux_x_local_v(k,j) = u_comp * & … … 2443 2478 DO k = nzb_v_inner(j,i)+1, nzt 2444 2479 u_comp = u(k,j-1,i) + u(k,j,i) - gu 2445 swap_flux_x_local_v(k,j) = u_comp * ( &2446 37. * ( v(k,j,i) + v(k,j,i-1) )&2447 - 8. * ( v(k,j,i+1) + v(k,j,i-2) )&2448 + ( v(k,j,i+2) + v(k,j,i-3) ) )&2480 swap_flux_x_local_v(k,j) = u_comp * ( & 2481 37.0 * ( v(k,j,i) + v(k,j,i-1) ) & 2482 - 8.0 * ( v(k,j,i+1) + v(k,j,i-2) ) & 2483 + ( v(k,j,i+2) + v(k,j,i-3) ) )& 2449 2484 * adv_mom_5 2450 swap_diss_x_local_v(k,j) = - abs(u_comp) * (&2451 10. * ( v(k,j,i) - v(k,j,i-1) )&2452 - 5. * ( v(k,j,i+1) - v(k,j,i-2) )&2453 + ( v(k,j,i+2) - v(k,j,i-3) ) )&2485 swap_diss_x_local_v(k,j) = - ABS( u_comp ) * ( & 2486 10.0 * ( v(k,j,i) - v(k,j,i-1) ) & 2487 - 5.0 * ( v(k,j,i+1) - v(k,j,i-2) ) & 2488 + ( v(k,j,i+2) - v(k,j,i-3) ) )& 2454 2489 * adv_mom_5 2455 2490 ENDDO … … 2463 2498 IF ( boundary_flags(j,i) == 7 ) THEN 2464 2499 2465 DO k = nzb_v_inner(j,i) +1, nzt2500 DO k = nzb_v_inner(j,i)+1, nzt 2466 2501 v_comp(k) = v(k,j,i) + v(k,j-1,i) - gv 2467 2502 swap_flux_y_local_v(k) = v_comp(k) * & … … 2475 2510 ELSE 2476 2511 2477 DO k = nzb_v_inner(j,i) +1, nzt2512 DO k = nzb_v_inner(j,i)+1, nzt 2478 2513 v_comp(k) = v(k,j,i) + v(k,j-1,i) - gv 2479 swap_flux_y_local_v(k) = v_comp(k) * ( &2480 37. * ( v(k,j,i) + v(k,j-1,i) )&2481 - 8. * ( v(k,j+1,i) + v(k,j-2,i) )&2482 + ( v(k,j+2,i) + v(k,j-3,i) ) )&2514 swap_flux_y_local_v(k) = v_comp(k) * ( & 2515 37.0 * ( v(k,j,i) + v(k,j-1,i) ) & 2516 - 8.0 * ( v(k,j+1,i) + v(k,j-2,i) ) & 2517 + ( v(k,j+2,i) + v(k,j-3,i) ) ) & 2483 2518 * adv_mom_5 2484 swap_diss_y_local_v(k) = - abs(v_comp(k)) * (&2485 10. * ( v(k,j,i) - v(k,j-1,i) )&2486 - 5. * ( v(k,j+1,i) - v(k,j-2,i) )&2487 + ( v(k,j+2,i) - v(k,j-3,i) ) )&2519 swap_diss_y_local_v(k) = - ABS( v_comp(k) ) * ( & 2520 10.0 * ( v(k,j,i) - v(k,j-1,i) ) & 2521 - 5.0 * ( v(k,j+1,i) - v(k,j-2,i) ) & 2522 + ( v(k,j+2,i) - v(k,j-3,i) ) ) & 2488 2523 * adv_mom_5 2489 2524 ENDDO … … 2498 2533 2499 2534 CASE ( 1 ) 2500 DO k = nzb_v_inner(j,i) +1, nzt2535 DO k = nzb_v_inner(j,i)+1, nzt 2501 2536 u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu 2502 flux_r(k) = u_comp * ( & 2503 7. * (v(k,j,i+1) + v(k,j,i) ) & 2504 - ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3 2505 diss_r(k) = - abs(u_comp) * ( & 2506 3. * ( v(k,j,i+1) - v(k,j,i) ) & 2507 - ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3 2537 flux_r(k) = u_comp * ( & 2538 7.0 * (v(k,j,i+1) + v(k,j,i) ) & 2539 - ( v(k,j,i+2) + v(k,j,i-1) ) ) & 2540 * adv_mom_3 2541 diss_r(k) = - ABS( u_comp ) * ( & 2542 3.0 * ( v(k,j,i+1) - v(k,j,i) ) & 2543 - ( v(k,j,i+2) - v(k,j,i-1) ) ) & 2544 * adv_mom_3 2508 2545 ENDDO 2509 2546 2510 2547 CASE ( 2 ) 2511 DO k = nzb_v_inner(j,i) +1, nzt2548 DO k = nzb_v_inner(j,i)+1, nzt 2512 2549 u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu 2513 2550 flux_r(k) = u_comp * ( v(k,j,i+1) + v(k,j,i) ) * 0.25 … … 2518 2555 2519 2556 CASE ( 3 ) 2520 DO k = nzb_v_inner(j,i) +1, nzt2557 DO k = nzb_v_inner(j,i)+1, nzt 2521 2558 v_comp(k) = v(k,j+1,i) + v(k,j,i) 2522 2559 flux_n(k) = ( v_comp(k)- gv ) * ( & 2523 7. * ( v(k,j+1,i) + v(k,j,i) ) & 2524 - ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3 2525 diss_n(k) = - abs(v_comp(k) - gv) * ( & 2526 3. * ( v(k,j+1,i) - v(k,j,i) ) & 2527 - ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3 2560 7.0 * ( v(k,j+1,i) + v(k,j,i) ) & 2561 - ( v(k,j+2,i) + v(k,j-1,i) ) ) & 2562 * adv_mom_3 2563 diss_n(k) = - ABS(v_comp(k) - gv) * ( & 2564 3.0 * ( v(k,j+1,i) - v(k,j,i) ) & 2565 - ( v(k,j+2,i) - v(k,j-1,i) ) ) & 2566 * adv_mom_3 2528 2567 ENDDO 2529 2568 2530 2569 CASE ( 4 ) 2531 DO k = nzb_v_inner(j,i) +1, nzt2570 DO k = nzb_v_inner(j,i)+1, nzt 2532 2571 v_comp(k) = v(k,j+1,i) + v(k,j,i) 2533 2572 flux_n(k) = ( v_comp(k) - gv ) * & … … 2539 2578 2540 2579 CASE ( 5 ) 2541 DO k = nzb_w_inner(j,i) +1, nzt2580 DO k = nzb_w_inner(j,i)+1, nzt 2542 2581 u_comp = u(k,j-1,i) + u(k,j,i) - gu 2543 2582 flux_r(k) = u_comp * ( & 2544 7. * ( v(k,j,i+1) + v(k,j,i) ) & 2545 - ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3 2546 diss_r(k) = - abs(u_comp) * ( & 2547 3. * ( w(k,j,i+1) - w(k,j,i) ) & 2548 - ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3 2583 7.0 * ( v(k,j,i+1) + v(k,j,i) ) & 2584 - ( v(k,j,i+2) + v(k,j,i-1) ) ) & 2585 * adv_mom_3 2586 diss_r(k) = - ABS (u_comp ) * ( & 2587 3.0 * ( w(k,j,i+1) - w(k,j,i) ) & 2588 - ( v(k,j,i+2) - v(k,j,i-1) ) ) & 2589 * adv_mom_3 2549 2590 ENDDO 2550 2591 2551 2592 CASE ( 6 ) 2552 DO k = nzb_v_inner(j,i) +1, nzt2593 DO k = nzb_v_inner(j,i)+1, nzt 2553 2594 2554 2595 u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu 2555 2596 flux_r(k) = u_comp * ( & 2556 7. * ( v(k,j,i+1) + v(k,j,i) ) & 2557 - ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3 2558 diss_r(k) = - abs(u_comp) * ( & 2559 3. * ( v(k,j,i+1) - v(k,j,i) ) & 2560 - ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3 2597 7.0 * ( v(k,j,i+1) + v(k,j,i) ) & 2598 - ( v(k,j,i+2) + v(k,j,i-1) ) ) & 2599 * adv_mom_3 2600 diss_r(k) = - ABS( u_comp ) * ( & 2601 3.0 * ( v(k,j,i+1) - v(k,j,i) ) & 2602 - ( v(k,j,i+2) - v(k,j,i-1) ) ) & 2603 * adv_mom_3 2561 2604 ENDDO 2562 2605 2563 2606 CASE ( 7 ) 2564 DO k = nzb_v_inner(j,i) +1, nzt2607 DO k = nzb_v_inner(j,i)+1, nzt 2565 2608 v_comp(k) = v(k,j+1,i) + v(k,j,i) 2566 2609 flux_n(k) = ( v_comp(k) - gv ) * ( & 2567 7. * ( v(k,j+1,i) + v(k,j,i) ) & 2568 - ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3 2569 diss_n(k) = - abs(v_comp(k) - gv) * ( & 2570 3. * ( v(k,j+1,i) - v(k,j,i) ) & 2571 - ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3 2610 7.0 * ( v(k,j+1,i) + v(k,j,i) ) & 2611 - ( v(k,j+2,i) + v(k,j-1,i) ) ) & 2612 * adv_mom_3 2613 diss_n(k) = - ABS( v_comp(k) - gv ) * ( & 2614 3.0 * ( v(k,j+1,i) - v(k,j,i) ) & 2615 - ( v(k,j+2,i) - v(k,j-1,i) ) ) & 2616 * adv_mom_3 2572 2617 ENDDO 2573 2618 … … 2576 2621 END SELECT 2577 2622 ! 2578 !-- Compute the crosswise 5th order fluxes at the outflow 2579 IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 & 2580 .OR. boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 ) THEN 2623 !-- Compute the crosswise 5th order fluxes at the outflow 2624 IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 .OR.& 2625 boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 ) & 2626 THEN 2581 2627 2582 DO k = nzb_v_inner(j,i) +1, nzt2628 DO k = nzb_v_inner(j,i)+1, nzt 2583 2629 v_comp(k) = v(k,j+1,i) + v(k,j,i) 2584 2630 flux_n(k) = ( v_comp(k) - gv ) * ( & 2585 37. * ( v(k,j+1,i) + v(k,j,i) ) & 2586 - 8. * ( v(k,j+2,i) + v(k,j-1,i) ) & 2587 + ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5 2588 diss_n(k) = - abs(v_comp(k) - gv ) * ( & 2589 10. * ( v(k,j+1,i) - v(k,j,i) ) & 2590 - 5. * ( v(k,j+2,i) - v(k,j-1,i) ) & 2591 + ( v(k,j+3,i) - v(k,j-2,i) ) ) *adv_mom_5 2631 37.0 * ( v(k,j+1,i) + v(k,j,i) ) & 2632 - 8.0 * ( v(k,j+2,i) + v(k,j-1,i) ) & 2633 + ( v(k,j+3,i) + v(k,j-2,i) ) ) & 2634 * adv_mom_5 2635 diss_n(k) = - ABS( v_comp(k) - gv ) * ( & 2636 10.0 * ( v(k,j+1,i) - v(k,j,i) ) & 2637 - 5.0 * ( v(k,j+2,i) - v(k,j-1,i) ) & 2638 + ( v(k,j+3,i) - v(k,j-2,i) ) ) & 2639 * adv_mom_5 2592 2640 ENDDO 2593 2641 2594 2642 ELSE 2595 2643 2596 DO k = nzb_v_inner(j,i) +1, nzt2644 DO k = nzb_v_inner(j,i)+1, nzt 2597 2645 u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu 2598 2646 flux_r(k) = u_comp * ( & 2599 37. * ( v(k,j,i+1) + v(k,j,i) ) & 2600 - 8. * ( v(k,j,i+2) + v(k,j,i-1) ) & 2601 + ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5 2602 diss_r(k) = - abs(u_comp) * ( & 2603 10. * ( v(k,j,i+1) - v(k,j,i) ) & 2604 -5. * ( v(k,j,i+2) - v(k,j,i-1) ) & 2605 + ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5 2647 37.0 * ( v(k,j,i+1) + v(k,j,i) ) & 2648 - 8.0 * ( v(k,j,i+2) + v(k,j,i-1) ) & 2649 + ( v(k,j,i+3) + v(k,j,i-2) ) ) & 2650 * adv_mom_5 2651 diss_r(k) = - ABS( u_comp ) * ( & 2652 10.0 * ( v(k,j,i+1) - v(k,j,i) ) & 2653 - 5.0 * ( v(k,j,i+2) - v(k,j,i-1) ) & 2654 + ( v(k,j,i+3) - v(k,j,i-2) ) ) & 2655 * adv_mom_5 2606 2656 ENDDO 2607 2657 … … 2611 2661 ELSE 2612 2662 2613 DO k = nzb_v_inner(j,i) +1, nzt2663 DO k = nzb_v_inner(j,i)+1, nzt 2614 2664 u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu 2615 2665 flux_r(k) = u_comp * ( & 2616 37. * ( v(k,j,i+1) + v(k,j,i) ) & 2617 - 8. * ( v(k,j,i+2) + v(k,j,i-1) ) & 2618 + ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5 2619 diss_r(k) = - abs(u_comp) * ( & 2620 10. * ( v(k,j,i+1) - v(k,j,i) ) & 2621 -5. * ( v(k,j,i+2) - v(k,j,i-1) ) & 2622 + ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5 2666 37.0 * ( v(k,j,i+1) + v(k,j,i) ) & 2667 - 8.0 * ( v(k,j,i+2) + v(k,j,i-1) ) & 2668 + ( v(k,j,i+3) + v(k,j,i-2) ) ) & 2669 * adv_mom_5 2670 diss_r(k) = - ABS ( u_comp ) * ( & 2671 10.0 * ( v(k,j,i+1) - v(k,j,i) ) & 2672 - 5.0 * ( v(k,j,i+2) - v(k,j,i-1) ) & 2673 + ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5 2623 2674 2624 2675 v_comp(k) = v(k,j+1,i) + v(k,j,i) 2625 2676 flux_n(k) = ( v_comp(k) - gv ) * ( & 2626 37. * ( v(k,j+1,i) + v(k,j,i) )&2627 - 8. * ( v(k,j+2,i) + v(k,j-1,i) )&2628 + ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_52629 diss_n(k) = - abs(v_comp(k) - gv ) * (&2630 10. * ( v(k,j+1,i) - v(k,j,i) )&2631 - 5. * ( v(k,j+2,i) - v(k,j-1,i) )&2632 + ( v(k,j+3,i) - v(k,j-2,i) ) ) *adv_mom_52677 37.0 * ( v(k,j+1,i) + v(k,j,i) ) & 2678 - 8.0 * ( v(k,j+2,i) + v(k,j-1,i) ) & 2679 + ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5 2680 diss_n(k) = - ABS( v_comp(k) - gv ) * ( & 2681 10.0 * ( v(k,j+1,i) - v(k,j,i) ) & 2682 - 5.0 * ( v(k,j+2,i) - v(k,j-1,i) ) & 2683 + ( v(k,j+3,i) - v(k,j-2,i) ) ) *adv_mom_5 2633 2684 2634 2685 ENDDO 2635 2686 ENDIF 2636 2687 2637 DO k = nzb_v_inner(j,i) +1, nzt2688 DO k = nzb_v_inner(j,i)+1, nzt 2638 2689 tend(k,j,i) = tend(k,j,i) - ( & 2639 2690 ( flux_r(k) + diss_r(k) & 2640 - ( swap_flux_x_local_v(k,j) + swap_diss_x_local_v(k,j) )) * ddx &2691 - swap_flux_x_local_v(k,j) - swap_diss_x_local_v(k,j) ) * ddx & 2641 2692 + ( flux_n(k) + diss_n(k) & 2642 - ( swap_flux_y_local_v(k) + swap_diss_y_local_v(k) )) * ddy )2693 - swap_flux_y_local_v(k) - swap_diss_y_local_v(k) ) * ddy ) 2643 2694 2644 2695 swap_flux_x_local_v(k,j) = flux_r(k) … … 2647 2698 swap_diss_y_local_v(k) = diss_n(k) 2648 2699 2649 sums_vs2_ws_l(k,:) = sums_vs2_ws_l(k,:) &2650 + ( flux_n(k) *( v_comp(k) - 2. * hom(k,1,2,:) )&2651 / ( v_comp(k) - gv + 1.0E-20 ) &2652 + diss_n(k) * abs(v_comp(k) - 2. * hom(k,1,2,:) )&2653 / (abs(v_comp(k) - gv) + 1.0E-20 ) )&2654 * weight_substep(intermediate_timestep_count) * rmask(j,i,:)2700 sums_vs2_ws_l(k,:) = sums_vs2_ws_l(k,:) & 2701 + ( flux_n(k) * ( v_comp(k) - 2.0 * hom(k,1,2,:) ) & 2702 / ( v_comp(k) - gv + 1.0E-20 ) & 2703 + diss_n(k) * ABS( v_comp(k) - 2.0 * hom(k,1,2,:) ) & 2704 / ( ABS( v_comp(k) - gv ) + 1.0E-20 ) ) & 2705 * weight_substep(intermediate_timestep_count) * rmask(j,i,:) 2655 2706 ENDDO 2656 sums_vs2_ws_l(nzb_v_inner(j,i),:) = &2707 sums_vs2_ws_l(nzb_v_inner(j,i),:) = & 2657 2708 sums_vs2_ws_l(nzb_v_inner(j,i)+1,:) 2658 2709 ENDDO … … 2667 2718 !-- The fluxes flux_d and diss_d at the surface are 0. Due to static 2668 2719 !-- reasons the top flux at the surface should be 0. 2669 flux_t(nzb_v_inner(j,i)) = 0. 2670 diss_t(nzb_v_inner(j,i)) = 0. 2671 k = nzb_v_inner(j,i) +12720 flux_t(nzb_v_inner(j,i)) = 0.0 2721 diss_t(nzb_v_inner(j,i)) = 0.0 2722 k = nzb_v_inner(j,i)+1 2672 2723 flux_d = flux_t(k-1) 2673 2724 diss_d = diss_t(k-1) 2674 2725 ! 2675 !-- 2nd order scheme 2726 !-- 2nd order scheme (bottom) 2676 2727 w_comp = w(k,j-1,i) + w(k,j,i) 2677 2728 flux_t(k) = w_comp * ( v(k+1,j,i) + v(k,j,i) ) * 0.25 2678 2729 diss_t(k) = diss_2nd( v(k+2,j,i), v(k+1,j,i), v(k,j,i), & 2679 0. , 0., w_comp, 0.25, ddzw(k) )2730 0.0, 0.0, w_comp, 0.25, ddzw(k) ) 2680 2731 2681 2732 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 2682 - ( flux_d + diss_d )) * ddzw(k)2683 ! 2684 !-- WS3 as an intermediate step2733 - flux_d - diss_d ) * ddzw(k) 2734 ! 2735 !-- WS3 as an intermediate step (bottom) 2685 2736 k = nzb_v_inner(j,i) + 2 2686 2737 flux_d = flux_t(k-1) … … 2688 2739 w_comp = w(k,j-1,i) + w(k,j,i) 2689 2740 flux_t(k) = w_comp * ( & 2690 7. * ( v(k+1,j,i) + v(k,j,i) )&2691 -( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_32692 diss_t(k) = - abs(w_comp) * (&2693 3. * ( v(k+1,j,i) - v(k,j,i) )&2694 -( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_32741 7.0 * ( v(k+1,j,i) + v(k,j,i) ) & 2742 - ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3 2743 diss_t(k) = - ABS( w_comp ) * ( & 2744 3.0 * ( v(k+1,j,i) - v(k,j,i) ) & 2745 - ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3 2695 2746 2696 2747 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 2697 - ( flux_d + diss_d )) * ddzw(k)2748 - flux_d - diss_d ) * ddzw(k) 2698 2749 2699 2750 ! 2700 2751 !-- WS5 2701 DO k = nzb_v_inner(j,i) + 3, nzt - 32752 DO k = nzb_v_inner(j,i)+3, nzt-2 2702 2753 flux_d = flux_t(k-1) 2703 2754 diss_d = diss_t(k-1) 2704 2755 w_comp = w(k,j-1,i) + w(k,j,i) 2705 2756 flux_t(k) = w_comp * ( & 2706 37. 2707 - 8.* ( v(k+2,j,i) + v(k-1,j,i) ) &2708 +( v(k+3,j,i) + v(k-2,j,i) ) ) * adv_mom_52709 diss_t(k) = - abs(w_comp) * (&2710 10. * ( v(k+1,j,i) - v(k,j,i) )&2711 -5. * ( v(k+2,j,i) - v(k-1,j,i) )&2712 +( v(k+3,j,i) - v(k-2,j,i) ) ) * adv_mom_52757 37.0 * ( v(k+1,j,i) + v(k,j,i) ) & 2758 - 8.0 * ( v(k+2,j,i) + v(k-1,j,i) ) & 2759 + ( v(k+3,j,i) + v(k-2,j,i) ) ) * adv_mom_5 2760 diss_t(k) = - ABS( w_comp ) * ( & 2761 10.0 * ( v(k+1,j,i) - v(k,j,i) ) & 2762 - 5.0 * ( v(k+2,j,i) - v(k-1,j,i) ) & 2763 + ( v(k+3,j,i) - v(k-2,j,i) ) ) * adv_mom_5 2713 2764 2714 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) 2715 - ( flux_d + diss_d )) * ddzw(k)2765 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 2766 - flux_d - diss_d ) * ddzw(k) 2716 2767 ENDDO 2717 2768 ! 2718 !-- WS3 as an intermediate step 2719 DO k = nzt - 2, nzt -12720 flux_d= flux_t(k-1)2721 diss_d= diss_t(k-1)2722 w_comp= w(k,j-1,i) + w(k,j,i)2723 flux_t(k) = w_comp * (&2724 7. * ( v(k+1,j,i) + v(k,j,i)) &2725 -( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_32726 diss_t(k) = - abs(w_comp) * (&2727 3. * ( v(k+1,j,i) - v(k,j,i)) &2728 -( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_32769 !-- WS3 as an intermediate step (top) 2770 k = nzt-1 2771 flux_d = flux_t(k-1) 2772 diss_d = diss_t(k-1) 2773 w_comp = w(k,j-1,i) + w(k,j,i) 2774 flux_t(k) = w_comp * ( & 2775 7.0 * ( v(k+1,j,i) + v(k,j,i) ) & 2776 - ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3 2777 diss_t(k) = - ABS( w_comp ) * ( & 2778 3.0 * ( v(k+1,j,i) - v(k,j,i) ) & 2779 - ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3 2729 2780 2730 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 2731 - ( flux_d + diss_d ) ) * ddzw(k) 2732 ENDDO 2733 ! 2734 !-- 2nd order scheme 2781 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 2782 - flux_d - diss_d ) * ddzw(k) 2783 ! 2784 !-- 2nd order scheme (top) 2735 2785 k = nzt 2736 2786 flux_d = flux_t(k-1) … … 2738 2788 w_comp = w(k,j-1,i) + w(k,j,i) 2739 2789 flux_t(k) = w_comp * ( v(k+1,j,i) + v(k,j,i) ) * 0.25 2740 diss_t(k) = diss_2nd( v(nzt+1,j,i), v(nzt+1,j,i), v(k,j,i), &2741 v(k-1,j,i), v(k-2,j,i), w_comp, 0.25, ddzw(k))2742 2743 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) &2744 - ( flux_d + diss_d )) * ddzw(k)2745 ! 2746 !- at last vertical momentum flux is accumulated2790 diss_t(k) = diss_2nd( v(nzt+1,j,i), v(nzt+1,j,i), v(k,j,i), & 2791 v(k-1,j,i), v(k-2,j,i), w_comp, 0.25, ddzw(k)) 2792 2793 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 2794 - flux_d - diss_d ) * ddzw(k) 2795 ! 2796 !- At last vertical momentum flux is accumulated. 2747 2797 DO k = nzb_v_inner(j,i), nzt 2748 2798 sums_wsvs_ws_l(k,:) = sums_wsvs_ws_l(k,:) & 2749 + ( flux_t(k) + diss_t(k) )&2750 * weight_substep(intermediate_timestep_count)&2751 *rmask(j,i,:)2799 + ( flux_t(k) + diss_t(k) ) & 2800 * weight_substep(intermediate_timestep_count) & 2801 * rmask(j,i,:) 2752 2802 ENDDO 2753 2803 sums_vs2_ws_l(nzb_v_inner(j,i),:) = & 2754 sums_vs2_ws_l(nzb_v_inner(j,i)+1,:)2804 sums_vs2_ws_l(nzb_v_inner(j,i)+1,:) 2755 2805 ENDDO 2756 2806 ENDDO … … 2759 2809 2760 2810 2761 ! 2811 !------------------------------------------------------------------------------! 2762 2812 ! Advection of w - Call for all grid points 2763 2813 !------------------------------------------------------------------------------! … … 2790 2840 IF ( boundary_flags(j,i) == 6 ) THEN 2791 2841 2792 DO k = nzb_v_inner(j,i) +1, nzt2842 DO k = nzb_v_inner(j,i)+1, nzt 2793 2843 u_comp = u(k+1,j,i) + u(k,j,i) - gu 2794 2844 swap_flux_x_local_w(k,j) = u_comp * & … … 2802 2852 ELSE 2803 2853 2804 DO k = nzb_v_inner(j,i) +1, nzt2854 DO k = nzb_v_inner(j,i)+1, nzt 2805 2855 u_comp = u(k+1,j,i) + u(k,j,i) - gu 2806 swap_flux_x_local_w(k,j) = u_comp * ( &2807 37. * ( w(k,j,i) + w(k,j,i-1) )&2808 - 8. * ( w(k,j,i+1) + w(k,j,i-2) )&2809 + ( w(k,j,i+2) + w(k,j,i-3) ) )&2856 swap_flux_x_local_w(k,j) = u_comp * ( & 2857 37.0 * ( w(k,j,i) + w(k,j,i-1) ) & 2858 - 8.0 * ( w(k,j,i+1) + w(k,j,i-2) ) & 2859 + ( w(k,j,i+2) + w(k,j,i-3) ) )& 2810 2860 * adv_mom_5 2811 swap_diss_x_local_w(k,j) = - abs(u_comp) * (&2812 10. * ( w(k,j,i) - w(k,j,i-1) )&2813 - 5. * ( w(k,j,i+1) - w(k,j,i-2) )&2814 + ( w(k,j,i+2) - w(k,j,i-3) ) )&2861 swap_diss_x_local_w(k,j) = - ABS( u_comp ) * ( & 2862 10.0 * ( w(k,j,i) - w(k,j,i-1) ) & 2863 - 5.0 * ( w(k,j,i+1) - w(k,j,i-2) ) & 2864 + ( w(k,j,i+2) - w(k,j,i-3) ) )& 2815 2865 * adv_mom_5 2816 2866 ENDDO … … 2824 2874 IF ( boundary_flags(j,i) == 8 ) THEN 2825 2875 2826 DO k = nzb_v_inner(j,i) +1, nzt2876 DO k = nzb_v_inner(j,i)+1, nzt 2827 2877 v_comp = v(k+1,j,i) + v(k,j,i) - gv 2828 swap_flux_y_local_w(k) = v_comp * &2878 swap_flux_y_local_w(k) = v_comp * & 2829 2879 ( w(k,j,i) + w(k,j-1,i) ) * 0.25 2830 swap_diss_y_local_w(k) = diss_2nd( w(k,j+2,i), w(k,j+1,i), &2831 w(k,j,i), w(k,j-1,i), &2832 w(k,j-1,i), v_comp, 0.25, ddy)2880 swap_diss_y_local_w(k) = diss_2nd( w(k,j+2,i), w(k,j+1,i), & 2881 w(k,j,i), w(k,j-1,i), & 2882 w(k,j-1,i), v_comp, 0.25, ddy) 2833 2883 ENDDO 2834 2884 2835 2885 ELSE 2836 2886 2837 DO k = nzb_v_inner(j,i) +1, nzt2887 DO k = nzb_v_inner(j,i)+1, nzt 2838 2888 v_comp = v(k+1,j,i) + v(k,j,i) - gv 2839 swap_flux_y_local_w(k) = v_comp * ( &2840 37. * ( w(k,j,i) + w(k,j-1,i) )&2841 - 8. * ( w(k,j+1,i) +w(k,j-2,i) )&2842 + ( w(k,j+2,i) +w(k,j-3,i) ) )&2889 swap_flux_y_local_w(k) = v_comp * ( & 2890 37.0 * ( w(k,j,i) + w(k,j-1,i) ) & 2891 - 8.0 * ( w(k,j+1,i) +w(k,j-2,i) ) & 2892 + ( w(k,j+2,i) +w(k,j-3,i) ) ) & 2843 2893 * adv_mom_5 2844 swap_diss_y_local_w(k) = - abs(v_comp) * (&2845 10. * ( w(k,j,i) - w(k,j-1,i) )&2846 - 5. * ( w(k,j+1,i) - w(k,j-2,i) )&2847 + ( w(k,j+2,i) - w(k,j-3,i) ) )&2894 swap_diss_y_local_w(k) = - ABS( v_comp ) * ( & 2895 10.0 * ( w(k,j,i) - w(k,j-1,i) ) & 2896 - 5.0 * ( w(k,j+1,i) - w(k,j-2,i) ) & 2897 + ( w(k,j+2,i) - w(k,j-3,i) ) ) & 2848 2898 * adv_mom_5 2849 2899 ENDDO … … 2859 2909 2860 2910 CASE ( 1 ) 2861 DO k = nzb_w_inner(j,i) +1, nzt2911 DO k = nzb_w_inner(j,i)+1, nzt 2862 2912 u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu 2863 2913 flux_r(k) = u_comp * ( & 2864 7. * ( w(k,j,i+1) + w(k,j,i) ) & 2865 - ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3 2866 diss_r(k) = -abs(u_comp) * ( & 2867 3. * ( w(k,j,i+1) - w(k,j,i) ) & 2868 - ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3 2914 7.0 * ( w(k,j,i+1) + w(k,j,i) ) & 2915 - ( w(k,j,i+2) + w(k,j,i-1) ) ) & 2916 * adv_mom_3 2917 diss_r(k) = -ABS( u_comp ) * ( & 2918 3.0 * ( w(k,j,i+1) - w(k,j,i) ) & 2919 - ( w(k,j,i+2) - w(k,j,i-1) ) ) & 2920 * adv_mom_3 2869 2921 ENDDO 2870 2922 2871 2923 CASE ( 2 ) 2872 DO k = nzb_w_inner(j,i) +1, nzt2924 DO k = nzb_w_inner(j,i)+1, nzt 2873 2925 u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu 2874 2926 flux_r(k) = u_comp * ( w(k,j,i+1) + w(k,j,i) ) * 0.25 … … 2879 2931 2880 2932 CASE ( 3 ) 2881 DO k = nzb_w_inner(j,i) +1, nzt2933 DO k = nzb_w_inner(j,i)+1, nzt 2882 2934 v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv 2883 2935 flux_n(k) = v_comp * ( & 2884 7. * ( w(k,j+1,i) + w(k,j,i) ) & 2885 - ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3 2886 diss_n(k) = -abs(v_comp) * ( & 2887 3. * ( w(k,j+1,i) - w(k,j,i) ) & 2888 - ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3 2936 7.0 * ( w(k,j+1,i) + w(k,j,i) ) & 2937 - ( w(k,j+2,i) + w(k,j-1,i) ) ) & 2938 * adv_mom_3 2939 diss_n(k) = -ABS( v_comp ) * ( & 2940 3.0 * ( w(k,j+1,i) - w(k,j,i) ) & 2941 - ( w(k,j+2,i) - w(k,j-1,i) ) ) & 2942 * adv_mom_3 2889 2943 ENDDO 2890 2944 2891 2945 CASE ( 4 ) 2892 DO k = nzb_w_inner(j,i) +1, nzt2946 DO k = nzb_w_inner(j,i)+1, nzt 2893 2947 v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv 2894 2948 flux_n(k) = v_comp * ( w(k,j+1,i) + w(k,j,i) ) * 0.25 … … 2899 2953 2900 2954 CASE ( 5 ) 2901 DO k = nzb_w_inner(j,i) +1, nzt2955 DO k = nzb_w_inner(j,i)+1, nzt 2902 2956 u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu 2903 2957 flux_r(k) = u_comp * ( & 2904 7. * ( w(k,j,i+1) + w(k,j,i) ) & 2905 - ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3 2906 diss_r(k) = - abs(u_comp) * ( & 2907 3. * ( w(k,j,i+1) - w(k,j,i) ) & 2908 - ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3 2958 7.0 * ( w(k,j,i+1) + w(k,j,i) ) & 2959 - ( w(k,j,i+2) + w(k,j,i-1) ) ) & 2960 * adv_mom_3 2961 diss_r(k) = - ABS( u_comp ) * ( & 2962 3.0 * ( w(k,j,i+1) - w(k,j,i) ) & 2963 - ( w(k,j,i+2) - w(k,j,i-1) ) ) & 2964 * adv_mom_3 2909 2965 ENDDO 2910 2966 2911 2967 CASE ( 6 ) 2912 DO k = nzb_w_inner(j,i) +1, nzt2968 DO k = nzb_w_inner(j,i)+1, nzt 2913 2969 u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu 2914 2970 flux_r(k) = u_comp *( & 2915 7. * ( w(k,j,i+1) + w(k,j,i) ) & 2916 - ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3 2917 diss_r(k) = - abs(u_comp) * ( & 2918 3. * ( w(k,j,i+1) - w(k,j,i) ) & 2919 - ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3 2971 7.0 * ( w(k,j,i+1) + w(k,j,i) ) & 2972 - ( w(k,j,i+2) + w(k,j,i-1) ) ) & 2973 * adv_mom_3 2974 diss_r(k) = - ABS( u_comp ) * ( & 2975 3.0 * ( w(k,j,i+1) - w(k,j,i) ) & 2976 - ( w(k,j,i+2) - w(k,j,i-1) ) ) & 2977 * adv_mom_3 2920 2978 ENDDO 2921 2979 2922 2980 CASE ( 7 ) 2923 DO k = nzb_w_inner(j,i) +1, nzt2981 DO k = nzb_w_inner(j,i)+1, nzt 2924 2982 v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv 2925 2983 flux_n(k) = v_comp *( & 2926 7. * ( w(k,j+1,i) + w(k,j,i) ) & 2927 - ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3 2928 diss_n(k) = - abs(v_comp) * ( & 2929 3. * ( w(k,j+1,i) - w(k,j,i) ) & 2930 - ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3 2984 7.0 * ( w(k,j+1,i) + w(k,j,i) ) & 2985 - ( w(k,j+2,i) + w(k,j-1,i) ) ) & 2986 * adv_mom_3 2987 diss_n(k) = - ABS( v_comp ) * ( & 2988 3.0 * ( w(k,j+1,i) - w(k,j,i) ) & 2989 - ( w(k,j+2,i) - w(k,j-1,i) ) ) & 2990 * adv_mom_3 2931 2991 ENDDO 2932 2992 2933 2993 CASE ( 8 ) 2934 DO k = nzb_w_inner(j,i) +1, nzt2994 DO k = nzb_w_inner(j,i)+1, nzt 2935 2995 v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv 2936 2996 flux_n(k) = v_comp * ( & 2937 7. * ( w(k,j+1,i) + w(k,j,i) ) & 2938 - ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3 2939 diss_n(k) = - abs(v_comp) * ( & 2940 3. * ( w(k,j+1,i) - w(k,j,i) ) & 2941 - ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3 2997 7.0 * ( w(k,j+1,i) + w(k,j,i) ) & 2998 - ( w(k,j+2,i) + w(k,j-1,i) ) ) & 2999 * adv_mom_3 3000 diss_n(k) = - ABS( v_comp ) * ( & 3001 3.0 * ( w(k,j+1,i) - w(k,j,i) ) & 3002 - ( w(k,j+2,i) - w(k,j-1,i) ) ) & 3003 * adv_mom_3 2942 3004 2943 3005 ENDDO … … 2947 3009 END SELECT 2948 3010 ! 2949 !-- Compute the crosswise 5th order fluxes at the outflow 2950 IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 & 2951 .OR. boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 ) THEN 3011 !-- Compute the crosswise 5th order fluxes at the outflow 3012 IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 .OR.& 3013 boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 ) & 3014 THEN 2952 3015 2953 DO k = nzb_w_inner(j,i) +1, nzt3016 DO k = nzb_w_inner(j,i)+1, nzt 2954 3017 v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv 2955 3018 flux_n(k) = v_comp * ( & 2956 37. * ( w(k,j+1,i) + w(k,j,i) ) & 2957 - 8. * ( w(k,j+2,i) + w(k,j-1,i) ) & 2958 + ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5 2959 diss_n(k) = - abs(v_comp) * ( & 2960 10. * ( w(k,j+1,i) - w(k,j,i) ) & 2961 - 5. * ( w(k,j+2,i) - w(k,j-1,i) ) & 2962 + ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5 3019 37.0 * ( w(k,j+1,i) + w(k,j,i) ) & 3020 - 8.0 * ( w(k,j+2,i) + w(k,j-1,i) ) & 3021 + ( w(k,j+3,i) + w(k,j-2,i) ) ) & 3022 * adv_mom_5 3023 diss_n(k) = - ABS( v_comp ) * ( & 3024 10.0 * ( w(k,j+1,i) - w(k,j,i) ) & 3025 - 5.0 * ( w(k,j+2,i) - w(k,j-1,i) ) & 3026 + ( w(k,j+3,i) - w(k,j-2,i) ) ) & 3027 * adv_mom_5 2963 3028 ENDDO 2964 3029 2965 3030 ELSE 2966 3031 2967 DO k = nzb_w_inner(j,i) +1, nzt3032 DO k = nzb_w_inner(j,i)+1, nzt 2968 3033 u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu 2969 3034 flux_r(k) = u_comp * ( & 2970 37. * ( w(k,j,i+1) + w(k,j,i) ) & 2971 - 8. * ( w(k,j,i+2) + w(k,j,i-1) ) & 2972 + ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5 2973 diss_r(k) = - abs(u_comp) * ( & 2974 10. * ( w(k,j,i+1) - w(k,j,i) ) & 2975 - 5. * ( w(k,j,i+2) - w(k,j,i-1) ) & 2976 + ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5 3035 37.0 * ( w(k,j,i+1) + w(k,j,i) ) & 3036 - 8.0 * ( w(k,j,i+2) + w(k,j,i-1) ) & 3037 + ( w(k,j,i+3) + w(k,j,i-2) ) ) & 3038 * adv_mom_5 3039 diss_r(k) = - ABS( u_comp ) * ( & 3040 10.0 * ( w(k,j,i+1) - w(k,j,i) ) & 3041 - 5.0 * ( w(k,j,i+2) - w(k,j,i-1) ) & 3042 + ( w(k,j,i+3) - w(k,j,i-2) ) ) & 3043 * adv_mom_5 2977 3044 ENDDO 2978 3045 … … 2982 3049 ! 2983 3050 !-- Compute the fifth order fluxes for the interior of PE domain. 2984 DO k = nzb_w_inner(j,i) +1, nzt3051 DO k = nzb_w_inner(j,i)+1, nzt 2985 3052 u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu 2986 3053 flux_r(k) = u_comp * ( & 2987 37. * ( w(k,j,i+1) + w(k,j,i) ) & 2988 - 8. * ( w(k,j,i+2) + w(k,j,i-1) ) & 2989 + ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5 2990 diss_r(k) = - abs(u_comp) * ( & 2991 10. * ( w(k,j,i+1) - w(k,j,i) ) & 2992 - 5. * ( w(k,j,i+2) - w(k,j,i-1) ) & 2993 + ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5 3054 37.0 * ( w(k,j,i+1) + w(k,j,i) ) & 3055 - 8.0 * ( w(k,j,i+2) + w(k,j,i-1) ) & 3056 + ( w(k,j,i+3) + w(k,j,i-2) ) ) & 3057 * adv_mom_5 3058 diss_r(k) = - ABS( u_comp ) * ( & 3059 10.0 * ( w(k,j,i+1) - w(k,j,i) ) & 3060 - 5.0 * ( w(k,j,i+2) - w(k,j,i-1) ) & 3061 + ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5 2994 3062 2995 3063 v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv 2996 3064 flux_n(k) = v_comp * ( & 2997 37. * ( w(k,j+1,i) + w(k,j,i) )&2998 - 8. * ( w(k,j+2,i) + w(k,j-1,i) )&2999 +( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_53000 diss_n(k) = - abs(v_comp) * (&3001 10. * ( w(k,j+1,i) - w(k,j,i) )&3002 - 5. * ( w(k,j+2,i) - w(k,j-1,i) )&3003 +( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_53065 37.0 * ( w(k,j+1,i) + w(k,j,i) ) & 3066 - 8.0 * ( w(k,j+2,i) + w(k,j-1,i) ) & 3067 + ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5 3068 diss_n(k) = - ABS( v_comp ) * ( & 3069 10.0 * ( w(k,j+1,i) - w(k,j,i) ) & 3070 - 5.0 * ( w(k,j+2,i) - w(k,j-1,i) ) & 3071 + ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5 3004 3072 ENDDO 3005 3073 … … 3009 3077 tend(k,j,i) = tend(k,j,i) - ( & 3010 3078 ( flux_r(k) + diss_r(k) & 3011 - ( swap_flux_x_local_w(k,j) + swap_diss_x_local_w(k,j) ) ) * ddx&3079 - swap_flux_x_local_w(k,j) - swap_diss_x_local_w(k,j) ) * ddx & 3012 3080 + ( flux_n(k) + diss_n(k) & 3013 - ( swap_flux_y_local_w(k) + swap_diss_y_local_w(k) )) * ddy )3081 - swap_flux_y_local_w(k) - swap_diss_y_local_w(k) ) * ddy ) 3014 3082 3015 3083 swap_flux_x_local_w(k,j) = flux_r(k) … … 3027 3095 DO i = nxl, nxr 3028 3096 DO j = nys, nyn 3029 k = nzb_w_inner(j,i) +13097 k = nzb_w_inner(j,i)+1 3030 3098 flux_d = w(k-1,j,i) * ( w(k,j,i) + w(k-1,j,i)) 3031 3099 flux_t(k-1) = flux_d … … 3033 3101 diss_t(k-1) = diss_d 3034 3102 ! 3035 !-- 2nd order scheme 3103 !-- 2nd order scheme (bottom) 3036 3104 w_comp = w(k+1,j,i) + w(k,j,i) 3037 3105 flux_t(k) = w_comp * ( w(k+1,j,i) + w(k,j,i) ) * 0.25 3038 diss_t(k) = diss_2nd( w(k+2,j,i), w(k+1,j,i), w(k,j,i), 0. , 0.,&3106 diss_t(k) = diss_2nd( w(k+2,j,i), w(k+1,j,i), w(k,j,i), 0.0, 0.0, & 3039 3107 w_comp, 0.25, ddzu(k+1) ) 3040 3108 3041 3109 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 3042 - ( flux_d + diss_d )) * ddzu(k+1)3110 - flux_d - diss_d ) * ddzu(k+1) 3043 3111 ! 3044 !-- WS3 as an intermediate step 3112 !-- WS3 as an intermediate step (bottom) 3045 3113 k = nzb_w_inner(j,i) + 2 3114 flux_d = flux_t(k-1) 3115 diss_d = diss_t(k-1) 3116 w_comp = w(k+1,j,i) + w(k,j,i) 3117 flux_t(k) = w_comp * ( & 3118 7.0 * ( w(k+1,j,i) + w(k,j,i) ) & 3119 - ( w(k+2,j,i) + w(k-1,j,i) ) ) * adv_mom_3 3120 diss_t(k) = - ABS( w_comp ) * ( & 3121 3.0 * ( w(k+1,j,i) - w(k,j,i) ) & 3122 - ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3 3123 3124 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 3125 - flux_d - diss_d ) * ddzu(k+1) 3126 ! 3127 !-- WS5 3128 DO k = nzb_v_inner(j,i)+3, nzt-2 3129 flux_d = flux_t(k-1) 3130 diss_d = diss_t(k-1) 3131 3132 w_comp = w(k+1,j,i) + w(k,j,i) 3133 flux_t(k) = w_comp * ( & 3134 37.0 * ( w(k+1,j,i) + w(k,j,i) ) & 3135 - 8.0 * ( w(k+2,j,i) + w(k-1,j,i) ) & 3136 + ( w(k+3,j,i) + w(k-2,j,i) ) ) * adv_mom_5 3137 diss_t(k) = - ABS( w_comp ) * ( & 3138 10.0 * ( w(k+1,j,i) - w(k,j,i) ) & 3139 - 5.0 * ( w(k+2,j,i) - w(k-1,j,i) ) & 3140 + ( w(k+3,j,i) - w(k-2,j,i) ) ) * adv_mom_5 3141 3142 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 3143 - flux_d - diss_d ) * ddzu(k+1) 3144 ENDDO 3145 ! 3146 !-- WS3 as an intermediate step (top) 3147 k = nzt-1 3046 3148 flux_d = flux_t(k-1) 3047 3149 diss_d = diss_t(k-1) 3048 3150 w_comp = w(k+1,j,i)+w(k,j,i) 3049 3151 flux_t(k) = w_comp * ( & 3050 7. * ( w(k+1,j,i) + w(k,j,i) )&3051 -( w(k+2,j,i) + w(k-1,j,i) ) ) * adv_mom_33052 diss_t(k) = - abs(w_comp) * (&3053 3. * ( w(k+1,j,i) - w(k,j,i) )&3054 -( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_33152 7.0 * ( w(k+1,j,i) + w(k,j,i) ) & 3153 - ( w(k+2,j,i) + w(k-1,j,i) ) ) * adv_mom_3 3154 diss_t(k) = - ABS( w_comp ) * ( & 3155 3.0 * ( w(k+1,j,i) - w(k,j,i) ) & 3156 - ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3 3055 3157 3056 3158 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 3057 - ( flux_d + diss_d ) ) * ddzu(k+1) 3058 ! 3059 !-- WS5 3060 DO k = nzb_v_inner(j,i) + 3, nzt 3061 flux_d = flux_t(k-1) 3062 diss_d = diss_t(k-1) 3063 3064 w_comp = w(k+1,j,i) + w(k,j,i) 3065 flux_t(k) = w_comp * ( & 3066 37. * ( w(k+1,j,i) + w(k,j,i) ) & 3067 - 8. * ( w(k+2,j,i) + w(k-1,j,i) ) & 3068 + ( w(k+3,j,i) + w(k-2,j,i) ) ) * adv_mom_5 3069 diss_t(k) = - abs(w_comp) * ( & 3070 10. * ( w(k+1,j,i) - w(k,j,i) ) & 3071 - 5. * ( w(k+2,j,i) - w(k-1,j,i) ) & 3072 + ( w(k+3,j,i) - w(k-2,j,i) ) ) * adv_mom_5 3073 3074 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 3075 - ( flux_d + diss_d ) ) * ddzu(k+1) 3076 ENDDO 3077 ! 3078 !-- WS3 as an intermediate step 3079 DO k = nzt - 2, nzt - 1 3080 flux_d = flux_t(k-1) 3081 diss_d = diss_t(k-1) 3082 w_comp = w(k+1,j,i)+w(k,j,i) 3083 flux_t(k) = w_comp * ( & 3084 7. * ( w(k+1,j,i) + w(k,j,i) ) & 3085 - ( w(k+2,j,i) + w(k-1,j,i) ) ) * adv_mom_3 3086 diss_t(k) = - abs(w_comp) * ( & 3087 3. * ( w(k+1,j,i) - w(k,j,i) ) & 3088 - ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3 3089 3090 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 3091 - ( flux_d + diss_d ) ) * ddzu(k+1) 3092 ENDDO 3093 ! 3094 !-- 2nd order scheme 3159 - flux_d - diss_d ) * ddzu(k+1) 3160 ! 3161 !-- 2nd order scheme (top) 3095 3162 k = nzt 3096 3163 flux_d = flux_t(k-1) … … 3103 3170 3104 3171 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & 3105 - ( flux_d + diss_d )) * ddzu(k+1)3172 - flux_d - diss_d ) * ddzu(k+1) 3106 3173 ! 3107 3174 !-- at last vertical momentum flux is accumulated 3108 3175 DO k = nzb_w_inner(j,i), nzt 3109 3176 sums_ws2_ws_l(k,:) = sums_ws2_ws_l(k,:) & 3110 + ( flux_t(k) + diss_t(k) )&3111 * weight_substep(intermediate_timestep_count)&3112 *rmask(j,i,:)3177 + ( flux_t(k) + diss_t(k) ) & 3178 * weight_substep(intermediate_timestep_count) & 3179 * rmask(j,i,:) 3113 3180 ENDDO 3114 3181 … … 3131 3198 ! 3132 3199 !-- Computation of 2nd order dissipation. This numerical dissipation is 3133 !-- necessary to keep a stable numerical solution in region where the3134 !-- order of the schemes degraded.3200 !-- necessary to keep a stable numerical solution in regions where the 3201 !-- order of the schemes is degraded. 3135 3202 3136 3203 REAL FUNCTION diss_2nd( v2, v1, v0, vm1, vm2, vel_comp, factor, grid ) & 3137 RESULT( dissip)3204 RESULT( dissip ) 3138 3205 3139 3206 IMPLICIT NONE … … 3144 3211 value_min_p, value_max_p, diss_m, diss_0, diss_p 3145 3212 3146 value_min_m = MIN( v0, vm1, vm2)3147 value_max_m = MAX( v0, vm1, vm2)3148 value_min = MIN(v1, v0, vm2)3149 value_max = MAX(v1, v0, vm2)3150 value_min_p = MIN( v2, v1, v0)3151 value_max_p = MAX( v2, v1, v0)3213 value_min_m = MIN( v0, vm1, vm2 ) 3214 value_max_m = MAX( v0, vm1, vm2 ) 3215 value_min = MIN( v1, v0, vm2 ) 3216 value_max = MAX( v1, v0, vm2 ) 3217 value_min_p = MIN( v2, v1, v0 ) 3218 value_max_p = MAX( v2, v1, v0 ) 3152 3219 3153 diss_m = MAX( 0.,vm1 - value_max_m) + MIN(0.,vm1 - value_min_m)3154 diss_0 = MAX( 0.,v0 - value_max) + MIN(0.,v0 - value_min)3155 diss_p = MAX( 0.,v1 - value_max_p) + MIN(0.,v1 - value_min_p)3220 diss_m = MAX( 0.0, vm1 - value_max_m ) + MIN( 0.0, vm1 - value_min_m ) 3221 diss_0 = MAX( 0.0, v0 - value_max ) + MIN( 0.0, v0 - value_min ) 3222 diss_p = MAX( 0.0, v1 - value_max_p ) + MIN( 0.0, v1 - value_min_p ) 3156 3223 3157 dissip = abs(vel_comp) * (diss_p - 2.*diss_0 + diss_m)&3158 * grid**2 * factor3224 dissip = ABS( vel_comp ) * ( diss_p - 2.0 * diss_0 + diss_m ) & 3225 * grid**2 * factor 3159 3226 3160 3227 END FUNCTION diss_2nd -
palm/trunk/SOURCE/init_3d_model.f90
r710 r713 7 7 ! Current revisions: 8 8 ! ----------------- 9 ! 9 ! ! Reformulate weight_substep and weight_pres as broken numbers. 10 10 ! 11 11 ! Former revisions: … … 1353 1353 IF ( TRIM(timestep_scheme) == 'runge-kutta-3' ) THEN ! for RK3-method 1354 1354 1355 weight_substep(1) = 0.1666666666666661356 weight_substep(2) = 0.31357 weight_substep(3) = 0.5333333333333331358 1359 weight_pres(1) = 0.3333333333333331360 weight_pres(2) = 0.4166666666666661361 weight_pres(3) = 0.251355 weight_substep(1) = 1./6. 1356 weight_substep(2) = 3./10. 1357 weight_substep(3) = 8./15. 1358 1359 weight_pres(1) = 1./3. 1360 weight_pres(2) = 5./12. 1361 weight_pres(3) = 1./4. 1362 1362 1363 1363 ELSEIF ( TRIM(timestep_scheme) == 'runge-kutta-2' ) THEN ! for RK2-method 1364 1364 1365 weight_substep(1) = 0.51366 weight_substep(2) = 0.51365 weight_substep(1) = 1./2. 1366 weight_substep(2) = 1./2. 1367 1367 1368 weight_pres(1) = 0.51369 weight_pres(2) = 0.51368 weight_pres(1) = 1./2. 1369 weight_pres(2) = 1./2. 1370 1370 1371 1371 ELSE ! for Euler- and leapfrog-method
Note: See TracChangeset
for help on using the changeset viewer.
