Ignore:
Timestamp:
Dec 14, 2017 5:12:51 PM (4 years ago)
Author:
kanani
Message:

Merge of branch palm4u into trunk

Location:
palm/trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk

  • palm/trunk/SOURCE

  • palm/trunk/SOURCE/plant_canopy_model_mod.f90

    r2669 r2696  
    11!> @file plant_canopy_model_mod.f90
    22!------------------------------------------------------------------------------!
    3 ! This file is part of PALM.
     3! This file is part of the PALM model system.
    44!
    55! PALM is free software: you can redistribute it and/or modify it under the
     
    2525! -----------------
    2626! $Id$
     27! Bugfix for vertical loop index pch_index in case of Netcdf input
     28! Introduce 2D index array incorporate canopy top index
     29! Check if canopy on top of topography do not exceed vertical dimension
     30! Add check for canopy_mode in case of Netcdf input.
     31! Enable _FillValue output for 3d quantities
     32! Bugfix in reading of PIDS leaf area density (MS)
     33!
     34! 2669 2017-12-06 16:03:27Z raasch
    2735! coupling_char removed
    2836!
     
    151159    CHARACTER (LEN=20)   ::  canopy_mode = 'block' !< canopy coverage
    152160
    153     INTEGER(iwp) ::  pch_index = 0                 !< plant canopy height/top index
    154     INTEGER(iwp) ::                                                            &
    155        lad_vertical_gradient_level_ind(10) = -9999 !< lad-profile levels (index)
     161    INTEGER(iwp) ::  pch_index = 0                               !< plant canopy height/top index
     162    INTEGER(iwp) ::  lad_vertical_gradient_level_ind(10) = -9999 !< lad-profile levels (index)
     163
     164    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  pch_index_ji   !< local plant canopy top
    156165
    157166    LOGICAL ::  calc_beta_lad_profile = .FALSE. !< switch for calc. of lad from beta func.
    158167    LOGICAL ::  plant_canopy = .FALSE.          !< switch for use of canopy model
    159 
    160     REAL(wp) ::  alpha_lad = 9999999.9_wp   !< coefficient for lad calculation
    161     REAL(wp) ::  beta_lad = 9999999.9_wp    !< coefficient for lad calculation
    162     REAL(wp) ::  canopy_drag_coeff = 0.0_wp !< canopy drag coefficient (parameter)
    163     REAL(wp) ::  cdc = 0.0_wp               !< canopy drag coeff. (abbreviation used in equations)
    164     REAL(wp) ::  cthf = 0.0_wp              !< canopy top heat flux
    165     REAL(wp) ::  dt_plant_canopy = 0.0_wp   !< timestep account. for canopy drag
    166     REAL(wp) ::  ext_coef = 0.6_wp          !< extinction coefficient
    167     REAL(wp) ::  lad_surface = 0.0_wp       !< lad surface value
    168     REAL(wp) ::  lai_beta = 0.0_wp          !< leaf area index (lai) for lad calc.
    169     REAL(wp) ::                                                                &
    170        leaf_scalar_exch_coeff = 0.0_wp      !< canopy scalar exchange coeff.
    171     REAL(wp) ::                                                                &
    172        leaf_surface_conc = 0.0_wp           !< leaf surface concentration
    173     REAL(wp) ::  lsec = 0.0_wp              !< leaf scalar exchange coeff.
    174     REAL(wp) ::  lsc = 0.0_wp               !< leaf surface concentration
    175 
    176     REAL(wp) ::                                                                &
    177        lad_vertical_gradient(10) = 0.0_wp              !< lad gradient
    178     REAL(wp) ::                                                                &
    179        lad_vertical_gradient_level(10) = -9999999.9_wp !< lad-prof. levels (in m)
     168    LOGICAL ::  usm_lad_rma = .TRUE.            !< use MPI RMA to access LAD for raytracing (instead of global array)
     169
     170    REAL(wp) ::  alpha_lad = 9999999.9_wp         !< coefficient for lad calculation
     171    REAL(wp) ::  beta_lad = 9999999.9_wp          !< coefficient for lad calculation
     172    REAL(wp) ::  canopy_drag_coeff = 0.0_wp       !< canopy drag coefficient (parameter)
     173    REAL(wp) ::  cdc = 0.0_wp                     !< canopy drag coeff. (abbreviation used in equations)
     174    REAL(wp) ::  cthf = 0.0_wp                    !< canopy top heat flux
     175    REAL(wp) ::  dt_plant_canopy = 0.0_wp         !< timestep account. for canopy drag
     176    REAL(wp) ::  ext_coef = 0.6_wp                !< extinction coefficient
     177    REAL(wp) ::  lad_surface = 0.0_wp             !< lad surface value
     178    REAL(wp) ::  lai_beta = 0.0_wp                !< leaf area index (lai) for lad calc.
     179    REAL(wp) ::  leaf_scalar_exch_coeff = 0.0_wp  !< canopy scalar exchange coeff.
     180    REAL(wp) ::  leaf_surface_conc = 0.0_wp       !< leaf surface concentration
     181    REAL(wp) ::  lsec = 0.0_wp                    !< leaf scalar exchange coeff.
     182    REAL(wp) ::  lsc = 0.0_wp                     !< leaf surface concentration
     183    REAL(wp) ::  prototype_lad                    !< prototype leaf area density for computing effective optical depth
     184
     185    REAL(wp) ::  lad_vertical_gradient(10) = 0.0_wp              !< lad gradient
     186    REAL(wp) ::  lad_vertical_gradient_level(10) = -9999999.9_wp !< lad-prof. levels (in m)
    180187
    181188    REAL(wp), DIMENSION(:), ALLOCATABLE ::  lad            !< leaf area density
     
    201208!-- Public variables and constants
    202209    PUBLIC pc_heating_rate, canopy_mode, cthf, dt_plant_canopy, lad, lad_s,   &
    203            pch_index, plant_canopy
     210           pch_index, plant_canopy, prototype_lad, usm_lad_rma
    204211           
    205212
     
    288295
    289296       USE control_parameters,                                                 &
    290            ONLY: cloud_physics, message_string, microphysics_seifert
     297           ONLY: cloud_physics, coupling_char, message_string,                 &
     298                 microphysics_seifert
     299
     300       USE netcdf_data_input_mod,                                              &
     301           ONLY:  input_file_static, input_pids_static
    291302                 
    292303   
     
    331342          CALL message( 'check_parameters', 'PA0360', 1, 2, 0, 6, 0 )
    332343       ENDIF
     344!
     345!--    If dynamic input file is used, canopy need to be read from file
     346       IF ( input_pids_static  .AND.                                           &
     347            TRIM( canopy_mode ) /= 'read_from_file_3d' )  THEN
     348          message_string = 'Usage of dynamic input file ' //                   &
     349                           TRIM( input_file_static ) //                        &
     350                           TRIM( coupling_char ) // ' requires ' //            &
     351                           'canopy_mode = read_from_file_3d'
     352          CALL message( 'check_parameters', 'PA0999', 1, 2, 0, 6, 0 )
     353       ENDIF
    333354
    334355 
     
    342363!> Subroutine defining 3D output variables
    343364!------------------------------------------------------------------------------!
    344  SUBROUTINE pcm_data_output_3d( av, variable, found, local_pf )
     365 SUBROUTINE pcm_data_output_3d( av, variable, found, local_pf, fill_value )
    345366 
    346367    USE control_parameters,                                                    &
     
    356377    CHARACTER (LEN=*) ::  variable !<
    357378
    358     INTEGER(iwp) ::  av    !<
    359     INTEGER(iwp) ::  i     !<
    360     INTEGER(iwp) ::  j     !<
    361     INTEGER(iwp) ::  k     !<
     379    INTEGER(iwp) ::  av     !<
     380    INTEGER(iwp) ::  i      !<
     381    INTEGER(iwp) ::  j      !<
     382    INTEGER(iwp) ::  k      !<
     383    INTEGER(iwp) ::  k_topo !< topography top index
    362384
    363385    LOGICAL      ::  found !<
    364386
     387    REAL(wp)     ::  fill_value
    365388    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb:nz_do3d) ::  local_pf !<
    366389
     
    368391    found = .TRUE.
    369392
     393    local_pf = REAL( fill_value, KIND = 4 )
    370394
    371395    SELECT CASE ( TRIM( variable ) )
     
    375399            DO  i = nxl, nxr
    376400               DO  j = nys, nyn
    377                   DO  k = nzb_s_inner(j,i), nz_do3d
    378                      local_pf(i,j,k) = pc_heating_rate(k-nzb_s_inner(j,i),j,i)
    379                   ENDDO
     401                  IF ( pch_index_ji(j,i) /= 0 )  THEN
     402                     k_topo = get_topography_top_index( j, i, 's' )
     403                     DO  k = k_topo, k_topo + pch_index_ji(j,i)
     404                        local_pf(i,j,k) = pc_heating_rate(k-k_topo,j,i)
     405                     ENDDO
     406                  ENDIF
    380407               ENDDO
    381408            ENDDO
     
    387414            DO  i = nxl, nxr
    388415               DO  j = nys, nyn
    389                   DO  k = nzb_s_inner(j,i), nz_do3d
    390                      local_pf(i,j,k) = lad_s(k-nzb_s_inner(j,i),j,i)
    391                   ENDDO
     416                  IF ( pch_index_ji(j,i) /= 0 )  THEN
     417                     k_topo = get_topography_top_index( j, i, 's' )
     418                     DO  k = k_topo, k_topo + pch_index_ji(j,i)
     419                        local_pf(i,j,k) = lad_s(k-k_topo,j,i)
     420                     ENDDO
     421                  ENDIF
    392422               ENDDO
    393423            ENDDO
     
    570600                 passive_scalar, urban_surface
    571601
     602       USE netcdf_data_input_mod,                                              &
     603           ONLY:  leaf_area_density_f
     604
    572605       USE surface_mod,                                                        &
    573606           ONLY: surf_def_h, surf_lsm_h, surf_usm_h
     
    716749          CASE ( 'read_from_file_3d' )
    717750!
     751!--          Initialize LAD with data from file. If LAD is given in NetCDF file,
     752!--          use these values, else take LAD profiles from ASCII file.
     753!--          Please note, in NetCDF file LAD is only given up to the maximum
     754!--          canopy top, indicated by leaf_area_density_f%nz. 
     755             lad_s = 0.0_wp
     756             IF ( leaf_area_density_f%from_file )  THEN
     757!
     758!--             Set also pch_index, used to be the upper bound of the vertical
     759!--             loops. Therefore, use the global top of the canopy layer.
     760                pch_index = leaf_area_density_f%nz - 1
     761
     762                DO  i = nxl, nxr
     763                   DO  j = nys, nyn
     764                      DO  k = 0, leaf_area_density_f%nz - 1
     765                      IF ( leaf_area_density_f%var(k,j,i) /=                   &
     766                           leaf_area_density_f%fill )                          &
     767                         lad_s(k,j,i) = leaf_area_density_f%var(k,j,i)
     768                      ENDDO
     769                   ENDDO
     770                ENDDO
     771                CALL exchange_horiz( lad_s, nbgp )
     772!
     773!            ASCII file
    718774!--          Initialize canopy parameters cdc (canopy drag coefficient),
    719775!--          lsec (leaf scalar exchange coefficient), lsc (leaf surface concentration)
    720776!--          from file which contains complete 3D data (separate vertical profiles for
    721777!--          each location).
    722              CALL pcm_read_plant_canopy_3d
     778             ELSE
     779                CALL pcm_read_plant_canopy_3d
     780             ENDIF
    723781
    724782          CASE DEFAULT
     
    732790 
    733791       END SELECT
     792!
     793!--    Initialize 2D index array indicating canopy top index.
     794       ALLOCATE( pch_index_ji(nysg:nyng,nxlg:nxrg) )
     795       pch_index_ji = 0
     796
     797       DO  i = nxl, nxr
     798          DO  j = nys, nyn
     799             DO  k = 0, pch_index
     800                IF ( lad_s(k,j,i) /= 0 )  pch_index_ji(j,i) = k
     801             ENDDO
     802!
     803!--          Check whether topography and local vegetation on top exceed
     804!--          height of the model domain.
     805             k = get_topography_top_index( j, i, 's' )
     806             IF ( k + pch_index_ji(j,i) >= nzt + 1 )  THEN
     807                message_string =  'Local vegetation height on top of ' //      &
     808                                  'topography exceeds height of model domain.'
     809                CALL message( 'pcm_init', 'PA0999', 2, 2, 0, 6, 0 )
     810             ENDIF
     811
     812          ENDDO
     813       ENDDO
     814
     815       CALL exchange_horiz_2d_int( pch_index_ji, nys, nyn, nxl, nxr, nbgp )
    734816
    735817!
     
    757839          DO  i = nxlg, nxrg
    758840             DO  j = nysg, nyng
    759                 DO  k = pch_index-1, 0, -1
    760                    IF ( k == pch_index-1 )  THEN
     841                DO  k = pch_index_ji(j,i)-1, 0, -1
     842                   IF ( k == pch_index_ji(j,i)-1 )  THEN
    761843                      cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) +                &
    762844                         ( 0.5_wp * lad_s(k+1,j,i) *                           &
     
    815897          DO  i = nxlg, nxrg
    816898             DO  j = nysg, nyng
    817                 DO  k = 1, pch_index
     899                DO  k = 1, pch_index_ji(j,i)
    818900                   IF ( cum_lai_hf(0,j,i) /= 0.0_wp )  THEN
    819901                      pc_heating_rate(k,j,i) = cthf *                             &
     
    9451027                   
    9461028                CASE DEFAULT
    947                    write(message_string, '(a,i2,a)')   &
     1029                   WRITE(message_string, '(a,i2,a)')   &
    9481030                      'Unknown record type in file PLANT_CANOPY_DATA_3D: "', dtype, '"'
    9491031                   CALL message( 'pcm_read_plant_canopy_3d', 'PA0530', 1, 2, 0, 6, 0 )
     
    10291111!--                Determine topography-top index on u-grid
    10301112                   k_wall = get_topography_top_index( j, i, 'u' )
    1031                    DO  k = k_wall+1, k_wall+pch_index
     1113                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
    10321114
    10331115                      kk = k - k_wall   !- lad arrays are defined flat
     
    10941176                   k_wall = get_topography_top_index( j, i, 'v' )
    10951177
    1096                    DO  k = k_wall+1, k_wall+pch_index
     1178                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
    10971179
    10981180                      kk = k - k_wall   !- lad arrays are defined flat
     
    11591241                   k_wall = get_topography_top_index( j, i, 'w' )
    11601242
    1161                    DO  k = k_wall+1, k_wall+pch_index-1
     1243                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i) - 1
    11621244
    11631245                      kk = k - k_wall   !- lad arrays are defined flat
     
    12111293                   k_wall = get_topography_top_index( j, i, 's' )
    12121294
    1213                    DO  k = k_wall+1, k_wall+pch_index
     1295                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
    12141296
    12151297                      kk = k - k_wall   !- lad arrays are defined flat
     
    12281310                   k_wall = get_topography_top_index( j, i, 's' )
    12291311
    1230                    DO  k = k_wall+1, k_wall+pch_index
     1312                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
    12311313
    12321314                      kk = k - k_wall   !- lad arrays are defined flat
     
    12581340                   k_wall = get_topography_top_index( j, i, 's' )
    12591341
    1260                    DO  k = k_wall+1, k_wall+pch_index
     1342                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
    12611343
    12621344                      kk = k - k_wall   !- lad arrays are defined flat
     
    12871369                   k_wall = get_topography_top_index( j, i, 's' )
    12881370
    1289                    DO  k = k_wall+1, k_wall+pch_index
     1371                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
    12901372
    12911373                      kk = k - k_wall   !- lad arrays are defined flat
     
    13701452
    13711453       ddt_3d = 1.0_wp / dt_3d
    1372 
    13731454!
    13741455!--    Compute drag for the three velocity components and the SGS-TKE
     
    13811462!--          Determine topography-top index on u-grid
    13821463             k_wall = get_topography_top_index( j, i, 'u' )
    1383 
    1384              DO  k = k_wall+1, k_wall+pch_index
     1464             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
    13851465
    13861466                kk = k - k_wall  !- lad arrays are defined flat
     1467
    13871468!
    13881469!--             In order to create sharp boundaries of the plant canopy,
     
    14421523             k_wall = get_topography_top_index( j, i, 'v' )
    14431524
    1444              DO  k = k_wall+1, k_wall+pch_index
     1525             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
    14451526
    14461527                kk = k - k_wall  !- lad arrays are defined flat
     
    15021583             k_wall = get_topography_top_index( j, i, 'w' )
    15031584
    1504              DO  k = k_wall+1, k_wall+pch_index-1
     1585             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i) - 1
    15051586
    15061587                kk = k - k_wall  !- lad arrays are defined flat
     
    15491630             k_wall = get_topography_top_index( j, i, 's' )
    15501631
    1551              DO  k = k_wall+1, k_wall+pch_index
     1632             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
    15521633                kk = k - k_wall  !- lad arrays are defined flat
    15531634                tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i)
     
    15621643             k_wall = get_topography_top_index( j, i, 's' )
    15631644
    1564              DO  k = k_wall+1, k_wall+pch_index
     1645             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
    15651646
    15661647                kk = k - k_wall
     
    15881669             k_wall = get_topography_top_index( j, i, 's' )
    15891670
    1590              DO  k = k_wall+1, k_wall+pch_index
     1671             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
    15911672
    15921673                kk = k - k_wall
     
    16141695             k_wall = get_topography_top_index( j, i, 's' )
    16151696
    1616              DO  k = k_wall+1, k_wall+pch_index
     1697             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
    16171698
    16181699                kk = k - k_wall
Note: See TracChangeset for help on using the changeset viewer.