Index: palm/trunk/UTIL/combine_plot_fields.f90
===================================================================
--- palm/trunk/UTIL/combine_plot_fields.f90	(revision 72)
+++ palm/trunk/UTIL/combine_plot_fields.f90	(revision 108)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! 
+! loop for processing of output by coupled runs, id_string does not contain 
+! modus any longer
 !
 ! Former revisions:
@@ -59,4 +60,8 @@
 !    -Mpreprocess -D__netcdf -I /home/usr5/mkanda/netcdf-3.6.1/src/f90 -r8
 !    -fast -L/home/usr5/mkanda/netcdf-3.6.1/src/libsrc -lnetcdf
+! FIMM:
+!    ifort combine...f90 -o combine...x
+!    -axW -cpp -openmp -r8 -nbs -convert little_endian -D__netcdf 
+!    -I /local/netcdf/include -Vaxlib -L/local/netcdf/lib -lnetcdf
 !------------------------------------------------------------------------------!
 
@@ -69,6 +74,6 @@
 !
 !-- Local variables
-    CHARACTER (LEN=2)    ::  modus
-    CHARACTER (LEN=7)    ::  id_char
+    CHARACTER (LEN=2)    ::  modus, model_string
+    CHARACTER (LEN=4)    ::  id_string
     CHARACTER (LEN=10)   ::  dimname, var_name
     CHARACTER (LEN=40)   ::  filename
@@ -79,9 +84,10 @@
 
     INTEGER ::  av, danz, i, id,             &
-                j, nc_stat, nxa, nxag, nxe, nxeg, nya,   &
+                j, model, models, nc_stat, nxa, nxag, nxe, nxeg, nya,   &
                 nyag, nye, nyeg, nza, nzag, nze, nzeg, pos, time_step, xa, xe, &
                 ya, ye, za, ze
 
-    INTEGER, DIMENSION(0:1) ::  current_level, current_var, fanz, id_set, id_var_time, num_var
+    INTEGER, DIMENSION(0:1) ::  current_level, current_var, fanz, id_set, &
+         id_var_time, num_var
 
     INTEGER, DIMENSION(4) ::  id_dims_loc
@@ -99,40 +105,369 @@
     REAL(spk), DIMENSION(:,:,:), ALLOCATABLE ::  pf3d
 
-
-!
-!-- 2D-arrays for ISO2D
-!-- Main loop for the three different cross-sections, starting with xy-section
-    modus = 'XY'
     PRINT*, ''
