Changeset 318 for palm/trunk/SOURCE/timestep.f90

Timestamp:

May 15, 2009 4:25:33 PM (16 years ago)

Author:

letzel

Message:

• Modified timestep criterion in case of simulations with plant canopy (timestep; modules rolled back to revision 305)

File:

• palm/trunk/SOURCE/timestep.f90 (modified) (2 diffs)

palm/trunk/SOURCE/timestep.f90

@@ -47 +47 @@
     INTEGER ::  i, j, k
 
-    REAL ::  div, dt_diff, dt_diff_l, dt_plant_canopy,  & 
-             dt_u, dt_v, dt_w, lad_max, percent_change, &
+    REAL ::  div, dt_diff, dt_diff_l, dt_plant_canopy,                 &
+             dt_plant_canopy_l,                                        &
+             dt_plant_canopy_u, dt_plant_canopy_v, dt_plant_canopy_w,  & 
+             dt_u, dt_v, dt_w, lad_max, percent_change,                &
              u_gtrans_l, vabs_max, value, v_gtrans_l
 
@@ -182 +184 @@
 !--    it is not allowed to extract more than the available momentum
        IF ( plant_canopy ) THEN
-!
-!--       Determine the maxima of the leaf drag coefficient and the leaf 
-!--       area density.
-          CALL global_min_max( nzb, nzt+1, nys-1, nyn+1, nxl-1, nxr+1, cdc, &
-                               'max', cdc_max, cdc_max_ijk )
-          CALL global_min_max( nzb, nzt+1, nys-1, nyn+1, nxl-1, nxr+1, lad_u, &
-                               'max', lad_u_max, lad_u_max_ijk )
-          CALL global_min_max( nzb, nzt+1, nys-1, nyn+1, nxl-1, nxr+1, lad_v, &
-                               'max', lad_v_max, lad_v_max_ijk )
-          CALL global_min_max( nzb, nzt+1, nys-1, nyn+1, nxl-1, nxr+1, lad_w, &
-                               'max', lad_w_max, lad_w_max_ijk )
-!
-!--       Determine the potential maximum of the absolute value of the velocity
-          vabs_max = SQRT( u_max**2 + v_max**2 + w_max**2 )
-          lad_max  = MAX( lad_u_max, lad_v_max, lad_w_max )
-
-          dt_plant_canopy = 1.0 / ( cdc_max * lad_max * vabs_max )
+
+          dt_plant_canopy_l = 0.0
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                DO k = nzb+1, nzt
+                   dt_plant_canopy_u = cdc(k,j,i) * lad_u(k,j,i) *  &
+                                       SQRT(     u(k,j,i)**2     +  &
+                                             ( ( v(k,j,i-1)      +  &
+                                                 v(k,j,i)        +  &
+                                                 v(k,j+1,i)      +  &
+                                                 v(k,j+1,i-1) )     &
+                                               / 4.0 )**2        +  &
+                                             ( ( w(k-1,j,i-1)    +  &
+                                                 w(k-1,j,i)      +  &
+                                                 w(k,j,i-1)      +  &
+                                                 w(k,j,i) )         &
+                                                 / 4.0 )**2 ) 
+                   IF ( dt_plant_canopy_u > dt_plant_canopy_l ) THEN
+                      dt_plant_canopy_l = dt_plant_canopy_u  
+                   ENDIF
+                   dt_plant_canopy_v = cdc(k,j,i) * lad_v(k,j,i) *  &
+                                       SQRT( ( ( u(k,j-1,i)      +  &
+                                                 u(k,j-1,i+1)    +  &
+                                                 u(k,j,i)        +  &
+                                                 u(k,j,i+1) )       &
+                                               / 4.0 )**2        +  &
+                                                 v(k,j,i)**2     +  &
+                                             ( ( w(k-1,j-1,i)    +  &
+                                                 w(k-1,j,i)      +  &
+                                                 w(k,j-1,i)      +  &
+                                                 w(k,j,i) )         &
+                                                 / 4.0 )**2 ) 
+                   IF ( dt_plant_canopy_v > dt_plant_canopy_l ) THEN
+                      dt_plant_canopy_l = dt_plant_canopy_v
+                   ENDIF                   
+                   dt_plant_canopy_w = cdc(k,j,i) * lad_w(k,j,i) *  &
+                                       SQRT( ( ( u(k,j,i)        +  &
+                                                 u(k,j,i+1)      +  &
+                                                 u(k+1,j,i)      +  &
+                                                 u(k+1,j,i+1) )     &
+                                               / 4.0 )**2        +  &
+                                             ( ( v(k,j,i)        +  &
+                                                 v(k,j+1,i)      +  &
+                                                 v(k+1,j,i)      +  &
+                                                 v(k+1,j+1,i) )     &
+                                               / 4.0 )**2        +  &
+                                                 w(k,j,i)**2 )     
+                   IF ( dt_plant_canopy_w > dt_plant_canopy_l ) THEN
+                      dt_plant_canopy_l = dt_plant_canopy_w
+                   ENDIF
+                ENDDO
+             ENDDO
+          ENDDO 
+
+          IF ( dt_plant_canopy_l > 0.0 ) THEN
+             dt_plant_canopy_l = 0.1 / dt_plant_canopy_l
+          ENDIF
+!
+!--       Determine the global minumum
+#if defined( __parallel )
+          CALL MPI_ALLREDUCE( dt_plant_canopy_l, dt_plant_canopy, 1, MPI_REAL,  &
+                              MPI_MIN, comm2d, ierr )
+#else
+          dt_plant_canopy = dt_plant_canopy_l
+#endif
 
        ELSE