Index: palm/trunk/SOURCE/chem_emissions_mod.f90
===================================================================
--- palm/trunk/SOURCE/chem_emissions_mod.f90	(revision 4227)
+++ palm/trunk/SOURCE/chem_emissions_mod.f90	(revision 4230)
@@ -27,4 +27,8 @@
 ! -----------------
 ! $Id$
+! Bugfix, consider that time_since_reference_point can be also negative when 
+! time indices are determined.
+! 
+! 4227 2019-09-10 18:04:34Z gronemeier
 ! implement new palm_date_time_mod
 ! 
@@ -1091,7 +1095,9 @@
 
           CALL get_date_time( 0.0_wp, second_of_day=time_utc_init )
-          CALL get_date_time( time_since_reference_point, hour=hour_of_day )
-
-          days_since_reference_point = INT( FLOOR( ( time_utc_init + time_since_reference_point ) &
+          CALL get_date_time( MAX( 0.0_wp, time_since_reference_point ),       &
+                              hour=hour_of_day )
+
+          days_since_reference_point = INT( FLOOR( (                            &
+                    time_utc_init + MAX( 0.0_wp, time_since_reference_point ) ) &
                                                    / seconds_per_day ) )
 
@@ -1124,5 +1130,5 @@
 !-- Update time indices
 
-             CALL get_date_time( time_since_reference_point, &
+             CALL get_date_time( MAX( time_since_reference_point, 0.0_wp ),    &
                                  day_of_year=day_of_year, hour=hour_of_day )
              index_hh = ( day_of_year - 1_iwp ) * hour_of_day
@@ -1156,8 +1162,8 @@
 !
 !-- Update time indices
-             CALL get_date_time( time_since_reference_point, &
-                                 month=month_of_year,        &
-                                 day=day_of_month,           &
-                                 hour=hour_of_day,           &
+             CALL get_date_time( MAX( time_since_reference_point, 0.0_wp ), &
+                                 month=month_of_year,                       &
+                                 day=day_of_month,                          &
+                                 hour=hour_of_day,                          &
                                  day_of_week=day_of_week     )
              index_mm = month_of_year
@@ -1230,5 +1236,6 @@
 !
 !--       Get time-factor for specific hour
-          CALL get_date_time( time_since_reference_point, hour=hour_of_day )
+          CALL get_date_time( MAX( time_since_reference_point, 0.0_wp ),              &
+                              hour=hour_of_day )
 
           index_hh = hour_of_day
Index: palm/trunk/SOURCE/chemistry_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 4227)
+++ palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 4230)
@@ -27,4 +27,8 @@
 ! -----------------
 ! $Id$
+! Bugfix, initialize mean profiles also in restart runs. Also initialize 
+! array used for Runge-Kutta tendecies in restart runs.  
+! 
+! 4227 2019-09-10 18:04:34Z gronemeier
 ! implement new palm_date_time_mod
 ! 
@@ -1948,12 +1952,12 @@
           ENDDO
        ENDIF 
-!
-!--    Initiale old and new time levels. 
-       DO  lsp = 1, nvar
-          chem_species(lsp)%tconc_m = 0.0_wp                      
-          chem_species(lsp)%conc_p  = chem_species(lsp)%conc     
-       ENDDO
 
     ENDIF
+!
+!-- Initiale old and new time levels. Note, this has to be done also in restart runs
+    DO  lsp = 1, nvar
+       chem_species(lsp)%tconc_m = 0.0_wp                      
+       chem_species(lsp)%conc_p  = chem_species(lsp)%conc     
+    ENDDO
 
     DO  lsp = 1, nphot
@@ -1964,5 +1968,5 @@
 !--          WRITE(6,'(a,i4,3x,a)')  'Photolysis: ',lsp,TRIM( phot_names(lsp) )
 !--       ENDIF 
-          phot_frequen(lsp)%freq(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  =>  freq_1(:,:,:,lsp)
+       phot_frequen(lsp)%freq(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  =>  freq_1(:,:,:,lsp)
     ENDDO
 
@@ -1982,7 +1986,4 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE chem_init_profiles              
-!
-!-- SUBROUTINE is called from chem_init in case of TRIM( initializing_actions ) /= 'read_restart_data'
-!< We still need to see what has to be done in case of restart run
 
     USE chem_modules
@@ -2000,57 +2001,59 @@
 !-- "cs_heights" are prescribed, their values will!override the constant profile given by
 !-- "cs_surface".
-    IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN 
-       lsp_usr = 1
-       DO  WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )   !'novalue' is the default
-          DO  lsp = 1, nspec                                !
-!
-!--          create initial profile (conc_pr_init) for each chemical species
-             IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_usr) ) )  THEN   !
-                IF ( cs_profile(lsp_usr,1) == 9999999.9_wp )  THEN
-!
-!--               set a vertically constant profile based on the surface conc (cs_surface(lsp_usr)) of each species
-                   DO lpr_lev = 0, nzt+1
-                      chem_species(lsp)%conc_pr_init(lpr_lev) = cs_surface(lsp_usr)
-                   ENDDO
-                ELSE
-                   IF ( cs_heights(1,1) /= 0.0_wp )  THEN
-                      message_string = 'The surface value of cs_heights must be 0.0'
-                      CALL message( 'chem_check_parameters', 'CM0434', 1, 2, 0, 6, 0 )
+!     IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN 
+    lsp_usr = 1
+    DO  WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )   !'novalue' is the default
+       DO  lsp = 1, nspec                                !
+!
+!--       create initial profile (conc_pr_init) for each chemical species
+          IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_usr) ) )  THEN   !
+             IF ( cs_profile(lsp_usr,1) == 9999999.9_wp )  THEN
+!
+!--            set a vertically constant profile based on the surface conc (cs_surface(lsp_usr)) of each species
+                DO lpr_lev = 0, nzt+1
+                   chem_species(lsp)%conc_pr_init(lpr_lev) = cs_surface(lsp_usr)
+                ENDDO
+             ELSE
+                IF ( cs_heights(1,1) /= 0.0_wp )  THEN
+                   message_string = 'The surface value of cs_heights must be 0.0'
+                   CALL message( 'chem_check_parameters', 'CM0434', 1, 2, 0, 6, 0 )
+                ENDIF
+
+                use_prescribed_profile_data = .TRUE.
+
+                npr_lev = 1
+!                chem_species(lsp)%conc_pr_init(0) = 0.0_wp
+                DO  lpr_lev = 1, nz+1
+                   IF ( npr_lev < 100 )  THEN
+                      DO  WHILE ( cs_heights(lsp_usr, npr_lev+1) <= zu(lpr_lev) )
+                         npr_lev = npr_lev + 1
+                         IF ( npr_lev == 100 )  THEN
+                            message_string = 'number of chem spcs exceeding the limit'
+                            CALL message( 'chem_check_parameters', 'CM0435', 1, 2, 0, 6, 0 )               
+                            EXIT
+                         ENDIF
+                      ENDDO
                    ENDIF
