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.