-    DO  WHILE ( modus == 'XY'  .OR.  modus == 'XZ'  .OR.  modus == 'YZ' )
-!
-!--    Check, if file from PE0 exists. If it does not exist, PALM did not
-!--    create any output for this cross-section.
-       danz = 0
-       WRITE (id_char,'(A2,''_'',I4.4)')  modus, danz
-       INQUIRE ( FILE='PLOT2D_'//id_char, EXIST=found )
-!
-!--    Find out the number of files (equal to the number of PEs which
-!--    have been used in PALM) and open them
-       DO  WHILE ( found )
-
-          OPEN ( danz+110, FILE='PLOT2D_'//id_char, FORM='UNFORMATTED' )
-          danz = danz + 1
-          WRITE (id_char,'(A2,''_'',I4.4)')  modus, danz
-          INQUIRE ( FILE='PLOT2D_'//id_char, EXIST=found )
+    PRINT*, '*** combine_plot_fields ***'
+    INQUIRE( FILE='COUPLING_PORT_OPENED', EXIST=found )
+    IF ( found )  THEN
+       models = 2
+       PRINT*, '    coupled run'
+    ELSE
+       models = 1
+       PRINT*, '    uncoupled run'
+    ENDIF
+    DO model = 1, models
+       IF ( models == 2 )  THEN
+          PRINT*, ''
+          PRINT*, ''
+          PRINT*, '*** combine_plot_fields ***'
+          IF ( model == 2 )  THEN
+             model_string = "_O"
+             PRINT*, '    now combining ocean data'
+             PRINT*, '    ========================'
+          ELSE
+             model_string = ""
+             PRINT*, '    now combining atmosphere data'
+             PRINT*, '    ============================='
+          ENDIF
+       ENDIF
+!
+!--    2D-arrays for ISO2D
+!--    Main loop for the three different cross-sections, starting with 
+!--    xy-section
+       modus = 'XY'
+       PRINT*, ''
+       DO  WHILE ( modus == 'XY'  .OR.  modus == 'XZ'  .OR.  modus == 'YZ' )
+!
+!--       Check, if file from PE0 exists. If it does not exist, PALM did not
+!--       create any output for this cross-section.
+          danz = 0
+          WRITE (id_string,'(I4.4)')  danz
+          INQUIRE ( &
+               FILE='PLOT2D_'//modus//TRIM( model_string )//'_'//id_string, &
+               EXIST=found )
+!
+!--       Find out the number of files (equal to the number of PEs which
+!--       have been used in PALM) and open them
+          DO  WHILE ( found )
+
+             OPEN ( danz+110, &
+                  FILE='PLOT2D_'//modus//TRIM( model_string )//'_'//id_string, &
+                  FORM='UNFORMATTED' )
+             danz = danz + 1
+             WRITE (id_string,'(I4.4)')  danz
+             INQUIRE ( &
+                  FILE='PLOT2D_'//modus//TRIM( model_string )//'_'//id_string, &
+                  EXIST=found )
+
+          ENDDO
+
+!
+!--       Inquire whether an iso2d parameter file exists
+          INQUIRE( FILE='PLOT2D_'//modus//'_GLOBAL'//TRIM( model_string ), &
+               EXIST=iso2d_output )
+
+!
+!--       Inquire whether a NetCDF file exists
+          INQUIRE( FILE='DATA_2D_'//modus//'_NETCDF'//TRIM( model_string ), &
+               EXIST=netcdf_0 )
+
+!
+!--       Inquire whether a NetCDF file for time-averaged data exists
+          INQUIRE( FILE='DATA_2D_'//modus//'_AV_NETCDF'//TRIM( model_string ),&
+               EXIST=netcdf_1 )
+
+          IF ( netcdf_0  .OR.  netcdf_1 )  THEN
+             netcdf_output = .TRUE.
+          ELSE
+             netcdf_output = .FALSE.
+          ENDIF
+
+!
+!--       Info-output
+          PRINT*, ''
+          PRINT*, '*** combine_plot_fields ***'
+#if defined( __netcdf )
+          IF ( netcdf_output )  PRINT*, '    NetCDF output enabled'
+#else
+          IF ( netcdf_output )  THEN
+             PRINT*, '--- Sorry, no NetCDF support on this host'
+             netcdf_output = .FALSE.
+          ENDIF
+#endif
+          IF ( danz /= 0 )  THEN
+             PRINT*, '    ',modus,'-section:  ', danz, ' file(s) found'
+          ELSE
+             PRINT*, '    no ', modus, '-section data available'
+          ENDIF
+
+          IF ( netcdf_output  .AND.  danz /= 0 )  THEN
+#if defined( __netcdf )
+             DO  av = 0, 1
+
+                IF ( av == 0  .AND.  .NOT.  netcdf_0 )  CYCLE
+                IF ( av == 1  .AND.  .NOT.  netcdf_1 )  CYCLE
+
+!
+!--             Open NetCDF dataset
+                IF ( av == 0 )  THEN
+                   filename = 'DATA_2D_'//modus//'_NETCDF' &
+                        //TRIM( model_string )
+                ELSE
+                   filename = 'DATA_2D_'//modus//'_AV_NETCDF' &
+                        //TRIM( model_string )
+                ENDIF
+                nc_stat = NF90_OPEN( filename, NF90_WRITE, id_set(av) )
+                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 1 )
+
+!
+!--             Get the list of variables (order of variables corresponds with
+!--             the order of data on the binary file)
+                var_list(av) = ' '    ! GET_ATT does not assign trailing blanks
+                nc_stat = NF90_GET_ATT( id_set(av), NF90_GLOBAL, 'VAR_LIST', &
+                     var_list(av) )
+                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 2 )
+
+!
+!--             Inquire id of the time coordinate variable
+                nc_stat = NF90_INQ_VARID( id_set(av), 'time', id_var_time(av) )
+                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 3 )
+
+!
+!--             Count number of variables; there is one more semicolon in the
+!--             string than variable names
+                num_var(av) = -1
+                DO  i = 1, LEN( var_list(av) )
+                   IF ( var_list(av)(i:i) == ';' )  num_var(av) = num_var(av) +1
+                ENDDO
+
+!
+!--             Extract the variable names from the list and inquire their
+!--             NetCDF IDs
+                pos = INDEX( var_list(av), ';' )
+!
+!--             Loop over all variables
+                DO  i = 1, num_var(av)
+
+!
+!--                Extract variable name from list
+                   var_list(av) = var_list(av)(pos+1:)
+                   pos = INDEX( var_list(av), ';' )
+                   var_name = var_list(av)(1:pos-1)
+
+!
+!--                Get variable ID from name
+                   nc_stat = NF90_INQ_VARID( id_set(av), TRIM( var_name ), &
+                        id_var(av,i) )
+                   IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 4 )
+
+!
+!--                Get number of x/y/z levels for that variable
+                   nc_stat = NF90_INQUIRE_VARIABLE( id_set(av), id_var(av,i), &
+                        dimids = id_dims_loc )
+                   IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 5 )
+                   id_dims(av,:) = id_dims_loc
+
+!
+!--                Inquire dimension ID
+                   DO  j = 1, 4
+                      nc_stat = NF90_INQUIRE_DIMENSION( id_set(av), &
+                           id_dims(av,j), dimname, levels(av,i) )
+                      IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 6 )
+
+                      IF ( modus == 'XY' .AND. INDEX(dimname, 'z') /= 0 )  EXIT
+                      IF ( modus == 'XZ' .AND. INDEX(dimname, 'y') /= 0 )  EXIT
+                      IF ( modus == 'YZ' .AND. INDEX(dimname, 'x') /= 0 )  EXIT
+                   ENDDO
+
+                ENDDO
+
+             ENDDO   ! av = 0, 1
+
+          ENDIF
+#endif
+
+!
+!--       Read the arrays, as long as the end of the file is reached
+          fanz          =         0
+          current_level =         1
+          current_var   = 999999999
+
+          DO  WHILE ( danz /= 0 )
+
+!
+!--          Loop over all files (reading data of the subdomains)
+             DO  id = 0, danz-1
+!
+!--             File from PE0 contains special information at the beginning,
+!--             concerning the lower and upper indices of the total-domain used
+!--             in PALM (nxag, nxeg, nyag, nyeg) and the lower and upper indices
+!--             of the array to be writte by this routine (nxa, nxe, nya,
+!--             nye). Usually in the horizontal directions nxag=-1 and nxa=0
+!--             while all other variables have the same value (i.e. nxeg=nxe).
+!--             Allocate necessary arrays, open the output file and write
+!--             the coordinate informations needed by ISO2D.
+                IF ( id == 0  .AND.  fanz(0) == 0  .AND.  fanz(1) == 0 )  THEN
+                   READ ( id+110 )  nxag, nxeg, nyag, nyeg
+                   READ ( id+110 )  nxa, nxe, nya, nye
+                   ALLOCATE ( eta(nya:nye), ho(nxa:nxe), hu(nxa:nxe), &
+                        pf(nxag:nxeg,nyag:nyeg) )
+                   READ ( id+110 )  dx, eta, hu, ho
+
+                   IF ( iso2d_output )  THEN
+                      OPEN ( 2, FILE='PLOT2D_'//modus//TRIM( model_string ), &
+                           FORM='UNFORMATTED' )
+                      WRITE ( 2 )  dx, eta, hu, ho
+                   ENDIF
+                ENDIF
+!
+!--             Read output time
+                IF ( netcdf_output  .AND.  id == 0 )  THEN
+                   IF ( netcdf_1 )  THEN
+                      READ ( id+110, END=998 )  simulated_time, time_step, av
+                   ELSE
+!
+!--                   For compatibility with earlier PALM versions
+                      READ ( id+110, END=998 )  simulated_time, time_step
+                      av = 0
+                   ENDIF
+                ENDIF
+!
+!--             Read subdomain indices
+                READ ( id+110, END=998 )  xa, xe, ya, ye
+!
+!--             IF the PE made no output (in case that no part of the
+!--             cross-section is situated on this PE), indices have the
+!--             value -1
+                IF ( .NOT. ( xa == -1  .AND.  xe == -1  .AND. &
+                             ya == -1  .AND.  ye == -1 ) )  THEN
+!
+!--                Read the subdomain grid-point values
+                   READ ( id+110 )  pf(xa:xe,ya:ye)
+                ENDIF
+                IF ( id == 0 )  fanz(av) = fanz(av) + 1
+
+             ENDDO
+!
+!--          Write the data of the total domain cross-section
+             IF ( iso2d_output )  WRITE ( 2 )  pf(nxa:nxe,nya:nye)
+       
+!
+!--          Write same data in NetCDF format
+             IF ( netcdf_output )  THEN
+#if defined( __netcdf )
+!
+!--             Check if a new time step has begun; if yes write data to time 
+!--             axis
+                IF ( current_var(av) > num_var(av) )  THEN
+                   current_var(av) = 1
+                   nc_stat = NF90_PUT_VAR( id_set(av), id_var_time(av), &
+                        (/ simulated_time /),        &
+                        start = (/ time_step /),     &
+                        count = (/ 1 /) )
+                   IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 7 )
+                ENDIF
+
+!
+!--             Now write the data; this is mode dependent
+                SELECT CASE ( modus )
+
+                   CASE ( 'XY' )
+                      nc_stat = NF90_PUT_VAR( id_set(av),                      &
+                                           id_var(av,current_var(av)),         &
+                                           pf(nxa:nxe,nya:nye),                &
+                             start = (/ 1, 1, current_level(av), time_step /), &
+                                      count = (/ nxe-nxa+1, nye-nya+1, 1, 1 /) )
+                      IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 8)
+                  
+                   CASE ( 'XZ' )
+                      nc_stat = NF90_PUT_VAR( id_set(av),                      &
+                                           id_var(av,current_var(av)),         &
+                                           pf(nxa:nxe,nya:nye),                &
+                             start = (/ 1, current_level(av), 1, time_step /), &
+                                      count = (/ nxe-nxa+1, 1, nye-nya+1, 1 /) )
+                      IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 9)
+
+                   CASE ( 'YZ' )
+                      nc_stat = NF90_PUT_VAR( id_set(av),                      &
+                                           id_var(av,current_var(av)),         &
+                                           pf(nxa:nxe,nya:nye),                &
+                             start = (/ current_level(av), 1, 1, time_step /), &
+                                      count = (/ 1, nxe-nxa+1, nye-nya+1, 1 /) )
+                      IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error(10)
+
+                END SELECT
+
+!
+!--             Data is written, check if max level is reached
+                current_level(av) = current_level(av) + 1
+                IF ( current_level(av) > levels(av,current_var(av)) )  THEN
+                   current_level(av) = 1
+                   current_var(av)   = current_var(av) + 1
+                ENDIF
+
+             ENDIF
+#endif
+
+          ENDDO
+
+998       IF ( danz /= 0 )  THEN
+!
+!--          Print the number of the arrays processed
+             WRITE (*,'(16X,I4,A)')  fanz(0)+fanz(1), ' array(s) processed'
+             IF ( fanz(1) /= 0 )  THEN
+                WRITE (*,'(16X,I4,A)')  fanz(1), ' array(s) are time-averaged'
+             ENDIF
+
+!
+!--          Close all files and deallocate arrays
+             DO  id = 0, danz-1
+                CLOSE ( id+110 )
+             ENDDO
+             CLOSE ( 2 )
+             DEALLOCATE ( eta, ho, hu, pf )
+          ENDIF
+
+!
+!--       Close the NetCDF file
+          IF ( netcdf_output  .AND.  danz /= 0 )  THEN
+#if defined( __netcdf )
+             IF ( netcdf_0 )  THEN
+                nc_stat = NF90_CLOSE( id_set(0) )
+                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 11 )
+             ENDIF
+             IF ( netcdf_1 )  THEN
+                nc_stat = NF90_CLOSE( id_set(1) )
+                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 12 )
+             ENDIF
+#endif
+          ENDIF
+
+!
+!--       Choose the next cross-section
+          SELECT CASE ( modus )
+             CASE ( 'XY' )
+                modus = 'XZ'
+             CASE ( 'XZ' )
+                modus = 'YZ'
+             CASE ( 'YZ' )
+                modus = 'no'
+          END SELECT
 
        ENDDO
 