-
-                   use_prescribed_profile_data = .TRUE.
-
-                   npr_lev = 1
-!                   chem_species(lsp)%conc_pr_init(0) = 0.0_wp
-                   DO  lpr_lev = 1, nz+1
-                      IF ( npr_lev < 100 )  THEN
-                         DO  WHILE ( cs_heights(lsp_usr, npr_lev+1) <= zu(lpr_lev) )
-                            npr_lev = npr_lev + 1
-                            IF ( npr_lev == 100 )  THEN
-                               message_string = 'number of chem spcs exceeding the limit'
-                               CALL message( 'chem_check_parameters', 'CM0435', 1, 2, 0, 6, 0 )               
-                               EXIT
-                            ENDIF
-                         ENDDO
-                      ENDIF
-                      IF ( npr_lev < 100  .AND.  cs_heights(lsp_usr,npr_lev+1) /= 9999999.9_wp )  THEN
-                         chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev) +       &
-                              ( zu(lpr_lev) - cs_heights(lsp_usr, npr_lev) ) /                          &
-                              ( cs_heights(lsp_usr, (npr_lev + 1)) - cs_heights(lsp_usr, npr_lev ) ) *  &
-                              ( cs_profile(lsp_usr, (npr_lev + 1)) - cs_profile(lsp_usr, npr_lev ) )
-                      ELSE
-                         chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev)
-                      ENDIF
-                   ENDDO
-                ENDIF
-!
-!--          If a profile is prescribed explicity using cs_profiles and cs_heights, then  
-!--          chem_species(lsp)%conc_pr_init is populated with the specific "lsp" based
-!--          on the cs_profiles(lsp_usr,:)  and cs_heights(lsp_usr,:). 
+                   IF ( npr_lev < 100  .AND.  cs_heights(lsp_usr,npr_lev+1) /= 9999999.9_wp )  THEN
+                      chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev) +       &
+                           ( zu(lpr_lev) - cs_heights(lsp_usr, npr_lev) ) /                          &
+                           ( cs_heights(lsp_usr, (npr_lev + 1)) - cs_heights(lsp_usr, npr_lev ) ) *  &
+                           ( cs_profile(lsp_usr, (npr_lev + 1)) - cs_profile(lsp_usr, npr_lev ) )
+                   ELSE
+                      chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev)
+                   ENDIF
+                ENDDO
              ENDIF
