Index: palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- palm/trunk/SOURCE/flow_statistics.f90	(revision 4646)
+++ palm/trunk/SOURCE/flow_statistics.f90	(revision 4671)
@@ -24,4 +24,7 @@
 ! -----------------
 ! $Id$
+! Implementation of downward facing USM and LSM surfaces
+!
+! 4646 2020-08-24 16:02:40Z raasch
 ! file re-formatted to follow the PALM coding standard
 !
@@ -523,5 +526,5 @@
        !$ACC PRESENT(e, u, v, w, km, kh, p, pt) &
        !$ACC PRESENT(ddx, ddy, ddzu, ddzw) &
-       !$ACC PRESENT(surf_def_h(0), surf_lsm_h, surf_usm_h) &
+       !$ACC PRESENT(surf_def_h(0), surf_lsm_h(0), surf_usm_h(0)) &
        !$ACC PRESENT(sums_l)
        DO  i = nxl, nxr
@@ -685,26 +688,26 @@
                                             surf_def_h(0)%pt_surface(m) * rmask(j,i,sr)
              ENDIF
-             IF ( surf_lsm_h%end_index(j,i) >= surf_lsm_h%start_index(j,i) )  THEN
-                m = surf_lsm_h%start_index(j,i)
+             IF ( surf_lsm_h(0)%end_index(j,i) >= surf_lsm_h(0)%start_index(j,i) )  THEN
+                m = surf_lsm_h(0)%start_index(j,i)
                 !$ACC ATOMIC
                 sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +                                &
-                                           surf_lsm_h%us(m)   * rmask(j,i,sr)
+                                           surf_lsm_h(0)%us(m)   * rmask(j,i,sr)
                 !$ACC ATOMIC
                 sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) +                              &
-                                           surf_lsm_h%usws(m) * rmask(j,i,sr)
+                                           surf_lsm_h(0)%usws(m) * rmask(j,i,sr)
                 !$ACC ATOMIC
                 sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) +                              &
-                                           surf_lsm_h%vsws(m) * rmask(j,i,sr)
+                                           surf_lsm_h(0)%vsws(m) * rmask(j,i,sr)
                 !$ACC ATOMIC
                 sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) +                              &
-                                           surf_lsm_h%ts(m)   * rmask(j,i,sr)
+                                           surf_lsm_h(0)%ts(m)   * rmask(j,i,sr)
 #ifndef _OPENACC
                 IF ( humidity )  THEN
                    sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +                         &
-                                               surf_lsm_h%qs(m)   * rmask(j,i,sr)
+                                               surf_lsm_h(0)%qs(m)   * rmask(j,i,sr)
                 ENDIF
                 IF ( passive_scalar )  THEN
                    sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) +                         &
-                                               surf_lsm_h%ss(m)   * rmask(j,i,sr)
+                                               surf_lsm_h(0)%ss(m)   * rmask(j,i,sr)
                 ENDIF
 #endif
@@ -713,28 +716,28 @@
                 !$ACC ATOMIC
                 sums_l(nzb+14,pr_palm,tn) = sums_l(nzb+14,pr_palm,tn) +                            &
-                                            surf_lsm_h%pt_surface(m) * rmask(j,i,sr)
+                                            surf_lsm_h(0)%pt_surface(m) * rmask(j,i,sr)
              ENDIF
-             IF ( surf_usm_h%end_index(j,i) >= surf_usm_h%start_index(j,i) )  THEN
-                m = surf_usm_h%start_index(j,i)
+             IF ( surf_usm_h(0)%end_index(j,i) >= surf_usm_h(0)%start_index(j,i) )  THEN
+                m = surf_usm_h(0)%start_index(j,i)
                 !$ACC ATOMIC
                 sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +                                &
-                                           surf_usm_h%us(m)   * rmask(j,i,sr)
+                                           surf_usm_h(0)%us(m)   * rmask(j,i,sr)
                 !$ACC ATOMIC
                 sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) +                              &
-                                           surf_usm_h%usws(m) * rmask(j,i,sr)
+                                           surf_usm_h(0)%usws(m) * rmask(j,i,sr)
                 !$ACC ATOMIC
                 sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) +                              &
