- Timestamp:
- Aug 17, 2011 2:13:26 PM (13 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 3 edited
-
advec_ws.f90 (modified) (48 diffs)
-
modules.f90 (modified) (3 diffs)
-
prognostic_equations.f90 (modified) (20 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/advec_ws.f90
r714 r736 3 3 !-----------------------------------------------------------------------------! 4 4 ! Current revisions: 5 ! ----------------- 5 ! ------------------ 6 ! Bugfix concerning OpenMP parallelization. i_omp introduced, because first 7 ! index where fluxes on left side have to be calculated explicitly is 8 ! different on each thread. Furthermore the swapping of fluxes is now 9 ! thread-safe by adding an additional dimension. 6 10 ! 7 11 ! Former revisions: … … 43 47 ! A further routine local_diss_ij is available (next weeks) and is used if a 44 48 ! control of dissipative fluxes is desired. 45 ! In case of vertical grid stretching the order of advection scheme is46 ! degraded. This is also realized for the ocean where the stretching is47 ! downwards.48 49 ! 49 50 ! OUTSTANDING: - Dirichlet and Radiation boundary conditions for … … 102 103 USE control_parameters 103 104 USE indices 105 USE pegrid 104 106 USE statistics 105 107 … … 166 168 IF ( ws_scheme_mom ) THEN 167 169 168 ALLOCATE( flux_s_u(nzb+1:nzt), flux_s_v(nzb+1:nzt), & 169 flux_s_w(nzb+1:nzt), diss_s_u(nzb+1:nzt), & 170 diss_s_v(nzb+1:nzt), diss_s_w(nzb+1:nzt) ) 171 ALLOCATE( flux_l_u(nzb+1:nzt,nys:nyn), & 172 flux_l_v(nzb+1:nzt,nys:nyn), & 173 flux_l_w(nzb+1:nzt,nys:nyn), & 174 diss_l_u(nzb+1:nzt,nys:nyn), & 175 diss_l_v(nzb+1:nzt,nys:nyn), & 176 diss_l_w(nzb+1:nzt,nys:nyn) ) 170 ALLOCATE( flux_s_u(nzb+1:nzt,0:threads_per_task-1), & 171 flux_s_v(nzb+1:nzt,0:threads_per_task-1), & 172 flux_s_w(nzb+1:nzt,0:threads_per_task-1), & 173 diss_s_u(nzb+1:nzt,0:threads_per_task-1), & 174 diss_s_v(nzb+1:nzt,0:threads_per_task-1), & 175 diss_s_w(nzb+1:nzt,0:threads_per_task-1) ) 176 ALLOCATE( flux_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & 177 flux_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & 178 flux_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & 179 diss_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & 180 diss_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & 181 diss_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ) 177 182 178 183 ENDIF … … 180 185 IF ( ws_scheme_sca ) THEN 181 186 182 ALLOCATE( flux_s_pt(nzb+1:nzt), flux_s_e(nzb+1:nzt), & 183 diss_s_pt(nzb+1:nzt), diss_s_e(nzb+1:nzt) ) 184 ALLOCATE( flux_l_pt(nzb+1:nzt,nys:nyn), & 185 flux_l_e(nzb+1:nzt,nys:nyn), & 186 diss_l_pt(nzb+1:nzt,nys:nyn), & 187 diss_l_e(nzb+1:nzt,nys:nyn) ) 187 ALLOCATE( flux_s_pt(nzb+1:nzt,0:threads_per_task-1), & 188 flux_s_e(nzb+1:nzt,0:threads_per_task-1), & 189 diss_s_pt(nzb+1:nzt,0:threads_per_task-1), & 190 diss_s_e(nzb+1:nzt,0:threads_per_task-1) ) 191 ALLOCATE( flux_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & 192 flux_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & 193 diss_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & 194 diss_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ) 188 195 189 196 IF ( humidity .OR. passive_scalar ) THEN 190 ALLOCATE( flux_s_q(nzb+1:nzt), diss_s_q(nzb+1:nzt) ) 191 ALLOCATE( flux_l_q(nzb+1:nzt,nys:nyn), & 192 diss_l_q(nzb+1:nzt,nys:nyn) ) 197 ALLOCATE( flux_s_q(nzb+1:nzt,0:threads_per_task-1), & 198 diss_s_q(nzb+1:nzt,0:threads_per_task-1) ) 199 ALLOCATE( flux_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & 200 diss_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ) 193 201 ENDIF 194 202 195 203 IF ( ocean ) THEN 196 ALLOCATE( flux_s_sa(nzb+1:nzt), diss_s_sa(nzb+1:nzt) ) 197 ALLOCATE( flux_l_sa(nzb+1:nzt,nys:nyn), & 198 diss_l_sa(nzb+1:nzt,nys:nyn) ) 204 ALLOCATE( flux_s_sa(nzb+1:nzt,0:threads_per_task-1), & 205 diss_s_sa(nzb+1:nzt,0:threads_per_task-1) ) 206 ALLOCATE( flux_l_sa(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & 207 diss_l_sa(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ) 199 208 ENDIF 200 209 … … 267 276 SUBROUTINE advec_s_ws_ij( i, j, sk, sk_char,swap_flux_y_local, & 268 277 swap_diss_y_local, swap_flux_x_local, & 269 swap_diss_x_local 278 swap_diss_x_local, i_omp, tn ) 270 279 271 280 USE arrays_3d … … 274 283 USE grid_variables 275 284 USE indices 285 USE pegrid 276 286 USE statistics 277 287 278 288 IMPLICIT NONE 279 289 280 INTEGER :: i, j, k290 INTEGER :: i, i_omp, j, k, tn 281 291 LOGICAL :: degraded_l, degraded_s 282 292 REAL :: flux_d, diss_d, u_comp, v_comp … … 284 294 REAL, DIMENSION(nzb:nzt+1) :: flux_t, diss_t, flux_r, diss_r, & 285 295 flux_n, diss_n 286 REAL, DIMENSION(nzb+1:nzt) :: swap_flux_y_local, & 287 swap_diss_y_local 288 REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local, & 289 swap_diss_x_local 296 REAL, DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: swap_flux_y_local, & 297 swap_diss_y_local 298 REAL, DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: & 299 swap_flux_x_local, & 300 swap_diss_x_local 290 301 CHARACTER (LEN = *), INTENT(IN) :: sk_char 291 302 … … 368 379 !-- Compute leftside fluxes for the left boundary of PE domain 369 380 u_comp = u(k,j,i) - u_gtrans 370 swap_flux_x_local(k,j ) = u_comp * ( &381 swap_flux_x_local(k,j,tn) = u_comp * ( & 371 382 sk(k,j,i) + sk(k,j,i-1) ) * 0.5 372 swap_diss_x_local(k,j ) = diss_2nd( sk(k,j,i+2),sk(k,j,i+1), &383 swap_diss_x_local(k,j,tn) = diss_2nd( sk(k,j,i+2),sk(k,j,i+1), & 373 384 sk(k,j,i), sk(k,j,i-1), & 374 385 sk(k,j,i-1), u_comp, & … … 402 413 !-- Compute southside fluxes for the south boundary of PE domain 403 414 v_comp = v(k,j,i) - v_gtrans 404 swap_flux_y_local(k ) = v_comp * &415 swap_flux_y_local(k,tn) = v_comp * & 405 416 ( sk(k,j,i) + sk(k,j-1,i) ) * 0.5 406 swap_diss_y_local(k ) = diss_2nd( sk(k,j+2,i), sk(k,j+1,i), &417 swap_diss_y_local(k,tn) = diss_2nd( sk(k,j+2,i), sk(k,j+1,i), & 407 418 sk(k,j,i), sk(k,j-1,i), & 408 419 sk(k,j-1,i), v_comp, & … … 477 488 ! 478 489 !-- Compute left- and southside fluxes of the respective PE bounds. 479 IF ( j == nys .AND. ( .NOT. degraded_s )) THEN490 IF ( j == nys .AND. .NOT. degraded_s ) THEN 480 491 481 492 DO k = nzb_s_inner(j,i)+1, nzt 482 493 v_comp = v(k,j,i) - v_gtrans 483 swap_flux_y_local(k ) = v_comp * ( &494 swap_flux_y_local(k,tn) = v_comp * ( & 484 495 37.0 * ( sk(k,j,i) + sk(k,j-1,i) ) & 485 496 - 8.0 * ( sk(k,j+1,i) + sk(k,j-2,i) ) & 486 497 + ( sk(k,j+2,i) + sk(k,j-3,i) ) & 487 498 ) * adv_sca_5 488 swap_diss_y_local(k ) = -ABS( v_comp ) * ( &499 swap_diss_y_local(k,tn) = -ABS( v_comp ) * ( & 489 500 10.0 * ( sk(k,j,i) - sk(k,j-1,i) ) & 490 501 - 5.0 * ( sk(k,j+1,i) - sk(k,j-2,i) ) & … … 495 506 ENDIF 496 507 497 IF ( i == nxl.AND. .NOT. degraded_l ) THEN508 IF ( i == i_omp .AND. .NOT. degraded_l ) THEN 498 509 499 510 DO k = nzb_s_inner(j,i)+1, nzt 500 511 u_comp = u(k,j,i) - u_gtrans 501 swap_flux_x_local(k,j ) = u_comp * ( &512 swap_flux_x_local(k,j,tn) = u_comp * ( & 502 513 37.0 * ( sk(k,j,i) + sk(k,j,i-1) ) & 503 514 - 8.0 * ( sk(k,j,i+1) + sk(k,j,i-2) ) & 504 515 + ( sk(k,j,i+2) + sk(k,j,i-3) ) & 505 516 ) * adv_sca_5 506 swap_diss_x_local(k,j ) = -ABS( u_comp ) * ( &517 swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * ( & 507 518 10.0 * ( sk(k,j,i) - sk(k,j,i-1) ) & 508 519 - 5.0 * ( sk(k,j,i+1) - sk(k,j,i-2) ) & … … 518 529 519 530 tend(k,j,i) = tend(k,j,i) - ( & 520 ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j ) - &521 swap_diss_x_local(k,j ) ) * ddx &522 + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k ) - &523 swap_diss_y_local(k ) ) * ddy &531 ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - & 532 swap_diss_x_local(k,j,tn) ) * ddx & 533 + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn) - & 534 swap_diss_y_local(k,tn) ) * ddy & 524 535 ) 525 536 526 swap_flux_y_local(k ) = flux_n(k)527 swap_diss_y_local(k ) = diss_n(k)528 swap_flux_x_local(k,j ) = flux_r(k)529 swap_diss_x_local(k,j ) = diss_r(k)537 swap_flux_y_local(k,tn) = flux_n(k) 538 swap_diss_y_local(k,tn) = diss_n(k) 539 swap_flux_x_local(k,j,tn) = flux_r(k) 540 swap_diss_x_local(k,j,tn) = diss_r(k) 530 541 531 542 ENDDO … … 661 672 ! Advection of u-component - Call for grid point i,j 662 673 !------------------------------------------------------------------------------! 663 SUBROUTINE advec_u_ws_ij( i, j )674 SUBROUTINE advec_u_ws_ij( i, j, i_omp, tn ) 664 675 665 676 USE arrays_3d … … 672 683 IMPLICIT NONE 673 684 674 INTEGER :: i, j, k685 INTEGER :: i, i_omp, j, k, tn 675 686 LOGICAL :: degraded_l, degraded_s 676 687 REAL :: gu, gv, flux_d, diss_d, u_comp_l, v_comp, w_comp 677 REAL, DIMENSION(nzb:nzt+1) 678 688 REAL, DIMENSION(nzb:nzt+1) :: flux_t, diss_t, flux_r, diss_r, & 689 flux_n, diss_n, u_comp 679 690 680 691 degraded_l = .FALSE. … … 696 707 - ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3 697 708 diss_r(k) = - ABS( u_comp(k) - gu ) * ( & 698 3.0 * ( u(k,j,i+1) - u(k,j,i) ) & 709 3.0 * ( u(k,j,i+1) - u(k,j,i) ) & 699 710 - ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3 700 711 ENDDO … … 717 728 - ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3 718 729 diss_n(k) = - ABS( v_comp ) * ( & 719 3.0 * ( u(k,j+1,i) - u(k,j,i) ) & 730 3.0 * ( u(k,j+1,i) - u(k,j,i) ) & 720 731 - ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3 721 732 ENDDO … … 735 746 !-- Compute leftside fluxes for the left boundary of PE domain 736 747 u_comp(k) = u(k,j,i) + u(k,j,i-1) - gu 737 flux_l_u(k,j ) = u_comp(k) * ( u(k,j,i) + u(k,j,i-1) ) * 0.25738 diss_l_u(k,j ) = diss_2nd( u(k,j,i+2), u(k,j,i+1), u(k,j,i), &748 flux_l_u(k,j,tn) = u_comp(k) * ( u(k,j,i) + u(k,j,i-1) ) * 0.25 749 diss_l_u(k,j,tn) = diss_2nd( u(k,j,i+2), u(k,j,i+1), u(k,j,i), & 739 750 u(k,j,i-1), u(k,j,i-1), u_comp(k),& 740 751 0.25, ddx ) … … 745 756 - ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3 746 757 diss_r(k) = - ABS( u_comp(k) -gu ) * ( & 747 3.0 * ( u(k,j,i+1) - u(k,j,i) ) & 758 3.0 * ( u(k,j,i+1) - u(k,j,i) ) & 748 759 - ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3 749 760 ENDDO … … 757 768 - ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3 758 769 diss_n(k) = - ABS( v_comp ) * ( & 759 3.0 * ( u(k,j+1,i) - u(k,j,i) ) & 770 3.0 * ( u(k,j+1,i) - u(k,j,i) ) & 760 771 - ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3 761 772 ENDDO … … 764 775 DO k = nzb_u_inner(j,i)+1, nzt 765 776 ! 766 !-- Compute southside fluxes for the south boundary of PE domain 777 !-- Compute southside fluxes for the south boundary of PE domain 767 778 v_comp = v(k,j,i) + v(k,j,i-1) - gv 768 flux_s_u(k ) = v_comp * ( u(k,j,i) + u(k,j-1,i) ) * 0.25769 diss_s_u(k ) = diss_2nd( u(k,j+2,i), u(k,j+1,i), u(k,j,i), &779 flux_s_u(k,tn) = v_comp * ( u(k,j,i) + u(k,j-1,i) ) * 0.25 780 diss_s_u(k,tn) = diss_2nd( u(k,j+2,i), u(k,j+1,i), u(k,j,i), & 770 781 u(k,j-1,i), u(k,j-1,i), v_comp, & 771 782 0.25, ddy ) … … 776 787 - ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3 777 788 diss_n(k) = - ABS( v_comp ) * ( & 778 3.0 * ( u(k,j+1,i) - u(k,j,i) ) & 789 3.0 * ( u(k,j+1,i) - u(k,j,i) ) & 779 790 - ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3 780 791 ENDDO … … 819 830 ELSE 820 831 ! 821 !-- Compute the fifth order fluxes for the interior of PE domain. 832 !-- Compute the fifth order fluxes for the interior of PE domain. 822 833 DO k = nzb_u_inner(j,i)+1, nzt 823 834 u_comp(k) = u(k,j,i+1) + u(k,j,i) … … 844 855 ENDIF 845 856 ! 846 !-- Compute left- and southside fluxes for the respective boundary of PE 857 !-- Compute left- and southside fluxes for the respective boundary of PE 847 858 IF ( j == nys .AND. ( .NOT. degraded_s ) ) THEN 848 859 849 860 DO k = nzb_u_inner(j,i)+1, nzt 850 861 v_comp = v(k,j,i) + v(k,j,i-1) - gv 851 flux_s_u(k ) = v_comp * ( &862 flux_s_u(k,tn) = v_comp * ( & 852 863 37.0 * ( u(k,j,i) + u(k,j-1,i) ) & 853 864 - 8.0 * ( u(k,j+1,i) + u(k,j-2,i) ) & 854 865 + ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5 855 diss_s_u(k ) = - ABS(v_comp) * ( &866 diss_s_u(k,tn) = - ABS(v_comp) * ( & 856 867 10.0 * ( u(k,j,i) - u(k,j-1,i) ) & 857 868 - 5.0 * ( u(k,j+1,i) - u(k,j-2,i) ) & … … 861 872 ENDIF 862 873 863 IF ( i == nxlu.AND. ( .NOT. degraded_l ) ) THEN874 IF ( i == i_omp .AND. ( .NOT. degraded_l ) ) THEN 864 875 865 876 DO k = nzb_u_inner(j,i)+1, nzt 866 877 u_comp_l = u(k,j,i)+u(k,j,i-1)-gu 867 flux_l_u(k,j ) = u_comp_l * ( &878 flux_l_u(k,j,tn) = u_comp_l * ( & 868 879 37.0 * ( u(k,j,i) + u(k,j,i-1) ) & 869 880 - 8.0 * ( u(k,j,i+1) + u(k,j,i-2) ) & 870 881 + ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5 871 diss_l_u(k,j ) = - ABS(u_comp_l) * ( &882 diss_l_u(k,j,tn) = - ABS(u_comp_l) * ( & 872 883 10.0 * ( u(k,j,i) - u(k,j,i-1) ) & 873 884 - 5.0 * ( u(k,j,i+1) - u(k,j,i-2) ) & … … 881 892 tend(k,j,i) = tend(k,j,i) - ( & 882 893 ( flux_r(k) + diss_r(k) & 883 - flux_l_u(k,j ) - diss_l_u(k,j) ) * ddx&894 - flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx & 884 895 + ( flux_n(k) + diss_n(k) & 885 - flux_s_u(k ) - diss_s_u(k) ) * ddy )886 887 flux_l_u(k,j ) = flux_r(k)888 diss_l_u(k,j ) = diss_r(k)889 flux_s_u(k ) = flux_n(k)890 diss_s_u(k ) = diss_n(k)896 - flux_s_u(k,tn) - diss_s_u(k,tn) ) * ddy ) 897 898 flux_l_u(k,j,tn) = flux_r(k) 899 diss_l_u(k,j,tn) = diss_r(k) 900 flux_s_u(k,tn) = flux_n(k) 901 diss_s_u(k,tn) = diss_n(k) 891 902 ! 892 903 !-- Statistical Evaluation of u'u'. The factor has to be applied for … … 897 908 / ( u_comp(k) - gu + 1.0E-20 ) & 898 909 + diss_r(k) * & 899 ABS( u_comp(k) - 2.0 * hom(k,1,1,:) ) & 910 ABS( u_comp(k) - 2.0 * hom(k,1,1,:) ) & 900 911 / ( ABS( u_comp(k) - gu ) + 1.0E-20 ) ) & 901 912 * weight_substep(intermediate_timestep_count) * rmask(j,i,:) … … 998 1009 999 1010 1000 !----------------------------------------------------------------------------- -!1011 !-----------------------------------------------------------------------------! 1001 1012 ! Advection of v-component - Call for grid point i,j 1002 !----------------------------------------------------------------------------- -!1003 SUBROUTINE advec_v_ws_ij( i, j )1013 !-----------------------------------------------------------------------------! 1014 SUBROUTINE advec_v_ws_ij( i, j, i_omp, tn ) 1004 1015 1005 1016 USE arrays_3d … … 1012 1023 IMPLICIT NONE 1013 1024 1014 INTEGER :: i, j, k1025 INTEGER :: i, i_omp, j, k, tn 1015 1026 LOGICAL :: degraded_l, degraded_s 1016 1027 REAL :: gu, gv, flux_d, diss_d, u_comp, v_comp_l, w_comp … … 1036 1047 - ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3 1037 1048 diss_r(k) = - ABS( u_comp ) * ( & 1038 3.0 * ( v(k,j,i+1) - v(k,j,i) ) & 1049 3.0 * ( v(k,j,i+1) - v(k,j,i) ) & 1039 1050 - ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3 1040 1051 ENDDO … … 1056 1067 - ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3 1057 1068 diss_n(k) = - ABS( v_comp(k) - gv ) * ( & 1058 3.0 * ( v(k,j+1,i) - v(k,j,i) ) & 1069 3.0 * ( v(k,j+1,i) - v(k,j,i) ) & 1059 1070 - ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3 1060 1071 ENDDO … … 1077 1088 - ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3 1078 1089 diss_r(k) = - ABS( u_comp ) * ( & 1079 3.0 * ( w(k,j,i+1) - w(k,j,i) ) & 1090 3.0 * ( w(k,j,i+1) - w(k,j,i) ) & 1080 1091 - ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3 1081 1092 ENDDO … … 1084 1095 DO k = nzb_v_inner(j,i)+1, nzt 1085 1096 u_comp = u(k,j-1,i) + u(k,j,i) - gu 1086 flux_l_v(k,j ) = u_comp * ( v(k,j,i) + v(k,j,i-1) ) * 0.251087 diss_l_v(k,j ) = diss_2nd( v(k,j,i+2), v(k,j,i+1), v(k,j,i),&1097 flux_l_v(k,j,tn) = u_comp * ( v(k,j,i) + v(k,j,i-1) ) * 0.25 1098 diss_l_v(k,j,tn) = diss_2nd( v(k,j,i+2), v(k,j,i+1), v(k,j,i),& 1088 1099 v(k,j,i-1), v(k,j,i-1), u_comp, & 1089 1100 0.25, ddx ) … … 1094 1105 - ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3 1095 1106 diss_r(k) = - ABS( u_comp ) * ( & 1096 3.0 * ( v(k,j,i+1) - v(k,j,i) ) & 1107 3.0 * ( v(k,j,i+1) - v(k,j,i) ) & 1097 1108 - ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3 1098 1109 ENDDO … … 1102 1113 DO k = nzb_v_inner(j,i)+1, nzt 1103 1114 v_comp(k) = v(k,j,i) + v(k,j-1,i) - gv 1104 flux_s_v(k ) = v_comp(k) * ( v(k,j,i) + v(k,j-1,i) ) * 0.251105 diss_s_v(k ) = diss_2nd( v(k,j+2,i), v(k,j+1,i), v(k,j,i), &1115 flux_s_v(k,tn) = v_comp(k) * ( v(k,j,i) + v(k,j-1,i) ) * 0.25 1116 diss_s_v(k,tn) = diss_2nd( v(k,j+2,i), v(k,j+1,i), v(k,j,i), & 1106 1117 v(k,j-1,i), v(k,j-1,i), v_comp(k), & 1107 1118 0.25, ddy ) … … 1112 1123 - ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3 1113 1124 diss_n(k) = - ABS( v_comp(k) - gv ) * ( & 1114 3.0 * ( v(k,j+1,i) - v(k,j,i) ) & 1125 3.0 * ( v(k,j+1,i) - v(k,j,i) ) & 1115 1126 - ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3 1116 1127 ENDDO … … 1181 1192 ! 1182 1193 !-- Compute left- and southside fluxes for the respective boundary 1183 IF ( i == nxl .AND. ( .NOT. degraded_l ) ) THEN 1194 IF ( i == i_omp .AND. .NOT. degraded_l ) THEN 1195 1184 1196 DO k = nzb_v_inner(j,i)+1, nzt 1185 1197 u_comp = u(k,j-1,i) + u(k,j,i) - gu 1186 flux_l_v(k,j ) = u_comp * ( &1198 flux_l_v(k,j,tn) = u_comp * ( & 1187 1199 37.0 * ( v(k,j,i) + v(k,j,i-1) ) & 1188 1200 - 8.0 * ( v(k,j,i+1) + v(k,j,i-2) ) & 1189 1201 + ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5 1190 diss_l_v(k,j ) = - ABS( u_comp ) * ( &1202 diss_l_v(k,j,tn) = - ABS( u_comp ) * ( & 1191 1203 10.0 * ( v(k,j,i) - v(k,j,i-1) ) & 1192 1204 - 5.0 * ( v(k,j,i+1) - v(k,j,i-2) ) & … … 1196 1208 ENDIF 1197 1209 1198 IF ( j == nysv .AND. ( .NOT. degraded_s )) THEN1210 IF ( j == nysv .AND. .NOT. degraded_s ) THEN 1199 1211 1200 1212 DO k = nzb_v_inner(j,i)+1, nzt 1201 1213 v_comp_l = v(k,j,i) + v(k,j-1,i) - gv 1202 flux_s_v(k ) = v_comp_l * ( &1214 flux_s_v(k,tn) = v_comp_l * ( & 1203 1215 37.0 * ( v(k,j,i) + v(k,j-1,i) ) & 1204 1216 - 8.0 * ( v(k,j+1,i) + v(k,j-2,i) ) & 1205 1217 + ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5 1206 diss_s_v(k ) = - ABS( v_comp_l ) * ( &1218 diss_s_v(k,tn) = - ABS( v_comp_l ) * ( & 1207 1219 10.0 * ( v(k,j,i) - v(k,j-1,i) ) & 1208 1220 - 5.0 * ( v(k,j+1,i) - v(k,j-2,i) ) & … … 1216 1228 tend(k,j,i) = tend(k,j,i) - ( & 1217 1229 ( flux_r(k) + diss_r(k) & 1218 - flux_l_v(k,j ) - diss_l_v(k,j) ) * ddx &1230 - flux_l_v(k,j,tn) - diss_l_v(k,j,tn) ) * ddx & 1219 1231 + ( flux_n(k) + diss_n(k) & 1220 - flux_s_v(k ) - diss_s_v(k) ) * ddy )1232 - flux_s_v(k,tn) - diss_s_v(k,tn) ) * ddy ) 1221 1233 1222 flux_l_v(k,j ) = flux_r(k)1223 diss_l_v(k,j ) = diss_r(k)1224 flux_s_v(k ) = flux_n(k)1225 diss_s_v(k ) = diss_n(k)1234 flux_l_v(k,j,tn) = flux_r(k) 1235 diss_l_v(k,j,tn) = diss_r(k) 1236 flux_s_v(k,tn) = flux_n(k) 1237 diss_s_v(k,tn) = diss_n(k) 1226 1238 1227 1239 ! … … 1333 1345 ! Advection of w-component - Call for grid point i,j 1334 1346 !------------------------------------------------------------------------------! 1335 SUBROUTINE advec_w_ws_ij( i, j )1347 SUBROUTINE advec_w_ws_ij( i, j, i_omp, tn ) 1336 1348 1337 1349 USE arrays_3d … … 1344 1356 IMPLICIT NONE 1345 1357 1346 INTEGER :: i, j, k1358 INTEGER :: i, i_omp, j, k, tn 1347 1359 LOGICAL :: degraded_l, degraded_s 1348 1360 REAL :: gu, gv, flux_d, diss_d, u_comp, v_comp, w_comp … … 1355 1367 gu = 2.0 * u_gtrans 1356 1368 gv = 2.0 * v_gtrans 1357 1369 1370 1358 1371 IF ( boundary_flags(j,i) /= 0 ) THEN 1359 1372 ! … … 1425 1438 1426 1439 u_comp = u(k+1,j,i) + u(k,j,i) - gu 1427 flux_l_w(k,j ) = u_comp * ( w(k,j,i) + w(k,j,i-1) ) * 0.251428 diss_l_w(k,j ) = diss_2nd( w(k,j,i+2), w(k,j,i+1), w(k,j,i), &1440 flux_l_w(k,j,tn) = u_comp * ( w(k,j,i) + w(k,j,i-1) ) * 0.25 1441 diss_l_w(k,j,tn) = diss_2nd( w(k,j,i+2), w(k,j,i+1), w(k,j,i), & 1429 1442 w(k,j,i-1), w(k,j,i-1), u_comp, & 1430 1443 0.25, ddx ) … … 1455 1468 !-- Compute southside fluxes for the south boundary of PE domain 1456 1469 v_comp = v(k+1,j,i) + v(k,j,i) - gv 1457 flux_s_w(k ) = v_comp * ( w(k,j,i) + w(k,j-1,i) ) * 0.251458 diss_s_w(k ) = diss_2nd( w(k,j+2,i), w(k,j+1,i), w(k,j,i), &1470 flux_s_w(k,tn) = v_comp * ( w(k,j,i) + w(k,j-1,i) ) * 0.25 1471 diss_s_w(k,tn) = diss_2nd( w(k,j+2,i), w(k,j+1,i), w(k,j,i), & 1459 1472 w(k,j-1,i), w(k,j-1,i), v_comp, & 1460 1473 0.25, ddy ) … … 1526 1539 ! 1527 1540 !-- Compute left- and southside fluxes for the respective boundary 1528 IF ( j == nys .AND. ( .NOT. degraded_s )) THEN1541 IF ( j == nys .AND. .NOT. degraded_s ) THEN 1529 1542 1530 1543 DO k = nzb_w_inner(j,i)+1, nzt 1531 1544 v_comp = v(k+1,j,i) + v(k,j,i) - gv 1532 flux_s_w(k ) = v_comp * ( &1545 flux_s_w(k,tn) = v_comp * ( & 1533 1546 37.0 * ( w(k,j,i) + w(k,j-1,i) ) & 1534 1547 - 8.0 * ( w(k,j+1,i) +w(k,j-2,i) ) & 1535 1548 + ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5 1536 diss_s_w(k ) = - ABS( v_comp ) * ( &1549 diss_s_w(k,tn) = - ABS( v_comp ) * ( & 1537 1550 10.0 * ( w(k,j,i) - w(k,j-1,i) ) & 1538 1551 - 5.0 * ( w(k,j+1,i) - w(k,j-2,i) ) & … … 1542 1555 ENDIF 1543 1556 1544 IF ( i == nxl .AND. ( .NOT. degraded_l )) THEN1545 1557 IF ( i == i_omp .AND. .NOT. degraded_l ) THEN 1558 1546 1559 DO k = nzb_w_inner(j,i)+1, nzt 1547 1560 u_comp = u(k+1,j,i) + u(k,j,i) - gu 1548 flux_l_w(k,j ) = u_comp * ( &1561 flux_l_w(k,j,tn) = u_comp * ( & 1549 1562 37.0 * ( w(k,j,i) + w(k,j,i-1) ) & 1550 1563 - 8.0 * ( w(k,j,i+1) + w(k,j,i-2) ) & 1551 1564 + ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5 1552 diss_l_w(k,j ) = - ABS( u_comp ) * ( &1565 diss_l_w(k,j,tn) = - ABS( u_comp ) * ( & 1553 1566 10.0 * ( w(k,j,i) - w(k,j,i-1) ) & 1554 1567 - 5.0 * ( w(k,j,i+1) - w(k,j,i-2) ) & … … 1557 1570 1558 1571 ENDIF 1572 1559 1573 ! 1560 1574 !-- Now compute the tendency terms for the horizontal parts. … … 1562 1576 tend(k,j,i) = tend(k,j,i) - ( & 1563 1577 ( flux_r(k) + diss_r(k) & 1564 - flux_l_w(k,j ) - diss_l_w(k,j) ) * ddx &1578 - flux_l_w(k,j,tn) - diss_l_w(k,j,tn) ) * ddx & 1565 1579 + ( flux_n(k) + diss_n(k) & 1566 - flux_s_w(k ) - diss_s_w(k) ) * ddy )1567 1568 flux_l_w(k,j ) = flux_r(k)1569 diss_l_w(k,j ) = diss_r(k)1570 flux_s_w(k ) = flux_n(k)1571 diss_s_w(k ) = diss_n(k)1580 - flux_s_w(k,tn) - diss_s_w(k,tn) ) * ddy ) 1581 1582 flux_l_w(k,j,tn) = flux_r(k) 1583 diss_l_w(k,j,tn) = diss_r(k) 1584 flux_s_w(k,tn) = flux_n(k) 1585 diss_s_w(k,tn) = diss_n(k) 1572 1586 ENDDO 1573 1587 -
palm/trunk/SOURCE/modules.f90
r723 r736 5 5 ! Current revisions: 6 6 ! ----------------- 7 ! 7 ! Dimension of fluxes needed for WS-scheme increased. 8 8 ! 9 9 ! Former revisions: … … 246 246 km_damp_x, km_damp_y, lad, l_grid, pt_init, q_init, rdf, sa_init, & 247 247 ug, u_init, u_nzb_p1_for_vfc, vg, v_init, v_nzb_p1_for_vfc, w_subs, & 248 zu, zw, flux_s_u, flux_s_v, flux_s_w, diss_s_u, diss_s_v, diss_s_w, & 249 flux_s_pt, diss_s_pt, flux_s_e, diss_s_e, flux_s_q, diss_s_q, & 250 flux_s_sa, diss_s_sa 248 zu, zw 251 249 252 250 REAL, DIMENSION(:,:), ALLOCATABLE :: & 253 251 c_u, c_v, c_w, dzu_mg, dzw_mg, f1_mg, f2_mg, f3_mg, & 254 252 mean_inflow_profiles, pt_slope_ref, qs, qswst_remote, ts, us, z0, & 255 flux_ l_u, flux_l_v, flux_l_w, diss_l_u, diss_l_v, diss_l_w,&256 flux_ l_pt, diss_l_pt, flux_l_e, diss_l_e, flux_l_q, diss_l_q, &257 flux_l_sa, diss_l_sa, total_2d_o, total_2d_a253 flux_s_pt, diss_s_pt, flux_s_e, diss_s_e, flux_s_q, diss_s_q, & 254 flux_s_sa, diss_s_sa, flux_s_u, flux_s_v, flux_s_w, diss_s_u, & 255 diss_s_v, diss_s_w, total_2d_o, total_2d_a 258 256 259 257 REAL, DIMENSION(:,:), ALLOCATABLE, TARGET :: & … … 270 268 canopy_heat_flux, cdc, d, diss, lad_s, lad_u, lad_v, lad_w, lai, & 271 269 l_wall, p_loc, sec, sls, tend, u_m_l, u_m_n, u_m_r, u_m_s, v_m_l, & 272 v_m_n, v_m_r, v_m_s, w_m_l, w_m_n, w_m_r, w_m_s 270 v_m_n, v_m_r, v_m_s, w_m_l, w_m_n, w_m_r, w_m_s, flux_l_pt, & 271 diss_l_pt, flux_l_e, diss_l_e, flux_l_q, diss_l_q, flux_l_sa, & 272 diss_l_sa, flux_l_u, flux_l_v, flux_l_w, diss_l_u, diss_l_v, diss_l_w 273 273 274 274 REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET :: & -
palm/trunk/SOURCE/prognostic_equations.f90
r710 r736 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! Bugfix: determination of first thread index i for WS-scheme 7 7 ! 8 8 ! Former revisions: … … 147 147 148 148 CHARACTER (LEN=9) :: time_to_string 149 INTEGER :: i, j, k149 INTEGER :: i, i_omp_start, j, k, tn = 0 150 150 REAL :: sat, sbt 151 151 … … 172 172 ! 173 173 !-- u-tendency terms with no communication 174 i_omp_start = nxlu 174 175 DO i = nxlu, nxr 175 176 DO j = nys, nyn … … 179 180 tend(:,j,i) = 0.0 180 181 IF ( ws_scheme_mom ) THEN 181 CALL advec_u_ws( i, j )182 CALL advec_u_ws( i, j, i_omp_start, tn ) 182 183 ELSE 183 184 CALL advec_u_pw( i, j ) … … 257 258 ! 258 259 !-- v-tendency terms with no communication 260 i_omp_start = nxl 259 261 DO i = nxl, nxr 260 262 DO j = nysv, nyn … … 264 266 tend(:,j,i) = 0.0 265 267 IF ( ws_scheme_mom ) THEN 266 CALL advec_v_ws( i, j )268 CALL advec_v_ws( i, j, i_omp_start, tn ) 267 269 ELSE 268 270 CALL advec_v_pw( i, j ) … … 348 350 tend(:,j,i) = 0.0 349 351 IF ( ws_scheme_mom ) THEN 350 CALL advec_w_ws( i, j )352 CALL advec_w_ws( i, j, i_omp_start, tn ) 351 353 ELSE 352 354 CALL advec_w_pw( i, j ) … … 453 455 IF ( ws_scheme_sca ) THEN 454 456 CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, & 455 diss_s_pt, flux_l_pt, diss_l_pt )457 diss_s_pt, flux_l_pt, diss_l_pt, i_omp_start, tn ) 456 458 ELSE 457 459 CALL advec_s_pw( i, j, pt ) … … 574 576 IF ( ws_scheme_sca ) THEN 575 577 CALL advec_s_ws( i, j, sa, 'sa', flux_s_sa, & 576 diss_s_sa, flux_l_sa, diss_l_sa )578 diss_s_sa, flux_l_sa, diss_l_sa, i_omp_start, tn ) 577 579 ELSE 578 580 CALL advec_s_pw( i, j, sa ) … … 673 675 IF ( ws_scheme_sca ) THEN 674 676 CALL advec_s_ws( i, j, q, 'q', flux_s_q, & 675 diss_s_q, flux_l_q, diss_l_q )677 diss_s_q, flux_l_q, diss_l_q, i_omp_start, tn ) 676 678 ELSE 677 679 CALL advec_s_pw( i, j, q ) … … 813 815 IF ( ws_scheme_sca ) THEN 814 816 CALL advec_s_ws( i, j, e, 'e', flux_s_e, & 815 diss_s_e, flux_l_e, diss_l_e )817 diss_s_e, flux_l_e, diss_l_e, i_omp_start, tn ) 816 818 ELSE 817 819 CALL advec_s_pw( i, j, e ) … … 916 918 917 919 CHARACTER (LEN=9) :: time_to_string 918 INTEGER :: i, j, k 920 INTEGER :: i, i_omp_start, j, k, omp_get_thread_num, tn = 0 921 LOGICAL :: loop_start 919 922 920 923 … … 934 937 intermediate_timestep_count == 1 ) CALL ws_statistics 935 938 936 937 939 ! 938 940 !-- Loop over all prognostic equations 939 !$OMP PARALLEL private (i,j,k) 941 !$OMP PARALLEL private (i,i_omp_start,j,k,loop_start,tn) 942 943 !$ tn = omp_get_thread_num() 944 loop_start = .TRUE. 940 945 !$OMP DO 941 946 DO i = nxl, nxr 947 948 ! 949 !-- Store the first loop index. It differs for each thread and is required 950 !-- later in advec_ws 951 IF ( loop_start ) THEN 952 loop_start = .FALSE. 953 i_omp_start = i 954 ENDIF 955 942 956 DO j = nys, nyn 943 957 ! … … 948 962 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN 949 963 IF ( ws_scheme_mom ) THEN 950 ! CALL local_diss( i, j, u) ! dissipation control 951 CALL advec_u_ws( i, j ) 964 IF ( outflow_l .AND. i_omp_start == nxl ) THEN 965 ! CALL local_diss( i, j, u) ! dissipation control 966 CALL advec_u_ws( i, j, i_omp_start + 1, tn ) 967 ELSE 968 CALL advec_u_ws( i, j, i_omp_start, tn ) 969 ENDIF 952 970 ELSE 953 971 CALL advec_u_pw( i, j ) 954 972 ENDIF 955 ELSE973 ELSE 956 974 CALL advec_u_up( i, j ) 957 975 ENDIF … … 1017 1035 IF ( ws_scheme_mom ) THEN 1018 1036 ! CALL local_diss( i, j, v) 1019 CALL advec_v_ws( i, j )1037 CALL advec_v_ws( i, j, i_omp_start, tn ) 1020 1038 ELSE 1021 1039 CALL advec_v_pw( i, j ) … … 1081 1099 IF ( ws_scheme_mom ) THEN 1082 1100 ! CALL local_diss( i, j, w) 1083 CALL advec_w_ws( i, j )1101 CALL advec_w_ws( i, j, i_omp_start, tn ) 1084 1102 ELSE 1085 1103 CALL advec_w_pw( i, j ) … … 1145 1163 ! CALL local_diss( i, j, pt ) 1146 1164 CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, & 1147 diss_s_pt, flux_l_pt, diss_l_pt )1165 diss_s_pt, flux_l_pt, diss_l_pt, i_omp_start, tn ) 1148 1166 ELSE 1149 1167 CALL advec_s_pw( i, j, pt ) … … 1227 1245 ! CALL local_diss( i, j, sa ) 1228 1246 CALL advec_s_ws( i, j, sa, 'sa', flux_s_sa, & 1229 diss_s_sa, flux_l_sa, diss_l_sa )1247 diss_s_sa, flux_l_sa, diss_l_sa, i_omp_start, tn ) 1230 1248 ELSE 1231 1249 CALL advec_s_pw( i, j, sa ) … … 1285 1303 ! CALL local_diss( i, j, q ) 1286 1304 CALL advec_s_ws( i, j, q, 'q', flux_s_q, & 1287 diss_s_q, flux_l_q, diss_l_q )1305 diss_s_q, flux_l_q, diss_l_q, i_omp_start, tn ) 1288 1306 ELSE 1289 1307 CALL advec_s_pw( i, j, q ) … … 1361 1379 ! CALL local_diss( i, j, e ) 1362 1380 CALL advec_s_ws( i, j, e, 'e', flux_s_e, & 1363 diss_s_e, flux_l_e, diss_l_e )1381 diss_s_e, flux_l_e, diss_l_e , i_omp_start, tn ) 1364 1382 ELSE 1365 1383 CALL advec_s_pw( i, j, e )
Note: See TracChangeset
for help on using the changeset viewer.
