Index: palm/trunk/SOURCE/plant_canopy_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 2001)
+++ palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 2007)
@@ -20,5 +20,9 @@
 ! Current revisions:
 ! -----------------
-! 
+! Added SUBROUTINE pcm_read_plant_canopy_3d for reading 3d plant canopy data
+! from file (new case canopy_mode=read_from_file_3d) in the course of
+! introduction of urban surface model,
+! introduced variable ext_coef,
+! resorted SUBROUTINEs to alphabetical order
 ! 
 ! Former revisions:
@@ -120,4 +124,5 @@
     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.
@@ -154,5 +159,7 @@
 !
 !-- Public variables and constants
-    PUBLIC canopy_mode, cthf, dt_plant_canopy, lad, lad_s, plant_canopy
+    PUBLIC canopy_heat_flux, canopy_mode, cthf, dt_plant_canopy, lad, lad_s,   &
+           pch_index, plant_canopy
+           
 
 
@@ -171,5 +178,9 @@
     INTERFACE pcm_parin
        MODULE PROCEDURE pcm_parin
-    END INTERFACE pcm_parin  
+    END INTERFACE pcm_parin
+
+    INTERFACE pcm_read_plant_canopy_3d
+       MODULE PROCEDURE pcm_read_plant_canopy_3d
+    END INTERFACE pcm_read_plant_canopy_3d
     
     INTERFACE pcm_tendency
@@ -363,18 +374,23 @@
 
        USE control_parameters,                                                 &
-           ONLY: dz, ocean, passive_scalar
+           ONLY: coupling_char, dz, humidity, io_blocks, io_group,             &
+                 message_string, ocean, passive_scalar 
 
 
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  i !< running index
-       INTEGER(iwp) ::  j !< running index
-       INTEGER(iwp) ::  k !< running index
-
-       REAL(wp) ::  int_bpdf      !< vertical integral for lad-profile construction
-       REAL(wp) ::  dzh           !< vertical grid spacing in units of canopy height
-       REAL(wp) ::  gradient      !< gradient for lad-profile construction
-       REAL(wp) ::  canopy_height !< canopy height for lad-profile construction
-
+       CHARACTER(10) :: pct
+       
+       INTEGER(iwp) ::  i   !< running index
+       INTEGER(iwp) ::  ii  !< index       
+       INTEGER(iwp) ::  j   !< running index
+       INTEGER(iwp) ::  k   !< running index
+
+       REAL(wp) ::  int_bpdf        !< vertical integral for lad-profile construction
+       REAL(wp) ::  dzh             !< vertical grid spacing in units of canopy height
+       REAL(wp) ::  gradient        !< gradient for lad-profile construction
+       REAL(wp) ::  canopy_height   !< canopy height for lad-profile construction
+       REAL(wp) ::  pcv(nzb:nzt+1)  !<
+       
 !
 !--    Allocate one-dimensional arrays for the computation of the 
@@ -509,13 +525,20 @@
              ENDDO
 
+          CASE ( 'read_from_file_3d' )
+!
+!--          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
+
           CASE DEFAULT
-
-!
-!--       The DEFAULT case is reached either if the parameter 
-!--       canopy mode contains a wrong character string or if the 
-!--       user has coded a special case in the user interface. 
-!--       There, the subroutine user_init_plant_canopy checks 
-!--       which of these two conditions applies.
-          CALL user_init_plant_canopy
+!
+!--          The DEFAULT case is reached either if the parameter 
+!--          canopy mode contains a wrong character string or if the 
+!--          user has coded a special case in the user interface. 
+!--          There, the subroutine user_init_plant_canopy checks 
+!--          which of these two conditions applies.
+             CALL user_init_plant_canopy
  
        END SELECT
@@ -561,5 +584,5 @@
                 DO  k = 0, pch_index
                    canopy_heat_flux(k,j,i) = cthf *                            &
-                                             exp( -0.6_wp * cum_lai_hf(k,j,i) )
+                                             exp( -ext_coef * cum_lai_hf(k,j,i) )
                 ENDDO
              ENDDO
