Index: palm/trunk/SOURCE/plant_canopy_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 2669)
+++ palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 2696)
@@ -1,5 +1,5 @@
 !> @file plant_canopy_model_mod.f90
 !------------------------------------------------------------------------------!
-! This file is part of PALM.
+! This file is part of the PALM model system.
 !
 ! PALM is free software: you can redistribute it and/or modify it under the
@@ -25,4 +25,12 @@
 ! -----------------
 ! $Id$
+! Bugfix for vertical loop index pch_index in case of Netcdf input
+! Introduce 2D index array incorporate canopy top index
+! Check if canopy on top of topography do not exceed vertical dimension 
+! Add check for canopy_mode in case of Netcdf input.
+! Enable _FillValue output for 3d quantities
+! Bugfix in reading of PIDS leaf area density (MS)
+! 
+! 2669 2017-12-06 16:03:27Z raasch
 ! coupling_char removed
 ! 
@@ -151,31 +159,30 @@
     CHARACTER (LEN=20)   ::  canopy_mode = 'block' !< canopy coverage
 
-    INTEGER(iwp) ::  pch_index = 0                 !< plant canopy height/top index
-    INTEGER(iwp) ::                                                            &
-       lad_vertical_gradient_level_ind(10) = -9999 !< lad-profile levels (index)
+    INTEGER(iwp) ::  pch_index = 0                               !< plant canopy height/top index
+    INTEGER(iwp) ::  lad_vertical_gradient_level_ind(10) = -9999 !< lad-profile levels (index)
+
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  pch_index_ji   !< local plant canopy top
 
     LOGICAL ::  calc_beta_lad_profile = .FALSE. !< switch for calc. of lad from beta func.
     LOGICAL ::  plant_canopy = .FALSE.          !< switch for use of canopy model
-
-    REAL(wp) ::  alpha_lad = 9999999.9_wp   !< coefficient for lad calculation
-    REAL(wp) ::  beta_lad = 9999999.9_wp    !< coefficient for lad calculation
-    REAL(wp) ::  canopy_drag_coeff = 0.0_wp !< canopy drag coefficient (parameter)
-    REAL(wp) ::  cdc = 0.0_wp               !< canopy drag coeff. (abbreviation used in equations)
-    REAL(wp) ::  cthf = 0.0_wp              !< canopy top heat flux
-    REAL(wp) ::  dt_plant_canopy = 0.0_wp   !< timestep account. for canopy drag
-    REAL(wp) ::  ext_coef = 0.6_wp          !< extinction coefficient
-    REAL(wp) ::  lad_surface = 0.0_wp       !< lad surface value
-    REAL(wp) ::  lai_beta = 0.0_wp          !< leaf area index (lai) for lad calc.
-    REAL(wp) ::                                                                &
-       leaf_scalar_exch_coeff = 0.0_wp      !< canopy scalar exchange coeff.
-    REAL(wp) ::                                                                &
-       leaf_surface_conc = 0.0_wp           !< leaf surface concentration
-    REAL(wp) ::  lsec = 0.0_wp              !< leaf scalar exchange coeff.
-    REAL(wp) ::  lsc = 0.0_wp               !< leaf surface concentration
-
-    REAL(wp) ::                                                                &
-       lad_vertical_gradient(10) = 0.0_wp              !< lad gradient
-    REAL(wp) ::                                                                &
-       lad_vertical_gradient_level(10) = -9999999.9_wp !< lad-prof. levels (in m)
+    LOGICAL ::  usm_lad_rma = .TRUE.            !< use MPI RMA to access LAD for raytracing (instead of global array)
+
+    REAL(wp) ::  alpha_lad = 9999999.9_wp         !< coefficient for lad calculation
+    REAL(wp) ::  beta_lad = 9999999.9_wp          !< coefficient for lad calculation
+    REAL(wp) ::  canopy_drag_coeff = 0.0_wp       !< canopy drag coefficient (parameter)
+    REAL(wp) ::  cdc = 0.0_wp                     !< canopy drag coeff. (abbreviation used in equations)
+    REAL(wp) ::  cthf = 0.0_wp                    !< canopy top heat flux
+    REAL(wp) ::  dt_plant_canopy = 0.0_wp         !< timestep account. for canopy drag
+    REAL(wp) ::  ext_coef = 0.6_wp                !< extinction coefficient
+    REAL(wp) ::  lad_surface = 0.0_wp             !< lad surface value
+    REAL(wp) ::  lai_beta = 0.0_wp                !< leaf area index (lai) for lad calc.
+    REAL(wp) ::  leaf_scalar_exch_coeff = 0.0_wp  !< canopy scalar exchange coeff.
+    REAL(wp) ::  leaf_surface_conc = 0.0_wp       !< leaf surface concentration
+    REAL(wp) ::  lsec = 0.0_wp                    !< leaf scalar exchange coeff.
+    REAL(wp) ::  lsc = 0.0_wp                     !< leaf surface concentration
+    REAL(wp) ::  prototype_lad                    !< prototype leaf area density for computing effective optical depth
+
+    REAL(wp) ::  lad_vertical_gradient(10) = 0.0_wp              !< lad gradient
+    REAL(wp) ::  lad_vertical_gradient_level(10) = -9999999.9_wp !< lad-prof. levels (in m)
 
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  lad            !< leaf area density
@@ -201,5 +208,5 @@
 !-- Public variables and constants
     PUBLIC pc_heating_rate, canopy_mode, cthf, dt_plant_canopy, lad, lad_s,   &