-                                           surf_usm_h%vsws(m) * rmask(j,i,sr)
+                                           surf_usm_h(0)%vsws(m) * rmask(j,i,sr)
                 !$ACC ATOMIC
                 sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) +                              &
-                                           surf_usm_h%ts(m)   * rmask(j,i,sr)
+                                           surf_usm_h(0)%ts(m)   * rmask(j,i,sr)
 #ifndef _OPENACC
                 IF ( humidity )  THEN
                    sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +                         &
-                                               surf_usm_h%qs(m)   * rmask(j,i,sr)
+                                               surf_usm_h(0)%qs(m)   * rmask(j,i,sr)
                 ENDIF
                 IF ( passive_scalar )  THEN
                    sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) +                         &
-                                               surf_usm_h%ss(m) * rmask(j,i,sr)
+                                               surf_usm_h(0)%ss(m) * rmask(j,i,sr)
                 ENDIF
 #endif
@@ -743,5 +746,5 @@
                 !$ACC ATOMIC
                 sums_l(nzb+14,pr_palm,tn) = sums_l(nzb+14,pr_palm,tn) +                            &
-                                            surf_usm_h%pt_surface(m)  * rmask(j,i,sr)
+                                            surf_usm_h(0)%pt_surface(m)  * rmask(j,i,sr)
              ENDIF
           ENDDO !j-loop
@@ -816,5 +819,5 @@
        !$ACC PRESENT(wall_flags_total_0, rmask, ddzu, rho_air_zw, hom(:,1,1:4,sr)) &
        !$ACC PRESENT(heatflux_output_conversion, momentumflux_output_conversion) &
