Changeset 2232 for palm/trunk/SOURCE/flow_statistics.f90
- Timestamp:
- May 30, 2017 5:47:52 PM (7 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/flow_statistics.f90 (modified) (48 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/flow_statistics.f90
r2119 r2232 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Adjustments to new topography and surface concept 23 23 ! 24 24 ! Former revisions: … … 224 224 USE arrays_3d, & 225 225 ONLY: ddzu, ddzw, e, heatflux_output_conversion, hyp, km, kh, & 226 momentumflux_output_conversion, nr, ol, p, prho, prr, pt, q, & 227 qc, ql, qr, qs, qsws, qswst, rho_air, rho_air_zw, rho_ocean, s, & 228 sa, ss, ssws, sswst, saswsb, saswst, shf, td_lsa_lpt, td_lsa_q, & 229 td_sub_lpt, td_sub_q, time_vert, ts, tswst, u, ug, us, usws, & 230 uswst, vsws, v, vg, vpt, vswst, w, w_subs, & 231 waterflux_output_conversion, zw 226 momentumflux_output_conversion, nr, p, prho, prr, pt, q, & 227 qc, ql, qr, rho_air, rho_air_zw, rho_ocean, s, & 228 sa, td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q, time_vert, u, & 229 ug, v, vg, vpt, w, w_subs, waterflux_output_conversion, zw 232 230 233 231 USE cloud_parameters, & … … 236 234 USE control_parameters, & 237 235 ONLY: average_count_pr, cloud_droplets, cloud_physics, do_sum, & 238 dt_3d, g, humidity, kappa, la rge_scale_forcing,&236 dt_3d, g, humidity, kappa, land_surface, large_scale_forcing, & 239 237 large_scale_subsidence, max_pr_user, message_string, neutral, & 240 238 microphysics_seifert, ocean, passive_scalar, simulated_time, & … … 250 248 USE indices, & 251 249 ONLY: ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, ngp_sums, & 252 ngp_sums_ls, nxl, nxr, nyn, nys, nzb, nzb_diff_s_inner, & 253 nzb_s_inner, nzt, nzt_diff 250 ngp_sums_ls, nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0 254 251 255 252 USE kinds 256 253 257 254 USE land_surface_model_mod, & 258 ONLY: ghf_eb, land_surface, m_soil, nzb_soil, nzt_soil, & 259 qsws_eb, qsws_liq_eb, qsws_soil_eb, qsws_veg_eb, r_a, r_s, & 260 shf_eb, t_soil 255 ONLY: m_soil_h, nzb_soil, nzt_soil, t_soil_h 261 256 262 257 USE netcdf_interface, & … … 277 272 USE statistics 278 273 274 USE surface_mod, & 275 ONLY : surf_def_h, surf_lsm_h, surf_usm_h 276 279 277 280 278 IMPLICIT NONE … … 283 281 INTEGER(iwp) :: j !< 284 282 INTEGER(iwp) :: k !< 283 INTEGER(iwp) :: ki !< 285 284 INTEGER(iwp) :: k_surface_level !< 285 INTEGER(iwp) :: m !< loop variable over all horizontal wall elements 286 INTEGER(iwp) :: l !< loop variable over surface facing -- up- or downward-facing 286 287 INTEGER(iwp) :: nt !< 287 288 INTEGER(iwp) :: omp_get_thread_num !< 288 289 INTEGER(iwp) :: sr !< 289 290 INTEGER(iwp) :: tn !< 290 291 291 292 LOGICAL :: first !< 292 293 293 294 REAL(wp) :: dptdz_threshold !< 294 295 REAL(wp) :: fac !< 296 REAL(wp) :: flag !< 295 297 REAL(wp) :: height !< 296 298 REAL(wp) :: pts !< … … 400 402 !-- for other horizontal averages. 401 403 tn = 0 402 403 !$OMP PARALLEL PRIVATE( i, j, k, tn ) 404 !$ tn = omp_get_thread_num() 405 404 !$OMP PARALLEL PRIVATE( i, j, k, tn, flag ) 405 !$ tn = omp_get_thread_num() 406 406 !$OMP DO 407 407 DO i = nxl, nxr 408 408 DO j = nys, nyn 409 DO k = nzb_s_inner(j,i), nzt+1 410 sums_l(k,1,tn) = sums_l(k,1,tn) + u(k,j,i) * rmask(j,i,sr) 411 sums_l(k,2,tn) = sums_l(k,2,tn) + v(k,j,i) * rmask(j,i,sr) 412 sums_l(k,4,tn) = sums_l(k,4,tn) + pt(k,j,i) * rmask(j,i,sr) 409 DO k = nzb, nzt+1 410 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) ) 411 sums_l(k,1,tn) = sums_l(k,1,tn) + u(k,j,i) * rmask(j,i,sr) & 412 * flag 413 sums_l(k,2,tn) = sums_l(k,2,tn) + v(k,j,i) * rmask(j,i,sr) & 414 * flag 415 sums_l(k,4,tn) = sums_l(k,4,tn) + pt(k,j,i) * rmask(j,i,sr) & 416 * flag 413 417 ENDDO 414 418 ENDDO … … 421 425 DO i = nxl, nxr 422 426 DO j = nys, nyn 423 DO k = nzb_s_inner(j,i), nzt+1 424 sums_l(k,23,tn) = sums_l(k,23,tn) + & 425 sa(k,j,i) * rmask(j,i,sr) 427 DO k = nzb, nzt+1 428 sums_l(k,23,tn) = sums_l(k,23,tn) + sa(k,j,i) & 429 * rmask(j,i,sr) & 430 * MERGE( 1.0_wp, 0.0_wp, & 431 BTEST( wall_flags_0(k,j,i), 22 ) ) 426 432 ENDDO 427 433 ENDDO … … 437 443 DO i = nxl, nxr 438 444 DO j = nys, nyn 439 DO k = nzb_s_inner(j,i), nzt+1 440 sums_l(k,44,tn) = sums_l(k,44,tn) + & 441 vpt(k,j,i) * rmask(j,i,sr) 442 sums_l(k,41,tn) = sums_l(k,41,tn) + & 443 q(k,j,i) * rmask(j,i,sr) 445 DO k = nzb, nzt+1 446 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) ) 447 sums_l(k,44,tn) = sums_l(k,44,tn) + & 448 vpt(k,j,i) * rmask(j,i,sr) * flag 449 sums_l(k,41,tn) = sums_l(k,41,tn) + & 450 q(k,j,i) * rmask(j,i,sr) * flag 444 451 ENDDO 445 452 ENDDO … … 449 456 DO i = nxl, nxr 450 457 DO j = nys, nyn 451 DO k = nzb_s_inner(j,i), nzt+1 452 sums_l(k,42,tn) = sums_l(k,42,tn) + & 453 ( q(k,j,i) - ql(k,j,i) ) * rmask(j,i,sr) 454 sums_l(k,43,tn) = sums_l(k,43,tn) + ( & 458 DO k = nzb, nzt+1 459 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) ) 460 sums_l(k,42,tn) = sums_l(k,42,tn) + & 461 ( q(k,j,i) - ql(k,j,i) ) * rmask(j,i,sr) & 462 * flag 463 sums_l(k,43,tn) = sums_l(k,43,tn) + ( & 455 464 pt(k,j,i) + l_d_cp*pt_d_t(k) * ql(k,j,i) & 456 ) * rmask(j,i,sr) 465 ) * rmask(j,i,sr) & 466 * flag 457 467 ENDDO 458 468 ENDDO … … 467 477 DO i = nxl, nxr 468 478 DO j = nys, nyn 469 DO k = nzb_s_inner(j,i), nzt+1 470 sums_l(k,117,tn) = sums_l(k,117,tn) + s(k,j,i) * rmask(j,i,sr) 479 DO k = nzb, nzt+1 480 sums_l(k,117,tn) = sums_l(k,117,tn) + s(k,j,i) & 481 * rmask(j,i,sr) & 482 * MERGE( 1.0_wp, 0.0_wp, & 483 BTEST( wall_flags_0(k,j,i), 22 ) ) 471 484 ENDDO 472 485 ENDDO … … 603 616 !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, & 604 617 !$OMP sums_l_etot, tn, ust, ust2, u2, vst, vst2, v2, & 605 !$OMP w2 ) 606 !$ tn = omp_get_thread_num() 607 618 !$OMP w2, flag, m, ki, l ) 619 !$ tn = omp_get_thread_num() 608 620 !$OMP DO 609 621 DO i = nxl, nxr 610 622 DO j = nys, nyn 611 623 sums_l_etot = 0.0_wp 612 DO k = nzb_s_inner(j,i), nzt+1 624 DO k = nzb, nzt+1 625 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) ) 613 626 ! 614 627 !-- Prognostic and diagnostic variables 615 sums_l(k,3,tn) = sums_l(k,3,tn) + w(k,j,i) * rmask(j,i,sr) 616 sums_l(k,8,tn) = sums_l(k,8,tn) + e(k,j,i) * rmask(j,i,sr) 617 sums_l(k,9,tn) = sums_l(k,9,tn) + km(k,j,i) * rmask(j,i,sr) 618 sums_l(k,10,tn) = sums_l(k,10,tn) + kh(k,j,i) * rmask(j,i,sr) 619 sums_l(k,40,tn) = sums_l(k,40,tn) + p(k,j,i) 628 sums_l(k,3,tn) = sums_l(k,3,tn) + w(k,j,i) * rmask(j,i,sr) & 629 * flag 630 sums_l(k,8,tn) = sums_l(k,8,tn) + e(k,j,i) * rmask(j,i,sr) & 631 * flag 632 sums_l(k,9,tn) = sums_l(k,9,tn) + km(k,j,i) * rmask(j,i,sr) & 633 * flag 634 sums_l(k,10,tn) = sums_l(k,10,tn) + kh(k,j,i) * rmask(j,i,sr) & 635 * flag 636 sums_l(k,40,tn) = sums_l(k,40,tn) + p(k,j,i) * flag 620 637 621 638 sums_l(k,33,tn) = sums_l(k,33,tn) + & 622 ( pt(k,j,i)-hom(k,1,4,sr) )**2 * rmask(j,i,sr) 639 ( pt(k,j,i)-hom(k,1,4,sr) )**2 * rmask(j,i,sr)& 640 * flag 623 641 624 642 IF ( humidity ) THEN 625 643 sums_l(k,70,tn) = sums_l(k,70,tn) + & 626 ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr) 644 ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr)& 645 * flag 627 646 ENDIF 628 647 IF ( passive_scalar ) THEN 629 648 sums_l(k,118,tn) = sums_l(k,118,tn) + & 630 ( s(k,j,i)-hom(k,1,117,sr) )**2 * rmask(j,i,sr) 649 ( s(k,j,i)-hom(k,1,117,sr) )**2 * rmask(j,i,sr)& 650 * flag 631 651 ENDIF 632 652 ! 633 653 !-- Higher moments 634 654 !-- (Computation of the skewness of w further below) 635 sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i)**3 * rmask(j,i,sr) 655 sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i)**3 * rmask(j,i,sr) & 656 * flag 636 657 637 658 sums_l_etot = sums_l_etot + & 638 0.5_wp * ( u(k,j,i)**2 + v(k,j,i)**2 + & 639 w(k,j,i)**2 ) * rmask(j,i,sr) 659 0.5_wp * ( u(k,j,i)**2 + v(k,j,i)**2 + & 660 w(k,j,i)**2 ) * rmask(j,i,sr)& 661 * flag 640 662 ENDDO 641 663 ! … … 647 669 ! 648 670 !-- 2D-arrays (being collected in the last column of sums_l) 649 sums_l(nzb,pr_palm,tn) = sums_l(nzb,pr_palm,tn) + & 650 us(j,i) * rmask(j,i,sr) 651 sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) + & 652 usws(j,i) * rmask(j,i,sr) 653 sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) + & 654 vsws(j,i) * rmask(j,i,sr) 655 sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) + & 656 ts(j,i) * rmask(j,i,sr) 657 IF ( humidity ) THEN 658 sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) + & 659 qs(j,i) * rmask(j,i,sr) 671 IF ( surf_def_h(0)%ns >= 1 ) THEN 672 m = surf_def_h(0)%start_index(j,i) 673 sums_l(nzb,pr_palm,tn) = sums_l(nzb,pr_palm,tn) + & 674 surf_def_h(0)%us(m) * rmask(j,i,sr) 675 sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) + & 676 surf_def_h(0)%usws(m) * rmask(j,i,sr) 677 sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) + & 678 surf_def_h(0)%vsws(m) * rmask(j,i,sr) 679 sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) + & 680 surf_def_h(0)%ts(m) * rmask(j,i,sr) 681 IF ( humidity ) THEN 682 sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) + & 683 surf_def_h(0)%qs(m) * rmask(j,i,sr) 684 ENDIF 685 IF ( passive_scalar ) THEN 686 sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) + & 687 surf_def_h(0)%ss(m) * rmask(j,i,sr) 688 ENDIF 660 689 ENDIF 661 IF ( passive_scalar ) THEN 662 sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) + & 663 ss(j,i) * rmask(j,i,sr) 690 IF ( surf_lsm_h%ns >= 1 ) THEN 691 m = surf_lsm_h%start_index(j,i) 692 sums_l(nzb,pr_palm,tn) = sums_l(nzb,pr_palm,tn) + & 693 surf_lsm_h%us(m) * rmask(j,i,sr) 694 sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) + & 695 surf_lsm_h%usws(m) * rmask(j,i,sr) 696 sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) + & 697 surf_lsm_h%vsws(m) * rmask(j,i,sr) 698 sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) + & 699 surf_lsm_h%ts(m) * rmask(j,i,sr) 700 IF ( humidity ) THEN 701 sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) + & 702 surf_lsm_h%qs(m) * rmask(j,i,sr) 703 ENDIF 704 IF ( passive_scalar ) THEN 705 sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) + & 706 surf_lsm_h%ss(m) * rmask(j,i,sr) 707 ENDIF 708 ENDIF 709 IF ( surf_usm_h%ns >= 1 ) THEN 710 m = surf_lsm_h%start_index(j,i) 711 sums_l(nzb,pr_palm,tn) = sums_l(nzb,pr_palm,tn) + & 712 surf_usm_h%us(m) * rmask(j,i,sr) 713 sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) + & 714 surf_usm_h%usws(m) * rmask(j,i,sr) 715 sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) + & 716 surf_usm_h%vsws(m) * rmask(j,i,sr) 717 sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) + & 718 surf_usm_h%ts(m) * rmask(j,i,sr) 719 IF ( humidity ) THEN 720 sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) + & 721 surf_usm_h%qs(m) * rmask(j,i,sr) 722 ENDIF 723 IF ( passive_scalar ) THEN 724 sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) + & 725 surf_usm_h%ss(m) * rmask(j,i,sr) 726 ENDIF 664 727 ENDIF 665 728 ENDDO … … 674 737 DO i = nxl, nxr 675 738 DO j = nys, nyn 676 DO k = nzb_s_inner(j,i), nzt+1 739 DO k = nzb, nzt+1 740 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) ) 741 677 742 u2 = u(k,j,i)**2 678 743 v2 = v(k,j,i)**2 … … 681 746 vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2 682 747 683 sums_l(k,30,tn) = sums_l(k,30,tn) + ust2 * rmask(j,i,sr) 684 sums_l(k,31,tn) = sums_l(k,31,tn) + vst2 * rmask(j,i,sr) 685 sums_l(k,32,tn) = sums_l(k,32,tn) + w2 * rmask(j,i,sr) 748 sums_l(k,30,tn) = sums_l(k,30,tn) + ust2 * rmask(j,i,sr) & 749 * flag 750 sums_l(k,31,tn) = sums_l(k,31,tn) + vst2 * rmask(j,i,sr) & 751 * flag 752 sums_l(k,32,tn) = sums_l(k,32,tn) + w2 * rmask(j,i,sr) & 753 * flag 686 754 ! 687 755 !-- Perturbation energy 688 756 689 757 sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5_wp * & 690 ( ust2 + vst2 + w2 ) * rmask(j,i,sr) 758 ( ust2 + vst2 + w2 ) * rmask(j,i,sr) & 759 * flag 691 760 ENDDO 692 761 ENDDO … … 702 771 DO i = nxl, nxr 703 772 DO j = nys, nyn 704 DO k = nzb_s_inner(j,i), nzt+1 773 DO k = nzb, nzt+1 774 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) ) 775 705 776 w2 = w(k,j,i)**2 706 777 ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2 … … 709 780 710 781 sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn) & 711 + 0.5_wp * ( ust2 + vst2 + w2 ) * rmask(j,i,sr) 782 + 0.5_wp * ( ust2 + vst2 + w2 ) & 783 * rmask(j,i,sr) & 784 * flag 712 785 ENDDO 713 786 ENDDO … … 728 801 !-- However, this implies no error since staggered velocity 729 802 !-- components are zero at the walls and inside buildings. 730 731 DO k = nzb_diff_s_inner(j,i)-1, nzt_diff 803 !-- Flag 23 is used to mask surface fluxes as well as model-top fluxes, 804 !-- which are added further below. 805 DO k = nzb, nzt 806 flag = MERGE( 1.0_wp, 0.0_wp, & 807 BTEST( wall_flags_0(k,j,i), 23 ) ) * & 808 MERGE( 1.0_wp, 0.0_wp, & 809 BTEST( wall_flags_0(k,j,i), 9 ) ) 732 810 ! 733 811 !-- Momentum flux w"u" … … 739 817 ) * rmask(j,i,sr) & 740 818 * rho_air_zw(k) & 741 * momentumflux_output_conversion(k) 819 * momentumflux_output_conversion(k) & 820 * flag 742 821 ! 743 822 !-- Momentum flux w"v" … … 749 828 ) * rmask(j,i,sr) & 750 829 * rho_air_zw(k) & 751 * momentumflux_output_conversion(k) 830 * momentumflux_output_conversion(k) & 831 * flag 752 832 ! 753 833 !-- Heat flux w"pt" … … 757 837 * rho_air_zw(k) & 758 838 * heatflux_output_conversion(k) & 759 * ddzu(k+1) * rmask(j,i,sr) 839 * ddzu(k+1) * rmask(j,i,sr) & 840 * flag 760 841 761 842 … … 766 847 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )& 767 848 * ( sa(k+1,j,i) - sa(k,j,i) ) & 768 * ddzu(k+1) * rmask(j,i,sr) 849 * ddzu(k+1) * rmask(j,i,sr) & 850 * flag 769 851 ENDIF 770 852 … … 777 859 * rho_air_zw(k) & 778 860 * heatflux_output_conversion(k) & 779 * ddzu(k+1) * rmask(j,i,sr) 861 * ddzu(k+1) * rmask(j,i,sr) * flag 780 862 sums_l(k,48,tn) = sums_l(k,48,tn) & 781 863 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )& … … 783 865 * rho_air_zw(k) & 784 866 * waterflux_output_conversion(k)& 785 * ddzu(k+1) * rmask(j,i,sr) 867 * ddzu(k+1) * rmask(j,i,sr) * flag 786 868 787 869 IF ( cloud_physics ) THEN … … 792 874 * rho_air_zw(k) & 793 875 * waterflux_output_conversion(k)& 794 * ddzu(k+1) * rmask(j,i,sr) 876 * ddzu(k+1) * rmask(j,i,sr) * flag 795 877 ENDIF 796 878 ENDIF … … 802 884 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )& 803 885 * ( s(k+1,j,i) - s(k,j,i) ) & 804 * ddzu(k+1) * rmask(j,i,sr) 886 * ddzu(k+1) * rmask(j,i,sr) & 887 * flag 805 888 ENDIF 806 889 … … 810 893 !-- Subgridscale fluxes in the Prandtl layer 811 894 IF ( use_surface_fluxes ) THEN 812 sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + & 895 DO l = 0, 1 896 ki = MERGE( -1, 0, l == 0 ) 897 IF ( surf_def_h(l)%ns >= 1 ) THEN 898 DO m = surf_def_h(l)%start_index(j,i), surf_def_h(l)%end_index(j,i) 899 k = surf_def_h(l)%k(m) 900 901 sums_l(k+ki,12,tn) = sums_l(k+ki,12,tn) + & 902 momentumflux_output_conversion(k+ki) * & 903 surf_def_h(l)%usws(m) * rmask(j,i,sr) ! w"u" 904 sums_l(k+ki,14,tn) = sums_l(k+ki,14,tn) + & 905 momentumflux_output_conversion(k+ki) * & 906 surf_def_h(l)%vsws(m) * rmask(j,i,sr) ! w"v" 907 sums_l(k+ki,16,tn) = sums_l(k+ki,16,tn) + & 908 heatflux_output_conversion(k+ki) * & 909 surf_def_h(l)%shf(m) * rmask(j,i,sr) ! w"pt" 910 sums_l(k+ki,58,tn) = sums_l(k+ki,58,tn) + & 911 0.0_wp * rmask(j,i,sr) ! u"pt" 912 sums_l(k+ki,61,tn) = sums_l(k+ki,61,tn) + & 913 0.0_wp * rmask(j,i,sr) ! v"pt" 914 IF ( ocean ) THEN 915 sums_l(k+ki,65,tn) = sums_l(k+ki,65,tn) + & 916 surf_def_h(l)%sasws(m) * rmask(j,i,sr) ! w"sa" 917 ENDIF 918 IF ( humidity ) THEN 919 sums_l(k+ki,48,tn) = sums_l(k+ki,48,tn) + & 920 waterflux_output_conversion(k+ki) * & 921 surf_def_h(l)%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv") 922 sums_l(k+ki,45,tn) = sums_l(k+ki,45,tn) + ( & 923 ( 1.0_wp + 0.61_wp * q(k+ki,j,i) ) * & 924 surf_def_h(l)%shf(m) + 0.61_wp * pt(k+ki,j,i) * & 925 surf_def_h(l)%qsws(m) ) & 926 * heatflux_output_conversion(k+ki) 927 IF ( cloud_droplets ) THEN 928 sums_l(k+ki,45,tn) = sums_l(k+ki,45,tn) + ( & 929 ( 1.0_wp + 0.61_wp * q(k+ki,j,i) - & 930 ql(k+ki,j,i) ) * surf_def_h(l)%shf(m) + & 931 0.61_wp * pt(k+ki,j,i) * surf_def_h(l)%qsws(m) ) & 932 * heatflux_output_conversion(k+ki) 933 ENDIF 934 IF ( cloud_physics ) THEN 935 ! 936 !-- Formula does not work if ql(k+ki) /= 0.0 937 sums_l(k+ki,51,tn) = sums_l(k+ki,51,tn) + & 938 waterflux_output_conversion(k+ki) * & 939 surf_def_h(l)%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv") 940 ENDIF 941 ENDIF 942 IF ( passive_scalar ) THEN 943 sums_l(k+ki,119,tn) = sums_l(k+ki,119,tn) + & 944 surf_def_h(l)%ssws(m) * rmask(j,i,sr) ! w"s" 945 ENDIF 946 947 ENDDO 948 949 ENDIF 950 ENDDO 951 IF ( surf_lsm_h%ns >= 1 ) THEN 952 m = surf_lsm_h%start_index(j,i) 953 sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + & 813 954 momentumflux_output_conversion(nzb) * & 814 usws(j,i) * rmask(j,i,sr) ! w"u"815 sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &955 surf_lsm_h%usws(m) * rmask(j,i,sr) ! w"u" 956 sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + & 816 957 momentumflux_output_conversion(nzb) * & 817 vsws(j,i) * rmask(j,i,sr) ! w"v"818 sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &958 surf_lsm_h%vsws(m) * rmask(j,i,sr) ! w"v" 959 sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + & 819 960 heatflux_output_conversion(nzb) * & 820 s hf(j,i) * rmask(j,i,sr) ! w"pt"821 sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &961 surf_lsm_h%shf(m) * rmask(j,i,sr) ! w"pt" 962 sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + & 822 963 0.0_wp * rmask(j,i,sr) ! u"pt" 823 sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &964 sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + & 824 965 0.0_wp * rmask(j,i,sr) ! v"pt" 825 IF ( ocean ) THEN826 sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &827 s aswsb(j,i) * rmask(j,i,sr) ! w"sa"828 ENDIF829 IF ( humidity ) THEN830 sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &966 IF ( ocean ) THEN 967 sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + & 968 surf_lsm_h%sasws(m) * rmask(j,i,sr) ! w"sa" 969 ENDIF 970 IF ( humidity ) THEN 971 sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + & 831 972 waterflux_output_conversion(nzb) * & 832 qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv")833 sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + ( &973 surf_lsm_h%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv") 974 sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + ( & 834 975 ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) * & 835 s hf(j,i) + 0.61_wp * pt(nzb,j,i) * &836 qsws(j,i) ) &976 surf_lsm_h%shf(m) + 0.61_wp * pt(nzb,j,i) * & 977 surf_lsm_h%qsws(m) ) & 837 978 * heatflux_output_conversion(nzb) 838 IF ( cloud_droplets ) THEN839 sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + ( &979 IF ( cloud_droplets ) THEN 980 sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + ( & 840 981 ( 1.0_wp + 0.61_wp * q(nzb,j,i) - & 841 ql(nzb,j,i) ) * s hf(j,i) + &842 0.61_wp * pt(nzb,j,i) * qsws(j,i) ) &982 ql(nzb,j,i) ) * surf_lsm_h%shf(m) + & 983 0.61_wp * pt(nzb,j,i) * surf_lsm_h%qsws(m) ) & 843 984 * heatflux_output_conversion(nzb) 985 ENDIF 986 IF ( cloud_physics ) THEN 987 ! 988 !-- Formula does not work if ql(nzb) /= 0.0 989 sums_l(nzb,51,tn) = sums_l(nzb,51,tn) + & 990 waterflux_output_conversion(nzb) * & 991 surf_lsm_h%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv") 992 ENDIF 844 993 ENDIF 845 IF ( cloud_physics ) THEN 846 ! 847 !-- Formula does not work if ql(nzb) /= 0.0 848 sums_l(nzb,51,tn) = sums_l(nzb,51,tn) + & 994 IF ( passive_scalar ) THEN 995 sums_l(nzb,119,tn) = sums_l(nzb,119,tn) + & 996 surf_lsm_h%ssws(m) * rmask(j,i,sr) ! w"s" 997 ENDIF 998 999 1000 ENDIF 1001 IF ( surf_usm_h%ns >= 1 ) THEN 1002 m = surf_usm_h%start_index(j,i) 1003 sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + & 1004 momentumflux_output_conversion(nzb) * & 1005 surf_usm_h%usws(m) * rmask(j,i,sr) ! w"u" 1006 sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + & 1007 momentumflux_output_conversion(nzb) * & 1008 surf_usm_h%vsws(m) * rmask(j,i,sr) ! w"v" 1009 sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + & 1010 heatflux_output_conversion(nzb) * & 1011 surf_usm_h%shf(m) * rmask(j,i,sr) ! w"pt" 1012 sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + & 1013 0.0_wp * rmask(j,i,sr) ! u"pt" 1014 sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + & 1015 0.0_wp * rmask(j,i,sr) ! v"pt" 1016 IF ( ocean ) THEN 1017 sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + & 1018 surf_usm_h%sasws(m) * rmask(j,i,sr) ! w"sa" 1019 ENDIF 1020 IF ( humidity ) THEN 1021 sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + & 1022 waterflux_output_conversion(nzb) * & 1023 surf_usm_h%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv") 1024 sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + ( & 1025 ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) * & 1026 surf_usm_h%shf(m) + 0.61_wp * pt(nzb,j,i) * & 1027 surf_usm_h%qsws(m) ) & 1028 * heatflux_output_conversion(nzb) 1029 IF ( cloud_droplets ) THEN 1030 sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + ( & 1031 ( 1.0_wp + 0.61_wp * q(nzb,j,i) - & 1032 ql(nzb,j,i) ) * surf_usm_h%shf(m) + & 1033 0.61_wp * pt(nzb,j,i) * surf_usm_h%qsws(m) ) & 1034 * heatflux_output_conversion(nzb) 1035 ENDIF 1036 IF ( cloud_physics ) THEN 1037 ! 1038 !-- Formula does not work if ql(nzb) /= 0.0 1039 sums_l(nzb,51,tn) = sums_l(nzb,51,tn) + & 849 1040 waterflux_output_conversion(nzb) * & 850 qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv") 1041 surf_usm_h%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv") 1042 ENDIF 851 1043 ENDIF 852 ENDIF 853 IF ( passive_scalar ) THEN 854 sums_l(nzb,119,tn) = sums_l(nzb,119,tn) + & 855 ssws(j,i) * rmask(j,i,sr) ! w"s" 856 ENDIF 1044 IF ( passive_scalar ) THEN 1045 sums_l(nzb,119,tn) = sums_l(nzb,119,tn) + & 1046 surf_usm_h%ssws(m) * rmask(j,i,sr) ! w"s" 1047 ENDIF 1048 1049 1050 ENDIF 1051 857 1052 ENDIF 858 1053 859 1054 IF ( .NOT. neutral ) THEN 860 sums_l(nzb,114,tn) = sums_l(nzb,114,tn) + & 861 ol(j,i) * rmask(j,i,sr) ! L 862 ENDIF 863 864 865 IF ( land_surface ) THEN 866 sums_l(nzb,93,tn) = sums_l(nzb,93,tn) + ghf_eb(j,i) 867 sums_l(nzb,94,tn) = sums_l(nzb,94,tn) + shf_eb(j,i) 868 sums_l(nzb,95,tn) = sums_l(nzb,95,tn) + qsws_eb(j,i) 869 sums_l(nzb,96,tn) = sums_l(nzb,96,tn) + qsws_liq_eb(j,i) 870 sums_l(nzb,97,tn) = sums_l(nzb,97,tn) + qsws_soil_eb(j,i) 871 sums_l(nzb,98,tn) = sums_l(nzb,98,tn) + qsws_veg_eb(j,i) 872 sums_l(nzb,99,tn) = sums_l(nzb,99,tn) + r_a(j,i) 873 sums_l(nzb,100,tn) = sums_l(nzb,100,tn)+ r_s(j,i) 1055 IF ( surf_def_h(0)%ns >= 1 ) THEN 1056 m = surf_def_h(0)%start_index(j,i) 1057 sums_l(nzb,114,tn) = sums_l(nzb,114,tn) + & 1058 surf_def_h(0)%ol(m) * rmask(j,i,sr) ! L 1059 ENDIF 1060 IF ( surf_lsm_h%ns >= 1 ) THEN 1061 m = surf_lsm_h%start_index(j,i) 1062 sums_l(nzb,114,tn) = sums_l(nzb,114,tn) + & 1063 surf_lsm_h%ol(m) * rmask(j,i,sr) ! L 1064 ENDIF 1065 IF ( surf_usm_h%ns >= 1 ) THEN 1066 m = surf_usm_h%start_index(j,i) 1067 sums_l(nzb,114,tn) = sums_l(nzb,114,tn) + & 1068 surf_usm_h%ol(m) * rmask(j,i,sr) ! L 1069 ENDIF 874 1070 ENDIF 875 1071 … … 893 1089 !-- Subgridscale fluxes at the top surface 894 1090 IF ( use_top_fluxes ) THEN 1091 m = surf_def_h(2)%start_index(j,i) 895 1092 sums_l(nzt:nzt+1,12,tn) = sums_l(nzt:nzt+1,12,tn) + & 896 1093 momentumflux_output_conversion(nzt:nzt+1) * & 897 uswst(j,i) * rmask(j,i,sr) ! w"u"1094 surf_def_h(2)%usws(m) * rmask(j,i,sr) ! w"u" 898 1095 sums_l(nzt:nzt+1,14,tn) = sums_l(nzt:nzt+1,14,tn) + & 899 1096 momentumflux_output_conversion(nzt:nzt+1) * & 900 vswst(j,i) * rmask(j,i,sr) ! w"v"1097 surf_def_h(2)%vsws(m) * rmask(j,i,sr) ! w"v" 901 1098 sums_l(nzt:nzt+1,16,tn) = sums_l(nzt:nzt+1,16,tn) + & 902 1099 heatflux_output_conversion(nzt:nzt+1) * & 903 tswst(j,i) * rmask(j,i,sr) ! w"pt"1100 surf_def_h(2)%shf(m) * rmask(j,i,sr) ! w"pt" 904 1101 sums_l(nzt:nzt+1,58,tn) = sums_l(nzt:nzt+1,58,tn) + & 905 1102 0.0_wp * rmask(j,i,sr) ! u"pt" … … 909 1106 IF ( ocean ) THEN 910 1107 sums_l(nzt,65,tn) = sums_l(nzt,65,tn) + & 911 s aswst(j,i) * rmask(j,i,sr) ! w"sa"1108 surf_def_h(2)%sasws(m) * rmask(j,i,sr) ! w"sa" 912 1109 ENDIF 913 1110 IF ( humidity ) THEN 914 1111 sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + & 915 1112 waterflux_output_conversion(nzt) * & 916 qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")1113 surf_def_h(2)%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv") 917 1114 sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + ( & 918 1115 ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) * & 919 tswst(j,i) + 0.61_wp * pt(nzt,j,i) * & 920 qswst(j,i) ) & 1116 surf_def_h(2)%shf(m) + & 1117 0.61_wp * pt(nzt,j,i) * & 1118 surf_def_h(2)%qsws(m) ) & 921 1119 * heatflux_output_conversion(nzt) 922 1120 IF ( cloud_droplets ) THEN 923 1121 sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + ( & 924 1122 ( 1.0_wp + 0.61_wp * q(nzt,j,i) - & 925 ql(nzt,j,i) ) * tswst(j,i) + & 926 0.61_wp * pt(nzt,j,i) * qswst(j,i) )& 1123 ql(nzt,j,i) ) * & 1124 surf_def_h(2)%shf(m) + & 1125 0.61_wp * pt(nzt,j,i) * & 1126 surf_def_h(2)%qsws(m) )& 927 1127 * heatflux_output_conversion(nzt) 928 1128 ENDIF … … 932 1132 sums_l(nzt,51,tn) = sums_l(nzt,51,tn) + & ! w"q" (w"qv") 933 1133 waterflux_output_conversion(nzt) * & 934 qswst(j,i) * rmask(j,i,sr)1134 surf_def_h(2)%qsws(m) * rmask(j,i,sr) 935 1135 ENDIF 936 1136 ENDIF 937 1137 IF ( passive_scalar ) THEN 938 1138 sums_l(nzt,119,tn) = sums_l(nzt,119,tn) + & 939 s swst(j,i) * rmask(j,i,sr) ! w"s"1139 surf_def_h(2)%ssws(m) * rmask(j,i,sr) ! w"s" 940 1140 ENDIF 941 1141 ENDIF … … 946 1146 !-- ---- speaking the following k-loop would have to be split up and 947 1147 !-- rearranged according to the staggered grid. 948 DO k = nzb_s_inner(j,i), nzt 1148 DO k = nzb, nzt 1149 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) ) 949 1150 ust = 0.5_wp * ( u(k,j,i) - hom(k,1,1,sr) + & 950 1151 u(k+1,j,i) - hom(k+1,1,1,sr) ) … … 956 1157 !-- Higher moments 957 1158 sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)**2 * & 958 rmask(j,i,sr) 1159 rmask(j,i,sr) * flag 959 1160 sums_l(k,36,tn) = sums_l(k,36,tn) + pts**2 * w(k,j,i) * & 960 rmask(j,i,sr) 1161 rmask(j,i,sr) * flag 961 1162 962 1163 ! … … 968 1169 sa(k+1,j,i) - hom(k+1,1,23,sr) ) 969 1170 sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) * & 970 rmask(j,i,sr) 1171 rmask(j,i,sr) * flag 971 1172 ENDIF 972 sums_l(k,64,tn) = sums_l(k,64,tn) + rho_ocean(k,j,i) * 973 rmask(j,i,sr) 1173 sums_l(k,64,tn) = sums_l(k,64,tn) + rho_ocean(k,j,i) * & 1174 rmask(j,i,sr) * flag 974 1175 sums_l(k,71,tn) = sums_l(k,71,tn) + prho(k,j,i) * & 975 rmask(j,i,sr) 1176 rmask(j,i,sr) * flag 976 1177 ENDIF 977 1178 … … 985 1186 sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * & 986 1187 heatflux_output_conversion(k) * & 987 rmask(j,i,sr) 988 sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * rmask(j,i,sr) 1188 rmask(j,i,sr) * flag 1189 sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * rmask(j,i,sr) & 1190 * flag 989 1191 990 1192 IF ( .NOT. cloud_droplets ) THEN … … 996 1198 sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) * & 997 1199 waterflux_output_conversion(k) * & 998 rmask(j,i,sr) 1200 rmask(j,i,sr) * & 1201 flag 999 1202 sums_l(k,75,tn) = sums_l(k,75,tn) + qc(k,j,i) * & 1000 rmask(j,i,sr) 1203 rmask(j,i,sr) * & 1204 flag 1001 1205 sums_l(k,76,tn) = sums_l(k,76,tn) + prr(k,j,i) * & 1002 rmask(j,i,sr) 1206 rmask(j,i,sr) * & 1207 flag 1003 1208 IF ( microphysics_seifert ) THEN 1004 1209 sums_l(k,73,tn) = sums_l(k,73,tn) + nr(k,j,i) * & 1005 rmask(j,i,sr) 1210 rmask(j,i,sr) *& 1211 flag 1006 1212 sums_l(k,74,tn) = sums_l(k,74,tn) + qr(k,j,i) * & 1007 rmask(j,i,sr) 1213 rmask(j,i,sr) *& 1214 flag 1008 1215 ENDIF 1009 1216 ENDIF … … 1015 1222 sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * & 1016 1223 heatflux_output_conversion(k) * & 1017 rmask(j,i,sr) 1224 rmask(j,i,sr) * & 1225 flag 1018 1226 ELSE IF ( ws_scheme_sca .AND. sr == 0 ) THEN 1019 1227 sums_l(k,46,tn) = ( ( 1.0_wp + 0.61_wp * & … … 1022 1230 0.61_wp * hom(k,1,4,sr) * & 1023 1231 sums_l(k,49,tn) & 1024 ) * heatflux_output_conversion(k) 1232 ) * heatflux_output_conversion(k) * & 1233 flag 1025 1234 END IF 1026 1235 END IF … … 1033 1242 s(k+1,j,i) - hom(k+1,1,117,sr) ) 1034 1243 sums_l(k,116,tn) = sums_l(k,116,tn) + pts * w(k,j,i) * & 1035 rmask(j,i,sr) 1244 rmask(j,i,sr) * flag 1036 1245 ENDIF 1037 1246 … … 1042 1251 ( ust**2 + vst**2 + w(k,j,i)**2 ) & 1043 1252 * momentumflux_output_conversion(k) & 1044 * rmask(j,i,sr)1253 * rmask(j,i,sr) * flag 1045 1254 ENDDO 1046 1255 ENDDO 1047 1256 ENDDO 1257 !$OMP END PARALLEL 1258 ! 1259 !-- Treat land-surface quantities according to new wall model structure. 1260 IF ( land_surface ) THEN 1261 tn = 0 1262 !$OMP PARALLEL PRIVATE( i, j, m, tn ) 1263 !$ tn = omp_get_thread_num() 1264 !$OMP DO 1265 DO m = 1, surf_lsm_h%ns 1266 i = surf_lsm_h%i(m) 1267 j = surf_lsm_h%j(m) 1268 1269 IF ( i >= nxl .AND. i <= nxr .AND. & 1270 j >= nys .AND. j <= nyn ) THEN 1271 sums_l(nzb,93,tn) = sums_l(nzb,93,tn) + surf_lsm_h%ghf_eb(m) 1272 sums_l(nzb,94,tn) = sums_l(nzb,94,tn) + surf_lsm_h%shf_eb(m) 1273 sums_l(nzb,95,tn) = sums_l(nzb,95,tn) + surf_lsm_h%qsws_eb(m) 1274 sums_l(nzb,96,tn) = sums_l(nzb,96,tn) + surf_lsm_h%qsws_liq_eb(m) 1275 sums_l(nzb,97,tn) = sums_l(nzb,97,tn) + surf_lsm_h%qsws_soil_eb(m) 1276 sums_l(nzb,98,tn) = sums_l(nzb,98,tn) + surf_lsm_h%qsws_veg_eb(m) 1277 sums_l(nzb,99,tn) = sums_l(nzb,99,tn) + surf_lsm_h%r_a(m) 1278 sums_l(nzb,100,tn) = sums_l(nzb,100,tn)+ surf_lsm_h%r_s(m) 1279 ENDIF 1280 ENDDO 1281 !$OMP END PARALLEL 1282 1283 tn = 0 1284 !$OMP PARALLEL PRIVATE( i, j, k, m, tn ) 1285 !$ tn = omp_get_thread_num() 1286 !$OMP DO 1287 DO m = 1, surf_lsm_h%ns 1288 1289 i = surf_lsm_h%i(m) 1290 j = surf_lsm_h%j(m) 1291 1292 IF ( i >= nxl .AND. i <= nxr .AND. & 1293 j >= nys .AND. j <= nyn ) THEN 1294 1295 DO k = nzb_soil, nzt_soil 1296 sums_l(k,89,tn) = sums_l(k,89,tn) + t_soil_h%var_2d(k,m) & 1297 * rmask(j,i,sr) 1298 sums_l(k,91,tn) = sums_l(k,91,tn) + m_soil_h%var_2d(k,m) & 1299 * rmask(j,i,sr) 1300 ENDDO 1301 ENDIF 1302 ENDDO 1303 !$OMP END PARALLEL 1304 ENDIF 1048 1305 ! 1049 1306 !-- For speed optimization fluxes which have been computed in part directly 1050 1307 !-- inside the WS advection routines are treated seperatly 1051 1308 !-- Momentum fluxes first: 1309 1310 tn = 0 1311 !$OMP PARALLEL PRIVATE( i, j, k, tn, flag ) 1312 !$ tn = omp_get_thread_num() 1052 1313 IF ( .NOT. ws_scheme_mom .OR. sr /= 0 ) THEN 1053 !$OMP DO 1054 DO i = nxl, nxr 1055 DO j = nys, nyn 1056 DO k = nzb_diff_s_inner(j,i)-1, nzt_diff 1057 ust = 0.5_wp * ( u(k,j,i) - hom(k,1,1,sr) + & 1058 u(k+1,j,i) - hom(k+1,1,1,sr) ) 1059 vst = 0.5_wp * ( v(k,j,i) - hom(k,1,2,sr) + & 1060 v(k+1,j,i) - hom(k+1,1,2,sr) ) 1061 ! 1062 !-- Momentum flux w*u* 1063 sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5_wp * & 1064 ( w(k,j,i-1) + w(k,j,i) ) & 1065 * momentumflux_output_conversion(k) & 1066 * ust * rmask(j,i,sr) 1067 ! 1068 !-- Momentum flux w*v* 1069 sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5_wp * & 1070 ( w(k,j-1,i) + w(k,j,i) ) & 1071 * momentumflux_output_conversion(k) & 1072 * vst * rmask(j,i,sr) 1073 ENDDO 1074 ENDDO 1075 ENDDO 1314 !$OMP DO 1315 DO i = nxl, nxr 1316 DO j = nys, nyn 1317 DO k = nzb, nzt 1318 ! 1319 !-- Flag 23 is used to mask surface fluxes as well as model-top 1320 !-- fluxes, which are added further below. 1321 flag = MERGE( 1.0_wp, 0.0_wp, & 1322 BTEST( wall_flags_0(k,j,i), 23 ) ) * & 1323 MERGE( 1.0_wp, 0.0_wp, & 1324 BTEST( wall_flags_0(k,j,i), 9 ) ) 1325 1326 ust = 0.5_wp * ( u(k,j,i) - hom(k,1,1,sr) + & 1327 u(k+1,j,i) - hom(k+1,1,1,sr) ) 1328 vst = 0.5_wp * ( v(k,j,i) - hom(k,1,2,sr) + & 1329 v(k+1,j,i) - hom(k+1,1,2,sr) ) 1330 ! 1331 !-- Momentum flux w*u* 1332 sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5_wp * & 1333 ( w(k,j,i-1) + w(k,j,i) ) & 1334 * momentumflux_output_conversion(k) & 1335 * ust * rmask(j,i,sr) & 1336 * flag 1337 ! 1338 !-- Momentum flux w*v* 1339 sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5_wp * & 1340 ( w(k,j-1,i) + w(k,j,i) ) & 1341 * momentumflux_output_conversion(k) & 1342 * vst * rmask(j,i,sr) & 1343 * flag 1344 ENDDO 1345 ENDDO 1346 ENDDO 1076 1347 1077 1348 ENDIF 1078 1349 IF ( .NOT. ws_scheme_sca .OR. sr /= 0 ) THEN 1079 !$OMP DO 1080 DO i = nxl, nxr 1081 DO j = nys, nyn 1082 DO k = nzb_diff_s_inner(j,i)-1, nzt_diff 1083 ! 1084 !-- Vertical heat flux 1085 sums_l(k,17,tn) = sums_l(k,17,tn) + 0.5_wp * & 1350 !$OMP DO 1351 DO i = nxl, nxr 1352 DO j = nys, nyn 1353 DO k = nzb, nzt 1354 flag = MERGE( 1.0_wp, 0.0_wp, & 1355 BTEST( wall_flags_0(k,j,i), 23 ) ) * & 1356 MERGE( 1.0_wp, 0.0_wp, & 1357 BTEST( wall_flags_0(k,j,i), 9 ) ) 1358 ! 1359 !-- Vertical heat flux 1360 sums_l(k,17,tn) = sums_l(k,17,tn) + 0.5_wp * & 1086 1361 ( pt(k,j,i) - hom(k,1,4,sr) + & 1087 1362 pt(k+1,j,i) - hom(k+1,1,4,sr) ) & 1088 1363 * heatflux_output_conversion(k) & 1089 * w(k,j,i) * rmask(j,i,sr) 1090 IF ( humidity ) THEN1091 pts = 0.5_wp * ( q(k,j,i) - hom(k,1,41,sr) +&1364 * w(k,j,i) * rmask(j,i,sr) * flag 1365 IF ( humidity ) THEN 1366 pts = 0.5_wp * ( q(k,j,i) - hom(k,1,41,sr) + & 1092 1367 q(k+1,j,i) - hom(k+1,1,41,sr) ) 1093 sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) *&1368 sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * & 1094 1369 waterflux_output_conversion(k) * & 1095 rmask(j,i,sr) 1096 ENDIF1097 IF ( passive_scalar ) THEN1098 pts = 0.5_wp * ( s(k,j,i) - hom(k,1,117,sr) +&1370 rmask(j,i,sr) * flag 1371 ENDIF 1372 IF ( passive_scalar ) THEN 1373 pts = 0.5_wp * ( s(k,j,i) - hom(k,1,117,sr) + & 1099 1374 s(k+1,j,i) - hom(k+1,1,117,sr) ) 1100 sums_l(k,116,tn) = sums_l(k,116,tn) + pts * w(k,j,i) *&1101 rmask(j,i,sr) 1102 ENDIF1103 ENDDO1104 ENDDO1105 ENDDO1375 sums_l(k,116,tn) = sums_l(k,116,tn) + pts * w(k,j,i) * & 1376 rmask(j,i,sr) * flag 1377 ENDIF 1378 ENDDO 1379 ENDDO 1380 ENDDO 1106 1381 1107 1382 ENDIF … … 1126 1401 DO i = nxl, nxr 1127 1402 DO j = nys, nyn 1128 DO k = nzb_s_inner(j,i)+1, nzt 1403 DO k = nzb+1, nzt 1404 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 1129 1405 1130 1406 sums_ll(k,1) = sums_ll(k,1) + 0.5_wp * w(k,j,i) * ( & … … 1133 1409 + ( 0.25_wp * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i) ) & 1134 1410 - 0.5_wp * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) ) )**2& 1135 + w(k,j,i)**2 ) 1411 + w(k,j,i)**2 ) * flag 1136 1412 1137 1413 sums_ll(k,2) = sums_ll(k,2) + 0.5_wp * w(k,j,i) & 1138 * ( p(k,j,i) + p(k+1,j,i) ) 1414 * ( p(k,j,i) + p(k+1,j,i) ) * flag 1139 1415 1140 1416 ENDDO … … 1164 1440 DO i = nxl, nxr 1165 1441 DO j = nys, nyn 1166 DO k = nzb_s_inner(j,i)+1, nzt 1442 DO k = nzb+1, nzt 1443 1444 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 1167 1445 1168 1446 sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5_wp * ( & 1169 1447 (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) & 1170 1448 - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k) & 1171 ) * ddzw(k) 1449 ) * ddzw(k) & 1450 * flag 1172 1451 1173 1452 sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5_wp * ( & 1174 1453 (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) & 1175 ) 1454 ) * flag 1176 1455 1177 1456 ENDDO … … 1191 1470 DO i = nxl, nxr 1192 1471 DO j = nys, nyn 1193 DO k = nzb_s_inner(j,i)+1, nzt 1472 DO k = nzb+1, nzt 1473 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 1194 1474 ! 1195 1475 !-- Subgrid horizontal heat fluxes u"pt", v"pt" … … 1199 1479 * rho_air_zw(k) & 1200 1480 * heatflux_output_conversion(k) & 1201 * ddx * rmask(j,i,sr) 1481 * ddx * rmask(j,i,sr) * flag 1202 1482 sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp * & 1203 1483 ( kh(k,j,i) + kh(k,j-1,i) ) & … … 1205 1485 * rho_air_zw(k) & 1206 1486 * heatflux_output_conversion(k) & 1207 * ddy * rmask(j,i,sr) 1487 * ddy * rmask(j,i,sr) * flag 1208 1488 ! 1209 1489 !-- Resolved horizontal heat fluxes u*pt*, v*pt* … … 1212 1492 * 0.5_wp * ( pt(k,j,i-1) - hom(k,1,4,sr) + & 1213 1493 pt(k,j,i) - hom(k,1,4,sr) ) & 1214 * heatflux_output_conversion(k) 1494 * heatflux_output_conversion(k) & 1495 * flag 1215 1496 pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) + & 1216 1497 pt(k,j,i) - hom(k,1,4,sr) ) … … 1219 1500 * 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) + & 1220 1501 pt(k,j,i) - hom(k,1,4,sr) ) & 1221 * heatflux_output_conversion(k) 1502 * heatflux_output_conversion(k) & 1503 * flag 1222 1504 ENDDO 1223 1505 ENDDO … … 1279 1561 ENDIF 1280 1562 1563 tn = 0 1281 1564 !$OMP PARALLEL PRIVATE( i, j, k, tn ) 1282 !$ tn = omp_get_thread_num() 1283 IF ( land_surface ) THEN 1284 !$OMP DO 1285 DO i = nxl, nxr 1286 DO j = nys, nyn 1287 DO k = nzb_soil, nzt_soil 1288 sums_l(k,89,tn) = sums_l(k,89,tn) + t_soil(k,j,i) & 1289 * rmask(j,i,sr) 1290 sums_l(k,91,tn) = sums_l(k,91,tn) + m_soil(k,j,i) & 1291 * rmask(j,i,sr) 1292 ENDDO 1293 ENDDO 1294 ENDDO 1295 ENDIF 1296 1565 !$ tn = omp_get_thread_num() 1297 1566 IF ( radiation .AND. radiation_scheme == 'rrtmg' ) THEN 1298 1567 !$OMP DO 1299 1568 DO i = nxl, nxr 1300 1569 DO j = nys, nyn 1301 DO k = nzb_s_inner(j,i)+1, nzt+1 1570 DO k = nzb+1, nzt+1 1571 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 1572 1302 1573 sums_l(k,102,tn) = sums_l(k,102,tn) + rad_lw_in(k,j,i) & 1303 * rmask(j,i,sr) 1574 * rmask(j,i,sr) * flag 1304 1575 sums_l(k,103,tn) = sums_l(k,103,tn) + rad_lw_out(k,j,i) & 1305 * rmask(j,i,sr) 1576 * rmask(j,i,sr) * flag 1306 1577 sums_l(k,104,tn) = sums_l(k,104,tn) + rad_sw_in(k,j,i) & 1307 * rmask(j,i,sr) 1578 * rmask(j,i,sr) * flag 1308 1579 sums_l(k,105,tn) = sums_l(k,105,tn) + rad_sw_out(k,j,i) & 1309 * rmask(j,i,sr) 1580 * rmask(j,i,sr) * flag 1310 1581 sums_l(k,106,tn) = sums_l(k,106,tn) + rad_lw_cs_hr(k,j,i) & 1311 * rmask(j,i,sr) 1582 * rmask(j,i,sr) * flag 1312 1583 sums_l(k,107,tn) = sums_l(k,107,tn) + rad_lw_hr(k,j,i) & 1313 * rmask(j,i,sr) 1584 * rmask(j,i,sr) * flag 1314 1585 sums_l(k,108,tn) = sums_l(k,108,tn) + rad_sw_cs_hr(k,j,i) & 1315 * rmask(j,i,sr) 1586 * rmask(j,i,sr) * flag 1316 1587 sums_l(k,109,tn) = sums_l(k,109,tn) + rad_sw_hr(k,j,i) & 1317 * rmask(j,i,sr) 1588 * rmask(j,i,sr) * flag 1318 1589 ENDDO 1319 1590 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.
