Index: palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
===================================================================
--- palm/trunk/SOURCE/surface_layer_fluxes_mod.f90	(revision 2547)
+++ palm/trunk/SOURCE/surface_layer_fluxes_mod.f90	(revision 2696)
@@ -1,5 +1,5 @@
 !> @file surface_layer_fluxes_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
@@ -18,4 +18,5 @@
 !
 !------------------------------------------------------------------------------!
+!
 ! Current revisions:
 ! ------------------
@@ -25,4 +26,10 @@
 ! -----------------
 ! $Id$
+! - Implementation of chemistry module (FK)
+! - Added calculation of pt1 and qv1 for all surface types. Added calculation of
+!   pt_surface for default-type surfaces (BM)
+! - Add flag to disable computation of qsws in case of urban surface (MS)
+! 
+! 2547 2017-10-16 12:41:56Z schwenkel
 ! extended by cloud_droplets option
 ! 
@@ -188,4 +195,5 @@
 !> @todo (re)move large_scale_forcing actions
 !> @todo check/optimize OpenMP directives
+!> @todo simplify if conditions (which flux need to be computed in which case)
 !------------------------------------------------------------------------------!
  MODULE surface_layer_fluxes_mod
@@ -195,4 +203,9 @@
                drho_air_zw, rho_air_zw
 
+#if defined( __chem )
+    USE chem_modules,                                                          &
+        ONLY:  constant_csflux, nvar
+#endif
+
     USE cloud_parameters,                                                      &
         ONLY:  l_d_cp, pt_d_t
@@ -204,7 +217,8 @@
 
     USE control_parameters,                                                    &
-        ONLY:  cloud_droplets, cloud_physics, constant_heatflux,               &
-               constant_scalarflux, constant_waterflux, coupling_mode, g,      &
-               humidity, ibc_e_b, ibc_pt_b, initializing_actions, kappa,       &
+        ONLY:  air_chemistry, cloud_droplets, cloud_physics,                   &
+               constant_heatflux, constant_scalarflux,                         &     
+               constant_waterflux, coupling_mode, g, humidity, ibc_e_b,        &
+               ibc_pt_b, initializing_actions, kappa,                          &
                intermediate_timestep_count, intermediate_timestep_count_max,   &
                land_surface, large_scale_forcing, lsf_surf,                    &
@@ -291,30 +305,45 @@
        downward      = .FALSE.
 !
-!--    In case cloud physics is used, it is required to derive potential
-!--    temperature and specific humidity at first grid level from the fields pt 
-!--    and q
-       IF ( cloud_physics  .OR.  cloud_droplets )  THEN
-!
-!--       First call for horizontal default-type surfaces (l=0 - upward facing, 
-!--       l=1 - downward facing)
-          DO  l = 0, 1
-             IF ( surf_def_h(l)%ns >= 1  )  THEN
-                surf => surf_def_h(l)
-                CALL calc_pt_q
-             ENDIF
-          ENDDO
-!
-!--       Call for natural-type surfaces
-          IF ( surf_lsm_h%ns >= 1  )  THEN
-             surf => surf_lsm_h
+!--    Derive potential temperature and specific humidity at first grid level 
+!--    from the fields pt and q
+!
+!--    First call for horizontal default-type surfaces (l=0 - upward facing, 
+!--    l=1 - downward facing)
+       DO  l = 0, 1
+          IF ( surf_def_h(l)%ns >= 1  )  THEN
+             surf => surf_def_h(l)
              CALL calc_pt_q
-          ENDIF
-!
-!--       Call for urban-type surfaces
-          IF ( surf_usm_h%ns >= 1  )  THEN
-             surf => surf_usm_h
+             IF ( .NOT. neutral )  CALL calc_pt_surface
+          ENDIF
+       ENDDO
+!
+!--    Call for natural-type horizontal surfaces
+       IF ( surf_lsm_h%ns >= 1  )  THEN
+          surf => surf_lsm_h
+          CALL calc_pt_q
+       ENDIF
+
+!
+!--    Call for urban-type horizontal surfaces
+       IF ( surf_usm_h%ns >= 1  )  THEN
+          surf => surf_usm_h
+          CALL calc_pt_q
+       ENDIF
+
+!
+!--    Call for natural-type vertical surfaces
+       DO  l = 0, 3
+          IF ( surf_lsm_v(l)%ns >= 1  )  THEN
+             surf => surf_lsm_v(l)
              CALL calc_pt_q
           ENDIF