-       !$ACC PRESENT(surf_def_h(0:2), surf_lsm_h, surf_usm_h) &
+       !$ACC PRESENT(surf_def_h(0:2), surf_lsm_h(0:1), surf_usm_h(0:1) &
        !$ACC PRESENT(sums_l)
        DO  i = nxl, nxr
@@ -988,115 +991,114 @@
 
                    ENDIF
-                ENDDO
-                IF ( surf_lsm_h%end_index(j,i) >= surf_lsm_h%start_index(j,i) )  THEN
-                   m = surf_lsm_h%start_index(j,i)
-                   !$ACC ATOMIC
-                   sums_l(nzb,12,tn) = sums_l(nzb,12,tn) +                                         &
-                                       momentumflux_output_conversion(nzb) *                       &
-                                       surf_lsm_h%usws(m) * rmask(j,i,sr)     ! w"u"
-                   !$ACC ATOMIC
-                   sums_l(nzb,14,tn) = sums_l(nzb,14,tn) +                                         &
-                                       momentumflux_output_conversion(nzb) *                       &
-                                       surf_lsm_h%vsws(m) * rmask(j,i,sr)     ! w"v"
-                   !$ACC ATOMIC
-                   sums_l(nzb,16,tn) = sums_l(nzb,16,tn) +                                         &
-                                       heatflux_output_conversion(nzb) *                           &
-                                       surf_lsm_h%shf(m)  * rmask(j,i,sr)     ! w"pt"
+                   IF ( surf_lsm_h(l)%end_index(j,i) >= surf_lsm_h(l)%start_index(j,i) )  THEN
+                      m = surf_lsm_h(l)%start_index(j,i)
+                      !$ACC ATOMIC
+                      sums_l(nzb,12,tn) = sums_l(nzb,12,tn) +                                         &
+                                          momentumflux_output_conversion(nzb) *                       &
+                                          surf_lsm_h(l)%usws(m) * rmask(j,i,sr)     ! w"u"
+                      !$ACC ATOMIC
+                      sums_l(nzb,14,tn) = sums_l(nzb,14,tn) +                                         &
+                                          momentumflux_output_conversion(nzb) *                       &
+                                          surf_lsm_h(l)%vsws(m) * rmask(j,i,sr)     ! w"v"
+                      !$ACC ATOMIC
+                      sums_l(nzb,16,tn) = sums_l(nzb,16,tn) +                                         &
+                                          heatflux_output_conversion(nzb) *                           &
+                                          surf_lsm_h(l)%shf(m)  * rmask(j,i,sr)     ! w"pt"
 #if 0
-                   sums_l(nzb,58,tn) = sums_l(nzb,58,tn) +                                         &
-                                       0.0_wp * rmask(j,i,sr)        ! u"pt"
-                   sums_l(nzb,61,tn) = sums_l(nzb,61,tn) +                                         &
-                                       0.0_wp * rmask(j,i,sr)        ! v"pt"
+                      sums_l(nzb,58,tn) = sums_l(nzb,58,tn) +                                         &
+                                          0.0_wp * rmask(j,i,sr)        ! u"pt"
+                      sums_l(nzb,61,tn) = sums_l(nzb,61,tn) +                                         &
+                                          0.0_wp * rmask(j,i,sr)        ! v"pt"
 #endif
 #ifndef _OPENACC
-                   IF ( ocean_mode )  THEN
-                      sums_l(nzb,65,tn) = sums_l(nzb,65,tn) +                                      &
-                                          surf_lsm_h%sasws(m) * rmask(j,i,sr)  ! w"sa"
+                      IF ( ocean_mode )  THEN
+                         sums_l(nzb,65,tn) = sums_l(nzb,65,tn) +                                      &
+                                             surf_lsm_h(l)%sasws(m) * rmask(j,i,sr)  ! w"sa"
+                      ENDIF
+                      IF ( humidity )  THEN
+                         sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                                      &
+                                             waterflux_output_conversion(nzb) *                       &
+                                             surf_lsm_h(l)%qsws(m) * rmask(j,i,sr)  ! w"q" (w"qv")
+                         sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                                    &
+                                             ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) * surf_lsm_h(l)%shf(m) +  &
+                                               0.61_wp * pt(nzb,j,i) * surf_lsm_h(l)%qsws(m) )           &
+                                             * heatflux_output_conversion(nzb)
+                         IF ( cloud_droplets )  THEN
+                            sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                                 &
+                                                ( 1.0_wp + 0.61_wp * q(nzb,j,i) -                     &
+                                                  ql(nzb,j,i) ) * surf_lsm_h(l)%shf(m) +                 &
+                                                  0.61_wp * pt(nzb,j,i) * surf_lsm_h(l)%qsws(m) )        &
+                                                * heatflux_output_conversion(nzb)
+                         ENDIF
+                         IF ( bulk_cloud_model )  THEN
+!
+!--                         Formula does not work if ql(nzb) /= 0.0
+                            sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                                   &
+                                                waterflux_output_conversion(nzb) *                    &
+                                                surf_lsm_h(l)%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
+                         ENDIF
+                      ENDIF
+                      IF ( passive_scalar )  THEN
+                         sums_l(nzb,117,tn) = sums_l(nzb,117,tn) +                                    &
+                                              surf_lsm_h(l)%ssws(m) * rmask(j,i,sr) ! w"s"
+                      ENDIF
+#endif
+
                    ENDIF