-          ENDDO
-          lsp_usr = lsp_usr + 1
+!
+!--       If a profile is prescribed explicity using cs_profiles and cs_heights, then  
+!--       chem_species(lsp)%conc_pr_init is populated with the specific "lsp" based
+!--       on the cs_profiles(lsp_usr,:)  and cs_heights(lsp_usr,:). 
+          ENDIF
+
        ENDDO
-    ENDIF
+
+       lsp_usr = lsp_usr + 1
+    ENDDO
+!     ENDIF
 
  END SUBROUTINE chem_init_profiles
Index: palm/trunk/SOURCE/nesting_offl_mod.f90
===================================================================
--- palm/trunk/SOURCE/nesting_offl_mod.f90	(revision 4227)
+++ palm/trunk/SOURCE/nesting_offl_mod.f90	(revision 4230)
@@ -20,9 +20,14 @@
 ! Current revisions:
 ! ------------------
-! implement new palm_date_time_mod
+! 
 ! 
 ! Former revisions:
 ! -----------------
 ! $Id$
+! Update mean chemistry profiles. These are also used for rayleigh damping. 
+! 
+! 4227 2019-09-10 18:04:34Z gronemeier
+! implement new palm_date_time_mod
+!
 ! - Data input moved into nesting_offl_mod
 ! - check rephrased
@@ -403,5 +408,5 @@
 !--    Check if dynamic driver data input is required.
        IF ( nest_offl%time(nest_offl%tind_p) <=                                &
-            MAX( time_since_reference_point, 0.0_wp) + time_utc_init  .OR.                   &
+            MAX( time_since_reference_point, 0.0_wp) + time_utc_init  .OR.     &
             .NOT.  nest_offl%init )  THEN
           CONTINUE
@@ -922,11 +927,13 @@
        REAL(wp), DIMENSION(nzb:nzt+1) ::  pt_ref   !< reference profile for potential temperature
        REAL(wp), DIMENSION(nzb:nzt+1) ::  pt_ref_l !< reference profile for potential temperature on subdomain
-       REAL(wp), DIMENSION(nzb:nzt+1) ::  q_ref    !< reference profile for mixing ratio on subdomain
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  q_ref    !< reference profile for mixing ratio
        REAL(wp), DIMENSION(nzb:nzt+1) ::  q_ref_l  !< reference profile for mixing ratio on subdomain
-       REAL(wp), DIMENSION(nzb:nzt+1) ::  u_ref    !< reference profile for u-component on subdomain
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  u_ref    !< reference profile for u-component
        REAL(wp), DIMENSION(nzb:nzt+1) ::  u_ref_l  !< reference profile for u-component on subdomain
-       REAL(wp), DIMENSION(nzb:nzt+1) ::  v_ref    !< reference profile for v-component on subdomain
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  v_ref    !< reference profile for v-component
        REAL(wp), DIMENSION(nzb:nzt+1) ::  v_ref_l  !< reference profile for v-component on subdomain
-       
+
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ref_chem   !< reference profile for chemical species
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ref_chem_l !< reference profile for chemical species on subdomain
 
        IF ( debug_output_timestep )  CALL debug_message( 'nesting_offl_bc', 'start' )
@@ -934,5 +941,5 @@
        CALL  cpu_log( log_point(58), 'offline nesting', 'start' )      
 !
-!--    Set mean profiles, derived from boundary data, to zero
+!--    Initialize mean profiles, derived from boundary data, to zero
        pt_ref   = 0.0_wp
        q_ref    = 0.0_wp
