Changeset 4275 for palm

Timestamp:

Oct 28, 2019 3:34:55 PM (5 years ago)

Author:

schwenkel

Message:

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

Location:

palm/trunk/SOURCE

Files:

• lagrangian_particle_model_mod.f90 (modified) (14 diffs)
• time_integration.f90 (modified) (4 diffs)

palm/trunk/SOURCE/lagrangian_particle_model_mod.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Change call of simple predictor corrector method, i.e. two divergence free
+! velocitiy fields are now used.
+!
+! 4232 2019-09-20 09:34:22Z knoop
 ! Removed INCLUDE "mpif.h", as it is not needed because of USE pegrid
 ! 
@@ -104 +108 @@
         ONLY:  de_dx, de_dy, de_dz, dzw, zu, zw,  ql_c, ql_v, ql_vp, hyp,      &
                pt, q, exner, ql, diss, e, u, v, w, km, ql_1, ql_2, pt_p, q_p,  &
-               d_exner, u_p, v_p, w_p
+               d_exner
  
     USE averaging,                                                             &
@@ -291 +295 @@
 
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ckernel !<
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  u_t   !< u value of old timelevel t
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  v_t   !< v value of old timelevel t
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  w_t   !< w value of old timelevel t
+
 
     INTEGER(iwp), PARAMETER         ::  PHASE_INIT    = 1  !<
@@ -769 +777 @@
           message_string = 'unknown collision kernel: collision_kernel = "' // &
                            TRIM( collision_kernel ) // '"'
-          CALL message( 'check_parameters', 'PA0350', 1, 2, 0, 6, 0 )
+          CALL message( 'lpm_check_parameters', 'PA0350', 1, 2, 0, 6, 0 )
 
     END SELECT
@@ -784 +792 @@
                            'simple interpolation method is not '            // &
                            'implemented'
-          CALL message( 'check_parameters', 'PA0659', 1, 2, 0, 6, 0 )
+          CALL message( 'lpm_check_parameters', 'PA0659', 1, 2, 0, 6, 0 )
     ENDIF
 
- END SUBROUTINE
+    IF ( nested_run  .AND.  cloud_droplets )  THEN
+       message_string = 'nested runs in combination with cloud droplets ' // &
+                        'is not implemented'
+          CALL message( 'lpm_check_parameters', 'PA0687', 1, 2, 0, 6, 0 )
+    ENDIF
+
+
+ END SUBROUTINE lpm_check_parameters
  
 !------------------------------------------------------------------------------!