@@ -585,4 +608,135 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Parin for &canopy_par for plant canopy model
+!------------------------------------------------------------------------------!
+    SUBROUTINE pcm_parin
+
+
+       IMPLICIT NONE
+
+       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file 
+       
+       NAMELIST /canopy_par/      alpha_lad, beta_lad, canopy_drag_coeff,      &
+                                  canopy_mode, cthf,                           &
+                                  lad_surface,                                 & 
+                                  lad_vertical_gradient,                       &
+                                  lad_vertical_gradient_level,                 &
+                                  lai_beta,                                    &
+                                  leaf_scalar_exch_coeff,                      &
+                                  leaf_surface_conc, pch_index
+       
+       line = ' '
+       
+!
+!--    Try to find radiation model package
+       REWIND ( 11 )
+       line = ' '
+       DO   WHILE ( INDEX( line, '&canopy_par' ) == 0 )
+          READ ( 11, '(A)', END=10 )  line
+       ENDDO
+       BACKSPACE ( 11 )
+
+!
+!--    Read user-defined namelist
+       READ ( 11, canopy_par )
+
+!
+!--    Set flag that indicates that the radiation model is switched on
+       plant_canopy = .TRUE.
+
+ 10    CONTINUE
+       
+
+    END SUBROUTINE pcm_parin
+
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!
+!> Loads 3D plant canopy data from file. File format is as follows:
+!>
+!> num_levels
+!> dtype,x,y,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
+!> dtype,x,y,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
+!> dtype,x,y,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
+!> ...
+!>
+!> i.e. first line determines number of levels and further lines represent plant
+!> canopy data, one line per column and variable. In each data line,
+!> dtype represents variable to be set:
+!>
+!> dtype=1: leaf area density (lad_s)
+!> dtype=2: canopy drag coefficient (cdc)
+!> dtype=3: leaf scalar exchange coefficient (lsec)
+!> dtype=4: leaf surface concentration (lsc)
+!>
+!> Zeros are added automatically above num_levels until top of domain.  Any
+!> non-specified (x,y) columns have zero values as default.
+!------------------------------------------------------------------------------!
+    SUBROUTINE pcm_read_plant_canopy_3d
+        USE control_parameters, &
+            ONLY: passive_scalar, message_string
+        IMPLICIT NONE
+
+        INTEGER(iwp)                            :: i, j, dtype, nzp, nzpltop, nzpl, kk
+        REAL(wp), DIMENSION(:), ALLOCATABLE     :: col
+       
+        lad_s = 0.0_wp
+!         cdc = 0.0_wp
+!         if ( passive_scalar )  then
+!            lsc = 0.0_wp
+!            lsec = 0.0_wp
+!         endif
+
+        WRITE(9,*) 'Reading PLANT_CANOPY_DATA_3D', nzt
+        FLUSH(9)
+        OPEN(152, file='PLANT_CANOPY_DATA_3D', access='SEQUENTIAL', &
+                action='READ', status='OLD', form='FORMATTED', err=515)
+        READ(152, *, err=516, end=517) nzp   !< read first line = number of vertical layers
+        ALLOCATE(col(1:nzp))
+        nzpltop = MIN(nzt+1, nzb+nzp-1)
+        nzpl = nzpltop - nzb + 1    !< no. of layers to assign
+
+        DO
+            READ(152, *, err=516, end=517) dtype, i, j, col(:)
+            IF ( i < nxlg .or. i > nxrg .or. j < nysg .or. j > nyng ) CYCLE
+            WRITE(9,*) 'Read ', i,j,nzb_s_inner(j,i),col(:)
+            FLUSH(9)
+
+            SELECT CASE (dtype)
+              CASE( 1 ) !< leaf area density
+                !-- only lad_s has flat z-coordinate, others have regular
+                kk = nzb_s_inner(j, i)
+                lad_s(nzb:nzpltop-kk, j, i) = col(1+kk:nzpl)
+!               CASE( 2 ) !< canopy drag coefficient
+!                 cdc(nzb:nzpltop, j, i) = col(1:nzpl)
+!               CASE( 3 ) !< leaf scalar exchange coefficient
+!                 lsec(nzb:nzpltop, j, i) = col(1:nzpl)
+!               CASE( 4 ) !< leaf surface concentration
+!                 lsc(nzb:nzpltop, j, i) = col(1:nzpl)
+              CASE DEFAULT
+                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 )
+            END SELECT
+        ENDDO
+
+515     message_string = 'error opening file PLANT_CANOPY_DATA_3D'
+        CALL message( 'pcm_read_plant_canopy_3d', 'PA0531', 1, 2, 0, 6, 0 )
+
+516     message_string = 'error reading file PLANT_CANOPY_DATA_3D'
+        CALL message( 'pcm_read_plant_canopy_3d', 'PA0532', 1, 2, 0, 6, 0 )
+
+517     CLOSE(152)
+        DEALLOCATE(col)
+        
+    END SUBROUTINE pcm_read_plant_canopy_3d
+    
+    
 
 !------------------------------------------------------------------------------!
@@ -913,51 +1067,4 @@
 ! Description:
 ! ------------
-!> Parin for &canopy_par for plant canopy model
-!------------------------------------------------------------------------------!
-    SUBROUTINE pcm_parin
-
-
-       IMPLICIT NONE
-
-       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file 
-       
-        NAMELIST /canopy_par/     alpha_lad, beta_lad, canopy_drag_coeff,      &
-                                  canopy_mode, cthf,                           &
-                                  lad_surface,                                 & 
-                                  lad_vertical_gradient,                       &
-                                  lad_vertical_gradient_level,                 &
-                                  lai_beta,                                    &
-                                  leaf_scalar_exch_coeff,                      &
-                                  leaf_surface_conc, pch_index
-       
-       line = ' '
-       
-!
-!--    Try to find radiation model package
-       REWIND ( 11 )
-       line = ' '
-       DO   WHILE ( INDEX( line, '&canopy_par' ) == 0 )
-          READ ( 11, '(A)', END=10 )  line
-       ENDDO
-       BACKSPACE ( 11 )
-
-!
-!--    Read user-defined namelist
-       READ ( 11, canopy_par )
-
-!
-!--    Set flag that indicates that the radiation model is switched on
-       plant_canopy = .TRUE.
-
- 10    CONTINUE
-       
-
-    END SUBROUTINE pcm_parin    
-    
-
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
 !> Calculation of the tendency terms, accounting for the effect of the plant
 !> canopy on momentum and scalar quantities.
@@ -1249,3 +1356,5 @@
     END SUBROUTINE pcm_tendency_ij
 
+
+
  END MODULE plant_canopy_model_mod