-           pch_index, plant_canopy
+           pch_index, plant_canopy, prototype_lad, usm_lad_rma
            
 
@@ -288,5 +295,9 @@
 
        USE control_parameters,                                                 &
-           ONLY: cloud_physics, message_string, microphysics_seifert
+           ONLY: cloud_physics, coupling_char, message_string,                 &
+                 microphysics_seifert
+
+       USE netcdf_data_input_mod,                                              &
+           ONLY:  input_file_static, input_pids_static
                  
     
@@ -331,4 +342,14 @@
           CALL message( 'check_parameters', 'PA0360', 1, 2, 0, 6, 0 )
        ENDIF
+!
+!--    If dynamic input file is used, canopy need to be read from file
+       IF ( input_pids_static  .AND.                                           &
+            TRIM( canopy_mode ) /= 'read_from_file_3d' )  THEN
+          message_string = 'Usage of dynamic input file ' //                   &
+                           TRIM( input_file_static ) //                        &
+                           TRIM( coupling_char ) // ' requires ' //            &
+                           'canopy_mode = read_from_file_3d'
+          CALL message( 'check_parameters', 'PA0999', 1, 2, 0, 6, 0 )
+       ENDIF
 
  
@@ -342,5 +363,5 @@
 !> Subroutine defining 3D output variables
 !------------------------------------------------------------------------------!
- SUBROUTINE pcm_data_output_3d( av, variable, found, local_pf )
+ SUBROUTINE pcm_data_output_3d( av, variable, found, local_pf, fill_value )
  
     USE control_parameters,                                                    &
@@ -356,11 +377,13 @@
     CHARACTER (LEN=*) ::  variable !< 
 
-    INTEGER(iwp) ::  av    !< 
-    INTEGER(iwp) ::  i     !< 
-    INTEGER(iwp) ::  j     !< 
-    INTEGER(iwp) ::  k     !< 
+    INTEGER(iwp) ::  av     !< 
+    INTEGER(iwp) ::  i      !< 
+    INTEGER(iwp) ::  j      !< 
+    INTEGER(iwp) ::  k      !< 
+    INTEGER(iwp) ::  k_topo !< topography top index
 
     LOGICAL      ::  found !< 
 
+    REAL(wp)     ::  fill_value
     REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb:nz_do3d) ::  local_pf !< 
 
@@ -368,4 +391,5 @@
     found = .TRUE.
 
+    local_pf = REAL( fill_value, KIND = 4 )
 
     SELECT CASE ( TRIM( variable ) )
@@ -375,7 +399,10 @@
             DO  i = nxl, nxr
                DO  j = nys, nyn
-                  DO  k = nzb_s_inner(j,i), nz_do3d
-                     local_pf(i,j,k) = pc_heating_rate(k-nzb_s_inner(j,i),j,i)
-                  ENDDO
+                  IF ( pch_index_ji(j,i) /= 0 )  THEN
+                     k_topo = get_topography_top_index( j, i, 's' )
+                     DO  k = k_topo, k_topo + pch_index_ji(j,i)
+                        local_pf(i,j,k) = pc_heating_rate(k-k_topo,j,i)
+                     ENDDO
+                  ENDIF
                ENDDO
             ENDDO
