Ignore:
Timestamp:
Oct 28, 2019 3:34:55 PM (20 months ago)
Author:
schwenkel
Message:

move call of lpm at the end of intermediate timeloop and improve simple corrector particle interpolation method

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/lagrangian_particle_model_mod.f90

    r4232 r4275  
    2525! -----------------
    2626! $Id$
     27! Change call of simple predictor corrector method, i.e. two divergence free
     28! velocitiy fields are now used.
     29!
     30! 4232 2019-09-20 09:34:22Z knoop
    2731! Removed INCLUDE "mpif.h", as it is not needed because of USE pegrid
    2832!
     
    104108        ONLY:  de_dx, de_dy, de_dz, dzw, zu, zw,  ql_c, ql_v, ql_vp, hyp,      &
    105109               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_p
     110               d_exner
    107111 
    108112    USE averaging,                                                             &
     
    291295
    292296    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
    293301
    294302    INTEGER(iwp), PARAMETER         ::  PHASE_INIT    = 1  !<
     
    769777          message_string = 'unknown collision kernel: collision_kernel = "' // &
    770778                           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 )
    772780
    773781    END SELECT
     
    784792                           'simple interpolation method is not '            // &
    785793                           'implemented'
    786           CALL message( 'check_parameters', 'PA0659', 1, 2, 0, 6, 0 )
     794          CALL message( 'lpm_check_parameters', 'PA0659', 1, 2, 0, 6, 0 )
    787795    ENDIF
    788796
    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
    790805 
    791806!------------------------------------------------------------------------------!
     
    806821                     ql_vp(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    807822    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
    811835    IF ( cloud_droplets )  THEN
    812836       ql   => ql_1
    813837       ql_c => ql_2
    814838    ENDIF
    815    
     839
    816840 END SUBROUTINE lpm_init_arrays
    817841 
     
    19741998    SELECT CASE ( location )
    19751999
    1976        CASE ( 'after_prognostic_equations' )
     2000       CASE ( 'after_pressure_solver' )
    19772001
    19782002          CALL cpu_log( log_point(25), 'lpm', 'start' )
     
    22502274!              ENDIF
    22512275!           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
    22522282
    22532283       CASE DEFAULT
     
    32603290!
    32613291!-- This case uses a simple interpolation method for the particle velocites,
    3262 !-- and applying a predictor-corrector 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 predictor-corrector 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).
    32663299    IF ( interpolation_simple_corrector )  THEN
    32673300!
     
    32713304
    32723305          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) * alpha
     3306          u_int(n) = u_t(kp,jp,ip) * ( 1.0_wp - alpha ) + u_t(kp,jp,ip+1) * alpha
    32743307
    32753308          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) * beta
     3309          v_int(n) = v_t(kp,jp,ip) * ( 1.0_wp - beta ) + v_t(kp,jp+1,ip) * beta
    32773310
    32783311          gamma = MAX( MIN( ( particles(n)%z - zw(kkw) ) /                   &
    32793312                            ( 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) * gamma
     3313          w_int(n) = w_t(kkw,jp,ip) * ( 1.0_wp - gamma ) + w_t(kkw+1,jp,ip) * gamma
    32813314
    32823315       ENDDO
     
    33093342!
    33103343!--          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) * alpha
    3313 
    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  ) * beta
    3316 
    3317              wnext = w_p(k_next,   j_next, i_next) * ( 1.0_wp - gamma ) +    &
    3318                      w_p(k_next+1, j_next, i_next  ) * gamma
     3344             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
    33193352
    33203353!
     
    42804313    v_sgs = v_sgs * fac + term1 + term2 + term3
    42814314
    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 'simple-corrector'
     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
    42854337!------------------------------------------------------------------------------! 
    42864338! Description:
     
    51145166             flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    51155167
    5116              q_p(k,j,i)  = q_p(k,j,i)  - ql_c(k,j,i) * flag
    5117              pt_p(k,j,i) = pt_p(k,j,i) + lv_d_cp * ql_c(k,j,i) * d_exner(k) &
    5118                                                      * flag
     5168             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
    51195171          ENDDO
    51205172       ENDDO
     
    51445196       flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    51455197
    5146        q_p(k,j,i)  = q_p(k,j,i)  - ql_c(k,j,i) * flag
    5147        pt_p(k,j,i) = pt_p(k,j,i) + lv_d_cp * ql_c(k,j,i) * d_exner(k) * flag
     5198       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
    51485200    ENDDO
    51495201
Note: See TracChangeset for help on using the changeset viewer.