-                   IF ( humidity )  THEN
-                      sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                                      &
-                                          waterflux_output_conversion(nzb) *                       &
-                                          surf_lsm_h%qsws(m) * rmask(j,i,sr)  ! w"q" (w"qv")
-                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                                    &
-                                          ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) * surf_lsm_h%shf(m) +  &
-                                            0.61_wp * pt(nzb,j,i) * surf_lsm_h%qsws(m) )           &
-                                          * heatflux_output_conversion(nzb)
-                      IF ( cloud_droplets )  THEN
-                         sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                                 &
-                                             ( 1.0_wp + 0.61_wp * q(nzb,j,i) -                     &
-                                               ql(nzb,j,i) ) * surf_lsm_h%shf(m) +                 &
-                                               0.61_wp * pt(nzb,j,i) * surf_lsm_h%qsws(m) )        &
-                                             * heatflux_output_conversion(nzb)
-                      ENDIF
-                      IF ( bulk_cloud_model )  THEN
-!
-!--                      Formula does not work if ql(nzb) /= 0.0
-                         sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                                   &
-                                             waterflux_output_conversion(nzb) *                    &
-                                             surf_lsm_h%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
-                      ENDIF
-                   ENDIF
-                   IF ( passive_scalar )  THEN
-                      sums_l(nzb,117,tn) = sums_l(nzb,117,tn) +                                    &
-                                           surf_lsm_h%ssws(m) * rmask(j,i,sr) ! w"s"
-                   ENDIF
-#endif
-
-                ENDIF
-                IF ( surf_usm_h%end_index(j,i) >= surf_usm_h%start_index(j,i) )  THEN
-                   m = surf_usm_h%start_index(j,i)
-                   !$ACC ATOMIC
-                   sums_l(nzb,12,tn) = sums_l(nzb,12,tn) +                                         &
-                                       momentumflux_output_conversion(nzb) *                       &
-                                       surf_usm_h%usws(m) * rmask(j,i,sr)                ! w"u"
-                   !$ACC ATOMIC
-                   sums_l(nzb,14,tn) = sums_l(nzb,14,tn) +                                         &
-                                       momentumflux_output_conversion(nzb) *                       &
-                                       surf_usm_h%vsws(m) * rmask(j,i,sr)                ! w"v"
-                   !$ACC ATOMIC
-                   sums_l(nzb,16,tn) = sums_l(nzb,16,tn) +                                         &
-                                       heatflux_output_conversion(nzb) *                           &
-                                       surf_usm_h%shf(m)  * rmask(j,i,sr)                ! w"pt"
+                   IF ( surf_usm_h(l)%end_index(j,i) >= surf_usm_h(l)%start_index(j,i) )  THEN
+                      m = surf_usm_h(l)%start_index(j,i)
+                      !$ACC ATOMIC
+                      sums_l(nzb,12,tn) = sums_l(nzb,12,tn) +                                         &
+                                          momentumflux_output_conversion(nzb) *                       &
+                                          surf_usm_h(l)%usws(m) * rmask(j,i,sr)                ! w"u"
+                      !$ACC ATOMIC
+                      sums_l(nzb,14,tn) = sums_l(nzb,14,tn) +                                         &
+                                          momentumflux_output_conversion(nzb) *                       &
+                                          surf_usm_h(l)%vsws(m) * rmask(j,i,sr)                ! w"v"
+                      !$ACC ATOMIC
+                      sums_l(nzb,16,tn) = sums_l(nzb,16,tn) +                                         &
+                                          heatflux_output_conversion(nzb) *                           &
+                                          surf_usm_h(l)%shf(m)  * rmask(j,i,sr)                ! w"pt"
 #if 0
-                   sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + 0.0_wp * rmask(j,i,sr)        ! u"pt"
-                   sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + 0.0_wp * rmask(j,i,sr)        ! v"pt"
+                      sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + 0.0_wp * rmask(j,i,sr)        ! u"pt"
+                      sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + 0.0_wp * rmask(j,i,sr)        ! v"pt"
 #endif
 #ifndef _OPENACC
-                   IF ( ocean_mode )  THEN
-                      sums_l(nzb,65,tn) = sums_l(nzb,65,tn) +                                      &
-                                          surf_usm_h%sasws(m) * rmask(j,i,sr)            ! w"sa"
+                      IF ( ocean_mode )  THEN
+                         sums_l(nzb,65,tn) = sums_l(nzb,65,tn) +                                      &
+                                             surf_usm_h(l)%sasws(m) * rmask(j,i,sr)            ! w"sa"
+                      ENDIF
+                      IF ( humidity )  THEN
+                         sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                                      &
+                                             waterflux_output_conversion(nzb) *                       &
+                                             surf_usm_h(l)%qsws(m) * rmask(j,i,sr)             ! w"q" (w"qv")
+                         sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                                    &
+                                             ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) *                      &
+                                             surf_usm_h(l)%shf(m) + 0.61_wp * pt(nzb,j,i) *              &
+                                             surf_usm_h(l)%qsws(m)  )                                    &
+                                             * heatflux_output_conversion(nzb)
+                         IF ( cloud_droplets )  THEN
+                            sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                                 &
+                                                ( 1.0_wp + 0.61_wp * q(nzb,j,i) -                     &
+                                                 ql(nzb,j,i) ) * surf_usm_h(l)%shf(m) +                  &
+                                                 0.61_wp * pt(nzb,j,i) * surf_usm_h(l)%qsws(m) )         &
+                                                * heatflux_output_conversion(nzb)
+                         ENDIF
+                         IF ( bulk_cloud_model )  THEN
+!
+!--                         Formula does not work if ql(nzb) /= 0.0
+                            sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                                   &
+                                                waterflux_output_conversion(nzb) *                    &
+                                                surf_usm_h(l)%qsws(m) * rmask(j,i,sr)          ! w"q" (w"qv")
+                         ENDIF
+                      ENDIF
+                      IF ( passive_scalar )  THEN
+                         sums_l(nzb,117,tn) = sums_l(nzb,117,tn) +                                    &
+                                              surf_usm_h(l)%ssws(m) * rmask(j,i,sr)            ! w"s"
+                      ENDIF
+#endif
+
                    ENDIF