-       ENDIF
+
+!
+!--    Call for urban-type vertical surfaces
+          IF ( surf_usm_v(l)%ns >= 1  )  THEN
+             surf => surf_usm_v(l)
+             CALL calc_pt_q
+          ENDIF
+       ENDDO
 
 !
@@ -923,5 +952,4 @@
 !--       However, this do not apply for downward-facing walls. To mask this, 
 !--       use ibit, which checks for upward/downward-facing surfaces. 
-
           surf%uvw_abs(m) = SQRT(                                              &
                               ( 0.5_wp * (   u(k,j,i)   + u(k,j,i+1)           &
@@ -1538,4 +1566,11 @@
              surf%pt1(m) = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
              surf%qv1(m) = q(k,j,i) 
+          ELSE
+             surf%pt1(m) = pt(k,j,i)
+             IF ( humidity )  THEN
+                surf%qv1(m) = q(k,j,i)
+             ELSE
+                surf%qv1(m) = 0.0_wp
+             ENDIF
           ENDIF 
 
@@ -1543,4 +1578,29 @@
 
     END SUBROUTINE calc_pt_q
+
+
+!
+!-- Calculate potential temperature and specific humidity at first grid level
+!-- ( only for upward-facing surfs )
+    SUBROUTINE calc_pt_surface
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  koff    !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls) 
+       INTEGER(iwp) ::  m       !< loop variable over all horizontal surf elements 
+       
+       koff = surf%koff
+       !$OMP PARALLEL DO PRIVATE( i, j, k )
+       DO  m = 1, surf%ns 
+
+          i   = surf%i(m)            
+          j   = surf%j(m)
+          k   = surf%k(m)
+
+          surf%pt_surface(m) = pt(k+koff,j,i)
+
+       ENDDO
+
+    END SUBROUTINE calc_pt_surface
 
 !
@@ -1552,5 +1612,5 @@
 
        INTEGER(iwp)  ::  m       !< loop variable over all horizontal surf elements 
-
+       INTEGER(iwp)  ::  lsp     !< running index for chemical species
 ! 
 !--    Compute theta* at horizontal surfaces
@@ -1592,36 +1652,16 @@
           ENDIF
 
-          IF ( cloud_physics  .OR.  cloud_droplets )  THEN
-             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
-             DO  m = 1, surf%ns   
-
-                i   = surf%i(m)            
-                j   = surf%j(m)
-                k   = surf%k(m)
-
-                z_mo = surf%z_mo(m)
-
-                surf%ts(m) = kappa * ( surf%pt1(m) - pt(k-1,j,i) )             &
-                                     / ( LOG( z_mo / surf%z0h(m) )             &
-                                         - psi_h( z_mo / surf%ol(m) )          &
-                                         + psi_h( surf%z0h(m) / surf%ol(m) ) )
-
-             ENDDO
-          ELSE
-             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
-             DO  m = 1, surf%ns   
-
-                i   = surf%i(m)            
-                j   = surf%j(m)
-                k   = surf%k(m)
-
-                z_mo = surf%z_mo(m)
-
-                surf%ts(m) = kappa * ( pt(k,j,i) - pt(k-1,j,i) )               &
-                                     / ( LOG( z_mo / surf%z0h(m) )             &
-                                         - psi_h( z_mo / surf%ol(m) )          &
-                                         + psi_h( surf%z0h(m) / surf%ol(m) ) )
-             ENDDO
-          ENDIF
+          !$OMP PARALLEL DO PRIVATE( z_mo )
+          DO  m = 1, surf%ns   
+
+             z_mo = surf%z_mo(m)
+
+             surf%ts(m) = kappa * ( surf%pt1(m) - surf%pt_surface(m) )      &
+                                  / ( LOG( z_mo / surf%z0h(m) )             &
+                                      - psi_h( z_mo / surf%ol(m) )          &
+                                      + psi_h( surf%z0h(m) / surf%ol(m) ) )
+
+          ENDDO
+
        ENDIF
 ! 
@@ -1769,4 +1809,36 @@
 
 !
+!--    If required compute cs* (chemical species)
+#if defined( __chem )
+       IF ( air_chemistry  )  THEN  
+!
+!--       At horizontal surfaces                             
+          DO  lsp = 1, nvar
+             IF ( constant_csflux(lsp)  .AND.  .NOT.  surf_vertical )  THEN
+!--             For a given chemical species' flux in the surface layer
+                !$OMP PARALLEL DO PRIVATE( i, j )
+                DO  m = 1, surf%ns  
+                   i   = surf%i(m)            
+                   j   = surf%j(m)
+                   surf%css(lsp,m) = -surf%cssws(lsp,m) / ( surf%us(m) + 1E-30_wp )
+                ENDDO
+             ENDIF
+          ENDDO
+!
+!--       At vertical surfaces
+          IF ( surf_vertical )  THEN
+             DO  lsp = 1, nvar
+                !$OMP PARALLEL DO PRIVATE( i, j )
+                DO  m = 1, surf%ns  
+                   i   = surf%i(m)            
+                   j   = surf%j(m)
+                   surf%css(lsp,m) = -surf%cssws(lsp,m) / ( surf%us(m) + 1E-30_wp )
+                ENDDO
+             ENDDO
+          ENDIF
+       ENDIF
+#endif
+
+!
 !--    If required compute qc* and nc*
        IF ( cloud_physics  .AND.  microphysics_morrison  .AND.                 &
@@ -1830,5 +1902,6 @@
        IMPLICIT NONE
 
-       INTEGER(iwp)  ::  m       !< loop variable over all horizontal surf elements 
+       INTEGER(iwp)  ::  m       !< loop variable over all horizontal surf elements
+       INTEGER(iwp)  ::  lsp     !< running index for chemical species
 
        REAL(wp)                            ::  dum     !< dummy to precalculate logarithm
@@ -1883,6 +1956,4 @@
 
                 surf%usws(m) = kappa * u(k,j,i) / LOG( z_mo / surf%z0(m) )
-!
-!--        To Do: Is the sign correct???
                 surf%usws(m) = surf%usws(m) * surf%us(m) * rho_air_zw(k)
 
@@ -1929,5 +2000,4 @@
 
                 surf%vsws(m) = kappa * v(k,j,i) / LOG( z_mo / surf%z0(m) )
- 
                 surf%vsws(m) = surf%vsws(m) * surf%us(m) * rho_air_zw(k)
              ENDDO
@@ -1937,5 +2007,5 @@
           IF (  .NOT.  constant_heatflux  .AND.  ( ( time_since_reference_point&
                <=  skip_time_do_lsm  .AND. simulated_time > 0.0_wp ) .OR.      &
-               .NOT.  land_surface )  .AND.  .NOT.  urban_surface  .AND.       &
+               .NOT.  land_surface )  .AND.  .NOT. urban_surface  .AND.        &
                .NOT. downward )  THEN
              !$OMP PARALLEL DO PRIVATE( i, j, k )
@@ -1951,6 +2021,6 @@
           IF (  .NOT.  constant_waterflux  .AND.  humidity  .AND.              &
                ( ( time_since_reference_point <= skip_time_do_lsm  .AND.       &
-               simulated_time > 0.0_wp ) .OR.  .NOT.  land_surface )  .AND.    &
-               .NOT. downward )  THEN
+               simulated_time > 0.0_wp ) .OR.  .NOT.  land_surface  )  .AND.   &
+               .NOT.  urban_surface  .AND.  .NOT. downward )  THEN
              !$OMP PARALLEL DO PRIVATE( i, j, k )
              DO  m = 1, surf%ns 
@@ -1974,4 +2044,20 @@
              ENDDO
           ENDIF   
+!
+!--       Compute the vertical chemical species' flux
+#if defined( __chem )
+          DO  lsp = 1, nvar
+             IF (  .NOT.  constant_csflux(lsp)  .AND.  air_chemistry  .AND.    &
+                   .NOT.  downward )  THEN
+                !$OMP PARALLEL DO PRIVATE( i, j )
+                DO  m = 1, surf%ns   
+                   i     = surf%i(m)            
+                   j     = surf%j(m)
+                   surf%cssws(lsp,m) = -surf%css(lsp,m) * surf%us(m)
+                ENDDO
+             ENDIF   
+          ENDDO
+#endif
+
 !
 !--       Compute (turbulent) fluxes of cloud water content and cloud drop conc.
