Changeset 4275 for palm/trunk/SOURCE/lagrangian_particle_model_mod.f90
 Timestamp:
 Oct 28, 2019 3:34:55 PM (20 months ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SOURCE/lagrangian_particle_model_mod.f90
r4232 r4275 25 25 !  26 26 ! $Id$ 27 ! Change call of simple predictor corrector method, i.e. two divergence free 28 ! velocitiy fields are now used. 29 ! 30 ! 4232 20190920 09:34:22Z knoop 27 31 ! Removed INCLUDE "mpif.h", as it is not needed because of USE pegrid 28 32 ! … … 104 108 ONLY: de_dx, de_dy, de_dz, dzw, zu, zw, ql_c, ql_v, ql_vp, hyp, & 105 109 pt, q, exner, ql, diss, e, u, v, w, km, ql_1, ql_2, pt_p, q_p, & 106 d_exner , u_p, v_p, w_p110 d_exner 107 111 108 112 USE averaging, & … … 291 295 292 296 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ckernel !< 297 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: u_t !< u value of old timelevel t 298 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: v_t !< v value of old timelevel t 299 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: w_t !< w value of old timelevel t 300 293 301 294 302 INTEGER(iwp), PARAMETER :: PHASE_INIT = 1 !< … … 769 777 message_string = 'unknown collision kernel: collision_kernel = "' // & 770 778 TRIM( collision_kernel ) // '"' 771 CALL message( ' check_parameters', 'PA0350', 1, 2, 0, 6, 0 )779 CALL message( 'lpm_check_parameters', 'PA0350', 1, 2, 0, 6, 0 ) 772 780 773 781 END SELECT … … 784 792 'simple interpolation method is not ' // & 785 793 'implemented' 786 CALL message( ' check_parameters', 'PA0659', 1, 2, 0, 6, 0 )794 CALL message( 'lpm_check_parameters', 'PA0659', 1, 2, 0, 6, 0 ) 787 795 ENDIF 788 796 789 END SUBROUTINE 797 IF ( nested_run .AND. cloud_droplets ) THEN 798 message_string = 'nested runs in combination with cloud droplets ' // & 799 'is not implemented' 800 CALL message( 'lpm_check_parameters', 'PA0687', 1, 2, 0, 6, 0 ) 801 ENDIF 802 803 804 END SUBROUTINE lpm_check_parameters 790 805 791 806 !! … … 806 821 ql_vp(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 807 822 ENDIF 808 809 ! 810 ! Initial assignment of the pointers 823 824 825 ALLOCATE( u_t(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 826 v_t(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 827 w_t(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 828 ! 829 ! Initialize values with current time step 830 u_t = u 831 v_t = v 832 w_t = w 833 ! 834 ! Initial assignment of the pointers 811 835 IF ( cloud_droplets ) THEN 812 836 ql => ql_1 813 837 ql_c => ql_2 814 838 ENDIF 815 839 816 840 END SUBROUTINE lpm_init_arrays 817 841 … … 1974 1998 SELECT CASE ( location ) 1975 1999 1976 CASE ( 'after_pr ognostic_equations' )2000 CASE ( 'after_pressure_solver' ) 1977 2001 1978 2002 CALL cpu_log( log_point(25), 'lpm', 'start' ) … … 2250 2274 ! ENDIF 2251 2275 ! ENDIF 2276 2277 CASE ( 'after_integration' ) 2278 ! 2279 ! Call at the end of timestep routine to save particle velocities fields 2280 ! for the next timestep 2281 CALL lpm_swap_timelevel_for_particle_advection 2252 2282 2253 2283 CASE DEFAULT … … 3260 3290 ! 3261 3291 ! This case uses a simple interpolation method for the particle velocites, 3262 ! and applying a predictorcorrector method. @attention: for the corrector 3263 ! step the velocities of t(n+1) are required. However, at this moment of 3264 ! the time integration they are not free of divergence. This interpolation 3265 ! method is described in more detail in Grabowski et al., 2018 (GMD). 3292 ! and applying a predictorcorrector method. @note the current time divergence 3293 ! free time step is denoted with u_t etc.; the velocities of the time level of 3294 ! t+1 wit u,v, and w, as the model is called after swap timelevel 3295 ! @attention: for the corrector step the velocities of t(n+1) are required. 3296 ! Therefore the particle code is executed at the end of the time intermediate 3297 ! timestep routine. This interpolation method is described in more detail 3298 ! in Grabowski et al., 2018 (GMD). 3266 3299 IF ( interpolation_simple_corrector ) THEN 3267 3300 ! … … 3271 3304 3272 3305 alpha = MAX( MIN( ( particles(n)%x  ip * dx ) * ddx, 1.0_wp ), 0.0_wp ) 3273 u_int(n) = u (kp,jp,ip) * ( 1.0_wp  alpha ) + u(kp,jp,ip+1) * alpha3306 u_int(n) = u_t(kp,jp,ip) * ( 1.0_wp  alpha ) + u_t(kp,jp,ip+1) * alpha 3274 3307 3275 3308 beta = MAX( MIN( ( particles(n)%y  jp * dy ) * ddy, 1.0_wp ), 0.0_wp ) 3276 v_int(n) = v (kp,jp,ip) * ( 1.0_wp  beta ) + v(kp,jp+1,ip) * beta3309 v_int(n) = v_t(kp,jp,ip) * ( 1.0_wp  beta ) + v_t(kp,jp+1,ip) * beta 3277 3310 3278 3311 gamma = MAX( MIN( ( particles(n)%z  zw(kkw) ) / & 3279 3312 ( zw(kkw+1)  zw(kkw) ), 1.0_wp ), 0.0_wp ) 3280 w_int(n) = w (kkw,jp,ip) * ( 1.0_wp  gamma ) + w(kkw+1,jp,ip) * gamma3313 w_int(n) = w_t(kkw,jp,ip) * ( 1.0_wp  gamma ) + w_t(kkw+1,jp,ip) * gamma 3281 3314 3282 3315 ENDDO … … 3309 3342 ! 3310 3343 ! Calculate part of the corrector step 3311 unext = u _p(k_next+1, j_next, i_next) * ( 1.0_wp  alpha ) + &3312 u _p(k_next+1, j_next, i_next+1) * alpha3313 3314 vnext = v _p(k_next+1, j_next, i_next) * ( 1.0_wp  beta ) + &3315 v _p(k_next+1, j_next+1, i_next ) * beta3316 3317 wnext = w _p(k_next, j_next, i_next) * ( 1.0_wp  gamma ) + &3318 w _p(k_next+1, j_next, i_next ) * gamma3344 unext = u(k_next+1, j_next, i_next) * ( 1.0_wp  alpha ) + & 3345 u(k_next+1, j_next, i_next+1) * alpha 3346 3347 vnext = v(k_next+1, j_next, i_next) * ( 1.0_wp  beta ) + & 3348 v(k_next+1, j_next+1, i_next ) * beta 3349 3350 wnext = w(k_next, j_next, i_next) * ( 1.0_wp  gamma ) + & 3351 w(k_next+1, j_next, i_next ) * gamma 3319 3352 3320 3353 ! … … 4280 4313 v_sgs = v_sgs * fac + term1 + term2 + term3 4281 4314 4282 END SUBROUTINE weil_stochastic_eq 4283 4284 4315 END SUBROUTINE weil_stochastic_eq 4316 4317 4318 !! 4319 ! Description: 4320 !  4321 !> swap timelevel in case of particle advection interpolation 'simplecorrector' 4322 !> This routine is called at the end of one timestep, the velocities are then 4323 !> used for the next timestep 4324 !! 4325 SUBROUTINE lpm_swap_timelevel_for_particle_advection 4326 4327 ! 4328 ! save the divergence free velocites of t+1 to use them at the end of the 4329 ! next time step 4330 u_t = u 4331 v_t = v 4332 w_t = w 4333 4334 END SUBROUTINE lpm_swap_timelevel_for_particle_advection 4335 4336 4285 4337 !! 4286 4338 ! Description: … … 5114 5166 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 5115 5167 5116 q _p(k,j,i) = q_p(k,j,i)  ql_c(k,j,i) * flag5117 pt _p(k,j,i) = pt_p(k,j,i) + lv_d_cp * ql_c(k,j,i) * d_exner(k) &5118 * flag5168 q(k,j,i) = q_p(k,j,i)  ql_c(k,j,i) * flag 5169 pt(k,j,i) = pt(k,j,i) + lv_d_cp * ql_c(k,j,i) * d_exner(k) & 5170 * flag 5119 5171 ENDDO 5120 5172 ENDDO … … 5144 5196 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 5145 5197 5146 q _p(k,j,i) = q_p(k,j,i)  ql_c(k,j,i) * flag5147 pt _p(k,j,i) = pt_p(k,j,i) + lv_d_cp * ql_c(k,j,i) * d_exner(k) * flag5198 q(k,j,i) = q(k,j,i)  ql_c(k,j,i) * flag 5199 pt(k,j,i) = pt(k,j,i) + lv_d_cp * ql_c(k,j,i) * d_exner(k) * flag 5148 5200 ENDDO 5149 5201
Note: See TracChangeset
for help on using the changeset viewer.