Index: palm/trunk/SOURCE/chemistry_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 4637)
+++ palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 4671)
@@ -26,4 +26,7 @@
 ! -----------------
 ! $Id$
+! Implementation of downward facing USM and LSM surfaces
+!
+! 4637 2020-08-07 07:49:33Z suehring
 ! Avoid usage of omp_lib, instead declare omp_get_thread_num explicitly
 !
@@ -684,6 +687,6 @@
                    DO  j = nys, nyn
                       match_def = surf_def_h(0)%start_index(j,i) <= surf_def_h(0)%end_index(j,i)
-                      match_lsm = surf_lsm_h%start_index(j,i)    <= surf_lsm_h%end_index(j,i)
-                      match_usm = surf_usm_h%start_index(j,i)    <= surf_usm_h%end_index(j,i)
+                      match_lsm = surf_lsm_h(0)%start_index(j,i) <= surf_lsm_h(0)%end_index(j,i)
+                      match_usm = surf_usm_h(0)%start_index(j,i) <= surf_usm_h(0)%end_index(j,i)
 
                       IF ( match_def )  THEN
@@ -692,11 +695,11 @@
                                                            surf_def_h(0)%cssws(lsp,m)
                       ELSEIF ( match_lsm  .AND.  .NOT. match_usm )  THEN
-                         m = surf_lsm_h%end_index(j,i)
+                         m = surf_lsm_h(0)%end_index(j,i)
                          chem_species(lsp)%cssws_av(j,i) = chem_species(lsp)%cssws_av(j,i) +       &
-                                                           surf_lsm_h%cssws(lsp,m)
+                                                           surf_lsm_h(0)%cssws(lsp,m)
                       ELSEIF ( match_usm )  THEN
-                         m = surf_usm_h%end_index(j,i)
+                         m = surf_usm_h(0)%end_index(j,i)
                          chem_species(lsp)%cssws_av(j,i) = chem_species(lsp)%cssws_av(j,i) +       &
-                                                           surf_usm_h%cssws(lsp,m)
+                                                           surf_usm_h(0)%cssws(lsp,m)
                       ENDIF
                    ENDDO
@@ -1468,13 +1471,15 @@
 !
 !--          No average output for now.
-             DO  m = 1, surf_lsm_h%ns
-                local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),nzb+1) =                                  &
-                                                   local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),nzb+1) &
-                                                   + surf_lsm_h%cssws(lsp,m)
+!--          !!! IT NEEDS TO RETHINK - with fully 3D structure, only lower (upper)
+!--          !!! upward facing horizontal surfaces should be taken into account here
+             DO  m = 1, surf_lsm_h(0)%ns
+                local_pf(surf_lsm_h(0)%i(m),surf_lsm_h(0)%j(m),nzb+1) =                                  &
+                                                   local_pf(surf_lsm_h(0)%i(m),surf_lsm_h(0)%j(m),nzb+1) &
+                                                   + surf_lsm_h(0)%cssws(lsp,m)
              ENDDO
-             DO  m = 1, surf_usm_h%ns
-                   local_pf(surf_usm_h%i(m),surf_usm_h%j(m),nzb+1) =                               &
-                                                   local_pf(surf_usm_h%i(m),surf_usm_h%j(m),nzb+1) &
-                                                   + surf_usm_h%cssws(lsp,m)
+             DO  m = 1, surf_usm_h(0)%ns
+                   local_pf(surf_usm_h(0)%i(m),surf_usm_h(0)%j(m),nzb+1) =                               &
+                                                   local_pf(surf_usm_h(0)%i(m),surf_usm_h(0)%j(m),nzb+1) &
+                                                   + surf_usm_h(0)%cssws(lsp,m)
              ENDDO
              grid = 'zu'
@@ -1577,13 +1582,15 @@
 !
 !--          no average for now