@@ -944,4 +951,12 @@
        u_ref_l  = 0.0_wp
        v_ref_l  = 0.0_wp
+!
+!--    If required, allocate temporary arrays to compute chemistry mean profiles
+       IF ( air_chemistry )  THEN
+          ALLOCATE( ref_chem(nzb:nzt+1,1:UBOUND( chem_species, 1 ) )   )
+          ALLOCATE( ref_chem_l(nzb:nzt+1,1:UBOUND( chem_species, 1 ) ) )
+          ref_chem   = 0.0_wp
+          ref_chem_l = 0.0_wp
+       ENDIF
 !
 !--    Set boundary conditions of u-, v-, w-component, as well as q, and pt.
@@ -1024,4 +1039,6 @@
                                                   fac_dt                   )
                       ENDDO
+                      ref_chem_l(nzb+1:nzt,n) = ref_chem_l(nzb+1:nzt,n)        &
+                                         + chem_species(n)%conc(nzb+1:nzt,j,-1)
                    ENDDO
                 ENDIF
@@ -1100,4 +1117,6 @@
                                                  fac_dt                       )
                       ENDDO
+                      ref_chem_l(nzb+1:nzt,n) = ref_chem_l(nzb+1:nzt,n)        &
+                                       + chem_species(n)%conc(nzb+1:nzt,j,nxr+1)
                    ENDDO
                 ENDIF
@@ -1179,4 +1198,6 @@
                                                  fac_dt                    )
                       ENDDO
+                      ref_chem_l(nzb+1:nzt,n) = ref_chem_l(nzb+1:nzt,n)        &
+                                       + chem_species(n)%conc(nzb+1:nzt,-1,i)
                    ENDDO
                 ENDIF
@@ -1257,4 +1278,6 @@
                                                  fac_dt                       )
                       ENDDO
+                      ref_chem_l(nzb+1:nzt,n) = ref_chem_l(nzb+1:nzt,n)        &
+                                       + chem_species(n)%conc(nzb+1:nzt,nyn+1,i)
                    ENDDO
                 ENDIF
@@ -1337,7 +1360,9 @@
                    DO  j = nys, nyn
                       chem_species(n)%conc(nzt+1,j,i) = interpolate_in_time(   &
-                                              nest_offl%chem_top(0,j,i,n),   &
-                                              nest_offl%chem_top(1,j,i,n),   &
+                                              nest_offl%chem_top(0,j,i,n),     &
+                                              nest_offl%chem_top(1,j,i,n),     &
                                               fac_dt                       )
+                      ref_chem_l(nzt+1,n) = ref_chem_l(nzt+1,n) +              &
+                                            chem_species(n)%conc(nzt+1,j,i)
                    ENDDO
                 ENDDO
@@ -1404,9 +1429,15 @@
                               comm2d, ierr )
        ENDIF
+       IF ( air_chemistry )  THEN
+          CALL MPI_ALLREDUCE( ref_chem_l, ref_chem,                            &
+                              ( nzt+1-nzb+1 ) * SIZE( ref_chem(nzb,:) ),       &
+                              MPI_REAL, MPI_SUM, comm2d, ierr )
+       ENDIF
 #else
        u_ref  = u_ref_l
        v_ref  = v_ref_l
-       IF ( humidity )       q_ref  = q_ref_l
-       IF ( .NOT. neutral )  pt_ref = pt_ref_l
+       IF ( humidity )       q_ref    = q_ref_l
+       IF ( .NOT. neutral )  pt_ref   = pt_ref_l
+       IF ( air_chemistry )  ref_chem = ref_chem_l
 #endif
 !
@@ -1427,4 +1458,8 @@
                                                           ( ny + 1 + nx + 1 ), &
                                               KIND = wp )
+       IF ( air_chemistry )                                                    &
+          ref_chem(nzb:nzt,:) = ref_chem(nzb:nzt,:) / REAL( 2.0_wp *           &
+                                                          ( ny + 1 + nx + 1 ), &
+                                                            KIND = wp )
 !
 !--    Derived from top boundary.    
@@ -1437,14 +1472,34 @@
           pt_ref(nzt+1) = pt_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx + 1 ),       &
                                                 KIND = wp )
