Changeset 2698 for palm/trunk/SOURCE
- Timestamp:
- Dec 14, 2017 6:46:24 PM (6 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 14 edited
-
date_and_time_mod.f90 (modified) (2 diffs)
-
header.f90 (modified) (3 diffs)
-
init_3d_model.f90 (modified) (4 diffs)
-
init_grid.f90 (modified) (4 diffs)
-
lpm.f90 (modified) (3 diffs)
-
lpm_advec.f90 (modified) (10 diffs)
-
lpm_boundary_conds.f90 (modified) (22 diffs)
-
lpm_init.f90 (modified) (4 diffs)
-
microphysics_mod.f90 (modified) (5 diffs)
-
plant_canopy_model_mod.f90 (modified) (19 diffs)
-
pmc_interface_mod.f90 (modified) (39 diffs)
-
radiation_model_mod.f90 (modified) (17 diffs)
-
surface_mod.f90 (modified) (2 diffs)
-
turbulence_closure_mod.f90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/date_and_time_mod.f90
r2696 r2698 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Bugfix in definition of d_seconds_year. 23 23 ! 24 24 ! Former revisions: … … 56 56 REAL(wp) :: time_utc_init = 43200.0_wp !< UTC time at model start 57 57 58 REAL(wp), PARAMETER :: d_hours_day = 0.0416666666667_wp!< inverse of hours per day (1/24)59 REAL(wp), PARAMETER :: d_seconds_hour = 0.000277777777778_wp!< inverse of seconds per hour (1/3600)60 REAL(wp), PARAMETER :: d_seconds_year = 31536000.0_wp!< inverse of the seconds per year (1/(365*86400))58 REAL(wp), PARAMETER :: d_hours_day = 1.0_wp / 24.0_wp !< inverse of hours per day (1/24) 59 REAL(wp), PARAMETER :: d_seconds_hour = 1.0_wp / 3600.0_wp !< inverse of seconds per hour (1/3600) 60 REAL(wp), PARAMETER :: d_seconds_year = 1.0_wp / 31536000.0_wp !< inverse of the seconds per year (1/(365*86400)) 61 61 62 62 SAVE -
palm/trunk/SOURCE/header.f90
r2696 r2698 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Bugfix in get_topography_top_index 23 23 ! 24 24 ! Former revisions: … … 410 410 411 411 USE surface_mod, & 412 ONLY: surf_def_h, get_topography_top_index 412 ONLY: surf_def_h, get_topography_top_index_ji 413 413 414 414 USE synthetic_turbulence_generator_mod, & … … 993 993 WRITE( io, 280 ) 994 994 IF ( turbulent_inflow ) THEN 995 WRITE( io, 281 ) zu( get_topography_top_index ( 0, 0, 's' ) )995 WRITE( io, 281 ) zu( get_topography_top_index_ji( 0, 0, 's' ) ) 996 996 ENDIF 997 997 IF ( TRIM( initializing_actions ) == 'cyclic_fill' ) THEN -
palm/trunk/SOURCE/init_3d_model.f90
r2696 r2698 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Bugfix in get_topography_top_index 23 23 ! 24 24 ! Former revisions: … … 466 466 USE surface_mod, & 467 467 ONLY : init_surface_arrays, init_surfaces, surf_def_h, surf_lsm_h, & 468 surf_usm_h, get_topography_top_index 468 surf_usm_h, get_topography_top_index_ji 469 469 470 470 USE transpose_indices … … 1580 1580 DO i = nxlg, nxrg 1581 1581 DO j = nysg, nyng 1582 nz_u_shift = get_topography_top_index ( j, i, 'u' )1583 nz_v_shift = get_topography_top_index ( j, i, 'v' )1584 nz_w_shift = get_topography_top_index ( j, i, 'w' )1585 nz_s_shift = get_topography_top_index ( j, i, 's' )1582 nz_u_shift = get_topography_top_index_ji( j, i, 'u' ) 1583 nz_v_shift = get_topography_top_index_ji( j, i, 'v' ) 1584 nz_w_shift = get_topography_top_index_ji( j, i, 'w' ) 1585 nz_s_shift = get_topography_top_index_ji( j, i, 's' ) 1586 1586 1587 1587 u(nz_u_shift:nzt+1,j,i) = u(0:nzt+1-nz_u_shift,j,i) … … 1625 1625 IF ( complex_terrain ) THEN 1626 1626 IF ( nxlg <= 0 .AND. nxrg >= 0 .AND. nysg <= 0 .AND. nyng >= 0 ) THEN 1627 nz_u_shift_l = get_topography_top_index ( 0, 0, 'u' )1628 nz_v_shift_l = get_topography_top_index ( 0, 0, 'v' )1629 nz_w_shift_l = get_topography_top_index ( 0, 0, 'w' )1630 nz_s_shift_l = get_topography_top_index ( 0, 0, 's' )1627 nz_u_shift_l = get_topography_top_index_ji( 0, 0, 'u' ) 1628 nz_v_shift_l = get_topography_top_index_ji( 0, 0, 'v' ) 1629 nz_w_shift_l = get_topography_top_index_ji( 0, 0, 'w' ) 1630 nz_s_shift_l = get_topography_top_index_ji( 0, 0, 's' ) 1631 1631 ELSE 1632 1632 nz_u_shift_l = 0 -
palm/trunk/SOURCE/init_grid.f90
r2696 r2698 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Bugfix in get_topography_top_index 23 23 ! 24 24 ! Former revisions: … … 298 298 299 299 USE surface_mod, & 300 ONLY: get_topography_top_index, init_bc300 ONLY: get_topography_top_index, get_topography_top_index_ji, init_bc 301 301 302 302 USE vertical_nesting_mod, & … … 552 552 !-- Topography height on scalar grid. Therefore, determine index of 553 553 !-- upward-facing surface element on scalar grid. 554 zu_s_inner(i,j) = zu( get_topography_top_index ( j, i, 's' ) )554 zu_s_inner(i,j) = zu( get_topography_top_index_ji( j, i, 's' ) ) 555 555 ! 556 556 !-- Topography height on w grid. Therefore, determine index of 557 557 !-- upward-facing surface element on w grid. 558 zw_w_inner(i,j) = zw( get_topography_top_index ( j, i, 's' ) )558 zw_w_inner(i,j) = zw( get_topography_top_index_ji( j, i, 's' ) ) 559 559 ENDDO 560 560 ENDDO … … 1471 1471 1472 1472 USE surface_mod, & 1473 ONLY: get_topography_top_index 1473 ONLY: get_topography_top_index, get_topography_top_index_ji 1474 1474 1475 1475 IMPLICIT NONE -
palm/trunk/SOURCE/lpm.f90
r2696 r2698 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Grid indices passed to lpm_boundary_conds. (responsible Philipp Thiele) 23 23 ! 24 24 ! Former revisions: … … 343 343 !-- optimization is still possible.) 344 344 IF ( topography /= 'flat' .AND. k < nzb_max + 2 ) THEN 345 CALL lpm_boundary_conds( 'walls')345 CALL lpm_boundary_conds( 'walls', i, j, k ) 346 346 ENDIF 347 347 ! … … 353 353 !-- the top or bottom boundary and delete those particles, which are 354 354 !-- older than allowed 355 CALL lpm_boundary_conds( 'bottom/top' )355 CALL lpm_boundary_conds( 'bottom/top', i, j, k ) 356 356 ! 357 357 !--- If not all particles of the actual grid cell have reached the -
palm/trunk/SOURCE/lpm_advec.f90
r2696 r2698 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Particle interpolations at walls in case of SGS velocities revised and not 23 ! required parts are removed. (responsible Philipp Thiele) 24 ! 25 ! Bugfix in get_topography_top_index 23 26 ! 24 27 ! Former revisions: … … 149 152 150 153 USE indices, & 151 ONLY: nzb, nzt 154 ONLY: nzb, nzt, wall_flags_0 152 155 153 156 USE kinds … … 163 166 164 167 USE surface_mod, & 165 ONLY: get_topography_top_index , surf_def_h, surf_lsm_h, surf_usm_h168 ONLY: get_topography_top_index_ji, surf_def_h, surf_lsm_h, surf_usm_h 166 169 167 170 IMPLICIT NONE 171 172 LOGICAL :: subbox_at_wall !< flag to see if the current subgridbox is adjacent to a wall 168 173 169 174 INTEGER(iwp) :: agp !< loop variable … … 302 307 !-- above topography (Prandtl-layer height) 303 308 !-- Determine vertical index of topography top 304 k_wall = get_topography_top_index ( jp, ip, 's' )309 k_wall = get_topography_top_index_ji( jp, ip, 's' ) 305 310 306 311 IF ( constant_flux_layer .AND. zv(n) - zw(k_wall) < z_p ) THEN … … 397 402 ! 398 403 !-- Determine vertical index of topography top 399 k_wall = get_topography_top_index ( jp,ip, 's' )404 k_wall = get_topography_top_index_ji( jp,ip, 's' ) 400 405 401 406 IF ( constant_flux_layer .AND. zv(n) - zw(k_wall) < z_p ) THEN … … 528 533 !-- velocities 529 534 IF ( use_sgs_for_particles .AND. .NOT. cloud_droplets ) THEN 530 531 IF ( topography == 'flat' ) THEN 532 533 DO nb = 0,7 534 535 536 DO nb = 0,7 537 538 subbox_at_wall = .FALSE. 539 ! 540 !-- In case of topography check if subbox is adjacent to a wall 541 IF ( .NOT. topography == 'flat' ) THEN 542 i = ip + MERGE( -1_iwp , 1_iwp, BTEST( nb, 2 ) ) 543 j = jp + MERGE( -1_iwp , 1_iwp, BTEST( nb, 1 ) ) 544 k = kp + MERGE( -1_iwp , 1_iwp, BTEST( nb, 0 ) ) 545 IF ( .NOT. BTEST(wall_flags_0(k, jp, ip), 0) .OR. & 546 .NOT. BTEST(wall_flags_0(kp, j, ip), 0) .OR. & 547 .NOT. BTEST(wall_flags_0(kp, jp, i ), 0) ) & 548 THEN 549 subbox_at_wall = .TRUE. 550 ENDIF 551 ENDIF 552 IF ( subbox_at_wall ) THEN 553 e_int(start_index(nb):end_index(nb)) = e(kp,jp,ip) 554 diss_int(start_index(nb):end_index(nb)) = diss(kp,jp,ip) 555 de_dx_int(start_index(nb):end_index(nb)) = de_dx(kp,jp,ip) 556 de_dy_int(start_index(nb):end_index(nb)) = de_dy(kp,jp,ip) 557 de_dz_int(start_index(nb):end_index(nb)) = de_dz(kp,jp,ip) 558 ! 559 !-- Set flag for stochastic equation. 560 term_1_2(start_index(nb):end_index(nb)) = 0.0_wp 561 ELSE 535 562 i = ip + block_offset(nb)%i_off 536 563 j = jp + block_offset(nb)%j_off … … 600 627 ELSE 601 628 de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i) + & 602 ( gg - bb ) * de_dy(k+1,j,i+1) + &603 ( gg - cc ) * de_dy(k+1,j+1,i) + &604 ( gg - dd ) * de_dy(k+1,j+1,i+1) &629 ( gg - bb ) * de_dy(k+1,j,i+1) + & 630 ( gg - cc ) * de_dy(k+1,j+1,i) + & 631 ( gg - dd ) * de_dy(k+1,j+1,i+1) & 605 632 ) / ( 3.0_wp * gg ) 606 de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / dz * &607 ( de_dy_int_u - de_dy_int_l )633 de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / dz * & 634 ( de_dy_int_u - de_dy_int_l ) 608 635 ENDIF 609 636 … … 638 665 ( gg - cc ) * diss(k,j+1,i) + & 639 666 ( gg - dd ) * diss(k,j+1,i+1) & 640 ) / ( 3.0_wp * gg )667 ) / ( 3.0_wp * gg ) 641 668 642 669 IF ( k == nzt ) THEN … … 649 676 ) / ( 3.0_wp * gg ) 650 677 diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dz * & 651 ( diss_int_u - diss_int_l )678 ( diss_int_u - diss_int_l ) 652 679 ENDIF 653 680 … … 655 682 !-- Set flag for stochastic equation. 656 683 term_1_2(n) = 1.0_wp 657 658 684 ENDDO 659 ENDDO 660 661 ELSE ! non-flat topography, e.g., buildings 662 663 DO nb = 0, 7 664 665 DO n = start_index(nb), end_index(nb) 666 667 i = particles(n)%x * ddx 668 j = particles(n)%y * ddy 669 k = ( zv(n) + dz * atmos_ocean_sign ) / dz & 670 + offset_ocean_nzt ! only exact if eq.dist 671 ! 672 !-- In case that there are buildings it has to be determined 673 !-- how many of the gridpoints defining the particle box are 674 !-- situated within a building 675 !-- gp_outside_of_building(1): i,j,k 676 !-- gp_outside_of_building(2): i,j+1,k 677 !-- gp_outside_of_building(3): i,j,k+1 678 !-- gp_outside_of_building(4): i,j+1,k+1 679 !-- gp_outside_of_building(5): i+1,j,k 680 !-- gp_outside_of_building(6): i+1,j+1,k 681 !-- gp_outside_of_building(7): i+1,j,k+1 682 !-- gp_outside_of_building(8): i+1,j+1,k+1 683 684 gp_outside_of_building = 0 685 location = 0.0_wp 686 num_gp = 0 687 688 ! 689 !-- Determine vertical index of topography top at (j,i) 690 k_wall = get_topography_top_index( j, i, 's' ) 691 ! 692 !-- To do: Reconsider order of computations in order to avoid 693 !-- unnecessary index calculations. 694 IF ( k > k_wall .OR. k_wall == 0 ) THEN 695 num_gp = num_gp + 1 696 gp_outside_of_building(1) = 1 697 location(num_gp,1) = i * dx 698 location(num_gp,2) = j * dy 699 location(num_gp,3) = k * dz - 0.5_wp * dz 700 ei(num_gp) = e(k,j,i) 701 dissi(num_gp) = diss(k,j,i) 702 de_dxi(num_gp) = de_dx(k,j,i) 703 de_dyi(num_gp) = de_dy(k,j,i) 704 de_dzi(num_gp) = de_dz(k,j,i) 705 ENDIF 706 707 ! 708 !-- Determine vertical index of topography top at (j+1,i) 709 k_wall = get_topography_top_index( j+1, i, 's' ) 710 711 IF ( k > k_wall .OR. k_wall == 0 ) THEN 712 num_gp = num_gp + 1 713 gp_outside_of_building(2) = 1 714 location(num_gp,1) = i * dx 715 location(num_gp,2) = (j+1) * dy 716 location(num_gp,3) = k * dz - 0.5_wp * dz 717 ei(num_gp) = e(k,j+1,i) 718 dissi(num_gp) = diss(k,j+1,i) 719 de_dxi(num_gp) = de_dx(k,j+1,i) 720 de_dyi(num_gp) = de_dy(k,j+1,i) 721 de_dzi(num_gp) = de_dz(k,j+1,i) 722 ENDIF 723 724 ! 725 !-- Determine vertical index of topography top at (j,i) 726 k_wall = get_topography_top_index( j, i, 's' ) 727 728 IF ( k+1 > k_wall .OR. k_wall == 0 ) THEN 729 num_gp = num_gp + 1 730 gp_outside_of_building(3) = 1 731 location(num_gp,1) = i * dx 732 location(num_gp,2) = j * dy 733 location(num_gp,3) = (k+1) * dz - 0.5_wp * dz 734 ei(num_gp) = e(k+1,j,i) 735 dissi(num_gp) = diss(k+1,j,i) 736 de_dxi(num_gp) = de_dx(k+1,j,i) 737 de_dyi(num_gp) = de_dy(k+1,j,i) 738 de_dzi(num_gp) = de_dz(k+1,j,i) 739 ENDIF 740 741 ! 742 !-- Determine vertical index of topography top at (j+1,i) 743 k_wall = get_topography_top_index( j+1, i, 's' ) 744 IF ( k+1 > k_wall .OR. k_wall == 0 ) THEN 745 num_gp = num_gp + 1 746 gp_outside_of_building(4) = 1 747 location(num_gp,1) = i * dx 748 location(num_gp,2) = (j+1) * dy 749 location(num_gp,3) = (k+1) * dz - 0.5_wp * dz 750 ei(num_gp) = e(k+1,j+1,i) 751 dissi(num_gp) = diss(k+1,j+1,i) 752 de_dxi(num_gp) = de_dx(k+1,j+1,i) 753 de_dyi(num_gp) = de_dy(k+1,j+1,i) 754 de_dzi(num_gp) = de_dz(k+1,j+1,i) 755 ENDIF 756 757 ! 758 !-- Determine vertical index of topography top at (j,i+1) 759 k_wall = get_topography_top_index( j, i+1, 's' ) 760 IF ( k > k_wall .OR. k_wall == 0 ) THEN 761 num_gp = num_gp + 1 762 gp_outside_of_building(5) = 1 763 location(num_gp,1) = (i+1) * dx 764 location(num_gp,2) = j * dy 765 location(num_gp,3) = k * dz - 0.5_wp * dz 766 ei(num_gp) = e(k,j,i+1) 767 dissi(num_gp) = diss(k,j,i+1) 768 de_dxi(num_gp) = de_dx(k,j,i+1) 769 de_dyi(num_gp) = de_dy(k,j,i+1) 770 de_dzi(num_gp) = de_dz(k,j,i+1) 771 ENDIF 772 773 ! 774 !-- Determine vertical index of topography top at (j+1,i+1) 775 k_wall = get_topography_top_index( j+1, i+1, 's' ) 776 777 IF ( k > k_wall .OR. k_wall == 0 ) THEN 778 num_gp = num_gp + 1 779 gp_outside_of_building(6) = 1 780 location(num_gp,1) = (i+1) * dx 781 location(num_gp,2) = (j+1) * dy 782 location(num_gp,3) = k * dz - 0.5_wp * dz 783 ei(num_gp) = e(k,j+1,i+1) 784 dissi(num_gp) = diss(k,j+1,i+1) 785 de_dxi(num_gp) = de_dx(k,j+1,i+1) 786 de_dyi(num_gp) = de_dy(k,j+1,i+1) 787 de_dzi(num_gp) = de_dz(k,j+1,i+1) 788 ENDIF 789 790 ! 791 !-- Determine vertical index of topography top at (j,i+1) 792 k_wall = get_topography_top_index( j, i+1, 's' ) 793 794 IF ( k+1 > k_wall .OR. k_wall == 0 ) THEN 795 num_gp = num_gp + 1 796 gp_outside_of_building(7) = 1 797 location(num_gp,1) = (i+1) * dx 798 location(num_gp,2) = j * dy 799 location(num_gp,3) = (k+1) * dz - 0.5_wp * dz 800 ei(num_gp) = e(k+1,j,i+1) 801 dissi(num_gp) = diss(k+1,j,i+1) 802 de_dxi(num_gp) = de_dx(k+1,j,i+1) 803 de_dyi(num_gp) = de_dy(k+1,j,i+1) 804 de_dzi(num_gp) = de_dz(k+1,j,i+1) 805 ENDIF 806 807 ! 808 !-- Determine vertical index of topography top at (j+1,i+1) 809 k_wall = get_topography_top_index( j+1, i+1, 's' ) 810 811 IF ( k+1 > k_wall .OR. k_wall == 0) THEN 812 num_gp = num_gp + 1 813 gp_outside_of_building(8) = 1 814 location(num_gp,1) = (i+1) * dx 815 location(num_gp,2) = (j+1) * dy 816 location(num_gp,3) = (k+1) * dz - 0.5_wp * dz 817 ei(num_gp) = e(k+1,j+1,i+1) 818 dissi(num_gp) = diss(k+1,j+1,i+1) 819 de_dxi(num_gp) = de_dx(k+1,j+1,i+1) 820 de_dyi(num_gp) = de_dy(k+1,j+1,i+1) 821 de_dzi(num_gp) = de_dz(k+1,j+1,i+1) 822 ENDIF 823 ! 824 !-- If all gridpoints are situated outside of a building, then the 825 !-- ordinary interpolation scheme can be used. 826 IF ( num_gp == 8 ) THEN 827 828 x = particles(n)%x - i * dx 829 y = particles(n)%y - j * dy 830 aa = x**2 + y**2 831 bb = ( dx - x )**2 + y**2 832 cc = x**2 + ( dy - y )**2 833 dd = ( dx - x )**2 + ( dy - y )**2 834 gg = aa + bb + cc + dd 835 836 e_int_l = ( ( gg - aa ) * e(k,j,i) + ( gg - bb ) * e(k,j,i+1) & 837 + ( gg - cc ) * e(k,j+1,i) + ( gg - dd ) * e(k,j+1,i+1) & 838 ) / ( 3.0_wp * gg ) 839 840 IF ( k == nzt ) THEN 841 e_int(n) = e_int_l 842 ELSE 843 e_int_u = ( ( gg - aa ) * e(k+1,j,i) + & 844 ( gg - bb ) * e(k+1,j,i+1) + & 845 ( gg - cc ) * e(k+1,j+1,i) + & 846 ( gg - dd ) * e(k+1,j+1,i+1) & 847 ) / ( 3.0_wp * gg ) 848 e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dz * & 849 ( e_int_u - e_int_l ) 850 ENDIF 851 ! 852 !-- Needed to avoid NaN particle velocities (this might not be 853 !-- required any more) 854 IF ( e_int(n) <= 0.0_wp ) THEN 855 e_int(n) = 1.0E-20_wp 856 ENDIF 857 ! 858 !-- Interpolate the TKE gradient along x (adopt incides i,j,k 859 !-- and all position variables from above (TKE)) 860 de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i) + & 861 ( gg - bb ) * de_dx(k,j,i+1) + & 862 ( gg - cc ) * de_dx(k,j+1,i) + & 863 ( gg - dd ) * de_dx(k,j+1,i+1) & 864 ) / ( 3.0_wp * gg ) 865 866 IF ( ( k == nzt ) .OR. ( k == nzb ) ) THEN 867 de_dx_int(n) = de_dx_int_l 868 ELSE 869 de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i) + & 870 ( gg - bb ) * de_dx(k+1,j,i+1) + & 871 ( gg - cc ) * de_dx(k+1,j+1,i) + & 872 ( gg - dd ) * de_dx(k+1,j+1,i+1) & 873 ) / ( 3.0_wp * gg ) 874 de_dx_int(n) = de_dx_int_l + ( zv(n) - zu(k) ) / & 875 dz * ( de_dx_int_u - de_dx_int_l ) 876 ENDIF 877 878 ! 879 !-- Interpolate the TKE gradient along y 880 de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i) + & 881 ( gg - bb ) * de_dy(k,j,i+1) + & 882 ( gg - cc ) * de_dy(k,j+1,i) + & 883 ( gg - dd ) * de_dy(k,j+1,i+1) & 884 ) / ( 3.0_wp * gg ) 885 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 886 de_dy_int(n) = de_dy_int_l 887 ELSE 888 de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i) + & 889 ( gg - bb ) * de_dy(k+1,j,i+1) + & 890 ( gg - cc ) * de_dy(k+1,j+1,i) + & 891 ( gg - dd ) * de_dy(k+1,j+1,i+1) & 892 ) / ( 3.0_wp * gg ) 893 de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / & 894 dz * ( de_dy_int_u - de_dy_int_l ) 895 ENDIF 896 897 ! 898 !-- Interpolate the TKE gradient along z 899 IF ( zv(n) < 0.5_wp * dz ) THEN 900 de_dz_int(n) = 0.0_wp 901 ELSE 902 de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i) + & 903 ( gg - bb ) * de_dz(k,j,i+1) + & 904 ( gg - cc ) * de_dz(k,j+1,i) + & 905 ( gg - dd ) * de_dz(k,j+1,i+1) & 906 ) / ( 3.0_wp * gg ) 907 908 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 909 de_dz_int(n) = de_dz_int_l 910 ELSE 911 de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i) + & 912 ( gg - bb ) * de_dz(k+1,j,i+1) + & 913 ( gg - cc ) * de_dz(k+1,j+1,i) + & 914 ( gg - dd ) * de_dz(k+1,j+1,i+1) & 915 ) / ( 3.0_wp * gg ) 916 de_dz_int(n) = de_dz_int_l + ( zv(n) - zu(k) ) /& 917 dz * ( de_dz_int_u - de_dz_int_l ) 918 ENDIF 919 ENDIF 920 921 ! 922 !-- Interpolate the dissipation of TKE 923 diss_int_l = ( ( gg - aa ) * diss(k,j,i) + & 924 ( gg - bb ) * diss(k,j,i+1) + & 925 ( gg - cc ) * diss(k,j+1,i) + & 926 ( gg - dd ) * diss(k,j+1,i+1) & 927 ) / ( 3.0_wp * gg ) 928 929 IF ( k == nzt ) THEN 930 diss_int(n) = diss_int_l 931 ELSE 932 diss_int_u = ( ( gg - aa ) * diss(k+1,j,i) + & 933 ( gg - bb ) * diss(k+1,j,i+1) + & 934 ( gg - cc ) * diss(k+1,j+1,i) + & 935 ( gg - dd ) * diss(k+1,j+1,i+1) & 936 ) / ( 3.0_wp * gg ) 937 diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dz *& 938 ( diss_int_u - diss_int_l ) 939 ENDIF 940 ! 941 !-- Set flag for stochastic equation. 942 term_1_2(n) = 1.0_wp 943 944 ELSE 945 946 ! 947 !-- If wall between gridpoint 1 and gridpoint 5, then 948 !-- Neumann boundary condition has to be applied 949 IF ( gp_outside_of_building(1) == 1 .AND. & 950 gp_outside_of_building(5) == 0 ) THEN 951 num_gp = num_gp + 1 952 location(num_gp,1) = i * dx + 0.5_wp * dx 953 location(num_gp,2) = j * dy 954 location(num_gp,3) = k * dz - 0.5_wp * dz 955 ei(num_gp) = e(k,j,i) 956 dissi(num_gp) = diss(k,j,i) 957 de_dxi(num_gp) = 0.0_wp 958 de_dyi(num_gp) = de_dy(k,j,i) 959 de_dzi(num_gp) = de_dz(k,j,i) 960 ENDIF 961 962 IF ( gp_outside_of_building(5) == 1 .AND. & 963 gp_outside_of_building(1) == 0 ) THEN 964 num_gp = num_gp + 1 965 location(num_gp,1) = i * dx + 0.5_wp * dx 966 location(num_gp,2) = j * dy 967 location(num_gp,3) = k * dz - 0.5_wp * dz 968 ei(num_gp) = e(k,j,i+1) 969 dissi(num_gp) = diss(k,j,i+1) 970 de_dxi(num_gp) = 0.0_wp 971 de_dyi(num_gp) = de_dy(k,j,i+1) 972 de_dzi(num_gp) = de_dz(k,j,i+1) 973 ENDIF 974 975 ! 976 !-- If wall between gridpoint 5 and gridpoint 6, then 977 !-- then Neumann boundary condition has to be applied 978 IF ( gp_outside_of_building(5) == 1 .AND. & 979 gp_outside_of_building(6) == 0 ) THEN 980 num_gp = num_gp + 1 981 location(num_gp,1) = (i+1) * dx 982 location(num_gp,2) = j * dy + 0.5_wp * dy 983 location(num_gp,3) = k * dz - 0.5_wp * dz 984 ei(num_gp) = e(k,j,i+1) 985 dissi(num_gp) = diss(k,j,i+1) 986 de_dxi(num_gp) = de_dx(k,j,i+1) 987 de_dyi(num_gp) = 0.0_wp 988 de_dzi(num_gp) = de_dz(k,j,i+1) 989 ENDIF 990 991 IF ( gp_outside_of_building(6) == 1 .AND. & 992 gp_outside_of_building(5) == 0 ) THEN 993 num_gp = num_gp + 1 994 location(num_gp,1) = (i+1) * dx 995 location(num_gp,2) = j * dy + 0.5_wp * dy 996 location(num_gp,3) = k * dz - 0.5_wp * dz 997 ei(num_gp) = e(k,j+1,i+1) 998 dissi(num_gp) = diss(k,j+1,i+1) 999 de_dxi(num_gp) = de_dx(k,j+1,i+1) 1000 de_dyi(num_gp) = 0.0_wp 1001 de_dzi(num_gp) = de_dz(k,j+1,i+1) 1002 ENDIF 1003 1004 ! 1005 !-- If wall between gridpoint 2 and gridpoint 6, then 1006 !-- Neumann boundary condition has to be applied 1007 IF ( gp_outside_of_building(2) == 1 .AND. & 1008 gp_outside_of_building(6) == 0 ) THEN 1009 num_gp = num_gp + 1 1010 location(num_gp,1) = i * dx + 0.5_wp * dx 1011 location(num_gp,2) = (j+1) * dy 1012 location(num_gp,3) = k * dz - 0.5_wp * dz 1013 ei(num_gp) = e(k,j+1,i) 1014 dissi(num_gp) = diss(k,j+1,i) 1015 de_dxi(num_gp) = 0.0_wp 1016 de_dyi(num_gp) = de_dy(k,j+1,i) 1017 de_dzi(num_gp) = de_dz(k,j+1,i) 1018 ENDIF 1019 1020 IF ( gp_outside_of_building(6) == 1 .AND. & 1021 gp_outside_of_building(2) == 0 ) THEN 1022 num_gp = num_gp + 1 1023 location(num_gp,1) = i * dx + 0.5_wp * dx 1024 location(num_gp,2) = (j+1) * dy 1025 location(num_gp,3) = k * dz - 0.5_wp * dz 1026 ei(num_gp) = e(k,j+1,i+1) 1027 dissi(num_gp) = diss(k,j+1,i+1) 1028 de_dxi(num_gp) = 0.0_wp 1029 de_dyi(num_gp) = de_dy(k,j+1,i+1) 1030 de_dzi(num_gp) = de_dz(k,j+1,i+1) 1031 ENDIF 1032 1033 ! 1034 !-- If wall between gridpoint 1 and gridpoint 2, then 1035 !-- Neumann boundary condition has to be applied 1036 IF ( gp_outside_of_building(1) == 1 .AND. & 1037 gp_outside_of_building(2) == 0 ) THEN 1038 num_gp = num_gp + 1 1039 location(num_gp,1) = i * dx 1040 location(num_gp,2) = j * dy + 0.5_wp * dy 1041 location(num_gp,3) = k * dz - 0.5_wp * dz 1042 ei(num_gp) = e(k,j,i) 1043 dissi(num_gp) = diss(k,j,i) 1044 de_dxi(num_gp) = de_dx(k,j,i) 1045 de_dyi(num_gp) = 0.0_wp 1046 de_dzi(num_gp) = de_dz(k,j,i) 1047 ENDIF 1048 1049 IF ( gp_outside_of_building(2) == 1 .AND. & 1050 gp_outside_of_building(1) == 0 ) THEN 1051 num_gp = num_gp + 1 1052 location(num_gp,1) = i * dx 1053 location(num_gp,2) = j * dy + 0.5_wp * dy 1054 location(num_gp,3) = k * dz - 0.5_wp * dz 1055 ei(num_gp) = e(k,j+1,i) 1056 dissi(num_gp) = diss(k,j+1,i) 1057 de_dxi(num_gp) = de_dx(k,j+1,i) 1058 de_dyi(num_gp) = 0.0_wp 1059 de_dzi(num_gp) = de_dz(k,j+1,i) 1060 ENDIF 1061 1062 ! 1063 !-- If wall between gridpoint 3 and gridpoint 7, then 1064 !-- Neumann boundary condition has to be applied 1065 IF ( gp_outside_of_building(3) == 1 .AND. & 1066 gp_outside_of_building(7) == 0 ) THEN 1067 num_gp = num_gp + 1 1068 location(num_gp,1) = i * dx + 0.5_wp * dx 1069 location(num_gp,2) = j * dy 1070 location(num_gp,3) = k * dz + 0.5_wp * dz 1071 ei(num_gp) = e(k+1,j,i) 1072 dissi(num_gp) = diss(k+1,j,i) 1073 de_dxi(num_gp) = 0.0_wp 1074 de_dyi(num_gp) = de_dy(k+1,j,i) 1075 de_dzi(num_gp) = de_dz(k+1,j,i) 1076 ENDIF 1077 1078 IF ( gp_outside_of_building(7) == 1 .AND. & 1079 gp_outside_of_building(3) == 0 ) THEN 1080 num_gp = num_gp + 1 1081 location(num_gp,1) = i * dx + 0.5_wp * dx 1082 location(num_gp,2) = j * dy 1083 location(num_gp,3) = k * dz + 0.5_wp * dz 1084 ei(num_gp) = e(k+1,j,i+1) 1085 dissi(num_gp) = diss(k+1,j,i+1) 1086 de_dxi(num_gp) = 0.0_wp 1087 de_dyi(num_gp) = de_dy(k+1,j,i+1) 1088 de_dzi(num_gp) = de_dz(k+1,j,i+1) 1089 ENDIF 1090 1091 ! 1092 !-- If wall between gridpoint 7 and gridpoint 8, then 1093 !-- Neumann boundary condition has to be applied 1094 IF ( gp_outside_of_building(7) == 1 .AND. & 1095 gp_outside_of_building(8) == 0 ) THEN 1096 num_gp = num_gp + 1 1097 location(num_gp,1) = (i+1) * dx 1098 location(num_gp,2) = j * dy + 0.5_wp * dy 1099 location(num_gp,3) = k * dz + 0.5_wp * dz 1100 ei(num_gp) = e(k+1,j,i+1) 1101 dissi(num_gp) = diss(k+1,j,i+1) 1102 de_dxi(num_gp) = de_dx(k+1,j,i+1) 1103 de_dyi(num_gp) = 0.0_wp 1104 de_dzi(num_gp) = de_dz(k+1,j,i+1) 1105 ENDIF 1106 1107 IF ( gp_outside_of_building(8) == 1 .AND. & 1108 gp_outside_of_building(7) == 0 ) THEN 1109 num_gp = num_gp + 1 1110 location(num_gp,1) = (i+1) * dx 1111 location(num_gp,2) = j * dy + 0.5_wp * dy 1112 location(num_gp,3) = k * dz + 0.5_wp * dz 1113 ei(num_gp) = e(k+1,j+1,i+1) 1114 dissi(num_gp) = diss(k+1,j+1,i+1) 1115 de_dxi(num_gp) = de_dx(k+1,j+1,i+1) 1116 de_dyi(num_gp) = 0.0_wp 1117 de_dzi(num_gp) = de_dz(k+1,j+1,i+1) 1118 ENDIF 1119 1120 ! 1121 !-- If wall between gridpoint 4 and gridpoint 8, then 1122 !-- Neumann boundary condition has to be applied 1123 IF ( gp_outside_of_building(4) == 1 .AND. & 1124 gp_outside_of_building(8) == 0 ) THEN 1125 num_gp = num_gp + 1 1126 location(num_gp,1) = i * dx + 0.5_wp * dx 1127 location(num_gp,2) = (j+1) * dy 1128 location(num_gp,3) = k * dz + 0.5_wp * dz 1129 ei(num_gp) = e(k+1,j+1,i) 1130 dissi(num_gp) = diss(k+1,j+1,i) 1131 de_dxi(num_gp) = 0.0_wp 1132 de_dyi(num_gp) = de_dy(k+1,j+1,i) 1133 de_dzi(num_gp) = de_dz(k+1,j+1,i) 1134 ENDIF 1135 1136 IF ( gp_outside_of_building(8) == 1 .AND. & 1137 gp_outside_of_building(4) == 0 ) THEN 1138 num_gp = num_gp + 1 1139 location(num_gp,1) = i * dx + 0.5_wp * dx 1140 location(num_gp,2) = (j+1) * dy 1141 location(num_gp,3) = k * dz + 0.5_wp * dz 1142 ei(num_gp) = e(k+1,j+1,i+1) 1143 dissi(num_gp) = diss(k+1,j+1,i+1) 1144 de_dxi(num_gp) = 0.0_wp 1145 de_dyi(num_gp) = de_dy(k+1,j+1,i+1) 1146 de_dzi(num_gp) = de_dz(k+1,j+1,i+1) 1147 ENDIF 1148 1149 ! 1150 !-- If wall between gridpoint 3 and gridpoint 4, then 1151 !-- Neumann boundary condition has to be applied 1152 IF ( gp_outside_of_building(3) == 1 .AND. & 1153 gp_outside_of_building(4) == 0 ) THEN 1154 num_gp = num_gp + 1 1155 location(num_gp,1) = i * dx 1156 location(num_gp,2) = j * dy + 0.5_wp * dy 1157 location(num_gp,3) = k * dz + 0.5_wp * dz 1158 ei(num_gp) = e(k+1,j,i) 1159 dissi(num_gp) = diss(k+1,j,i) 1160 de_dxi(num_gp) = de_dx(k+1,j,i) 1161 de_dyi(num_gp) = 0.0_wp 1162 de_dzi(num_gp) = de_dz(k+1,j,i) 1163 ENDIF 1164 1165 IF ( gp_outside_of_building(4) == 1 .AND. & 1166 gp_outside_of_building(3) == 0 ) THEN 1167 num_gp = num_gp + 1 1168 location(num_gp,1) = i * dx 1169 location(num_gp,2) = j * dy + 0.5_wp * dy 1170 location(num_gp,3) = k * dz + 0.5_wp * dz 1171 ei(num_gp) = e(k+1,j+1,i) 1172 dissi(num_gp) = diss(k+1,j+1,i) 1173 de_dxi(num_gp) = de_dx(k+1,j+1,i) 1174 de_dyi(num_gp) = 0.0_wp 1175 de_dzi(num_gp) = de_dz(k+1,j+1,i) 1176 ENDIF 1177 1178 ! 1179 !-- If wall between gridpoint 1 and gridpoint 3, then 1180 !-- Neumann boundary condition has to be applied 1181 !-- (only one case as only building beneath is possible) 1182 IF ( gp_outside_of_building(1) == 0 .AND. & 1183 gp_outside_of_building(3) == 1 ) THEN 1184 num_gp = num_gp + 1 1185 location(num_gp,1) = i * dx 1186 location(num_gp,2) = j * dy 1187 location(num_gp,3) = k * dz 1188 ei(num_gp) = e(k+1,j,i) 1189 dissi(num_gp) = diss(k+1,j,i) 1190 de_dxi(num_gp) = de_dx(k+1,j,i) 1191 de_dyi(num_gp) = de_dy(k+1,j,i) 1192 de_dzi(num_gp) = 0.0_wp 1193 ENDIF 1194 1195 ! 1196 !-- If wall between gridpoint 5 and gridpoint 7, then 1197 !-- Neumann boundary condition has to be applied 1198 !-- (only one case as only building beneath is possible) 1199 IF ( gp_outside_of_building(5) == 0 .AND. & 1200 gp_outside_of_building(7) == 1 ) THEN 1201 num_gp = num_gp + 1 1202 location(num_gp,1) = (i+1) * dx 1203 location(num_gp,2) = j * dy 1204 location(num_gp,3) = k * dz 1205 ei(num_gp) = e(k+1,j,i+1) 1206 dissi(num_gp) = diss(k+1,j,i+1) 1207 de_dxi(num_gp) = de_dx(k+1,j,i+1) 1208 de_dyi(num_gp) = de_dy(k+1,j,i+1) 1209 de_dzi(num_gp) = 0.0_wp 1210 ENDIF 1211 1212 ! 1213 !-- If wall between gridpoint 2 and gridpoint 4, then 1214 !-- Neumann boundary condition has to be applied 1215 !-- (only one case as only building beneath is possible) 1216 IF ( gp_outside_of_building(2) == 0 .AND. & 1217 gp_outside_of_building(4) == 1 ) THEN 1218 num_gp = num_gp + 1 1219 location(num_gp,1) = i * dx 1220 location(num_gp,2) = (j+1) * dy 1221 location(num_gp,3) = k * dz 1222 ei(num_gp) = e(k+1,j+1,i) 1223 dissi(num_gp) = diss(k+1,j+1,i) 1224 de_dxi(num_gp) = de_dx(k+1,j+1,i) 1225 de_dyi(num_gp) = de_dy(k+1,j+1,i) 1226 de_dzi(num_gp) = 0.0_wp 1227 ENDIF 1228 1229 ! 1230 !-- If wall between gridpoint 6 and gridpoint 8, then 1231 !-- Neumann boundary condition has to be applied 1232 !-- (only one case as only building beneath is possible) 1233 IF ( gp_outside_of_building(6) == 0 .AND. & 1234 gp_outside_of_building(8) == 1 ) THEN 1235 num_gp = num_gp + 1 1236 location(num_gp,1) = (i+1) * dx 1237 location(num_gp,2) = (j+1) * dy 1238 location(num_gp,3) = k * dz 1239 ei(num_gp) = e(k+1,j+1,i+1) 1240 dissi(num_gp) = diss(k+1,j+1,i+1) 1241 de_dxi(num_gp) = de_dx(k+1,j+1,i+1) 1242 de_dyi(num_gp) = de_dy(k+1,j+1,i+1) 1243 de_dzi(num_gp) = 0.0_wp 1244 ENDIF 1245 1246 ! 1247 !-- Carry out the interpolation 1248 IF ( num_gp == 1 ) THEN 1249 ! 1250 !-- If only one of the gridpoints is situated outside of the 1251 !-- building, it follows that the values at the particle 1252 !-- location are the same as the gridpoint values 1253 e_int(n) = ei(num_gp) 1254 diss_int(n) = dissi(num_gp) 1255 de_dx_int(n) = de_dxi(num_gp) 1256 de_dy_int(n) = de_dyi(num_gp) 1257 de_dz_int(n) = de_dzi(num_gp) 1258 ! 1259 !-- Set flag for stochastic equation which disables calculation 1260 !-- of drift and memory term near topography. 1261 term_1_2(n) = 0.0_wp 1262 ELSE IF ( num_gp > 1 ) THEN 1263 1264 d_sum = 0.0_wp 1265 ! 1266 !-- Evaluation of the distances between the gridpoints 1267 !-- contributing to the interpolated values, and the particle 1268 !-- location 1269 DO agp = 1, num_gp 1270 d_gp_pl(agp) = ( particles(n)%x-location(agp,1) )**2 & 1271 + ( particles(n)%y-location(agp,2) )**2 & 1272 + ( zv(n)-location(agp,3) )**2 1273 d_sum = d_sum + d_gp_pl(agp) 1274 ENDDO 1275 1276 ! 1277 !-- Finally the interpolation can be carried out 1278 e_int(n) = 0.0_wp 1279 diss_int(n) = 0.0_wp 1280 de_dx_int(n) = 0.0_wp 1281 de_dy_int(n) = 0.0_wp 1282 de_dz_int(n) = 0.0_wp 1283 DO agp = 1, num_gp 1284 e_int(n) = e_int(n) + ( d_sum - d_gp_pl(agp) ) * & 1285 ei(agp) / ( (num_gp-1) * d_sum ) 1286 diss_int(n) = diss_int(n) + ( d_sum - d_gp_pl(agp) ) * & 1287 dissi(agp) / ( (num_gp-1) * d_sum ) 1288 de_dx_int(n) = de_dx_int(n) + ( d_sum - d_gp_pl(agp) ) * & 1289 de_dxi(agp) / ( (num_gp-1) * d_sum ) 1290 de_dy_int(n) = de_dy_int(n) + ( d_sum - d_gp_pl(agp) ) * & 1291 de_dyi(agp) / ( (num_gp-1) * d_sum ) 1292 de_dz_int(n) = de_dz_int(n) + ( d_sum - d_gp_pl(agp) ) * & 1293 de_dzi(agp) / ( (num_gp-1) * d_sum ) 1294 ENDDO 1295 1296 ENDIF 1297 e_int(n) = MAX( 1E-20_wp, e_int(n) ) 1298 diss_int(n) = MAX( 1E-20_wp, diss_int(n) ) 1299 de_dx_int(n) = MAX( 1E-20_wp, de_dx_int(n) ) 1300 de_dy_int(n) = MAX( 1E-20_wp, de_dy_int(n) ) 1301 de_dz_int(n) = MAX( 1E-20_wp, de_dz_int(n) ) 1302 ! 1303 !-- Set flag for stochastic equation which disables calculation 1304 !-- of drift and memory term near topography. 1305 term_1_2(n) = 0.0_wp 1306 ENDIF 1307 ENDDO 1308 ENDDO 1309 ENDIF 685 ENDIF 686 ENDDO 1310 687 1311 688 DO nb = 0,7 -
palm/trunk/SOURCE/lpm_boundary_conds.f90
r2696 r2698 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Particle reflections at downward-facing walls implemented. Moreover, 23 ! reflections are adjusted to revised particle grid box location. 24 ! (responsible Philipp Thiele) 23 25 ! 24 26 ! Former revisions: … … 92 94 !> (see offset_ocean_*) 93 95 !------------------------------------------------------------------------------! 94 SUBROUTINE lpm_boundary_conds( location )96 SUBROUTINE lpm_boundary_conds( location , i, j, k ) 95 97 96 98 … … 108 110 109 111 USE indices, & 110 ONLY: nxl, nxr, nyn, nys, nz, nzb 112 ONLY: nxl, nxr, nyn, nys, nz, nzb, wall_flags_0,nyng,nysg 111 113 112 114 USE kinds … … 119 121 USE pegrid 120 122 121 USE surface_mod, &122 ONLY: get_topography_top_index123 124 123 IMPLICIT NONE 125 124 126 125 CHARACTER (LEN=*) :: location !< 126 127 INTEGER(iwp), INTENT(IN) :: i !< 128 INTEGER(iwp), INTENT(IN) :: j !< 129 INTEGER(iwp), INTENT(IN) :: k !< 127 130 128 131 INTEGER(iwp) :: inc !< dummy for sorting algorithmus … … 133 136 INTEGER(iwp) :: jr !< dummy for sorting algorithmus 134 137 INTEGER(iwp) :: j1 !< grid index (y) of old particle position 135 INTEGER(iwp) :: j2 !< grid index (x) of current particle position 136 INTEGER(iwp) :: j3 !< grid index (x) of intermediate particle position 137 INTEGER(iwp) :: k_wall !< vertical index of topography top 138 INTEGER(iwp) :: j2 !< grid index (y) of current particle position 139 INTEGER(iwp) :: j3 !< grid index (y) of intermediate particle position 140 INTEGER(iwp) :: k1 !< grid index (z) of old particle position 141 INTEGER(iwp) :: k2 !< grid index (z) of current particle position 142 INTEGER(iwp) :: k3 !< grid index (z) of intermediate particle position 138 143 INTEGER(iwp) :: n !< particle number 139 144 INTEGER(iwp) :: t_index !< running index for intermediate particle timesteps in reflection algorithmus 140 145 INTEGER(iwp) :: t_index_number !< number of intermediate particle timesteps in reflection algorithmus 141 INTEGER(iwp) :: tmp_x !< dummy for sorting algorithmus 142 INTEGER(iwp) :: tmp_y !< dummy for sorting algorithmus 146 INTEGER(iwp) :: tmp_x !< dummy for sorting algorithm 147 INTEGER(iwp) :: tmp_y !< dummy for sorting algorithm 148 INTEGER(iwp) :: tmp_z !< dummy for sorting algorithm 143 149 144 150 INTEGER(iwp), DIMENSION(0:10) :: x_ind(0:10) = 0 !< index array (x) of intermediate particle positions 145 INTEGER(iwp), DIMENSION(0:10) :: y_ind(0:10) = 0 !< index array (x) of intermediate particle positions 151 INTEGER(iwp), DIMENSION(0:10) :: y_ind(0:10) = 0 !< index array (y) of intermediate particle positions 152 INTEGER(iwp), DIMENSION(0:10) :: z_ind(0:10) = 0 !< index array (z) of intermediate particle positions 146 153 147 154 LOGICAL :: cross_wall_x !< flag to check if particle reflection along x is necessary 148 155 LOGICAL :: cross_wall_y !< flag to check if particle reflection along y is necessary 149 LOGICAL :: downwards !< flag to check if particle reflection along z is necessary (only if particle move downwards)156 LOGICAL :: cross_wall_z !< flag to check if particle reflection along z is necessary 150 157 LOGICAL :: reflect_x !< flag to check if particle is already reflected along x 151 158 LOGICAL :: reflect_y !< flag to check if particle is already reflected along y … … 156 163 LOGICAL :: x_wall_reached !< flag to check if particle has already reached wall 157 164 LOGICAL :: y_wall_reached !< flag to check if particle has already reached wall 165 LOGICAL :: z_wall_reached !< flag to check if particle has already reached wall 158 166 159 167 LOGICAL, DIMENSION(0:10) :: reach_x !< flag to check if particle is at a yz-wall … … 177 185 REAL(wp) :: xwall !< location of wall in x 178 186 REAL(wp) :: ywall !< location of wall in y 179 REAL(wp) :: zwall1 !< location of wall in z (old grid box) 180 REAL(wp) :: zwall2 !< location of wall in z (current grid box) 181 REAL(wp) :: zwall3 !< location of wall in z (old y, current x) 182 REAL(wp) :: zwall4 !< location of wall in z (current y, old x) 187 REAL(wp) :: zwall !< location of wall in z 183 188 184 189 REAL(wp), DIMENSION(0:10) :: t !< reflection time … … 263 268 i2 = particles(n)%x * ddx 264 269 j2 = particles(n)%y * ddy 270 IF (zw(k) < particles(n)%z ) k2 = k + 1 271 IF (zw(k-1) > particles(n)%z ) k2 = k - 1 265 272 ! 266 273 !-- Save current particle positions … … 275 282 ! 276 283 !-- Obtain x/y indices for old particle positions 277 i1 = pos_x_old * ddx 278 j1 = pos_y_old * ddy 284 i1 = i 285 j1 = j 286 k1 = k 279 287 ! 280 288 !-- Determine horizontal as well as vertical walls at which particle can … … 283 291 !-- Wall to the right 284 292 IF ( prt_x > pos_x_old ) THEN 285 xwall = ( i1 + 0.5_wp) * dx293 xwall = ( i1 + 1 ) * dx 286 294 ! 287 295 !-- Wall to the left 288 296 ELSE 289 xwall = ( i1 - 0.5_wp )* dx297 xwall = i1 * dx 290 298 ENDIF 291 299 ! … … 293 301 !-- Wall to the north 294 302 IF ( prt_y > pos_y_old ) THEN 295 ywall = ( j1 + 0.5_wp) * dy303 ywall = ( j1 +1 ) * dy 296 304 !-- Wall to the south 297 305 ELSE 298 ywall = ( j1 - 0.5_wp ) * dy 299 ENDIF 300 ! 301 !-- Walls aligned in xy layer at which particle can be possiblly reflected. 302 !-- The construct of MERGE and BTEST is used to determine the topography- 303 !-- top index (former nzb_s_inner). 304 zwall1 = zw( get_topography_top_index( j2, i2, 's' ) ) 305 zwall2 = zw( get_topography_top_index( j1, i1, 's' ) ) 306 zwall3 = zw( get_topography_top_index( j1, i2, 's' ) ) 307 zwall4 = zw( get_topography_top_index( j2, i1, 's' ) ) 306 ywall = j1 * dy 307 ENDIF 308 309 IF ( prt_z > pos_z_old ) THEN 310 zwall = zw(k) 311 ELSE 312 zwall = zw(k-1) 313 ENDIF 308 314 ! 309 315 !-- Initialize flags to check if particle reflection is necessary 310 downwards = .FALSE.311 316 cross_wall_x = .FALSE. 312 317 cross_wall_y = .FALSE. 318 cross_wall_z = .FALSE. 313 319 ! 314 320 !-- Initialize flags to check if a wall is reached … … 318 324 ! 319 325 !-- Initialize flags to check if a particle was already reflected 320 reflect_x = .FALSE.321 reflect_y = .FALSE.322 reflect_z = .FALSE.323 ! 324 !-- Initialize flags to check if a verticalwall is already crossed.326 reflect_x = .FALSE. 327 reflect_y = .FALSE. 328 reflect_z = .FALSE. 329 ! 330 !-- Initialize flags to check if a wall is already crossed. 325 331 !-- ( Required to obtain correct indices. ) 326 332 x_wall_reached = .FALSE. 327 333 y_wall_reached = .FALSE. 334 z_wall_reached = .FALSE. 328 335 ! 329 336 !-- Initialize time array … … 342 349 x_ind(t_index) = i2 343 350 y_ind(t_index) = j1 351 z_ind(t_index) = k1 344 352 reach_x(t_index) = .TRUE. 345 353 reach_y(t_index) = .FALSE. … … 360 368 x_ind(t_index) = i1 361 369 y_ind(t_index) = j2 370 z_ind(t_index) = k1 362 371 reach_x(t_index) = .FALSE. 363 372 reach_y(t_index) = .TRUE. … … 369 378 ! 370 379 !-- z-direction 371 !-- At first, check if particle moves downwards. Only in this case a 372 !-- particle can be reflected vertically. 373 IF ( prt_z < pos_z_old ) THEN 374 375 downwards = .TRUE. 376 dum = 1.0_wp / MERGE( MAX( prt_z - pos_z_old, 1E-30_wp ), & 377 MIN( prt_z - pos_z_old, -1E-30_wp ), & 378 prt_z > pos_z_old ) 379 380 t(t_index) = ( zwall1 - pos_z_old ) * dum 381 x_ind(t_index) = i2 382 y_ind(t_index) = j2 383 reach_x(t_index) = .FALSE. 384 reach_y(t_index) = .FALSE. 385 reach_z(t_index) = .TRUE. 386 IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) & 387 t_index = t_index + 1 388 389 reach_x(t_index) = .FALSE. 390 reach_y(t_index) = .FALSE. 391 reach_z(t_index) = .TRUE. 392 t(t_index) = ( zwall2 - pos_z_old ) * dum 393 x_ind(t_index) = i1 394 y_ind(t_index) = j1 395 IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) & 396 t_index = t_index + 1 397 398 reach_x(t_index) = .FALSE. 399 reach_y(t_index) = .FALSE. 400 reach_z(t_index) = .TRUE. 401 t(t_index) = ( zwall3 - pos_z_old ) * dum 402 x_ind(t_index) = i2 403 y_ind(t_index) = j1 404 IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) & 405 t_index = t_index + 1 406 407 reach_x(t_index) = .FALSE. 408 reach_y(t_index) = .FALSE. 409 reach_z(t_index) = .TRUE. 410 t(t_index) = ( zwall4 - pos_z_old ) * dum 411 x_ind(t_index) = i1 412 y_ind(t_index) = j2 413 IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) & 414 t_index = t_index + 1 415 416 END IF 380 t(t_index) = (zwall - pos_z_old ) & 381 / MERGE( MAX( prt_z - pos_z_old, 1E-30_wp ), & 382 MIN( prt_z - pos_z_old, -1E-30_wp ), & 383 prt_z > pos_z_old ) 384 385 x_ind(t_index) = i1 386 y_ind(t_index) = j1 387 z_ind(t_index) = k2 388 reach_x(t_index) = .FALSE. 389 reach_y(t_index) = .FALSE. 390 reach_z(t_index) = .TRUE. 391 IF( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp) THEN 392 t_index = t_index + 1 393 cross_wall_z = .TRUE. 394 ENDIF 395 417 396 t_index_number = t_index - 1 418 397 ! 419 398 !-- Carry out reflection only if particle reaches any wall 420 IF ( cross_wall_x .OR. cross_wall_y .OR. downwards) THEN399 IF ( cross_wall_x .OR. cross_wall_y .OR. cross_wall_z ) THEN 421 400 ! 422 401 !-- Sort fractional timesteps in ascending order. Also sort the … … 435 414 tmp_x = x_ind(ir) 436 415 tmp_y = y_ind(ir) 416 tmp_z = z_ind(ir) 437 417 tmp_reach_x = reach_x(ir) 438 418 tmp_reach_y = reach_y(ir) … … 443 423 x_ind(jr) = x_ind(jr-inc) 444 424 y_ind(jr) = y_ind(jr-inc) 425 z_ind(jr) = z_ind(jr-inc) 445 426 reach_x(jr) = reach_x(jr-inc) 446 427 reach_y(jr) = reach_y(jr-inc) … … 452 433 x_ind(jr) = tmp_x 453 434 y_ind(jr) = tmp_y 435 z_ind(jr) = tmp_z 454 436 reach_x(jr) = tmp_reach_x 455 437 reach_y(jr) = tmp_reach_y … … 480 462 i3 = x_ind(t_index) 481 463 j3 = y_ind(t_index) 464 k3 = z_ind(t_index) 482 465 ! 483 466 !-- Check which wall is already reached 484 467 IF ( .NOT. x_wall_reached ) x_wall_reached = reach_x(t_index) 485 IF ( .NOT. y_wall_reached ) y_wall_reached = reach_y(t_index) 468 IF ( .NOT. y_wall_reached ) y_wall_reached = reach_y(t_index) 469 IF ( .NOT. z_wall_reached ) z_wall_reached = reach_z(t_index) 486 470 ! 487 471 !-- Check if a particle needs to be reflected at any yz-wall. If 488 472 !-- necessary, carry out reflection. Please note, a security 489 !-- constant is required, as the particle position do not473 !-- constant is required, as the particle position does not 490 474 !-- necessarily exactly match the wall location due to rounding 491 !-- errors. At first, determine index of topography top at (j3,i3) 492 k_wall = get_topography_top_index( j3, i3, 's' ) 493 IF ( ABS( pos_x - xwall ) < eps .AND. & 494 pos_z <= zw(k_wall) .AND. & 495 reach_x(t_index) .AND. & 475 !-- errors. 476 IF ( reach_x(t_index) .AND. & 477 ABS( pos_x - xwall ) < eps .AND. & 478 .NOT. BTEST(wall_flags_0(k3,j3,i3),0) .AND. & 496 479 .NOT. reflect_x ) THEN 497 ! 480 ! 481 ! 498 482 !-- Reflection in x-direction. 499 483 !-- Ensure correct reflection by MIN/MAX functions, depending on 500 484 !-- direction of particle transport. 501 !-- Due to rounding errors pos_x do not exactly matchesthe wall485 !-- Due to rounding errors pos_x does not exactly match the wall 502 486 !-- location, leading to erroneous reflection. 503 487 pos_x = MERGE( MIN( 2.0_wp * xwall - pos_x, xwall ), & … … 508 492 particles(n)%speed_x = - particles(n)%speed_x 509 493 ! 510 !-- Change alsosign of subgrid-scale particle speed494 !-- Also change sign of subgrid-scale particle speed 511 495 particles(n)%rvar1 = - particles(n)%rvar1 512 496 ! … … 514 498 reflect_x = .TRUE. 515 499 ! 516 !-- As particle do not crosses any further yz-wall during500 !-- As the particle does not cross any further yz-wall during 517 501 !-- this timestep, set further x-indices to the current one. 518 502 x_ind(t_index:t_index_number) = i1 … … 522 506 ELSEIF ( x_wall_reached .AND. .NOT. reflect_x ) THEN 523 507 x_ind(t_index:t_index_number) = i2 524 ENDIF 508 ENDIF !particle reflection in x direction done 509 525 510 ! 526 511 !-- Check if a particle needs to be reflected at any xz-wall. If 527 !-- necessary, carry out reflection. At first, determine index of 528 !-- topography top at (j3,i3) 529 k_wall = get_topography_top_index( j3, i3, 's' ) 530 IF ( ABS( pos_y - ywall ) < eps .AND. & 531 pos_z <= zw(k_wall) .AND. & 532 reach_y(t_index) .AND. & 533 .NOT. reflect_y ) THEN 534 512 !-- necessary, carry out reflection. Please note, a security 513 !-- constant is required, as the particle position does not 514 !-- necessarily exactly match the wall location due to rounding 515 !-- errors. 516 WRITE(9,*) 'wall_test_y nysg:',nysg ,' y: ',particles(n)%y,' j3: ',j3,'nyng:', nyng 517 FLUSH(9) 518 WRITE(9,*) BTEST(wall_flags_0(k3,j3,i3),0) 519 FLUSH(9) 520 IF ( reach_y(t_index) .AND. & 521 ABS( pos_y - ywall ) < eps .AND. & 522 .NOT. BTEST(wall_flags_0(k3,j3,i3),0) .AND. & 523 .NOT. reflect_y ) THEN 524 ! 525 ! 526 !-- Reflection in y-direction. 527 !-- Ensure correct reflection by MIN/MAX functions, depending on 528 !-- direction of particle transport. 529 !-- Due to rounding errors pos_y does not exactly match the wall 530 !-- location, leading to erroneous reflection. 535 531 pos_y = MERGE( MIN( 2.0_wp * ywall - pos_y, ywall ), & 536 532 MAX( 2.0_wp * ywall - pos_y, ywall ), & 537 particles(n)%y > ywall ) 538 539 particles(n)%speed_y = - particles(n)%speed_y 540 particles(n)%rvar2 = - particles(n)%rvar2 541 542 reflect_y = .TRUE. 533 particles(n)%y > ywall ) 534 ! 535 !-- Change sign of particle speed 536 particles(n)%speed_y = - particles(n)%speed_y 537 ! 538 !-- Also change sign of subgrid-scale particle speed 539 particles(n)%rvar2 = - particles(n)%rvar2 540 ! 541 !-- Set flag that reflection along y is already done 542 reflect_y = .TRUE. 543 ! 544 !-- As the particle does not cross any further xz-wall during 545 !-- this timestep, set further y-indices to the current one. 543 546 y_ind(t_index:t_index_number) = j1 544 547 ! 548 !-- If particle already reached the wall but was not reflected, 549 !-- set further y-indices to the new one. 545 550 ELSEIF ( y_wall_reached .AND. .NOT. reflect_y ) THEN 546 y_ind(t_index:t_index_number) = j2 547 ENDIF 551 y_ind(t_index:t_index_number) = j2 552 ENDIF !particle reflection in y direction done 553 548 554 ! 549 555 !-- Check if a particle needs to be reflected at any xy-wall. If 550 !-- necessary, carry out reflection. 551 IF ( downwards .AND. reach_z(t_index) .AND. & 556 !-- necessary, carry out reflection. Please note, a security 557 !-- constant is required, as the particle position does not 558 !-- necessarily exactly match the wall location due to rounding 559 !-- errors. 560 IF ( reach_z(t_index) .AND. & 561 ABS( pos_z - zwall ) < eps .AND. & 562 .NOT. BTEST(wall_flags_0(k3,j3,i3),0) .AND. & 552 563 .NOT. reflect_z ) THEN 553 ! 554 !-- Determine index of topography top at (j3,i3) and chick if 555 !-- particle is below. 556 k_wall = get_topography_top_index( j3, i3, 's' ) 557 IF ( pos_z - zw(k_wall) < eps ) THEN 558 559 pos_z = MAX( 2.0_wp * zw(k_wall) - pos_z, & 560 zw(k_wall) ) 561 562 particles(n)%speed_z = - particles(n)%speed_z 563 particles(n)%rvar3 = - particles(n)%rvar3 564 565 reflect_z = .TRUE. 566 567 ENDIF 568 569 ENDIF 564 ! 565 ! 566 !-- Reflection in z-direction. 567 !-- Ensure correct reflection by MIN/MAX functions, depending on 568 !-- direction of particle transport. 569 !-- Due to rounding errors pos_z does not exactly match the wall 570 !-- location, leading to erroneous reflection. 571 pos_z = MERGE( MIN( 2.0_wp * zwall - pos_z, zwall ), & 572 MAX( 2.0_wp * zwall - pos_z, zwall ), & 573 particles(n)%z > zwall ) 574 ! 575 !-- Change sign of particle speed 576 particles(n)%speed_z = - particles(n)%speed_z 577 ! 578 !-- Also change sign of subgrid-scale particle speed 579 particles(n)%rvar3 = - particles(n)%rvar3 580 ! 581 !-- Set flag that reflection along z is already done 582 reflect_z = .TRUE. 583 ! 584 !-- As the particle does not cross any further xy-wall during 585 !-- this timestep, set further z-indices to the current one. 586 z_ind(t_index:t_index_number) = k1 587 ! 588 !-- If particle already reached the wall but was not reflected, 589 !-- set further z-indices to the new one. 590 ELSEIF ( z_wall_reached .AND. .NOT. reflect_z ) THEN 591 z_ind(t_index:t_index_number) = k2 592 ENDIF !particle reflection in z direction done 593 570 594 ! 571 595 !-- Swap time -
palm/trunk/SOURCE/lpm_init.f90
r2696 r2698 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Grid indices passed to lpm_boundary_conds. (responsible Philipp Thiele) 23 23 ! 24 24 ! Former revisions: … … 236 236 237 237 USE surface_mod, & 238 ONLY: get_topography_top_index , surf_def_h, surf_lsm_h, surf_usm_h238 ONLY: get_topography_top_index_ji, surf_def_h, surf_lsm_h, surf_usm_h 239 239 240 240 IMPLICIT NONE … … 757 757 !-- Determine surface level. Therefore, check for 758 758 !-- upward-facing wall on w-grid. 759 k_surf = get_topography_top_index ( jp, ip, 'w' )759 k_surf = get_topography_top_index_ji( jp, ip, 'w' ) 760 760 761 761 IF ( seed_follows_topography ) THEN … … 930 930 !-- Identify particles located outside the model domain and reflect 931 931 !-- or absorb them if necessary. 932 CALL lpm_boundary_conds( 'bottom/top' )932 CALL lpm_boundary_conds( 'bottom/top', i, j, k ) 933 933 ! 934 934 !-- Furthermore, remove particles located in topography. Note, as -
palm/trunk/SOURCE/microphysics_mod.f90
r2696 r2698 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Bugfix in get_topography_top_index 23 23 ! 24 24 ! Former revisions: … … 928 928 929 929 USE surface_mod, & 930 ONLY: get_topography_top_index 930 ONLY: get_topography_top_index_ji 931 931 932 932 … … 945 945 ! 946 946 !-- Determine vertical index of topography top 947 k_wall = get_topography_top_index ( j, i, 's' )947 k_wall = get_topography_top_index_ji( j, i, 's' ) 948 948 DO k = nzb+1, nzt 949 949 ! … … 2319 2319 2320 2320 USE surface_mod, & 2321 ONLY: get_topography_top_index 2321 ONLY: get_topography_top_index_ji 2322 2322 2323 2323 … … 2334 2334 ! 2335 2335 !-- Determine vertical index of topography top 2336 k_wall = get_topography_top_index ( j, i, 's' )2336 k_wall = get_topography_top_index_ji( j, i, 's' ) 2337 2337 DO k = nzb+1, nzt 2338 2338 ! -
palm/trunk/SOURCE/plant_canopy_model_mod.f90
r2696 r2698 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Bugfix in get_topography_top_index 23 23 ! 24 24 ! Former revisions: … … 151 151 152 152 USE surface_mod, & 153 ONLY: get_topography_top_index 153 ONLY: get_topography_top_index_ji 154 154 155 155 … … 400 400 DO j = nys, nyn 401 401 IF ( pch_index_ji(j,i) /= 0 ) THEN 402 k_topo = get_topography_top_index ( j, i, 's' )402 k_topo = get_topography_top_index_ji( j, i, 's' ) 403 403 DO k = k_topo, k_topo + pch_index_ji(j,i) 404 404 local_pf(i,j,k) = pc_heating_rate(k-k_topo,j,i) … … 415 415 DO j = nys, nyn 416 416 IF ( pch_index_ji(j,i) /= 0 ) THEN 417 k_topo = get_topography_top_index ( j, i, 's' )417 k_topo = get_topography_top_index_ji( j, i, 's' ) 418 418 DO k = k_topo, k_topo + pch_index_ji(j,i) 419 419 local_pf(i,j,k) = lad_s(k-k_topo,j,i) … … 803 803 !-- Check whether topography and local vegetation on top exceed 804 804 !-- height of the model domain. 805 k = get_topography_top_index ( j, i, 's' )805 k = get_topography_top_index_ji( j, i, 's' ) 806 806 IF ( k + pch_index_ji(j,i) >= nzt + 1 ) THEN 807 807 message_string = 'Local vegetation height on top of ' // & … … 1110 1110 ! 1111 1111 !-- Determine topography-top index on u-grid 1112 k_wall = get_topography_top_index ( j, i, 'u' )1112 k_wall = get_topography_top_index_ji( j, i, 'u' ) 1113 1113 DO k = k_wall+1, k_wall + pch_index_ji(j,i) 1114 1114 … … 1174 1174 ! 1175 1175 !-- Determine topography-top index on v-grid 1176 k_wall = get_topography_top_index ( j, i, 'v' )1176 k_wall = get_topography_top_index_ji( j, i, 'v' ) 1177 1177 1178 1178 DO k = k_wall+1, k_wall + pch_index_ji(j,i) … … 1239 1239 ! 1240 1240 !-- Determine topography-top index on w-grid 1241 k_wall = get_topography_top_index ( j, i, 'w' )1241 k_wall = get_topography_top_index_ji( j, i, 'w' ) 1242 1242 1243 1243 DO k = k_wall+1, k_wall + pch_index_ji(j,i) - 1 … … 1291 1291 ! 1292 1292 !-- Determine topography-top index on scalar-grid 1293 k_wall = get_topography_top_index ( j, i, 's' )1293 k_wall = get_topography_top_index_ji( j, i, 's' ) 1294 1294 1295 1295 DO k = k_wall+1, k_wall + pch_index_ji(j,i) … … 1308 1308 ! 1309 1309 !-- Determine topography-top index on scalar-grid 1310 k_wall = get_topography_top_index ( j, i, 's' )1310 k_wall = get_topography_top_index_ji( j, i, 's' ) 1311 1311 1312 1312 DO k = k_wall+1, k_wall + pch_index_ji(j,i) … … 1338 1338 ! 1339 1339 !-- Determine topography-top index on scalar-grid 1340 k_wall = get_topography_top_index ( j, i, 's' )1340 k_wall = get_topography_top_index_ji( j, i, 's' ) 1341 1341 1342 1342 DO k = k_wall+1, k_wall + pch_index_ji(j,i) … … 1367 1367 ! 1368 1368 !-- Determine topography-top index on scalar-grid 1369 k_wall = get_topography_top_index ( j, i, 's' )1369 k_wall = get_topography_top_index_ji( j, i, 's' ) 1370 1370 1371 1371 DO k = k_wall+1, k_wall + pch_index_ji(j,i) … … 1461 1461 ! 1462 1462 !-- Determine topography-top index on u-grid 1463 k_wall = get_topography_top_index ( j, i, 'u' )1463 k_wall = get_topography_top_index_ji( j, i, 'u' ) 1464 1464 DO k = k_wall + 1, k_wall + pch_index_ji(j,i) 1465 1465 … … 1521 1521 ! 1522 1522 !-- Determine topography-top index on v-grid 1523 k_wall = get_topography_top_index ( j, i, 'v' )1523 k_wall = get_topography_top_index_ji( j, i, 'v' ) 1524 1524 1525 1525 DO k = k_wall + 1, k_wall + pch_index_ji(j,i) … … 1581 1581 ! 1582 1582 !-- Determine topography-top index on w-grid 1583 k_wall = get_topography_top_index ( j, i, 'w' )1583 k_wall = get_topography_top_index_ji( j, i, 'w' ) 1584 1584 1585 1585 DO k = k_wall + 1, k_wall + pch_index_ji(j,i) - 1 … … 1628 1628 ! 1629 1629 !-- Determine topography-top index on scalar grid 1630 k_wall = get_topography_top_index ( j, i, 's' )1630 k_wall = get_topography_top_index_ji( j, i, 's' ) 1631 1631 1632 1632 DO k = k_wall + 1, k_wall + pch_index_ji(j,i) … … 1641 1641 ! 1642 1642 !-- Determine topography-top index on scalar grid 1643 k_wall = get_topography_top_index ( j, i, 's' )1643 k_wall = get_topography_top_index_ji( j, i, 's' ) 1644 1644 1645 1645 DO k = k_wall + 1, k_wall + pch_index_ji(j,i) … … 1667 1667 ! 1668 1668 !-- Determine topography-top index on scalar grid 1669 k_wall = get_topography_top_index ( j, i, 's' )1669 k_wall = get_topography_top_index_ji( j, i, 's' ) 1670 1670 1671 1671 DO k = k_wall + 1, k_wall + pch_index_ji(j,i) … … 1693 1693 ! 1694 1694 !-- Determine topography-top index on scalar grid 1695 k_wall = get_topography_top_index ( j, i, 's' )1695 k_wall = get_topography_top_index_ji( j, i, 's' ) 1696 1696 1697 1697 DO k = k_wall + 1, k_wall + pch_index_ji(j,i) -
palm/trunk/SOURCE/pmc_interface_mod.f90
r2696 r2698 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Bugfix in get_topography_top_index 23 23 ! 24 24 ! Former revisions: … … 265 265 266 266 USE surface_mod, & 267 ONLY: get_topography_top_index , surf_def_h, surf_lsm_h, surf_usm_h267 ONLY: get_topography_top_index_ji, surf_def_h, surf_lsm_h, surf_usm_h 268 268 269 269 IMPLICIT NONE … … 1410 1410 !-- Determine largest topography index on scalar grid 1411 1411 nzt_topo_nestbc_l = MAX( nzt_topo_nestbc_l, & 1412 get_topography_top_index ( j, i, 's' ) )1412 get_topography_top_index_ji( j, i, 's' ) ) 1413 1413 ! 1414 1414 !-- Determine largest topography index on u grid 1415 1415 nzt_topo_nestbc_l = MAX( nzt_topo_nestbc_l, & 1416 get_topography_top_index ( j, i, 'u' ) )1416 get_topography_top_index_ji( j, i, 'u' ) ) 1417 1417 ! 1418 1418 !-- Determine largest topography index on v grid 1419 1419 nzt_topo_nestbc_l = MAX( nzt_topo_nestbc_l, & 1420 get_topography_top_index ( j, i, 'v' ) )1420 get_topography_top_index_ji( j, i, 'v' ) ) 1421 1421 ! 1422 1422 !-- Determine largest topography index on w grid 1423 1423 nzt_topo_nestbc_l = MAX( nzt_topo_nestbc_l, & 1424 get_topography_top_index ( j, i, 'w' ) )1424 get_topography_top_index_ji( j, i, 'w' ) ) 1425 1425 ENDDO 1426 1426 ENDDO … … 1436 1436 !-- Determine largest topography index on scalar grid 1437 1437 nzt_topo_nestbc_r = MAX( nzt_topo_nestbc_r, & 1438 get_topography_top_index ( j, i, 's' ) )1438 get_topography_top_index_ji( j, i, 's' ) ) 1439 1439 ! 1440 1440 !-- Determine largest topography index on u grid 1441 1441 nzt_topo_nestbc_r = MAX( nzt_topo_nestbc_r, & 1442 get_topography_top_index ( j, i, 'u' ) )1442 get_topography_top_index_ji( j, i, 'u' ) ) 1443 1443 ! 1444 1444 !-- Determine largest topography index on v grid 1445 1445 nzt_topo_nestbc_r = MAX( nzt_topo_nestbc_r, & 1446 get_topography_top_index ( j, i, 'v' ) )1446 get_topography_top_index_ji( j, i, 'v' ) ) 1447 1447 ! 1448 1448 !-- Determine largest topography index on w grid 1449 1449 nzt_topo_nestbc_r = MAX( nzt_topo_nestbc_r, & 1450 get_topography_top_index ( j, i, 'w' ) )1450 get_topography_top_index_ji( j, i, 'w' ) ) 1451 1451 ENDDO 1452 1452 nzt_topo_nestbc_r = nzt_topo_nestbc_r + 1 … … 1461 1461 !-- Determine largest topography index on scalar grid 1462 1462 nzt_topo_nestbc_s = MAX( nzt_topo_nestbc_s, & 1463 get_topography_top_index ( j, i, 's' ) )1463 get_topography_top_index_ji( j, i, 's' ) ) 1464 1464 ! 1465 1465 !-- Determine largest topography index on u grid 1466 1466 nzt_topo_nestbc_s = MAX( nzt_topo_nestbc_s, & 1467 get_topography_top_index ( j, i, 'u' ) )1467 get_topography_top_index_ji( j, i, 'u' ) ) 1468 1468 ! 1469 1469 !-- Determine largest topography index on v grid 1470 1470 nzt_topo_nestbc_s = MAX( nzt_topo_nestbc_s, & 1471 get_topography_top_index ( j, i, 'v' ) )1471 get_topography_top_index_ji( j, i, 'v' ) ) 1472 1472 ! 1473 1473 !-- Determine largest topography index on w grid 1474 1474 nzt_topo_nestbc_s = MAX( nzt_topo_nestbc_s, & 1475 get_topography_top_index ( j, i, 'w' ) )1475 get_topography_top_index_ji( j, i, 'w' ) ) 1476 1476 ENDDO 1477 1477 ENDDO … … 1487 1487 !-- Determine largest topography index on scalar grid 1488 1488 nzt_topo_nestbc_n = MAX( nzt_topo_nestbc_n, & 1489 get_topography_top_index ( j, i, 's' ) )1489 get_topography_top_index_ji( j, i, 's' ) ) 1490 1490 ! 1491 1491 !-- Determine largest topography index on u grid 1492 1492 nzt_topo_nestbc_n = MAX( nzt_topo_nestbc_n, & 1493 get_topography_top_index ( j, i, 'u' ) )1493 get_topography_top_index_ji( j, i, 'u' ) ) 1494 1494 ! 1495 1495 !-- Determine largest topography index on v grid 1496 1496 nzt_topo_nestbc_n = MAX( nzt_topo_nestbc_n, & 1497 get_topography_top_index ( j, i, 'v' ) )1497 get_topography_top_index_ji( j, i, 'v' ) ) 1498 1498 ! 1499 1499 !-- Determine largest topography index on w grid 1500 1500 nzt_topo_nestbc_n = MAX( nzt_topo_nestbc_n, & 1501 get_topography_top_index ( j, i, 'w' ) )1501 get_topography_top_index_ji( j, i, 'w' ) ) 1502 1502 ENDDO 1503 1503 nzt_topo_nestbc_n = nzt_topo_nestbc_n + 1 … … 1552 1552 !-- is part of the surfacetypes now. Set default roughness instead. 1553 1553 !-- Determine topography top index on u-grid 1554 kb = get_topography_top_index ( j, i, 'u' )1554 kb = get_topography_top_index_ji( j, i, 'u' ) 1555 1555 k = kb + 1 1556 1556 wall_index = kb … … 1568 1568 ! 1569 1569 !-- Determine topography top index on v-grid 1570 kb = get_topography_top_index ( j, i, 'v' )1570 kb = get_topography_top_index_ji( j, i, 'v' ) 1571 1571 k = kb + 1 1572 1572 wall_index = kb … … 1611 1611 !-- to the present surface tpye. 1612 1612 !-- Determine topography top index on u-grid 1613 kb = get_topography_top_index ( j, i, 'u' )1613 kb = get_topography_top_index_ji( j, i, 'u' ) 1614 1614 k = kb + 1 1615 1615 wall_index = kb … … 1627 1627 ! 1628 1628 !-- Determine topography top index on v-grid 1629 kb = get_topography_top_index ( j, i, 'v' )1629 kb = get_topography_top_index_ji( j, i, 'v' ) 1630 1630 k = kb + 1 1631 1631 wall_index = kb … … 1667 1667 ! 1668 1668 !-- Determine topography top index on u-grid 1669 kb = get_topography_top_index ( j, i, 'u' )1669 kb = get_topography_top_index_ji( j, i, 'u' ) 1670 1670 k = kb + 1 1671 1671 wall_index = kb … … 1683 1683 ! 1684 1684 !-- Determine topography top index on v-grid 1685 kb = get_topography_top_index ( j, i, 'v' )1685 kb = get_topography_top_index_ji( j, i, 'v' ) 1686 1686 k = kb + 1 1687 1687 wall_index = kb … … 1723 1723 ! 1724 1724 !-- Determine topography top index on u-grid 1725 kb = get_topography_top_index ( j, i, 'u' )1725 kb = get_topography_top_index_ji( j, i, 'u' ) 1726 1726 k = kb + 1 1727 1727 wall_index = kb … … 1739 1739 ! 1740 1740 !-- Determine topography top index on v-grid 1741 kb = get_topography_top_index ( j, i, 'v' )1741 kb = get_topography_top_index_ji( j, i, 'v' ) 1742 1742 k = kb + 1 1743 1743 wall_index = kb … … 1779 1779 ! 1780 1780 !-- Determine lowest grid on outer grids for u and w. 1781 k_wall_u_ji = get_topography_top_index ( j, 0, 'u_out' )1782 k_wall_u_ji_p = get_topography_top_index ( j+1, 0, 'u_out' )1783 k_wall_u_ji_m = get_topography_top_index ( j-1, 0, 'u_out' )1784 1785 k_wall_w_ji = get_topography_top_index ( j, -1, 'w_out' )1786 k_wall_w_ji_p = get_topography_top_index ( j+1, -1, 'w_out' )1787 k_wall_w_ji_m = get_topography_top_index ( j-1, -1, 'w_out' )1781 k_wall_u_ji = get_topography_top_index_ji( j, 0, 'u_out' ) 1782 k_wall_u_ji_p = get_topography_top_index_ji( j+1, 0, 'u_out' ) 1783 k_wall_u_ji_m = get_topography_top_index_ji( j-1, 0, 'u_out' ) 1784 1785 k_wall_w_ji = get_topography_top_index_ji( j, -1, 'w_out' ) 1786 k_wall_w_ji_p = get_topography_top_index_ji( j+1, -1, 'w_out' ) 1787 k_wall_w_ji_m = get_topography_top_index_ji( j-1, -1, 'w_out' ) 1788 1788 1789 1789 DO k = nzb, nzt_topo_nestbc_l … … 1869 1869 ! 1870 1870 !-- Determine lowest grid on outer grids for u and w. 1871 k_wall_u_ji = get_topography_top_index ( j, i, 'u_out' )1872 k_wall_u_ji_p = get_topography_top_index ( j+1, i, 'u_out' )1873 k_wall_u_ji_m = get_topography_top_index ( j-1, i, 'u_out' )1874 1875 k_wall_w_ji = get_topography_top_index ( j, i, 'w_out' )1876 k_wall_w_ji_p = get_topography_top_index ( j+1, i, 'w_out' )1877 k_wall_w_ji_m = get_topography_top_index ( j-1, i, 'w_out' )1871 k_wall_u_ji = get_topography_top_index_ji( j, i, 'u_out' ) 1872 k_wall_u_ji_p = get_topography_top_index_ji( j+1, i, 'u_out' ) 1873 k_wall_u_ji_m = get_topography_top_index_ji( j-1, i, 'u_out' ) 1874 1875 k_wall_w_ji = get_topography_top_index_ji( j, i, 'w_out' ) 1876 k_wall_w_ji_p = get_topography_top_index_ji( j+1, i, 'w_out' ) 1877 k_wall_w_ji_m = get_topography_top_index_ji( j-1, i, 'w_out' ) 1878 1878 1879 1879 DO k = nzb, nzt_topo_nestbc_r … … 1952 1952 ! 1953 1953 !-- Determine lowest grid on outer grids for v and w. 1954 k_wall_v_ji = get_topography_top_index ( 0, i, 'v_out' )1955 k_wall_v_ji_p = get_topography_top_index ( 0, i+1, 'v_out' )1956 k_wall_v_ji_m = get_topography_top_index ( 0, i-1, 'v_out' )1957 1958 k_wall_w_ji = get_topography_top_index ( -1, i, 'w_out' )1959 k_wall_w_ji_p = get_topography_top_index ( -1, i+1, 'w_out' )1960 k_wall_w_ji_m = get_topography_top_index ( -1, i-1, 'w_out' )1954 k_wall_v_ji = get_topography_top_index_ji( 0, i, 'v_out' ) 1955 k_wall_v_ji_p = get_topography_top_index_ji( 0, i+1, 'v_out' ) 1956 k_wall_v_ji_m = get_topography_top_index_ji( 0, i-1, 'v_out' ) 1957 1958 k_wall_w_ji = get_topography_top_index_ji( -1, i, 'w_out' ) 1959 k_wall_w_ji_p = get_topography_top_index_ji( -1, i+1, 'w_out' ) 1960 k_wall_w_ji_m = get_topography_top_index_ji( -1, i-1, 'w_out' ) 1961 1961 1962 1962 DO k = nzb, nzt_topo_nestbc_s … … 2041 2041 ! 2042 2042 !-- Determine lowest grid on outer grids for v and w. 2043 k_wall_v_ji = get_topography_top_index ( j, i, 'v_out' )2044 k_wall_v_ji_p = get_topography_top_index ( j, i+1, 'v_out' )2045 k_wall_v_ji_m = get_topography_top_index ( j, i-1, 'v_out' )2046 2047 k_wall_w_ji = get_topography_top_index ( j, i, 'w_out' )2048 k_wall_w_ji_p = get_topography_top_index ( j, i+1, 'w_out' )2049 k_wall_w_ji_m = get_topography_top_index ( j, i-1, 'w_out' )2043 k_wall_v_ji = get_topography_top_index_ji( j, i, 'v_out' ) 2044 k_wall_v_ji_p = get_topography_top_index_ji( j, i+1, 'v_out' ) 2045 k_wall_v_ji_m = get_topography_top_index_ji( j, i-1, 'v_out' ) 2046 2047 k_wall_w_ji = get_topography_top_index_ji( j, i, 'w_out' ) 2048 k_wall_w_ji_p = get_topography_top_index_ji( j, i+1, 'w_out' ) 2049 k_wall_w_ji_m = get_topography_top_index_ji( j, i-1, 'w_out' ) 2050 2050 2051 2051 DO k = nzb, nzt_topo_nestbc_n … … 2643 2643 i = nxl - 1 2644 2644 DO j = nysg, nyng 2645 k_wall = get_topography_top_index ( j, i, 's' )2645 k_wall = get_topography_top_index_ji( j, i, 's' ) 2646 2646 2647 2647 DO k = k_wall + 1, nzt … … 2664 2664 i = nxr + 1 2665 2665 DO j = nysg, nyng 2666 k_wall = get_topography_top_index ( j, i, 's' )2666 k_wall = get_topography_top_index_ji( j, i, 's' ) 2667 2667 2668 2668 DO k = k_wall + 1, nzt … … 2685 2685 j = nys - 1 2686 2686 DO i = nxlg, nxrg 2687 k_wall = get_topography_top_index ( j, i, 's' )2687 k_wall = get_topography_top_index_ji( j, i, 's' ) 2688 2688 2689 2689 DO k = k_wall + 1, nzt … … 2706 2706 j = nyn + 1 2707 2707 DO i = nxlg, nxrg 2708 k_wall = get_topography_top_index ( j, i, 's' )2708 k_wall = get_topography_top_index_ji( j, i, 's' ) 2709 2709 2710 2710 DO k = k_wall + 1, nzt … … 2729 2729 ! 2730 2730 !-- Determine vertical index for local topography top 2731 k_wall = get_topography_top_index ( j, i, 's' )2731 k_wall = get_topography_top_index_ji( j, i, 's' ) 2732 2732 2733 2733 kc = kco(k) + 1 … … 3196 3196 ! 3197 3197 !-- Determine vertical index of topography top at grid point (j,i) 3198 k_wall = get_topography_top_index ( j, i, TRIM ( var ) )3198 k_wall = get_topography_top_index_ji( j, i, TRIM ( var ) ) 3199 3199 ! 3200 3200 !-- kbc is the first coarse-grid point above the surface … … 3227 3227 ! 3228 3228 !-- Determine vertical index of topography top at grid point (j,i) 3229 k_wall = get_topography_top_index ( j, i, 'w' )3229 k_wall = get_topography_top_index_ji( j, i, 'w' ) 3230 3230 3231 3231 f(k_wall,j,i) = 0.0_wp … … 4369 4369 ! 4370 4370 !-- Determine vertical index of topography top at grid point (j,i) 4371 k_wall = get_topography_top_index ( j, i, TRIM( var ) )4371 k_wall = get_topography_top_index_ji( j, i, TRIM( var ) ) 4372 4372 4373 4373 DO k = k_wall, nzt+1 … … 4395 4395 ! 4396 4396 !-- Determine vertical index of topography top at grid point (j,i) 4397 k_wall = get_topography_top_index ( j, i, TRIM ( var ) )4397 k_wall = get_topography_top_index_ji( j, i, TRIM ( var ) ) 4398 4398 4399 4399 k = k_wall+1 … … 4420 4420 ! 4421 4421 !-- Determine vertical index of topography top at grid point (j,i) 4422 k_wall = get_topography_top_index ( j, i, TRIM ( var ) )4422 k_wall = get_topography_top_index_ji( j, i, TRIM ( var ) ) 4423 4423 4424 4424 DO k = k_wall+1, nzt_topo_nestbc … … 4445 4445 ! 4446 4446 !-- Determine vertical index of topography top at grid point (j,i) 4447 k_wall = get_topography_top_index ( j, i, TRIM ( var ) )4447 k_wall = get_topography_top_index_ji( j, i, TRIM ( var ) ) 4448 4448 4449 4449 k = k_wall + 1 … … 4476 4476 ! 4477 4477 !-- Determine vertical index of topography top at grid point (j,i) 4478 k_wall = get_topography_top_index ( j, i, 's' )4478 k_wall = get_topography_top_index_ji( j, i, 's' ) 4479 4479 4480 4480 DO k = k_wall, nzt + 1 … … 4486 4486 ! 4487 4487 !-- Determine vertical index of topography top at grid point (j,i) 4488 k_wall = get_topography_top_index ( j, i, 's' )4488 k_wall = get_topography_top_index_ji( j, i, 's' ) 4489 4489 4490 4490 DO k = k_wall, nzt+1 … … 4593 4593 ! 4594 4594 !-- Determine vertical index of topography top at grid point (j,i) 4595 k_wall = get_topography_top_index ( j, i, TRIM( var ) )4595 k_wall = get_topography_top_index_ji( j, i, TRIM( var ) ) 4596 4596 4597 4597 DO k = k_wall, nzt+1 … … 4619 4619 ! 4620 4620 !-- Determine vertical index of topography top at grid point (j,i) 4621 k_wall = get_topography_top_index ( j, i, TRIM( var ) )4621 k_wall = get_topography_top_index_ji( j, i, TRIM( var ) ) 4622 4622 4623 4623 k = k_wall + 1 … … 4643 4643 ! 4644 4644 !-- Determine vertical index of topography top at grid point (j,i) 4645 k_wall = get_topography_top_index ( j, i, TRIM( var ) )4645 k_wall = get_topography_top_index_ji( j, i, TRIM( var ) ) 4646 4646 4647 4647 DO k = k_wall, nzt_topo_nestbc … … 4670 4670 ! 4671 4671 !-- Determine vertical index of topography top at grid point (j,i) 4672 k_wall = get_topography_top_index ( j, i, TRIM( var ) )4672 k_wall = get_topography_top_index_ji( j, i, TRIM( var ) ) 4673 4673 4674 4674 k = k_wall + 1 … … 4701 4701 ! 4702 4702 !-- Determine vertical index of topography top at grid point (j,i) 4703 k_wall = get_topography_top_index ( j, i, 's' )4703 k_wall = get_topography_top_index_ji( j, i, 's' ) 4704 4704 DO k = k_wall, nzt+1 4705 4705 f(k,j,i) = tkefactor_s(k,i) * f(k,j,i) … … 4710 4710 ! 4711 4711 !-- Determine vertical index of topography top at grid point (j,i) 4712 k_wall = get_topography_top_index ( j, i, 's' )4712 k_wall = get_topography_top_index_ji( j, i, 's' ) 4713 4713 DO k = k_wall, nzt+1 4714 4714 f(k,j,i) = tkefactor_n(k,i) * f(k,j,i) … … 4891 4891 ! 4892 4892 !-- Determine vertical index of topography top at grid point (j,i) 4893 k_wall = get_topography_top_index ( j, i, TRIM( var ) )4893 k_wall = get_topography_top_index_ji( j, i, TRIM( var ) ) 4894 4894 4895 4895 DO k = k_wall, nzt+1 … … 4970 4970 ! 4971 4971 !-- Determine vertical index of topography top at grid point (j,i) 4972 k_wall = get_topography_top_index ( j, i, TRIM( var ) )4972 k_wall = get_topography_top_index_ji( j, i, TRIM( var ) ) 4973 4973 4974 4974 DO k = k_wall, nzt+1 -
palm/trunk/SOURCE/radiation_model_mod.f90
r2696 r2698 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Bugfix in get_topography_top_index 23 23 ! 24 24 ! Former revisions: … … 272 272 273 273 USE surface_mod, & 274 ONLY: get_topography_top_index, surf_def_h, surf_def_v, surf_lsm_h, & 274 ONLY: get_topography_top_index, get_topography_top_index_ji, & 275 surf_def_h, surf_def_v, surf_lsm_h, & 275 276 surf_lsm_v, surf_type, surf_usm_h, surf_usm_v 276 277 … … 3244 3245 ! 3245 3246 !-- Obtain topography top index (lower bound of RRTMG) 3246 k_topo = get_topography_top_index ( j, i, 's' )3247 k_topo = get_topography_top_index_ji( j, i, 's' ) 3247 3248 3248 3249 IF ( lw_radiation ) THEN … … 4411 4412 pc_heating_rate, lad_s, prototype_lad, usm_lad_rma 4412 4413 4413 USE surface_mod, &4414 ONLY: get_topography_top_index, surf_lsm_h, surf_lsm_v, surf_usm_h,&4415 surf_usm_v4416 4417 4414 IMPLICIT NONE 4418 4415 … … 4448 4445 ! 4449 4446 !-- Find topography top index 4450 k_topo = get_topography_top_index ( j, i, 's' )4447 k_topo = get_topography_top_index_ji( j, i, 's' ) 4451 4448 4452 4449 DO k = nzt+1, 0, -1 … … 4553 4550 DO j = ijdb(3,ids), ijdb(4,ids) 4554 4551 4555 k_topo = get_topography_top_index ( j, i, 's' )4556 k_topo2 = get_topography_top_index ( j-jdir(ids), i-idir(ids), 's' )4552 k_topo = get_topography_top_index_ji( j, i, 's' ) 4553 k_topo2 = get_topography_top_index_ji( j-jdir(ids), i-idir(ids), 's' ) 4557 4554 4558 4555 … … 4572 4569 DO j = nys, nyn 4573 4570 !-- Find topography top index 4574 k_topo = get_topography_top_index ( j, i, 's' )4571 k_topo = get_topography_top_index_ji( j, i, 's' ) 4575 4572 k = nzut - k_topo 4576 4573 nsurfl = nsurfl + 6 * k … … 4593 4590 ! 4594 4591 !-- Find topography top index 4595 k_topo = get_topography_top_index ( j, i, 's' )4592 k_topo = get_topography_top_index_ji( j, i, 's' ) 4596 4593 4597 4594 DO k = k_topo + 1, pct(j,i) … … 4709 4706 DO i = ijdb(1,ids), ijdb(2,ids) 4710 4707 DO j = ijdb(3,ids), ijdb(4,ids) 4711 k_topo = get_topography_top_index ( j, i, 's' )4712 k_topo2 = get_topography_top_index ( j-jdir(ids), i-idir(ids), 's' )4708 k_topo = get_topography_top_index_ji( j, i, 's' ) 4709 k_topo2 = get_topography_top_index_ji( j-jdir(ids), i-idir(ids), 's' ) 4713 4710 4714 4711 DO k = MAX(k_topo,k_topo2)+1, nzut … … 4727 4724 DO i = nxl, nxr 4728 4725 DO j = nys, nyn 4729 k_topo = get_topography_top_index ( j, i, 's' )4726 k_topo = get_topography_top_index_ji( j, i, 's' ) 4730 4727 4731 4728 !-- add upward surface … … 4746 4743 DO i = nxl, nxr 4747 4744 DO j = nys, nyn 4748 k_topo = get_topography_top_index ( j, i, 's' )4745 k_topo = get_topography_top_index_ji( j, i, 's' ) 4749 4746 !-- north 4750 4747 IF ( j /= ny ) THEN … … 4752 4749 jr = min(max(j-jdir(ids),0),ny) 4753 4750 ir = min(max(i-idir(ids),0),nx) 4754 k_topo2 = get_topography_top_index ( jr, ir, 's' )4751 k_topo2 = get_topography_top_index_ji( jr, ir, 's' ) 4755 4752 DO k = MAX(k_topo,k_topo2)+1, nzut 4756 4753 isurf = isurf + 1 … … 4763 4760 jr = min(max(j-jdir(ids),0),ny) 4764 4761 ir = min(max(i-idir(ids),0),nx) 4765 k_topo2 = get_topography_top_index ( jr, ir, 's' )4762 k_topo2 = get_topography_top_index_ji( jr, ir, 's' ) 4766 4763 4767 4764 DO k = MAX(k_topo,k_topo2)+1, nzut … … 4775 4772 jr = min(max(j-jdir(ids),0),ny) 4776 4773 ir = min(max(i-idir(ids),0),nx) 4777 k_topo2 = get_topography_top_index ( jr, ir, 's' )4774 k_topo2 = get_topography_top_index_ji( jr, ir, 's' ) 4778 4775 4779 4776 DO k = MAX(k_topo,k_topo2)+1, nzut … … 4787 4784 jr = min(max(j-jdir(ids),0),ny) 4788 4785 ir = min(max(i-idir(ids),0),nx) 4789 k_topo2 = get_topography_top_index ( jr, ir, 's' )4786 k_topo2 = get_topography_top_index_ji( jr, ir, 's' ) 4790 4787 4791 4788 DO k = MAX(k_topo,k_topo2)+1, nzut … … 5255 5252 ! 5256 5253 !-- Following expression equals former kk = k - nzb_s_inner(j,i) 5257 kk = k - get_topography_top_index ( j, i, 's' ) !- lad arrays are defined flat5254 kk = k - get_topography_top_index_ji( j, i, 's' ) !- lad arrays are defined flat 5258 5255 pc_heating_rate(kk, j, i) = (pcbinsw(ipcgb) + pcbinlw(ipcgb)) & 5259 5256 * pchf_prep(k) * pt(k, j, i) !-- = dT/dt … … 5856 5853 DO i = nxl, nxr 5857 5854 DO j = nys, nyn 5858 k = get_topography_top_index ( j, i, 's' )5855 k = get_topography_top_index_ji( j, i, 's' ) 5859 5856 5860 5857 usm_lad(k:nzut, j, i) = lad_s(0:nzut-k, j, i) -
palm/trunk/SOURCE/surface_mod.f90
r2696 r2698 391 391 PRIVATE 392 392 393 INTERFACE get_topography_top_index394 MODULE PROCEDURE get_topography_top_index395 MODULE PROCEDURE get_topography_top_index_ji396 END INTERFACE get_topography_top_index397 398 393 INTERFACE init_bc 399 394 MODULE PROCEDURE init_bc … … 431 426 ! 432 427 !-- Public subroutines and functions 433 PUBLIC get_topography_top_index, init_bc, init_surfaces, & 428 PUBLIC get_topography_top_index, get_topography_top_index_ji, init_bc, & 429 init_surfaces, & 434 430 init_surface_arrays, surface_read_restart_data, & 435 431 surface_restore_elements, surface_write_restart_data, & -
palm/trunk/SOURCE/turbulence_closure_mod.f90
r2696 r2698 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Bugfix in get_topography_top_index 23 23 ! 24 24 ! Former revisions: … … 870 870 871 871 USE surface_mod, & 872 ONLY: get_topography_top_index 872 ONLY: get_topography_top_index_ji 873 873 874 874 IMPLICIT NONE … … 972 972 DO i = nxlg, nxrg 973 973 DO j = nysg, nyng 974 nz_s_shift = get_topography_top_index ( j, i, 's' )974 nz_s_shift = get_topography_top_index_ji( j, i, 's' ) 975 975 976 976 e(nz_s_shift:nzt+1,j,i) = e(0:nzt+1-nz_s_shift,j,i) … … 991 991 IF ( complex_terrain ) THEN 992 992 IF ( nxlg <= 0 .AND. nxrg >= 0 .AND. nysg <= 0 .AND. nyng >= 0 ) THEN 993 nz_s_shift_l = get_topography_top_index ( 0, 0, 's' )993 nz_s_shift_l = get_topography_top_index_ji( 0, 0, 's' ) 994 994 ELSE 995 995 nz_s_shift_l = 0
Note: See TracChangeset
for help on using the changeset viewer.