-             DO  m = 1, surf_usm_h%ns
-                local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) =                        &
-                                         local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) &
-                                         + surf_usm_h%cssws(lsp,m)
-             ENDDO
-             DO  m = 1, surf_lsm_h%ns
-                local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),surf_lsm_h%k(m)) =                        &
-                                         local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),surf_lsm_h%k(m)) &
-                                         + surf_lsm_h%cssws(lsp,m)
+             DO l = 0, 1
+                DO  m = 1, surf_usm_h(l)%ns
+                   local_pf(surf_usm_h(l)%i(m),surf_usm_h(l)%j(m),surf_usm_h(l)%k(m)) =            &
+                                local_pf(surf_usm_h(l)%i(m),surf_usm_h(l)%j(m),surf_usm_h(l)%k(m)) &
+                                + surf_usm_h(l)%cssws(lsp,m)
+                ENDDO
+                DO  m = 1, surf_lsm_h(l)%ns
+                   local_pf(surf_lsm_h(l)%i(m),surf_lsm_h(l)%j(m),surf_lsm_h(l)%k(m)) =            &
+                                local_pf(surf_lsm_h(l)%i(m),surf_lsm_h(l)%j(m),surf_lsm_h(l)%k(m)) &
+                                + surf_lsm_h(l)%cssws(lsp,m)
+                ENDDO
              ENDDO
              DO  l = 0, 3
@@ -2994,6 +3001,8 @@
             surf_def_h(1)%cssws(ilsp,:),                                                           &
             surf_def_h(2)%cssws(ilsp,:),                                                           &
-            surf_lsm_h%cssws(ilsp,:),                                                              &
-            surf_usm_h%cssws(ilsp,:),                                                              &
+            surf_lsm_h(0)%cssws(ilsp,:),                                                           &
+            surf_lsm_h(1)%cssws(ilsp,:),                                                           &
+            surf_usm_h(0)%cssws(ilsp,:),                                                           &
+            surf_usm_h(1)%cssws(ilsp,:),                                                           &
             surf_def_v(0)%cssws(ilsp,:),                                                           &
             surf_def_v(1)%cssws(ilsp,:),                                                           &
@@ -3115,5 +3124,6 @@
             surf_def_h(0)%cssws(ilsp,:), surf_def_h(1)%cssws(ilsp,:),                              &
             surf_def_h(2)%cssws(ilsp,:),                                                           &
-            surf_lsm_h%cssws(ilsp,:), surf_usm_h%cssws(ilsp,:),                                    &
+            surf_lsm_h(0)%cssws(ilsp,:), surf_lsm_h(1)%cssws(ilsp,:),                              &
+            surf_usm_h(0)%cssws(ilsp,:), surf_usm_h(1)%cssws(ilsp,:),                              &
             surf_def_v(0)%cssws(ilsp,:), surf_def_v(1)%cssws(ilsp,:),                              &
             surf_def_v(2)%cssws(ilsp,:), surf_def_v(3)%cssws(ilsp,:),                              &
@@ -3273,4 +3283,5 @@
     INTEGER(iwp) ::  j                   !< running index on y-axis
     INTEGER(iwp) ::  k                   !< vertical index counter
+    INTEGER(iwp) ::  l                   !< direction index counter
     INTEGER(iwp) ::  lpr                 !< running index chem spcs
     INTEGER(iwp) ::  m                   !< running index for surface elements
@@ -3328,28 +3339,28 @@
                 ENDDO
 !
-!--             Add surface fluxes
-                surf_s = surf_def_h(0)%start_index(j,i)
-                surf_e = surf_def_h(0)%end_index(j,i)
-                DO  m = surf_s, surf_e
-                   k   = surf_def_h(0)%k(m) + surf_def_h(0)%koff
-                   sums_tmp(k,tn) = sums_tmp(k,tn) + surf_def_h(0)%cssws(cs_pr_index_fl_sgs(lpr),m)
+!--             Add surface fluxes (?Is the order mandatory or could it be done in one cycle?)
+                DO l = 0, 1
+                   surf_s = surf_def_h(0)%start_index(j,i)
+                   surf_e = surf_def_h(0)%end_index(j,i)
+                   DO  m = surf_s, surf_e
+                      k   = surf_def_h(0)%k(m) + surf_def_h(0)%koff
+                      sums_tmp(k,tn) = sums_tmp(k,tn) + surf_def_h(0)%cssws(cs_pr_index_fl_sgs(lpr),m)
+                   ENDDO
                 ENDDO