-!
-!--    Inquire whether an iso2d parameter file exists
-       INQUIRE( FILE='PLOT2D_' // modus // '_GLOBAL', EXIST=iso2d_output )
+
+!
+!--    Combine the 3D-arrays
+
+!
+!--    Inquire whether an avs fld file exists
+       INQUIRE( FILE='PLOT3D_FLD'//TRIM( model_string ), EXIST=avs_output )
 
 !
 !--    Inquire whether a NetCDF file exists
-       INQUIRE( FILE='DATA_2D_' // modus // '_NETCDF', EXIST=netcdf_0 )
+       INQUIRE( FILE='DATA_3D_NETCDF'//TRIM( model_string ), EXIST=netcdf_0 )
 
 !
 !--    Inquire whether a NetCDF file for time-averaged data exists
-       INQUIRE( FILE='DATA_2D_' // modus // '_AV_NETCDF', EXIST=netcdf_1 )
+       INQUIRE( FILE='DATA_3D_AV_NETCDF'//TRIM( model_string ), EXIST=netcdf_1 )
 
        IF ( netcdf_0  .OR.  netcdf_1 )  THEN
@@ -143,6 +478,35 @@
 
 !
+!--    Check, if file from PE0 exists
+       danz = 0
+       WRITE (id_string,'(I4.4)')  danz
+       INQUIRE ( &
+            FILE='PLOT3D_DATA'//TRIM( model_string )//'_'//TRIM( id_string ),  &
+            EXIST=found )
+
+!
+!--    Combination only works, if data are not compressed. In that case, 
+!--    PALM created a flag file (PLOT3D_COMPRESSED)
+       INQUIRE ( FILE='PLOT3D_COMPRESSED'//TRIM( model_string ), &
+            EXIST=compressed )
+
+!
+!--    Find out the number of files and open them
+       DO  WHILE ( found  .AND.  .NOT. compressed )
+
+          OPEN ( danz+110, &
+               FILE='PLOT3D_DATA'//TRIM( model_string )//'_'//TRIM(id_string), &
+               FORM='UNFORMATTED')
+          danz = danz + 1
+          WRITE (id_string,'(I4.4)')  danz
+          INQUIRE ( &
+               FILE='PLOT3D_DATA'//TRIM( model_string )//'_'//TRIM(id_string), &
+               EXIST=found )
+
+       ENDDO
+
+!
 !--    Info-output
-       PRINT*, ''
+       PRINT*, ' '
        PRINT*, '*** combine_plot_fields ***'
 #if defined( __netcdf )
@@ -155,7 +519,11 @@
 #endif
        IF ( danz /= 0 )  THEN
-          PRINT*, '    ',modus,'-section:  ', danz, ' file(s) found'
+          PRINT*, '    3D-data:     ', danz, ' file(s) found'
        ELSE
-          PRINT*, '    no ', modus, '-section data available'
+          IF ( found .AND. compressed )  THEN
+             PRINT*, '+++ no 3D-data processing, since data are compressed'
+          ELSE
+             PRINT*, '    no 3D-data file available'
+          ENDIF
        ENDIF
 
@@ -170,10 +538,11 @@
 !--          Open NetCDF dataset
              IF ( av == 0 )  THEN
-                filename = 'DATA_2D_' // modus // '_NETCDF'
+                filename = 'DATA_3D_NETCDF'//TRIM( model_string )
              ELSE
-                filename = 'DATA_2D_' // modus // '_AV_NETCDF'
+                filename = 'DATA_3D_AV_NETCDF'//TRIM( model_string )
              ENDIF
              nc_stat = NF90_OPEN( filename, NF90_WRITE, id_set(av) )
-             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 1 )
+             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 13 )
+
 
 !
@@ -182,11 +551,11 @@
              var_list(av) = ' '    ! GET_ATT does not assign trailing blanks
              nc_stat = NF90_GET_ATT( id_set(av), NF90_GLOBAL, 'VAR_LIST', &
-                                     var_list(av) )
-             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 2 )
+                  var_list(av) )
+             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 14 )
 
 !
 !--          Inquire id of the time coordinate variable
              nc_stat = NF90_INQ_VARID( id_set(av), 'time', id_var_time(av) )
-             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 3 )
+             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 15 )
 
 !