@@ -806 +821 @@
                      ql_vp(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     ENDIF
-    
-!
-!--    Initial assignment of the pointers    
+
+
+    ALLOCATE( u_t(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
+              v_t(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
+              w_t(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+!
+!-- Initialize values with current time step
+    u_t = u
+    v_t = v
+    w_t = w
+!
+!--    Initial assignment of the pointers
     IF ( cloud_droplets )  THEN
        ql   => ql_1
        ql_c => ql_2
     ENDIF
-    
+
  END SUBROUTINE lpm_init_arrays
  
@@ -1974 +1998 @@
     SELECT CASE ( location )
 
-       CASE ( 'after_prognostic_equations' )
+       CASE ( 'after_pressure_solver' )
 
           CALL cpu_log( log_point(25), 'lpm', 'start' )
@@ -2250 +2274 @@
 !              ENDIF
 !           ENDIF
+
+       CASE ( 'after_integration' )
+!
+!--       Call at the end of timestep routine to save particle velocities fields
+!--       for the next timestep
+          CALL lpm_swap_timelevel_for_particle_advection
 
        CASE DEFAULT
@@ -3260 +3290 @@
 !
 !-- This case uses a simple interpolation method for the particle velocites,
-!-- and applying a predictor-corrector method. @attention: for the corrector
-!-- step the velocities of t(n+1) are required. However, at this moment of
-!-- the time integration they are not free of divergence. This interpolation
-!-- method is described in more detail in Grabowski et al., 2018 (GMD).
+!-- and applying a predictor-corrector method. @note the current time divergence
+!-- free time step is denoted with u_t etc.; the velocities of the time level of
+!-- t+1 wit u,v, and w, as the model is called after swap timelevel
+!-- @attention: for the corrector step the velocities of t(n+1) are required.
+!-- Therefore the particle code is executed at the end of the time intermediate
+!-- timestep routine. This interpolation method is described in more detail
+!-- in Grabowski et al., 2018 (GMD).
     IF ( interpolation_simple_corrector )  THEN
 !
@@ -3271 +3304 @@
 
           alpha = MAX( MIN( ( particles(n)%x - ip * dx ) * ddx, 1.0_wp ), 0.0_wp )
-          u_int(n) = u(kp,jp,ip) * ( 1.0_wp - alpha ) + u(kp,jp,ip+1) * alpha
+          u_int(n) = u_t(kp,jp,ip) * ( 1.0_wp - alpha ) + u_t(kp,jp,ip+1) * alpha
 
           beta  = MAX( MIN( ( particles(n)%y - jp * dy ) * ddy, 1.0_wp ), 0.0_wp )
-          v_int(n) = v(kp,jp,ip) * ( 1.0_wp - beta ) + v(kp,jp+1,ip) * beta
+          v_int(n) = v_t(kp,jp,ip) * ( 1.0_wp - beta ) + v_t(kp,jp+1,ip) * beta
 
           gamma = MAX( MIN( ( particles(n)%z - zw(kkw) ) /                   &
                             ( zw(kkw+1) - zw(kkw) ), 1.0_wp ), 0.0_wp )
-          w_int(n) = w(kkw,jp,ip) * ( 1.0_wp - gamma ) + w(kkw+1,jp,ip) * gamma
+          w_int(n) = w_t(kkw,jp,ip) * ( 1.0_wp - gamma ) + w_t(kkw+1,jp,ip) * gamma
 
        ENDDO
@@ -3309 +3342 @@
 !
 !--          Calculate part of the corrector step
-             unext = u_p(k_next+1, j_next, i_next) * ( 1.0_wp - alpha ) +    &
-                     u_p(k_next+1, j_next,   i_next+1) * alpha
-
-             vnext = v_p(k_next+1, j_next, i_next) * ( 1.0_wp - beta  ) +    &
-                     v_p(k_next+1, j_next+1, i_next  ) * beta
-
-             wnext = w_p(k_next,   j_next, i_next) * ( 1.0_wp - gamma ) +    &
-                     w_p(k_next+1, j_next, i_next  ) * gamma
+             unext = u(k_next+1, j_next, i_next) * ( 1.0_wp - alpha ) +    &
+                     u(k_next+1, j_next,   i_next+1) * alpha
+
+             vnext = v(k_next+1, j_next, i_next) * ( 1.0_wp - beta  ) +    &
+                     v(k_next+1, j_next+1, i_next  ) * beta
+
+             wnext = w(k_next,   j_next, i_next) * ( 1.0_wp - gamma ) +    &
+                     w(k_next+1, j_next, i_next  ) * gamma
 
 !
@@ -4280 +4313 @@
     v_sgs = v_sgs * fac + term1 + term2 + term3
 
- END SUBROUTINE weil_stochastic_eq 
- 
- 
+ END SUBROUTINE weil_stochastic_eq
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> swap timelevel in case of particle advection interpolation 'simple-corrector'
+!> This routine is called at the end of one timestep, the velocities are then
+!> used for the next timestep
+!------------------------------------------------------------------------------!
+ SUBROUTINE lpm_swap_timelevel_for_particle_advection
+
+!
+!-- save the divergence free velocites of t+1 to use them at the end of the
+!-- next time step
+    u_t = u
+    v_t = v
+    w_t = w
+
+ END SUBROUTINE lpm_swap_timelevel_for_particle_advection
+
+
 !------------------------------------------------------------------------------!  
 ! Description:
@@ -5114 +5166 @@
              flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
 
-             q_p(k,j,i)  = q_p(k,j,i)  - ql_c(k,j,i) * flag
-             pt_p(k,j,i) = pt_p(k,j,i) + lv_d_cp * ql_c(k,j,i) * d_exner(k) &
-                                                     * flag
+             q(k,j,i)  = q_p(k,j,i)  - ql_c(k,j,i) * flag
+             pt(k,j,i) = pt(k,j,i) + lv_d_cp * ql_c(k,j,i) * d_exner(k) &
+                                                           * flag
           ENDDO
        ENDDO
@@ -5144 +5196 @@
        flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
 
-       q_p(k,j,i)  = q_p(k,j,i)  - ql_c(k,j,i) * flag
-       pt_p(k,j,i) = pt_p(k,j,i) + lv_d_cp * ql_c(k,j,i) * d_exner(k) * flag
+       q(k,j,i)  = q(k,j,i)  - ql_c(k,j,i) * flag
+       pt(k,j,i) = pt(k,j,i) + lv_d_cp * ql_c(k,j,i) * d_exner(k) * flag
     ENDDO
 

palm/trunk/SOURCE/time_integration.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Move call oft lpm to the end of intermediate timestep loop
+! 
+! 4268 2019-10-17 11:29:38Z schwenkel
 ! Removing module specific boundary conditions an put them into their modules
 ! 
@@ -656 +659 @@
              CALL prognostic_equations_vector
           ENDIF
-
-!
-!--       Particle transport/physics with the Lagrangian particle model
-!--       (only once during intermediate steps, because it uses an Euler-step)
-!--       ### particle model should be moved before prognostic_equations, in order
-!--       to regard droplet interactions directly
-          IF ( particle_advection  .AND.  time_since_reference_point >= particle_advection_start   &
-               .AND.  intermediate_timestep_count == 1 )                                           &
-          THEN
-             CALL module_interface_actions( 'after_prognostic_equations' )
-             first_call_lpm = .FALSE.
-          ENDIF
-
-!
-!--       Interaction of droplets with temperature and mixing ratio.
-!--       Droplet condensation and evaporation is calculated within
-!--       advec_particles.
-          IF ( cloud_droplets  .AND.  intermediate_timestep_count == intermediate_timestep_count_max ) &
-          THEN
-             CALL lpm_interaction_droplets_ptq
-          ENDIF
-
 !
 !--       Movement of agents in multi agent system
@@ -720 +701 @@
                 CALL exchange_horiz( nr_p, nbgp )
              ENDIF
-          ENDIF
-          IF ( cloud_droplets )  THEN
-             CALL exchange_horiz( ql, nbgp )
-             CALL exchange_horiz( ql_c, nbgp )
-             CALL exchange_horiz( ql_v, nbgp )
-             CALL exchange_horiz( ql_vp, nbgp )
           ENDIF
           IF ( passive_scalar )  CALL exchange_horiz( s_p, nbgp )
@@ -995 +970 @@
 
           ENDIF
-
+!
+!--       Particle transport/physics with the Lagrangian particle model
+!--       (only once during intermediate steps, because it uses an Euler-step)
+!--       ### particle model should be moved before prognostic_equations, in order
+!--       to regard droplet interactions directly
+          IF ( particle_advection  .AND.  time_since_reference_point >= particle_advection_start   &
+               .AND.  intermediate_timestep_count == intermediate_timestep_count_max )             &
+          THEN
+             CALL module_interface_actions( 'after_pressure_solver' )
+             first_call_lpm = .FALSE.
+          ENDIF
+
+          IF ( cloud_droplets )  THEN
+             CALL exchange_horiz( ql, nbgp )
+             CALL exchange_horiz( ql_c, nbgp )
+             CALL exchange_horiz( ql_v, nbgp )
+             CALL exchange_horiz( ql_vp, nbgp )
+          ENDIF
+!
+!--       Interaction of droplets with temperature and mixing ratio.
+!--       Droplet condensation and evaporation is calculated within
+!--       advec_particles.
+          IF ( cloud_droplets  .AND.  intermediate_timestep_count == intermediate_timestep_count_max ) &
+          THEN
+             CALL lpm_interaction_droplets_ptq
+             CALL exchange_horiz( pt, nbgp )
+             CALL exchange_horiz( q, nbgp )
+          ENDIF
 !
 !--       If required, compute liquid water content