Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 4131)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 4150)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Input of additional surface variables independent on land- or urban-surface
+! model
+! 
+! 4131 2019-08-02 11:06:18Z monakurppa
 ! Allocate sums and sums_l to allow profile output for salsa variables.
 ! 
@@ -661,7 +665,18 @@
 
     USE netcdf_data_input_mod,                                                 &
-        ONLY:  init_3d,                              &
-               netcdf_data_input_init_3d, netcdf_data_input_interpolate
-
+        ONLY:  char_fill,                                                      &
+               check_existence,                                                &
+               close_input_file,                                               &
+               get_attribute,                                                  &
+               get_variable,                                                   &
+               init_3d,                                                        &
+               input_pids_static,                                              &
+               inquire_num_variables,                                          &
+               inquire_variable_names,                                         &
+               input_file_static,                                              &
+               netcdf_data_input_init_3d,                                      &
+               open_read_file,                                                 &
+               real_2d
+               
     USE nesting_offl_mod,                                                      &
         ONLY:  nesting_offl_init
@@ -695,6 +710,12 @@
 
     USE surface_mod,                                                           &
-        ONLY :  init_surface_arrays, init_surfaces, surf_def_h, surf_lsm_h,    &
-                surf_usm_h, get_topography_top_index_ji
+        ONLY :  init_single_surface_properties,                                &
+                init_surface_arrays,                                           &
+                init_surfaces,                                                 &
+                surf_def_h,                                                    &
+                surf_def_v,                                                    &
+                surf_lsm_h,                                                    &
+                surf_usm_h,                                                    &
+                get_topography_top_index_ji
 
 #if defined( _OPENACC )
@@ -710,11 +731,25 @@
     IMPLICIT NONE
 
-    INTEGER(iwp) ::  i             !< grid index in x direction
-    INTEGER(iwp) ::  ind_array(1)  !< dummy used to determine start index for external pressure forcing
-    INTEGER(iwp) ::  j             !< grid index in y direction
-    INTEGER(iwp) ::  k             !< grid index in z direction
-    INTEGER(iwp) ::  k_surf        !< surface level index
-    INTEGER(iwp) ::  m             !< index of surface element in surface data type
-    INTEGER(iwp) ::  sr            !< index of statistic region
+    CHARACTER(LEN=50), DIMENSION(:), ALLOCATABLE ::  vars_pids_static  !< variable names in static input file
+    
+    INTEGER(iwp) ::  i                    !< grid index in x direction
+    INTEGER(iwp) ::  ind_array(1)         !< dummy used to determine start index for external pressure forcing
+    INTEGER(iwp) ::  j                    !< grid index in y direction
+    INTEGER(iwp) ::  k                    !< grid index in z direction
+    INTEGER(iwp) ::  k_surf               !< surface level index
+    INTEGER(iwp) ::  l                    !< running index over surface orientation   
+    INTEGER(iwp) ::  m                    !< index of surface element in surface data type    
+    INTEGER(iwp) ::  num_vars_pids_static !< number of variables in static input file
+    INTEGER(iwp) ::  nz_u_shift           !< topography-top index on u-grid, used to vertically shift initial profiles
+    INTEGER(iwp) ::  nz_v_shift           !< topography-top index on v-grid, used to vertically shift initial profiles
+    INTEGER(iwp) ::  nz_w_shift           !< topography-top index on w-grid, used to vertically shift initial profiles
+    INTEGER(iwp) ::  nz_s_shift           !< topography-top index on scalar-grid, used to vertically shift initial profiles
+    INTEGER(iwp) ::  nz_u_shift_l         !< topography-top index on u-grid, used to vertically shift initial profiles
+    INTEGER(iwp) ::  nz_v_shift_l         !< topography-top index on v-grid, used to vertically shift initial profiles
+    INTEGER(iwp) ::  nz_w_shift_l         !< topography-top index on w-grid, used to vertically shift initial profiles
+    INTEGER(iwp) ::  nz_s_shift_l         !< topography-top index on scalar-grid, used to vertically shift initial profiles
+    INTEGER(iwp) ::  nzt_l                !< index of top PE boundary for multigrid level
+    INTEGER(iwp) ::  pids_static_id       !< file id for static input file
+    INTEGER(iwp) ::  sr                   !< index of statistic region
 
     INTEGER(iwp), DIMENSION(:), ALLOCATABLE   ::  ngp_2dh_l  !< toal number of horizontal grid points in statistical region on subdomain
@@ -723,9 +758,9 @@
     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  ngp_2dh_s_inner_l  !< number of horizontal non-topography grid points on subdomain
 
+
+    
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  init_l        !< dummy array used for averaging 3D data to obtain inital profiles 
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  p_hydrostatic !< hydrostatic pressure
 
