Changeset 114 for palm/trunk/SOURCE
- Timestamp:
- Oct 10, 2007 12:03:15 AM (17 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 9 edited
-
CURRENT_MODIFICATIONS (modified) (2 diffs)
-
advec_particles.f90 (modified) (3 diffs)
-
check_parameters.f90 (modified) (2 diffs)
-
init_grid.f90 (modified) (9 diffs)
-
init_pegrid.f90 (modified) (2 diffs)
-
modules.f90 (modified) (2 diffs)
-
palm.f90 (modified) (1 diff)
-
poismg.f90 (modified) (20 diffs)
-
user_interface.f90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/CURRENT_MODIFICATIONS
r110 r114 1 1 New: 2 2 --- 3 4 Pressure boundary conditions for vertical walls added to the multigrid solver. 5 They are applied using new wall flag arrays (wall_flags_..) which are defined 6 for each grid level. New argument gls added to routine user_init_grid 7 (user_interface). 8 9 init_grid, init_pegrid, modules, user_interface 3 10 4 11 … … 9 16 Errors: 10 17 ------ 18 19 Bugfix: pleft/pright changed to pnorth/psouth in sendrecv of particle tail 20 numbers along y (advec_particles) 21 22 Bugfix: model_string needed a default value (combine_plot_fields) 23 24 advec_particles, combine_plot_fields -
palm/trunk/SOURCE/advec_particles.f90
r110 r114 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! Bugfix: pleft/pright changed to pnorth/psouth in sendrecv of particle tail 7 ! numbers along y 6 8 ! TEST: PRINT statements on unit 9 (commented out) 7 9 ! … … 2773 2775 IF ( use_particle_tails ) THEN 2774 2776 2775 CALL MPI_SENDRECV( trspt_count, 1, MPI_INTEGER, p left,0, &2776 trnpt_count_recv, 1, MPI_INTEGER, p right, 0, &2777 CALL MPI_SENDRECV( trspt_count, 1, MPI_INTEGER, psouth, 0, & 2778 trnpt_count_recv, 1, MPI_INTEGER, pnorth, 0, & 2777 2779 comm2d, status, ierr ) 2778 2780 … … 2850 2852 IF ( use_particle_tails ) THEN 2851 2853 2852 CALL MPI_SENDRECV( trnpt_count, 1, MPI_INTEGER, p right, 0, &2853 trspt_count_recv, 1, MPI_INTEGER, p left,0, &2854 CALL MPI_SENDRECV( trnpt_count, 1, MPI_INTEGER, pnorth, 0, & 2855 trspt_count_recv, 1, MPI_INTEGER, psouth, 0, & 2854 2856 comm2d, status, ierr ) 2855 2857 -
palm/trunk/SOURCE/check_parameters.f90
r110 r114 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! Multigrid solver allows topography 7 7 ! 8 8 ! Former revisions: … … 286 286 WRITE( action, '(A,A)' ) 'timestep_scheme = ', timestep_scheme 287 287 ENDIF 288 IF ( psolver == ' multigrid' .OR. psolver == 'sor' ) THEN288 IF ( psolver == 'sor' ) THEN 289 289 WRITE( action, '(A,A)' ) 'psolver = ', psolver 290 290 ENDIF -
palm/trunk/SOURCE/init_grid.f90
r98 r114 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! Calculation of wall flag arrays 7 7 ! 8 8 ! Former revisions: … … 44 44 IMPLICIT NONE 45 45 46 INTEGER :: bh, blx, bly, bxl, bxr, byn, bys, i, i_center, j, j_center, k, &47 l, nzb_si, nzt_l, vi46 INTEGER :: bh, blx, bly, bxl, bxr, byn, bys, gls, i, inc, i_center, j, & 47 j_center, k, l, nxl_l, nxr_l, nyn_l, nys_l, nzb_si, nzt_l, vi 48 48 49 49 INTEGER, DIMENSION(:), ALLOCATABLE :: vertical_influence … … 217 217 ! 218 218 !-- Allocate outer and inner index arrays for topography and set 219 !-- defaults 220 ALLOCATE( corner_nl(nys:nyn,nxl:nxr), corner_nr(nys:nyn,nxl:nxr), & 221 corner_sl(nys:nyn,nxl:nxr), corner_sr(nys:nyn,nxl:nxr), & 222 nzb_local(-1:ny+1,-1:nx+1), nzb_tmp(-1:ny+1,-1:nx+1), & 223 wall_l(nys:nyn,nxl:nxr), wall_n(nys:nyn,nxl:nxr), & 219 !-- defaults. 220 !-- nzb_local has to contain additional layers of ghost points for calculating 221 !-- the flag arrays needed for the multigrid method 222 gls = 2**( maximum_grid_level ) 223 ALLOCATE( corner_nl(nys:nyn,nxl:nxr), corner_nr(nys:nyn,nxl:nxr), & 224 corner_sl(nys:nyn,nxl:nxr), corner_sr(nys:nyn,nxl:nxr), & 225 nzb_local(-gls:ny+gls,-gls:nx+gls), nzb_tmp(-1:ny+1,-1:nx+1), & 226 wall_l(nys:nyn,nxl:nxr), wall_n(nys:nyn,nxl:nxr), & 224 227 wall_r(nys:nyn,nxl:nxr), wall_s(nys:nyn,nxl:nxr) ) 225 228 ALLOCATE( fwxm(nys-1:nyn+1,nxl-1:nxr+1), fwxp(nys-1:nyn+1,nxl-1:nxr+1), & … … 388 391 ENDDO 389 392 ENDDO 390 nzb_local(-1,0:nx) = nzb_local(ny,0:nx) 391 nzb_local(ny+1,0:nx) = nzb_local(0,0:nx) 392 nzb_local(:,-1) = nzb_local(:,nx) 393 nzb_local(:,nx+1) = nzb_local(:,0) 393 ! 394 !-- Add cyclic boundaries (additional layers are for calculating flag 395 !-- arrays needed for the multigrid sover) 396 nzb_local(-gls:-1,0:nx) = nzb_local(ny-gls+1:ny,0:nx) 397 nzb_local(ny+1:ny+gls,0:nx) = nzb_local(0:gls-1,0:nx) 398 nzb_local(:,-gls:-1) = nzb_local(:,nx-gls+1:nx) 399 nzb_local(:,nx+1:nx+gls) = nzb_local(:,0:gls-1) 394 400 395 401 GOTO 12 … … 413 419 !-- case in the user interface. There, the subroutine user_init_grid 414 420 !-- checks which of these two conditions applies. 415 CALL user_init_grid( nzb_local )421 CALL user_init_grid( gls, nzb_local ) 416 422 417 423 END SELECT 424 425 ! 426 !-- Test output of nzb_local -1:ny+1,-1:nx+1 427 WRITE (9,*) '*** nzb_local ***' 428 DO j = ny+1, -1, -1 429 WRITE (9,'(194(1X,I2))') ( nzb_local(j,i), i = -1, nx+1 ) 430 ENDDO 418 431 419 432 ! … … 478 491 !-- nzb_s_outer: 479 492 !-- extend nzb_local east-/westwards first, then north-/southwards 480 nzb_tmp = nzb_local 493 nzb_tmp = nzb_local(-1:ny+1,-1:nx+1) 481 494 DO j = -1, ny + 1 482 495 DO i = 0, nx … … 510 523 !-- nzb_u_inner: 511 524 !-- extend nzb_local rightwards only 512 nzb_tmp = nzb_local 525 nzb_tmp = nzb_local(-1:ny+1,-1:nx+1) 513 526 DO j = -1, ny + 1 514 527 DO i = 0, nx + 1 … … 542 555 !-- nzb_v_inner: 543 556 !-- extend nzb_local northwards only 544 nzb_tmp = nzb_local 557 nzb_tmp = nzb_local(-1:ny+1,-1:nx+1) 545 558 DO i = -1, nx + 1 546 559 DO j = 0, ny + 1 … … 728 741 729 742 ! 743 !-- Calculate wall flag arrays for the multigrid method 744 IF ( psolver == 'multigrid' ) THEN 745 ! 746 !-- Gridpoint increment of the current level 747 inc = 1 748 749 DO l = maximum_grid_level, 1 , -1 750 751 nxl_l = nxl_mg(l) 752 nxr_l = nxr_mg(l) 753 nys_l = nys_mg(l) 754 nyn_l = nyn_mg(l) 755 nzt_l = nzt_mg(l) 756 757 ! 758 !-- Assign the flag level to be calculated 759 SELECT CASE ( l ) 760 CASE ( 1 ) 761 flags => wall_flags_1 762 CASE ( 2 ) 763 flags => wall_flags_2 764 CASE ( 3 ) 765 flags => wall_flags_3 766 CASE ( 4 ) 767 flags => wall_flags_4 768 CASE ( 5 ) 769 flags => wall_flags_5 770 CASE ( 6 ) 771 flags => wall_flags_6 772 CASE ( 7 ) 773 flags => wall_flags_7 774 CASE ( 8 ) 775 flags => wall_flags_8 776 CASE ( 9 ) 777 flags => wall_flags_9 778 CASE ( 10 ) 779 flags => wall_flags_10 780 END SELECT 781 782 ! 783 !-- Depending on the grid level, set the respective bits in case of 784 !-- neighbouring walls 785 !-- Bit 0: wall to the bottom 786 !-- Bit 1: wall to the top (not realized in remaining PALM code so far) 787 !-- Bit 2: wall to the south 788 !-- Bit 3: wall to the north 789 !-- Bit 4: wall to the left 790 !-- Bit 5: wall to the right 791 !-- Bit 6: inside / outside building (1/0) 792 793 flags = 0 794 795 DO i = nxl_l-1, nxr_l+1 796 DO j = nys_l-1, nyn_l+1 797 DO k = nzb, nzt_l+1 798 799 ! 800 !-- Inside/outside building (inside building does not need 801 !-- further tests for walls) 802 IF ( k*inc <= nzb_local(j*inc,i*inc) ) THEN 803 804 flags(k,j,i) = IBSET( flags(k,j,i), 6 ) 805 806 ELSE 807 ! 808 !-- Bottom wall 809 IF ( (k-1)*inc <= nzb_local(j*inc,i*inc) ) THEN 810 flags(k,j,i) = IBSET( flags(k,j,i), 0 ) 811 ENDIF 812 ! 813 !-- South wall 814 IF ( k*inc <= nzb_local((j-1)*inc,i*inc) ) THEN 815 flags(k,j,i) = IBSET( flags(k,j,i), 2 ) 816 ENDIF 817 ! 818 !-- North wall 819 IF ( k*inc <= nzb_local((j+1)*inc,i*inc) ) THEN 820 flags(k,j,i) = IBSET( flags(k,j,i), 3 ) 821 ENDIF 822 ! 823 !-- Left wall 824 IF ( k*inc <= nzb_local(j*inc,(i-1)*inc) ) THEN 825 flags(k,j,i) = IBSET( flags(k,j,i), 4 ) 826 ENDIF 827 ! 828 !-- Right wall 829 IF ( k*inc <= nzb_local(j*inc,(i+1)*inc) ) THEN 830 flags(k,j,i) = IBSET( flags(k,j,i), 5 ) 831 ENDIF 832 833 ENDIF 834 835 ENDDO 836 ENDDO 837 ENDDO 838 839 ! 840 !-- Test output of flag arrays 841 i = nxl_l 842 WRITE (9,*) ' ' 843 WRITE (9,*) '*** mg level ', l, ' ***', mg_switch_to_pe0_level 844 WRITE (9,*) ' inc=', inc, ' i =', nxl_l 845 WRITE (9,*) ' nxl_l',nxl_l,' nxr_l=',nxr_l,' nys_l=',nys_l,' nyn_l=',nyn_l 846 DO k = nzt_l+1, nzb, -1 847 WRITE (9,'(194(1X,I2))') ( flags(k,j,i), j = nys_l-1, nyn_l+1 ) 848 ENDDO 849 850 inc = inc * 2 851 852 ENDDO 853 854 ENDIF 855 856 ! 730 857 !-- In case of topography: limit near-wall mixing length l_wall further: 731 858 !-- Go through all points of the subdomain one by one and look for the closest -
palm/trunk/SOURCE/init_pegrid.f90
r110 r114 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! Allocation of wall flag arrays for multigrid solver 6 7 ! TEST OUTPUT (TO BE REMOVED) logging mpi2 ierr values 7 8 ! … … 915 916 ENDIF 916 917 918 ! 919 !-- Allocate wall flag arrays used in the multigrid solver 920 IF ( psolver == 'multigrid' ) THEN 921 922 DO i = maximum_grid_level, 1, -1 923 924 SELECT CASE ( i ) 925 926 CASE ( 1 ) 927 ALLOCATE( wall_flags_1(nzb:nzt_mg(i)+1, & 928 nys_mg(i)-1:nyn_mg(i)+1, & 929 nxl_mg(i)-1:nxr_mg(i)+1) ) 930 931 CASE ( 2 ) 932 ALLOCATE( wall_flags_2(nzb:nzt_mg(i)+1, & 933 nys_mg(i)-1:nyn_mg(i)+1, & 934 nxl_mg(i)-1:nxr_mg(i)+1) ) 935 936 CASE ( 3 ) 937 ALLOCATE( wall_flags_3(nzb:nzt_mg(i)+1, & 938 nys_mg(i)-1:nyn_mg(i)+1, & 939 nxl_mg(i)-1:nxr_mg(i)+1) ) 940 941 CASE ( 4 ) 942 ALLOCATE( wall_flags_4(nzb:nzt_mg(i)+1, & 943 nys_mg(i)-1:nyn_mg(i)+1, & 944 nxl_mg(i)-1:nxr_mg(i)+1) ) 945 946 CASE ( 5 ) 947 ALLOCATE( wall_flags_5(nzb:nzt_mg(i)+1, & 948 nys_mg(i)-1:nyn_mg(i)+1, & 949 nxl_mg(i)-1:nxr_mg(i)+1) ) 950 951 CASE ( 6 ) 952 ALLOCATE( wall_flags_6(nzb:nzt_mg(i)+1, & 953 nys_mg(i)-1:nyn_mg(i)+1, & 954 nxl_mg(i)-1:nxr_mg(i)+1) ) 955 956 CASE ( 7 ) 957 ALLOCATE( wall_flags_7(nzb:nzt_mg(i)+1, & 958 nys_mg(i)-1:nyn_mg(i)+1, & 959 nxl_mg(i)-1:nxr_mg(i)+1) ) 960 961 CASE ( 8 ) 962 ALLOCATE( wall_flags_8(nzb:nzt_mg(i)+1, & 963 nys_mg(i)-1:nyn_mg(i)+1, & 964 nxl_mg(i)-1:nxr_mg(i)+1) ) 965 966 CASE ( 9 ) 967 ALLOCATE( wall_flags_9(nzb:nzt_mg(i)+1, & 968 nys_mg(i)-1:nyn_mg(i)+1, & 969 nxl_mg(i)-1:nxr_mg(i)+1) ) 970 971 CASE ( 10 ) 972 ALLOCATE( wall_flags_10(nzb:nzt_mg(i)+1, & 973 nys_mg(i)-1:nyn_mg(i)+1, & 974 nxl_mg(i)-1:nxr_mg(i)+1) ) 975 976 CASE DEFAULT 977 IF ( myid == 0 ) PRINT*, '+++ init_pegrid: more than 10 ', & 978 ' multigrid levels' 979 CALL local_stop 980 981 END SELECT 982 983 ENDDO 984 985 ENDIF 986 917 987 END SUBROUTINE init_pegrid -
palm/trunk/SOURCE/modules.f90
r110 r114 5 5 ! Actual revisions: 6 6 ! ----------------- 7 ! 7 ! +flags, wall_flags_1..10 8 8 ! 9 9 ! Former revisions: … … 589 589 nzb_v_outer, nzb_w_inner, nzb_w_outer, nzb_2d 590 590 591 INTEGER, DIMENSION(:,:,:), POINTER :: flags 592 593 INTEGER, DIMENSION(:,:,:), ALLOCATABLE, TARGET :: wall_flags_1, & 594 wall_flags_2, wall_flags_3, wall_flags_4, wall_flags_5, & 595 wall_flags_6, wall_flags_7, wall_flags_8, wall_flags_9, & 596 wall_flags_10 597 591 598 SAVE 592 599 -
palm/trunk/SOURCE/palm.f90
r110 r114 66 66 INTEGER :: i, run_description_header_i(80) 67 67 68 version = 'PALM 3.4 '68 version = 'PALM 3.4a' 69 69 70 70 #if defined( __parallel ) -
palm/trunk/SOURCE/poismg.f90
r77 r114 8 8 ! Actual revisions: 9 9 ! ----------------- 10 ! 10 ! Boundary conditions at walls are implicitly set using flag arrays. Only 11 ! Neumann BC is allowed. Upper walls are still not realized. 12 ! Bottom and top BCs for array f_mg in restrict removed because boundary 13 ! values are not needed (right hand side of SOR iteration) 11 14 ! 12 15 ! Former revisions: … … 150 153 l = grid_level 151 154 155 ! 156 !-- Choose flag array of this level 157 SELECT CASE ( l ) 158 CASE ( 1 ) 159 flags => wall_flags_1 160 CASE ( 2 ) 161 flags => wall_flags_2 162 CASE ( 3 ) 163 flags => wall_flags_3 164 CASE ( 4 ) 165 flags => wall_flags_4 166 CASE ( 5 ) 167 flags => wall_flags_5 168 CASE ( 6 ) 169 flags => wall_flags_6 170 CASE ( 7 ) 171 flags => wall_flags_7 172 CASE ( 8 ) 173 flags => wall_flags_8 174 CASE ( 9 ) 175 flags => wall_flags_9 176 CASE ( 10 ) 177 flags => wall_flags_10 178 END SELECT 179 152 180 !$OMP PARALLEL PRIVATE (i,j,k) 153 181 !$OMP DO … … 155 183 DO j = nys_mg(l), nyn_mg(l) 156 184 DO k = nzb+1, nzt_mg(l) 157 r(k,j,i) = f_mg(k,j,i) & 158 - ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 159 - ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 160 - f2_mg(k,l) * p_mg(k+1,j,i) & 161 - f3_mg(k,l) * p_mg(k-1,j,i) & 185 r(k,j,i) = f_mg(k,j,i) & 186 - ddx2_mg(l) * & 187 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 188 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 189 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 190 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 191 - ddy2_mg(l) * & 192 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 193 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & 194 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * & 195 ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) & 196 - f2_mg(k,l) * p_mg(k+1,j,i) & 197 - f3_mg(k,l) * & 198 ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * & 199 ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) & 162 200 + f1_mg(k,l) * p_mg(k,j,i) 201 ! 202 !-- Residual within topography should be zero 203 r(k,j,i) = r(k,j,i) * ( 1.0 - IBITS( flags(k,j,i), 6, 1 ) ) 163 204 ENDDO 164 205 ENDDO … … 181 222 182 223 ! 183 !-- Bottom and top boundary conditions 184 IF ( ibc_p_b == 1 ) THEN 185 r(nzb,:,: ) = r(nzb+1,:,:) 186 ELSE 187 r(nzb,:,: ) = 0.0 188 ENDIF 189 224 !-- Top boundary condition 225 !-- A Neumann boundary condition for r is implicitly set in routine restrict 190 226 IF ( ibc_p_t == 1 ) THEN 191 227 r(nzt_mg(l)+1,:,: ) = r(nzt_mg(l),:,:) … … 217 253 INTEGER :: i, ic, j, jc, k, kc, l 218 254 255 REAL :: rkjim, rkjip, rkjmi, rkjmim, rkjmip, rkjpi, rkjpim, rkjpip, & 256 rkmji, rkmjim, rkmjip, rkmjmi, rkmjmim, rkmjmip, rkmjpi, rkmjpim, & 257 rkmjpip 258 219 259 REAL, DIMENSION(nzb:nzt_mg(grid_level)+1, & 220 260 nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & … … 229 269 l = grid_level 230 270 271 ! 272 !-- Choose flag array of the upper level 273 SELECT CASE ( l ) 274 CASE ( 1 ) 275 flags => wall_flags_1 276 CASE ( 2 ) 277 flags => wall_flags_2 278 CASE ( 3 ) 279 flags => wall_flags_3 280 CASE ( 4 ) 281 flags => wall_flags_4 282 CASE ( 5 ) 283 flags => wall_flags_5 284 CASE ( 6 ) 285 flags => wall_flags_6 286 CASE ( 7 ) 287 flags => wall_flags_7 288 CASE ( 8 ) 289 flags => wall_flags_8 290 CASE ( 9 ) 291 flags => wall_flags_9 292 CASE ( 10 ) 293 flags => wall_flags_10 294 END SELECT 295 231 296 !$OMP PARALLEL PRIVATE (i,j,k,ic,jc,kc) 232 297 !$OMP DO … … 237 302 DO kc = nzb+1, nzt_mg(l) 238 303 k = 2*kc-1 304 ! 305 !-- Use implicit Neumann BCs if the respective gridpoint is inside 306 !-- the building 307 rkjim = r(k,j,i-1) + IBITS( flags(k,j,i-1), 6, 1 ) * & 308 ( r(k,j,i) - r(k,j,i-1) ) 309 rkjip = r(k,j,i+1) + IBITS( flags(k,j,i+1), 6, 1 ) * & 310 ( r(k,j,i) - r(k,j,i+1) ) 311 rkjpi = r(k,j+1,i) + IBITS( flags(k,j+1,i), 6, 1 ) * & 312 ( r(k,j,i) - r(k,j+1,i) ) 313 rkjmi = r(k,j-1,i) + IBITS( flags(k,j-1,i), 6, 1 ) * & 314 ( r(k,j,i) - r(k,j-1,i) ) 315 rkjmim = r(k,j-1,i-1) + IBITS( flags(k,j-1,i-1), 6, 1 ) * & 316 ( r(k,j,i) - r(k,j-1,i-1) ) 317 rkjpim = r(k,j+1,i-1) + IBITS( flags(k,j+1,i-1), 6, 1 ) * & 318 ( r(k,j,i) - r(k,j+1,i-1) ) 319 rkjmip = r(k,j-1,i+1) + IBITS( flags(k,j-1,i+1), 6, 1 ) * & 320 ( r(k,j,i) - r(k,j-1,i+1) ) 321 rkjpip = r(k,j+1,i+1) + IBITS( flags(k,j+1,i+1), 6, 1 ) * & 322 ( r(k,j,i) - r(k,j+1,i+1) ) 323 rkmji = r(k-1,j,i) + IBITS( flags(k-1,j,i), 6, 1 ) * & 324 ( r(k,j,i) - r(k-1,j,i) ) 325 rkmjim = r(k-1,j,i-1) + IBITS( flags(k-1,j,i-1), 6, 1 ) * & 326 ( r(k,j,i) - r(k-1,j,i-1) ) 327 rkmjip = r(k-1,j,i+1) + IBITS( flags(k-1,j,i+1), 6, 1 ) * & 328 ( r(k,j,i) - r(k-1,j,i+1) ) 329 rkmjpi = r(k-1,j+1,i) + IBITS( flags(k-1,j+1,i), 6, 1 ) * & 330 ( r(k,j,i) - r(k-1,j+1,i) ) 331 rkmjmi = r(k-1,j-1,i) + IBITS( flags(k-1,j-1,i), 6, 1 ) * & 332 ( r(k,j,i) - r(k-1,j-1,i) ) 333 rkmjmim = r(k-1,j-1,i-1) + IBITS( flags(k-1,j-1,i-1), 6, 1 ) * & 334 ( r(k,j,i) - r(k-1,j-1,i-1) ) 335 rkmjpim = r(k-1,j+1,i-1) + IBITS( flags(k-1,j+1,i-1), 6, 1 ) * & 336 ( r(k,j,i) - r(k-1,j+1,i-1) ) 337 rkmjmip = r(k-1,j-1,i+1) + IBITS( flags(k-1,j-1,i+1), 6, 1 ) * & 338 ( r(k,j,i) - r(k-1,j-1,i+1) ) 339 rkmjpip = r(k-1,j+1,i+1) + IBITS( flags(k-1,j+1,i+1), 6, 1 ) * & 340 ( r(k,j,i) - r(k-1,j+1,i+1) ) 341 239 342 f_mg(kc,jc,ic) = 1.0 / 64.0 * ( & 240 343 8.0 * r(k,j,i) & 241 + 4.0 * ( r (k,j,i-1) + r(k,j,i+1) +&242 r (k,j+1,i) + r(k,j-1,i) )&243 + 2.0 * ( r (k,j-1,i-1) + r(k,j+1,i-1) +&244 r (k,j-1,i+1) + r(k,j+1,i+1) )&245 + 4.0 * r (k-1,j,i)&246 + 2.0 * ( r (k-1,j,i-1) + r(k-1,j,i+1) +&247 r (k-1,j+1,i) + r(k-1,j-1,i) )&248 + ( r (k-1,j-1,i-1) + r(k-1,j+1,i-1) +&249 r (k-1,j-1,i+1) + r(k-1,j+1,i+1) )&344 + 4.0 * ( rkjim + rkjip + & 345 rkjpi + rkjmi ) & 346 + 2.0 * ( rkjmim + rkjpim + & 347 rkjmip + rkjpip ) & 348 + 4.0 * rkmji & 349 + 2.0 * ( rkmjim + rkmjim + & 350 rkmjpi + rkmjmi ) & 351 + ( rkmjmim + rkmjpim + & 352 rkmjmip + rkmjpip ) & 250 353 + 4.0 * r(k+1,j,i) & 251 354 + 2.0 * ( r(k+1,j,i-1) + r(k+1,j,i+1) + & … … 254 357 r(k+1,j-1,i+1) + r(k+1,j+1,i+1) ) & 255 358 ) 359 360 ! f_mg(kc,jc,ic) = 1.0 / 64.0 * ( & 361 ! 8.0 * r(k,j,i) & 362 ! + 4.0 * ( r(k,j,i-1) + r(k,j,i+1) + & 363 ! r(k,j+1,i) + r(k,j-1,i) ) & 364 ! + 2.0 * ( r(k,j-1,i-1) + r(k,j+1,i-1) + & 365 ! r(k,j-1,i+1) + r(k,j+1,i+1) ) & 366 ! + 4.0 * r(k-1,j,i) & 367 ! + 2.0 * ( r(k-1,j,i-1) + r(k-1,j,i+1) + & 368 ! r(k-1,j+1,i) + r(k-1,j-1,i) ) & 369 ! + ( r(k-1,j-1,i-1) + r(k-1,j+1,i-1) + & 370 ! r(k-1,j-1,i+1) + r(k-1,j+1,i+1) ) & 371 ! + 4.0 * r(k+1,j,i) & 372 ! + 2.0 * ( r(k+1,j,i-1) + r(k+1,j,i+1) + & 373 ! r(k+1,j+1,i) + r(k+1,j-1,i) ) & 374 ! + ( r(k+1,j-1,i-1) + r(k+1,j+1,i-1) + & 375 ! r(k+1,j-1,i+1) + r(k+1,j+1,i+1) ) & 376 ! ) 256 377 ENDDO 257 378 ENDDO … … 275 396 ! 276 397 !-- Bottom and top boundary conditions 277 IF ( ibc_p_b == 1 ) THEN278 f_mg(nzb,:,: ) = f_mg(nzb+1,:,:)279 ELSE280 f_mg(nzb,:,: ) = 0.0281 ENDIF282 283 IF ( ibc_p_t == 1 ) THEN284 f_mg(nzt_mg(l)+1,:,: ) = f_mg(nzt_mg(l),:,:)285 ELSE286 f_mg(nzt_mg(l)+1,:,: ) = 0.0287 ENDIF398 ! IF ( ibc_p_b == 1 ) THEN 399 ! f_mg(nzb,:,: ) = f_mg(nzb+1,:,:) 400 ! ELSE 401 ! f_mg(nzb,:,: ) = 0.0 402 ! ENDIF 403 ! 404 ! IF ( ibc_p_t == 1 ) THEN 405 ! f_mg(nzt_mg(l)+1,:,: ) = f_mg(nzt_mg(l),:,:) 406 ! ELSE 407 ! f_mg(nzt_mg(l)+1,:,: ) = 0.0 408 ! ENDIF 288 409 289 410 … … 412 533 LOGICAL :: unroll 413 534 535 REAL :: wall_left, wall_north, wall_right, wall_south, wall_total, wall_top 536 414 537 REAL, DIMENSION(nzb:nzt_mg(grid_level)+1, & 415 538 nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & … … 418 541 419 542 l = grid_level 543 544 ! 545 !-- Choose flag array of this level 546 SELECT CASE ( l ) 547 CASE ( 1 ) 548 flags => wall_flags_1 549 CASE ( 2 ) 550 flags => wall_flags_2 551 CASE ( 3 ) 552 flags => wall_flags_3 553 CASE ( 4 ) 554 flags => wall_flags_4 555 CASE ( 5 ) 556 flags => wall_flags_5 557 CASE ( 6 ) 558 flags => wall_flags_6 559 CASE ( 7 ) 560 flags => wall_flags_7 561 CASE ( 8 ) 562 flags => wall_flags_8 563 CASE ( 9 ) 564 flags => wall_flags_9 565 CASE ( 10 ) 566 flags => wall_flags_10 567 END SELECT 420 568 421 569 unroll = ( MOD( nyn_mg(l)-nys_mg(l)+1, 4 ) == 0 .AND. & … … 434 582 DO j = nys_mg(l) + 2 - colour, nyn_mg(l), 2 435 583 DO k = nzb+1, nzt_mg(l), 2 584 ! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 585 ! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 586 ! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 587 ! + f2_mg(k,l) * p_mg(k+1,j,i) & 588 ! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 589 ! ) 590 436 591 p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 437 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 438 + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 439 + f2_mg(k,l) * p_mg(k+1,j,i) & 440 + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 441 ) 592 ddx2_mg(l) * & 593 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 594 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 595 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 596 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 597 + ddy2_mg(l) * & 598 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 599 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & 600 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * & 601 ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) & 602 + f2_mg(k,l) * p_mg(k+1,j,i) & 603 + f3_mg(k,l) * & 604 ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * & 605 ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) & 606 - f_mg(k,j,i) ) 442 607 ENDDO 443 608 ENDDO … … 448 613 DO k = nzb+1, nzt_mg(l), 2 449 614 p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 450 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 451 + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 452 + f2_mg(k,l) * p_mg(k+1,j,i) & 453 + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 454 ) 615 ddx2_mg(l) * & 616 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 617 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 618 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 619 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 620 + ddy2_mg(l) * & 621 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 622 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & 623 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * & 624 ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) & 625 + f2_mg(k,l) * p_mg(k+1,j,i) & 626 + f3_mg(k,l) * & 627 ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * & 628 ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) & 629 - f_mg(k,j,i) ) 455 630 ENDDO 456 631 ENDDO … … 461 636 DO k = nzb+2, nzt_mg(l), 2 462 637 p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 463 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 464 + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 465 + f2_mg(k,l) * p_mg(k+1,j,i) & 466 + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 467 ) 638 ddx2_mg(l) * & 639 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 640 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 641 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 642 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 643 + ddy2_mg(l) * & 644 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 645 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & 646 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * & 647 ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) & 648 + f2_mg(k,l) * p_mg(k+1,j,i) & 649 + f3_mg(k,l) * & 650 ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * & 651 ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) & 652 - f_mg(k,j,i) ) 468 653 ENDDO 469 654 ENDDO … … 474 659 DO k = nzb+2, nzt_mg(l), 2 475 660 p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 476 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 477 + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 478 + f2_mg(k,l) * p_mg(k+1,j,i) & 479 + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 480 ) 661 ddx2_mg(l) * & 662 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 663 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 664 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 665 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 666 + ddy2_mg(l) * & 667 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 668 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & 669 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * & 670 ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) & 671 + f2_mg(k,l) * p_mg(k+1,j,i) & 672 + f3_mg(k,l) * & 673 ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * & 674 ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) & 675 - f_mg(k,j,i) ) 481 676 ENDDO 482 677 ENDDO … … 496 691 j = jj 497 692 p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 498 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 499 + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 500 + f2_mg(k,l) * p_mg(k+1,j,i) & 501 + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 502 ) 693 ddx2_mg(l) * & 694 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 695 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 696 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 697 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 698 + ddy2_mg(l) * & 699 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 700 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & 701 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * & 702 ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) & 703 + f2_mg(k,l) * p_mg(k+1,j,i) & 704 + f3_mg(k,l) * & 705 ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * & 706 ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) & 707 - f_mg(k,j,i) ) 503 708 j = jj+2 504 709 p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 505 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 506 + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 507 + f2_mg(k,l) * p_mg(k+1,j,i) & 508 + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 509 ) 510 ! j = jj+4 511 ! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 512 ! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 513 ! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 514 ! + f2_mg(k,l) * p_mg(k+1,j,i) & 515 ! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 516 ! ) 517 ! j = jj+6 518 ! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 519 ! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 520 ! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 521 ! + f2_mg(k,l) * p_mg(k+1,j,i) & 522 ! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 523 ! ) 710 ddx2_mg(l) * & 711 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 712 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 713 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 714 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 715 + ddy2_mg(l) * & 716 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 717 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & 718 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * & 719 ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) & 720 + f2_mg(k,l) * p_mg(k+1,j,i) & 721 + f3_mg(k,l) * & 722 ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * & 723 ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) & 724 - f_mg(k,j,i) ) 524 725 ENDDO 525 726 … … 529 730 j =jj 530 731 p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 531 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 532 + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 533 + f2_mg(k,l) * p_mg(k+1,j,i) & 534 + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 535 ) 732 ddx2_mg(l) * & 733 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 734 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 735 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 736 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 737 + ddy2_mg(l) * & 738 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 739 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & 740 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * & 741 ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) & 742 + f2_mg(k,l) * p_mg(k+1,j,i) & 743 + f3_mg(k,l) * & 744 ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * & 745 ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) & 746 - f_mg(k,j,i) ) 536 747 j = jj+2 537 748 p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 538 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 539 + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 540 + f2_mg(k,l) * p_mg(k+1,j,i) & 541 + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 542 ) 543 ! j = jj+4 544 ! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 545 ! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 546 ! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 547 ! + f2_mg(k,l) * p_mg(k+1,j,i) & 548 ! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 549 ! ) 550 ! j = jj+6 551 ! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 552 ! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 553 ! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 554 ! + f2_mg(k,l) * p_mg(k+1,j,i) & 555 ! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 556 ! ) 749 ddx2_mg(l) * & 750 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 751 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 752 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 753 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 754 + ddy2_mg(l) * & 755 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 756 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & 757 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * & 758 ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) & 759 + f2_mg(k,l) * p_mg(k+1,j,i) & 760 + f3_mg(k,l) * & 761 ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * & 762 ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) & 763 - f_mg(k,j,i) ) 557 764 ENDDO 558 765 … … 562 769 j =jj 563 770 p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 564 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 565 + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 566 + f2_mg(k,l) * p_mg(k+1,j,i) & 567 + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 568 ) 771 ddx2_mg(l) * & 772 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 773 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 774 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 775 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 776 + ddy2_mg(l) * & 777 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 778 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & 779 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * & 780 ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) & 781 + f2_mg(k,l) * p_mg(k+1,j,i) & 782 + f3_mg(k,l) * & 783 ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * & 784 ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) & 785 - f_mg(k,j,i) ) 569 786 j = jj+2 570 787 p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 571 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 572 + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 573 + f2_mg(k,l) * p_mg(k+1,j,i) & 574 + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 575 ) 576 ! j = jj+4 577 ! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 578 ! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 579 ! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 580 ! + f2_mg(k,l) * p_mg(k+1,j,i) & 581 ! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 582 ! ) 583 ! j = jj+6 584 ! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 585 ! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 586 ! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 587 ! + f2_mg(k,l) * p_mg(k+1,j,i) & 588 ! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 589 ! ) 788 ddx2_mg(l) * & 789 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 790 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 791 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 792 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 793 + ddy2_mg(l) * & 794 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 795 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & 796 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * & 797 ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) & 798 + f2_mg(k,l) * p_mg(k+1,j,i) & 799 + f3_mg(k,l) * & 800 ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * & 801 ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) & 802 - f_mg(k,j,i) ) 590 803 ENDDO 591 804 … … 595 808 j =jj 596 809 p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 597 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 598 + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 599 + f2_mg(k,l) * p_mg(k+1,j,i) & 600 + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 601 ) 810 ddx2_mg(l) * & 811 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 812 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 813 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 814 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 815 + ddy2_mg(l) * & 816 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 817 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & 818 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * & 819 ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) & 820 + f2_mg(k,l) * p_mg(k+1,j,i) & 821 + f3_mg(k,l) * & 822 ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * & 823 ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) & 824 - f_mg(k,j,i) ) 602 825 j = jj+2 603 826 p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 604 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 605 + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 606 + f2_mg(k,l) * p_mg(k+1,j,i) & 607 + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 608 ) 609 ! j = jj+4 610 ! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 611 ! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 612 ! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 613 ! + f2_mg(k,l) * p_mg(k+1,j,i) & 614 ! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 615 ! ) 616 ! j = jj+6 617 ! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 618 ! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 619 ! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 620 ! + f2_mg(k,l) * p_mg(k+1,j,i) & 621 ! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 622 ! ) 827 ddx2_mg(l) * & 828 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 829 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 830 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 831 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 832 + ddy2_mg(l) * & 833 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 834 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & 835 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * & 836 ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) & 837 + f2_mg(k,l) * p_mg(k+1,j,i) & 838 + f3_mg(k,l) * & 839 ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * & 840 ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) & 841 - f_mg(k,j,i) ) 623 842 ENDDO 624 843 … … 669 888 ENDDO 670 889 890 ! 891 !-- Set pressure within topography and at the topography surfaces 892 !$OMP PARALLEL PRIVATE (i,j,k,wall_left,wall_north,wall_right,wall_south,wall_top,wall_total) 893 !$OMP DO 894 DO i = nxl_mg(l), nxr_mg(l) 895 DO j = nys_mg(l), nyn_mg(l) 896 DO k = nzb, nzt_mg(l) 897 ! 898 !-- First, set pressure inside topography to zero 899 p_mg(k,j,i) = p_mg(k,j,i) * ( 1.0 - IBITS( flags(k,j,i), 6, 1 ) ) 900 ! 901 !-- Second, determine if the gridpoint inside topography is adjacent 902 !-- to a wall and set its value to a value given by the average of 903 !-- those values obtained from Neumann boundary condition 904 wall_left = IBITS( flags(k,j,i-1), 5, 1 ) 905 wall_right = IBITS( flags(k,j,i+1), 4, 1 ) 906 wall_south = IBITS( flags(k,j-1,i), 3, 1 ) 907 wall_north = IBITS( flags(k,j+1,i), 2, 1 ) 908 wall_top = IBITS( flags(k+1,j,i), 0, 1 ) 909 wall_total = wall_left + wall_right + wall_south + wall_north + & 910 wall_top 911 912 IF ( wall_total > 0.0 ) THEN 913 p_mg(k,j,i) = 1.0 / wall_total * & 914 ( wall_left * p_mg(k,j,i-1) + & 915 wall_right * p_mg(k,j,i+1) + & 916 wall_south * p_mg(k,j-1,i) + & 917 wall_north * p_mg(k,j+1,i) + & 918 wall_top * p_mg(k+1,j,i) ) 919 ENDIF 920 ENDDO 921 ENDDO 922 ENDDO 923 !$OMP END PARALLEL 924 925 ! 926 !-- One more time horizontal boundary conditions 927 CALL exchange_horiz( p_mg ) 671 928 672 929 END SUBROUTINE redblack -
palm/trunk/SOURCE/user_interface.f90
r110 r114 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! +argument gls in user_init_grid 7 7 ! 8 8 ! Former revisions: … … 238 238 239 239 240 SUBROUTINE user_init_grid( nzb_local )240 SUBROUTINE user_init_grid( gls, nzb_local ) 241 241 242 242 !------------------------------------------------------------------------------! … … 245 245 ! ------------ 246 246 ! Execution of user-defined grid initializing actions 247 ! First argument gls contains the number of ghost layers, which is > 1 if the 248 ! multigrid method for the pressure solver is used 247 249 !------------------------------------------------------------------------------! 248 250 … … 253 255 IMPLICIT NONE 254 256 255 INTEGER, DIMENSION(-1:ny+1,-1:nx+1) :: nzb_local 257 INTEGER :: gls 258 259 INTEGER, DIMENSION(-gls:ny+gls,-gls:nx+gls) :: nzb_local 256 260 257 261 !
Note: See TracChangeset
for help on using the changeset viewer.
