- Timestamp:
- Nov 20, 2017 12:40:38 PM (7 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 5 edited
-
check_parameters.f90 (modified) (2 diffs)
-
lpm_advec.f90 (modified) (8 diffs)
-
lpm_exchange_horiz.f90 (modified) (8 diffs)
-
lpm_init.f90 (modified) (4 diffs)
-
lpm_pack_arrays.f90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/check_parameters.f90
r2575 r2628 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Enabled particle advection with grid stretching -> Removed parameter check 28 ! 29 ! 2575 2017-10-24 09:57:58Z maronga 27 30 ! Renamed phi --> latitude 28 31 ! … … 880 883 881 884 ! 882 !-- Check if vertical grid stretching is used together with particles883 IF ( dz_stretch_level < 100000.0_wp .AND. particle_advection ) THEN884 message_string = 'Vertical grid stretching is not allowed together ' // &885 'with particle advection.'886 CALL message( 'check_parameters', 'PA0017', 1, 2, 0, 6, 0 )887 ENDIF888 889 !890 885 !-- Check topography setting (check for illegal parameter combinations) 891 886 IF ( topography /= 'flat' ) THEN -
palm/trunk/SOURCE/lpm_advec.f90
r2610 r2628 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Removed indices ilog and jlog which are no longer needed since particle box 23 ! locations are identical to scalar boxes and topography. 24 23 25 ! 24 26 ! Former revisions: … … 167 169 INTEGER(iwp) :: i !< index variable along x 168 170 INTEGER(iwp) :: ip !< index variable along x 169 INTEGER(iwp) :: ilog !< index variable along x170 171 INTEGER(iwp) :: j !< index variable along y 171 172 INTEGER(iwp) :: jp !< index variable along y 172 INTEGER(iwp) :: jlog !< index variable along y173 173 INTEGER(iwp) :: k !< index variable along z 174 174 INTEGER(iwp) :: k_wall !< vertical index of topography top … … 299 299 !-- First, check if particle is located below first vertical grid level 300 300 !-- above topography (Prandtl-layer height) 301 ilog = particles(n)%x * ddx302 jlog = particles(n)%y * ddy303 !304 301 !-- Determine vertical index of topography top 305 k_wall = get_topography_top_index( j log, ilog, 's' )302 k_wall = get_topography_top_index( jp, ip, 's' ) 306 303 307 304 IF ( constant_flux_layer .AND. zv(n) - zw(k_wall) < z_p ) THEN … … 327 324 !-- Get friction velocity and momentum flux from new surface data 328 325 !-- types. 329 IF ( surf_def_h(0)%start_index(j log,ilog) <= &330 surf_def_h(0)%end_index(j log,ilog) ) THEN331 surf_start = surf_def_h(0)%start_index(j log,ilog)326 IF ( surf_def_h(0)%start_index(jp,ip) <= & 327 surf_def_h(0)%end_index(jp,ip) ) THEN 328 surf_start = surf_def_h(0)%start_index(jp,ip) 332 329 !-- Limit friction velocity. In narrow canyons or holes the 333 330 !-- friction velocity can become very small, resulting in a too … … 335 332 us_int = MAX( surf_def_h(0)%us(surf_start), 0.01_wp ) 336 333 usws_int = surf_def_h(0)%usws(surf_start) 337 ELSEIF ( surf_lsm_h%start_index(j log,ilog) <= &338 surf_lsm_h%end_index(j log,ilog) ) THEN339 surf_start = surf_lsm_h%start_index(j log,ilog)334 ELSEIF ( surf_lsm_h%start_index(jp,ip) <= & 335 surf_lsm_h%end_index(jp,ip) ) THEN 336 surf_start = surf_lsm_h%start_index(jp,ip) 340 337 us_int = MAX( surf_lsm_h%us(surf_start), 0.01_wp ) 341 338 usws_int = surf_lsm_h%usws(surf_start) 342 ELSEIF ( surf_usm_h%start_index(j log,ilog) <= &343 surf_usm_h%end_index(j log,ilog) ) THEN344 surf_start = surf_usm_h%start_index(j log,ilog)339 ELSEIF ( surf_usm_h%start_index(jp,ip) <= & 340 surf_usm_h%end_index(jp,ip) ) THEN 341 surf_start = surf_usm_h%start_index(jp,ip) 345 342 us_int = MAX( surf_usm_h%us(surf_start), 0.01_wp ) 346 343 usws_int = surf_usm_h%usws(surf_start) … … 396 393 DO n = start_index(nb), end_index(nb) 397 394 398 ilog = particles(n)%x * ddx399 jlog = particles(n)%y * ddy400 395 ! 401 396 !-- Determine vertical index of topography top 402 k_wall = get_topography_top_index( j log,ilog, 's' )397 k_wall = get_topography_top_index( jp,ip, 's' ) 403 398 404 399 IF ( constant_flux_layer .AND. zv(n) - zw(k_wall) < z_p ) THEN … … 427 422 !-- Get friction velocity and momentum flux from new surface data 428 423 !-- types. 429 IF ( surf_def_h(0)%start_index(j log,ilog) <= &430 surf_def_h(0)%end_index(j log,ilog) ) THEN431 surf_start = surf_def_h(0)%start_index(j log,ilog)424 IF ( surf_def_h(0)%start_index(jp,ip) <= & 425 surf_def_h(0)%end_index(jp,ip) ) THEN 426 surf_start = surf_def_h(0)%start_index(jp,ip) 432 427 !-- Limit friction velocity. In narrow canyons or holes the 433 428 !-- friction velocity can become very small, resulting in a too … … 435 430 us_int = MAX( surf_def_h(0)%us(surf_start), 0.01_wp ) 436 431 vsws_int = surf_def_h(0)%vsws(surf_start) 437 ELSEIF ( surf_lsm_h%start_index(j log,ilog) <= &438 surf_lsm_h%end_index(j log,ilog) ) THEN439 surf_start = surf_lsm_h%start_index(j log,ilog)432 ELSEIF ( surf_lsm_h%start_index(jp,ip) <= & 433 surf_lsm_h%end_index(jp,ip) ) THEN 434 surf_start = surf_lsm_h%start_index(jp,ip) 440 435 us_int = MAX( surf_lsm_h%us(surf_start), 0.01_wp ) 441 436 vsws_int = surf_lsm_h%vsws(surf_start) 442 ELSEIF ( surf_usm_h%start_index(j log,ilog) <= &443 surf_usm_h%end_index(j log,ilog) ) THEN444 surf_start = surf_usm_h%start_index(j log,ilog)437 ELSEIF ( surf_usm_h%start_index(jp,ip) <= & 438 surf_usm_h%end_index(jp,ip) ) THEN 439 surf_start = surf_usm_h%start_index(jp,ip) 445 440 us_int = MAX( surf_usm_h%us(surf_start), 0.01_wp ) 446 441 vsws_int = surf_usm_h%vsws(surf_start) -
palm/trunk/SOURCE/lpm_exchange_horiz.f90
r2606 r2628 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Enabled particle advection with grid stretching. Furthermore, the CFL- 28 ! criterion is checked for every particle at every time step. 29 ! 30 ! 2606 2017-11-10 10:36:31Z schwenkel 27 31 ! Changed particle box locations: center of particle box now coincides 28 32 ! with scalar grid point of same index. … … 124 128 125 129 USE, INTRINSIC :: ISO_C_BINDING 126 130 131 USE arrays_3d, & 132 ONLY: zw 133 127 134 USE control_parameters, & 128 135 ONLY: dz, message_string, simulated_time … … 887 894 jp = particle_array(n)%y * ddy 888 895 kp = particle_array(n)%z / dz + 1 + offset_ocean_nzt 896 ! 897 !-- In case of grid stretching a particle might be above or below the 898 !-- previously calculated particle grid box (indices). 899 DO WHILE( zw(kp) < particle_array(n)%z ) 900 kp = kp + 1 901 ENDDO 902 903 DO WHILE( zw(kp-1) > particle_array(n)%z ) 904 kp = kp - 1 905 ENDDO 889 906 890 907 IF ( ip >= nxl .AND. ip <= nxr .AND. jp >= nys .AND. jp <= nyn & … … 1022 1039 !------------------------------------------------------------------------------! 1023 1040 SUBROUTINE lpm_move_particle 1024 1041 1025 1042 IMPLICIT NONE 1026 1043 … … 1038 1055 1039 1056 CALL cpu_log( log_point_s(41), 'lpm_move_particle', 'start' ) 1040 1057 CALL lpm_check_cfl 1041 1058 DO ip = nxl, nxr 1042 1059 DO jp = nys, nyn … … 1050 1067 i = particles_before_move(n)%x * ddx 1051 1068 j = particles_before_move(n)%y * ddy 1052 k = particles_before_move(n)%z / dz + 1 + offset_ocean_nzt 1053 1069 k = kp 1070 ! 1071 !-- Find correct vertical particle grid box (necessary in case of grid stretching) 1072 !-- Due to the CFL limitations only the neighbouring grid boxes are considered. 1073 IF( zw(k) < particles_before_move(n)%z ) k = k + 1 1074 IF( zw(k-1) > particles_before_move(n)%z ) k = k - 1 1075 1054 1076 !-- For lpm_exchange_horiz to work properly particles need to be moved to the outermost gridboxes 1055 1077 !-- of the respective processor. If the particle index is inside the processor the following lines … … 1059 1081 j = MIN ( j , nyn ) 1060 1082 j = MAX ( j , nys ) 1083 1061 1084 k = MIN ( k , nzt ) 1062 1085 k = MAX ( k , nzb+1 ) 1086 1063 1087 ! 1064 1088 !-- Check, if particle has moved to another grid cell. … … 1092 1116 1093 1117 END SUBROUTINE lpm_move_particle 1094 1118 1119 !------------------------------------------------------------------------------! 1120 ! Description: 1121 ! ------------ 1122 !> Check CFL-criterion for each particle. If one particle violated the 1123 !> criterion the particle will be deleted and a warning message is given. 1124 !------------------------------------------------------------------------------! 1125 SUBROUTINE lpm_check_cfl 1126 1127 IMPLICIT NONE 1128 1129 INTEGER(iwp) :: i 1130 INTEGER(iwp) :: j 1131 INTEGER(iwp) :: k 1132 INTEGER(iwp) :: n 1133 1134 DO i = nxl, nxr 1135 DO j = nys, nyn 1136 DO k = nzb+1, nzt 1137 number_of_particles = prt_count(k,j,i) 1138 IF ( number_of_particles <= 0 ) CYCLE 1139 particles => grid_particles(k,j,i)%particles(1:number_of_particles) 1140 DO n = 1, number_of_particles 1141 1142 IF(ABS(particles(n)%speed_x) > & 1143 (dx/(particles(n)%age-particles(n)%age_m)) .OR. & 1144 ABS(particles(n)%speed_y) > & 1145 (dx/(particles(n)%age-particles(n)%age_m)) .OR. & 1146 ABS(particles(n)%speed_z) > & 1147 ((zw(k)-zw(k-1))/(particles(n)%age-particles(n)%age_m))) THEN 1148 WRITE( message_string, * ) 'PARTICLE VIOLATED CFL CRITERION'& 1149 ': particle with id ',particles(n)%id,' will be deleted!' 1150 CALL message( 'lpm_check_cfl', 'PA0500', 0, 1, 0, 6, 0 ) 1151 particles(n)%particle_mask= .FALSE. 1152 ENDIF 1153 ENDDO 1154 ENDDO 1155 ENDDO 1156 ENDDO 1157 1158 END SUBROUTINE lpm_check_cfl 1159 1095 1160 !------------------------------------------------------------------------------! 1096 1161 ! Description: -
palm/trunk/SOURCE/lpm_init.f90
r2608 r2628 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Enabled particle advection with grid stretching. 28 ! 29 ! 2608 2017-11-13 14:04:26Z schwenkel 27 30 ! Calculation of magnus equation in external module (diagnostic_quantities_mod). 28 31 ! … … 624 627 !------------------------------------------------------------------------------! 625 628 SUBROUTINE lpm_create_particle (phase) 626 629 630 USE arrays_3d, & 631 ONLY: zw 627 632 USE lpm_exchange_horiz_mod, & 628 633 ONLY: lpm_exchange_horiz, lpm_move_particle, realloc_particles_array 629 634 630 USE lpm_pack_and_sort_mod, 635 USE lpm_pack_and_sort_mod, & 631 636 ONLY: lpm_sort_in_subboxes 632 637 … … 742 747 ip = tmp_particle%x * ddx 743 748 jp = tmp_particle%y * ddy 744 kp = tmp_particle%z / dz + 1 + offset_ocean_nzt 749 kp = tmp_particle%z / dz + 1 + offset_ocean_nzt 750 DO WHILE( zw(kp) < tmp_particle%z ) 751 kp = kp + 1 752 ENDDO 753 DO WHILE( zw(kp-1) > tmp_particle%z ) 754 kp = kp - 1 755 ENDDO 745 756 ! 746 757 !-- Determine surface level. Therefore, check for … … 929 940 i = particles(n)%x * ddx 930 941 j = particles(n)%y * ddy 931 k = particles(n)%z / dz + 1 + offset_ocean_nzt 942 k = particles(n)%z / dz + 1 + offset_ocean_nzt 943 DO WHILE( zw(k) < particles(n)%z ) 944 k = k + 1 945 ENDDO 946 DO WHILE( zw(k-1) > particles(n)%z ) 947 k = k - 1 948 ENDDO 932 949 ! 933 950 !-- Check if particle is within topography -
palm/trunk/SOURCE/lpm_pack_arrays.f90
r2609 r2628 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Enabled particle advection with grid stretching. 28 ! 29 ! 2609 2017-11-14 14:14:44Z schwenkel 27 30 ! Integrated subroutine pack_and_sort into lpm_sort_in_subboxes 28 31 ! … … 112 115 113 116 USE cpulog, & 114 117 ONLY: cpu_log, log_point_s 115 118 116 119 USE indices, & 117 120 ONLY: nxl, nxr, nys, nyn, nzb, nzt 118 121 119 122 USE kinds 120 123 121 124 USE control_parameters, & 122 125 ONLY: dz 123 126 124 127 USE grid_variables, & 125 ONLY: dx,dy,ddx, ddy 126 128 ONLY: dx,dy,ddx, ddy 129 130 USE arrays_3d, & 131 ONLY: zu 127 132 IMPLICIT NONE 128 133 … … 167 172 i = (particles(n)%x + 0.5_wp * dx) * ddx 168 173 j = (particles(n)%y + 0.5_wp * dy) * ddy 169 k = ( particles(n)%z+ 0.5_wp *dz ) / dz+1 + offset_ocean_nzt170 174 171 175 IF ( i == ip ) sort_index = sort_index + 4 172 IF ( j == jp ) sort_index = sort_index + 2 173 IF ( k == kp )sort_index = sort_index + 1176 IF ( j == jp ) sort_index = sort_index + 2 177 IF ( zu(kp) > particles(n)%z ) sort_index = sort_index + 1 174 178 175 179 sort_count(sort_index) = sort_count(sort_index) + 1 … … 343 347 344 348 345 END modulelpm_pack_and_sort_mod349 END MODULE lpm_pack_and_sort_mod
Note: See TracChangeset
for help on using the changeset viewer.