-                surf_s = surf_def_h(1)%start_index(j,i)
-                surf_e = surf_def_h(1)%end_index(j,i)
-                DO  m = surf_s, surf_e
-                   k   = surf_def_h(1)%k(m)  + surf_def_h(1)%koff
-                   sums_tmp(k,tn) = sums_tmp(k,tn) + surf_def_h(1)%cssws(cs_pr_index_fl_sgs(lpr),m)
+                DO l = 0, 1
+                   surf_s = surf_lsm_h(0)%start_index(j,i)
+                   surf_e = surf_lsm_h(0)%end_index(j,i)
+                   DO  m = surf_s, surf_e
+                      k   = surf_lsm_h(0)%k(m) + surf_lsm_h(0)%koff
+                      sums_tmp(k,tn) = sums_tmp(k,tn) + surf_lsm_h(0)%cssws(cs_pr_index_fl_sgs(lpr),m)
+                   ENDDO
                 ENDDO
-                surf_s = surf_lsm_h%start_index(j,i)
-                surf_e = surf_lsm_h%end_index(j,i)
-                DO  m = surf_s, surf_e
-                   k   = surf_lsm_h%k(m) + surf_lsm_h%koff
-                   sums_tmp(k,tn) = sums_tmp(k,tn) + surf_lsm_h%cssws(cs_pr_index_fl_sgs(lpr),m)
-                ENDDO
-                surf_s = surf_usm_h%start_index(j,i)
-                surf_e = surf_usm_h%end_index(j,i)
-                DO  m = surf_s, surf_e
-                   k   = surf_usm_h%k(m) + surf_usm_h%koff
-                   sums_tmp(k,tn) = sums_tmp(k,tn) + surf_usm_h%cssws(cs_pr_index_fl_sgs(lpr),m)
+                DO l = 0, 1
+                   surf_s = surf_usm_h(0)%start_index(j,i)
+                   surf_e = surf_usm_h(0)%end_index(j,i)
+                   DO  m = surf_s, surf_e
+                      k   = surf_usm_h(0)%k(m) + surf_usm_h(0)%koff
+                      sums_tmp(k,tn) = sums_tmp(k,tn) + surf_usm_h(0)%cssws(cs_pr_index_fl_sgs(lpr),m)
+                   ENDDO
                 ENDDO
              ENDDO
@@ -3450,7 +3461,7 @@
 ! ------------
 !> Subroutine to calculate the deposition of gases and PMs. For now deposition only takes place on
-!> lsm and usm horizontal surfaces. Default surfaces are NOT considered. The deposition of particles
-!> is derived following Zhang et al., 2001, gases are deposited using the DEPAC module
-!> (van Zanten et al., 2010).
+!> lsm and usm horizontal upward faceing surfaces. Default surfaces are NOT considered.
+!> The deposition of particlesis derived following Zhang et al., 2001, gases are deposited using
+!>  the DEPAC module (van Zanten et al., 2010).
 !>
 !> @TODO: Consider deposition on vertical surfaces
@@ -3768,8 +3779,8 @@
     k = 0
 !
-!-- LSM or USM surface present at i,j:
+!-- LSM or USM horizintal upward facing surface present at i,j:
 !-- Default surfaces are NOT considered for deposition
-    match_lsm = surf_lsm_h%start_index(j,i) <= surf_lsm_h%end_index(j,i)
-    match_usm = surf_usm_h%start_index(j,i) <= surf_usm_h%end_index(j,i)
+    match_lsm = surf_lsm_h(0)%start_index(j,i) <= surf_lsm_h(0)%end_index(j,i)
+    match_usm = surf_usm_h(0)%start_index(j,i) <= surf_usm_h(0)%end_index(j,i)
 !
 !--For LSM surfaces
@@ -3777,14 +3788,14 @@
 !
 !--    Get surface element information at i,j:
-       m = surf_lsm_h%start_index(j,i)
-       k = surf_lsm_h%k(m)
+       m = surf_lsm_h(0)%start_index(j,i)
+       k = surf_lsm_h(0)%k(m)
 !
 !--    Get needed variables for surface element m