@@ -199,6 +568,6 @@
 
 !
-!--          Extract the variable names from the list and inquire their
-!--          NetCDF IDs
+!--          Extract the variable names from the list and inquire their NetCDF
+!--          IDs
              pos = INDEX( var_list(av), ';' )
 !
@@ -214,30 +583,13 @@
 !
 !--             Get variable ID from name
+!                print*, '*** find id for "',TRIM( var_name ),'" begin'
                 nc_stat = NF90_INQ_VARID( id_set(av), TRIM( var_name ), &
-                                          id_var(av,i) )
-                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 4 )
-
-!
-!--             Get number of x/y/z levels for that variable
-                nc_stat = NF90_INQUIRE_VARIABLE( id_set(av), id_var(av,i), &
-                                                 dimids = id_dims_loc )
-                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 5 )
-                id_dims(av,:) = id_dims_loc
-
-!
-!--             Inquire dimension ID
-                DO  j = 1, 4
-                   nc_stat = NF90_INQUIRE_DIMENSION( id_set(av), id_dims(av,j),&
-                                                     dimname, levels(av,i) )
-                   IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 6 )
-
-                   IF ( modus == 'XY' .AND. INDEX(dimname, 'z') /= 0 )  EXIT
-                   IF ( modus == 'XZ' .AND. INDEX(dimname, 'y') /= 0 )  EXIT
-                   IF ( modus == 'YZ' .AND. INDEX(dimname, 'x') /= 0 )  EXIT
-                ENDDO
+                     id_var(av,i) )
+                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 16 )
+!                print*, '*** find id for "',TRIM( var_name ),'" end'
 
              ENDDO
 
