Ignore:
Timestamp:
May 15, 2009 4:25:33 PM (15 years ago)
Author:
letzel
Message:
  • Modified timestep criterion in case of simulations with plant canopy (timestep; modules rolled back to revision 305)
File:
1 edited

Legend:

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

    r316 r318  
    4747    INTEGER ::  i, j, k
    4848
    49     REAL ::  div, dt_diff, dt_diff_l, dt_plant_canopy,  &
    50              dt_u, dt_v, dt_w, lad_max, percent_change, &
     49    REAL ::  div, dt_diff, dt_diff_l, dt_plant_canopy,                 &
     50             dt_plant_canopy_l,                                        &
     51             dt_plant_canopy_u, dt_plant_canopy_v, dt_plant_canopy_w,  &
     52             dt_u, dt_v, dt_w, lad_max, percent_change,                &
    5153             u_gtrans_l, vabs_max, value, v_gtrans_l
    5254
     
    182184!--    it is not allowed to extract more than the available momentum
    183185       IF ( plant_canopy ) THEN
    184 !
    185 !--       Determine the maxima of the leaf drag coefficient and the leaf
    186 !--       area density.
    187           CALL global_min_max( nzb, nzt+1, nys-1, nyn+1, nxl-1, nxr+1, cdc, &
    188                                'max', cdc_max, cdc_max_ijk )
    189           CALL global_min_max( nzb, nzt+1, nys-1, nyn+1, nxl-1, nxr+1, lad_u, &
    190                                'max', lad_u_max, lad_u_max_ijk )
    191           CALL global_min_max( nzb, nzt+1, nys-1, nyn+1, nxl-1, nxr+1, lad_v, &
    192                                'max', lad_v_max, lad_v_max_ijk )
    193           CALL global_min_max( nzb, nzt+1, nys-1, nyn+1, nxl-1, nxr+1, lad_w, &
    194                                'max', lad_w_max, lad_w_max_ijk )
    195 !
    196 !--       Determine the potential maximum of the absolute value of the velocity
    197           vabs_max = SQRT( u_max**2 + v_max**2 + w_max**2 )
    198           lad_max  = MAX( lad_u_max, lad_v_max, lad_w_max )
    199 
    200           dt_plant_canopy = 1.0 / ( cdc_max * lad_max * vabs_max )
     186
     187          dt_plant_canopy_l = 0.0
     188          DO  i = nxl, nxr
     189             DO  j = nys, nyn
     190                DO k = nzb+1, nzt
     191                   dt_plant_canopy_u = cdc(k,j,i) * lad_u(k,j,i) *  &
     192                                       SQRT(     u(k,j,i)**2     +  &
     193                                             ( ( v(k,j,i-1)      +  &
     194                                                 v(k,j,i)        +  &
     195                                                 v(k,j+1,i)      +  &
     196                                                 v(k,j+1,i-1) )     &
     197                                               / 4.0 )**2        +  &
     198                                             ( ( w(k-1,j,i-1)    +  &
     199                                                 w(k-1,j,i)      +  &
     200                                                 w(k,j,i-1)      +  &
     201                                                 w(k,j,i) )         &
     202                                                 / 4.0 )**2 )
     203                   IF ( dt_plant_canopy_u > dt_plant_canopy_l ) THEN
     204                      dt_plant_canopy_l = dt_plant_canopy_u 
     205                   ENDIF
     206                   dt_plant_canopy_v = cdc(k,j,i) * lad_v(k,j,i) *  &
     207                                       SQRT( ( ( u(k,j-1,i)      +  &
     208                                                 u(k,j-1,i+1)    +  &
     209                                                 u(k,j,i)        +  &
     210                                                 u(k,j,i+1) )       &
     211                                               / 4.0 )**2        +  &
     212                                                 v(k,j,i)**2     +  &
     213                                             ( ( w(k-1,j-1,i)    +  &
     214                                                 w(k-1,j,i)      +  &
     215                                                 w(k,j-1,i)      +  &
     216                                                 w(k,j,i) )         &
     217                                                 / 4.0 )**2 )
     218                   IF ( dt_plant_canopy_v > dt_plant_canopy_l ) THEN
     219                      dt_plant_canopy_l = dt_plant_canopy_v
     220                   ENDIF                   
     221                   dt_plant_canopy_w = cdc(k,j,i) * lad_w(k,j,i) *  &
     222                                       SQRT( ( ( u(k,j,i)        +  &
     223                                                 u(k,j,i+1)      +  &
     224                                                 u(k+1,j,i)      +  &
     225                                                 u(k+1,j,i+1) )     &
     226                                               / 4.0 )**2        +  &
     227                                             ( ( v(k,j,i)        +  &
     228                                                 v(k,j+1,i)      +  &
     229                                                 v(k+1,j,i)      +  &
     230                                                 v(k+1,j+1,i) )     &
     231                                               / 4.0 )**2        +  &
     232                                                 w(k,j,i)**2 )     
     233                   IF ( dt_plant_canopy_w > dt_plant_canopy_l ) THEN
     234                      dt_plant_canopy_l = dt_plant_canopy_w
     235                   ENDIF
     236                ENDDO
     237             ENDDO
     238          ENDDO
     239
     240          IF ( dt_plant_canopy_l > 0.0 ) THEN
     241             dt_plant_canopy_l = 0.1 / dt_plant_canopy_l
     242          ENDIF
     243!
     244!--       Determine the global minumum
     245#if defined( __parallel )
     246          CALL MPI_ALLREDUCE( dt_plant_canopy_l, dt_plant_canopy, 1, MPI_REAL,  &
     247                              MPI_MIN, comm2d, ierr )
     248#else
     249          dt_plant_canopy = dt_plant_canopy_l
     250#endif
    201251
    202252       ELSE
Note: See TracChangeset for help on using the changeset viewer.