-                   IF ( humidity )  THEN
-                      sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                                      &
-                                          waterflux_output_conversion(nzb) *                       &
-                                          surf_usm_h%qsws(m) * rmask(j,i,sr)             ! w"q" (w"qv")
-                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                                    &
-                                          ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) *                      &
-                                          surf_usm_h%shf(m) + 0.61_wp * pt(nzb,j,i) *              &
-                                          surf_usm_h%qsws(m)  )                                    &
-                                          * heatflux_output_conversion(nzb)
-                      IF ( cloud_droplets )  THEN
-                         sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                                 &
-                                             ( 1.0_wp + 0.61_wp * q(nzb,j,i) -                     &
-                                              ql(nzb,j,i) ) * surf_usm_h%shf(m) +                  &
-                                              0.61_wp * pt(nzb,j,i) * surf_usm_h%qsws(m) )         &
-                                             * heatflux_output_conversion(nzb)
-                      ENDIF
-                      IF ( bulk_cloud_model )  THEN
-!
-!--                      Formula does not work if ql(nzb) /= 0.0
-                         sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                                   &
-                                             waterflux_output_conversion(nzb) *                    &
-                                             surf_usm_h%qsws(m) * rmask(j,i,sr)          ! w"q" (w"qv")
-                      ENDIF
-                   ENDIF
-                   IF ( passive_scalar )  THEN
-                      sums_l(nzb,117,tn) = sums_l(nzb,117,tn) +                                    &
-                                           surf_usm_h%ssws(m) * rmask(j,i,sr)            ! w"s"
-                   ENDIF
-#endif
-
-                ENDIF
-
+                ENDDO
              ENDIF
 
@@ -1107,11 +1109,11 @@
                    sums_l(nzb,112,tn) = sums_l(nzb,112,tn) + surf_def_h(0)%ol(m) * rmask(j,i,sr) ! L
                 ENDIF
-                IF ( surf_lsm_h%end_index(j,i) >= surf_lsm_h%start_index(j,i) )  THEN
-                   m = surf_lsm_h%start_index(j,i)
-                   sums_l(nzb,112,tn) = sums_l(nzb,112,tn) + surf_lsm_h%ol(m)    * rmask(j,i,sr) ! L
-                ENDIF
-                IF ( surf_usm_h%end_index(j,i) >= surf_usm_h%start_index(j,i) )  THEN
-                   m = surf_usm_h%start_index(j,i)
-                   sums_l(nzb,112,tn) = sums_l(nzb,112,tn) + surf_usm_h%ol(m)    * rmask(j,i,sr) ! L
+                IF ( surf_lsm_h(0)%end_index(j,i) >= surf_lsm_h(0)%start_index(j,i) )  THEN
+                   m = surf_lsm_h(0)%start_index(j,i)
+                   sums_l(nzb,112,tn) = sums_l(nzb,112,tn) + surf_lsm_h(0)%ol(m)    * rmask(j,i,sr) ! L
+                ENDIF
+                IF ( surf_usm_h(0)%end_index(j,i) >= surf_usm_h(0)%start_index(j,i) )  THEN
+                   m = surf_usm_h(0)%start_index(j,i)
+                   sums_l(nzb,112,tn) = sums_l(nzb,112,tn) + surf_usm_h(0)%ol(m)    * rmask(j,i,sr) ! L
                 ENDIF
              ENDIF