-    INTEGER(iwp) ::  l       !< loop variable
-    INTEGER(iwp) ::  nzt_l   !< index of top PE boundary for multigrid level
     REAL(wp) ::  dx_l !< grid spacing along x on different multigrid level
     REAL(wp) ::  dy_l !< grid spacing along y on different multigrid level
@@ -738,14 +773,6 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  ngp_3d_inner_tmp  !< total number of non-topography grid points
 
-    INTEGER(iwp) ::  nz_u_shift   !< topography-top index on u-grid, used to vertically shift initial profiles
-    INTEGER(iwp) ::  nz_v_shift   !< topography-top index on v-grid, used to vertically shift initial profiles
-    INTEGER(iwp) ::  nz_w_shift   !< topography-top index on w-grid, used to vertically shift initial profiles
-    INTEGER(iwp) ::  nz_s_shift   !< topography-top index on scalar-grid, used to vertically shift initial profiles
-    INTEGER(iwp) ::  nz_u_shift_l !< topography-top index on u-grid, used to vertically shift initial profiles
-    INTEGER(iwp) ::  nz_v_shift_l !< topography-top index on v-grid, used to vertically shift initial profiles
-    INTEGER(iwp) ::  nz_w_shift_l !< topography-top index on w-grid, used to vertically shift initial profiles
-    INTEGER(iwp) ::  nz_s_shift_l !< topography-top index on scalar-grid, used to vertically shift initial profiles
-
-
+    TYPE(real_2d) ::  tmp_2d !< temporary variable to input additional surface-data from static file
+    
     CALL location_message( 'model initialization', 'start' )
 
@@ -1619,5 +1646,4 @@
           ENDDO
        ENDIF
-
 !
 !--    Initialization of the turbulence recycling method
@@ -2005,4 +2031,73 @@
        ENDIF
 
+    ENDIF
+!
+!-- In the following, surface properties can be further initialized with
+!-- input from static driver file.
+!-- At the moment this affects only default surfaces. For example, 
+!-- roughness length or sensible / latent heat fluxes can be initialized
+!-- heterogeneously for default surfaces. Therefore, a generic routine
+!-- from netcdf_data_input_mod is called to read a 2D array.
+    IF ( input_pids_static )  THEN
+!
+!--    Open the static input file
+       CALL open_read_file( TRIM( input_file_static ) //                    &
+                            TRIM( coupling_char ),                          &
+                            pids_static_id )
+                            
+       CALL inquire_num_variables( pids_static_id, num_vars_pids_static )
+!
+!--    Allocate memory to store variable names and read them
+       ALLOCATE( vars_pids_static(1:num_vars_pids_static) )
+       CALL inquire_variable_names( pids_static_id, vars_pids_static )
+                            
+       ALLOCATE( tmp_2d%var(nys:nyn,nxl:nxr) )
+       tmp_2d%var = 0.0_wp
+!
+!--    Input roughness length. 
+       IF ( check_existence( vars_pids_static, 'z0' ) )  THEN
+!
+!--       Read _FillValue attribute
+          CALL get_attribute( pids_static_id, char_fill, tmp_2d%fill,          &
+                              .FALSE., 'z0' )                                  
+!                                                                              
+!--       Read variable                                                        
+          CALL get_variable( pids_static_id, 'z0', tmp_2d%var,                 &
+                             nxl, nxr, nys, nyn )                              
+!                                                                              
+!--       Initialize roughness length. Note, z0 will be only initialized at    
+!--       default-type surfaces. At natural or urban z0 is implicitly          
+!--       initialized bythe respective parameter lists.                        
+!--       Initialize horizontal surface elements.                              
+          CALL init_single_surface_properties( surf_def_h(0)%z0,               &
+                                               tmp_2d%var,                     &
+                                               surf_def_h(0)%ns,               &
+                                               tmp_2d%fill,                    &
+                                               surf_def_h(0)%i,                &
+                                               surf_def_h(0)%j )               
+!                                                                              
+!--       Initialize roughness also at vertical surface elements.              
+!--       Note, the actual 2D input arrays are only defined on the             
+!--       subdomain. Therefore, pass the index arrays with their respective    
+!--       offset values.                                                       
+          DO  l = 0, 3                                                         
+             CALL init_single_surface_properties(                              &
+                                         surf_def_v(l)%z0,                     &
+                                         tmp_2d%var,                           &
+                                         surf_def_v(l)%ns,                     &
+                                         tmp_2d%fill,                          &
+                                         surf_def_v(l)%i + surf_def_v(l)%ioff, &
+                                         surf_def_v(l)%j + surf_def_v(l)%joff )
+          ENDDO
+          
+       ENDIF
+!
+!--    Additional variables, e.g. shf, qsws, etc, can be initialized the 
+!--    same way.
+                           
+       DEALLOCATE( tmp_2d%var )
+!
+!--    Finally, close the input file.
+       CALL close_input_file( pids_static_id )
     ENDIF
 !
