Changeset 1484 for palm


Ignore:
Timestamp:
Oct 21, 2014 10:53:05 AM (10 years ago)
Author:
kanani
Message:

New:
---
Subroutine init_plant_canopy added to module plant_canopy_model_mod. (plant_canopy_model)
Alternative method for lad-profile construction added, also, new parameters added.
(header, package_parin, plant_canopy_model, read_var_list, write_var_list)
plant_canopy_model-dependency added to several subroutines. (Makefile)
New package/namelist canopy_par for canopy-related parameters added. (package_parin)

Changed:
---
Code structure of the plant canopy model changed, all canopy-model related code
combined to module plant_canopy_model_mod. (check_parameters, init_3d_model,
modules, timestep)
Module plant_canopy_model_mod added in USE-lists of some subroutines. (check_parameters,
header, init_3d_model, package_parin, read_var_list, user_init_plant_canopy, write_var_list)
Canopy initialization moved to new subroutine init_plant_canopy. (check_parameters,
init_3d_model, plant_canopy_model)
Calculation of canopy timestep-criterion removed, instead, the canopy
drag is now directly limited in the calculation of the canopy tendency terms.
(plant_canopy_model, timestep)
Some parameters renamed. (check_parameters, header, init_plant_canopy,
plant_canopy_model, read_var_list, write_var_list)
Unnecessary 3d-arrays removed. (init_plant_canopy, plant_canopy_model, user_init_plant_canopy)
Parameter checks regarding canopy initialization added. (check_parameters)
All canopy steering parameters moved from namelist inipar to canopy_par. (package_parin, parin)
Some redundant MPI communication removed. (init_plant_canopy)

Bugfix:
---
Missing KIND-attribute for REAL constant added. (check_parameters)
DO-WHILE-loop for lad-profile output restricted. (header)
Removed double-listing of use_upstream_for_tke in ONLY-list of module
control_parameters. (prognostic_equations)

