Changeset 1353 for palm/trunk/SOURCE/init_grid.f90
 Timestamp:
 Apr 8, 2014 3:21:23 PM (10 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SOURCE/init_grid.f90
r1323 r1353 20 20 ! Current revisions: 21 21 !  22 ! 22 ! REAL constants provided with KINDattribute 23 23 ! 24 24 ! Former revisions: … … 105 105 f2_mg, f3_mg, l_grid, l_wall, zu, zw 106 106 107 USE control_parameters, 107 USE control_parameters, & 108 108 ONLY: bc_lr, bc_ns, building_height, building_length_x, & 109 109 building_length_y, building_wall_left, building_wall_south, & … … 201 201 ! 202 202 ! Allocate grid arrays 203 ALLOCATE( ddzu(1:nzt+1), ddzw(1:nzt+1), dd2zu(1:nzt), dzu(1:nzt+1), &203 ALLOCATE( ddzu(1:nzt+1), ddzw(1:nzt+1), dd2zu(1:nzt), dzu(1:nzt+1), & 204 204 dzw(1:nzt+1), l_grid(1:nzt), zu(nzb:nzt+1), zw(nzb:nzt+1) ) 205 205 206 206 ! 207 207 ! Compute height of ulevels from constant grid length and dz stretch factors 208 IF ( dz == 1.0 ) THEN208 IF ( dz == 1.0_wp ) THEN 209 209 message_string = 'missing dz' 210 210 CALL message( 'init_grid', 'PA0200', 1, 2, 0, 6, 0 ) 211 ELSEIF ( dz <= 0.0 ) THEN211 ELSEIF ( dz <= 0.0_wp ) THEN 212 212 WRITE( message_string, * ) 'dz=',dz,' <= 0.0' 213 213 CALL message( 'init_grid', 'PA0201', 1, 2, 0, 6, 0 ) … … 223 223 224 224 IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2 ) THEN 225 zu(0) = 0.0 226 ! zu(0) =  dz * 0.5 225 zu(0) = 0.0_wp 226 ! zu(0) =  dz * 0.5_wp 227 227 ELSE 228 zu(0) =  dz * 0.5 228 zu(0) =  dz * 0.5_wp 229 229 ENDIF 230 zu(1) = dz * 0.5 230 zu(1) = dz * 0.5_wp 231 231 232 232 dz_stretch_level_index = nzt+1 … … 246 246 ! ground the first u and wlevel (k=0) are defined at same height (z=0). 247 247 ! The top wlevel is extrapolated linearly. 248 zw(0) = 0.0 248 zw(0) = 0.0_wp 249 249 DO k = 1, nzt 250 zw(k) = ( zu(k) + zu(k+1) ) * 0.5 251 ENDDO 252 zw(nzt+1) = zw(nzt) + 2.0 * ( zu(nzt+1)  zw(nzt) )250 zw(k) = ( zu(k) + zu(k+1) ) * 0.5_wp 251 ENDDO 252 zw(nzt+1) = zw(nzt) + 2.0_wp * ( zu(nzt+1)  zw(nzt) ) 253 253 254 254 ELSE … … 259 259 ! wlevel are defined at same height, but staggered from the second level. 260 260 ! The second ulevel (k=1) corresponds to the top of the Prandtllayer. 261 zu(nzt+1) = dz * 0.5 262 zu(nzt) =  dz * 0.5 261 zu(nzt+1) = dz * 0.5_wp 262 zu(nzt) =  dz * 0.5_wp 263 263 264 264 dz_stretch_level_index = 0 … … 281 281 ! consistency, since w and all scalar variables are defined up tp nzt+1. 282 282 zw(nzt+1) = dz 283 zw(nzt) = 0.0 283 zw(nzt) = 0.0_wp 284 284 DO k = 0, nzt 285 zw(k) = ( zu(k) + zu(k+1) ) * 0.5 285 zw(k) = ( zu(k) + zu(k+1) ) * 0.5_wp 286 286 ENDDO 287 287 … … 299 299 DO k = 1, nzt+1 300 300 dzu(k) = zu(k)  zu(k1) 301 ddzu(k) = 1.0 / dzu(k)301 ddzu(k) = 1.0_wp / dzu(k) 302 302 dzw(k) = zw(k)  zw(k1) 303 ddzw(k) = 1.0 / dzw(k)303 ddzw(k) = 1.0_wp / dzw(k) 304 304 ENDDO 305 305 306 306 DO k = 1, nzt 307 dd2zu(k) = 1.0 / ( dzu(k) + dzu(k+1) )307 dd2zu(k) = 1.0_wp / ( dzu(k) + dzu(k+1) ) 308 308 ENDDO 309 309 … … 340 340 nzt_l = nzt 341 341 DO l = maximum_grid_level1, 1, 1 342 dzu_mg(nzb+1,l) = 2.0 * dzu_mg(nzb+1,l+1)343 dzw_mg(nzb+1,l) = 2.0 * dzw_mg(nzb+1,l+1)342 dzu_mg(nzb+1,l) = 2.0_wp * dzu_mg(nzb+1,l+1) 343 dzw_mg(nzb+1,l) = 2.0_wp * dzw_mg(nzb+1,l+1) 344 344 nzt_l = nzt_l / 2 345 345 DO k = 2, nzt_l+1 … … 353 353 dy_l = dy 354 354 DO l = maximum_grid_level, 1, 1 355 ddx2_mg(l) = 1.0 / dx_l**2356 ddy2_mg(l) = 1.0 / dy_l**2355 ddx2_mg(l) = 1.0_wp / dx_l**2 356 ddy2_mg(l) = 1.0_wp / dy_l**2 357 357 DO k = nzb+1, nzt_l 358 f2_mg(k,l) = 1.0 / ( dzu_mg(k+1,l) * dzw_mg(k,l) )359 f3_mg(k,l) = 1.0 / ( dzu_mg(k,l) * dzw_mg(k,l) )360 f1_mg(k,l) = 2.0 * ( ddx2_mg(l) + ddy2_mg(l) ) + &358 f2_mg(k,l) = 1.0_wp / ( dzu_mg(k+1,l) * dzw_mg(k,l) ) 359 f3_mg(k,l) = 1.0_wp / ( dzu_mg(k,l) * dzw_mg(k,l) ) 360 f1_mg(k,l) = 2.0_wp * ( ddx2_mg(l) + ddy2_mg(l) ) + & 361 361 f2_mg(k,l) + f3_mg(k,l) 362 362 ENDDO 363 363 nzt_l = nzt_l / 2 364 dx_l = dx_l * 2.0 365 dy_l = dy_l * 2.0 364 dx_l = dx_l * 2.0_wp 365 dy_l = dy_l * 2.0_wp 366 366 ENDDO 367 367 … … 370 370 ! 371 371 ! Compute the reciprocal values of the horizontal grid lengths. 372 ddx = 1.0 / dx373 ddy = 1.0 / dy372 ddx = 1.0_wp / dx 373 ddy = 1.0_wp / dy 374 374 dx2 = dx * dx 375 375 dy2 = dy * dy 376 ddx2 = 1.0 / dx2377 ddy2 = 1.0 / dy2376 ddx2 = 1.0_wp / dx2 377 ddy2 = 1.0_wp / dy2 378 378 379 379 ! … … 433 433 nzb_w_inner = nzb; nzb_w_outer = nzb 434 434 435 rflags_s_inner = 1.0 436 rflags_invers = 1.0 435 rflags_s_inner = 1.0_wp 436 rflags_invers = 1.0_wp 437 437 438 438 ! … … 453 453 nzb_diff_u = nzb_diff; nzb_diff_v = nzb_diff 454 454 455 wall_e_x = 0.0 ; wall_e_y = 0.0; wall_u = 0.0; wall_v = 0.0456 wall_w_x = 0.0 ; wall_w_y = 0.0457 fwxp = 1.0 ; fwxm = 1.0; fwyp = 1.0; fwym = 1.0458 fxp = 1.0 ; fxm = 1.0; fyp = 1.0; fym = 1.0455 wall_e_x = 0.0_wp; wall_e_y = 0.0_wp; wall_u = 0.0_wp; wall_v = 0.0_wp 456 wall_w_x = 0.0_wp; wall_w_y = 0.0_wp 457 fwxp = 1.0_wp; fwxm = 1.0_wp; fwyp = 1.0_wp; fwym = 1.0_wp 458 fxp = 1.0_wp; fxm = 1.0_wp; fyp = 1.0_wp; fym = 1.0_wp 459 459 460 460 ! … … 471 471 DO k = 1, nzt 472 472 vertical_influence(k) = MIN ( INT( l_grid(k) / & 473 ( wall_adjustment_factor * dzw(k) ) + 0.5 ), nzt  k )473 ( wall_adjustment_factor * dzw(k) ) + 0.5_wp ), nzt  k ) 474 474 ENDDO 475 475 476 476 DO k = 1, MAXVAL( nzb_s_inner ) 477 IF ( l_grid(k) > 1.5 * dx * wall_adjustment_factor .OR. &478 l_grid(k) > 1.5 * dy * wall_adjustment_factor ) THEN477 IF ( l_grid(k) > 1.5_wp * dx * wall_adjustment_factor .OR. & 478 l_grid(k) > 1.5_wp * dy * wall_adjustment_factor ) THEN 479 479 WRITE( message_string, * ) 'grid anisotropy exceeds ', & 480 480 'threshold given by only local', & … … 582 582 ! Street canyon size has to meet some requirements 583 583 IF ( canyon_width_x /= 9999999.9_wp ) THEN 584 IF ( ( cxl < 1 ) .OR. ( cxr > nx1 ) .OR. ( cwx < 3 ) .OR. &584 IF ( ( cxl < 1 ) .OR. ( cxr > nx1 ) .OR. ( cwx < 3 ) .OR. & 585 585 ( ch < 3 ) ) THEN 586 WRITE( message_string, * ) 'inconsistent canyon parameters:', &587 '&cxl=', cxl, 'cxr=', cxr, &588 'cwx=', cwx, &586 WRITE( message_string, * ) 'inconsistent canyon parameters:', & 587 '&cxl=', cxl, 'cxr=', cxr, & 588 'cwx=', cwx, & 589 589 'ch=', ch, 'nx=', nx, 'ny=', ny 590 590 CALL message( 'init_grid', 'PA0205', 1, 2, 0, 6, 0 ) 591 591 ENDIF 592 592 ELSEIF ( canyon_width_y /= 9999999.9_wp ) THEN 593 IF ( ( cys < 1 ) .OR. ( cyn > ny1 ) .OR. ( cwy < 3 ) .OR. &593 IF ( ( cys < 1 ) .OR. ( cyn > ny1 ) .OR. ( cwy < 3 ) .OR. & 594 594 ( ch < 3 ) ) THEN 595 WRITE( message_string, * ) 'inconsistent canyon parameters:', &596 '&cys=', cys, 'cyn=', cyn, &597 'cwy=', cwy, &595 WRITE( message_string, * ) 'inconsistent canyon parameters:', & 596 '&cys=', cys, 'cyn=', cyn, & 597 'cwy=', cwy, & 598 598 'ch=', ch, 'nx=', nx, 'ny=', ny 599 599 CALL message( 'init_grid', 'PA0206', 1, 2, 0, 6, 0 ) 600 600 ENDIF 601 601 ENDIF 602 IF ( canyon_width_x /= 9999999.9_wp .AND. canyon_width_y /= 9999999.9_wp ) &603 THEN604 message_string = 'inconsistent canyon parameters:' // &605 '&street canyon can only be oriented' // &602 IF ( canyon_width_x /= 9999999.9_wp .AND. & 603 canyon_width_y /= 9999999.9_wp ) THEN 604 message_string = 'inconsistent canyon parameters:' // & 605 '&street canyon can only be oriented' // & 606 606 '&either in x or in ydirection' 607 607 CALL message( 'init_grid', 'PA0207', 1, 2, 0, 6, 0 ) … … 679 679 ! scheme have to be reduced up to nzt (required at the lateral boundaries). 680 680 nzb_max = MAXVAL( nzb_local ) 681 IF ( inflow_l .OR. outflow_l .OR. inflow_r .OR. outflow_r .OR. &681 IF ( inflow_l .OR. outflow_l .OR. inflow_r .OR. outflow_r .OR. & 682 682 inflow_n .OR. outflow_n .OR. inflow_s .OR. outflow_s ) THEN 683 683 nzb_max = nzt … … 693 693 ! Consistency checks 694 694 IF ( MINVAL( nzb_local ) < 0 .OR. MAXVAL( nzb_local ) > nz + 1 ) THEN 695 WRITE( message_string, * ) 'nzb_local values are outside the', &696 'model domain', &697 '&MINVAL( nzb_local ) = ', MINVAL(nzb_local), &695 WRITE( message_string, * ) 'nzb_local values are outside the', & 696 'model domain', & 697 '&MINVAL( nzb_local ) = ', MINVAL(nzb_local), & 698 698 '&MAXVAL( nzb_local ) = ', MAXVAL(nzb_local) 699 699 CALL message( 'init_grid', 'PA0210', 1, 2, 0, 6, 0 ) … … 701 701 702 702 IF ( bc_lr == 'cyclic' ) THEN 703 IF ( ANY( nzb_local(:,1) /= nzb_local(:,nx) ) .OR. &703 IF ( ANY( nzb_local(:,1) /= nzb_local(:,nx) ) .OR. & 704 704 ANY( nzb_local(:,0) /= nzb_local(:,nx+1) ) ) THEN 705 message_string = 'nzb_local does not fulfill cyclic' // &705 message_string = 'nzb_local does not fulfill cyclic' // & 706 706 ' boundary condition in xdirection' 707 707 CALL message( 'init_grid', 'PA0211', 1, 2, 0, 6, 0 ) … … 709 709 ENDIF 710 710 IF ( bc_ns == 'cyclic' ) THEN 711 IF ( ANY( nzb_local(1,:) /= nzb_local(ny,:) ) .OR. &711 IF ( ANY( nzb_local(1,:) /= nzb_local(ny,:) ) .OR. & 712 712 ANY( nzb_local(0,:) /= nzb_local(ny+1,:) ) ) THEN 713 message_string = 'nzb_local does not fulfill cyclic' // &713 message_string = 'nzb_local does not fulfill cyclic' // & 714 714 ' boundary condition in ydirection' 715 715 CALL message( 'init_grid', 'PA0212', 1, 2, 0, 6, 0 ) … … 770 770 DO j = 1, ny + 1 771 771 DO i = 0, nx 772 nzb_tmp(j,i) = MAX( nzb_local(j,i1), nzb_local(j,i), &772 nzb_tmp(j,i) = MAX( nzb_local(j,i1), nzb_local(j,i), & 773 773 nzb_local(j,i+1) ) 774 774 ENDDO … … 776 776 DO i = nxl, nxr 777 777 DO j = nys, nyn 778 nzb_s_outer(j,i) = MAX( nzb_tmp(j1,i), nzb_tmp(j,i), &778 nzb_s_outer(j,i) = MAX( nzb_tmp(j1,i), nzb_tmp(j,i), & 779 779 nzb_tmp(j+1,i) ) 780 780 ENDDO … … 812 812 DO i = nxl, nxr 813 813 DO j = nys, nyn 814 nzb_u_outer(j,i) = MAX( nzb_tmp(j1,i), nzb_tmp(j,i), &814 nzb_u_outer(j,i) = MAX( nzb_tmp(j1,i), nzb_tmp(j,i), & 815 815 nzb_tmp(j+1,i) ) 816 816 ENDDO … … 844 844 DO j = nys, nyn 845 845 DO i = nxl, nxr 846 nzb_v_outer(j,i) = MAX( nzb_tmp(j,i1), nzb_tmp(j,i), &846 nzb_v_outer(j,i) = MAX( nzb_tmp(j,i1), nzb_tmp(j,i), & 847 847 nzb_tmp(j,i+1) ) 848 848 ENDDO … … 891 891 DO j = nysg, nyng 892 892 DO k = nzb, nzt+1 893 IF ( k <= nzb_s_inner(j,i) ) rflags_s_inner(k,j,i) = 0.0 894 IF ( k <= nzb_s_inner(j,i) ) rflags_invers(j,i,k) = 0.0 893 IF ( k <= nzb_s_inner(j,i) ) rflags_s_inner(k,j,i) = 0.0_wp 894 IF ( k <= nzb_s_inner(j,i) ) rflags_invers(j,i,k) = 0.0_wp 895 895 ENDDO 896 896 ENDDO … … 938 938 ! ucomponent 939 939 IF ( nzb_u_outer(j,i) > nzb_u_outer(j+1,i) ) THEN 940 wall_u(j,i) = 1.0 ! north wall (location of adjacent fluid)941 fym(j,i) = 0.0 942 fyp(j,i) = 1.0 940 wall_u(j,i) = 1.0_wp ! north wall (location of adjacent fluid) 941 fym(j,i) = 0.0_wp 942 fyp(j,i) = 1.0_wp 943 943 ELSEIF ( nzb_u_outer(j,i) > nzb_u_outer(j1,i) ) THEN 944 wall_u(j,i) = 1.0 ! south wall (location of adjacent fluid)945 fym(j,i) = 1.0 946 fyp(j,i) = 0.0 944 wall_u(j,i) = 1.0_wp ! south wall (location of adjacent fluid) 945 fym(j,i) = 1.0_wp 946 fyp(j,i) = 0.0_wp 947 947 ENDIF 948 948 ! 949 949 ! vcomponent 950 950 IF ( nzb_v_outer(j,i) > nzb_v_outer(j,i+1) ) THEN 951 wall_v(j,i) = 1.0 ! rigth wall (location of adjacent fluid)952 fxm(j,i) = 0.0 953 fxp(j,i) = 1.0 951 wall_v(j,i) = 1.0_wp ! rigth wall (location of adjacent fluid) 952 fxm(j,i) = 0.0_wp 953 fxp(j,i) = 1.0_wp 954 954 ELSEIF ( nzb_v_outer(j,i) > nzb_v_outer(j,i1) ) THEN 955 wall_v(j,i) = 1.0 ! left wall (location of adjacent fluid)956 fxm(j,i) = 1.0 957 fxp(j,i) = 0.0 955 wall_v(j,i) = 1.0_wp ! left wall (location of adjacent fluid) 956 fxm(j,i) = 1.0_wp 957 fxp(j,i) = 0.0_wp 958 958 ENDIF 959 959 ! … … 961 961 ! production of tke 962 962 IF ( nzb_w_outer(j,i) > nzb_w_outer(j+1,i) ) THEN 963 wall_e_y(j,i) = 1.0 ! north wall (location of adjacent fluid)964 wall_w_y(j,i) = 1.0 965 fwym(j,i) = 0.0 966 fwyp(j,i) = 1.0 963 wall_e_y(j,i) = 1.0_wp ! north wall (location of adjacent fluid) 964 wall_w_y(j,i) = 1.0_wp 965 fwym(j,i) = 0.0_wp 966 fwyp(j,i) = 1.0_wp 967 967 ELSEIF ( nzb_w_outer(j,i) > nzb_w_outer(j1,i) ) THEN 968 wall_e_y(j,i) = 1.0 ! south wall (location of adjacent fluid)969 wall_w_y(j,i) = 1.0 970 fwym(j,i) = 1.0 971 fwyp(j,i) = 0.0 968 wall_e_y(j,i) = 1.0_wp ! south wall (location of adjacent fluid) 969 wall_w_y(j,i) = 1.0_wp 970 fwym(j,i) = 1.0_wp 971 fwyp(j,i) = 0.0_wp 972 972 ENDIF 973 973 IF ( nzb_w_outer(j,i) > nzb_w_outer(j,i+1) ) THEN 974 wall_e_x(j,i) = 1.0 ! right wall (location of adjacent fluid)975 wall_w_x(j,i) = 1.0 976 fwxm(j,i) = 0.0 977 fwxp(j,i) = 1.0 974 wall_e_x(j,i) = 1.0_wp ! right wall (location of adjacent fluid) 975 wall_w_x(j,i) = 1.0_wp 976 fwxm(j,i) = 0.0_wp 977 fwxp(j,i) = 1.0_wp 978 978 ELSEIF ( nzb_w_outer(j,i) > nzb_w_outer(j,i1) ) THEN 979 wall_e_x(j,i) = 1.0 ! left wall (location of adjacent fluid)980 wall_w_x(j,i) = 1.0 981 fwxm(j,i) = 1.0 982 fwxp(j,i) = 0.0 979 wall_e_x(j,i) = 1.0_wp ! left wall (location of adjacent fluid) 980 wall_w_x(j,i) = 1.0_wp 981 fwxm(j,i) = 1.0_wp 982 fwxp(j,i) = 0.0_wp 983 983 ENDIF 984 984 ! … … 1399 1399 ! North wall (y distance) 1400 1400 DO k = wall_n(j,i), nzb_si 1401 l_wall(k,j+1,i) = MIN( l_wall(k,j+1,i), 0.5 * dy )1401 l_wall(k,j+1,i) = MIN( l_wall(k,j+1,i), 0.5_wp * dy ) 1402 1402 ENDDO 1403 1403 ! 1404 1404 ! Above North wall (yz distance) 1405 1405 DO k = nzb_si + 1, nzb_si + vi 1406 l_wall(k,j+1,i) = MIN( l_wall(k,j+1,i), &1407 SQRT( 0.25 * dy**2 +&1406 l_wall(k,j+1,i) = MIN( l_wall(k,j+1,i), & 1407 SQRT( 0.25_wp * dy**2 + & 1408 1408 ( zu(k)  zw(nzb_si) )**2 ) ) 1409 1409 ENDDO … … 1413 1413 DO k = corner_nl(j,i), nzb_si 1414 1414 l_wall(k,j+1,i1) = MIN( l_wall(k,j+1,i1), & 1415 0.5 * SQRT( dx**2 + dy**2 ) )1415 0.5_wp * SQRT( dx**2 + dy**2 ) ) 1416 1416 ENDDO 1417 1417 ! 1418 1418 ! Above Northleft corner (xyz distance) 1419 1419 DO k = nzb_si + 1, nzb_si + vi 1420 l_wall(k,j+1,i1) = MIN( l_wall(k,j+1,i1), &1421 SQRT( 0.25 * (dx**2 + dy**2) +&1422 1420 l_wall(k,j+1,i1) = MIN( l_wall(k,j+1,i1), & 1421 SQRT( 0.25_wp * (dx**2 + dy**2) + & 1422 ( zu(k)  zw(nzb_si) )**2 ) ) 1423 1423 ENDDO 1424 1424 ENDIF … … 1427 1427 IF ( corner_nr(j,i) > 0 ) THEN 1428 1428 DO k = corner_nr(j,i), nzb_si 1429 l_wall(k,j+1,i+1) = MIN( l_wall(k,j+1,i+1), &1430 0.5 * SQRT( dx**2 + dy**2 ) )1429 l_wall(k,j+1,i+1) = MIN( l_wall(k,j+1,i+1), & 1430 0.5_wp * SQRT( dx**2 + dy**2 ) ) 1431 1431 ENDDO 1432 1432 ! 1433 1433 ! Above northright corner (xyz distance) 1434 1434 DO k = nzb_si + 1, nzb_si + vi 1435 l_wall(k,j+1,i+1) = MIN( l_wall(k,j+1,i+1), &1436 SQRT( 0.25 * (dx**2 + dy**2) +&1437 1435 l_wall(k,j+1,i+1) = MIN( l_wall(k,j+1,i+1), & 1436 SQRT( 0.25_wp * (dx**2 + dy**2) + & 1437 ( zu(k)  zw(nzb_si) )**2 ) ) 1438 1438 ENDDO 1439 1439 ENDIF … … 1444 1444 ! South wall (y distance) 1445 1445 DO k = wall_s(j,i), nzb_si 1446 l_wall(k,j1,i) = MIN( l_wall(k,j1,i), 0.5 * dy )1446 l_wall(k,j1,i) = MIN( l_wall(k,j1,i), 0.5_wp * dy ) 1447 1447 ENDDO 1448 1448 ! 1449 1449 ! Above south wall (yz distance) 1450 DO k = nzb_si + 1, & 1451 nzb_si + vi 1452 l_wall(k,j1,i) = MIN( l_wall(k,j1,i), & 1453 SQRT( 0.25 * dy**2 + & 1450 DO k = nzb_si + 1, nzb_si + vi 1451 l_wall(k,j1,i) = MIN( l_wall(k,j1,i), & 1452 SQRT( 0.25_wp * dy**2 + & 1454 1453 ( zu(k)  zw(nzb_si) )**2 ) ) 1455 1454 ENDDO … … 1458 1457 IF ( corner_sl(j,i) > 0 ) THEN 1459 1458 DO k = corner_sl(j,i), nzb_si 1460 l_wall(k,j1,i1) = MIN( l_wall(k,j1,i1), &1461 0.5 * SQRT( dx**2 + dy**2 ) )1459 l_wall(k,j1,i1) = MIN( l_wall(k,j1,i1), & 1460 0.5_wp * SQRT( dx**2 + dy**2 ) ) 1462 1461 ENDDO 1463 1462 ! 1464 1463 ! Above southleft corner (xyz distance) 1465 1464 DO k = nzb_si + 1, nzb_si + vi 1466 l_wall(k,j1,i1) = MIN( l_wall(k,j1,i1), &1467 SQRT( 0.25 * (dx**2 + dy**2) +&1468 1465 l_wall(k,j1,i1) = MIN( l_wall(k,j1,i1), & 1466 SQRT( 0.25_wp * (dx**2 + dy**2) + & 1467 ( zu(k)  zw(nzb_si) )**2 ) ) 1469 1468 ENDDO 1470 1469 ENDIF … … 1473 1472 IF ( corner_sr(j,i) > 0 ) THEN 1474 1473 DO k = corner_sr(j,i), nzb_si 1475 l_wall(k,j1,i+1) = MIN( l_wall(k,j1,i+1), &1476 0.5 * SQRT( dx**2 + dy**2 ) )1474 l_wall(k,j1,i+1) = MIN( l_wall(k,j1,i+1), & 1475 0.5_wp * SQRT( dx**2 + dy**2 ) ) 1477 1476 ENDDO 1478 1477 ! 1479 1478 ! Above southright corner (xyz distance) 1480 1479 DO k = nzb_si + 1, nzb_si + vi 1481 l_wall(k,j1,i+1) = MIN( l_wall(k,j1,i+1), &1482 SQRT( 0.25 * (dx**2 + dy**2) +&1483 1480 l_wall(k,j1,i+1) = MIN( l_wall(k,j1,i+1), & 1481 SQRT( 0.25_wp * (dx**2 + dy**2) + & 1482 ( zu(k)  zw(nzb_si) )**2 ) ) 1484 1483 ENDDO 1485 1484 ENDIF … … 1491 1490 ! Left wall (x distance) 1492 1491 DO k = wall_l(j,i), nzb_si 1493 l_wall(k,j,i1) = MIN( l_wall(k,j,i1), 0.5 * dx )1492 l_wall(k,j,i1) = MIN( l_wall(k,j,i1), 0.5_wp * dx ) 1494 1493 ENDDO 1495 1494 ! 1496 1495 ! Above left wall (xz distance) 1497 1496 DO k = nzb_si + 1, nzb_si + vi 1498 l_wall(k,j,i1) = MIN( l_wall(k,j,i1), & 1499 SQRT( 0.25 * dx**2 + & 1497 l_wall(k,j,i1) = MIN( l_wall(k,j,i1), & 1498 SQRT( 0.25_wp * dx**2 + & 1499 ( zu(k)  zw(nzb_si) )**2 ) ) 1500 ENDDO 1501 ENDIF 1502 1503 IF ( wall_r(j,i) > 0 ) THEN 1504 ! 1505 ! Right wall (x distance) 1506 DO k = wall_r(j,i), nzb_si 1507 l_wall(k,j,i+1) = MIN( l_wall(k,j,i+1), 0.5_wp * dx ) 1508 ENDDO 1509 ! 1510 ! Above right wall (xz distance) 1511 DO k = nzb_si + 1, nzb_si + vi 1512 l_wall(k,j,i+1) = MIN( l_wall(k,j,i+1), & 1513 SQRT( 0.25_wp * dx**2 + & 1500 1514 ( zu(k)  zw(nzb_si) )**2 ) ) 1501 1515 ENDDO 1502 ENDIF1503 1504 IF ( wall_r(j,i) > 0 ) THEN1505 !1506 ! Right wall (x distance)1507 DO k = wall_r(j,i), nzb_si1508 l_wall(k,j,i+1) = MIN( l_wall(k,j,i+1), 0.5 * dx )1509 ENDDO1510 !1511 ! Above right wall (xz distance)1512 DO k = nzb_si + 1, nzb_si + vi1513 l_wall(k,j,i+1) = MIN( l_wall(k,j,i+1), &1514 SQRT( 0.25 * dx**2 + &1515 ( zu(k)  zw(nzb_si) )**2 ) )1516 ENDDO1517 1516 1518 1517 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.