-          ENDDO   ! av = 0, 1
+          ENDDO    ! av=0,1
 
        ENDIF
@@ -245,66 +597,54 @@
 
 !
-!--    Read the arrays, as long as the end of the file is reached
-       fanz          =         0
-       current_level =         1
-       current_var   = 999999999
-
+!--    Read arrays, until the end of the file is reached
+       current_var = 999999999
+       fanz = 0
        DO  WHILE ( danz /= 0 )
 
 !
-!--       Loop over all files (reading data of the subdomains)
+!--       Loop over all files
           DO  id = 0, danz-1
 !
 !--          File from PE0 contains special information at the beginning,
-!--          concerning the lower and upper indices of the total-domain used
-!--          in PALM (nxag, nxeg, nyag, nyeg) and the lower and upper
-!--          indices of the array to be writte by this routine (nxa, nxe, nya,
-!--          nye). Usually in the horizontal directions nxag=-1 and nxa=0
-!--          while all other variables have the same value (i.e. nxeg=nxe).
-!--          Allocate necessary arrays, open the output file and write
-!--          the coordinate informations needed by ISO2D.
+!--          concerning the lower and upper indices of the total-domain used in
+!--          PALM (nxag, nxeg, nyag, nyeg, nzag, nzeg) and the lower and upper
+!--          indices of the array to be written by this routine (nxa, nxe, nya,
+!--          nye, nza, nze). Usually nxag=-1 and nxa=0, nyag=-1 and nya=0,
+!--          nzeg=nz and nze=nz_plot3d.
+!--          Allocate necessary array and open the output file.
              IF ( id == 0  .AND.  fanz(0) == 0  .AND.  fanz(1) == 0 )  THEN