Location:
palm/trunk/SOURCE
Files:
13 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/Makefile

    r1445 r1484  
    2020# Current revisions:
    2121# ------------------
    22 #
     22# plant_canopy_model-dependency added for check_parameters, header, init_3d_model,
     23# package_parin, read_var_list, user_init_plant_canopy, write_var_list
    2324#
    2425# Former revisions:
     
    288289check_for_restart.o: modules.o mod_kinds.o
    289290check_open.o: modules.o mod_kinds.o mod_particle_attributes.o
    290 check_parameters.o: modules.o mod_kinds.o subsidence.o
     291check_parameters.o: modules.o mod_kinds.o subsidence.o plant_canopy_model.o
    291292close_file.o: modules.o mod_kinds.o
    292293compute_vpt.o: modules.o mod_kinds.o
     
    318319flow_statistics.o: modules.o cpulog.o mod_kinds.o
    319320global_min_max.o: modules.o mod_kinds.o
    320 header.o: modules.o cpulog.o mod_kinds.o subsidence.o
     321header.o: modules.o cpulog.o mod_kinds.o plant_canopy_model.o subsidence.o
    321322impact_of_latent_heat.o: modules.o mod_kinds.o
    322323inflow_turbulence.o: modules.o cpulog.o mod_kinds.o
    323324init_1d_model.o: modules.o mod_kinds.o
    324325init_3d_model.o: modules.o cpulog.o mod_kinds.o random_function.o advec_ws.o \
    325         ls_forcing.o lpm_init.o random_generator_parallel.o
     326        ls_forcing.o lpm_init.o plant_canopy_model.o random_generator_parallel.o
    326327init_advec.o: modules.o mod_kinds.o
    327328init_cloud_physics.o: modules.o mod_kinds.o
     
    379380netcdf.o: modules.o mod_kinds.o
    380381nudging.o: modules.o cpulog.o mod_kinds.o
    381 package_parin.o: modules.o mod_kinds.o
     382package_parin.o: modules.o mod_kinds.o plant_canopy_model.o
    382383palm.o: modules.o cpulog.o ls_forcing.o mod_kinds.o nudging.o
    383384parin.o: modules.o cpulog.o mod_kinds.o progress_bar.o
     
    400401random_generator_parallel.o: mod_kinds.o
    401402read_3d_binary.o: modules.o cpulog.o mod_kinds.o random_function.o random_generator_parallel.o
    402 read_var_list.o: modules.o mod_kinds.o
     403read_var_list.o: modules.o mod_kinds.o plant_canopy_model.o
    403404run_control.o: modules.o cpulog.o mod_kinds.o
    404405set_slicer_attributes_dvrp.o: modules.o mod_kinds.o
     
    434435user_init_3d_model.o: modules.o mod_kinds.o user_module.o
    435436user_init_grid.o: modules.o mod_kinds.o user_module.o
    436 user_init_plant_canopy.o: modules.o mod_kinds.o user_module.o
     437user_init_plant_canopy.o: modules.o mod_kinds.o user_module.o plant_canopy_model.o
    437438user_last_actions.o: modules.o mod_kinds.o user_module.o
    438439user_lpm_advec.o: modules.o mod_kinds.o user_module.o
     
    446447wall_fluxes.o: modules.o mod_kinds.o
    447448write_3d_binary.o: modules.o cpulog.o mod_kinds.o random_function.o random_generator_parallel.o
    448 write_var_list.o: modules.o mod_kinds.o
     449write_var_list.o: modules.o mod_kinds.o plant_canopy_model.o
  • palm/trunk/SOURCE/check_parameters.f90

    r1456 r1484  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Changes due to new module structure of the plant canopy model:
     23!   module plant_canopy_model_mod added,
     24!   checks regarding usage of new method for leaf area density profile
     25!   construction added,
     26!   lad-profile construction moved to new subroutine init_plant_canopy within
     27!   the module plant_canopy_model_mod,
     28!   drag_coefficient renamed to canopy_drag_coeff.
     29! Missing KIND-attribute for REAL constant added
    2330!
    2431! Former revisions:
     
    246253    USE particle_attributes
    247254    USE pegrid
     255    USE plant_canopy_model_mod
    248256    USE profil_parameter
    249257    USE spectrum
     
    891899    ENDIF
    892900
    893     IF ( plant_canopy .AND. ( drag_coefficient == 0.0_wp ) ) THEN
    894        message_string = 'plant_canopy = .TRUE. requires a non-zero drag ' // &
    895                         'coefficient & given value is drag_coefficient = 0.0'
    896        CALL message( 'check_parameters', 'PA0041', 1, 2, 0, 6, 0 )
    897     ENDIF
    898 
    899     IF ( plant_canopy  .AND.  cloud_physics  .AND.  icloud_scheme == 0 ) THEN
    900        message_string = 'plant_canopy = .TRUE. requires cloud_scheme /=' //  &
    901                         ' seifert_beheng'
    902        CALL message( 'check_parameters', 'PA0360', 1, 2, 0, 6, 0 )
    903     ENDIF
     901    IF ( plant_canopy )  THEN
     902   
     903       IF ( canopy_drag_coeff == 0.0_wp )  THEN
     904          message_string = 'plant_canopy = .TRUE. requires a non-zero drag '// &
     905                           'coefficient & given value is canopy_drag_coeff = 0.0'
     906          CALL message( 'check_parameters', 'PA0041', 1, 2, 0, 6, 0 )
     907       ENDIF
     908   
     909       IF ( ( alpha_lad /= 9999999.9_wp  .AND.  beta_lad == 9999999.9_wp )  .OR.&
     910              beta_lad /= 9999999.9_wp   .AND.  alpha_lad == 9999999.9_wp )  THEN
     911          message_string = 'using the beta function for the construction ' //  &
     912                           'of the leaf area density profile requires '    //  &
     913                           'both alpha_lad and beta_lad to be /= 9999999.9'
     914          CALL message( 'check_parameters', 'PA0118', 1, 2, 0, 6, 0 )
     915       ENDIF
     916   
     917       IF ( calc_beta_lad_profile  .AND.  lai_beta == 0.0_wp )  THEN
     918          message_string = 'using the beta function for the construction ' //  &
     919                           'of the leaf area density profile requires '    //  &
     920                           'a non-zero lai_beta, but given value is '      //  &
     921                           'lai_beta = 0.0'
     922          CALL message( 'check_parameters', 'PA0119', 1, 2, 0, 6, 0 )
     923       ENDIF
     924
     925       IF ( calc_beta_lad_profile  .AND.  lad_surface /= 0.0_wp )  THEN
     926          message_string = 'simultaneous setting of alpha_lad /= 9999999.9' // &
     927                           'and lad_surface /= 0.0 is not possible, '       // &
     928                           'use either vertical gradients or the beta '     // &
     929                           'function for the construction of the leaf area '// &
     930                           'density profile'
     931          CALL message( 'check_parameters', 'PA0120', 1, 2, 0, 6, 0 )
     932       ENDIF
     933
     934       IF ( cloud_physics  .AND.  icloud_scheme == 0 )  THEN
     935          message_string = 'plant_canopy = .TRUE. requires cloud_scheme /=' // &
     936                           ' seifert_beheng'
     937          CALL message( 'check_parameters', 'PA0360', 1, 2, 0, 6, 0 )
     938       ENDIF
     939
     940    ENDIF
    904941
    905942    IF ( .NOT. ( loop_optimization == 'cache'  .OR.                            &
     
    924961       IF ( ocean )           sa_init = sa_surface
    925962       IF ( passive_scalar )  q_init  = s_surface
    926        IF ( plant_canopy )    lad = 0.0_wp
    927963
    928964!
     
    12641300       ENDIF
    12651301
    1266 !
    1267 !--    If required compute the profile of leaf area density used in the plant
    1268 !--    canopy model
    1269        IF ( plant_canopy ) THEN
    1270        
    1271           i = 1
    1272           gradient = 0.0_wp
    1273 
    1274           IF ( .NOT. ocean ) THEN
    1275 
    1276              lad(0) = lad_surface
    1277  
    1278              lad_vertical_gradient_level_ind(1) = 0
    1279              DO k = 1, pch_index
    1280                 IF ( i < 11 ) THEN
    1281                    IF ( lad_vertical_gradient_level(i) < zu(k) .AND.  &
    1282                         lad_vertical_gradient_level(i) >= 0.0_wp ) THEN
    1283                       gradient = lad_vertical_gradient(i)
    1284                       lad_vertical_gradient_level_ind(i) = k - 1
    1285                       i = i + 1
    1286                    ENDIF
    1287                 ENDIF
    1288                 IF ( gradient /= 0.0_wp ) THEN
    1289                    IF ( k /= 1 ) THEN
    1290                       lad(k) = lad(k-1) + dzu(k) * gradient
    1291                    ELSE
    1292                       lad(k) = lad_surface + dzu(k) *gradient
    1293                    ENDIF
    1294                 ELSE
    1295                    lad(k) = lad(k-1)
    1296                 ENDIF
    1297              ENDDO
    1298 
    1299           ENDIF
    1300 
    1301 !
    1302 !--       In case of no given leaf area density gradients, choose a vanishing
    1303 !--       gradient
    1304           IF ( lad_vertical_gradient_level(1) == -9999999.9_wp ) THEN
    1305              lad_vertical_gradient_level(1) = 0.0_wp
    1306           ENDIF
    1307 
    1308        ENDIF
    13091302         
    13101303    ENDIF
     
    16311624       ENDIF
    16321625
    1633        IF ( top_heatflux /= 0.0  .AND.  top_heatflux /= 9999999.9_wp ) &
     1626       IF ( top_heatflux /= 0.0_wp  .AND.  top_heatflux /= 9999999.9_wp ) &
    16341627       THEN
    16351628          message_string = 'heatflux must not be set for pure neutral flow'
     
    18281821       rayleigh_damping_factor = 0.0_wp
    18291822    ELSE
    1830        IF ( rayleigh_damping_factor < 0.0 .OR. rayleigh_damping_factor > 1.0_wp ) &
     1823       IF ( rayleigh_damping_factor < 0.0_wp .OR. rayleigh_damping_factor > 1.0_wp ) &
    18311824       THEN
    18321825          WRITE( message_string, * )  'rayleigh_damping_factor = ', &
  • palm/trunk/SOURCE/header.f90

    r1483 r1484  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Changes due to new module structure of the plant canopy model:
     23!   module plant_canopy_model_mod and output for new canopy model parameters
     24!   (alpha_lad, beta_lad, lai_beta,...) added,
     25!   drag_coefficient, leaf_surface_concentration and scalar_exchange_coefficient
     26!   renamed to canopy_drag_coeff, leaf_surface_conc and leaf_scalar_exch_coeff,
     27!   learde renamed leaf_area_density.
     28! Bugfix: DO-WHILE-loop for lad header information additionally restricted
     29! by maximum number of gradient levels (currently 10)
    2330!
    2431! Former revisions:
     
    165172
    166173    USE arrays_3d,                                                             &
    167         ONLY:  lad, pt_init, qsws, q_init, sa_init, shf, ug, vg, w_subs, zu
     174        ONLY:  pt_init, qsws, q_init, sa_init, shf, ug, vg, w_subs, zu
    168175       
    169176    USE control_parameters
     
    205212       
    206213    USE pegrid
     214
     215    USE plant_canopy_model_mod,                                                &
     216        ONLY:  alpha_lad, beta_lad, calc_beta_lad_profile, canopy_drag_coeff,  &
     217               canopy_mode, cthf, lad, lad_surface, lad_vertical_gradient,     &
     218               lad_vertical_gradient_level, lad_vertical_gradient_level_ind,   &
     219               lai_beta, leaf_scalar_exch_coeff, leaf_surface_conc, pch_index, &
     220               plant_canopy
    207221   
    208222    USE spectrum,                                                              &
     
    242256    CHARACTER (LEN=86) ::  coordinates         !:
    243257    CHARACTER (LEN=86) ::  gradients           !:
    244     CHARACTER (LEN=86) ::  learde              !:
     258    CHARACTER (LEN=86) ::  leaf_area_density   !:
    245259    CHARACTER (LEN=86) ::  slices              !:
    246260    CHARACTER (LEN=86) ::  temperatures        !:
     
    270284    INTEGER(iwp) ::  io        !:
    271285    INTEGER(iwp) ::  j         !:
     286    INTEGER(iwp) ::  k         !:
    272287    INTEGER(iwp) ::  l         !:
    273288    INTEGER(iwp) ::  ll        !:
    274289    INTEGER(iwp) ::  mpi_type  !:
    275290   
     291    REAL(wp) ::  canopy_height                    !: canopy height (in m)
    276292    REAL(wp) ::  cpuseconds_per_simulated_second  !:
    277293
     
    760776
    761777    IF ( plant_canopy )  THEN
    762 
    763        WRITE ( io, 280 ) canopy_mode, pch_index, drag_coefficient
     778   
     779       canopy_height = pch_index * dz
     780
     781       WRITE ( io, 280 )  canopy_mode, canopy_height, pch_index,               &
     782                          canopy_drag_coeff
    764783       IF ( passive_scalar )  THEN
    765           WRITE ( io, 281 ) scalar_exchange_coefficient,   &
    766                             leaf_surface_concentration
     784          WRITE ( io, 281 )  leaf_scalar_exch_coeff,                           &
     785                             leaf_surface_conc
    767786       ENDIF
    768787
    769788!
    770789!--    Heat flux at the top of vegetation
    771        WRITE ( io, 282 ) cthf
    772 
    773 !
    774 !--    Leaf area density profile
    775 !--    Building output strings, starting with surface value
    776        WRITE ( learde, '(F6.4)' )  lad_surface
    777        gradients = '------'
    778        slices = '     0'
    779        coordinates = '   0.0'
    780        i = 1
    781        DO  WHILE ( lad_vertical_gradient_level_ind(i) /= -9999 )
    782 
    783           WRITE (coor_chr,'(F7.2)')  lad(lad_vertical_gradient_level_ind(i))
    784           learde = TRIM( learde ) // ' ' // TRIM( coor_chr )
    785 
    786           WRITE (coor_chr,'(F7.2)')  lad_vertical_gradient(i)
    787           gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr )
    788 
    789           WRITE (coor_chr,'(I7)')  lad_vertical_gradient_level_ind(i)
    790           slices = TRIM( slices ) // ' ' // TRIM( coor_chr )
    791 
    792           WRITE (coor_chr,'(F7.1)')  lad_vertical_gradient_level(i)
    793           coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
    794 
    795           i = i + 1
    796        ENDDO
    797 
    798        WRITE ( io, 283 )  TRIM( coordinates ), TRIM( learde ), &
    799                           TRIM( gradients ), TRIM( slices )
    800 
    801     ENDIF
     790       WRITE ( io, 282 )  cthf
     791
     792!
     793!--    Leaf area density profile, calculated either from given vertical
     794!--    gradients or from beta probability density function.
     795       IF (  .NOT.  calc_beta_lad_profile )  THEN
     796
     797!--       Building output strings, starting with surface value
     798          WRITE ( leaf_area_density, '(F6.4)' )  lad_surface
     799          gradients = '------'
     800          slices = '     0'
     801          coordinates = '   0.0'
     802          i = 1
     803          DO  WHILE ( i < 11  .AND.  lad_vertical_gradient_level_ind(i) /= -9999 )
     804
     805             WRITE (coor_chr,'(F7.2)')  lad(lad_vertical_gradient_level_ind(i))
     806             leaf_area_density = TRIM( leaf_area_density ) // ' ' // TRIM( coor_chr )
     807 
     808             WRITE (coor_chr,'(F7.2)')  lad_vertical_gradient(i)
     809             gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr )
     810
     811             WRITE (coor_chr,'(I7)')  lad_vertical_gradient_level_ind(i)
     812             slices = TRIM( slices ) // ' ' // TRIM( coor_chr )
     813
     814             WRITE (coor_chr,'(F7.1)')  lad_vertical_gradient_level(i)
     815             coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
     816
     817             i = i + 1
     818          ENDDO
     819
     820          WRITE ( io, 283 )  TRIM( coordinates ), TRIM( leaf_area_density ),              &
     821                             TRIM( gradients ), TRIM( slices )
     822
     823       ELSE
     824       
     825          WRITE ( leaf_area_density, '(F6.4)' )  lad_surface
     826          coordinates = '   0.0'
     827         
     828          DO  k = 1, pch_index
     829
     830             WRITE (coor_chr,'(F7.2)')  lad(k)
     831             leaf_area_density = TRIM( leaf_area_density ) // ' ' // TRIM( coor_chr )
     832 
     833             WRITE (coor_chr,'(F7.1)')  zu(k)
     834             coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
     835
     836          ENDDO       
     837
     838          WRITE ( io, 284 ) TRIM( coordinates ), TRIM( leaf_area_density ), alpha_lad,    &
     839                            beta_lad, lai_beta
     840
     841       ENDIF 
     842
     843    ENDIF
     844
    802845
    803846!
     
    18361879              ' ------------------------------'// &
    18371880              ' Canopy mode: ', A / &
    1838               ' Canopy top: ',I4 / &
     1881              ' Canopy height: ',F6.2,'m (',I4,' grid points)' / &
    18391882              ' Leaf drag coefficient: ',F6.2 /)
    1840 281 FORMAT (/ ' Scalar_exchange_coefficient: ',F6.2 / &
     1883281 FORMAT (/ ' Scalar exchange coefficient: ',F6.2 / &
    18411884              ' Scalar concentration at leaf surfaces in kg/m**3: ',F6.2 /)
    18421885282 FORMAT (' Predefined constant heatflux at the top of the vegetation: ',F6.2,' K m/s')
     
    18461889              ' Gradient:            ',A,'  m**2/m**4'/ &
    18471890              ' Gridpoint:           ',A)
     1891284 FORMAT (//' Characteristic levels of the leaf area density and coefficients:'// &
     1892              ' Height:              ',A,'  m'/ &
     1893              ' Leaf area density:   ',A,'  m**2/m**3'/ &
     1894              ' Coefficient alpha: ',F6.2 / &
     1895              ' Coefficient beta: ',F6.2 / &
     1896              ' Leaf area index: ',F6.2,'  m**2/m**2' /)
    18481897               
    18491898300 FORMAT (//' Boundary conditions:'/ &
  • palm/trunk/SOURCE/init_3d_model.f90

    r1432 r1484  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Changes due to new module structure of the plant canopy model:
     23!   canopy-related initialization (e.g. lad and canopy_heat_flux) moved to new
     24!   subroutine init_plant_canopy within the module plant_canopy_model_mod,
     25!   call of subroutine init_plant_canopy added.
    2326!
    2427! Former revisions:
     
    217220    USE pegrid
    218221   
     222    USE plant_canopy_model_mod,                                                &
     223        ONLY:  init_plant_canopy, plant_canopy
     224   
    219225    USE random_function_mod
    220226   
     
    517523   
    518524!
    519 !-- 3D-arrays for the leaf area density and the canopy drag coefficient
    520     IF ( plant_canopy ) THEN
    521        ALLOCATE ( lad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
    522                   lad_u(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
    523                   lad_v(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
    524                   lad_w(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
    525                   cdc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    526 
    527        IF ( passive_scalar ) THEN
    528           ALLOCATE ( sls(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
    529                      sec(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    530        ENDIF
    531 
    532        IF ( cthf /= 0.0_wp ) THEN
    533           ALLOCATE ( lai(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
    534                      canopy_heat_flux(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    535        ENDIF
    536 
    537     ENDIF
    538 
    539 !
    540525!-- 4D-array for storing the Rif-values at vertical walls
    541526    IF ( topography /= 'flat' )  THEN
     
    15711556
    15721557!
    1573 !-- Initialization of the leaf area density
    1574     IF ( plant_canopy )  THEN
    1575  
    1576        SELECT CASE ( TRIM( canopy_mode ) )
    1577 
    1578           CASE( 'block' )
    1579 
    1580              DO  i = nxlg, nxrg
    1581                 DO  j = nysg, nyng
    1582                    lad_s(:,j,i) = lad(:)
    1583                    cdc(:,j,i)   = drag_coefficient
    1584                    IF ( passive_scalar )  THEN
    1585                       sls(:,j,i) = leaf_surface_concentration
    1586                       sec(:,j,i) = scalar_exchange_coefficient
    1587                    ENDIF
    1588                 ENDDO
    1589              ENDDO
    1590 
    1591           CASE DEFAULT
    1592 
    1593 !
    1594 !--          The DEFAULT case is reached either if the parameter
    1595 !--          canopy mode contains a wrong character string or if the
    1596 !--          user has coded a special case in the user interface.
    1597 !--          There, the subroutine user_init_plant_canopy checks
    1598 !--          which of these two conditions applies.
    1599              CALL user_init_plant_canopy
    1600  
    1601           END SELECT
    1602 
    1603        CALL exchange_horiz( lad_s, nbgp )
    1604        CALL exchange_horiz( cdc, nbgp )
    1605 
    1606        IF ( passive_scalar )  THEN
    1607           CALL exchange_horiz( sls, nbgp )
    1608           CALL exchange_horiz( sec, nbgp )
    1609        ENDIF
    1610 
    1611 !
    1612 !--    Sharp boundaries of the plant canopy in horizontal directions
    1613 !--    In vertical direction the interpolation is retained, as the leaf
    1614 !--    area density is initialised by prescribing a vertical profile
    1615 !--    consisting of piecewise linear segments. The upper boundary
    1616 !--    of the plant canopy is now defined by lad_w(pch_index,:,:) = 0.0.
    1617 
    1618        DO  i = nxl, nxr
    1619           DO  j = nys, nyn
    1620              DO  k = nzb, nzt+1
    1621                 IF ( lad_s(k,j,i) > 0.0_wp )  THEN
    1622                    lad_u(k,j,i)   = lad_s(k,j,i)
    1623                    lad_u(k,j,i+1) = lad_s(k,j,i)
    1624                    lad_v(k,j,i)   = lad_s(k,j,i)
    1625                    lad_v(k,j+1,i) = lad_s(k,j,i)
    1626                 ENDIF
    1627              ENDDO
    1628              DO  k = nzb, nzt
    1629                 lad_w(k,j,i) = 0.5_wp * ( lad_s(k+1,j,i) + lad_s(k,j,i) )
    1630              ENDDO
    1631           ENDDO
    1632        ENDDO
    1633 
    1634        lad_w(pch_index,:,:) = 0.0_wp
    1635        lad_w(nzt+1,:,:)     = lad_w(nzt,:,:)
    1636 
    1637        CALL exchange_horiz( lad_u, nbgp )
    1638        CALL exchange_horiz( lad_v, nbgp )
    1639        CALL exchange_horiz( lad_w, nbgp )
    1640 
    1641 !
    1642 !--    Initialisation of the canopy heat source distribution
    1643        IF ( cthf /= 0.0_wp )  THEN
    1644 !
    1645 !--       Piecewise evaluation of the leaf area index by
    1646 !--       integration of the leaf area density
    1647           lai(:,:,:) = 0.0_wp
    1648           DO  i = nxlg, nxrg
    1649              DO  j = nysg, nyng
    1650                 DO  k = pch_index-1, 0, -1
    1651                    lai(k,j,i) = lai(k+1,j,i) +                   &
    1652                                 ( 0.5_wp * ( lad_w(k+1,j,i) +    &
    1653                                           lad_s(k+1,j,i) ) *     &
    1654                                   ( zw(k+1) - zu(k+1) ) )  +     &
    1655                                 ( 0.5_wp * ( lad_w(k,j,i)   +    &
    1656                                           lad_s(k+1,j,i) ) *     &
    1657                                   ( zu(k+1) - zw(k) ) )
    1658                 ENDDO
    1659              ENDDO
    1660           ENDDO
    1661 
    1662 !
    1663 !--       Evaluation of the upward kinematic vertical heat flux within the
    1664 !--       canopy
    1665           DO  i = nxlg, nxrg
    1666              DO  j = nysg, nyng
    1667                 DO  k = 0, pch_index
    1668                    canopy_heat_flux(k,j,i) = cthf *                    &
    1669                                              exp( -0.6_wp * lai(k,j,i) )
    1670                 ENDDO
    1671              ENDDO
    1672           ENDDO
    1673 
    1674 !
    1675 !--       The near surface heat flux is derived from the heat flux
    1676 !--       distribution within the canopy
    1677           shf(:,:) = canopy_heat_flux(0,:,:)
    1678 
    1679        ENDIF
    1680 
    1681     ENDIF
     1558!-- If required, initialize quantities needed for the plant canopy model
     1559    IF ( plant_canopy )  CALL init_plant_canopy
    16821560
    16831561!
  • palm/trunk/SOURCE/modules.f90

    r1469 r1484  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Changes due to new module structure of the plant canopy model:
     23!   canopy-model related parameters/variables moved to module
     24!   plant_canopy_model_mod
    2325!
    2426! Former revisions:
     
    3840! 1429 2014-07-15 12:53:45Z knoop
    3941! +ensemble_member_nr to prepare the random_generator for ensemble runs
    40 ! 
     42!
    4143! 1382 2014-04-30 12:15:41Z boeske
    4244! Renamed variables which store large scale forcing tendencies
     
    294296    REAL(wp), DIMENSION(:), ALLOCATABLE ::                                     &
    295297          c_u_m, c_u_m_l, c_v_m, c_v_m_l, c_w_m, c_w_m_l, ddzu, ddzu_pres,     &
    296           dd2zu, dzu, ddzw, dzw, hyp, inflow_damping_factor, lad, l_grid,      &
     298          dd2zu, dzu, ddzw, dzw, hyp, inflow_damping_factor, l_grid,           &
    297299          nc_1d, nr_1d, ptdf_x, ptdf_y, p_surf, pt_surf, pt_1d, pt_init,       &
    298300          qsws_surf, q_1d, q_init, q_surf, qc_1d, qr_1d, rdf, rdf_sc,          &
     
    313315
    314316    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
    315           canopy_heat_flux, cdc, d, de_dx, de_dy, de_dz, diss, diss_l_e,       &
     317          d, de_dx, de_dy, de_dz, diss, diss_l_e,                              &
    316318          diss_l_nr, diss_l_pt, diss_l_q, diss_l_qr, diss_l_sa, diss_l_u,      &
    317319          diss_l_v, diss_l_w, flux_l_e, flux_l_nr, flux_l_pt, flux_l_q,        &
    318           flux_l_qr, flux_l_sa, flux_l_u, flux_l_v, flux_l_w, kh, km, lad_s,   &
    319           lad_u, lad_v, lad_w, lai, l_wall, p_loc, sec, sls, tend, tric,       &
     320          flux_l_qr, flux_l_sa, flux_l_u, flux_l_v, flux_l_w, kh, km,          &
     321          l_wall, p_loc, tend, tric,                                           &
    320322          u_m_l, u_m_n, u_m_r, u_m_s, v_m_l, v_m_n, v_m_r, v_m_s, w_m_l,       &
    321323          w_m_n, w_m_r, w_m_s
     
    533535                             bc_sa_t = 'neumann', &
    534536                             bc_uv_b = 'dirichlet', bc_uv_t = 'dirichlet', &
    535                              canopy_mode = 'block', &
    536537                             cloud_scheme = 'seifert_beheng', &
    537538                             coupling_mode = 'uncoupled', &
     
    586587                     nr_timesteps_this_run = 0, &
    587588                     nsor = 20, nsor_ini = 100, n_sor, normalizing_region = 0, &
    588                      nz_do3d = -9999, pch_index = 0, prt_time_count = 0, &
     589                     nz_do3d = -9999, prt_time_count = 0, &
    589590                     recycling_plane, runnr = 0, &
    590591                     skip_do_avs = 0, subdomain_size, terminate_coupled = 0, &
     
    596597                     do3d_no(0:1) = 0, do3d_time_count(0:1), &
    597598                     domask_no(max_masks,0:1) = 0, domask_time_count(max_masks,0:1),&
    598                      lad_vertical_gradient_level_ind(10) = -9999, &
    599599                     mask_size(max_masks,3) = -1, mask_size_l(max_masks,3) = -1, &
    600600                     mask_start_l(max_masks,3) = -1, &
     
    652652                outflow_l = .FALSE., outflow_n = .FALSE., outflow_r = .FALSE., &
    653653                outflow_s = .FALSE., passive_scalar = .FALSE., &
    654                 plant_canopy = .FALSE., prandtl_layer = .TRUE., &
     654                prandtl_layer = .TRUE., &
    655655                precipitation = .FALSE., &
    656656                radiation = .FALSE., random_heatflux = .FALSE., &
     
    682682                 canyon_width_x = 9999999.9_wp, canyon_width_y = 9999999.9_wp, &
    683683                 canyon_wall_left = 9999999.9_wp, canyon_wall_south = 9999999.9_wp, &
    684                  cthf = 0.0_wp, cfl_factor = -1.0_wp, cos_alpha_surface, &
     684                 cfl_factor = -1.0_wp, cos_alpha_surface, &
    685685                 coupling_start_time = 0.0_wp, disturbance_amplitude = 0.25_wp, &
    686686                 disturbance_energy_limit = 0.01_wp, &
    687687                 disturbance_level_b = -9999999.9_wp, &
    688688                 disturbance_level_t = -9999999.9_wp, &
    689                  dp_level_b = 0.0_wp, drag_coefficient = 0.0_wp, &
     689                 dp_level_b = 0.0_wp, &
    690690                 dt = -1.0_wp, dt_averaging_input = 0.0_wp, &
    691691                 dt_averaging_input_pr = 9999999.9_wp, dt_coupling = 9999999.9_wp, &
     
    703703                 f = 0.0_wp, fs = 0.0_wp, g = 9.81_wp, inflow_damping_height = 9999999.9_wp, &
    704704                 inflow_damping_width = 9999999.9_wp, kappa = 0.4_wp, km_constant = -1.0_wp,&
    705                  lad_surface = 0.0_wp, leaf_surface_concentration = 0.0_wp, &
    706705                 mask_scale_x = 1.0_wp, mask_scale_y = 1.0_wp, mask_scale_z = 1.0_wp, &
    707706                 maximum_cpu_time_allowed = 0.0_wp,  &
     
    719718                 restart_time = 9999999.9_wp, rho_reference, rho_surface, &
    720719                 rif_max = 1.0_wp, rif_min = -5.0_wp, roughness_length = 0.1_wp, &
    721                  sa_surface = 35.0_wp, scalar_exchange_coefficient = 0.0_wp, &
     720                 sa_surface = 35.0_wp, &
    722721                 simulated_time = 0.0_wp, simulated_time_at_begin, sin_alpha_surface, &
    723722                 skip_time_data_output = 0.0_wp, skip_time_data_output_av = 9999999.9_wp,&
     
    745744    REAL(wp) ::  do2d_xy_last_time(0:1) = -1.0_wp, do2d_xz_last_time(0:1) = -1.0_wp, &
    746745                 do2d_yz_last_time(0:1) = -1.0_wp, dpdxy(1:2) = 0.0_wp, &
    747                  dt_domask(max_masks) = 9999999.9_wp, lad_vertical_gradient(10) = 0.0_wp,&
    748                  lad_vertical_gradient_level(10) = -9999999.9_wp, &
     746                 dt_domask(max_masks) = 9999999.9_wp, &
    749747                 mask_scale(3), &
    750748                 pt_vertical_gradient(10) = 0.0_wp, &
  • palm/trunk/SOURCE/package_parin.f90

    r1368 r1484  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Changes due to new module structure of the plant canopy model:
     23!   module plant_canopy_model_mod added,
     24!   new package/namelist canopy_par added, i.e. the canopy model is no longer
     25!   steered over the inipar-namelist,
     26!   drag_coefficient, leaf_surface_concentration and scalar_exchange_coefficient
     27!   renamed to canopy_drag_coeff, leaf_surface_conc and leaf_scalar_exch_coeff.
     28! Changed statement tags in CONTINUE-statement
    2329!
    2430! Former revisions:
     
    103109               write_particle_statistics
    104110
     111    USE plant_canopy_model_mod,                                                &
     112        ONLY:  alpha_lad, beta_lad, calc_beta_lad_profile, canopy_drag_coeff,  &
     113               canopy_mode, cthf, lad_surface,                                 &
     114               lad_vertical_gradient, lad_vertical_gradient_level, lai_beta,   &
     115               leaf_scalar_exch_coeff, leaf_surface_conc, pch_index,           &
     116               plant_canopy
     117
    105118    USE spectrum,                                                              &
    106119        ONLY:  comp_spectra_level, data_output_sp, plot_spectra_level,         &
     
    110123
    111124    CHARACTER (LEN=80) ::  line  !:
     125
     126    NAMELIST /canopy_par/         alpha_lad, beta_lad, canopy_drag_coeff,      &
     127                                  canopy_mode, cthf,                           &
     128                                  lad_surface,                                 &
     129                                  lad_vertical_gradient,                       &
     130                                  lad_vertical_gradient_level,                 &
     131                                  lai_beta,                                    &
     132                                  leaf_scalar_exch_coeff,                      &
     133                                  leaf_surface_conc, pch_index
    112134
    113135    NAMELIST /dvrp_graphics_par/  clip_dvrp_l, clip_dvrp_n, clip_dvrp_r,       &
     
    159181    line = ' '
    160182
     183!
     184!-- Try to find canopy package
     185    REWIND ( 11 )
     186    line = ' '
     187    DO   WHILE ( INDEX( line, '&canopy_par' ) == 0 )
     188       READ ( 11, '(A)', END=10 )  line
     189    ENDDO
     190    BACKSPACE ( 11 )
     191
     192!
     193!-- Read user-defined namelist
     194    READ ( 11, canopy_par )
     195
     196!
     197!-- Set flag that indicates that canopy model is switched on
     198    plant_canopy = .TRUE.
     199
     200!
     201!-- Set flag that indicates that the lad-profile shall be calculated by using
     202!-- a beta probability density function
     203    IF ( alpha_lad /= 9999999.9_wp  .AND.  beta_lad /= 9999999.9_wp )          &
     204       calc_beta_lad_profile = .TRUE.
     205
     206 10 CONTINUE
     207
     208
    161209#if defined( __dvrp_graphics )
    162210    REWIND ( 11 )
    163211    line = ' '
    164212    DO   WHILE ( INDEX( line, '&dvrp_graphics_par' ) == 0 )
    165        READ ( 11, '(A)', END=10 )  line
     213       READ ( 11, '(A)', END=20 )  line
    166214    ENDDO
    167215    BACKSPACE ( 11 )
     
    171219    READ ( 11, dvrp_graphics_par )
    172220
    173  10 CONTINUE
     221 20 CONTINUE
    174222#endif
    175223
     
    179227    line = ' '
    180228    DO   WHILE ( INDEX( line, '&particles_par' ) == 0 )
    181        READ ( 11, '(A)', END=20 )  line
     229       READ ( 11, '(A)', END=30 )  line
    182230    ENDDO
    183231    BACKSPACE ( 11 )
     
    191239    particle_advection = .TRUE.
    192240
    193  20 CONTINUE
     241 30 CONTINUE
    194242
    195243
     
    198246    line = ' '
    199247    DO   WHILE ( INDEX( line, '&spectra_par' ) == 0 )
    200        READ ( 11, '(A)', END=30 )  line
     248       READ ( 11, '(A)', END=40 )  line
    201249    ENDDO
    202250    BACKSPACE ( 11 )
     
    211259    IF ( dt_dosp == 9999999.9_wp )  dt_dosp = dt_data_output
    212260
    213  30 CONTINUE
     261 40 CONTINUE
    214262#endif
    215263
  • palm/trunk/SOURCE/parin.f90

    r1430 r1484  
    1919!
    2020! Current revisions:
    21 ! ------------------
    22 !
     21! -----------------
     22! Changes due to new module structure of the plant canopy model:
     23!   canopy-model related parameters moved to new package canopy_par in
     24!   subroutine package_parin
    2325!
    2426! Former revisions:
     
    151153
    152154    USE arrays_3d,                                                             &
    153         ONLY:  lad, pt_init, q_init, ref_state, sa_init, ug, u_init, v_init,   &
     155        ONLY:  pt_init, q_init, ref_state, sa_init, ug, u_init, v_init,        &
    154156               vg
    155157
     
    165167               building_length_y, building_wall_left, building_wall_south,     &
    166168               call_microphysics_at_all_substeps, call_psolver_at_all_substeps,&
    167                canopy_mode, canyon_height,                                     &
     169               canyon_height,                                                  &
    168170               canyon_width_x, canyon_width_y, canyon_wall_left,               &
    169171               canyon_wall_south, cfl_factor,                                  &
    170172               cloud_droplets, cloud_physics, cloud_scheme,                    &
    171173               conserve_volume_flow, conserve_volume_flow_mode,                &
    172                coupling_start_time, create_disturbances, cthf, cycle_mg,       &
     174               coupling_start_time, create_disturbances, cycle_mg,             &
    173175               data_output, data_output_masks,                                 &
    174176               data_output_pr, data_output_2d_on_each_pe,                      &
     
    176178               disturbance_level_b, disturbance_level_t, dissipation_1d,       &
    177179               do2d_at_begin, do3d_at_begin,                                   &
    178                dp_external, dp_level_b, dp_smooth, dpdxy, drag_coefficient,    &
     180               dp_external, dp_level_b, dp_smooth, dpdxy,                      &
    179181               drizzle, dt, dt_averaging_input, dt_averaging_input_pr,         &
    180182               dt_coupling, dt_data_output, dt_data_output_av, dt_disturb,     &
     
    187189               inflow_damping_width, inflow_disturbance_begin,                 &
    188190               inflow_disturbance_end,  initializing_actions, io_blocks,       &
    189                io_group, km_constant, lad_surface, lad_vertical_gradient,      &
    190                lad_vertical_gradient_level, large_scale_forcing,               &
    191                large_scale_subsidence, leaf_surface_concentration,             &
     191               io_group, km_constant,                                          &
     192               large_scale_forcing,                                            &
     193               large_scale_subsidence,                                         &
    192194               loop_optimization, masking_method, mask_scale_x, mask_scale_y,  &
    193195               mask_scale_z, mask_x, mask_y, mask_z, mask_x_loop,              &
     
    197199               momentum_advec, netcdf_data_format, netcdf_precision, neutral,  &
    198200               ngsrb, normalizing_region, nsor, nsor_ini, nudging, ocean,      &
    199                omega, omega_sor, passive_scalar, pch_index, phi, nz_do3d,      &
    200                plant_canopy, prandtl_layer, prandtl_number, precipitation,     &
     201               omega, omega_sor, passive_scalar, phi, nz_do3d,                 &
     202               prandtl_layer, prandtl_number, precipitation,                   &
    201203               precipitation_amount_interval, psolver, pt_damping_factor,      &
    202204               pt_damping_width, pt_reference, pt_surface,                     &
     
    211213               run_identifier, sa_surface, sa_vertical_gradient,               &
    212214               sa_vertical_gradient_level, scalar_advec,                       &
    213                scalar_exchange_coefficient, scalar_rayleigh_damping,           &
     215               scalar_rayleigh_damping,                                        &
    214216               section_xy, section_xz, section_yz, skip_time_data_output,      &
    215217               skip_time_data_output_av, skip_time_dopr, skip_time_do2d_xy,    &
     
    270272             building_length_y, building_wall_left, building_wall_south,       &
    271273             call_psolver_at_all_substeps, call_microphysics_at_all_substeps,  &
    272              canopy_mode, canyon_height,                                       &
     274             canyon_height,                                                    &
    273275             canyon_width_x, canyon_width_y, canyon_wall_left,                 &
    274276             canyon_wall_south, c_sedimentation, cfl_factor, cloud_droplets,   &
    275277             cloud_physics, cloud_scheme, collective_wait,                     &
    276278             conserve_volume_flow,                                             &
    277              conserve_volume_flow_mode, coupling_start_time, cthf,             &
     279             conserve_volume_flow_mode, coupling_start_time,                   &
    278280             curvature_solution_effects, cycle_mg, damp_level_1d,              &
    279281             dissipation_1d,                                                   &
    280              dp_external, dp_level_b, dp_smooth, dpdxy, drag_coefficient,      &
     282             dp_external, dp_level_b, dp_smooth, dpdxy,                        &
    281283             drizzle, dt, dt_pr_1d, dt_run_control_1d, dx, dy, dz, dz_max,     &
    282284             dz_stretch_factor, dz_stretch_level, end_time_1d,                 &
     
    285287             inflow_damping_height, inflow_damping_width,                      &
    286288             inflow_disturbance_begin, inflow_disturbance_end,                 &
    287              initializing_actions, km_constant, lad_surface,                   &
    288              lad_vertical_gradient, lad_vertical_gradient_level,               &
     289             initializing_actions, km_constant,                                &
    289290             large_scale_forcing, large_scale_subsidence,                      &
    290              leaf_surface_concentration, limiter_sedimentation,                &
     291             limiter_sedimentation, &
    291292             loop_optimization, masking_method, mg_cycles,                     &
    292293             mg_switch_to_pe0_level, mixing_length_1d, momentum_advec,         &
    293294             nc_const, netcdf_precision, neutral, ngsrb,                       &
    294295             nsor, nsor_ini, nudging, nx, ny, nz, ocean, omega, omega_sor,     &
    295              passive_scalar, pch_index, phi, plant_canopy, prandtl_layer,      &
     296             passive_scalar, phi, prandtl_layer, &
    296297             prandtl_number, precipitation, psolver, pt_damping_factor,        &
    297298             pt_damping_width, pt_reference, pt_surface,                       &
     
    304305             rif_max, rif_min, roughness_length, sa_surface,                   &
    305306             sa_vertical_gradient, sa_vertical_gradient_level, scalar_advec,   &
    306              scalar_exchange_coefficient, scalar_rayleigh_damping,             &
     307             scalar_rayleigh_damping,                                          &
    307308             statistic_regions, subs_vertical_gradient,                        &
    308309             subs_vertical_gradient_level, surface_heatflux, surface_pressure, &
     
    478479!--          ATTENTION: in case of changes to the following statement please
    479480!--                  also check the allocate statement in routine read_var_list
    480              ALLOCATE( lad(0:nz+1), pt_init(0:nz+1), q_init(0:nz+1),           &
     481             ALLOCATE( pt_init(0:nz+1), q_init(0:nz+1),                       &
    481482                       ref_state(0:nz+1), sa_init(0:nz+1), ug(0:nz+1),         &
    482483                       u_init(0:nz+1), v_init(0:nz+1), vg(0:nz+1),             &
  • palm/trunk/SOURCE/plant_canopy_model.f90

    r1341 r1484  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Changes due to new module structure of the plant canopy model:
     23!   module plant_canopy_model_mod now contains a subroutine for the
     24!   initialization of the canopy model (init_plant_canopy),
     25!   limitation of the canopy drag (previously accounted for by calculation of
     26!   a limiting canopy timestep for the determination of the maximum LES timestep
     27!   in subroutine timestep) is now realized by the calculation of pre-tendencies
     28!   and preliminary velocities in subroutine plant_canopy_model,
     29!   some redundant MPI communication removed in subroutine init_plant_canopy
     30!   (was previously in init_3d_model),
     31!   unnecessary 3d-arrays lad_u, lad_v, lad_w removed - lad information on the
     32!   respective grid is now provided only by lad_s (e.g. in the calculation of
     33!   the tendency terms or of cum_lai_hf),
     34!   drag_coefficient, lai, leaf_surface_concentration,
     35!   scalar_exchange_coefficient, sec and sls renamed to canopy_drag_coeff,
     36!   cum_lai_hf, leaf_surface_conc, leaf_scalar_exch_coeff, lsec and lsc,
     37!   respectively,
     38!   unnecessary 3d-arrays cdc, lsc and lsec now defined as single-value constants,
     39!   USE-statements and ONLY-lists modified accordingly
    2340!
    2441! Former revisions:
     
    4663! Description:
    4764! ------------
    48 ! Evaluation of sinks and sources of momentum, heat and scalar concentration
    49 ! due to canopy elements
     65! 1) Initialization of the canopy model, e.g. construction of leaf area density
     66! profile (subroutine init_plant_canopy).
     67! 2) Calculation of sinks and sources of momentum, heat and scalar concentration
     68! due to canopy elements (subroutine plant_canopy_model).
    5069!------------------------------------------------------------------------------!
     70    USE arrays_3d,                                                             &
     71        ONLY:  dzu, dzw, e, q, shf, tend, u, v, w, zu, zw
     72
     73    USE indices,                                                               &
     74        ONLY:  nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv,   &
     75               nz, nzb, nzb_s_inner, nzb_u_inner, nzb_v_inner, nzb_w_inner, nzt
     76
     77    USE kinds
     78
     79
     80    IMPLICIT NONE
     81
     82
     83    CHARACTER (LEN=20)   ::  canopy_mode = 'block' !: canopy coverage
     84
     85    INTEGER(iwp) ::  pch_index = 0                 !: plant canopy height/top index
     86    INTEGER(iwp) ::                                                            &
     87       lad_vertical_gradient_level_ind(10) = -9999 !: lad-profile levels (index)
     88
     89    LOGICAL ::  calc_beta_lad_profile = .FALSE. !: switch for calc. of lad from beta func.
     90    LOGICAL ::  plant_canopy = .FALSE.          !: switch for use of canopy model
     91
     92    REAL(wp) ::  alpha_lad = 9999999.9_wp   !: coefficient for lad calculation
     93    REAL(wp) ::  beta_lad = 9999999.9_wp    !: coefficient for lad calculation
     94    REAL(wp) ::  canopy_drag_coeff = 0.0_wp !: canopy drag coefficient (parameter)
     95    REAL(wp) ::  cdc = 0.0_wp               !: canopy drag coeff. (abbreviation used in equations)
     96    REAL(wp) ::  cthf = 0.0_wp              !: canopy top heat flux
     97    REAL(wp) ::  dt_plant_canopy = 0.0_wp   !: timestep account. for canopy drag
     98    REAL(wp) ::  lad_surface = 0.0_wp       !: lad surface value
     99    REAL(wp) ::  lai_beta = 0.0_wp          !: leaf area index (lai) for lad calc.
     100    REAL(wp) ::                                                                &
     101       leaf_scalar_exch_coeff = 0.0_wp      !: canopy scalar exchange coeff.
     102    REAL(wp) ::                                                                &
     103       leaf_surface_conc = 0.0_wp           !: leaf surface concentration
     104    REAL(wp) ::  lsec = 0.0_wp              !: leaf scalar exchange coeff.
     105    REAL(wp) ::  lsc = 0.0_wp               !: leaf surface concentration
     106
     107    REAL(wp) ::                                                                &
     108       lad_vertical_gradient(10) = 0.0_wp              !: lad gradient
     109    REAL(wp) ::                                                                &
     110       lad_vertical_gradient_level(10) = -9999999.9_wp !: lad-prof. levels (in m)
     111
     112    REAL(wp), DIMENSION(:), ALLOCATABLE ::  lad            !: leaf area density
     113    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pre_lad        !: preliminary lad
     114   
     115    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
     116       canopy_heat_flux                                    !: canopy heat flux
     117    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  cum_lai_hf !: cumulative lai for heatflux calc.
     118    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  lad_s      !: lad on scalar-grid
     119
     120
     121    SAVE
     122
    51123
    52124    PRIVATE
    53     PUBLIC plant_canopy_model
     125    PUBLIC alpha_lad, beta_lad, calc_beta_lad_profile, canopy_drag_coeff,      &
     126           canopy_mode, cdc, cthf, dt_plant_canopy, init_plant_canopy, lad,    &
     127           lad_s, lad_surface, lad_vertical_gradient,                          &
     128           lad_vertical_gradient_level, lad_vertical_gradient_level_ind,       &
     129           lai_beta, leaf_scalar_exch_coeff, leaf_surface_conc, lsc, lsec,     &
     130           pch_index, plant_canopy, plant_canopy_model
     131
     132
     133    INTERFACE init_plant_canopy
     134       MODULE PROCEDURE init_plant_canopy
     135    END INTERFACE init_plant_canopy
    54136
    55137    INTERFACE plant_canopy_model
     
    58140    END INTERFACE plant_canopy_model
    59141
     142
     143
     144
    60145 CONTAINS
    61146
    62147
    63148!------------------------------------------------------------------------------!
    64 ! Call for all grid points
     149! Description:
     150! ------------
     151!-- Initialization of the plant canopy model
     152!------------------------------------------------------------------------------!
     153    SUBROUTINE init_plant_canopy
     154   
     155
     156       USE control_parameters,                                                 &
     157           ONLY: dz, ocean, passive_scalar
     158
     159
     160       IMPLICIT NONE
     161
     162       INTEGER(iwp) ::  i !: running index
     163       INTEGER(iwp) ::  j !: running index
     164       INTEGER(iwp) ::  k !: running index
     165
     166       REAL(wp) ::  int_bpdf      !: vertical integral for lad-profile construction
     167       REAL(wp) ::  dzh           !: vertical grid spacing in units of canopy height
     168       REAL(wp) ::  gradient      !: gradient for lad-profile construction
     169       REAL(wp) ::  canopy_height !: canopy height for lad-profile construction
     170
     171!
     172!--    Allocate one-dimensional arrays for the computation of the
     173!--    leaf area density (lad) profile
     174       ALLOCATE( lad(0:nz+1), pre_lad(0:nz+1) )
     175       lad = 0.0_wp
     176       pre_lad = 0.0_wp
     177
     178!
     179!--    Compute the profile of leaf area density used in the plant
     180!--    canopy model. The profile can either be constructed from
     181!--    prescribed vertical gradients of the leaf area density or by
     182!--    using a beta probability density function (see e.g. Markkanen et al.,
     183!--    2003: Boundary-Layer Meteorology, 106, 437-459)
     184       IF (  .NOT.  calc_beta_lad_profile )  THEN   
     185
     186!
     187!--       Use vertical gradients for lad-profile construction   
     188          i = 1
     189          gradient = 0.0_wp
     190
     191          IF (  .NOT.  ocean )  THEN
     192
     193             lad(0) = lad_surface
     194             lad_vertical_gradient_level_ind(1) = 0
     195 
     196             DO k = 1, pch_index
     197                IF ( i < 11 )  THEN
     198                   IF ( lad_vertical_gradient_level(i) < zu(k)  .AND.          &
     199                        lad_vertical_gradient_level(i) >= 0.0_wp )  THEN
     200                      gradient = lad_vertical_gradient(i)
     201                      lad_vertical_gradient_level_ind(i) = k - 1
     202                      i = i + 1
     203                   ENDIF
     204                ENDIF
     205                IF ( gradient /= 0.0_wp )  THEN
     206                   IF ( k /= 1 )  THEN
     207                      lad(k) = lad(k-1) + dzu(k) * gradient
     208                   ELSE
     209                      lad(k) = lad_surface + dzu(k) * gradient
     210                   ENDIF
     211                ELSE
     212                   lad(k) = lad(k-1)
     213                ENDIF
     214             ENDDO
     215
     216          ENDIF
     217
     218!
     219!--       In case of no given leaf area density gradients, choose a vanishing
     220!--       gradient. This information is used for the HEADER and the RUN_CONTROL
     221!--       file.
     222          IF ( lad_vertical_gradient_level(1) == -9999999.9_wp )  THEN
     223             lad_vertical_gradient_level(1) = 0.0_wp
     224          ENDIF
     225
     226       ELSE
     227
     228!
     229!--       Use beta function for lad-profile construction
     230          int_bpdf = 0.0_wp
     231          canopy_height = pch_index * dz
     232
     233          DO k = nzb, pch_index
     234             int_bpdf = int_bpdf +                                             &
     235                           ( ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) ) * &
     236                           ( ( 1.0_wp - ( zw(k) / canopy_height ) )**( beta_lad-1.0_wp ) ) *   &
     237                           ( ( zw(k+1)-zw(k) ) / canopy_height ) )
     238          ENDDO
     239
     240!
     241!--       Preliminary lad profile (defined on w-grid)
     242          DO k = nzb, pch_index
     243             pre_lad(k) =  lai_beta *                                              &
     244                           (   ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) ) *          &
     245                               ( ( 1.0_wp - ( zw(k) / canopy_height ) )**( beta_lad-1.0_wp ) ) /   &
     246                               int_bpdf                                            &
     247                           ) / canopy_height
     248          ENDDO
     249
     250!
     251!--       Final lad profile (defined on scalar-grid level, since most prognostic
     252!--       quantities are defined there, hence, less interpolation is required
     253!--       when calculating the canopy tendencies)
     254          lad(0) = pre_lad(0)
     255          DO k = nzb+1, pch_index
     256             lad(k) = 0.5 * ( pre_lad(k-1) + pre_lad(k) )
     257          ENDDO         
     258
     259       ENDIF
     260
     261!
     262!--    Allocate 3D-array for the leaf area density (lad_s). In case of a
     263!--    prescribed canopy-top heat flux (cthf), allocate 3D-arrays for
     264!--    the cumulative leaf area index (cum_lai_hf) and the canopy heat flux.
     265       ALLOCATE( lad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     266
     267       IF ( cthf /= 0.0_wp )  THEN
     268          ALLOCATE( cum_lai_hf(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                 &
     269                    canopy_heat_flux(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     270       ENDIF
     271
     272!
     273!--    Initialize canopy parameters cdc (canopy drag coefficient),
     274!--    lsec (leaf scalar exchange coefficient), lsc (leaf surface concentration)
     275!--    with the prescribed values
     276       cdc = canopy_drag_coeff
     277       lsec = leaf_scalar_exch_coeff
     278       lsc = leaf_surface_conc
     279
     280!
     281!--    Initialization of the canopy coverage in the model domain:
     282!--    Setting the parameter canopy_mode = 'block' initializes a canopy, which
     283!--    fully covers the domain surface
     284       SELECT CASE ( TRIM( canopy_mode ) )
     285
     286          CASE( 'block' )
     287
     288             DO  i = nxlg, nxrg
     289                DO  j = nysg, nyng
     290                   lad_s(:,j,i) = lad(:)
     291                ENDDO
     292             ENDDO
     293
     294          CASE DEFAULT
     295
     296!
     297!--       The DEFAULT case is reached either if the parameter
     298!--       canopy mode contains a wrong character string or if the
     299!--       user has coded a special case in the user interface.
     300!--       There, the subroutine user_init_plant_canopy checks
     301!--       which of these two conditions applies.
     302          CALL user_init_plant_canopy
     303 
     304       END SELECT
     305
     306!
     307!--    Initialization of the canopy heat source distribution
     308       IF ( cthf /= 0.0_wp )  THEN
     309!
     310!--       Piecewise calculation of the leaf area index by vertical
     311!--       integration of the leaf area density
     312          cum_lai_hf(:,:,:) = 0.0_wp
     313          DO  i = nxlg, nxrg
     314             DO  j = nysg, nyng
     315                DO  k = pch_index-1, 0, -1
     316                   IF ( k == pch_index-1 )  THEN
     317                      cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) +                &
     318                         ( 0.5_wp * lad_s(k+1,j,i) *                           &
     319                           ( zw(k+1) - zu(k+1) ) )  +                          &
     320                         ( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) +              &
     321                                                 lad_s(k,j,i) ) +              &
     322                                      lad_s(k+1,j,i) ) *                       &
     323                           ( zu(k+1) - zw(k) ) )
     324                   ELSE
     325                      cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) +                &
     326                         ( 0.5_wp * ( 0.5_wp * ( lad_s(k+2,j,i) +              &
     327                                                 lad_s(k+1,j,i) ) +            &
     328                                      lad_s(k+1,j,i) ) *                       &
     329                           ( zw(k+1) - zu(k+1) ) )  +                          &
     330                         ( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) +              &
     331                                                 lad_s(k,j,i) ) +              &
     332                                      lad_s(k+1,j,i) ) *                       &
     333                           ( zu(k+1) - zw(k) ) )
     334                   ENDIF
     335                ENDDO
     336             ENDDO
     337          ENDDO
     338
     339!
     340!--       Calculation of the upward kinematic vertical heat flux within the
     341!--       canopy
     342          DO  i = nxlg, nxrg
     343             DO  j = nysg, nyng
     344                DO  k = 0, pch_index
     345                   canopy_heat_flux(k,j,i) = cthf *                            &
     346                                             exp( -0.6_wp * cum_lai_hf(k,j,i) )
     347                ENDDO
     348             ENDDO
     349          ENDDO
     350
     351!
     352!--       The surface heat flux is set to the surface value of the calculated
     353!--       in-canopy heat flux distribution 
     354          shf(:,:) = canopy_heat_flux(0,:,:)
     355
     356       ENDIF
     357
     358
     359
     360    END SUBROUTINE init_plant_canopy
     361
     362
     363
     364!------------------------------------------------------------------------------!
     365! Description:
     366! ------------
     367!-- Calculation of the tendency terms, accounting for the effect of the plant
     368!-- canopy on momentum and scalar quantities.
     369!--
     370!-- The canopy is located where the leaf area density lad_s(k,j,i) > 0.0
     371!-- (defined on scalar grid), as initialized in subroutine init_plant_canopy.
     372!-- The lad on the w-grid is vertically interpolated from the surrounding
     373!-- lad_s. The upper boundary of the canopy is defined on the w-grid at
     374!-- k = pch_index. Here, the lad is zero.
     375!--
     376!-- The canopy drag must be limited (previously accounted for by calculation of
     377!-- a limiting canopy timestep for the determination of the maximum LES timestep
     378!-- in subroutine timestep), since it is physically impossible that the canopy
     379!-- drag alone can locally change the sign of a velocity component. This
     380!-- limitation is realized by calculating preliminary tendencies and velocities.
     381!-- It is subsequently checked if the preliminary new velocity has a different
     382!-- sign than the current velocity. If so, the tendency is limited in a way that
     383!-- the velocity can at maximum be reduced to zero by the canopy drag.
     384!--
     385!--
     386!-- Call for all grid points
    65387!------------------------------------------------------------------------------!
    66388    SUBROUTINE plant_canopy_model( component )
    67389
    68        USE arrays_3d,                                                          &
    69            ONLY:  canopy_heat_flux, cdc, dzw, e, lad_s, lad_u, lad_v, lad_w,   &
    70                   q, sec, sls, tend, u, v, w
    71390
    72391       USE control_parameters,                                                 &
    73            ONLY:  pch_index, message_string
    74 
    75        USE indices,                                                            &
    76            ONLY:  nxl, nxlu, nxr, nys, nysv, nyn, nzb_s_inner, nzb_u_inner,    &
    77                   nzb_v_inner, nzb_w_inner
     392           ONLY:  dt_3d, message_string
    78393
    79394       USE kinds
     
    81396       IMPLICIT NONE
    82397
    83        INTEGER(iwp) ::  component  !:
    84        INTEGER(iwp) ::  i          !:
    85        INTEGER(iwp) ::  j          !:
    86        INTEGER(iwp) ::  k          !:
     398       INTEGER(iwp) ::  component !: prognostic variable (u,v,w,pt,q,e)
     399       INTEGER(iwp) ::  i         !: running index
     400       INTEGER(iwp) ::  j         !: running index
     401       INTEGER(iwp) ::  k         !: running index
     402
     403       REAL(wp) ::  ddt_3d    !: inverse of the LES timestep (dt_3d)
     404       REAL(wp) ::  lad_local !: local lad value
     405       REAL(wp) ::  pre_tend  !: preliminary tendency
     406       REAL(wp) ::  pre_u     !: preliminary u-value
     407       REAL(wp) ::  pre_v     !: preliminary v-value
     408       REAL(wp) ::  pre_w     !: preliminary w-value
     409
     410
     411       ddt_3d = 1.0_wp / dt_3d
    87412 
    88413!
    89 !--    Compute drag for the three velocity components and the SGS-TKE
     414!--    Compute drag for the three velocity components and the SGS-TKE:
    90415       SELECT CASE ( component )
    91416
     
    96421                DO  j = nys, nyn
    97422                   DO  k = nzb_u_inner(j,i)+1, pch_index
    98                       tend(k,j,i) = tend(k,j,i) -                &
    99                                     cdc(k,j,i) * lad_u(k,j,i) *  &
    100                                     SQRT(     u(k,j,i)**2     +  &
    101                                           ( ( v(k,j,i-1)      +  &
    102                                               v(k,j,i)        +  &
    103                                               v(k,j+1,i)      +  &
    104                                               v(k,j+1,i-1) )     &
    105                                             / 4.0_wp )**2     +  &
    106                                           ( ( w(k-1,j,i-1)    +  &
    107                                               w(k-1,j,i)      +  &
    108                                               w(k,j,i-1)      +  &
    109                                               w(k,j,i) )         &
    110                                             / 4.0_wp )**2 )   *  &
    111                                     u(k,j,i)
     423
     424!
     425!--                   In order to create sharp boundaries of the plant canopy,
     426!--                   the lad on the u-grid at index (k,j,i) is equal to
     427!--                   lad_s(k,j,i), rather than being interpolated from the
     428!--                   surrounding lad_s, because this would yield smaller lad
     429!--                   at the canopy boundaries than inside of the canopy.
     430!--                   For the same reason, the lad at the rightmost(i+1)canopy
     431!--                   boundary on the u-grid equals lad_s(k,j,i).
     432                      lad_local = lad_s(k,j,i)
     433                      IF ( lad_local == 0.0_wp  .AND.                          &
     434                           lad_s(k,j,i-1) > 0.0_wp )  THEN
     435                         lad_local = lad_s(k,j,i-1)
     436                      ENDIF
     437
     438                      pre_tend = 0.0_wp
     439                      pre_u = 0.0_wp
     440!
     441!--                   Calculate preliminary value (pre_tend) of the tendency
     442                      pre_tend = - cdc *                                       &
     443                                   lad_local *                                 &
     444                                   SQRT( u(k,j,i)**2 +                         &
     445                                         ( 0.25_wp * ( v(k,j,i-1) +            &
     446                                                       v(k,j,i)   +            &
     447                                                       v(k,j+1,i) +            &
     448                                                       v(k,j+1,i-1) )          &
     449                                         )**2 +                                &
     450                                         ( 0.25_wp * ( w(k-1,j,i-1) +          &
     451                                                       w(k-1,j,i)   +          &
     452                                                       w(k,j,i-1)   +          &
     453                                                       w(k,j,i) )              &
     454                                         )**2                                  &
     455                                       ) *                                     &
     456                                   u(k,j,i)
     457
     458!
     459!--                   Calculate preliminary new velocity, based on pre_tend
     460                      pre_u = u(k,j,i) + dt_3d * pre_tend
     461!
     462!--                   Compare sign of old velocity and new preliminary velocity,
     463!--                   and in case the signs are different, limit the tendency
     464                      IF ( SIGN(pre_u,u(k,j,i)) /= pre_u )  THEN
     465                         pre_tend = - u(k,j,i) * ddt_3d
     466                      ELSE
     467                         pre_tend = pre_tend
     468                      ENDIF
     469!
     470!--                   Calculate final tendency
     471                      tend(k,j,i) = tend(k,j,i) + pre_tend
     472
    112473                   ENDDO
    113474                ENDDO
     
    120481                DO  j = nysv, nyn
    121482                   DO  k = nzb_v_inner(j,i)+1, pch_index
    122                       tend(k,j,i) = tend(k,j,i) -                &
    123                                     cdc(k,j,i) * lad_v(k,j,i) *  &
    124                                     SQRT( ( ( u(k,j-1,i)      +  &
    125                                               u(k,j-1,i+1)    +  &
    126                                               u(k,j,i)        +  &
    127                                               u(k,j,i+1) )       &
    128                                             / 4.0_wp )**2     +  &
    129                                               v(k,j,i)**2     +  &
    130                                           ( ( w(k-1,j-1,i)    +  &
    131                                               w(k-1,j,i)      +  &
    132                                               w(k,j-1,i)      +  &
    133                                               w(k,j,i) )         &
    134                                             / 4.0_wp )**2 )   *  &
    135                                     v(k,j,i) 
     483
     484!
     485!--                   In order to create sharp boundaries of the plant canopy,
     486!--                   the lad on the v-grid at index (k,j,i) is equal to
     487!--                   lad_s(k,j,i), rather than being interpolated from the
     488!--                   surrounding lad_s, because this would yield smaller lad
     489!--                   at the canopy boundaries than inside of the canopy.
     490!--                   For the same reason, the lad at the northmost(j+1) canopy
     491!--                   boundary on the v-grid equals lad_s(k,j,i).
     492                      lad_local = lad_s(k,j,i)
     493                      IF ( lad_local == 0.0_wp  .AND.                          &
     494                           lad_s(k,j-1,i) > 0.0_wp )  THEN
     495                         lad_local = lad_s(k,j-1,i)
     496                      ENDIF
     497
     498                      pre_tend = 0.0_wp
     499                      pre_v = 0.0_wp
     500!
     501!--                   Calculate preliminary value (pre_tend) of the tendency
     502                      pre_tend = - cdc *                                       &
     503                                   lad_local *                                 &
     504                                   SQRT( ( 0.25_wp * ( u(k,j-1,i)   +          &
     505                                                       u(k,j-1,i+1) +          &
     506                                                       u(k,j,i)     +          &
     507                                                       u(k,j,i+1) )            &
     508                                         )**2 +                                &
     509                                         v(k,j,i)**2 +                         &
     510                                         ( 0.25_wp * ( w(k-1,j-1,i) +          &
     511                                                       w(k-1,j,i)   +          &
     512                                                       w(k,j-1,i)   +          &
     513                                                       w(k,j,i) )              &
     514                                         )**2                                  &
     515                                       ) *                                     &
     516                                   v(k,j,i)
     517
     518!
     519!--                   Calculate preliminary new velocity, based on pre_tend
     520                      pre_v = v(k,j,i) + dt_3d * pre_tend
     521!
     522!--                   Compare sign of old velocity and new preliminary velocity,
     523!--                   and in case the signs are different, limit the tendency
     524                      IF ( SIGN(pre_v,v(k,j,i)) /= pre_v )  THEN
     525                         pre_tend = - v(k,j,i) * ddt_3d
     526                      ELSE
     527                         pre_tend = pre_tend
     528                      ENDIF
     529!
     530!--                   Calculate final tendency
     531                      tend(k,j,i) = tend(k,j,i) + pre_tend
     532
    136533                   ENDDO
    137534                ENDDO
     
    143540             DO  i = nxl, nxr
    144541                DO  j = nys, nyn
    145                    DO  k = nzb_w_inner(j,i)+1, pch_index
    146                       tend(k,j,i) = tend(k,j,i) -                &
    147                                     cdc(k,j,i) * lad_w(k,j,i) *  &
    148                                     SQRT( ( ( u(k,j,i)        +  &
    149                                               u(k,j,i+1)      +  &
    150                                               u(k+1,j,i)      +  &
    151                                               u(k+1,j,i+1) )     &
    152                                             / 4.0_wp )**2     +  &
    153                                           ( ( v(k,j,i)        +  &
    154                                               v(k,j+1,i)      +  &
    155                                               v(k+1,j,i)      +  &
    156                                               v(k+1,j+1,i) )     &
    157                                             / 4.0_wp )**2     +  &
    158                                               w(k,j,i)**2 )   *  &
    159                                     w(k,j,i)
     542                   DO  k = nzb_w_inner(j,i)+1, pch_index-1
     543
     544                      pre_tend = 0.0_wp
     545                      pre_w = 0.0_wp
     546!
     547!--                   Calculate preliminary value (pre_tend) of the tendency
     548                      pre_tend = - cdc *                                       &
     549                                   (0.5_wp *                                   &
     550                                      ( lad_s(k+1,j,i) + lad_s(k,j,i) )) *     &
     551                                   SQRT( ( 0.25_wp * ( u(k,j,i)   +            &
     552                                                       u(k,j,i+1) +            &
     553                                                       u(k+1,j,i) +            &
     554                                                       u(k+1,j,i+1) )          &
     555                                         )**2 +                                &
     556                                         ( 0.25_wp * ( v(k,j,i)   +            &
     557                                                       v(k,j+1,i) +            &
     558                                                       v(k+1,j,i) +            &
     559                                                       v(k+1,j+1,i) )          &
     560                                         )**2 +                                &
     561                                         w(k,j,i)**2                           &
     562                                       ) *                                     &
     563                                   w(k,j,i)
     564!
     565!--                   Calculate preliminary new velocity, based on pre_tend
     566                      pre_w = w(k,j,i) + dt_3d * pre_tend
     567!
     568!--                   Compare sign of old velocity and new preliminary velocity,
     569!--                   and in case the signs are different, limit the tendency
     570                      IF ( SIGN(pre_w,w(k,j,i)) /= pre_w )  THEN
     571                         pre_tend = - w(k,j,i) * ddt_3d
     572                      ELSE
     573                         pre_tend = pre_tend
     574                      ENDIF
     575!
     576!--                   Calculate final tendency
     577                      tend(k,j,i) = tend(k,j,i) + pre_tend
     578
    160579                   ENDDO
    161580                ENDDO
     
    168587                DO  j = nys, nyn
    169588                   DO  k = nzb_s_inner(j,i)+1, pch_index
    170                       tend(k,j,i) = tend(k,j,i) +                   &
    171                                     ( canopy_heat_flux(k,j,i) -     &
    172                                       canopy_heat_flux(k-1,j,i) ) / &
    173                                       dzw(k)
     589                      tend(k,j,i) = tend(k,j,i) +                              &
     590                                       ( canopy_heat_flux(k,j,i) -             &
     591                                         canopy_heat_flux(k-1,j,i) ) / dzw(k)
    174592                   ENDDO
    175593                ENDDO
     
    182600                DO  j = nys, nyn
    183601                   DO  k = nzb_s_inner(j,i)+1, pch_index
    184                       tend(k,j,i) = tend(k,j,i) -                     &
    185                                     sec(k,j,i) * lad_s(k,j,i) *       &
    186                                     SQRT( ( ( u(k,j,i)        +       &
    187                                               u(k,j,i+1) )            &
    188                                             / 2.0_wp )**2     +       &
    189                                           ( ( v(k,j,i)        +       &
    190                                               v(k,j+1,i) )            &
    191                                             / 2.0_wp )**2     +       &
    192                                           ( ( w(k-1,j,i)      +       &
    193                                               w(k,j,i) )              &
    194                                             / 2.0_wp )**2 )   *       &
    195                                     ( q(k,j,i) - sls(k,j,i) )
     602                      tend(k,j,i) = tend(k,j,i) -                              &
     603                                       lsec *                                  &
     604                                       lad_s(k,j,i) *                          &
     605                                       SQRT( ( 0.5_wp * ( u(k,j,i) +           &
     606                                                          u(k,j,i+1) )         &
     607                                             )**2 +                            &
     608                                             ( 0.5_wp * ( v(k,j,i) +           &
     609                                                          v(k,j+1,i) )         &
     610                                             )**2 +                            &
     611                                             ( 0.5_wp * ( w(k-1,j,i) +         &
     612                                                          w(k,j,i) )           &
     613                                             )**2                              &
     614                                           ) *                                 &
     615                                       ( q(k,j,i) - lsc )
    196616                   ENDDO
    197617                ENDDO
     
    204624                DO  j = nys, nyn
    205625                   DO  k = nzb_s_inner(j,i)+1, pch_index
    206                       tend(k,j,i) = tend(k,j,i) -                     &
    207                                     2.0_wp * cdc(k,j,i) * lad_s(k,j,i) * &
    208                                     SQRT( ( ( u(k,j,i)              + &
    209                                               u(k,j,i+1) )            &
    210                                             / 2.0_wp )**2           + &
    211                                           ( ( v(k,j,i)              + &
    212                                               v(k,j+1,i) )            &
    213                                             / 2.0_wp )**2           + &
    214                                           ( ( w(k,j,i)              + &
    215                                               w(k+1,j,i) )            &
    216                                             / 2.0_wp )**2 )         * &
    217                                     e(k,j,i)
     626                      tend(k,j,i) = tend(k,j,i) -                              &
     627                                       2.0_wp * cdc *                          &
     628                                       lad_s(k,j,i) *                          &
     629                                       SQRT( ( 0.5_wp * ( u(k,j,i) +           &
     630                                                          u(k,j,i+1) )         &
     631                                             )**2 +                            &
     632                                             ( 0.5_wp * ( v(k,j,i) +           &
     633                                                          v(k,j+1,i) )         &
     634                                             )**2 +                            &
     635                                             ( 0.5_wp * ( w(k,j,i) +           &
     636                                                          w(k+1,j,i) )         &
     637                                             )**2                              &
     638                                           ) *                                 &
     639                                       e(k,j,i)
    218640                   ENDDO
    219641                ENDDO
    220642             ENDDO
    221                        
     643
     644
    222645          CASE DEFAULT
    223646
     
    231654
    232655!------------------------------------------------------------------------------!
    233 ! Call for grid point i,j
     656! Description:
     657! ------------
     658!-- Calculation of the tendency terms, accounting for the effect of the plant
     659!-- canopy on momentum and scalar quantities.
     660!--
     661!-- The canopy is located where the leaf area density lad_s(k,j,i) > 0.0
     662!-- (defined on scalar grid), as initialized in subroutine init_plant_canopy.
     663!-- The lad on the w-grid is vertically interpolated from the surrounding
     664!-- lad_s. The upper boundary of the canopy is defined on the w-grid at
     665!-- k = pch_index. Here, the lad is zero.
     666!--
     667!-- The canopy drag must be limited (previously accounted for by calculation of
     668!-- a limiting canopy timestep for the determination of the maximum LES timestep
     669!-- in subroutine timestep), since it is physically impossible that the canopy
     670!-- drag alone can locally change the sign of a velocity component. This
     671!-- limitation is realized by calculating preliminary tendencies and velocities.
     672!-- It is subsequently checked if the preliminary new velocity has a different
     673!-- sign than the current velocity. If so, the tendency is limited in a way that
     674!-- the velocity can at maximum be reduced to zero by the canopy drag.
     675!--
     676!--
     677!-- Call for grid point i,j
    234678!------------------------------------------------------------------------------!
    235679    SUBROUTINE plant_canopy_model_ij( i, j, component )
    236680
    237        USE arrays_3d,                                                          &
    238            ONLY:  canopy_heat_flux, cdc, dzw, e, lad_s, lad_u, lad_v, lad_w,   &
    239                   q, sec, sls, tend, u, v, w
    240681
    241682       USE control_parameters,                                                 &
    242            ONLY:  pch_index, message_string
    243 
    244        USE indices,                                                            &
    245            ONLY:  nxl, nxlu, nxr, nys, nysv, nyn, nzb_s_inner, nzb_u_inner,    &
    246                   nzb_v_inner, nzb_w_inner
     683           ONLY:  dt_3d, message_string
    247684
    248685       USE kinds
     
    250687       IMPLICIT NONE
    251688
    252        INTEGER(iwp) ::  component  !:
    253        INTEGER(iwp) ::  i          !:
    254        INTEGER(iwp) ::  j          !:
    255        INTEGER(iwp) ::  k          !:
    256 
    257 !
    258 !--    Compute drag for the three velocity components
     689       INTEGER(iwp) ::  component !: prognostic variable (u,v,w,pt,q,e)
     690       INTEGER(iwp) ::  i         !: running index
     691       INTEGER(iwp) ::  j         !: running index
     692       INTEGER(iwp) ::  k         !: running index
     693
     694       REAL(wp) ::  ddt_3d    !: inverse of the LES timestep (dt_3d)
     695       REAL(wp) ::  lad_local !: local lad value
     696       REAL(wp) ::  pre_tend  !: preliminary tendency
     697       REAL(wp) ::  pre_u     !: preliminary u-value
     698       REAL(wp) ::  pre_v     !: preliminary v-value
     699       REAL(wp) ::  pre_w     !: preliminary w-value
     700
     701
     702       ddt_3d = 1.0_wp / dt_3d
     703
     704!
     705!--    Compute drag for the three velocity components and the SGS-TKE
    259706       SELECT CASE ( component )
    260707
    261708!
    262709!--       u-component
    263        CASE ( 1 )
    264           DO  k = nzb_u_inner(j,i)+1, pch_index
    265              tend(k,j,i) = tend(k,j,i) -                     &
    266                               cdc(k,j,i) * lad_u(k,j,i) *    &   
    267                               SQRT(     u(k,j,i)**2 +        &
    268                                     ( ( v(k,j,i-1)  +        &
    269                                         v(k,j,i)    +        &
    270                                         v(k,j+1,i)  +        &
    271                                         v(k,j+1,i-1) )       &
    272                                       / 4.0_wp )**2 +        &
    273                                     ( ( w(k-1,j,i-1) +       &
    274                                         w(k-1,j,i)   +       &
    275                                         w(k,j,i-1)   +       &
    276                                         w(k,j,i) )           &
    277                                       / 4.0_wp )**2 ) *      &
    278                               u(k,j,i)
    279           ENDDO
     710          CASE ( 1 )
     711             DO  k = nzb_u_inner(j,i)+1, pch_index
     712
     713!
     714!--             In order to create sharp boundaries of the plant canopy,
     715!--             the lad on the u-grid at index (k,j,i) is equal to lad_s(k,j,i),
     716!--             rather than being interpolated from the surrounding lad_s,
     717!--             because this would yield smaller lad at the canopy boundaries
     718!--             than inside of the canopy.
     719!--             For the same reason, the lad at the rightmost(i+1)canopy
     720!--             boundary on the u-grid equals lad_s(k,j,i).
     721                lad_local = lad_s(k,j,i)
     722                IF ( lad_local == 0.0_wp  .AND.                                &
     723                     lad_s(k,j,i-1) > 0.0_wp )  THEN
     724                   lad_local = lad_s(k,j,i-1)
     725                ENDIF
     726
     727                pre_tend = 0.0_wp
     728                pre_u = 0.0_wp
     729!
     730!--             Calculate preliminary value (pre_tend) of the tendency
     731                pre_tend = - cdc *                                             &
     732                             lad_local *                                       &   
     733                             SQRT( u(k,j,i)**2 +                               &
     734                                   ( 0.25_wp * ( v(k,j,i-1)  +                 &
     735                                                 v(k,j,i)    +                 &
     736                                                 v(k,j+1,i)  +                 &
     737                                                 v(k,j+1,i-1) )                &
     738                                   )**2 +                                      &
     739                                   ( 0.25_wp * ( w(k-1,j,i-1) +                &
     740                                                 w(k-1,j,i)   +                &
     741                                                 w(k,j,i-1)   +                &
     742                                                 w(k,j,i) )                    &
     743                                   )**2                                        &
     744                                 ) *                                           &
     745                             u(k,j,i)
     746
     747!
     748!--             Calculate preliminary new velocity, based on pre_tend
     749                pre_u = u(k,j,i) + dt_3d * pre_tend
     750!
     751!--             Compare sign of old velocity and new preliminary velocity,
     752!--             and in case the signs are different, limit the tendency
     753                IF ( SIGN(pre_u,u(k,j,i)) /= pre_u )  THEN
     754                   pre_tend = - u(k,j,i) * ddt_3d
     755                ELSE
     756                   pre_tend = pre_tend
     757                ENDIF
     758!
     759!--             Calculate final tendency
     760                tend(k,j,i) = tend(k,j,i) + pre_tend
     761             ENDDO
     762
    280763
    281764!
    282765!--       v-component
    283        CASE ( 2 )
    284           DO  k = nzb_v_inner(j,i)+1, pch_index
    285              tend(k,j,i) = tend(k,j,i) -                     &
    286                               cdc(k,j,i) * lad_v(k,j,i) *    &
    287                               SQRT( ( ( u(k,j-1,i)   +       &
    288                                         u(k,j-1,i+1) +       &
    289                                         u(k,j,i)     +       &
    290                                         u(k,j,i+1) )         &
    291                                       / 4.0_wp )**2  +       &
    292                                         v(k,j,i)**2  +       &
    293                                     ( ( w(k-1,j-1,i) +       &
    294                                         w(k-1,j,i)   +       &
    295                                         w(k,j-1,i)   +       &
    296                                         w(k,j,i) )           &
    297                                       / 4.0_wp )**2 ) *      &
    298                               v(k,j,i)
    299           ENDDO
     766          CASE ( 2 )
     767             DO  k = nzb_v_inner(j,i)+1, pch_index
     768
     769!
     770!--             In order to create sharp boundaries of the plant canopy,
     771!--             the lad on the v-grid at index (k,j,i) is equal to lad_s(k,j,i),
     772!--             rather than being interpolated from the surrounding lad_s,
     773!--             because this would yield smaller lad at the canopy boundaries
     774!--             than inside of the canopy.
     775!--             For the same reason, the lad at the northmost(j+1)canopy
     776!--             boundary on the v-grid equals lad_s(k,j,i).
     777                lad_local = lad_s(k,j,i)
     778                IF ( lad_local == 0.0_wp  .AND.                                &
     779                     lad_s(k,j-1,i) > 0.0_wp )  THEN
     780                   lad_local = lad_s(k,j-1,i)
     781                ENDIF
     782
     783                pre_tend = 0.0_wp
     784                pre_v = 0.0_wp
     785!
     786!--             Calculate preliminary value (pre_tend) of the tendency
     787                pre_tend = - cdc *                                             &
     788                             lad_local *                                       &
     789                             SQRT( ( 0.25_wp * ( u(k,j-1,i)   +                &
     790                                                 u(k,j-1,i+1) +                &
     791                                                 u(k,j,i)     +                &
     792                                                 u(k,j,i+1) )                  &
     793                                   )**2 +                                      &
     794                                   v(k,j,i)**2 +                               &
     795                                   ( 0.25_wp * ( w(k-1,j-1,i) +                &
     796                                                 w(k-1,j,i)   +                &
     797                                                 w(k,j-1,i)   +                &
     798                                                 w(k,j,i) )                    &
     799                                   )**2                                        &
     800                                 ) *                                           &
     801                             v(k,j,i)
     802
     803!
     804!--             Calculate preliminary new velocity, based on pre_tend
     805                pre_v = v(k,j,i) + dt_3d * pre_tend
     806!
     807!--             Compare sign of old velocity and new preliminary velocity,
     808!--             and in case the signs are different, limit the tendency
     809                IF ( SIGN(pre_v,v(k,j,i)) /= pre_v )  THEN
     810                   pre_tend = - v(k,j,i) * ddt_3d
     811                ELSE
     812                   pre_tend = pre_tend
     813                ENDIF
     814!
     815!--             Calculate final tendency
     816                tend(k,j,i) = tend(k,j,i) + pre_tend
     817             ENDDO
     818
    300819
    301820!
    302821!--       w-component
    303        CASE ( 3 )
    304           DO  k = nzb_w_inner(j,i)+1, pch_index
    305              tend(k,j,i) = tend(k,j,i) -                     &
    306                               cdc(k,j,i) * lad_w(k,j,i) *    &
    307                               SQRT( ( ( u(k,j,i)    +        & 
    308                                         u(k,j,i+1)  +        &
    309                                         u(k+1,j,i)  +        &
    310                                         u(k+1,j,i+1) )       &
    311                                       / 4.0_wp )**2 +        &
    312                                     ( ( v(k,j,i)    +        &
    313                                         v(k,j+1,i)  +        &
    314                                         v(k+1,j,i)  +        &
    315                                         v(k+1,j+1,i) )       &
    316                                       / 4.0_wp )**2 +        &
    317                                         w(k,j,i)**2 ) *      &
    318                               w(k,j,i)
    319    
    320           ENDDO
     822          CASE ( 3 )
     823             DO  k = nzb_w_inner(j,i)+1, pch_index-1
     824
     825                pre_tend = 0.0_wp
     826                pre_w = 0.0_wp
     827!
     828!--             Calculate preliminary value (pre_tend) of the tendency
     829                pre_tend = - cdc *                                             &
     830                             (0.5_wp *                                         &
     831                                ( lad_s(k+1,j,i) + lad_s(k,j,i) )) *           &
     832                             SQRT( ( 0.25_wp * ( u(k,j,i)    +                 & 
     833                                                 u(k,j,i+1)  +                 &
     834                                                 u(k+1,j,i)  +                 &
     835                                                 u(k+1,j,i+1) )                &
     836                                   )**2 +                                      &
     837                                   ( 0.25_wp * ( v(k,j,i)    +                 &
     838                                                 v(k,j+1,i)  +                 &
     839                                                 v(k+1,j,i)  +                 &
     840                                                 v(k+1,j+1,i) )                &
     841                                   )**2 +                                      &
     842                                   w(k,j,i)**2                                 &
     843                                 ) *                                           &
     844                             w(k,j,i)
     845!
     846!--             Calculate preliminary new velocity, based on pre_tend
     847                pre_w = w(k,j,i) + dt_3d * pre_tend
     848!
     849!--             Compare sign of old velocity and new preliminary velocity,
     850!--             and in case the signs are different, limit the tendency
     851                IF ( SIGN(pre_w,w(k,j,i)) /= pre_w )  THEN
     852                   pre_tend = - w(k,j,i) * ddt_3d
     853                ELSE
     854                   pre_tend = pre_tend
     855                ENDIF
     856!
     857!--             Calculate final tendency
     858                tend(k,j,i) = tend(k,j,i) + pre_tend
     859             ENDDO
    321860
    322861!
     
    324863          CASE ( 4 )
    325864             DO  k = nzb_s_inner(j,i)+1, pch_index
    326                 tend(k,j,i) = tend(k,j,i) +                   &
    327                               ( canopy_heat_flux(k,j,i) -     &
    328                                 canopy_heat_flux(k-1,j,i) ) / &
    329                                 dzw(k)
     865                tend(k,j,i) = tend(k,j,i) +                                    &
     866                                 ( canopy_heat_flux(k,j,i) -                   &
     867                                   canopy_heat_flux(k-1,j,i) ) / dzw(k)
    330868             ENDDO
    331869
     
    335873          CASE ( 5 )
    336874             DO  k = nzb_s_inner(j,i)+1, pch_index
    337                 tend(k,j,i) = tend(k,j,i) -                     &
    338                               sec(k,j,i) * lad_s(k,j,i) *       &
    339                               SQRT( ( ( u(k,j,i)        +       &
    340                                         u(k,j,i+1) )            &
    341                                       / 2.0_wp )**2     +       &
    342                                     ( ( v(k,j,i)        +       &
    343                                         v(k,j+1,i) )            &
    344                                       / 2.0_wp )**2     +       &
    345                                     ( ( w(k-1,j,i)      +       &
    346                                         w(k,j,i) )              &
    347                                       / 2.0_wp )**2 )   *       &
    348                               ( q(k,j,i) - sls(k,j,i) )
     875                tend(k,j,i) = tend(k,j,i) -                                    &
     876                                 lsec *                                        &
     877                                 lad_s(k,j,i) *                                &
     878                                 SQRT( ( 0.5_wp * ( u(k,j,i) +                 &
     879                                                    u(k,j,i+1) )               &
     880                                       )**2  +                                 &
     881                                       ( 0.5_wp * ( v(k,j,i) +                 &
     882                                                    v(k,j+1,i) )               &
     883                                       )**2 +                                  &
     884                                       ( 0.5_wp * ( w(k-1,j,i) +               &
     885                                                    w(k,j,i) )                 &
     886                                       )**2                                    &
     887                                     ) *                                       &
     888                                 ( q(k,j,i) - lsc )
    349889             ENDDO   
    350890
    351891!
    352892!--       sgs-tke
    353        CASE ( 6 )
    354           DO  k = nzb_s_inner(j,i)+1, pch_index   
    355              tend(k,j,i) = tend(k,j,i) -                        &
    356                               2.0_wp * cdc(k,j,i) * lad_s(k,j,i) * &
    357                               SQRT( ( ( u(k,j,i)           +    &
    358                                         u(k,j,i+1) )            &
    359                                       / 2.0_wp )**2        +    & 
    360                                     ( ( v(k,j,i)           +    &
    361                                         v(k,j+1,i) )            &
    362                                       / 2.0_wp )**2        +    &
    363                                     ( ( w(k,j,i)           +    &
    364                                         w(k+1,j,i) )            &
    365                                       / 2.0_wp )**2 )      *    &
    366                               e(k,j,i)
    367           ENDDO
     893          CASE ( 6 )
     894             DO  k = nzb_s_inner(j,i)+1, pch_index   
     895                tend(k,j,i) = tend(k,j,i) -                                    &
     896                                 2.0_wp * cdc *                                &
     897                                 lad_s(k,j,i) *                                &
     898                                 SQRT( ( 0.5_wp * ( u(k,j,i) +                 &
     899                                                    u(k,j,i+1) )               &
     900                                       )**2 +                                  & 
     901                                       ( 0.5_wp * ( v(k,j,i) +                 &
     902                                                    v(k,j+1,i) )               &
     903                                       )**2 +                                  &
     904                                       ( 0.5_wp * ( w(k,j,i) +                 &
     905                                                    w(k+1,j,i) )               &
     906                                       )**2                                    &
     907                                     ) *                                       &
     908                                 e(k,j,i)
     909             ENDDO
    368910
    369911       CASE DEFAULT
  • palm/trunk/SOURCE/prognostic_equations.f90

    r1410 r1484  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Changes due to new module structure of the plant canopy model:
     23!   parameters cthf and plant_canopy moved to module plant_canopy_model_mod.
     24! Removed double-listing of use_upstream_for_tke in ONLY-list of module
     25! control_parameters
    2326!
    2427! Former revisions:
     
    163166    USE control_parameters,                                                    &
    164167        ONLY:  call_microphysics_at_all_substeps, cloud_physics,               &
    165                constant_diffusion, cthf, dp_external,                          &
     168               constant_diffusion, dp_external,                                &
    166169               dp_level_ind_b, dp_smooth_factor, dpdxy, dt_3d, humidity,       &
    167170               icloud_scheme, inflow_l, intermediate_timestep_count,           &
    168171               intermediate_timestep_count_max, large_scale_forcing,           &
    169172               large_scale_subsidence, neutral, nudging, ocean, outflow_l,     &
    170                outflow_s, passive_scalar, plant_canopy, precipitation,         &
     173               outflow_s, passive_scalar, precipitation,                       &
    171174               prho_reference, prho_reference, prho_reference, pt_reference,   &
    172175               pt_reference, pt_reference, radiation, scalar_advec,            &
    173176               scalar_advec, simulated_time, sloping_surface, timestep_scheme, &
    174177               tsc, use_subsidence_tendencies, use_upstream_for_tke,           &
    175                use_upstream_for_tke, use_upstream_for_tke, wall_heatflux,      &
     178               wall_heatflux,                                                  &
    176179               wall_nrflux, wall_qflux, wall_qflux, wall_qflux, wall_qrflux,   &
    177180               wall_salinityflux, ws_scheme_mom, ws_scheme_sca
     
    257260
    258261    USE plant_canopy_model_mod,                                                &
    259         ONLY:  plant_canopy_model
     262        ONLY:  cthf, plant_canopy, plant_canopy_model
    260263
    261264    USE production_e_mod,                                                      &
  • palm/trunk/SOURCE/read_var_list.f90

    r1323 r1484  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Changes in the course of the canopy-model modularization:
     23!   parameters alpha_lad, beta_lad, lai_beta added,
     24!   module plant_canopy_model_mod added,
     25!   drag_coefficient, leaf_surface_concentration and scalar_exchange_coefficient
     26!   renamed to canopy_drag_coeff, leaf_surface_conc and leaf_scalar_exch_coeff
    2327!
    2428! Former revisions:
     
    117121
    118122    USE arrays_3d,                                                             &
    119         ONLY:  inflow_damping_factor, lad, mean_inflow_profiles, pt_init,      &
     123        ONLY:  inflow_damping_factor, mean_inflow_profiles, pt_init,           &
    120124               q_init, ref_state, sa_init, u_init, ug, v_init, vg
    121125
     
    141145
    142146    USE pegrid
     147
     148    USE plant_canopy_model_mod,                                                &
     149        ONLY:  alpha_lad, beta_lad, canopy_drag_coeff, canopy_mode, cthf, lad, &
     150               lad_surface, lad_vertical_gradient, lad_vertical_gradient_level,&
     151               lad_vertical_gradient_level_ind,                                &
     152               lai_beta, leaf_scalar_exch_coeff, leaf_surface_conc, pch_index, &
     153               plant_canopy
    143154
    144155    USE statistics,                                                            &
     
    242253          CASE ( 'advected_distance_y' )
    243254             READ ( 13 )  advected_distance_y
     255          CASE ( 'alpha_lad' )
     256             READ ( 13 )  alpha_lad
    244257          CASE ( 'alpha_surface' )
    245258             READ ( 13 )  alpha_surface
     
    282295          CASE ( 'bc_uv_t' )
    283296             READ ( 13 )  bc_uv_t
     297          CASE ( 'beta_lad' )
     298             READ ( 13 )  beta_lad
    284299          CASE ( 'bottom_salinityflux' )
    285300             READ ( 13 )  bottom_salinityflux
     
    296311          CASE ( 'call_psolver_at_all_substeps' )
    297312             READ ( 13 )  call_psolver_at_all_substeps
     313          CASE ( 'canopy_drag_coeff' )
     314             READ ( 13 )  canopy_drag_coeff
    298315          CASE ( 'canopy_mode' )
    299316             READ ( 13 )  canopy_mode
     
    354371          CASE ( 'dpdxy' )
    355372             READ ( 13 )  dpdxy
    356           CASE ( 'drag_coefficient' )
    357              READ ( 13 )  drag_coefficient
    358373          CASE ( 'drizzle' )
    359374             READ ( 13 )  drizzle
     
    419434          CASE ( 'lad_vertical_gradient_level_in' )
    420435             READ ( 13 )  lad_vertical_gradient_level_ind
     436          CASE ( 'lai_beta' )
     437             READ ( 13 )  lai_beta
    421438          CASE ( 'large_scale_forcing' )
    422439             READ ( 13 )  large_scale_forcing
    423440           CASE ( 'large_scale_subsidence' )
    424441             READ ( 13 )  large_scale_subsidence
    425           CASE ( 'leaf_surface_concentration' )
    426              READ ( 13 )  leaf_surface_concentration
     442          CASE ( 'leaf_scalar_exch_coeff' )
     443             READ ( 13 )  leaf_scalar_exch_coeff
     444          CASE ( 'leaf_surface_conc' )
     445             READ ( 13 )  leaf_surface_conc
    427446          CASE ( 'limiter_sedimentation' )
    428447             READ ( 13 )  limiter_sedimentation
     
    558577          CASE ( 'scalar_advec' )
    559578             READ ( 13 )  scalar_advec
    560           CASE ( 'scalar_exchange_coefficient' )
    561              READ ( 13 )  scalar_exchange_coefficient
    562579          CASE ( 'simulated_time' )
    563580             READ ( 13 )  simulated_time
     
    731748
    732749    USE arrays_3d,                                                             &
    733         ONLY:  inflow_damping_factor, lad, mean_inflow_profiles, pt_init,      &
     750        ONLY:  inflow_damping_factor, mean_inflow_profiles, pt_init,           &
    734751               q_init, ref_state, sa_init, u_init, ug, v_init, vg
    735752
  • palm/trunk/SOURCE/timestep.f90

    r1343 r1484  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Changes due to new module structure of the plant canopy model:
     23!   calculations and parameters related to the plant canopy model removed
     24!   (the limitation of the canopy drag, i.e. that the canopy drag itself should
     25!   not change the sign of the velocity components, is now assured for in the
     26!   calculation of the canopy tendency terms in subroutine plant_canopy_model)
    2327!
    2428! Former revisions:
     
    7680
    7781    USE arrays_3d,                                                             &
    78         ONLY:  cdc, dzu, dzw, kh, km, lad_u, lad_v, lad_w, u, v, w
     82        ONLY:  dzu, dzw, kh, km, u, v, w
    7983
    8084    USE cloud_parameters,                                                      &
     
    8387    USE control_parameters,                                                    &
    8488        ONLY:  cfl_factor, coupling_mode, dt_3d, dt_fixed, dt_max,             &
    85                galilei_transformation, old_dt, plant_canopy, message_string,   &
     89               galilei_transformation, old_dt, message_string,                 &
    8690               stop_dt, terminate_coupled, terminate_coupled_remote,           &
    8791               timestep_reason, u_gtrans, use_ug_for_galilei_tr, v_gtrans
     
    115119    REAL(wp) ::  dt_diff           !:
    116120    REAL(wp) ::  dt_diff_l         !:
    117     REAL(wp) ::  dt_plant_canopy   !:
    118     REAL(wp) ::  dt_plant_canopy_l !:
    119     REAL(wp) ::  dt_plant_canopy_u !:
    120     REAL(wp) ::  dt_plant_canopy_v !:
    121     REAL(wp) ::  dt_plant_canopy_w !:
    122121    REAL(wp) ::  dt_u              !:
    123122    REAL(wp) ::  dt_u_l            !:
     
    413412
    414413!
    415 !--    Additional timestep criterion with plant canopies:
    416 !--    it is not allowed to extract more than the available momentum
    417        IF ( plant_canopy ) THEN
    418 
    419           dt_plant_canopy_l = 0.0_wp
    420           DO  i = nxl, nxr
    421              DO  j = nys, nyn
    422                 DO k = nzb+1, nzt
    423                    dt_plant_canopy_u = cdc(k,j,i) * lad_u(k,j,i) *  &
    424                                        SQRT(     u(k,j,i)**2     +  &
    425                                              ( ( v(k,j,i-1)      +  &
    426                                                  v(k,j,i)        +  &
    427                                                  v(k,j+1,i)      +  &
    428                                                  v(k,j+1,i-1) )     &
    429                                                / 4.0_wp )**2     +  &
    430                                              ( ( w(k-1,j,i-1)    +  &
    431                                                  w(k-1,j,i)      +  &
    432                                                  w(k,j,i-1)      +  &
    433                                                  w(k,j,i) )         &
    434                                                  / 4.0_wp )**2 )
    435                    IF ( dt_plant_canopy_u > dt_plant_canopy_l ) THEN
    436                       dt_plant_canopy_l = dt_plant_canopy_u 
    437                    ENDIF
    438                    dt_plant_canopy_v = cdc(k,j,i) * lad_v(k,j,i) *  &
    439                                        SQRT( ( ( u(k,j-1,i)      +  &
    440                                                  u(k,j-1,i+1)    +  &
    441                                                  u(k,j,i)        +  &
    442                                                  u(k,j,i+1) )       &
    443                                                / 4.0_wp )**2     +  &
    444                                                  v(k,j,i)**2     +  &
    445                                              ( ( w(k-1,j-1,i)    +  &
    446                                                  w(k-1,j,i)      +  &
    447                                                  w(k,j-1,i)      +  &
    448                                                  w(k,j,i) )         &
    449                                                  / 4.0_wp )**2 )
    450                    IF ( dt_plant_canopy_v > dt_plant_canopy_l ) THEN
    451                       dt_plant_canopy_l = dt_plant_canopy_v
    452                    ENDIF                   
    453                    dt_plant_canopy_w = cdc(k,j,i) * lad_w(k,j,i) *  &
    454                                        SQRT( ( ( u(k,j,i)        +  &
    455                                                  u(k,j,i+1)      +  &
    456                                                  u(k+1,j,i)      +  &
    457                                                  u(k+1,j,i+1) )     &
    458                                                / 4.0_wp )**2     +  &
    459                                              ( ( v(k,j,i)        +  &
    460                                                  v(k,j+1,i)      +  &
    461                                                  v(k+1,j,i)      +  &
    462                                                  v(k+1,j+1,i) )     &
    463                                                / 4.0_wp )**2        +  &
    464                                                  w(k,j,i)**2 )     
    465                    IF ( dt_plant_canopy_w > dt_plant_canopy_l ) THEN
    466                       dt_plant_canopy_l = dt_plant_canopy_w
    467                    ENDIF
    468                 ENDDO
    469              ENDDO
    470           ENDDO
    471 
    472           IF ( dt_plant_canopy_l > 0.0_wp ) THEN
    473 !
    474 !--          Invert dt_plant_canopy_l and apply a security timestep factor 0.1
    475              dt_plant_canopy_l = 0.1_wp / dt_plant_canopy_l
    476           ELSE
    477 !
    478 !--          In case of inhomogeneous plant canopy, some processors may have no
    479 !--          canopy at all. Then use dt_max as dummy instead.
    480              dt_plant_canopy_l = dt_max
    481           ENDIF
    482 
    483 !
    484 !--       Determine the global minumum
    485 #if defined( __parallel )
    486           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    487           CALL MPI_ALLREDUCE( dt_plant_canopy_l, dt_plant_canopy, 1, MPI_REAL, &
    488                               MPI_MIN, comm2d, ierr )
    489 #else
    490           dt_plant_canopy = dt_plant_canopy_l
    491 #endif
    492 
    493        ELSE
    494 !
    495 !--       Use dt_diff as dummy value to avoid extra IF branches further below
    496           dt_plant_canopy = dt_diff
    497 
    498        ENDIF
    499 
    500 !
    501414!--    The time step is the minimum of the 3-4 components and the diffusion time
    502415!--    step minus a reduction (cfl_factor) to be on the safe side.
    503416!--    The time step must not exceed the maximum allowed value.
    504        dt_3d = cfl_factor * MIN( dt_diff, dt_plant_canopy, dt_u, dt_v, dt_w,   &
     417       dt_3d = cfl_factor * MIN( dt_diff, dt_u, dt_v, dt_w,   &
    505418                                 dt_precipitation )
    506419       dt_3d = MIN( dt_3d, dt_max )
     
    508421!
    509422!--    Remember the restricting time step criterion for later output.
    510        IF ( MIN( dt_u, dt_v, dt_w ) < MIN( dt_diff, dt_plant_canopy ) )  THEN
     423       IF ( MIN( dt_u, dt_v, dt_w ) < dt_diff )  THEN
    511424          timestep_reason = 'A'
    512        ELSEIF ( dt_plant_canopy < dt_diff )  THEN
    513           timestep_reason = 'C'
    514425       ELSE
    515426          timestep_reason = 'D'
     
    528439               '&dt_w            = ', dt_w, ' s',                             &
    529440               '&dt_diff         = ', dt_diff, ' s',                          &
    530                '&dt_plant_canopy = ', dt_plant_canopy, ' s',                  &
    531441               '&u_max           = ', u_max, ' m/s   k=', u_max_ijk(1),       &
    532442               '  j=', u_max_ijk(2), '  i=', u_max_ijk(3),                    &
  • palm/trunk/SOURCE/user_init_plant_canopy.f90

    r1354 r1484  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Changes in the course of the canopy-model modularization:
     23!   module plant_canopy_model_mod added,
     24!   definition of array cdc (canopy drag coefficient) removed, since it is now
     25!   defined purely as a single constant value (see module plant_canopy_model)
    2326!
    2427! Former revisions:
     
    5760   
    5861    USE kinds
     62
     63    USE plant_canopy_model_mod
    5964   
    6065    USE user
     
    6267    IMPLICIT NONE
    6368
    64     INTEGER(iwp) :: i   !:
    65     INTEGER(iwp) :: j   !:
     69    INTEGER(iwp) :: i   !: running index
     70    INTEGER(iwp) :: j   !: running index
    6671
    6772!
     
    6974
    7075!
    71 !-- Set the 3D-arrays lad_s and cdc for user defined canopies
     76!-- Set the 3D-array lad_s for user defined canopies
    7277    SELECT CASE ( TRIM( canopy_mode ) )
    7378
     
    7883       CASE ( 'user_defined_canopy_1' )
    7984!
    80 !--       Here the user can define his own topography.
    81 !--       The following lines contain an example in that the
    82 !--       plant canopy extends only over the second half of the
    83 !--       model domain along x
     85!--       Here the user can define his own forest topography.
     86!--       The following lines contain an example, where the plant canopy extends
     87!--       only over the second half of the model domain along x.
     88!--       Attention: DO-loops have to include the ghost points (nxlg-nxrg,
     89!--       nysg-nyng), because no exchange of ghost point information is intended,
     90!--       in order to minimize communication between CPUs
    8491!          DO  i = nxlg, nxrg
    8592!             IF ( i >= INT(nx+1/2) ) THEN
    8693!                DO  j = nysg, nyng
    8794!                   lad_s(:,j,i) = lad(:)
    88 !                   cdc(:,j,i)   = drag_coefficient
    8995!                ENDDO
    9096!             ELSE
    9197!                lad_s(:,:,i) = 0.0_wp
    92 !                cdc(:,:,i)   = 0.0_wp
    9398!             ENDIF
    9499!          ENDDO
     100!
    95101!--       After definition, please
    96102!--       remove the following three lines!
  • palm/trunk/SOURCE/write_var_list.f90

    r1353 r1484  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Changes in the course of the canopy-model modularization:
     23!   parameters alpha_lad, beta_lad, lai_beta added,
     24!   module plant_canopy_model_mod added,
     25!   drag_coefficient, leaf_surface_concentration and scalar_exchange_coefficient
     26!   renamed to canopy_drag_coeff, leaf_surface_conc and leaf_scalar_exch_coeff
    2327!
    2428! Former revisions:
     
    104108
    105109    USE arrays_3d,                                                             &
    106         ONLY:  inflow_damping_factor, lad, mean_inflow_profiles, pt_init,      &
     110        ONLY:  inflow_damping_factor, mean_inflow_profiles, pt_init,           &
    107111               q_init, ref_state, sa_init, u_init, ug, v_init, vg
    108112   
     
    128132   
    129133    USE pegrid
     134
     135    USE plant_canopy_model_mod,                                                &
     136        ONLY:  alpha_lad, beta_lad, canopy_drag_coeff, canopy_mode, cthf, lad, &
     137               lad_surface, lad_vertical_gradient, lad_vertical_gradient_level,&
     138               lad_vertical_gradient_level_ind,                                &
     139               lai_beta, leaf_scalar_exch_coeff, leaf_surface_conc, pch_index, &
     140               plant_canopy
    130141   
    131142    USE statistics,                                                            &
     
    163174    WRITE ( 14 )  'advected_distance_y           '
    164175    WRITE ( 14 )  advected_distance_y
     176    WRITE ( 14 )  'alpha_lad                     '
     177    WRITE ( 14 )  alpha_lad
    165178    WRITE ( 14 )  'alpha_surface                 '
    166179    WRITE ( 14 )  alpha_surface
     
    203216    WRITE ( 14 )  'bc_uv_t                       '
    204217    WRITE ( 14 )  bc_uv_t
     218    WRITE ( 14 )  'beta_lad                      '
     219    WRITE ( 14 )  beta_lad
    205220    WRITE ( 14 )  'bottom_salinityflux           '
    206221    WRITE ( 14 )  bottom_salinityflux
     
    217232    WRITE ( 14 )  'call_psolver_at_all_substeps  '
    218233    WRITE ( 14 )  call_psolver_at_all_substeps
     234    WRITE ( 14 )  'canopy_drag_coeff              '
     235    WRITE ( 14 )  canopy_drag_coeff
    219236    WRITE ( 14 )  'canopy_mode                   '
    220237    WRITE ( 14 )  canopy_mode
     
    275292    WRITE ( 14 )  'dpdxy                         '
    276293    WRITE ( 14 )  dpdxy
    277     WRITE ( 14 )  'drag_coefficient              '
    278     WRITE ( 14 )  drag_coefficient
    279294    WRITE ( 14 )  'drizzle                       '
    280295    WRITE ( 14 )  drizzle
     
    339354    WRITE ( 14 )  'lad_vertical_gradient_level_in'
    340355    WRITE ( 14 )  lad_vertical_gradient_level_ind
     356    WRITE ( 14 )  'lai_beta                      '
     357    WRITE ( 14 )  lai_beta
    341358    WRITE ( 14 )  'large_scale_forcing           '
    342359    WRITE ( 14 )  large_scale_forcing
    343360    WRITE ( 14 )  'large_scale_subsidence        '
    344361    WRITE ( 14 )  large_scale_subsidence
    345     WRITE ( 14 )  'leaf_surface_concentration    '
    346     WRITE ( 14 )  leaf_surface_concentration
     362    WRITE ( 14 )  'leaf_scalar_exch_coeff   '
     363    WRITE ( 14 )  leaf_scalar_exch_coeff
     364    WRITE ( 14 )  'leaf_surface_conc    '
     365    WRITE ( 14 )  leaf_surface_conc
    347366    WRITE ( 14 )  'limiter_sedimentation         '
    348367    WRITE ( 14 )  limiter_sedimentation
     
    475494    WRITE ( 14 )  'scalar_advec                  '
    476495    WRITE ( 14 )  scalar_advec
    477     WRITE ( 14 )  'scalar_exchange_coefficient   '
    478     WRITE ( 14 )  scalar_exchange_coefficient
    479496    WRITE ( 14 )  'simulated_time                '
    480497    WRITE ( 14 )  simulated_time
Note: See TracChangeset for help on using the changeset viewer.