Changeset 4770


Ignore:
Timestamp:
Nov 3, 2020 2:04:53 PM (4 years ago)
Author:
suehring
Message:

Consider basal-area density as an additional sink for momentum in the prognostic equations for momentum and SGS TKE

File:
1 edited

Legend:

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

    r4768 r4770  
    2727! -----------------
    2828! $Id$
     29! Consider basal-area density as an additional sink for momentum in the
     30! prognostic equations for momentum and SGS TKE
     31!
     32! 4768 2020-11-02 19:11:23Z suehring
    2933! Enable 3D data output also with 64-bit precision
    3034!
     
    293297    REAL(wp), DIMENSION(:), ALLOCATABLE ::  lad            !< leaf area density
    294298    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pre_lad        !< preliminary lad
    295    
     299
     300    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  bad_s                    !< basal-area density on scalar-grid
    296301    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  cum_lai_hf               !< cumulative lai for heatflux calc.
    297302    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  lad_s                    !< lad on scalar-grid
     
    546551          unit = 'K s-1'
    547552
    548        CASE ( 'pcm_lad' )
     553       CASE ( 'pcm_bad', 'pcm_lad' )
    549554          unit = 'm2 m-3'
    550555
     
    870875          ENDIF
    871876
     877       CASE ( 'pcm_bad' )
     878          IF ( av == 0 )  THEN
     879             DO  i = nxl, nxr
     880                DO  j = nys, nyn
     881                   DO  k = MAX( 1, nzb_do ), MIN( pch_index_ji(j,i), nzt_do )
     882                      local_pf(i,j,k) = bad_s(k,j,i)
     883                   ENDDO
     884                ENDDO
     885             ENDDO
     886          ENDIF
     887
    872888       CASE DEFAULT
    873889          found = .FALSE.
     
    898914     SELECT CASE ( TRIM( var ) )
    899915
    900         CASE ( 'pcm_heatrate', 'pcm_lad', 'pcm_transpirationrate', 'pcm_latentrate')
     916        CASE ( 'pcm_heatrate', 'pcm_bad', 'pcm_lad', 'pcm_transpirationrate', 'pcm_latentrate')
    901917           grid_x = 'x'
    902918           grid_y = 'y'
     
    10391055
    10401056       LOGICAL      ::  lad_on_top = .FALSE.  !< dummy flag to indicate that LAD is defined on a building roof
     1057       LOGICAL      ::  bad_on_top = .FALSE.  !< dummy flag to indicate that BAD is defined on a building roof
    10411058
    10421059       REAL(wp) ::  canopy_height   !< canopy height for lad-profile construction
     
    11421159
    11431160!
    1144 !--    Allocate 3D-array for the leaf area density (lad_s).
     1161!--    Allocate 3D-array for the leaf-area density (lad_s) as well as
     1162!--    for basal-area densitiy (bad_s). Note, by default bad_s is zero.
    11451163       ALLOCATE( lad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    1146 
     1164       ALLOCATE( bad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     1165       bad_s = 0.0_wp
    11471166!
    11481167!--    Initialization of the canopy coverage in the model domain:
     
    13031322#endif
    13041323                IF ( lad_on_top )  THEN
    1305                    WRITE( message_string, * )                                  &
    1306                                         'Resolved plant-canopy is ' //         &
    1307                                         'defined on top of an artificially '// &
    1308                                         'created building grid point(s) ' //   &
    1309                                         '(emerged from the filtering) - ' //   &
    1310                                         'LAD profile is omitted at this / ' // &
     1324                   WRITE( message_string, * )                                                      &
     1325                                        'Resolved plant-canopy is defined on top of an ' //        &
     1326                                        'artificially created building grid point(s) '//           &
     1327                                        'the filtering) - LAD/BAD profile is omitted at this / ' //&
    13111328                                        'these grid point(s).'
    13121329                   CALL message( 'pcm_init', 'PA0313', 0, 0, 0, 6, 0 )
     
    13211338             ELSE
    13221339                CALL pcm_read_plant_canopy_3d
     1340             ENDIF
     1341!
     1342!--          Initialize LAD with data from file. If LAD is given in NetCDF file,
     1343!--          use these values, else take LAD profiles from ASCII file.
     1344!--          Please note, in NetCDF file LAD is only given up to the maximum
     1345!--          canopy top, indicated by basal_area_density_f%nz. 
     1346             bad_s = 0.0_wp
     1347             IF ( basal_area_density_f%from_file )  THEN
     1348!
     1349!--             Set also pch_index, used to be the upper bound of the vertical
     1350!--             loops. Therefore, use the global top of the canopy layer.
     1351                pch_index = basal_area_density_f%nz - 1
     1352
     1353                DO  i = nxl, nxr
     1354                   DO  j = nys, nyn
     1355                      DO  k = 0, basal_area_density_f%nz - 1
     1356                         IF ( basal_area_density_f%var(k,j,i) /= basal_area_density_f%fill )       &
     1357                            bad_s(k,j,i) = basal_area_density_f%var(k,j,i)
     1358                      ENDDO
     1359!
     1360!--                   Check if resolved vegetation is mapped onto buildings.
     1361!--                   Please see comment for leaf_area density
     1362                      IF ( ANY( bad_s(:,j,i) /= 0.0_wp )  .AND.                &
     1363                           ANY( BTEST( wall_flags_total_0(:,j,i), 6 ) ) .AND.  &
     1364                           ANY( BTEST( wall_flags_total_0(:,j,i), 4 ) ) )  THEN
     1365                         bad_s(:,j,i) = 0.0_wp
     1366                         bad_on_top   = .TRUE.
     1367                      ENDIF
     1368                   ENDDO
     1369                ENDDO
     1370#if defined( __parallel )
     1371               CALL MPI_ALLREDUCE( MPI_IN_PLACE, bad_on_top, 1, MPI_LOGICAL,  &
     1372                                   MPI_LOR, comm2d, ierr)
     1373#endif
     1374                IF ( bad_on_top )  THEN
     1375                   WRITE( message_string, * )                                                      &
     1376                                        'Resolved plant-canopy is defined on top of an ' //        &
     1377                                        'artificially created building grid point(s) '//           &
     1378                                        'the filtering) - LAD/BAD profile is omitted at this / ' //&
     1379                                        'these grid point(s).'
     1380                   CALL message( 'pcm_init', 'PA0313', 0, 0, 0, 6, 0 )
     1381                ENDIF
     1382                CALL exchange_horiz( bad_s, nbgp )
    13231383             ENDIF
    13241384
     
    13551415          DO  j = nysg, nyng
    13561416             DO  k = 0, pch_index
    1357                 IF ( lad_s(k,j,i) /= 0 )  pch_index_ji(j,i) = k
     1417                IF ( lad_s(k,j,i) /= 0.0_wp  .OR.  bad_s(k,j,i) /= 0.0_wp )  pch_index_ji(j,i) = k
    13581418             ENDDO
    13591419!
     
    18621922       LOGICAL ::  building_edge_w !< control flag indicating a westward-facing building edge
    18631923
     1924       REAL(wp) ::  bad_local !< local bad value
    18641925       REAL(wp) ::  ddt_3d    !< inverse of the LES timestep (dt_3d)
    18651926       REAL(wp) ::  lad_local !< local lad value
     
    19061967!--                   this is not applied, here the LAD is equals the LAD at grid
    19071968!--                   point (k,j,i), in order to avoid that LAD is mistakenly mapped
    1908 !--                   on top of a roof where (usually) is no LAD is defined.
     1969!--                   on top of a roof where (usually) is no LAD is defined. The same is
     1970!--                   also valid for bad_s.
    19091971                      lad_local = lad_s(kk,j,i)
    1910                       IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp  &
     1972                      IF ( lad_local == 0.0_wp  .AND.  lad_s(kk,j,i-1) > 0.0_wp                    &
    19111973                           .AND.  .NOT. building_edge_w )  lad_local = lad_s(kk,j,i-1)
     1974
     1975                      bad_local = bad_s(kk,j,i)
     1976                      IF ( bad_local == 0.0_wp  .AND.  bad_s(kk,j,i-1) > 0.0_wp                    &
     1977                           .AND.  .NOT. building_edge_w )  bad_local = bad_s(kk,j,i-1)
    19121978!
    19131979!--                   In order to avoid that LAD is mistakenly considered at right-
     
    19151981!--                   u-component at index j,i is still on the building while the
    19161982!--                   topography top for the scalar isn't), LAD is taken from grid
    1917 !--                   point (j,i-1).
    1918                       IF ( lad_local > 0.0_wp .AND. lad_s(kk,j,i-1) == 0.0_wp  &
     1983!--                   point (j,i-1). The same is also valid for bad_s.
     1984                      IF ( lad_local > 0.0_wp  .AND.  lad_s(kk,j,i-1) == 0.0_wp                    &
    19191985                           .AND.  building_edge_e )  lad_local = lad_s(kk,j,i-1)
     1986                      IF ( bad_local > 0.0_wp  .AND.  bad_s(kk,j,i-1) == 0.0_wp                    &
     1987                           .AND.  building_edge_e )  bad_local = bad_s(kk,j,i-1)
    19201988
    19211989                      pre_tend = 0.0_wp
     
    19241992!--                   Calculate preliminary value (pre_tend) of the tendency
    19251993                      pre_tend = - canopy_drag_coeff *                         &
    1926                                    lad_local *                                 &
     1994                                   ( lad_local + bad_local ) *                 &
    19271995                                   SQRT( u(k,j,i)**2 +                         &
    19281996                                         ( 0.25_wp * ( v(k,j,i-1) +            &
     
    19872055!--                   point (k,j,i), in order to avoid that LAD is mistakenly mapped
    19882056!--                   on top of a roof where (usually) is no LAD is defined.
     2057!--                   The same is also valid for bad_s.
    19892058                      lad_local = lad_s(kk,j,i)
    1990                       IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp &
     2059                      IF ( lad_local == 0.0_wp  .AND.  lad_s(kk,j-1,i) > 0.0_wp                    &
    19912060                       .AND.  .NOT. building_edge_s )  lad_local = lad_s(kk,j-1,i)
     2061
     2062                      bad_local = bad_s(kk,j,i)
     2063                      IF ( bad_local == 0.0_wp  .AND.  bad_s(kk,j-1,i) > 0.0_wp                    &
     2064                       .AND.  .NOT. building_edge_s )  bad_local = bad_s(kk,j-1,i)
    19922065!
    19932066!--                   In order to avoid that LAD is mistakenly considered at right-
     
    19952068!--                   u-component at index j,i is still on the building while the
    19962069!--                   topography top for the scalar isn't), LAD is taken from grid
    1997 !--                   point (j,i-1).
    1998                       IF ( lad_local > 0.0_wp .AND. lad_s(kk,j-1,i) == 0.0_wp  &
     2070!--                   point (j,i-1). The same is also valid for bad_s.
     2071                      IF ( lad_local > 0.0_wp  .AND.  lad_s(kk,j-1,i) == 0.0_wp                    &
    19992072                       .AND.  building_edge_n )  lad_local = lad_s(kk,j-1,i)
     2073
     2074                      IF ( bad_local > 0.0_wp  .AND.  bad_s(kk,j-1,i) == 0.0_wp                    &
     2075                       .AND.  building_edge_n )  bad_local = bad_s(kk,j-1,i)
    20002076
    20012077                      pre_tend = 0.0_wp
     
    20042080!--                   Calculate preliminary value (pre_tend) of the tendency
    20052081                      pre_tend = - canopy_drag_coeff *                         &
    2006                                    lad_local *                                 &
     2082                                   ( lad_local + bad_local ) *                 &
    20072083                                   SQRT( ( 0.25_wp * ( u(k,j-1,i)   +          &
    20082084                                                       u(k,j-1,i+1) +          &
     
    20542130!--                   Calculate preliminary value (pre_tend) of the tendency
    20552131                      pre_tend = - canopy_drag_coeff *                         &
    2056                                    (0.5_wp *                                   &
    2057                                       ( lad_s(kk+1,j,i) + lad_s(kk,j,i) )) *   &
     2132                                   ( 0.5_wp *                                  &
     2133                                      ( lad_s(kk+1,j,i) + lad_s(kk,j,i) ) +    &
     2134                                     0.5_wp *                                  &
     2135                                      ( bad_s(kk+1,j,i) + bad_s(kk,j,i) ) ) *  &
    20582136                                   SQRT( ( 0.25_wp * ( u(k,j,i)   +            &
    20592137                                                       u(k,j,i+1) +            &
     
    21612239                      tend(k,j,i) = tend(k,j,i) -                              &
    21622240                                       2.0_wp * canopy_drag_coeff *            &
    2163                                        lad_s(kk,j,i) *                         &
     2241                                       ( lad_s(kk,j,i) + bad_s(kk,j,i) ) *     &
    21642242                                       SQRT( ( 0.5_wp * ( u(k,j,i) +           &
    21652243                                                          u(k,j,i+1) )         &
     
    22532331       LOGICAL ::  building_edge_w !< control flag indicating a westward-facing building edge
    22542332
     2333       REAL(wp) ::  bad_local !< local lad value
    22552334       REAL(wp) ::  ddt_3d    !< inverse of the LES timestep (dt_3d)
    22562335       REAL(wp) ::  lad_local !< local lad value
     
    22962375!--             point (k,j,i), in order to avoid that LAD is mistakenly mapped
    22972376!--             on top of a roof where (usually) is no LAD is defined.
     2377!--             The same is also valid for bad_s.
    22982378                lad_local = lad_s(kk,j,i)
    2299                 IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp  .AND. &
     2379                IF ( lad_local == 0.0_wp  .AND.  lad_s(kk,j,i-1) > 0.0_wp  .AND.                  &
    23002380                     .NOT. building_edge_w )  lad_local = lad_s(kk,j,i-1)
     2381
     2382                bad_local = bad_s(kk,j,i)
     2383                IF ( bad_local == 0.0_wp  .AND.  bad_s(kk,j,i-1) > 0.0_wp  .AND.                   &
     2384                     .NOT. building_edge_w )  bad_local = bad_s(kk,j,i-1)
    23012385!
    23022386!--             In order to avoid that LAD is mistakenly considered at right-
     
    23042388!--             u-component at index j,i is still on the building while the
    23052389!--             topography top for the scalar isn't), LAD is taken from grid
    2306 !--             point (j,i-1).
    2307                 IF ( lad_local > 0.0_wp .AND. lad_s(kk,j,i-1) == 0.0_wp  .AND. &
     2390!--             point (j,i-1). The same is also valid for bad_s.
     2391                IF ( lad_local > 0.0_wp  .AND.  lad_s(kk,j,i-1) == 0.0_wp  .AND.                  &
    23082392                     building_edge_e )  lad_local = lad_s(kk,j,i-1)
     2393                IF ( bad_local > 0.0_wp  .AND.  bad_s(kk,j,i-1) == 0.0_wp  .AND.                   &
     2394                     building_edge_e )  bad_local = bad_s(kk,j,i-1)
    23092395
    23102396                pre_tend = 0.0_wp
     
    23132399!--             Calculate preliminary value (pre_tend) of the tendency
    23142400                pre_tend = - canopy_drag_coeff *                               &
    2315                              lad_local *                                       &   
     2401                             ( lad_local + bad_local ) *                       &
    23162402                             SQRT( u(k,j,i)**2 +                               &
    23172403                                   ( 0.25_wp * ( v(k,j,i-1)  +                 &
     
    23742460!--             point (k,j,i), in order to avoid that LAD is mistakenly mapped
    23752461!--             on top of a roof where (usually) is no LAD is defined.
     2462!--             The same is also valid for bad_s.
    23762463                lad_local = lad_s(kk,j,i)
    2377                 IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp  .AND. &
     2464                IF ( lad_local == 0.0_wp  .AND.  lad_s(kk,j-1,i) > 0.0_wp  .AND.                  &
    23782465                     .NOT. building_edge_s )  lad_local = lad_s(kk,j-1,i)
     2466
     2467                bad_local = bad_s(kk,j,i)
     2468                IF ( bad_local == 0.0_wp  .AND.  bad_s(kk,j-1,i) > 0.0_wp  .AND.                   &
     2469                     .NOT. building_edge_s )  bad_local = bad_s(kk,j-1,i)
    23792470!
    23802471!--             In order to avoid that LAD is mistakenly considered at right-
     
    23822473!--             u-component at index j,i is still on the building while the
    23832474!--             topography top for the scalar isn't), LAD is taken from grid
    2384 !--             point (j,i-1).
    2385                 IF ( lad_local > 0.0_wp .AND. lad_s(kk,j-1,i) == 0.0_wp  .AND. &
     2475!--             point (j,i-1). The same is also valid for bad_s.
     2476                IF ( lad_local > 0.0_wp  .AND.  lad_s(kk,j-1,i) == 0.0_wp  .AND.                  &
    23862477                     building_edge_n )  lad_local = lad_s(kk,j-1,i)
     2478                IF ( bad_local > 0.0_wp  .AND.  bad_s(kk,j-1,i) == 0.0_wp  .AND.                   &
     2479                     building_edge_n )  bad_local = bad_s(kk,j-1,i)
    23872480
    23882481                pre_tend = 0.0_wp
     
    23912484!--             Calculate preliminary value (pre_tend) of the tendency
    23922485                pre_tend = - canopy_drag_coeff *                               &
    2393                              lad_local *                                       &
     2486                             ( lad_local + bad_local ) *                       &
    23942487                             SQRT( ( 0.25_wp * ( u(k,j-1,i)   +                &
    23952488                                                 u(k,j-1,i+1) +                &
     
    24372530!--             Calculate preliminary value (pre_tend) of the tendency
    24382531                pre_tend = - canopy_drag_coeff *                               &
    2439                              (0.5_wp *                                         &
    2440                                 ( lad_s(kk+1,j,i) + lad_s(kk,j,i) )) *         &
    2441                              SQRT( ( 0.25_wp * ( u(k,j,i)    +                 & 
     2532                             ( 0.5_wp *                                        &
     2533                                ( lad_s(kk+1,j,i) + lad_s(kk,j,i) ) +          &
     2534                               0.5_wp *                                        &
     2535                                ( bad_s(kk+1,j,i) + bad_s(kk,j,i) ) ) *        &
     2536                             SQRT( ( 0.25_wp * ( u(k,j,i)    +                 &
    24422537                                                 u(k,j,i+1)  +                 &
    24432538                                                 u(k+1,j,i)  +                 &
     
    25272622                tend(k,j,i) = tend(k,j,i) -                                    &
    25282623                                 2.0_wp * canopy_drag_coeff *                  &
    2529                                  lad_s(kk,j,i) *                               &
     2624                                 ( lad_s(kk,j,i) + bad_s(kk,j,i) ) *           &
    25302625                                 SQRT( ( 0.5_wp * ( u(k,j,i) +                 &
    25312626                                                    u(k,j,i+1) )               &
     
    25402635                                 e(k,j,i)
    25412636             ENDDO
    2542              
    25432637!
    25442638!--       scalar concentration
Note: See TracChangeset for help on using the changeset viewer.