-                READ ( id+110 )  nxag, nxeg, nyag, nyeg
-                READ ( id+110 )  nxa, nxe, nya, nye
-                ALLOCATE ( eta(nya:nye), ho(nxa:nxe), hu(nxa:nxe), &
-                           pf(nxag:nxeg,nyag:nyeg) )
-                READ ( id+110 )  dx, eta, hu, ho
-
-                IF ( iso2d_output )  THEN
-                   OPEN ( 2, FILE='PLOT2D_'//modus, FORM='UNFORMATTED' )
-                   WRITE ( 2 )  dx, eta, hu, ho
+                READ ( id+110 )  nxag, nxeg, nyag, nyeg, nzag, nzeg
+                READ ( id+110 )  nxa, nxe, nya, nye, nza, nze
+                ALLOCATE ( pf3d(nxag:nxeg,nyag:nyeg,nzag:nzeg) )
+                IF ( avs_output )  THEN
+                   OPEN ( 2, FILE='PLOT3D_DATA'//TRIM( model_string ), &
+                        FORM='UNFORMATTED' )
                 ENDIF
              ENDIF
+
 !
 !--          Read output time
              IF ( netcdf_output  .AND.  id == 0 )  THEN
                 IF ( netcdf_1 )  THEN
-                   READ ( id+110, END=998 )  simulated_time, time_step, av
+                   READ ( id+110, END=999 )  simulated_time, time_step, av
                 ELSE
 !
 !--                For compatibility with earlier PALM versions
-                   READ ( id+110, END=998 )  simulated_time, time_step
+                   READ ( id+110, END=999 )  simulated_time, time_step
                    av = 0
                 ENDIF
              ENDIF
-!
-!--          Read subdomain indices
-             READ ( id+110, END=998 )  xa, xe, ya, ye
-!
-!--          IF the PE made no output (in case that no part of the
-!--          cross-section is situated on this PE), indices have the
-!--          value -1
-             IF ( .NOT. ( xa == -1  .AND.  xe == -1  .AND. &
-                          ya == -1  .AND.  ye == -1 ) )  THEN
-!
-!--              Read the subdomain grid-point values
-                 READ ( id+110 )  pf(xa:xe,ya:ye)
-             ENDIF
+
+!
+!--          Read subdomain indices and grid point values
+             READ ( id+110, END=999 )  xa, xe, ya, ye, za, ze
+             READ ( id+110 )  pf3d(xa:xe,ya:ye,za:ze)
              IF ( id == 0 )  fanz(av) = fanz(av) + 1
 
           ENDDO
-!
-!--       Write the data of the total domain cross-section
-          IF ( iso2d_output )  WRITE ( 2 )  pf(nxa:nxe,nya:nye)
+
+!
+!--       Write data of the total domain
+          IF ( avs_output )  WRITE ( 2 )  pf3d(nxa:nxe,nya:nye,nza:nze)
        
 !
@@ -317,341 +657,55 @@
                 current_var(av) = 1
                 nc_stat = NF90_PUT_VAR( id_set(av), id_var_time(av), &
-                                        (/ simulated_time /),        &
-                                        start = (/ time_step /),     &
-                                        count = (/ 1 /) )
-                IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 7 )
-             ENDIF
-
-!
-!--          Now write the data; this is mode dependent
-             SELECT CASE ( modus )
-
-                CASE ( 'XY' )
-                   nc_stat = NF90_PUT_VAR( id_set(av),                         &
-                                           id_var(av,current_var(av)),         &
-                                           pf(nxa:nxe,nya:nye),                &
-                             start = (/ 1, 1, current_level(av), time_step /), &
-                                      count = (/ nxe-nxa+1, nye-nya+1, 1, 1 /) )
-                   IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 8 )
-                  
-                CASE ( 'XZ' )
-                   nc_stat = NF90_PUT_VAR( id_set(av),                         &
-                                           id_var(av,current_var(av)),         &
-                                           pf(nxa:nxe,nya:nye),                &
-                             start = (/ 1, current_level(av), 1, time_step /), &
-                                      count = (/ nxe-nxa+1, 1, nye-nya+1, 1 /) )
-                   IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 9 )
-
-                CASE ( 'YZ' )
-                   nc_stat = NF90_PUT_VAR( id_set(av),                         &
-                                           id_var(av,current_var(av)),         &
-                                           pf(nxa:nxe,nya:nye),                &
-                             start = (/ current_level(av), 1, 1, time_step /), &
-                                      count = (/ 1, nxe-nxa+1, nye-nya+1, 1 /) )
-                   IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 10 )
-
-             END SELECT
-
-!
-!--          Data is written, check if max level is reached
-             current_level(av) = current_level(av) + 1
-             IF ( current_level(av) > levels(av,current_var(av)) )  THEN
-                current_level(av) = 1
-                current_var(av)   = current_var(av) + 1
-             ENDIF
-
-          ENDIF
-#endif
+                                     (/ simulated_time /),&
+                                     start = (/ time_step /), count = (/ 1 /) )
+                IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 17 )
+             ENDIF
+
+!
+!--          Now write the data
+             nc_stat = NF90_PUT_VAR( id_set(av), id_var(av,current_var(av)), &
+                                  pf3d(nxa:nxe,nya:nye,nza:nze),      &
+                                  start = (/ 1, 1, 1, time_step /),   &
+                              count = (/ nxe-nxa+1, nye-nya+1, nze-nza+1, 1 /) )
+             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 18 )
+
+             current_var(av) = current_var(av) + 1
+
+#endif
+          ENDIF
 
        ENDDO
 
-998    IF ( danz /= 0 )  THEN
-!
-!--       Print the number of the arrays processed
+999    IF ( danz /= 0 )  THEN
+!
+!--       Print the number of arrays processed
           WRITE (*,'(16X,I4,A)')  fanz(0)+fanz(1), ' array(s) processed'
           IF ( fanz(1) /= 0 )  THEN
              WRITE (*,'(16X,I4,A)')  fanz(1), ' array(s) are time-averaged'
           ENDIF
-
-!
-!--       Close all files and deallocate arrays
+!
+!--       Close all files and deallocate array
           DO  id = 0, danz-1
              CLOSE ( id+110 )
           ENDDO
           CLOSE ( 2 )
-          DEALLOCATE ( eta, ho, hu, pf )
+          DEALLOCATE ( pf3d )
+!
+!--       Close the NetCDF file
+          IF ( netcdf_output )  THEN
+#if defined( __netcdf )
+             IF ( netcdf_0 )  THEN
+                nc_stat = NF90_CLOSE( id_set(0) )
+                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 19 )
+             ENDIF
+             IF ( netcdf_1 )  THEN
+                nc_stat = NF90_CLOSE( id_set(1) )
+                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 20 )
+             ENDIF
+#endif
+          ENDIF
        ENDIF
 