@@ -1131,29 +1133,29 @@
                                         surf_def_h(0)%rad_sw_out(m) * rmask(j,i,sr)
                 ENDIF
-                IF ( surf_lsm_h%end_index(j,i) >= surf_lsm_h%start_index(j,i) )  THEN
-                   m = surf_lsm_h%start_index(j,i)
+                IF ( surf_lsm_h(0)%end_index(j,i) >= surf_lsm_h(0)%start_index(j,i) )  THEN
+                   m = surf_lsm_h(0)%start_index(j,i)
                    sums_l(nzb,99,tn)  = sums_l(nzb,99,tn)   +                                      &
-                                        surf_lsm_h%rad_net(m)    * rmask(j,i,sr)
+                                        surf_lsm_h(0)%rad_net(m)    * rmask(j,i,sr)
                    sums_l(nzb,100,tn) = sums_l(nzb,100,tn)  +                                      &
-                                        surf_lsm_h%rad_lw_in(m)  * rmask(j,i,sr)
+                                        surf_lsm_h(0)%rad_lw_in(m)  * rmask(j,i,sr)
                    sums_l(nzb,101,tn) = sums_l(nzb,101,tn)  +                                      &
-                                        surf_lsm_h%rad_lw_out(m) * rmask(j,i,sr)
+                                        surf_lsm_h(0)%rad_lw_out(m) * rmask(j,i,sr)
                    sums_l(nzb,102,tn) = sums_l(nzb,102,tn)  +                                      &
-                                        surf_lsm_h%rad_sw_in(m)  * rmask(j,i,sr)
+                                        surf_lsm_h(0)%rad_sw_in(m)  * rmask(j,i,sr)
                    sums_l(nzb,103,tn) = sums_l(nzb,103,tn)  +                                      &
-                                        surf_lsm_h%rad_sw_out(m) * rmask(j,i,sr)
-                ENDIF
-                IF ( surf_usm_h%end_index(j,i) >= surf_usm_h%start_index(j,i) )  THEN
-                   m = surf_usm_h%start_index(j,i)
+                                        surf_lsm_h(0)%rad_sw_out(m) * rmask(j,i,sr)
+                ENDIF
+                IF ( surf_usm_h(0)%end_index(j,i) >= surf_usm_h(0)%start_index(j,i) )  THEN
+                   m = surf_usm_h(0)%start_index(j,i)
                    sums_l(nzb,99,tn)  = sums_l(nzb,99,tn)   +                                      &
-                                        surf_usm_h%rad_net(m)    * rmask(j,i,sr)
+                                        surf_usm_h(0)%rad_net(m)    * rmask(j,i,sr)
                    sums_l(nzb,100,tn) = sums_l(nzb,100,tn)  +                                      &
-                                        surf_usm_h%rad_lw_in(m)  * rmask(j,i,sr)
+                                        surf_usm_h(0)%rad_lw_in(m)  * rmask(j,i,sr)
                    sums_l(nzb,101,tn) = sums_l(nzb,101,tn)  +                                      &
-                                        surf_usm_h%rad_lw_out(m) * rmask(j,i,sr)
+                                        surf_usm_h(0)%rad_lw_out(m) * rmask(j,i,sr)
                    sums_l(nzb,102,tn) = sums_l(nzb,102,tn)  +                                      &
-                                        surf_usm_h%rad_sw_in(m)  * rmask(j,i,sr)
+                                        surf_usm_h(0)%rad_sw_in(m)  * rmask(j,i,sr)
                    sums_l(nzb,103,tn) = sums_l(nzb,103,tn)  +                                      &
-                                        surf_usm_h%rad_sw_out(m) * rmask(j,i,sr)
+                                        surf_usm_h(0)%rad_sw_out(m) * rmask(j,i,sr)
                 ENDIF
 
@@ -1172,33 +1174,33 @@
                                            surf_def_h(0)%rrtm_asdir(m,0) * rmask(j,i,sr)
                    ENDIF
-                   IF ( surf_lsm_h%end_index(j,i) >= surf_lsm_h%start_index(j,i) )  THEN
-                      m = surf_lsm_h%start_index(j,i)
+                   IF ( surf_lsm_h(0)%end_index(j,i) >= surf_lsm_h(0)%start_index(j,i) )  THEN
+                      m = surf_lsm_h(0)%start_index(j,i)
                       sums_l(nzb,108,tn)  = sums_l(nzb,108,tn)  +                                  &
-                                            SUM( surf_lsm_h%frac(m,:) *                            &
-                                                 surf_lsm_h%rrtm_aldif(m,:) ) * rmask(j,i,sr)
+                                            SUM( surf_lsm_h(0)%frac(m,:) *                            &
+                                                 surf_lsm_h(0)%rrtm_aldif(m,:) ) * rmask(j,i,sr)
                       sums_l(nzb,109,tn) = sums_l(nzb,109,tn)  +                                   &
-                                           SUM( surf_lsm_h%frac(m,:) *                             &
-                                                surf_lsm_h%rrtm_aldir(m,:) ) * rmask(j,i,sr)
+                                           SUM( surf_lsm_h(0)%frac(m,:) *                             &
+                                                surf_lsm_h(0)%rrtm_aldir(m,:) ) * rmask(j,i,sr)
                       sums_l(nzb,110,tn) = sums_l(nzb,110,tn)  +                                   &
-                                           SUM( surf_lsm_h%frac(m,:) *                             &
-                                                surf_lsm_h%rrtm_asdif(m,:) ) * rmask(j,i,sr)
+                                           SUM( surf_lsm_h(0)%frac(m,:) *                             &
+                                                surf_lsm_h(0)%rrtm_asdif(m,:) ) * rmask(j,i,sr)
                       sums_l(nzb,111,tn) = sums_l(nzb,111,tn)  +                                   &
-                                           SUM( surf_lsm_h%frac(m,:) *                             &
-                                                surf_lsm_h%rrtm_asdir(m,:) ) * rmask(j,i,sr)
+                                           SUM( surf_lsm_h(0)%frac(m,:) *                             &
+                                                surf_lsm_h(0)%rrtm_asdir(m,:) ) * rmask(j,i,sr)
                    ENDIF
-                   IF ( surf_usm_h%end_index(j,i) >= surf_usm_h%start_index(j,i) )  THEN
-                      m = surf_usm_h%start_index(j,i)
+                   IF ( surf_usm_h(0)%end_index(j,i) >= surf_usm_h(0)%start_index(j,i) )  THEN
+                      m = surf_usm_h(0)%start_index(j,i)
                       sums_l(nzb,108,tn)  = sums_l(nzb,108,tn)  +                                  &
-                                            SUM( surf_usm_h%frac(m,:) *                            &
-                                                 surf_usm_h%rrtm_aldif(m,:) ) * rmask(j,i,sr)
+                                            SUM( surf_usm_h(0)%frac(m,:) *                            &
+                                                 surf_usm_h(0)%rrtm_aldif(m,:) ) * rmask(j,i,sr)
                       sums_l(nzb,109,tn) = sums_l(nzb,109,tn)  +                                   &
-                                           SUM( surf_usm_h%frac(m,:) *                             &
-                                                surf_usm_h%rrtm_aldir(m,:) ) * rmask(j,i,sr)
+                                           SUM( surf_usm_h(0)%frac(m,:) *                             &
+                                                surf_usm_h(0)%rrtm_aldir(m,:) ) * rmask(j,i,sr)
                       sums_l(nzb,110,tn) = sums_l(nzb,110,tn)  +                                   &
-                                           SUM( surf_usm_h%frac(m,:) *                             &
-                                                surf_usm_h%rrtm_asdif(m,:) ) * rmask(j,i,sr)
+                                           SUM( surf_usm_h(0)%frac(m,:) *                             &
+                                                surf_usm_h(0)%rrtm_asdif(m,:) ) * rmask(j,i,sr)
                       sums_l(nzb,111,tn) = sums_l(nzb,111,tn)  +                                   &