@@ -387,7 +414,10 @@
             DO  i = nxl, nxr
                DO  j = nys, nyn
-                  DO  k = nzb_s_inner(j,i), nz_do3d
-                     local_pf(i,j,k) = lad_s(k-nzb_s_inner(j,i),j,i)
-                  ENDDO
+                  IF ( pch_index_ji(j,i) /= 0 )  THEN
+                     k_topo = get_topography_top_index( j, i, 's' )
+                     DO  k = k_topo, k_topo + pch_index_ji(j,i)
+                        local_pf(i,j,k) = lad_s(k-k_topo,j,i)
+                     ENDDO
+                  ENDIF
                ENDDO
             ENDDO
@@ -570,4 +600,7 @@
                  passive_scalar, urban_surface
 
+       USE netcdf_data_input_mod,                                              &
+           ONLY:  leaf_area_density_f
+
        USE surface_mod,                                                        &
            ONLY: surf_def_h, surf_lsm_h, surf_usm_h
@@ -716,9 +749,34 @@
           CASE ( 'read_from_file_3d' )
 !
+!--          Initialize LAD with data from file. If LAD is given in NetCDF file,
+!--          use these values, else take LAD profiles from ASCII file. 
+!--          Please note, in NetCDF file LAD is only given up to the maximum 
+!--          canopy top, indicated by leaf_area_density_f%nz.  
+             lad_s = 0.0_wp
+             IF ( leaf_area_density_f%from_file )  THEN
+!
+!--             Set also pch_index, used to be the upper bound of the vertical 
+!--             loops. Therefore, use the global top of the canopy layer. 
+                pch_index = leaf_area_density_f%nz - 1
+
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = 0, leaf_area_density_f%nz - 1
+                      IF ( leaf_area_density_f%var(k,j,i) /=                   &
+                           leaf_area_density_f%fill )                          &
+                         lad_s(k,j,i) = leaf_area_density_f%var(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+                CALL exchange_horiz( lad_s, nbgp )
+!
+!            ASCII file
 !--          Initialize canopy parameters cdc (canopy drag coefficient),
 !--          lsec (leaf scalar exchange coefficient), lsc (leaf surface concentration)
 !--          from file which contains complete 3D data (separate vertical profiles for
 !--          each location).
-             CALL pcm_read_plant_canopy_3d
+             ELSE
+                CALL pcm_read_plant_canopy_3d
+             ENDIF
 
           CASE DEFAULT
@@ -732,4 +790,28 @@
  
        END SELECT
+!
+!--    Initialize 2D index array indicating canopy top index. 
+       ALLOCATE( pch_index_ji(nysg:nyng,nxlg:nxrg) )
+       pch_index_ji = 0
+
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = 0, pch_index
+                IF ( lad_s(k,j,i) /= 0 )  pch_index_ji(j,i) = k
+             ENDDO
+!
+!--          Check whether topography and local vegetation on top exceed 
+!--          height of the model domain.
+             k = get_topography_top_index( j, i, 's' )
+             IF ( k + pch_index_ji(j,i) >= nzt + 1 )  THEN
+                message_string =  'Local vegetation height on top of ' //      &
+                                  'topography exceeds height of model domain.'
+                CALL message( 'pcm_init', 'PA0999', 2, 2, 0, 6, 0 )
+             ENDIF
+
+          ENDDO
+       ENDDO
+
+       CALL exchange_horiz_2d_int( pch_index_ji, nys, nyn, nxl, nxr, nbgp )
 
 !
@@ -757,6 +839,6 @@
           DO  i = nxlg, nxrg
              DO  j = nysg, nyng
-                DO  k = pch_index-1, 0, -1
-                   IF ( k == pch_index-1 )  THEN
+                DO  k = pch_index_ji(j,i)-1, 0, -1
+                   IF ( k == pch_index_ji(j,i)-1 )  THEN
                       cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) +                &
                          ( 0.5_wp * lad_s(k+1,j,i) *                           &
@@ -815,5 +897,5 @@
           DO  i = nxlg, nxrg
              DO  j = nysg, nyng
-                DO  k = 1, pch_index
+                DO  k = 1, pch_index_ji(j,i)
                    IF ( cum_lai_hf(0,j,i) /= 0.0_wp )  THEN
                       pc_heating_rate(k,j,i) = cthf *                             &
@@ -945,5 +1027,5 @@
                     
                 CASE DEFAULT
-                   write(message_string, '(a,i2,a)')   &
+                   WRITE(message_string, '(a,i2,a)')   &
                       'Unknown record type in file PLANT_CANOPY_DATA_3D: "', dtype, '"'
                    CALL message( 'pcm_read_plant_canopy_3d', 'PA0530', 1, 2, 0, 6, 0 )
@@ -1029,5 +1111,5 @@
 !--                Determine topography-top index on u-grid
                    k_wall = get_topography_top_index( j, i, 'u' )
-                   DO  k = k_wall+1, k_wall+pch_index
+                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
 
                       kk = k - k_wall   !- lad arrays are defined flat
@@ -1094,5 +1176,5 @@
                    k_wall = get_topography_top_index( j, i, 'v' )
 
-                   DO  k = k_wall+1, k_wall+pch_index
+                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
 
                       kk = k - k_wall   !- lad arrays are defined flat
@@ -1159,5 +1241,5 @@
                    k_wall = get_topography_top_index( j, i, 'w' )
 
-                   DO  k = k_wall+1, k_wall+pch_index-1
+                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i) - 1
 
                       kk = k - k_wall   !- lad arrays are defined flat
