Changeset 87 for palm/trunk/SOURCE/flow_statistics.f90
- Timestamp:
- May 22, 2007 3:46:47 PM (17 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/flow_statistics.f90
r83 r87 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! Two more arguments added to user_statistics, which is now also called for 7 ! user-defined profiles, 8 ! var_hom and var_sum renamed pr_palm 7 9 ! 8 10 ! Former revisions: … … 87 89 !-- WARNING: next line still has to be adjusted for OpenMP 88 90 sums_l(:,21,0) = sums_wsts_bc_l(:,sr) ! heat flux from advec_s_bc 89 sums_l(nzb+9, var_sum,0) = sums_divold_l(sr) ! old divergence from pres90 sums_l(nzb+10, var_sum,0) = sums_divnew_l(sr) ! new divergence from pres91 sums_l(nzb+9,pr_palm,0) = sums_divold_l(sr) ! old divergence from pres 92 sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr) ! new divergence from pres 91 93 !-- WARNING: next four lines still may have to be adjusted for OpenMP 92 sums_l(nzb:nzb+2, var_sum-1,0) = sums_up_fraction_l(1,1:3,sr)! upstream93 sums_l(nzb+3:nzb+5, var_sum-1,0) = sums_up_fraction_l(2,1:3,sr)! parts94 sums_l(nzb+6:nzb+8, var_sum-1,0) = sums_up_fraction_l(3,1:3,sr)! from95 sums_l(nzb+9:nzb+11, var_sum-1,0) = sums_up_fraction_l(4,1:3,sr)! spline94 sums_l(nzb:nzb+2,pr_palm-1,0) = sums_up_fraction_l(1,1:3,sr)! upstream 95 sums_l(nzb+3:nzb+5,pr_palm-1,0) = sums_up_fraction_l(2,1:3,sr)! parts 96 sums_l(nzb+6:nzb+8,pr_palm-1,0) = sums_up_fraction_l(3,1:3,sr)! from 97 sums_l(nzb+9:nzb+11,pr_palm-1,0) = sums_up_fraction_l(4,1:3,sr)! spline 96 98 97 99 ! … … 316 318 !-- quantities is seperated from the previous loop in order to 317 319 !-- allow vectorization of that loop. 318 sums_l(nzb+4, var_sum,tn) = sums_l(nzb+4,var_sum,tn) + sums_l_etot319 sums_l(nzb+5, var_sum,tn) = sums_l(nzb+5,var_sum,tn) + sums_l_eper320 sums_l(nzb+4,pr_palm,tn) = sums_l(nzb+4,pr_palm,tn) + sums_l_etot 321 sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn) + sums_l_eper 320 322 ! 321 323 !-- 2D-arrays (being collected in the last column of sums_l) 322 sums_l(nzb, var_sum,tn) = sums_l(nzb,var_sum,tn) + &324 sums_l(nzb,pr_palm,tn) = sums_l(nzb,pr_palm,tn) + & 323 325 us(j,i) * rmask(j,i,sr) 324 sums_l(nzb+1, var_sum,tn) = sums_l(nzb+1,var_sum,tn) + &326 sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) + & 325 327 usws(j,i) * rmask(j,i,sr) 326 sums_l(nzb+2, var_sum,tn) = sums_l(nzb+2,var_sum,tn) + &328 sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) + & 327 329 vsws(j,i) * rmask(j,i,sr) 328 sums_l(nzb+3, var_sum,tn) = sums_l(nzb+3,var_sum,tn) + &330 sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) + & 329 331 ts(j,i) * rmask(j,i,sr) 330 332 ENDDO … … 652 654 653 655 ENDIF 656 657 ! 658 !-- Calculate the user-defined profiles 659 CALL user_statistics( 'profiles', sr, tn ) 654 660 !$OMP END PARALLEL 655 661 … … 660 666 sums_l(:,3,0) = sums_l(:,3,0) + sums_l(:,3,i) 661 667 sums_l(:,4:40,0) = sums_l(:,4:40,0) + sums_l(:,4:40,i) 662 sums_l(:,45:var_sum,0) = sums_l(:,45:var_sum,0) + & 663 sums_l(:,45:var_sum,i) 668 sums_l(:,45:pr_palm,0) = sums_l(:,45:pr_palm,0) + & 669 sums_l(:,45:pr_palm,i) 670 IF ( max_pr_user > 0 ) THEN 671 sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) = & 672 sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + & 673 sums_l(:,pr_palm+1:pr_palm+max_pr_user,i) 674 ENDIF 664 675 ENDDO 665 676 ENDIF … … 679 690 !-- Profiles: 680 691 DO k = nzb, nzt+1 681 sums(k,: var_sum-2) = sums(k,:var_sum-2) / ngp_2dh_outer(k,sr)692 sums(k,:pr_palm-2) = sums(k,:pr_palm-2) / ngp_2dh_outer(k,sr) 682 693 ENDDO 683 694 !-- Upstream-parts 684 sums(nzb:nzb+11, var_sum-1) = sums(nzb:nzb+11,var_sum-1) / ngp_3d(sr)695 sums(nzb:nzb+11,pr_palm-1) = sums(nzb:nzb+11,pr_palm-1) / ngp_3d(sr) 685 696 !-- u* and so on 686 !-- As sums(nzb:nzb+3, var_sum) are full 2D arrays (us, usws, vsws, ts) whose697 !-- As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose 687 698 !-- size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer 688 699 !-- above the topography, they are being divided by ngp_2dh(sr) 689 sums(nzb:nzb+3, var_sum) = sums(nzb:nzb+3,var_sum) / &700 sums(nzb:nzb+3,pr_palm) = sums(nzb:nzb+3,pr_palm) / & 690 701 ngp_2dh(sr) 691 702 !-- eges, e* 692 sums(nzb+4:nzb+5, var_sum) = sums(nzb+4:nzb+5,var_sum) / &703 sums(nzb+4:nzb+5,pr_palm) = sums(nzb+4:nzb+5,pr_palm) / & 693 704 ngp_3d_inner(sr) 694 705 !-- Old and new divergence 695 sums(nzb+9:nzb+10, var_sum) = sums(nzb+9:nzb+10,var_sum) / &706 sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / & 696 707 ngp_3d_inner(sr) 708 709 !-- User-defined profiles 710 IF ( max_pr_user > 0 ) THEN 711 DO k = nzb, nzt+1 712 sums(k,pr_palm+1:pr_palm+max_pr_user) = & 713 sums(k,pr_palm+1:pr_palm+max_pr_user) / & 714 ngp_2dh_outer(k,sr) 715 ENDDO 716 ENDIF 697 717 698 718 ! … … 748 768 hom(:,1,63,sr) = sums(:,61) + sums(:,62) ! vpt_t 749 769 750 hom(:,1, var_hom-1,sr) = sums(:,var_sum-1)770 hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1) 751 771 ! upstream-parts u_x, u_y, u_z, v_x, 752 772 ! v_y, usw. (in last but one profile) 753 hom(:,1, var_hom,sr) = sums(:,var_sum)773 hom(:,1,pr_palm,sr) = sums(:,pr_palm) 754 774 ! u*, w'u', w'v', t* (in last profile) 775 776 IF ( max_pr_user > 0 ) THEN ! user-defined profiles 777 hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = & 778 sums(:,pr_palm+1:pr_palm+max_pr_user) 779 ENDIF 755 780 756 781 ! … … 792 817 ENDDO 793 818 794 hom(nzb+6,1, var_hom,sr) = z_i(1)795 hom(nzb+7,1, var_hom,sr) = z_i(2)819 hom(nzb+6,1,pr_palm,sr) = z_i(1) 820 hom(nzb+7,1,pr_palm,sr) = z_i(2) 796 821 797 822 ! … … 800 825 !-- The horizontal average at nzb+1 is input for the average temperature. 801 826 IF ( hom(nzb,1,18,sr) > 0.0 .AND. z_i(1) /= 0.0 ) THEN 802 hom(nzb+8,1, var_hom,sr) = ( g / hom(nzb+1,1,4,sr) * &827 hom(nzb+8,1,pr_palm,sr) = ( g / hom(nzb+1,1,4,sr) * & 803 828 hom(nzb,1,18,sr) * z_i(1) )**0.333333333 804 829 !-- so far this only works if Prandtl layer is used 805 hom(nzb+11,1, var_hom,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,var_hom,sr)830 hom(nzb+11,1,pr_palm,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,pr_palm,sr) 806 831 ELSE 807 hom(nzb+8,1, var_hom,sr) = 0.0808 hom(nzb+11,1, var_hom,sr) = 0.0832 hom(nzb+8,1,pr_palm,sr) = 0.0 833 hom(nzb+11,1,pr_palm,sr) = 0.0 809 834 ENDIF 810 835 811 836 ! 812 837 !-- Collect the time series quantities 813 ts_value(1,sr) = hom(nzb+4,1, var_hom,sr) ! E814 ts_value(2,sr) = hom(nzb+5,1, var_hom,sr) ! E*838 ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr) ! E 839 ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr) ! E* 815 840 ts_value(3,sr) = dt_3d 816 ts_value(4,sr) = hom(nzb,1, var_hom,sr) ! u*817 ts_value(5,sr) = hom(nzb+3,1, var_hom,sr) ! th*841 ts_value(4,sr) = hom(nzb,1,pr_palm,sr) ! u* 842 ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr) ! th* 818 843 ts_value(6,sr) = u_max 819 844 ts_value(7,sr) = v_max 820 845 ts_value(8,sr) = w_max 821 ts_value(9,sr) = hom(nzb+10,1, var_sum,sr) ! new divergence822 ts_value(10,sr) = hom(nzb+9,1, var_hom,sr) ! old Divergence823 ts_value(11,sr) = hom(nzb+6,1, var_hom,sr) ! z_i(1)824 ts_value(12,sr) = hom(nzb+7,1, var_hom,sr) ! z_i(2)825 ts_value(13,sr) = hom(nzb+8,1, var_hom,sr) ! w*846 ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr) ! new divergence 847 ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr) ! old Divergence 848 ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr) ! z_i(1) 849 ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr) ! z_i(2) 850 ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr) ! w* 826 851 ts_value(14,sr) = hom(nzb,1,16,sr) ! w'pt' at k=0 827 852 ts_value(15,sr) = hom(nzb+1,1,16,sr) ! w'pt' at k=1 … … 829 854 ts_value(17,sr) = hom(nzb,1,4,sr) ! pt(0) 830 855 ts_value(18,sr) = hom(nzb+1,1,4,sr) ! pt(zp) 831 ts_value(19,sr) = hom(nzb+9,1, var_hom-1,sr) ! splptx832 ts_value(20,sr) = hom(nzb+10,1, var_hom-1,sr) ! splpty833 ts_value(21,sr) = hom(nzb+11,1, var_hom-1,sr) ! splptz856 ts_value(19,sr) = hom(nzb+9,1,pr_palm-1,sr) ! splptx 857 ts_value(20,sr) = hom(nzb+10,1,pr_palm-1,sr) ! splpty 858 ts_value(21,sr) = hom(nzb+11,1,pr_palm-1,sr) ! splptz 834 859 IF ( ts_value(5,sr) /= 0.0 ) THEN 835 860 ts_value(22,sr) = ts_value(4,sr)**2 / & … … 841 866 ! 842 867 !-- Calculate additional statistics provided by the user interface 843 CALL user_statistics( sr)868 CALL user_statistics( 'time_series', sr, 0 ) 844 869 845 870 ENDDO ! loop of the subregions
Note: See TracChangeset
for help on using the changeset viewer.