-!
-!--    Close the NetCDF file
-       IF ( netcdf_output  .AND.  danz /= 0 )  THEN
-#if defined( __netcdf )
-          IF ( netcdf_0 )  THEN
-             nc_stat = NF90_CLOSE( id_set(0) )
-             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 11 )
-          ENDIF
-          IF ( netcdf_1 )  THEN
-             nc_stat = NF90_CLOSE( id_set(1) )
-             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 12 )
-          ENDIF
-#endif
-       ENDIF
-
-!
-!--    Choose the next cross-section
-       SELECT CASE ( modus )
-          CASE ( 'XY' )
-             modus = 'XZ'
-          CASE ( 'XZ' )
-             modus = 'YZ'
-          CASE ( 'YZ' )
-             modus = 'no'
-       END SELECT
-
-    ENDDO
-
-
-!
-!-- Combine the 3D-arrays
-
-!
-!-- Inquire whether an avs fld file exists
-    INQUIRE( FILE='PLOT3D_FLD', EXIST=avs_output )
-
-!
-!-- Inquire whether a NetCDF file exists
-    INQUIRE( FILE='DATA_3D_NETCDF', EXIST=netcdf_0 )
-
-!
-!-- Inquire whether a NetCDF file for time-averaged data exists
-    INQUIRE( FILE='DATA_3D_AV_NETCDF', EXIST=netcdf_1 )
-
-    IF ( netcdf_0  .OR.  netcdf_1 )  THEN
-       netcdf_output = .TRUE.
-    ELSE
-       netcdf_output = .FALSE.
-    ENDIF
-
-!
-!-- Check, if file from PE0 exists
-    danz = 0
-    WRITE (id_char,'(I4.4)')  danz
-    INQUIRE ( FILE='PLOT3D_DATA_'//TRIM( id_char ), EXIST=found )
-
-!
-!-- Combination only works, if data are not compressed. In that case, 
-!-- PALM created a flag file (PLOT3D_COMPRESSED)
-    INQUIRE ( FILE='PLOT3D_COMPRESSED', EXIST=compressed )
-
-!
-!-- Find out the number of files and open them
-    DO  WHILE ( found  .AND.  .NOT. compressed )
-
-       OPEN ( danz+110, FILE='PLOT3D_DATA_'//TRIM( id_char ), &
-              FORM='UNFORMATTED')
-       danz = danz + 1
-       WRITE (id_char,'(I4.4)')  danz
-       INQUIRE ( FILE='PLOT3D_DATA_'//TRIM( id_char ), EXIST=found )
-
-    ENDDO
-
-!
-!-- Info-output
-    PRINT*, ' '
-    PRINT*, '*** combine_plot_fields ***'
-#if defined( __netcdf )
-       IF ( netcdf_output )  PRINT*, '    NetCDF output enabled'
-#else
-       IF ( netcdf_output )  THEN
-          PRINT*, '--- Sorry, no NetCDF support on this host'
-          netcdf_output = .FALSE.
-       ENDIF
-#endif
-    IF ( danz /= 0 )  THEN
-       PRINT*, '    3D-data:     ', danz, ' file(s) found'
-    ELSE
-       IF ( found .AND. compressed )  THEN
-          PRINT*, '+++ no 3D-data processing, since data are compressed'
-       ELSE
-          PRINT*, '    no 3D-data file available'
-       ENDIF
-    ENDIF
-
-    IF ( netcdf_output  .AND.  danz /= 0 )  THEN
-#if defined( __netcdf )
-       DO  av = 0, 1
-
-          IF ( av == 0  .AND.  .NOT.  netcdf_0 )  CYCLE
-          IF ( av == 1  .AND.  .NOT.  netcdf_1 )  CYCLE
-
-!
-!--       Open NetCDF dataset
-          IF ( av == 0 )  THEN
-             filename = 'DATA_3D_NETCDF'
-          ELSE
-             filename = 'DATA_3D_AV_NETCDF'
-          ENDIF
-          nc_stat = NF90_OPEN( filename, NF90_WRITE, id_set(av) )
-          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 13 )
-
-
-!
-!--       Get the list of variables (order of variables corresponds with the
-!--       order of data on the binary file)
-          var_list(av) = ' '    ! GET_ATT does not assign trailing blanks
-          nc_stat = NF90_GET_ATT( id_set(av), NF90_GLOBAL, 'VAR_LIST', &
-                                  var_list(av) )
-          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 14 )
-
-!
-!--       Inquire id of the time coordinate variable
-          nc_stat = NF90_INQ_VARID( id_set(av), 'time', id_var_time(av) )
-          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 15 )
-
-!
-!--       Count number of variables; there is one more semicolon in the string
-!--       than variable names
-          num_var(av) = -1
-          DO  i = 1, LEN( var_list(av) )
-             IF ( var_list(av)(i:i) == ';' )  num_var(av) = num_var(av) + 1
-          ENDDO
-
-!
-!--       Extract the variable names from the list and inquire their NetCDF IDs
-          pos = INDEX( var_list(av), ';' )
-!
-!--       Loop over all variables
-          DO  i = 1, num_var(av)
-
-!
-!--          Extract variable name from list
-             var_list(av) = var_list(av)(pos+1:)
-             pos = INDEX( var_list(av), ';' )
-             var_name = var_list(av)(1:pos-1)
-
-!
-!--          Get variable ID from name
-!             print*, '*** find id for "',TRIM( var_name ),'" begin'
-             nc_stat = NF90_INQ_VARID( id_set(av), TRIM( var_name ), &
-                                       id_var(av,i) )
-             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 16 )
-!             print*, '*** find id for "',TRIM( var_name ),'" end'
-
-          ENDDO
-
-       ENDDO    ! av=0,1
-
-    ENDIF
-#endif
-
-!
-!-- Read arrays, until the end of the file is reached
-    current_var = 999999999
-    fanz = 0
-    DO  WHILE ( danz /= 0 )
-
-!
-!--    Loop over all files
-       DO  id = 0, danz-1
-!
-!--       File from PE0 contains special information at the beginning,
-!--       concerning the lower and upper indices of the total-domain used
-!--       in PALM (nxag, nxeg, nyag, nyeg, nzag, nzeg) and the lower and upper
-!--       indices of the array to be written by this routine (nxa, nxe, nya,
-!--       nye, nza, nze). Usually nxag=-1 and nxa=0, nyag=-1 and nya=0,
-!--       nzeg=nz and nze=nz_plot3d.
-!--       Allocate necessary array and open the output file.
-          IF ( id == 0  .AND.  fanz(0) == 0  .AND.  fanz(1) == 0 )  THEN
-             READ ( id+110 )  nxag, nxeg, nyag, nyeg, nzag, nzeg
-             READ ( id+110 )  nxa, nxe, nya, nye, nza, nze
-             ALLOCATE ( pf3d(nxag:nxeg,nyag:nyeg,nzag:nzeg) )
-             IF ( avs_output )  THEN
-                OPEN ( 2, FILE='PLOT3D_DATA', FORM='UNFORMATTED' )
-             ENDIF
-          ENDIF
-
-!
-!--       Read output time
-          IF ( netcdf_output  .AND.  id == 0 )  THEN
-             IF ( netcdf_1 )  THEN
-                READ ( id+110, END=999 )  simulated_time, time_step, av
-             ELSE
-!
-!--             For compatibility with earlier PALM versions
-                READ ( id+110, END=999 )  simulated_time, time_step
-                av = 0
-             ENDIF
-          ENDIF
-
-!
-!--       Read subdomain indices and grid point values
-          READ ( id+110, END=999 )  xa, xe, ya, ye, za, ze
-          READ ( id+110 )  pf3d(xa:xe,ya:ye,za:ze)
-          IF ( id == 0 )  fanz(av) = fanz(av) + 1
-
-       ENDDO
-
-!
-!--    Write data of the total domain
-       IF ( avs_output )  WRITE ( 2 )  pf3d(nxa:nxe,nya:nye,nza:nze)
-       
-!
-!--    Write same data in NetCDF format
-       IF ( netcdf_output )  THEN
-#if defined( __netcdf )
-!
-!--       Check if a new time step has begun; if yes write data to time axis
-          IF ( current_var(av) > num_var(av) )  THEN
-             current_var(av) = 1
-             nc_stat = NF90_PUT_VAR( id_set(av), id_var_time(av), &
-                                     (/ simulated_time /),&
-                                     start = (/ time_step /), count = (/ 1 /) )
-             IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 17 )
-          ENDIF
-
-!
-!--       Now write the data
-          nc_stat = NF90_PUT_VAR( id_set(av), id_var(av,current_var(av)), &
-                                  pf3d(nxa:nxe,nya:nye,nza:nze),      &
-                                  start = (/ 1, 1, 1, time_step /),   &
-                              count = (/ nxe-nxa+1, nye-nya+1, nze-nza+1, 1 /) )
-          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 18 )
-
-          current_var(av) = current_var(av) + 1
-
-#endif
-       ENDIF
-
-    ENDDO
-
-999 IF ( danz /= 0 )  THEN
-!
-!--    Print the number of arrays processed
-       WRITE (*,'(16X,I4,A)')  fanz(0)+fanz(1), ' array(s) processed'
-       IF ( fanz(1) /= 0 )  THEN
-          WRITE (*,'(16X,I4,A)')  fanz(1), ' array(s) are time-averaged'
-       ENDIF
-!
-!--    Close all files and deallocate array
-       DO  id = 0, danz-1
-          CLOSE ( id+110 )
-       ENDDO
-       CLOSE ( 2 )
-       DEALLOCATE ( pf3d )
-!
-!--    Close the NetCDF file
-       IF ( netcdf_output )  THEN
-#if defined( __netcdf )
-          IF ( netcdf_0 )  THEN
-             nc_stat = NF90_CLOSE( id_set(0) )
-             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 19 )
-          ENDIF
-          IF ( netcdf_1 )  THEN
-             nc_stat = NF90_CLOSE( id_set(1) )
-             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 20 )
-          ENDIF
-#endif
-       ENDIF
-    ENDIF
+    ENDDO  ! models
 
 