-                                           SUM( surf_usm_h%frac(m,:) *                             &
-                                                surf_usm_h%rrtm_asdir(m,:) ) * rmask(j,i,sr)
+                                           SUM( surf_usm_h(0)%frac(m,:) *                             &
+                                                surf_usm_h(0)%rrtm_asdir(m,:) ) * rmask(j,i,sr)
                    ENDIF
 
@@ -1405,15 +1407,15 @@
           !$ tn = omp_get_thread_num()
           !$OMP DO
-          DO  m = 1, surf_lsm_h%ns
-             i = surf_lsm_h%i(m)
-             j = surf_lsm_h%j(m)
+          DO  m = 1, surf_lsm_h(0)%ns
+             i = surf_lsm_h(0)%i(m)
+             j = surf_lsm_h(0)%j(m)
 
              IF ( i >= nxl  .AND.  i <= nxr  .AND.  j >= nys  .AND.  j <= nyn )  THEN
-                sums_l(nzb,93,tn)  = sums_l(nzb,93,tn) + surf_lsm_h%ghf(m)       * rmask(j,i,sr)
-                sums_l(nzb,94,tn)  = sums_l(nzb,94,tn) + surf_lsm_h%qsws_liq(m)  * rmask(j,i,sr)
-                sums_l(nzb,95,tn)  = sums_l(nzb,95,tn) + surf_lsm_h%qsws_soil(m) * rmask(j,i,sr)
-                sums_l(nzb,96,tn)  = sums_l(nzb,96,tn) + surf_lsm_h%qsws_veg(m)  * rmask(j,i,sr)
-                sums_l(nzb,97,tn)  = sums_l(nzb,97,tn) + surf_lsm_h%r_a(m)       * rmask(j,i,sr)
-                sums_l(nzb,98,tn)  = sums_l(nzb,98,tn) + surf_lsm_h%r_s(m)       * rmask(j,i,sr)
+                sums_l(nzb,93,tn)  = sums_l(nzb,93,tn) + surf_lsm_h(0)%ghf(m)       * rmask(j,i,sr)
+                sums_l(nzb,94,tn)  = sums_l(nzb,94,tn) + surf_lsm_h(0)%qsws_liq(m)  * rmask(j,i,sr)
+                sums_l(nzb,95,tn)  = sums_l(nzb,95,tn) + surf_lsm_h(0)%qsws_soil(m) * rmask(j,i,sr)
+                sums_l(nzb,96,tn)  = sums_l(nzb,96,tn) + surf_lsm_h(0)%qsws_veg(m)  * rmask(j,i,sr)
+                sums_l(nzb,97,tn)  = sums_l(nzb,97,tn) + surf_lsm_h(0)%r_a(m)       * rmask(j,i,sr)
+                sums_l(nzb,98,tn)  = sums_l(nzb,98,tn) + surf_lsm_h(0)%r_s(m)       * rmask(j,i,sr)
              ENDIF
           ENDDO
@@ -1424,14 +1426,14 @@
           !$ tn = omp_get_thread_num()
           !$OMP DO
-          DO  m = 1, surf_lsm_h%ns
-
-             i = surf_lsm_h%i(m)
-             j = surf_lsm_h%j(m)
+          DO  m = 1, surf_lsm_h(0)%ns
+
+             i = surf_lsm_h(0)%i(m)
+             j = surf_lsm_h(0)%j(m)
 
              IF ( i >= nxl  .AND.  i <= nxr  .AND.  j >= nys  .AND.  j <= nyn )  THEN
 
                 DO  k = nzb_soil, nzt_soil
-                   sums_l(k,89,tn)  = sums_l(k,89,tn)  + t_soil_h%var_2d(k,m) * rmask(j,i,sr)
-                   sums_l(k,91,tn)  = sums_l(k,91,tn)  + m_soil_h%var_2d(k,m) * rmask(j,i,sr)
+                   sums_l(k,89,tn)  = sums_l(k,89,tn)  + t_soil_h(0)%var_2d(k,m) * rmask(j,i,sr)
+                   sums_l(k,91,tn)  = sums_l(k,91,tn)  + m_soil_h(0)%var_2d(k,m) * rmask(j,i,sr)
                 ENDDO
              ENDIF