-!
-!--    Write onto init profiles, which are used for damping
+       IF ( air_chemistry )                                                    &
+          ref_chem(nzt+1,:) = ref_chem(nzt+1,:) /                              &
+                              REAL( ( ny + 1 ) * ( nx + 1 ),KIND = wp )
+!
+!--    Write onto init profiles, which are used for damping. Also set lower
+!--    boundary condition for scalars (not required for u and v as these are 
+!--    zero at k=nzb.
        u_init = u_ref
        v_init = v_ref
-       IF ( humidity      )  q_init  = q_ref
-       IF ( .NOT. neutral )  pt_init = pt_ref
-!
-!--    Set bottom boundary condition
-       IF ( humidity      )  q_init(nzb)  = q_init(nzb+1)
-       IF ( .NOT. neutral )  pt_init(nzb) = pt_init(nzb+1)
+       IF ( humidity      )  THEN 
+          q_init      = q_ref
+          q_init(nzb) = q_init(nzb+1)
+       ENDIF
+       IF ( .NOT. neutral )  THEN
+          pt_init      = pt_ref
+          pt_init(nzb) = pt_init(nzb+1)
+       ENDIF
+
+       IF ( air_chemistry )  THEN
+          DO  n = 1, UBOUND( chem_species, 1 )
+             IF ( nest_offl%chem_from_file_t(n) )  THEN 
+                chem_species(n)%conc_pr_init(:) = ref_chem(:,n)
+                chem_species(n)%conc_pr_init(nzb) =                            &
+                                            chem_species(n)%conc_pr_init(nzb+1)
+             ENDIF
+          ENDDO
+       ENDIF
+
+       DEALLOCATE( ref_chem   )
+       DEALLOCATE( ref_chem_l )     
 !
 !--    Further, adjust Rayleigh damping height in case of time-changing conditions.
Index: palm/trunk/SOURCE/urban_surface_mod.f90
===================================================================
--- palm/trunk/SOURCE/urban_surface_mod.f90	(revision 4227)
+++ palm/trunk/SOURCE/urban_surface_mod.f90	(revision 4230)
@@ -28,4 +28,8 @@
 ! -----------------
 ! $Id$
+! Bugfix, initialize canopy resistance. Even if no green fraction is set, 
+! r_canopy must be initialized for output purposes. 
+! 
+! 4227 2019-09-10 18:04:34Z gronemeier
 ! implement new palm_date_time_mod
 ! 
@@ -3555,6 +3559,7 @@
 !--    Map values onto horizontal elemements
        DO  m = 1, surf_usm_h%ns
-             surf_usm_h%r_canopy_min(m)     = 200.0_wp !< min_canopy_resistance
-             surf_usm_h%g_d(m)              = 0.0_wp   !< canopy_resistance_coefficient
+          surf_usm_h%r_canopy(m)     = 200.0_wp !< canopy_resistance
+          surf_usm_h%r_canopy_min(m) = 200.0_wp !< min_canopy_resistance
+          surf_usm_h%g_d(m)          = 0.0_wp   !< canopy_resistance_coefficient
        ENDDO
 !
@@ -3563,6 +3568,7 @@
        DO  l = 0, 3
           DO  m = 1, surf_usm_v(l)%ns
-                surf_usm_v(l)%r_canopy_min(m)     = 200.0_wp !< min_canopy_resistance
-                surf_usm_v(l)%g_d(m)              = 0.0_wp   !< canopy_resistance_coefficient
+             surf_usm_v(l)%r_canopy(m)     = 200.0_wp !< canopy_resistance
+             surf_usm_v(l)%r_canopy_min(m) = 200.0_wp !< min_canopy_resistance
+             surf_usm_v(l)%g_d(m)          = 0.0_wp   !< canopy_resistance_coefficient
           ENDDO
        ENDDO
@@ -7585,4 +7591,5 @@
 !--           vegetated surfaces and bare soils.
               m_liq_max = m_max_depth * ( surf_usm_h%lai(m) )
+
               surf_usm_h%c_liq(m) = MIN( 1.0_wp, ( m_liq_usm_h%var_usm_1d(m) / m_liq_max )**0.67 )
 !
@@ -7987,5 +7994,5 @@
  
 
-              IF ( surf_usm_v(l)%frac(1,m) > 0.0_wp ) THEN
+              IF ( surf_usm_v(l)%frac(ind_pav_green,m) > 0.0_wp ) THEN
 !
 !--             Adapted from LSM:
