Index: palm/trunk/SOURCE/data_output_2d.f90
===================================================================
--- palm/trunk/SOURCE/data_output_2d.f90	(revision 1360)
+++ palm/trunk/SOURCE/data_output_2d.f90	(revision 1551)
@@ -20,5 +20,8 @@
 ! Current revisions:
 ! -----------------
-! 
+! Added suppport for land surface model and radiation model output. In the course
+! of this action, the limits for vertical loops have been changed (from nzb and 
+! nzt+1 to nzb_do and nzt_do, respectively in order to allow soil model output).
+! Moreover, a new vertical grid zs was introduced.
 ! 
 ! Former revisions:
@@ -127,5 +130,13 @@
                
     USE kinds
-        
+    
+    USE land_surface_model_mod,                                                &
+        ONLY:  c_liq, c_liq_av, c_soil_av, c_veg, c_veg_av, ghf_eb,            &
+               ghf_eb_av, lai, lai_av, m_liq_eb, m_liq_eb_av, m_soil,          &
+               m_soil_av, nzb_soil, nzt_soil, qsws_eb, qsws_eb_av,             &
+               qsws_liq_eb, qsws_liq_eb_av, qsws_soil_eb, qsws_soil_eb_av,     &
+               qsws_veg_eb, qsws_veg_eb_av, shf_eb, shf_eb_av, t_soil,         &
+               t_soil_av, zs
+    
     USE netcdf_control
 
@@ -135,4 +146,7 @@
     
     USE pegrid
+
+    USE radiation_model_mod,                                                   &
+        ONLY:  rad_net, rad_net_av, rad_sw_in, rad_sw_in_av
 
     IMPLICIT NONE
@@ -157,4 +171,6 @@
     INTEGER(iwp) ::  n         !: 
     INTEGER(iwp) ::  ns        !: 
+    INTEGER(iwp) ::  nzb_do    !: lower limit of the data field (usually nzb)
+    INTEGER(iwp) ::  nzt_do    !: upper limit of the data field (usually nzt+1)
     INTEGER(iwp) ::  psi       !: 
     INTEGER(iwp) ::  s         !: 
@@ -343,4 +359,7 @@
 
        IF ( do2d_mode == mode )  THEN
+
+          nzb_do = nzb
+          nzt_do = nzt+1
 !
 !--       Store the array chosen on the temporary array.
@@ -356,4 +375,94 @@
                 IF ( mode == 'xy' )  level_z = zu
 
+             CASE ( 'c_liq*_xy' )        ! 2d-array
+                IF ( av == 0 )  THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) = c_liq(j,i) * c_veg(j,i)
+                      ENDDO
+                   ENDDO
+                ELSE
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) = c_liq_av(j,i)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                two_d = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)
+
+             CASE ( 'c_soil*_xy' )        ! 2d-array
+                IF ( av == 0 )  THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) = 1.0_wp - c_veg(j,i)
+                      ENDDO
+                   ENDDO
+                ELSE
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) = c_soil_av(j,i)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                two_d = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)
+
+             CASE ( 'c_veg*_xy' )        ! 2d-array
+                IF ( av == 0 )  THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) = c_veg(j,i)
+                      ENDDO
+                   ENDDO
+                ELSE
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) = c_veg_av(j,i)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                two_d = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)
+
+             CASE ( 'ghf_eb*_xy' )        ! 2d-array
+                IF ( av == 0 )  THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) = ghf_eb(j,i)
+                      ENDDO
+                   ENDDO
+                ELSE
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) = ghf_eb_av(j,i)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                two_d = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)
+
+             CASE ( 'lai*_xy' )        ! 2d-array
+                IF ( av == 0 )  THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) = lai(j,i)
+                      ENDDO
+                   ENDDO
+                ELSE
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) = lai_av(j,i)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                two_d = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)
+
              CASE ( 'lpt_xy', 'lpt_xz', 'lpt_yz' )
                 IF ( av == 0 )  THEN
@@ -382,4 +491,32 @@
                 two_d = .TRUE.
                 level_z(nzb+1) = zu(nzb+1)
+
+             CASE ( 'm_liq_eb*_xy' )        ! 2d-array
+                IF ( av == 0 )  THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) = m_liq_eb(j,i)
+                      ENDDO
+                   ENDDO
+                ELSE
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) = m_liq_eb_av(j,i)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                two_d = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)
+
+             CASE ( 'm_soil_xy', 'm_soil_xz', 'm_soil_yz' )
+                nzb_do = nzb_soil
+                nzt_do = nzt_soil
+                IF ( av == 0 )  THEN
+                   to_be_resorted => m_soil
+                ELSE
+                   to_be_resorted => m_soil_av
+                ENDIF
+                IF ( mode == 'xy' )  level_z = zs
 
              CASE ( 'nr_xy', 'nr_xz', 'nr_yz' )
@@ -665,4 +802,76 @@
                 level_z(nzb+1) = zu(nzb+1)
 
+             CASE ( 'qsws_eb*_xy' )        ! 2d-array
+                IF ( av == 0 ) THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) =  qsws_eb(j,i)
+                      ENDDO
+                   ENDDO
+                ELSE
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng 
+                         local_pf(i,j,nzb+1) =  qsws_eb_av(j,i)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                two_d = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)
+
+             CASE ( 'qsws_liq_eb*_xy' )        ! 2d-array
+                IF ( av == 0 ) THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) =  qsws_liq_eb(j,i)
+                      ENDDO
+                   ENDDO
+                ELSE
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng 
+                         local_pf(i,j,nzb+1) =  qsws_liq_eb_av(j,i)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                two_d = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)
+
+             CASE ( 'qsws_soil_eb*_xy' )        ! 2d-array
+                IF ( av == 0 ) THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) =  qsws_soil_eb(j,i)
+                      ENDDO
+                   ENDDO
+                ELSE
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng 
+                         local_pf(i,j,nzb+1) =  qsws_soil_eb_av(j,i)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                two_d = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)
+
+             CASE ( 'qsws_veg_eb*_xy' )        ! 2d-array
+                IF ( av == 0 ) THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) =  qsws_veg_eb(j,i)
+                      ENDDO
+                   ENDDO
+                ELSE
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng 
+                         local_pf(i,j,nzb+1) =  qsws_veg_eb_av(j,i)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                two_d = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)
+
              CASE ( 'qv_xy', 'qv_xz', 'qv_yz' )
                 IF ( av == 0 )  THEN
@@ -680,4 +889,40 @@
                 IF ( mode == 'xy' )  level_z = zu
 
+             CASE ( 'rad_net*_xy' )        ! 2d-array
+                IF ( av == 0 ) THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) =  rad_net(j,i)
+                      ENDDO
+                   ENDDO
+                ELSE
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng 
+                         local_pf(i,j,nzb+1) =  rad_net_av(j,i)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                two_d = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)
+
+             CASE ( 'rad_sw_in*_xy' )        ! 2d-array
+                IF ( av == 0 ) THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) =  rad_sw_in(j,i)
+                      ENDDO
+                   ENDDO
+                ELSE
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng 
+                         local_pf(i,j,nzb+1) =  rad_sw_in_av(j,i)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                two_d = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)
+
              CASE ( 'rho_xy', 'rho_xz', 'rho_yz' )
                 IF ( av == 0 )  THEN
@@ -719,4 +964,22 @@
                 level_z(nzb+1) = zu(nzb+1)
 
+             CASE ( 'shf_eb*_xy' )        ! 2d-array
+                IF ( av == 0 ) THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) =  shf_eb(j,i)
+                      ENDDO
+                   ENDDO
+                ELSE
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) =  shf_eb_av(j,i)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                two_d = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)
+
              CASE ( 't*_xy' )        ! 2d-array
                 IF ( av == 0 )  THEN
@@ -736,4 +999,14 @@
                 two_d = .TRUE.
                 level_z(nzb+1) = zu(nzb+1)
+
+             CASE ( 't_soil_xy', 't_soil_xz', 't_soil_yz' )
+                nzb_do = nzb_soil
+                nzt_do = nzt_soil
+                IF ( av == 0 )  THEN
+                   to_be_resorted => t_soil
+                ELSE
+                   to_be_resorted => t_soil_av
+                ENDIF
+                IF ( mode == 'xy' )  level_z = zs
 
              CASE ( 'u_xy', 'u_xz', 'u_yz' )
@@ -839,5 +1112,5 @@
 !--             User defined quantity
                 CALL user_data_output_2d( av, do2d(av,if), found, grid,        &
-                                          local_pf, two_d )
+                                          local_pf, two_d, nzb_do, nzt_do )
                 resorted = .TRUE.
 
@@ -848,4 +1121,6 @@
                 ELSEIF ( grid == 'zu1' ) THEN
                    IF ( mode == 'xy' )  level_z(nzb+1) = zu(nzb+1)
+                ELSEIF ( grid == 'zs' ) THEN
+                   IF ( mode == 'xy' )  level_z = zs
                 ENDIF
 
@@ -863,5 +1138,5 @@
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = to_be_resorted(k,j,i)
                    ENDDO
@@ -874,5 +1149,5 @@
 !--       section mode chosen.
           is = 1
-   loop1: DO  WHILE ( section(is,s) /= -9999  .OR.  two_d )
+   loop1: DO WHILE ( section(is,s) /= -9999  .OR.  two_d )
 
              SELECT CASE ( mode )
@@ -885,4 +1160,11 @@
                    ELSE
                       layer_xy = section(is,s)
+                   ENDIF
+
+!
+!--                Exit the loop for layers beyond the data output domain
+!--                (used for soil model)
+                   IF ( layer_xy .GT. nzt_do )  THEN 
+                      EXIT loop1
                    ENDIF
 
@@ -916,5 +1198,5 @@
 !
 !--                   Carry out the averaging (all data are on the PE)
-                      DO  k = nzb, nzt+1
+                      DO  k = nzb_do, nzt_do
                          DO  j = nysg, nyng
                             DO  i = nxlg, nxrg
@@ -924,5 +1206,5 @@
                       ENDDO
 
-                      local_2d = local_2d / ( nzt -nzb + 2.0_wp)
+                      local_2d = local_2d / ( nzt_do - nzb_do + 1.0_wp)
 
                    ELSE
@@ -967,5 +1249,5 @@
                          DO  i = 0, io_blocks-1
                             IF ( i == io_group )  THEN
-                               WRITE ( 21 )  nxlg, nxrg, nysg, nyng
+                               WRITE ( 21 )  nxlg, nxrg, nysg, nyng, nysg, nyng
                                WRITE ( 21 )  local_2d
                             ENDIF
@@ -1103,10 +1385,10 @@
                    IF ( section(is,s) == -1 )  THEN
 
-                      ALLOCATE( local_2d_l(nxlg:nxrg,nzb:nzt+1) )
+                      ALLOCATE( local_2d_l(nxlg:nxrg,nzb_do:nzt_do) )
                       local_2d_l = 0.0_wp
-                      ngp = ( nxrg-nxlg+1 ) * ( nzt-nzb+2 )
+                      ngp = ( nxrg-nxlg + 1 ) * ( nzt_do-nzb_do + 1 )
 !
 !--                   First local averaging on the PE
-                      DO  k = nzb, nzt+1
+                      DO  k = nzb_do, nzt_do
                          DO  j = nys, nyn
                             DO  i = nxlg, nxrg
@@ -1120,6 +1402,6 @@
 !--                   Now do the averaging over all PEs along y
                       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-                      CALL MPI_ALLREDUCE( local_2d_l(nxlg,nzb),                &
-                                          local_2d(nxlg,nzb), ngp, MPI_REAL,   &
+                      CALL MPI_ALLREDUCE( local_2d_l(nxlg,nzb_do),                &
+                                          local_2d(nxlg,nzb_do), ngp, MPI_REAL,   &
                                           MPI_SUM, comm1dy, ierr )
 #else
@@ -1136,5 +1418,5 @@
                       IF ( section(is,s) >= nys  .AND.  section(is,s) <= nyn ) &
                       THEN
-                         local_2d = local_pf(:,section(is,s),nzb:nzt+1)
+                         local_2d = local_pf(:,section(is,s),nzb_do:nzt_do)
                       ENDIF
 
@@ -1157,5 +1439,5 @@
 !--                      output file afterwards to increase the performance.
                          DO  i = nxlg, nxrg
-                            DO  k = nzb, nzt+1
+                            DO  k = nzb_do, nzt_do
                                local_2d_sections_l(i,is,k) = local_2d(i,k)
                             ENDDO
@@ -1184,8 +1466,8 @@
                                       nys-1 == -1 ) )                          &
                                THEN
-                                  WRITE (22)  nxlg, nxrg, nzb, nzt+1
+                                  WRITE (22)  nxlg, nxrg, nzb_do, nzt_do, nzb, nzt+1
                                   WRITE (22)  local_2d
                                ELSE
-                                  WRITE (22)  -1, -1, -1, -1
+                                  WRITE (22)  -1, -1, -1, -1, -1, -1
                                ENDIF
                             ENDIF
@@ -1203,5 +1485,5 @@
                          CALL MPI_BARRIER( comm2d, ierr )
 
-                         ngp = ( nxrg-nxlg+1 ) * ( nzt-nzb+2 )
+                         ngp = ( nxrg-nxlg + 1 ) * ( nzt_do-nzb_do + 1 )
                          IF ( myid == 0 )  THEN
 !
@@ -1211,5 +1493,5 @@
                                  ( section(is,s) == -1  .AND.  nys-1 == -1 ) ) &
                             THEN
-                               total_2d(nxlg:nxrg,nzb:nzt+1) = local_2d
+                               total_2d(nxlg:nxrg,nzb_do:nzt_do) = local_2d
                             ENDIF
 !
@@ -1240,12 +1522,12 @@
 !--                         Relocate the local array for the next loop increment
                             DEALLOCATE( local_2d )
-                            ALLOCATE( local_2d(nxlg:nxrg,nzb:nzt+1) )
+                            ALLOCATE( local_2d(nxlg:nxrg,nzb_do:nzt_do) )
 
 #if defined( __netcdf )
                             nc_stat = NF90_PUT_VAR( id_set_xz(av),          &
                                                  id_var_do2d(av,if),        &
-                                                 total_2d(0:nx+1,nzb:nzt+1),&
+                                                 total_2d(0:nx+1,nzb_do:nzt_do),&
                             start = (/ 1, is, 1, do2d_xz_time_count(av) /), &
-                                             count = (/ nx+2, 1, nz+2, 1 /) )
+                                             count = (/ nx+2, 1, nzt_do-nzb_do+1, 1 /) )
                             CALL handle_netcdf_error( 'data_output_2d', 58 )
 #endif
@@ -1260,5 +1542,5 @@
                             THEN
                                ind(1) = nxlg; ind(2) = nxrg
-                               ind(3) = nzb;   ind(4) = nzt+1
+                               ind(3) = nzb_do;   ind(4) = nzt_do
                             ELSE
                                ind(1) = -9999; ind(2) = -9999
@@ -1270,5 +1552,5 @@
 !--                         If applicable, send data to PE0.
                             IF ( ind(1) /= -9999 )  THEN
-                               CALL MPI_SEND( local_2d(nxlg,nzb), ngp,         &
+                               CALL MPI_SEND( local_2d(nxlg,nzb_do), ngp,         &
                                               MPI_REAL, 0, 1, comm2d, ierr )
                             ENDIF
@@ -1286,7 +1568,7 @@
                    nc_stat = NF90_PUT_VAR( id_set_xz(av),                   &
                                            id_var_do2d(av,if),              &
-                                           local_2d(nxl:nxr+1,nzb:nzt+1),   &
+                                           local_2d(nxl:nxr+1,nzb_do:nzt_do),   &
                             start = (/ 1, is, 1, do2d_xz_time_count(av) /), &
-                                           count = (/ nx+2, 1, nz+2, 1 /) )
+                                           count = (/ nx+2, 1, nzt_do-nzb_do+1, 1 /) )
                    CALL handle_netcdf_error( 'data_output_2d', 451 )
 #endif
@@ -1322,10 +1604,10 @@
                    IF ( section(is,s) == -1 )  THEN
 
-                      ALLOCATE( local_2d_l(nysg:nyng,nzb:nzt+1) )
+                      ALLOCATE( local_2d_l(nysg:nyng,nzb_do:nzt_do) )
                       local_2d_l = 0.0_wp
-                      ngp = ( nyng-nysg+1 ) * ( nzt-nzb+2 )
+                      ngp = ( nyng-nysg+1 ) * ( nzt_do-nzb_do+1 )
 !
 !--                   First local averaging on the PE
-                      DO  k = nzb, nzt+1
+                      DO  k = nzb_do, nzt_do
                          DO  j = nysg, nyng
                             DO  i = nxl, nxr
@@ -1339,6 +1621,6 @@
 !--                   Now do the averaging over all PEs along x
                       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-                      CALL MPI_ALLREDUCE( local_2d_l(nysg,nzb),                &
-                                          local_2d(nysg,nzb), ngp, MPI_REAL,   &
+                      CALL MPI_ALLREDUCE( local_2d_l(nysg,nzb_do),                &
+                                          local_2d(nysg,nzb_do), ngp, MPI_REAL,   &
                                           MPI_SUM, comm1dx, ierr )
 #else
@@ -1355,5 +1637,5 @@
                       IF ( section(is,s) >= nxl  .AND.  section(is,s) <= nxr ) &
                       THEN
-                         local_2d = local_pf(section(is,s),:,nzb:nzt+1)
+                         local_2d = local_pf(section(is,s),:,nzb_do:nzt_do)
                       ENDIF
 
@@ -1376,5 +1658,5 @@
 !--                      output file afterwards to increase the performance.
                          DO  j = nysg, nyng
-                            DO  k = nzb, nzt+1
+                            DO  k = nzb_do, nzt_do
                                local_2d_sections_l(is,j,k) = local_2d(j,k)
                             ENDDO
@@ -1403,8 +1685,8 @@
                                       nxl-1 == -1 ) )                          &
                                THEN
-                                  WRITE (23)  nysg, nyng, nzb, nzt+1
+                                  WRITE (23)  nysg, nyng, nzb_do, nzt_do, nzb, nzt+1
                                   WRITE (23)  local_2d
                                ELSE
-                                  WRITE (23)  -1, -1, -1, -1
+                                  WRITE (23)  -1, -1, -1, -1, -1, -1
                                ENDIF
                             ENDIF
@@ -1422,5 +1704,5 @@
                          CALL MPI_BARRIER( comm2d, ierr )
 
-                         ngp = ( nyng-nysg+1 ) * ( nzt-nzb+2 )
+                         ngp = ( nyng-nysg+1 ) * ( nzt_do-nzb_do+1 )
                          IF ( myid == 0 )  THEN
 !
@@ -1430,5 +1712,5 @@
                                  ( section(is,s) == -1  .AND.  nxl-1 == -1 ) ) &
                             THEN
-                               total_2d(nysg:nyng,nzb:nzt+1) = local_2d
+                               total_2d(nysg:nyng,nzb_do:nzt_do) = local_2d
                             ENDIF
 !
@@ -1459,12 +1741,12 @@
 !--                         Relocate the local array for the next loop increment
                             DEALLOCATE( local_2d )
-                            ALLOCATE( local_2d(nysg:nyng,nzb:nzt+1) )
+                            ALLOCATE( local_2d(nysg:nyng,nzb_do:nzt_do) )
 
 #if defined( __netcdf )
                             nc_stat = NF90_PUT_VAR( id_set_yz(av),          &
                                                  id_var_do2d(av,if),        &
-                                                 total_2d(0:ny+1,nzb:nzt+1),&
+                                                 total_2d(0:ny+1,nzb_do:nzt_do),&
                             start = (/ is, 1, 1, do2d_yz_time_count(av) /), &
-                                             count = (/ 1, ny+2, nz+2, 1 /) )
+                                             count = (/ 1, ny+2, nzt_do-nzb_do+1, 1 /) )
                             CALL handle_netcdf_error( 'data_output_2d', 61 )
 #endif
@@ -1479,5 +1761,5 @@
                             THEN
                                ind(1) = nysg; ind(2) = nyng
-                               ind(3) = nzb;   ind(4) = nzt+1
+                               ind(3) = nzb_do;   ind(4) = nzt_do
                             ELSE
                                ind(1) = -9999; ind(2) = -9999
@@ -1489,5 +1771,5 @@
 !--                         If applicable, send data to PE0.
                             IF ( ind(1) /= -9999 )  THEN
-                               CALL MPI_SEND( local_2d(nysg,nzb), ngp,         &
+                               CALL MPI_SEND( local_2d(nysg,nzb_do), ngp,         &
                                               MPI_REAL, 0, 1, comm2d, ierr )
                             ENDIF
@@ -1505,7 +1787,7 @@
                    nc_stat = NF90_PUT_VAR( id_set_yz(av),                   &
                                            id_var_do2d(av,if),              &
-                                           local_2d(nys:nyn+1,nzb:nzt+1),   &
+                                           local_2d(nys:nyn+1,nzb_do:nzt_do),   &
                             start = (/ is, 1, 1, do2d_xz_time_count(av) /), &
-                                           count = (/ 1, ny+2, nz+2, 1 /) )
+                                           count = (/ 1, ny+2, nzt_do-nzb_do+1, 1 /) )
                    CALL handle_netcdf_error( 'data_output_2d', 452 )
 #endif
@@ -1595,8 +1877,8 @@
 !
 !--                      Distribute data over all PEs along y
-                         ngp = ( nxrg-nxlg+1 ) * ( nzt-nzb+2 ) * ns
+                         ngp = ( nxrg-nxlg+1 ) * ( nzt_do-nzb_do+1 ) * ns
                          IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
-                         CALL MPI_ALLREDUCE( local_2d_sections_l(nxlg,1,nzb),  &
-                                             local_2d_sections(nxlg,1,nzb),    &
+                         CALL MPI_ALLREDUCE( local_2d_sections_l(nxlg,1,nzb_do),  &
+                                             local_2d_sections(nxlg,1,nzb_do),    &
                                              ngp, MPI_REAL, MPI_SUM, comm1dy,  &
                                              ierr )
@@ -1612,8 +1894,8 @@
                                              id_var_do2d(av,if),               & 
                                              local_2d_sections(nxl:nxr+1,1:ns, &
-                                                nzb:nzt+1),                    &
+                                                nzb_do:nzt_do),                &
                                              start = (/ nxl+1, 1, 1,           &
                                                 do2d_xz_time_count(av) /),     &
-                                             count = (/ nxr-nxl+2, ns, nzt+2,  &
+                                             count = (/ nxr-nxl+2, ns, nzt_do-nzb_do+1,  &
                                                         1 /) )
                       ELSE
@@ -1621,8 +1903,8 @@
                                              id_var_do2d(av,if),               &
                                              local_2d_sections(nxl:nxr,1:ns,   &
-                                                nzb:nzt+1),                    &
+                                                nzb_do:nzt_do),                &
                                              start = (/ nxl+1, 1, 1,           &
                                                 do2d_xz_time_count(av) /),     &
-                                             count = (/ nxr-nxl+1, ns, nzt+2,  &
+                                             count = (/ nxr-nxl+1, ns, nzt_do-nzb_do+1,  &
                                                 1 /) )
                       ENDIF
@@ -1647,6 +1929,6 @@
                          ngp = ( nyng-nysg+1 ) * ( nzt-nzb + 2 ) * ns
                          IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
-                         CALL MPI_ALLREDUCE( local_2d_sections_l(1,nysg,nzb),  &
-                                             local_2d_sections(1,nysg,nzb),    &
+                         CALL MPI_ALLREDUCE( local_2d_sections_l(1,nysg,nzb_do),  &
+                                             local_2d_sections(1,nysg,nzb_do),    &
                                              ngp, MPI_REAL, MPI_SUM, comm1dx,  &
                                              ierr )
@@ -1662,18 +1944,18 @@
                                              id_var_do2d(av,if),               &
                                              local_2d_sections(1:ns,           &
-                                                nys:nyn+1,nzb:nzt+1),          &
+                                                nys:nyn+1,nzb_do:nzt_do),      &
                                              start = (/ 1, nys+1, 1,           &
                                                 do2d_yz_time_count(av) /),     &
                                              count = (/ ns, nyn-nys+2,         &
-                                                        nzt+2, 1 /) )
+                                                        nzt_do-nzb_do+1, 1 /) )
                       ELSE
                          nc_stat = NF90_PUT_VAR( id_set_yz(av),                &
                                              id_var_do2d(av,if),               &
                                              local_2d_sections(1:ns,nys:nyn,   &
-                                                nzb:nzt+1),                    &
+                                                nzb_do:nzt_do),                &
                                              start = (/ 1, nys+1, 1,           &
                                                 do2d_yz_time_count(av) /),     &
                                              count = (/ ns, nyn-nys+1,         &
-                                                        nzt+2, 1 /) )
+                                                        nzt_do-nzb_do+1, 1 /) )
                       ENDIF
 