-       ustar_surf  = surf_lsm_h%us(m)
-       z0h_surf    = surf_lsm_h%z0h(m)
-       r_aero_surf = surf_lsm_h%r_a(m)
-       solar_rad   = surf_lsm_h%rad_sw_dir(m) + surf_lsm_h%rad_sw_dif(m)
-
-       lai = surf_lsm_h%lai(m)
+       ustar_surf  = surf_lsm_h(0)%us(m)
+       z0h_surf    = surf_lsm_h(0)%z0h(m)
+       r_aero_surf = surf_lsm_h(0)%r_a(m)
+       solar_rad   = surf_lsm_h(0)%rad_sw_dir(m) + surf_lsm_h(0)%rad_sw_dif(m)
+
+       lai = surf_lsm_h(0)%lai(m)
        sai = lai + 1
 !
@@ -3808,6 +3819,6 @@
 !
 !--    Get land use for i,j and assign to DEPAC lu
-       IF ( surf_lsm_h%frac(m,ind_veg_wall) > 0 )  THEN
-          luv_palm = surf_lsm_h%vegetation_type(m)
+       IF ( surf_lsm_h(0)%frac(m,ind_veg_wall) > 0 )  THEN
+          luv_palm = surf_lsm_h(0)%vegetation_type(m)
           IF ( luv_palm == ind_luv_user )  THEN
              message_string = 'No vegetation type defined. Please define vegetation type to enable deposition calculation'
@@ -3852,6 +3863,6 @@
        ENDIF
 
-       IF ( surf_lsm_h%frac(m,ind_pav_green) > 0 )  THEN
-          lup_palm = surf_lsm_h%pavement_type(m)
+       IF ( surf_lsm_h(0)%frac(m,ind_pav_green) > 0 )  THEN
+          lup_palm = surf_lsm_h(0)%pavement_type(m)
           IF ( lup_palm == ind_lup_user )  THEN
              message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation'
@@ -3890,6 +3901,6 @@
        ENDIF
 
-       IF ( surf_lsm_h%frac(m,ind_wat_win) > 0 )  THEN
-          luw_palm = surf_lsm_h%water_type(m)
+       IF ( surf_lsm_h(0)%frac(m,ind_wat_win) > 0 )  THEN
+          luw_palm = surf_lsm_h(0)%water_type(m)
           IF ( luw_palm == ind_luw_user )  THEN
              message_string = 'No water type defined. Please define water type to enable deposition calculation'
@@ -3909,5 +3920,5 @@
 !
 !--    Set wetness indicator to dry or wet for lsm vegetation or pavement
-       IF ( surf_lsm_h%c_liq(m) > 0 )  THEN
+       IF ( surf_lsm_h(0)%c_liq(m) > 0 )  THEN
           nwet = 1
        ELSE
@@ -3950,5 +3961,5 @@
 !
 !--    Vegetation
-       IF ( surf_lsm_h%frac(m,ind_veg_wall) > 0 )  THEN
+       IF ( surf_lsm_h(0)%frac(m,ind_veg_wall) > 0 )  THEN
 
 !
@@ -4084,5 +4095,5 @@
 !
 !--    Pavement
-       IF ( surf_lsm_h%frac(m,ind_pav_green) > 0 )  THEN
+       IF ( surf_lsm_h(0)%frac(m,ind_pav_green) > 0 )  THEN
 !
 !--       No vegetation on pavements:
@@ -4208,5 +4219,5 @@
 !
 !--    Water
-       IF ( surf_lsm_h%frac(m,ind_wat_win) > 0 )  THEN
+       IF ( surf_lsm_h(0)%frac(m,ind_wat_win) > 0 )  THEN
 !
 !--       No vegetation on water:
@@ -4339,7 +4350,7 @@
        DO  lsp = 1, nspec
 
-          bud(lsp) = surf_lsm_h%frac(m,ind_veg_wall)  * bud_luv(lsp) +                             &
-                     surf_lsm_h%frac(m,ind_pav_green) * bud_lup(lsp) +                             &
-                     surf_lsm_h%frac(m,ind_wat_win)   * bud_luw(lsp)
+          bud(lsp) = surf_lsm_h(0)%frac(m,ind_veg_wall)  * bud_luv(lsp) +                             &
+                     surf_lsm_h(0)%frac(m,ind_pav_green) * bud_lup(lsp) +                             &
+                     surf_lsm_h(0)%frac(m,ind_wat_win)   * bud_luw(lsp)
 !
 !--       Compute new concentration:
@@ -4356,13 +4367,13 @@
 !
 !--    Get surface element information at i,j:
-       m = surf_usm_h%start_index(j,i)
-       k = surf_usm_h%k(m)
+       m = surf_usm_h(0)%start_index(j,i)
+       k = surf_usm_h(0)%k(m)
 !
 !--    Get needed variables for surface element m
-       ustar_surf  = surf_usm_h%us(m)
-       z0h_surf    = surf_usm_h%z0h(m)
-       r_aero_surf = surf_usm_h%r_a(m)
-       solar_rad   = surf_usm_h%rad_sw_dir(m) + surf_usm_h%rad_sw_dif(m)
-       lai = surf_usm_h%lai(m)
+       ustar_surf  = surf_usm_h(0)%us(m)
+       z0h_surf    = surf_usm_h(0)%z0h(m)
+       r_aero_surf = surf_usm_h(0)%r_a(m)
+       solar_rad   = surf_usm_h(0)%rad_sw_dir(m) + surf_usm_h(0)%rad_sw_dif(m)
+       lai = surf_usm_h(0)%lai(m)
        sai = lai + 1
 !
@@ -4386,5 +4397,5 @@
 !
 !--    Get land use for i,j and assign to DEPAC lu
-       IF ( surf_usm_h%frac(m,ind_pav_green) > 0 )  THEN
+       IF ( surf_usm_h(0)%frac(m,ind_pav_green) > 0 )  THEN
 !
 !--       For green urban surfaces (e.g. green roofs assume LU short grass
@@ -4432,5 +4443,5 @@
        ENDIF
 
-       IF ( surf_usm_h%frac(m,ind_veg_wall) > 0 )  THEN
+       IF ( surf_usm_h(0)%frac(m,ind_veg_wall) > 0 )  THEN
 !
 !--       For walls in USM assume concrete walls/roofs,
@@ -4473,5 +4484,5 @@
        ENDIF
 
-       IF ( surf_usm_h%frac(m,ind_wat_win) > 0 )  THEN
+       IF ( surf_usm_h(0)%frac(m,ind_wat_win) > 0 )  THEN
 !
 !--       For windows in USM assume metal as this is as close as we get, assumed LU class desert as
@@ -4516,5 +4527,5 @@
 !--    @TODO: Activate these lines as soon as new ebsolver branch is merged:
 !--    Set wetness indicator to dry or wet for usm vegetation or pavement
-       !IF ( surf_usm_h%c_liq(m) > 0 )  THEN
+       !IF ( surf_usm_h(0)%c_liq(m) > 0 )  THEN
        !   nwet = 1
        !ELSE
@@ -4557,5 +4568,5 @@
 
 !--    Walls/roofs
-       IF ( surf_usm_h%frac(m,ind_veg_wall) > 0 )  THEN
+       IF ( surf_usm_h(0)%frac(m,ind_veg_wall) > 0 )  THEN
 !
 !--       No vegetation on non-green walls:
@@ -4684,5 +4695,5 @@
 !
 !--    Green usm surfaces
-       IF ( surf_usm_h%frac(m,ind_pav_green) > 0 )  THEN
+       IF ( surf_usm_h(0)%frac(m,ind_pav_green) > 0 )  THEN
 
 !
@@ -4817,5 +4828,5 @@
 !
 !--    Windows
-       IF ( surf_usm_h%frac(m,ind_wat_win) > 0 )  THEN
+       IF ( surf_usm_h(0)%frac(m,ind_wat_win) > 0 )  THEN
 !
 !--       No vegetation on windows:
@@ -4947,7 +4958,7 @@
 
 
-          bud(lsp) = surf_usm_h%frac(m,ind_veg_wall)  * bud_luu(lsp) +                             &
-                     surf_usm_h%frac(m,ind_pav_green) * bud_lug(lsp) +                             &
-                     surf_usm_h%frac(m,ind_wat_win)   * bud_lud(lsp)
+          bud(lsp) = surf_usm_h(0)%frac(m,ind_veg_wall)  * bud_luu(lsp) +                             &
+                     surf_usm_h(0)%frac(m,ind_pav_green) * bud_lug(lsp) +                             &
+                     surf_usm_h(0)%frac(m,ind_wat_win)   * bud_lud(lsp)
 !
 !--       Compute new concentration