@@ -1211,5 +1293,5 @@
                    k_wall = get_topography_top_index( j, i, 's' )
 
-                   DO  k = k_wall+1, k_wall+pch_index
+                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
 
                       kk = k - k_wall   !- lad arrays are defined flat
@@ -1228,5 +1310,5 @@
                    k_wall = get_topography_top_index( j, i, 's' )
 
-                   DO  k = k_wall+1, k_wall+pch_index
+                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
 
                       kk = k - k_wall   !- lad arrays are defined flat
@@ -1258,5 +1340,5 @@
                    k_wall = get_topography_top_index( j, i, 's' )
 
-                   DO  k = k_wall+1, k_wall+pch_index
+                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
 
                       kk = k - k_wall   !- lad arrays are defined flat
@@ -1287,5 +1369,5 @@
                    k_wall = get_topography_top_index( j, i, 's' )
 
-                   DO  k = k_wall+1, k_wall+pch_index
+                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
 
                       kk = k - k_wall   !- lad arrays are defined flat
@@ -1370,5 +1452,4 @@
 
        ddt_3d = 1.0_wp / dt_3d
-
 !
 !--    Compute drag for the three velocity components and the SGS-TKE
@@ -1381,8 +1462,8 @@
 !--          Determine topography-top index on u-grid
              k_wall = get_topography_top_index( j, i, 'u' )
-
-             DO  k = k_wall+1, k_wall+pch_index
+             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
 
                 kk = k - k_wall  !- lad arrays are defined flat
+
 !
 !--             In order to create sharp boundaries of the plant canopy, 
@@ -1442,5 +1523,5 @@
              k_wall = get_topography_top_index( j, i, 'v' )
 
-             DO  k = k_wall+1, k_wall+pch_index
+             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
 
                 kk = k - k_wall  !- lad arrays are defined flat
@@ -1502,5 +1583,5 @@
              k_wall = get_topography_top_index( j, i, 'w' )
 
-             DO  k = k_wall+1, k_wall+pch_index-1
+             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i) - 1
 
                 kk = k - k_wall  !- lad arrays are defined flat
@@ -1549,5 +1630,5 @@
              k_wall = get_topography_top_index( j, i, 's' )
 
-             DO  k = k_wall+1, k_wall+pch_index
+             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
                 kk = k - k_wall  !- lad arrays are defined flat
                 tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i)
@@ -1562,5 +1643,5 @@
              k_wall = get_topography_top_index( j, i, 's' )
 
-             DO  k = k_wall+1, k_wall+pch_index
+             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
 
                 kk = k - k_wall
@@ -1588,5 +1669,5 @@
              k_wall = get_topography_top_index( j, i, 's' )
 
-             DO  k = k_wall+1, k_wall+pch_index
+             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
 
                 kk = k - k_wall
@@ -1614,5 +1695,5 @@
              k_wall = get_topography_top_index( j, i, 's' )
 
-             DO  k = k_wall+1, k_wall+pch_index
+             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
 
                 kk = k - k_wall
