Ignore:
Timestamp:
Feb 8, 2019 2:52:10 PM (3 years ago)
Author:
gronemeier
Message:

New: surface output available in NetCDF format (Makefile, netcdf_interface_mod, surface_data_output_mod)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/surface_data_output_mod.f90

    r3691 r3727  
    2525! -----------------
    2626! $Id$
     27! Enable NetCDF output for surface data (suehring, gronemeier)
     28!
     29! 3691 2019-01-23 09:57:04Z suehring
    2730! Add output of surface-parallel flow speed
    2831!
     
    7780       
    7881   USE control_parameters,                                                     &
    79        ONLY:  surface_output
     82       ONLY:  coupling_char, data_output_during_spinup, end_time,              &
     83              message_string, run_description_header, simulated_time_at_begin, &
     84              spinup_time, surface_output
    8085       
    8186   USE grid_variables,                                                         &
     
    8590       ONLY: nxl, nxr, nys, nyn, nzb, nzt
    8691       
     92#if defined( __netcdf )
     93   USE NETCDF
     94#endif
     95
     96   USE netcdf_data_input_mod,                                                  &
     97       ONLY:  init_model
     98
     99   USE netcdf_interface,                                                       &
     100       ONLY:  netcdf_create_att, netcdf_create_dim, netcdf_create_file,        &
     101              netcdf_create_global_atts, netcdf_create_var, netcdf_handle_error
     102       
    87103   USE pegrid
    88104       
     
    100116      INTEGER(iwp) ::  npoints_total  !< total number of points / vertices which define a surface element
    101117     
     118      INTEGER(iwp), DIMENSION(:), ALLOCATABLE   ::  s        !< coordinate for NetCDF output, number of the surface element
     119     
    102120      REAL(wp) ::  fillvalue = -9999.9_wp !< fillvalue for surface elements which are not defined
    103121     
     122      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  azimuth  !< azimuth orientation coordinate for NetCDF output
     123      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  es_utm   !< E-UTM coordinate for NetCDF output
     124      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  ns_utm   !< E-UTM coordinate for NetCDF output
     125      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  xs       !< x-coordinate for NetCDF output
     126      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  ys       !< y-coordinate for NetCDF output
     127      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  zs       !< z-coordinate for NetCDF output
     128      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  zenith   !< zenith orientation coordinate for NetCDF output
    104129      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  var_out  !< output variables
    105130      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  var_av   !< variables used for averaging
     
    109134   
    110135   CHARACTER(LEN=100), DIMENSION(300)       ::  data_output_surf = ' '  !< namelist variable which describes the output variables
    111    CHARACTER(LEN=100), DIMENSION(0:1,300)   ::  dosurf = ' '            !< internal variable which describes the output variables and separates
    112                                                                         !< averaged from non-averaged output
     136   CHARACTER(LEN=100), DIMENSION(0:1,300)   ::  dosurf = ' '            !< internal variable which describes the output variables
     137                                                                        !<  and separates averaged from non-averaged output
     138   CHARACTER(LEN=100), DIMENSION(0:1,300)   ::  dosurf_unit = ' '       !< internal variable which holds the unit of the given output variable
    113139   
    114    INTEGER(iwp) ::  average_count_surf = 0                !< number of ensemble members used for averaging
    115    INTEGER(iwp) ::  dosurf_no(0:1) = 0                    !< number of surface output quantities
     140   INTEGER(iwp) ::  average_count_surf = 0   !< number of ensemble members used for averaging
     141   INTEGER(iwp) ::  dosurf_no(0:1) = 0       !< number of surface output quantities
     142
     143   INTEGER(iwp) :: nc_stat                   !< error code for netcdf routines
     144   INTEGER(iwp), SAVE ::  oldmode            !< save old set-fill-mode of netcdf file (not needed, but required for routine call)
     145
     146   INTEGER(iwp), DIMENSION(0:1) ::  id_dim_s_surf         !< netcdf ID for dimension s
     147   INTEGER(iwp), DIMENSION(0:1) ::  id_dim_time_surf      !< netcdf ID for dimension time
     148   INTEGER(iwp), DIMENSION(0:1) ::  id_set_surf           !< netcdf ID for file
     149   INTEGER(iwp), DIMENSION(0:1) ::  id_var_etum_surf      !< netcdf ID for variable Es_UTM
     150   INTEGER(iwp), DIMENSION(0:1) ::  id_var_nutm_surf      !< netcdf ID for variable Ns_UTM
     151   INTEGER(iwp), DIMENSION(0:1) ::  id_var_time_surf      !< netcdf ID for variable time
     152   INTEGER(iwp), DIMENSION(0:1) ::  id_var_s_surf         !< netcdf ID for variable s
     153   INTEGER(iwp), DIMENSION(0:1) ::  id_var_xs_surf        !< netcdf ID for variable xs
     154   INTEGER(iwp), DIMENSION(0:1) ::  id_var_ys_surf        !< netcdf ID for variable ys
     155   INTEGER(iwp), DIMENSION(0:1) ::  id_var_zs_surf        !< netcdf ID for variable zs
     156   INTEGER(iwp), DIMENSION(0:1) ::  dosurf_time_count = 0 !< count of output time steps
     157   INTEGER(iwp), DIMENSION(0:1) ::  ntdim_surf            !< number of output time steps
     158   
     159   INTEGER(iwp), DIMENSION(0:1,300) ::  id_var_dosurf     !< netcdf ID for output variables
    116160
    117161   LOGICAL :: first_output(0:1) = .FALSE.                 !< true if first output was already called
     162   LOGICAL :: to_netcdf = .TRUE.                          !< flag indicating parallel NetCDF output
     163   LOGICAL :: to_vtk = .TRUE.                             !< flag indicating binary surface-data output that can be further
     164                                                          !< processed to VTK format
    118165
    119166   REAL(wp) ::  averaging_interval_surf  = 9999999.9_wp   !< averaging interval
     
    125172   REAL(wp) ::  time_dosurf_av = 0.0_wp                   !< internal counter variable to check for averaged data output
    126173
    127    
    128174   TYPE(surf_out) ::  surfaces      !< variable which contains all required output information
    129175   
     
    178224      IMPLICIT NONE
    179225
     226      CHARACTER (LEN=100)  :: filename            !< name of output file
     227      CHARACTER (LEN=80)   ::  time_average_text  !< string written to file attribute time_avg
     228      CHARACTER (LEN=4000) ::  var_list           !< list of variables written to NetCDF file
     229
     230      INTEGER(iwp) ::  av                !< flag for averaged (=1) and non-averaged (=0) data
    180231      INTEGER(iwp) ::  i                 !< grid index in x-direction, also running variable for counting non-average data output
    181232      INTEGER(iwp) ::  j                 !< grid index in y-direction, also running variable for counting average data output
     
    183234      INTEGER(iwp) ::  l                 !< running index for surface-element orientation
    184235      INTEGER(iwp) ::  m                 !< running index for surface elements
     236      INTEGER(iwp) ::  mm                !< local counting variable for surface elements
    185237      INTEGER(iwp) ::  npg               !< counter variable for all surface elements ( or polygons )
    186238      INTEGER(iwp) ::  point_index_count !< local counter variable for point index
     239      INTEGER(iwp) ::  start_count       !< local start counter for the surface index
    187240     
    188241      INTEGER(iwp), DIMENSION(0:numprocs-1) :: num_points_on_pe   !< array which contains the number of points on all mpi ranks
     242      INTEGER(iwp), DIMENSION(0:numprocs-1) :: num_surfaces_on_pe !< array which contains the number of surfaces on all mpi ranks
    189243      INTEGER(iwp), ALLOCATABLE, DIMENSION(:,:,:) ::  point_index !< dummy array used to check where the reference points for surface polygons are located
    190244           
     245      REAL(wp) ::  az    !< azimuth angle, indicated the vertical orientation of a surface element
     246      REAL(wp) ::  off_x !< grid offset in x-direction between the stored grid index and the actual wall
     247      REAL(wp) ::  off_y !< grid offset in y-direction between the stored grid index and the actual wall
     248
     249      REAL(wp), DIMENSION(:), ALLOCATABLE ::  netcdf_data_1d  !< dummy array to output 1D data into netcdf file
     250
     251
    191252!
    192253!--   Determine the number of surface elements on subdomain
     
    198259                  + surf_def_v(3)%ns + surf_lsm_v(3)%ns + surf_usm_v(3)%ns       !eastward-facing   
    199260!
     261!--    Determine the total number of surfaces in the model domain
     262#if defined( __parallel )
     263       CALL MPI_ALLREDUCE( surfaces%ns, surfaces%ns_total, 1,                  &
     264                           MPI_INTEGER, MPI_SUM, comm2d, ierr )
     265#else
     266       surfaces%ns_total = surfaces%ns
     267#endif
     268!
    200269!--   Allocate output variable and set to _FillValue attribute
    201270      ALLOCATE ( surfaces%var_out(1:surfaces%ns) )
     
    209278      ENDIF
    210279!
    211 !--   Initialize points and polygon data.
     280!--   If output to VTK format is enabled, initialize point and polygon data.
    212281!--   In a first step, count the number of points which are defining
    213282!--   the surfaces and the polygons.
    214       ALLOCATE( point_index(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
    215       point_index = -1
    216 !
    217 !--   Horizontal default surfaces
    218       surfaces%npoints = 0
    219       DO  l = 0, 1
    220          DO  m = 1, surf_def_h(0)%ns
    221 !
    222 !--         Determine the indices of the respective grid cell inside the topography
    223             i = surf_def_h(0)%i(m) + surf_def_h(0)%ioff
    224             j = surf_def_h(0)%j(m) + surf_def_h(0)%joff
    225             k = surf_def_h(0)%k(m) + surf_def_h(0)%koff
    226 !
    227 !--         Check if the vertices that define the surface element are already
    228 !--         defined, if not, increment the counter.
     283      IF ( to_vtk )  THEN
     284         ALLOCATE( point_index(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
     285         point_index = -1
     286!
     287!--      Horizontal default surfaces
     288         surfaces%npoints = 0
     289         DO  l = 0, 1
     290            DO  m = 1, surf_def_h(0)%ns
     291!
     292!--            Determine the indices of the respective grid cell inside the topography
     293               i = surf_def_h(0)%i(m) + surf_def_h(0)%ioff
     294               j = surf_def_h(0)%j(m) + surf_def_h(0)%joff
     295               k = surf_def_h(0)%k(m) + surf_def_h(0)%koff
     296!
     297!--            Check if the vertices that define the surface element are already
     298!--            defined, if not, increment the counter.
     299               IF ( point_index(k,j,i) < 0 )  THEN
     300                  surfaces%npoints   = surfaces%npoints + 1   
     301                  point_index(k,j,i) = surfaces%npoints - 1
     302               ENDIF
     303               IF ( point_index(k,j,i+1) < 0 )  THEN
     304                  surfaces%npoints     = surfaces%npoints + 1   
     305                  point_index(k,j,i+1) = surfaces%npoints - 1
     306               ENDIF
     307               IF ( point_index(k,j+1,i+1) < 0 )  THEN
     308                  surfaces%npoints       = surfaces%npoints + 1   
     309                  point_index(k,j+1,i+1) = surfaces%npoints - 1
     310               ENDIF
     311               IF ( point_index(k,j+1,i) < 0 )  THEN
     312                  surfaces%npoints     = surfaces%npoints + 1   
     313                  point_index(k,j+1,i) = surfaces%npoints - 1
     314               ENDIF           
     315            ENDDO
     316         ENDDO
     317         DO  m = 1, surf_lsm_h%ns
     318            i = surf_lsm_h%i(m) + surf_lsm_h%ioff
     319            j = surf_lsm_h%j(m) + surf_lsm_h%joff
     320            k = surf_lsm_h%k(m) + surf_lsm_h%koff
     321
    229322            IF ( point_index(k,j,i) < 0 )  THEN
    230323               surfaces%npoints   = surfaces%npoints + 1   
     
    242335               surfaces%npoints     = surfaces%npoints + 1   
    243336               point_index(k,j+1,i) = surfaces%npoints - 1
    244             ENDIF           
     337            ENDIF   
    245338         ENDDO
    246       ENDDO
    247       DO  m = 1, surf_lsm_h%ns
    248          i = surf_lsm_h%i(m) + surf_lsm_h%ioff
    249          j = surf_lsm_h%j(m) + surf_lsm_h%joff
    250          k = surf_lsm_h%k(m) + surf_lsm_h%koff
    251 
    252          IF ( point_index(k,j,i) < 0 )  THEN
    253             surfaces%npoints   = surfaces%npoints + 1   
    254             point_index(k,j,i) = surfaces%npoints - 1
    255          ENDIF
    256          IF ( point_index(k,j,i+1) < 0 )  THEN
    257             surfaces%npoints     = surfaces%npoints + 1   
    258             point_index(k,j,i+1) = surfaces%npoints - 1
    259          ENDIF
    260          IF ( point_index(k,j+1,i+1) < 0 )  THEN
    261             surfaces%npoints       = surfaces%npoints + 1   
    262             point_index(k,j+1,i+1) = surfaces%npoints - 1
    263          ENDIF
    264          IF ( point_index(k,j+1,i) < 0 )  THEN
    265             surfaces%npoints     = surfaces%npoints + 1   
    266             point_index(k,j+1,i) = surfaces%npoints - 1
    267          ENDIF   
    268       ENDDO
    269       DO  m = 1, surf_usm_h%ns
    270          i = surf_usm_h%i(m) + surf_usm_h%ioff
    271          j = surf_usm_h%j(m) + surf_usm_h%joff
    272          k = surf_usm_h%k(m) + surf_usm_h%koff
    273 
    274          IF ( point_index(k,j,i) < 0 )  THEN
    275             surfaces%npoints   = surfaces%npoints + 1   
    276             point_index(k,j,i) = surfaces%npoints - 1
    277          ENDIF
    278          IF ( point_index(k,j,i+1) < 0 )  THEN
    279             surfaces%npoints     = surfaces%npoints + 1   
    280             point_index(k,j,i+1) = surfaces%npoints - 1
    281          ENDIF
    282          IF ( point_index(k,j+1,i+1) < 0 )  THEN
    283             surfaces%npoints       = surfaces%npoints + 1   
    284             point_index(k,j+1,i+1) = surfaces%npoints - 1
    285          ENDIF
    286          IF ( point_index(k,j+1,i) < 0 )  THEN
    287             surfaces%npoints     = surfaces%npoints + 1   
    288             point_index(k,j+1,i) = surfaces%npoints - 1
    289          ENDIF     
    290       ENDDO
    291 !
    292 !--   Vertical surfaces
    293       DO  l = 0, 3
    294          DO  m = 1, surf_def_v(l)%ns
    295 !
    296 !--         Determine the indices of the respective grid cell inside the topography.
    297 !--         Please note, j-index for north-facing surfaces ( l==0 ) is
    298 !--         identical to the reference j-index outside the grid box.
    299 !--         Equivalent for east-facing surfaces and i-index.
    300             i = surf_def_v(l)%i(m) + MERGE( surf_def_v(l)%ioff, 0, l == 3 )
    301             j = surf_def_v(l)%j(m) + MERGE( surf_def_v(l)%joff, 0, l == 1 )
    302             k = surf_def_v(l)%k(m) + surf_def_v(l)%koff
    303 !
    304 !--         Lower left /front vertex
     339         DO  m = 1, surf_usm_h%ns
     340            i = surf_usm_h%i(m) + surf_usm_h%ioff
     341            j = surf_usm_h%j(m) + surf_usm_h%joff
     342            k = surf_usm_h%k(m) + surf_usm_h%koff
     343
    305344            IF ( point_index(k,j,i) < 0 )  THEN
    306345               surfaces%npoints   = surfaces%npoints + 1   
    307346               point_index(k,j,i) = surfaces%npoints - 1
    308             ENDIF
    309 !
    310 !--         Upper / lower right index for north- and south-facing surfaces
    311             IF ( l == 0  .OR.  l == 1 )  THEN
     347            ENDIF
     348            IF ( point_index(k,j,i+1) < 0 )  THEN
     349               surfaces%npoints     = surfaces%npoints + 1   
     350               point_index(k,j,i+1) = surfaces%npoints - 1
     351            ENDIF
     352            IF ( point_index(k,j+1,i+1) < 0 )  THEN
     353               surfaces%npoints       = surfaces%npoints + 1   
     354               point_index(k,j+1,i+1) = surfaces%npoints - 1
     355            ENDIF
     356            IF ( point_index(k,j+1,i) < 0 )  THEN
     357               surfaces%npoints     = surfaces%npoints + 1   
     358               point_index(k,j+1,i) = surfaces%npoints - 1
     359            ENDIF     
     360         ENDDO
     361!
     362!--      Vertical surfaces
     363         DO  l = 0, 3
     364            DO  m = 1, surf_def_v(l)%ns
     365!
     366!--            Determine the indices of the respective grid cell inside the
     367!--            topography. Please note, j-index for north-facing surfaces
     368!--            ( l==0 ) is identical to the reference j-index outside the grid 
     369!--            box. Equivalent for east-facing surfaces and i-index.
     370               i = surf_def_v(l)%i(m) + MERGE( surf_def_v(l)%ioff, 0, l == 3 )
     371               j = surf_def_v(l)%j(m) + MERGE( surf_def_v(l)%joff, 0, l == 1 )
     372               k = surf_def_v(l)%k(m) + surf_def_v(l)%koff
     373!
     374!--            Lower left /front vertex
     375               IF ( point_index(k,j,i) < 0 )  THEN
     376                  surfaces%npoints   = surfaces%npoints + 1   
     377                  point_index(k,j,i) = surfaces%npoints - 1
     378               ENDIF
     379!
     380!--            Upper / lower right index for north- and south-facing surfaces
     381               IF ( l == 0  .OR.  l == 1 )  THEN
     382                  IF ( point_index(k,j,i+1) < 0 )  THEN
     383                     surfaces%npoints     = surfaces%npoints + 1   
     384                     point_index(k,j,i+1) = surfaces%npoints - 1
     385                  ENDIF
     386                  IF ( point_index(k+1,j,i+1) < 0 )  THEN
     387                     surfaces%npoints       = surfaces%npoints + 1   
     388                     point_index(k+1,j,i+1) = surfaces%npoints - 1
     389                  ENDIF
     390!
     391!--            Upper / lower front index for east- and west-facing surfaces
     392               ELSEIF ( l == 2  .OR.  l == 3 )  THEN
     393                  IF ( point_index(k,j+1,i) < 0 )  THEN
     394                     surfaces%npoints     = surfaces%npoints + 1   
     395                     point_index(k,j+1,i) = surfaces%npoints - 1
     396                  ENDIF
     397                  IF ( point_index(k+1,j+1,i) < 0 )  THEN
     398                     surfaces%npoints       = surfaces%npoints + 1   
     399                     point_index(k+1,j+1,i) = surfaces%npoints - 1
     400                  ENDIF
     401               ENDIF
     402!
     403!--            Upper left / front vertex
     404               IF ( point_index(k+1,j,i) < 0 )  THEN
     405                  surfaces%npoints     = surfaces%npoints + 1   
     406                  point_index(k+1,j,i) = surfaces%npoints - 1
     407               ENDIF
     408            ENDDO
     409            DO  m = 1, surf_lsm_v(l)%ns
     410               i = surf_lsm_v(l)%i(m) + MERGE( surf_lsm_v(l)%ioff, 0, l == 3 )
     411               j = surf_lsm_v(l)%j(m) + MERGE( surf_lsm_v(l)%joff, 0, l == 1 )
     412               k = surf_lsm_v(l)%k(m) + surf_lsm_v(l)%koff
     413!
     414!--            Lower left /front vertex
     415               IF ( point_index(k,j,i) < 0 )  THEN
     416                  surfaces%npoints   = surfaces%npoints + 1   
     417                  point_index(k,j,i) = surfaces%npoints - 1
     418               ENDIF
     419!
     420!--            Upper / lower right index for north- and south-facing surfaces
     421               IF ( l == 0  .OR.  l == 1 )  THEN
     422                  IF ( point_index(k,j,i+1) < 0 )  THEN
     423                     surfaces%npoints     = surfaces%npoints + 1   
     424                     point_index(k,j,i+1) = surfaces%npoints - 1
     425                  ENDIF
     426                  IF ( point_index(k+1,j,i+1) < 0 )  THEN
     427                     surfaces%npoints       = surfaces%npoints + 1   
     428                     point_index(k+1,j,i+1) = surfaces%npoints - 1
     429                  ENDIF
     430!
     431!--            Upper / lower front index for east- and west-facing surfaces
     432               ELSEIF ( l == 2  .OR.  l == 3 )  THEN
     433                  IF ( point_index(k,j+1,i) < 0 )  THEN
     434                     surfaces%npoints     = surfaces%npoints + 1   
     435                     point_index(k,j+1,i) = surfaces%npoints - 1
     436                  ENDIF
     437                  IF ( point_index(k+1,j+1,i) < 0 )  THEN
     438                     surfaces%npoints       = surfaces%npoints + 1   
     439                     point_index(k+1,j+1,i) = surfaces%npoints - 1
     440                  ENDIF
     441               ENDIF
     442!
     443!--            Upper left / front vertex
     444               IF ( point_index(k+1,j,i) < 0 )  THEN
     445                  surfaces%npoints     = surfaces%npoints + 1   
     446                  point_index(k+1,j,i) = surfaces%npoints - 1
     447               ENDIF
     448            ENDDO
     449
     450            DO  m = 1, surf_usm_v(l)%ns
     451               i = surf_usm_v(l)%i(m) + MERGE( surf_usm_v(l)%ioff, 0, l == 3 )
     452               j = surf_usm_v(l)%j(m) + MERGE( surf_usm_v(l)%joff, 0, l == 1 )
     453               k = surf_usm_v(l)%k(m) + surf_usm_v(l)%koff
     454!
     455!--            Lower left /front vertex
     456               IF ( point_index(k,j,i) < 0 )  THEN
     457                  surfaces%npoints   = surfaces%npoints + 1   
     458                  point_index(k,j,i) = surfaces%npoints - 1
     459               ENDIF
     460!
     461!--            Upper / lower right index for north- and south-facing surfaces
     462               IF ( l == 0  .OR.  l == 1 )  THEN
     463                  IF ( point_index(k,j,i+1) < 0 )  THEN
     464                     surfaces%npoints     = surfaces%npoints + 1   
     465                     point_index(k,j,i+1) = surfaces%npoints - 1
     466                  ENDIF
     467                  IF ( point_index(k+1,j,i+1) < 0 )  THEN
     468                     surfaces%npoints       = surfaces%npoints + 1   
     469                     point_index(k+1,j,i+1) = surfaces%npoints - 1
     470                  ENDIF
     471!
     472!--            Upper / lower front index for east- and west-facing surfaces
     473               ELSEIF ( l == 2  .OR.  l == 3 )  THEN
     474                  IF ( point_index(k,j+1,i) < 0 )  THEN
     475                     surfaces%npoints     = surfaces%npoints + 1   
     476                     point_index(k,j+1,i) = surfaces%npoints - 1
     477                  ENDIF
     478                  IF ( point_index(k+1,j+1,i) < 0 )  THEN
     479                     surfaces%npoints       = surfaces%npoints + 1   
     480                     point_index(k+1,j+1,i) = surfaces%npoints - 1
     481                  ENDIF
     482               ENDIF
     483!
     484!--            Upper left / front vertex
     485               IF ( point_index(k+1,j,i) < 0 )  THEN
     486                  surfaces%npoints     = surfaces%npoints + 1   
     487                  point_index(k+1,j,i) = surfaces%npoints - 1
     488               ENDIF
     489            ENDDO
     490
     491         ENDDO
     492!
     493!--      Allocate the number of points and polygons. Note, the number of polygons
     494!--      is identical to the number of surfaces elements, whereas the number
     495!--      of points (vertices), which define the polygons, can be larger.
     496         ALLOCATE( surfaces%points(3,1:surfaces%npoints) )
     497         ALLOCATE( surfaces%polygons(5,1:surfaces%ns)    )
     498!
     499!--      Note, PARAVIEW expects consecutively ordered points, in order to
     500!--      unambiguously identify surfaces.
     501!--      Hence, all PEs should know where they start counting, depending on the
     502!--      number of points on the other PE's with lower MPI rank.
     503#if defined( __parallel )
     504         CALL MPI_ALLGATHER( surfaces%npoints, 1, MPI_INTEGER,                 &
     505                             num_points_on_pe, 1, MPI_INTEGER, comm2d, ierr  )
     506#else
     507         num_points_on_pe = surfaces%npoints
     508#endif
     509
     510!
     511!--      After the number of vertices is counted, repeat the loops and define the
     512!--      vertices. Start with the horizontal default surfaces.
     513!--      First, however, determine the offset where couting of points should be
     514!--      started, which is the sum of points of all PE's with lower MPI rank.
     515         i                 = 0
     516         point_index_count = 0
     517         DO WHILE ( i < myid  .AND.  i <= SIZE( num_points_on_pe ) )
     518            point_index_count = point_index_count + num_points_on_pe(i)
     519            i                 = i + 1
     520         ENDDO
     521         
     522         surfaces%npoints = 0
     523         point_index      = -1
     524         npg              = 0
     525     
     526         DO  l = 0, 1
     527            DO  m = 1, surf_def_h(0)%ns
     528!
     529!--            Determine the indices of the respective grid cell inside the
     530!--            topography.
     531               i = surf_def_h(0)%i(m) + surf_def_h(0)%ioff
     532               j = surf_def_h(0)%j(m) + surf_def_h(0)%joff
     533               k = surf_def_h(0)%k(m) + surf_def_h(0)%koff
     534!
     535!--            Check if the vertices that define the surface element are 
     536!--            already defined, if not, increment the counter.
     537               IF ( point_index(k,j,i) < 0 )  THEN
     538                  surfaces%npoints   = surfaces%npoints + 1   
     539                  point_index(k,j,i) = point_index_count
     540                  point_index_count  = point_index_count + 1               
     541                  surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
     542                  surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
     543                  surfaces%points(3,surfaces%npoints) = zw(k)
     544               ENDIF
    312545               IF ( point_index(k,j,i+1) < 0 )  THEN
    313546                  surfaces%npoints     = surfaces%npoints + 1   
    314                   point_index(k,j,i+1) = surfaces%npoints - 1
    315                ENDIF
    316                IF ( point_index(k+1,j,i+1) < 0 )  THEN
     547                  point_index(k,j,i+1) = point_index_count
     548                  point_index_count    = point_index_count + 1
     549                  surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
     550                  surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
     551                  surfaces%points(3,surfaces%npoints) = zw(k)
     552               ENDIF
     553               IF ( point_index(k,j+1,i+1) < 0 )  THEN
    317554                  surfaces%npoints       = surfaces%npoints + 1   
    318                   point_index(k+1,j,i+1) = surfaces%npoints - 1
    319                ENDIF
    320 !
    321 !--         Upper / lower front index for east- and west-facing surfaces
    322             ELSEIF ( l == 2  .OR.  l == 3 )  THEN
     555                  point_index(k,j+1,i+1) = point_index_count
     556                  point_index_count      = point_index_count + 1
     557                  surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
     558                  surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
     559                  surfaces%points(3,surfaces%npoints) = zw(k)
     560               ENDIF
    323561               IF ( point_index(k,j+1,i) < 0 )  THEN
    324562                  surfaces%npoints     = surfaces%npoints + 1   
    325                   point_index(k,j+1,i) = surfaces%npoints - 1
    326                ENDIF
    327                IF ( point_index(k+1,j+1,i) < 0 )  THEN
    328                   surfaces%npoints       = surfaces%npoints + 1   
    329                   point_index(k+1,j+1,i) = surfaces%npoints - 1
    330                ENDIF
    331             ENDIF
    332 !
    333 !--         Upper left / front vertex
    334             IF ( point_index(k+1,j,i) < 0 )  THEN
    335                surfaces%npoints     = surfaces%npoints + 1   
    336                point_index(k+1,j,i) = surfaces%npoints - 1
    337             ENDIF
     563                  point_index(k,j+1,i) = point_index_count
     564                  point_index_count    = point_index_count + 1
     565                  surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
     566                  surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
     567                  surfaces%points(3,surfaces%npoints) = zw(k)
     568               ENDIF
     569               
     570               npg                        = npg + 1
     571               surfaces%polygons(1,npg)   = 4
     572               surfaces%polygons(2,npg)   = point_index(k,j,i)
     573               surfaces%polygons(3,npg)   = point_index(k,j,i+1)
     574               surfaces%polygons(4,npg)   = point_index(k,j+1,i+1)
     575               surfaces%polygons(5,npg)   = point_index(k,j+1,i)
     576            ENDDO
    338577         ENDDO
    339          DO  m = 1, surf_lsm_v(l)%ns
    340             i = surf_lsm_v(l)%i(m) + MERGE( surf_lsm_v(l)%ioff, 0, l == 3 )
    341             j = surf_lsm_v(l)%j(m) + MERGE( surf_lsm_v(l)%joff, 0, l == 1 )
    342             k = surf_lsm_v(l)%k(m) + surf_lsm_v(l)%koff
    343 !
    344 !--         Lower left /front vertex
    345             IF ( point_index(k,j,i) < 0 )  THEN
    346                surfaces%npoints   = surfaces%npoints + 1   
    347                point_index(k,j,i) = surfaces%npoints - 1
    348             ENDIF
    349 !
    350 !--         Upper / lower right index for north- and south-facing surfaces
    351             IF ( l == 0  .OR.  l == 1 )  THEN
    352                IF ( point_index(k,j,i+1) < 0 )  THEN
    353                   surfaces%npoints     = surfaces%npoints + 1   
    354                   point_index(k,j,i+1) = surfaces%npoints - 1
    355                ENDIF
    356                IF ( point_index(k+1,j,i+1) < 0 )  THEN
    357                   surfaces%npoints       = surfaces%npoints + 1   
    358                   point_index(k+1,j,i+1) = surfaces%npoints - 1
    359                ENDIF
    360 !
    361 !--         Upper / lower front index for east- and west-facing surfaces
    362             ELSEIF ( l == 2  .OR.  l == 3 )  THEN
    363                IF ( point_index(k,j+1,i) < 0 )  THEN
    364                   surfaces%npoints     = surfaces%npoints + 1   
    365                   point_index(k,j+1,i) = surfaces%npoints - 1
    366                ENDIF
    367                IF ( point_index(k+1,j+1,i) < 0 )  THEN
    368                   surfaces%npoints       = surfaces%npoints + 1   
    369                   point_index(k+1,j+1,i) = surfaces%npoints - 1
    370                ENDIF
    371             ENDIF
    372 !
    373 !--         Upper left / front vertex
    374             IF ( point_index(k+1,j,i) < 0 )  THEN
    375                surfaces%npoints     = surfaces%npoints + 1   
    376                point_index(k+1,j,i) = surfaces%npoints - 1
    377             ENDIF
    378          ENDDO
    379 
    380          DO  m = 1, surf_usm_v(l)%ns
    381             i = surf_usm_v(l)%i(m) + MERGE( surf_usm_v(l)%ioff, 0, l == 3 )
    382             j = surf_usm_v(l)%j(m) + MERGE( surf_usm_v(l)%joff, 0, l == 1 )
    383             k = surf_usm_v(l)%k(m) + surf_usm_v(l)%koff
    384 !
    385 !--         Lower left /front vertex
    386             IF ( point_index(k,j,i) < 0 )  THEN
    387                surfaces%npoints   = surfaces%npoints + 1   
    388                point_index(k,j,i) = surfaces%npoints - 1
    389             ENDIF
    390 !
    391 !--         Upper / lower right index for north- and south-facing surfaces
    392             IF ( l == 0  .OR.  l == 1 )  THEN
    393                IF ( point_index(k,j,i+1) < 0 )  THEN
    394                   surfaces%npoints     = surfaces%npoints + 1   
    395                   point_index(k,j,i+1) = surfaces%npoints - 1
    396                ENDIF
    397                IF ( point_index(k+1,j,i+1) < 0 )  THEN
    398                   surfaces%npoints       = surfaces%npoints + 1   
    399                   point_index(k+1,j,i+1) = surfaces%npoints - 1
    400                ENDIF
    401 !
    402 !--         Upper / lower front index for east- and west-facing surfaces
    403             ELSEIF ( l == 2  .OR.  l == 3 )  THEN
    404                IF ( point_index(k,j+1,i) < 0 )  THEN
    405                   surfaces%npoints     = surfaces%npoints + 1   
    406                   point_index(k,j+1,i) = surfaces%npoints - 1
    407                ENDIF
    408                IF ( point_index(k+1,j+1,i) < 0 )  THEN
    409                   surfaces%npoints       = surfaces%npoints + 1   
    410                   point_index(k+1,j+1,i) = surfaces%npoints - 1
    411                ENDIF
    412             ENDIF
    413 !
    414 !--         Upper left / front vertex
    415             IF ( point_index(k+1,j,i) < 0 )  THEN
    416                surfaces%npoints     = surfaces%npoints + 1   
    417                point_index(k+1,j,i) = surfaces%npoints - 1
    418             ENDIF
    419          ENDDO
    420 
    421       ENDDO
    422 !
    423 !--   Allocate the number of points and polygons. Note, the number of polygons
    424 !--   is identical to the number of surfaces elements, whereas the number
    425 !--   of points (vertices) which define the polygons can be larger.
    426       ALLOCATE( surfaces%points(3,1:surfaces%npoints) )
    427       ALLOCATE( surfaces%polygons(5,1:surfaces%ns)    )
    428 !
    429 !--   Note, PARAVIEW expects consecutively ordered points which can be
    430 !--   unambiguously identified.
    431 !--   Hence, all PEs know where they should start counting, depending on the
    432 !--   number of points on the other PE's with lower MPI rank.
    433 #if defined( __parallel )
    434       CALL MPI_ALLGATHER( surfaces%npoints, 1, MPI_INTEGER,                    &
    435                           num_points_on_pe, 1, MPI_INTEGER, comm2d, ierr  )
    436 #else
    437       num_points_on_pe = surfaces%npoints
    438 #endif
    439 
    440 !
    441 !--   After the number of vertices is counted, repeat the loops and define the
    442 !--   vertices. Start with the horizontal default surfaces
    443 !--   First, however, determine the offset where couting of points should be
    444 !--   started, which is the sum of points of all PE's with lower MPI rank.
    445       i                 = 0
    446       point_index_count = 0
    447       DO WHILE ( i < myid  .AND.  i <= SIZE( num_points_on_pe ) )
    448          point_index_count = point_index_count + num_points_on_pe(i)
    449          i                 = i + 1
    450       ENDDO
    451          
    452      
    453       surfaces%npoints = 0
    454       point_index      = -1
    455       npg              = 0
    456      
    457       DO  l = 0, 1
    458          DO  m = 1, surf_def_h(0)%ns
    459 !
    460 !--         Determine the indices of the respective grid cell inside the topography
    461             i = surf_def_h(0)%i(m) + surf_def_h(0)%ioff
    462             j = surf_def_h(0)%j(m) + surf_def_h(0)%joff
    463             k = surf_def_h(0)%k(m) + surf_def_h(0)%koff
    464 !
    465 !--         Check if the vertices that define the surface element are already
    466 !--         defined, if not, increment the counter.
     578         DO  m = 1, surf_lsm_h%ns
     579            i = surf_lsm_h%i(m) + surf_lsm_h%ioff
     580            j = surf_lsm_h%j(m) + surf_lsm_h%joff
     581            k = surf_lsm_h%k(m) + surf_lsm_h%koff
    467582            IF ( point_index(k,j,i) < 0 )  THEN
    468583               surfaces%npoints   = surfaces%npoints + 1   
    469584               point_index(k,j,i) = point_index_count
    470                point_index_count  = point_index_count + 1              
     585               point_index_count  = point_index_count + 1
    471586               surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
    472587               surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
     
    503618            surfaces%polygons(3,npg)   = point_index(k,j,i+1)
    504619            surfaces%polygons(4,npg)   = point_index(k,j+1,i+1)
    505             surfaces%polygons(5,npg)   = point_index(k,j+1,i)
     620            surfaces%polygons(5,npg)   = point_index(k,j+1,i) 
    506621         ENDDO
    507       ENDDO
    508       DO  m = 1, surf_lsm_h%ns
    509          i = surf_lsm_h%i(m) + surf_lsm_h%ioff
    510          j = surf_lsm_h%j(m) + surf_lsm_h%joff
    511          k = surf_lsm_h%k(m) + surf_lsm_h%koff
    512          IF ( point_index(k,j,i) < 0 )  THEN
    513             surfaces%npoints   = surfaces%npoints + 1   
    514             point_index(k,j,i) = point_index_count
    515             point_index_count  = point_index_count + 1
    516             surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
    517             surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
    518             surfaces%points(3,surfaces%npoints) = zw(k)
    519          ENDIF
    520          IF ( point_index(k,j,i+1) < 0 )  THEN
    521             surfaces%npoints     = surfaces%npoints + 1   
    522             point_index(k,j,i+1) = point_index_count
    523             point_index_count    = point_index_count + 1
    524             surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
    525             surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
    526             surfaces%points(3,surfaces%npoints) = zw(k)
    527          ENDIF
    528          IF ( point_index(k,j+1,i+1) < 0 )  THEN
    529             surfaces%npoints       = surfaces%npoints + 1   
    530             point_index(k,j+1,i+1) = point_index_count
    531             point_index_count      = point_index_count + 1
    532             surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
    533             surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
    534             surfaces%points(3,surfaces%npoints) = zw(k)
    535          ENDIF
    536          IF ( point_index(k,j+1,i) < 0 )  THEN
    537             surfaces%npoints     = surfaces%npoints + 1   
    538             point_index(k,j+1,i) = point_index_count
    539             point_index_count    = point_index_count + 1
    540             surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
    541             surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
    542             surfaces%points(3,surfaces%npoints) = zw(k)
    543          ENDIF
    544          
    545          npg                        = npg + 1
    546          surfaces%polygons(1,npg)   = 4
    547          surfaces%polygons(2,npg)   = point_index(k,j,i)
    548          surfaces%polygons(3,npg)   = point_index(k,j,i+1)
    549          surfaces%polygons(4,npg)   = point_index(k,j+1,i+1)
    550          surfaces%polygons(5,npg)   = point_index(k,j+1,i)
    551       ENDDO
    552 
    553       DO  m = 1, surf_usm_h%ns
    554          i = surf_usm_h%i(m) + surf_usm_h%ioff
    555          j = surf_usm_h%j(m) + surf_usm_h%joff
    556          k = surf_usm_h%k(m) + surf_usm_h%koff
    557 
    558          IF ( point_index(k,j,i) < 0 )  THEN
    559             surfaces%npoints   = surfaces%npoints + 1   
    560             point_index(k,j,i) = point_index_count
    561             point_index_count  = point_index_count + 1
    562             surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
    563             surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
    564             surfaces%points(3,surfaces%npoints) = zw(k)
    565          ENDIF
    566          IF ( point_index(k,j,i+1) < 0 )  THEN
    567             surfaces%npoints     = surfaces%npoints + 1   
    568             point_index(k,j,i+1) = point_index_count
    569             point_index_count    = point_index_count + 1
    570             surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
    571             surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
    572             surfaces%points(3,surfaces%npoints) = zw(k)
    573          ENDIF
    574          IF ( point_index(k,j+1,i+1) < 0 )  THEN
    575             surfaces%npoints       = surfaces%npoints + 1   
    576             point_index(k,j+1,i+1) = point_index_count
    577             point_index_count      = point_index_count + 1
    578             surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
    579             surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
    580             surfaces%points(3,surfaces%npoints) = zw(k)
    581          ENDIF
    582          IF ( point_index(k,j+1,i) < 0 )  THEN
    583             surfaces%npoints     = surfaces%npoints + 1   
    584             point_index(k,j+1,i) = point_index_count
    585             point_index_count    = point_index_count + 1
    586             surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
    587             surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
    588             surfaces%points(3,surfaces%npoints) = zw(k)
    589          ENDIF
    590          
    591          npg                        = npg + 1
    592          surfaces%polygons(1,npg)   = 4
    593          surfaces%polygons(2,npg)   = point_index(k,j,i)
    594          surfaces%polygons(3,npg)   = point_index(k,j,i+1)
    595          surfaces%polygons(4,npg)   = point_index(k,j+1,i+1)
    596          surfaces%polygons(5,npg)   = point_index(k,j+1,i)   
    597       ENDDO
    598 
    599       DO  l = 0, 3
    600          DO  m = 1, surf_def_v(l)%ns
    601 !
    602 !--         Determine the indices of the respective grid cell inside the topography.
    603 !--         Please note, j-index for north-facing surfaces ( l==0 ) is
    604 !--         identical to the reference j-index outside the grid box.
    605 !--         Equivalent for east-facing surfaces and i-index.
    606             i = surf_def_v(l)%i(m) + MERGE( surf_def_v(l)%ioff, 0, l == 3 )
    607             j = surf_def_v(l)%j(m) + MERGE( surf_def_v(l)%joff, 0, l == 1 )
    608             k = surf_def_v(l)%k(m) + surf_def_v(l)%koff
    609 !
    610 !--         Lower left /front vertex
     622     
     623         DO  m = 1, surf_usm_h%ns
     624            i = surf_usm_h%i(m) + surf_usm_h%ioff
     625            j = surf_usm_h%j(m) + surf_usm_h%joff
     626            k = surf_usm_h%k(m) + surf_usm_h%koff
     627     
    611628            IF ( point_index(k,j,i) < 0 )  THEN
    612629               surfaces%npoints   = surfaces%npoints + 1   
     
    615632               surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
    616633               surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
    617                surfaces%points(3,surfaces%npoints) = zw(k-1)
    618             ENDIF
    619 !
    620 !--         Upper / lower right index for north- and south-facing surfaces
    621             IF ( l == 0  .OR.  l == 1 )  THEN
    622                IF ( point_index(k,j,i+1) < 0 )  THEN
    623                   surfaces%npoints     = surfaces%npoints + 1   
    624                   point_index(k,j,i+1) = point_index_count
    625                   point_index_count    = point_index_count + 1
    626                   surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
    627                   surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
    628                   surfaces%points(3,surfaces%npoints) = zw(k-1)
    629                ENDIF
    630                IF ( point_index(k+1,j,i+1) < 0 )  THEN
    631                   surfaces%npoints       = surfaces%npoints + 1   
    632                   point_index(k+1,j,i+1) = point_index_count
    633                   point_index_count      = point_index_count + 1
    634                   surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
    635                   surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
    636                   surfaces%points(3,surfaces%npoints) = zw(k)
    637                ENDIF
    638 !
    639 !--         Upper / lower front index for east- and west-facing surfaces
    640             ELSEIF ( l == 2  .OR.  l == 3 )  THEN
    641                IF ( point_index(k,j+1,i) < 0 )  THEN
    642                   surfaces%npoints     = surfaces%npoints + 1   
    643                   point_index(k,j+1,i) = point_index_count
    644                   point_index_count    = point_index_count + 1
    645                   surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
    646                   surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
    647                   surfaces%points(3,surfaces%npoints) = zw(k-1)
    648                ENDIF
    649                IF ( point_index(k+1,j+1,i) < 0 )  THEN
    650                   surfaces%npoints       = surfaces%npoints + 1   
    651                   point_index(k+1,j+1,i) = point_index_count
    652                   point_index_count      = point_index_count + 1
    653                   surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
    654                   surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
    655                   surfaces%points(3,surfaces%npoints) = zw(k)
    656                ENDIF
     634               surfaces%points(3,surfaces%npoints) = zw(k)
    657635            ENDIF
    658 !
    659 !--         Upper left / front vertex
    660             IF ( point_index(k+1,j,i) < 0 )  THEN
     636            IF ( point_index(k,j,i+1) < 0 )  THEN
    661637               surfaces%npoints     = surfaces%npoints + 1   
    662                point_index(k+1,j,i) = point_index_count
    663                point_index_count    = point_index_count + 1
    664                surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
    665                surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
     638               point_index(k,j,i+1) = point_index_count
     639               point_index_count    = point_index_count + 1
     640               surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
     641               surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
     642               surfaces%points(3,surfaces%npoints) = zw(k)
     643            ENDIF
     644            IF ( point_index(k,j+1,i+1) < 0 )  THEN
     645               surfaces%npoints       = surfaces%npoints + 1   
     646               point_index(k,j+1,i+1) = point_index_count
     647               point_index_count      = point_index_count + 1
     648               surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
     649               surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
     650               surfaces%points(3,surfaces%npoints) = zw(k)
     651            ENDIF
     652            IF ( point_index(k,j+1,i) < 0 )  THEN
     653               surfaces%npoints     = surfaces%npoints + 1   
     654               point_index(k,j+1,i) = point_index_count
     655               point_index_count    = point_index_count + 1
     656               surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
     657               surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
    666658               surfaces%points(3,surfaces%npoints) = zw(k)
    667659            ENDIF
    668660           
    669             npg = npg + 1           
    670             IF ( l == 0  .OR.  l == 1 )  THEN
    671                surfaces%polygons(1,npg)   = 4
    672                surfaces%polygons(2,npg)   = point_index(k,j,i)
    673                surfaces%polygons(3,npg)   = point_index(k,j,i+1)
    674                surfaces%polygons(4,npg)   = point_index(k+1,j,i+1)
    675                surfaces%polygons(5,npg)   = point_index(k+1,j,i) 
    676             ELSE
    677                surfaces%polygons(1,npg)   = 4
    678                surfaces%polygons(2,npg)   = point_index(k,j,i)
    679                surfaces%polygons(3,npg)   = point_index(k,j+1,i)
    680                surfaces%polygons(4,npg)   = point_index(k+1,j+1,i)
    681                surfaces%polygons(5,npg)   = point_index(k+1,j,i)
    682             ENDIF
    683            
     661            npg                        = npg + 1
     662            surfaces%polygons(1,npg)   = 4
     663            surfaces%polygons(2,npg)   = point_index(k,j,i)
     664            surfaces%polygons(3,npg)   = point_index(k,j,i+1)
     665            surfaces%polygons(4,npg)   = point_index(k,j+1,i+1)
     666            surfaces%polygons(5,npg)   = point_index(k,j+1,i)   
    684667         ENDDO
    685668
    686          DO  m = 1, surf_lsm_v(l)%ns
    687             i = surf_lsm_v(l)%i(m) + MERGE( surf_lsm_v(l)%ioff, 0, l == 3 )
    688             j = surf_lsm_v(l)%j(m) + MERGE( surf_lsm_v(l)%joff, 0, l == 1 )
    689             k = surf_lsm_v(l)%k(m) + surf_lsm_v(l)%koff
    690 !
    691 !--         Lower left /front vertex
    692             IF ( point_index(k,j,i) < 0 )  THEN
    693                surfaces%npoints   = surfaces%npoints + 1   
    694                point_index(k,j,i) = point_index_count
    695                point_index_count  = point_index_count + 1
    696                surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
    697                surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
    698                surfaces%points(3,surfaces%npoints) = zw(k-1)
    699             ENDIF
    700 !
    701 !--         Upper / lower right index for north- and south-facing surfaces
    702             IF ( l == 0  .OR.  l == 1 )  THEN
    703                IF ( point_index(k,j,i+1) < 0 )  THEN
    704                   surfaces%npoints     = surfaces%npoints + 1   
    705                   point_index(k,j,i+1) = point_index_count
    706                   point_index_count    = point_index_count + 1
    707                   surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
    708                   surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
     669         DO  l = 0, 3
     670            DO  m = 1, surf_def_v(l)%ns
     671!       
     672!--            Determine the indices of the respective grid cell inside the topography.
     673!--            Please note, j-index for north-facing surfaces ( l==0 ) is
     674!--            identical to the reference j-index outside the grid box.
     675!--            Equivalent for east-facing surfaces and i-index.
     676               i = surf_def_v(l)%i(m) + MERGE( surf_def_v(l)%ioff, 0, l == 3 )
     677               j = surf_def_v(l)%j(m) + MERGE( surf_def_v(l)%joff, 0, l == 1 )
     678               k = surf_def_v(l)%k(m) + surf_def_v(l)%koff
     679!       
     680!--            Lower left /front vertex
     681               IF ( point_index(k,j,i) < 0 )  THEN
     682                  surfaces%npoints   = surfaces%npoints + 1   
     683                  point_index(k,j,i) = point_index_count
     684                  point_index_count  = point_index_count + 1
     685                  surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
     686                  surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
    709687                  surfaces%points(3,surfaces%npoints) = zw(k-1)
    710688               ENDIF
    711                IF ( point_index(k+1,j,i+1) < 0 )  THEN
    712                   surfaces%npoints       = surfaces%npoints + 1   
    713                   point_index(k+1,j,i+1) = point_index_count
    714                   point_index_count      = point_index_count + 1
    715                   surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
    716                   surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
     689!       
     690!--            Upper / lower right index for north- and south-facing surfaces
     691               IF ( l == 0  .OR.  l == 1 )  THEN
     692                  IF ( point_index(k,j,i+1) < 0 )  THEN
     693                     surfaces%npoints     = surfaces%npoints + 1   
     694                     point_index(k,j,i+1) = point_index_count
     695                     point_index_count    = point_index_count + 1
     696                     surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
     697                     surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
     698                     surfaces%points(3,surfaces%npoints) = zw(k-1)
     699                  ENDIF
     700                  IF ( point_index(k+1,j,i+1) < 0 )  THEN
     701                     surfaces%npoints       = surfaces%npoints + 1   
     702                     point_index(k+1,j,i+1) = point_index_count
     703                     point_index_count      = point_index_count + 1
     704                     surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
     705                     surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
     706                     surfaces%points(3,surfaces%npoints) = zw(k)
     707                  ENDIF
     708!       
     709!--            Upper / lower front index for east- and west-facing surfaces
     710               ELSEIF ( l == 2  .OR.  l == 3 )  THEN
     711                  IF ( point_index(k,j+1,i) < 0 )  THEN
     712                     surfaces%npoints     = surfaces%npoints + 1   
     713                     point_index(k,j+1,i) = point_index_count
     714                     point_index_count    = point_index_count + 1
     715                     surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
     716                     surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
     717                     surfaces%points(3,surfaces%npoints) = zw(k-1)
     718                  ENDIF
     719                  IF ( point_index(k+1,j+1,i) < 0 )  THEN
     720                     surfaces%npoints       = surfaces%npoints + 1   
     721                     point_index(k+1,j+1,i) = point_index_count
     722                     point_index_count      = point_index_count + 1
     723                     surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
     724                     surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
     725                     surfaces%points(3,surfaces%npoints) = zw(k)
     726                  ENDIF
     727               ENDIF
     728!       
     729!--            Upper left / front vertex
     730               IF ( point_index(k+1,j,i) < 0 )  THEN
     731                  surfaces%npoints     = surfaces%npoints + 1   
     732                  point_index(k+1,j,i) = point_index_count
     733                  point_index_count    = point_index_count + 1
     734                  surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
     735                  surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
    717736                  surfaces%points(3,surfaces%npoints) = zw(k)
    718737               ENDIF
    719 !
    720 !--         Upper / lower front index for east- and west-facing surfaces
    721             ELSEIF ( l == 2  .OR.  l == 3 )  THEN
    722                IF ( point_index(k,j+1,i) < 0 )  THEN
    723                   surfaces%npoints     = surfaces%npoints + 1   
    724                   point_index(k,j+1,i) = point_index_count
    725                   point_index_count    = point_index_count + 1
    726                   surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
    727                   surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
     738               
     739               npg = npg + 1           
     740               IF ( l == 0  .OR.  l == 1 )  THEN
     741                  surfaces%polygons(1,npg)   = 4
     742                  surfaces%polygons(2,npg)   = point_index(k,j,i)
     743                  surfaces%polygons(3,npg)   = point_index(k,j,i+1)
     744                  surfaces%polygons(4,npg)   = point_index(k+1,j,i+1)
     745                  surfaces%polygons(5,npg)   = point_index(k+1,j,i) 
     746               ELSE
     747                  surfaces%polygons(1,npg)   = 4
     748                  surfaces%polygons(2,npg)   = point_index(k,j,i)
     749                  surfaces%polygons(3,npg)   = point_index(k,j+1,i)
     750                  surfaces%polygons(4,npg)   = point_index(k+1,j+1,i)
     751                  surfaces%polygons(5,npg)   = point_index(k+1,j,i)
     752               ENDIF
     753               
     754            ENDDO
     755         
     756            DO  m = 1, surf_lsm_v(l)%ns
     757               i = surf_lsm_v(l)%i(m) + MERGE( surf_lsm_v(l)%ioff, 0, l == 3 )
     758               j = surf_lsm_v(l)%j(m) + MERGE( surf_lsm_v(l)%joff, 0, l == 1 )
     759               k = surf_lsm_v(l)%k(m) + surf_lsm_v(l)%koff
     760!       
     761!--            Lower left /front vertex
     762               IF ( point_index(k,j,i) < 0 )  THEN
     763                  surfaces%npoints   = surfaces%npoints + 1   
     764                  point_index(k,j,i) = point_index_count
     765                  point_index_count  = point_index_count + 1
     766                  surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
     767                  surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
    728768                  surfaces%points(3,surfaces%npoints) = zw(k-1)
    729769               ENDIF
    730                IF ( point_index(k+1,j+1,i) < 0 )  THEN
    731                   surfaces%npoints       = surfaces%npoints + 1   
    732                   point_index(k+1,j+1,i) = point_index_count
    733                   point_index_count      = point_index_count + 1
    734                   surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
    735                   surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
     770!       
     771!--            Upper / lower right index for north- and south-facing surfaces
     772               IF ( l == 0  .OR.  l == 1 )  THEN
     773                  IF ( point_index(k,j,i+1) < 0 )  THEN
     774                     surfaces%npoints     = surfaces%npoints + 1   
     775                     point_index(k,j,i+1) = point_index_count
     776                     point_index_count    = point_index_count + 1
     777                     surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
     778                     surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
     779                     surfaces%points(3,surfaces%npoints) = zw(k-1)
     780                  ENDIF
     781                  IF ( point_index(k+1,j,i+1) < 0 )  THEN
     782                     surfaces%npoints       = surfaces%npoints + 1   
     783                     point_index(k+1,j,i+1) = point_index_count
     784                     point_index_count      = point_index_count + 1
     785                     surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
     786                     surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
     787                     surfaces%points(3,surfaces%npoints) = zw(k)
     788                  ENDIF
     789!       
     790!--            Upper / lower front index for east- and west-facing surfaces
     791               ELSEIF ( l == 2  .OR.  l == 3 )  THEN
     792                  IF ( point_index(k,j+1,i) < 0 )  THEN
     793                     surfaces%npoints     = surfaces%npoints + 1   
     794                     point_index(k,j+1,i) = point_index_count
     795                     point_index_count    = point_index_count + 1
     796                     surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
     797                     surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
     798                     surfaces%points(3,surfaces%npoints) = zw(k-1)
     799                  ENDIF
     800                  IF ( point_index(k+1,j+1,i) < 0 )  THEN
     801                     surfaces%npoints       = surfaces%npoints + 1   
     802                     point_index(k+1,j+1,i) = point_index_count
     803                     point_index_count      = point_index_count + 1
     804                     surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
     805                     surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
     806                     surfaces%points(3,surfaces%npoints) = zw(k)
     807                  ENDIF
     808               ENDIF
     809!       
     810!--            Upper left / front vertex
     811               IF ( point_index(k+1,j,i) < 0 )  THEN
     812                  surfaces%npoints     = surfaces%npoints + 1   
     813                  point_index(k+1,j,i) = point_index_count
     814                  point_index_count    = point_index_count + 1
     815                  surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
     816                  surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
    736817                  surfaces%points(3,surfaces%npoints) = zw(k)
    737818               ENDIF
    738             ENDIF
    739 !
    740 !--         Upper left / front vertex
    741             IF ( point_index(k+1,j,i) < 0 )  THEN
    742                surfaces%npoints     = surfaces%npoints + 1   
    743                point_index(k+1,j,i) = point_index_count
    744                point_index_count    = point_index_count + 1
    745                surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
    746                surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
    747                surfaces%points(3,surfaces%npoints) = zw(k)
    748             ENDIF
    749            
    750             npg = npg + 1           
    751             IF ( l == 0  .OR.  l == 1 )  THEN
    752                surfaces%polygons(1,npg)   = 4
    753                surfaces%polygons(2,npg)   = point_index(k,j,i)
    754                surfaces%polygons(3,npg)   = point_index(k,j,i+1)
    755                surfaces%polygons(4,npg)   = point_index(k+1,j,i+1)
    756                surfaces%polygons(5,npg)   = point_index(k+1,j,i) 
    757             ELSE
    758                surfaces%polygons(1,npg)   = 4
    759                surfaces%polygons(2,npg)   = point_index(k,j,i)
    760                surfaces%polygons(3,npg)   = point_index(k,j+1,i)
    761                surfaces%polygons(4,npg)   = point_index(k+1,j+1,i)
    762                surfaces%polygons(5,npg)   = point_index(k+1,j,i)
    763             ENDIF
    764          ENDDO
    765          DO  m = 1, surf_usm_v(l)%ns
    766             i = surf_usm_v(l)%i(m) + MERGE( surf_usm_v(l)%ioff, 0, l == 3 )
    767             j = surf_usm_v(l)%j(m) + MERGE( surf_usm_v(l)%joff, 0, l == 1 )
    768             k = surf_usm_v(l)%k(m) + surf_usm_v(l)%koff
    769 !
    770 !--         Lower left /front vertex
    771             IF ( point_index(k,j,i) < 0 )  THEN
    772                surfaces%npoints   = surfaces%npoints + 1   
    773                point_index(k,j,i) = point_index_count
    774                point_index_count  = point_index_count + 1
    775                surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
    776                surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
    777                surfaces%points(3,surfaces%npoints) = zw(k-1)
    778             ENDIF
    779 !
    780 !--         Upper / lower right index for north- and south-facing surfaces
    781             IF ( l == 0  .OR.  l == 1 )  THEN
    782                IF ( point_index(k,j,i+1) < 0 )  THEN
    783                   surfaces%npoints     = surfaces%npoints + 1   
    784                   point_index(k,j,i+1) = point_index_count
    785                   point_index_count    = point_index_count + 1
    786                   surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
    787                   surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
     819               
     820               npg = npg + 1           
     821               IF ( l == 0  .OR.  l == 1 )  THEN
     822                  surfaces%polygons(1,npg)   = 4
     823                  surfaces%polygons(2,npg)   = point_index(k,j,i)
     824                  surfaces%polygons(3,npg)   = point_index(k,j,i+1)
     825                  surfaces%polygons(4,npg)   = point_index(k+1,j,i+1)
     826                  surfaces%polygons(5,npg)   = point_index(k+1,j,i) 
     827               ELSE
     828                  surfaces%polygons(1,npg)   = 4
     829                  surfaces%polygons(2,npg)   = point_index(k,j,i)
     830                  surfaces%polygons(3,npg)   = point_index(k,j+1,i)
     831                  surfaces%polygons(4,npg)   = point_index(k+1,j+1,i)
     832                  surfaces%polygons(5,npg)   = point_index(k+1,j,i)
     833               ENDIF
     834            ENDDO
     835            DO  m = 1, surf_usm_v(l)%ns
     836               i = surf_usm_v(l)%i(m) + MERGE( surf_usm_v(l)%ioff, 0, l == 3 )
     837               j = surf_usm_v(l)%j(m) + MERGE( surf_usm_v(l)%joff, 0, l == 1 )
     838               k = surf_usm_v(l)%k(m) + surf_usm_v(l)%koff
     839!       
     840!--            Lower left /front vertex
     841               IF ( point_index(k,j,i) < 0 )  THEN
     842                  surfaces%npoints   = surfaces%npoints + 1   
     843                  point_index(k,j,i) = point_index_count
     844                  point_index_count  = point_index_count + 1
     845                  surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
     846                  surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
    788847                  surfaces%points(3,surfaces%npoints) = zw(k-1)
    789848               ENDIF
    790                IF ( point_index(k+1,j,i+1) < 0 )  THEN
    791                   surfaces%npoints       = surfaces%npoints + 1   
    792                   point_index(k+1,j,i+1) = point_index_count
    793                   point_index_count      = point_index_count + 1
    794                   surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
    795                   surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
     849!       
     850!--            Upper / lower right index for north- and south-facing surfaces
     851               IF ( l == 0  .OR.  l == 1 )  THEN
     852                  IF ( point_index(k,j,i+1) < 0 )  THEN
     853                     surfaces%npoints     = surfaces%npoints + 1   
     854                     point_index(k,j,i+1) = point_index_count
     855                     point_index_count    = point_index_count + 1
     856                     surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
     857                     surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
     858                     surfaces%points(3,surfaces%npoints) = zw(k-1)
     859                  ENDIF
     860                  IF ( point_index(k+1,j,i+1) < 0 )  THEN
     861                     surfaces%npoints       = surfaces%npoints + 1   
     862                     point_index(k+1,j,i+1) = point_index_count
     863                     point_index_count      = point_index_count + 1
     864                     surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
     865                     surfaces%points(2,surfaces%npoints) = ( j     - 0.5_wp ) * dy
     866                     surfaces%points(3,surfaces%npoints) = zw(k)
     867                  ENDIF
     868!       
     869!--            Upper / lower front index for east- and west-facing surfaces
     870               ELSEIF ( l == 2  .OR.  l == 3 )  THEN
     871                  IF ( point_index(k,j+1,i) < 0 )  THEN
     872                     surfaces%npoints     = surfaces%npoints + 1   
     873                     point_index(k,j+1,i) = point_index_count
     874                     point_index_count    = point_index_count + 1
     875                     surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
     876                     surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
     877                     surfaces%points(3,surfaces%npoints) = zw(k-1)
     878                  ENDIF
     879                  IF ( point_index(k+1,j+1,i) < 0 )  THEN
     880                     surfaces%npoints       = surfaces%npoints + 1   
     881                     point_index(k+1,j+1,i) = point_index_count
     882                     point_index_count      = point_index_count + 1
     883                     surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
     884                     surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
     885                     surfaces%points(3,surfaces%npoints) = zw(k)
     886                  ENDIF
     887               ENDIF
     888!       
     889!--            Upper left / front vertex
     890               IF ( point_index(k+1,j,i) < 0 )  THEN
     891                  surfaces%npoints     = surfaces%npoints + 1   
     892                  point_index(k+1,j,i) = point_index_count
     893                  point_index_count    = point_index_count + 1
     894                  surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
     895                  surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
    796896                  surfaces%points(3,surfaces%npoints) = zw(k)
    797897               ENDIF
    798 !
    799 !--         Upper / lower front index for east- and west-facing surfaces
    800             ELSEIF ( l == 2  .OR.  l == 3 )  THEN
    801                IF ( point_index(k,j+1,i) < 0 )  THEN
    802                   surfaces%npoints     = surfaces%npoints + 1   
    803                   point_index(k,j+1,i) = point_index_count
    804                   point_index_count    = point_index_count + 1
    805                   surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
    806                   surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
    807                   surfaces%points(3,surfaces%npoints) = zw(k-1)
    808                ENDIF
    809                IF ( point_index(k+1,j+1,i) < 0 )  THEN
    810                   surfaces%npoints       = surfaces%npoints + 1   
    811                   point_index(k+1,j+1,i) = point_index_count
    812                   point_index_count      = point_index_count + 1
    813                   surfaces%points(1,surfaces%npoints) = ( i     - 0.5_wp ) * dx
    814                   surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
    815                   surfaces%points(3,surfaces%npoints) = zw(k)
    816                ENDIF
    817             ENDIF
    818 !
    819 !--         Upper left / front vertex
    820             IF ( point_index(k+1,j,i) < 0 )  THEN
    821                surfaces%npoints     = surfaces%npoints + 1   
    822                point_index(k+1,j,i) = point_index_count
    823                point_index_count    = point_index_count + 1
    824                surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
    825                surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
    826                surfaces%points(3,surfaces%npoints) = zw(k)
    827             ENDIF
    828            
    829             npg = npg + 1           
    830             IF ( l == 0  .OR.  l == 1 )  THEN
    831                surfaces%polygons(1,npg)   = 4
    832                surfaces%polygons(2,npg)   = point_index(k,j,i)
    833                surfaces%polygons(3,npg)   = point_index(k,j,i+1)
    834                surfaces%polygons(4,npg)   = point_index(k+1,j,i+1)
    835                surfaces%polygons(5,npg)   = point_index(k+1,j,i) 
    836             ELSE
    837                surfaces%polygons(1,npg)   = 4
    838                surfaces%polygons(2,npg)   = point_index(k,j,i)
    839                surfaces%polygons(3,npg)   = point_index(k,j+1,i)
    840                surfaces%polygons(4,npg)   = point_index(k+1,j+1,i)
    841                surfaces%polygons(5,npg)   = point_index(k+1,j,i)
    842             ENDIF
     898               
     899               npg = npg + 1           
     900               IF ( l == 0  .OR.  l == 1 )  THEN
     901                  surfaces%polygons(1,npg)   = 4
     902                  surfaces%polygons(2,npg)   = point_index(k,j,i)
     903                  surfaces%polygons(3,npg)   = point_index(k,j,i+1)
     904                  surfaces%polygons(4,npg)   = point_index(k+1,j,i+1)
     905                  surfaces%polygons(5,npg)   = point_index(k+1,j,i) 
     906               ELSE
     907                  surfaces%polygons(1,npg)   = 4
     908                  surfaces%polygons(2,npg)   = point_index(k,j,i)
     909                  surfaces%polygons(3,npg)   = point_index(k,j+1,i)
     910                  surfaces%polygons(4,npg)   = point_index(k+1,j+1,i)
     911                  surfaces%polygons(5,npg)   = point_index(k+1,j,i)
     912               ENDIF
     913            ENDDO
     914         
    843915         ENDDO
    844 
    845       ENDDO
    846 !
    847 !--   Deallocate temporary dummy variable
    848       DEALLOCATE ( point_index )
    849 !
    850 !--   Sum-up total number of surface elements and vertices on domain. This
    851 !--   will be needed for post-processing.
    852       surfaces%npoints_total = 0
     916!       
     917!--      Deallocate temporary dummy variable
     918         DEALLOCATE ( point_index )
     919!
     920!--      Sum-up total number of vertices on domain. This
     921!--      will be needed for post-processing.
     922         surfaces%npoints_total = 0
    853923#if defined( __parallel )
    854        CALL MPI_ALLREDUCE( surfaces%npoints, surfaces%npoints_total, 1,        &
    855                            MPI_INTEGER, MPI_SUM, comm2d, ierr )
    856        CALL MPI_ALLREDUCE( surfaces%ns, surfaces%ns_total, 1,                  &
    857                            MPI_INTEGER, MPI_SUM, comm2d, ierr )
     924          CALL MPI_ALLREDUCE( surfaces%npoints, surfaces%npoints_total, 1,     &
     925                              MPI_INTEGER, MPI_SUM, comm2d, ierr )
    858926#else
    859        surfaces%npoints_total = surfaces%npoints
    860        surfaces%ns_total      = surfaces%ns
     927          surfaces%npoints_total = surfaces%npoints
    861928#endif
     929       ENDIF
     930!
     931!--    If output to netcdf is enabled, set-up the coordinate arrays that
     932!--    unambiguously describe the position and orientation of each surface
     933!--    element.
     934       IF ( to_netcdf )  THEN
     935!
     936!--       Allocate local coordinate arrays
     937          ALLOCATE( surfaces%s(1:surfaces%ns)       )
     938          ALLOCATE( surfaces%xs(1:surfaces%ns)      )
     939          ALLOCATE( surfaces%ys(1:surfaces%ns)      )
     940          ALLOCATE( surfaces%zs(1:surfaces%ns)      )
     941          ALLOCATE( surfaces%azimuth(1:surfaces%ns) )
     942          ALLOCATE( surfaces%zenith(1:surfaces%ns)  )
     943          ALLOCATE( surfaces%es_utm(1:surfaces%ns)  )
     944          ALLOCATE( surfaces%ns_utm(1:surfaces%ns)  )
     945!
     946!--       Gather the number of surface on each processor, in order to number
     947!--       the surface elements in ascending order with respect to the total
     948!--       number of surfaces in the domain.
     949#if defined( __parallel )
     950          CALL MPI_ALLGATHER( surfaces%ns, 1, MPI_INTEGER,                     &
     951                              num_surfaces_on_pe, 1, MPI_INTEGER, comm2d, ierr  )
     952#else
     953          num_surfaces_on_pe = surfaces%ns
     954#endif
     955!
     956!--       First, however, determine the offset where couting of the surfaces
     957!--       should start (the sum of surfaces on all PE's with lower MPI rank).
     958          i           = 0
     959          start_count = 1
     960          DO WHILE ( i < myid  .AND.  i <= SIZE( num_surfaces_on_pe ) )
     961             start_count = start_count + num_surfaces_on_pe(i)
     962             i           = i + 1
     963          ENDDO
     964!       
     965!--       Set coordinate arrays. For horizontal surfaces, azimuth
     966!--       angles are not defined (fill value). Zenith angle is 0 (180) for
     967!--       upward (downward)-facing surfaces.
     968          i  = start_count
     969          mm = 1
     970          DO  m = 1, surf_def_h(0)%ns
     971             surfaces%s(mm)       = i
     972             surfaces%xs(mm)      = ( surf_def_h(0)%i(m) + 0.5_wp ) * dx
     973             surfaces%ys(mm)      = ( surf_def_h(0)%j(m) + 0.5_wp ) * dy
     974             surfaces%zs(mm)      = zw(surf_def_h(0)%k(m)+surf_def_h(0)%koff)
     975             surfaces%azimuth(mm) = surfaces%fillvalue
     976             surfaces%zenith(mm)  = 0.0
     977             i                    = i + 1
     978             mm                   = mm + 1
     979          ENDDO
     980          DO  m = 1, surf_lsm_h%ns
     981             surfaces%s(mm)       = i
     982             surfaces%xs(mm)      = ( surf_lsm_h%i(m) + 0.5_wp ) * dx
     983             surfaces%ys(mm)      = ( surf_lsm_h%j(m) + 0.5_wp ) * dy
     984             surfaces%zs(mm)      = zw(surf_lsm_h%k(m)+surf_lsm_h%koff)
     985             surfaces%azimuth(mm) = surfaces%fillvalue
     986             surfaces%zenith(mm)  = 0.0
     987             i                    = i + 1
     988             mm                   = mm + 1
     989          ENDDO
     990          DO  m = 1, surf_usm_h%ns
     991             surfaces%s(mm)       = i
     992             surfaces%xs(mm)      = ( surf_usm_h%i(m) + 0.5_wp ) * dx
     993             surfaces%ys(mm)      = ( surf_usm_h%j(m) + 0.5_wp ) * dy
     994             surfaces%zs(mm)      = zw(surf_usm_h%k(m)+surf_usm_h%koff)
     995             surfaces%azimuth(mm) = surfaces%fillvalue
     996             surfaces%zenith(mm)  = 0.0
     997             i                    = i + 1
     998             mm                   = mm + 1
     999          ENDDO
     1000          DO  m = 1, surf_def_h(1)%ns
     1001             surfaces%s(mm)       = i
     1002             surfaces%xs(mm)      = ( surf_def_h(1)%i(m) + 0.5_wp ) * dx
     1003             surfaces%ys(mm)      = ( surf_def_h(1)%j(m) + 0.5_wp ) * dy
     1004             surfaces%zs(mm)      = zw(surf_def_h(1)%k(m)+surf_def_h(1)%koff)
     1005             surfaces%azimuth(mm) = surfaces%fillvalue
     1006             surfaces%zenith(mm)  = 180.0
     1007             i                    = i + 1
     1008             mm                   = mm + 1
     1009          ENDDO
     1010!       
     1011!--       For vertical surfaces, zenith angles are not defined (fill value). 
     1012!--       Azimuth angle: eastward (0), westward (180), northward (270),
     1013!--       southward (90).
     1014          DO  l = 0, 3
     1015             IF ( l == 0 )  THEN
     1016                az    = 270.0_wp
     1017                off_x =   0.0_wp
     1018                off_y = - 0.5_wp
     1019             ELSEIF ( l == 1 )  THEN
     1020                az    = 90.0_wp
     1021                off_x =  0.0_wp
     1022                off_y =  0.5_wp
     1023             ELSEIF ( l == 2 )  THEN
     1024                az    =   0.0_wp
     1025                off_x = - 0.5_wp
     1026                off_y =   0.0_wp
     1027             ELSEIF ( l == 3 )  THEN
     1028                az   = 180.0_wp
     1029                off_x =  0.5_wp
     1030                off_y =  0.0_wp
     1031             ENDIF
     1032               
     1033             DO  m = 1, surf_def_v(l)%ns
     1034                surfaces%s(mm)       = i
     1035                surfaces%xs(mm)      = ( surf_def_v(l)%i(m) + off_x ) * dx
     1036                surfaces%ys(mm)      = ( surf_def_v(l)%j(m) + off_y ) * dy
     1037                surfaces%zs(mm)      = zu(surf_def_v(l)%k(m))
     1038                surfaces%azimuth(mm) = az
     1039                surfaces%zenith(mm)  = surfaces%fillvalue
     1040                i                    = i + 1
     1041                mm                   = mm + 1
     1042             ENDDO
     1043             DO  m = 1, surf_lsm_v(l)%ns
     1044                surfaces%s(mm)       = i
     1045                surfaces%xs(mm)      = ( surf_lsm_v(l)%i(m) + off_x ) * dx
     1046                surfaces%ys(mm)      = ( surf_lsm_v(l)%j(m) + off_y ) * dy
     1047                surfaces%zs(mm)      = zu(surf_lsm_v(l)%k(m))
     1048                surfaces%azimuth(mm) = az
     1049                surfaces%zenith(mm)  = surfaces%fillvalue
     1050                i                    = i + 1
     1051                mm                   = mm + 1
     1052             ENDDO
     1053             DO  m = 1, surf_usm_v(l)%ns
     1054                surfaces%s(mm)       = i
     1055                surfaces%xs(mm)      = ( surf_usm_v(l)%i(m) + off_x ) * dx
     1056                surfaces%ys(mm)      = ( surf_usm_v(l)%j(m) + off_y ) * dy
     1057                surfaces%zs(mm)      = zu(surf_usm_v(l)%k(m))
     1058                surfaces%azimuth(mm) = az
     1059                surfaces%zenith(mm)  = surfaces%fillvalue
     1060                i                    = i + 1
     1061                mm                   = mm + 1
     1062             ENDDO
     1063          ENDDO
     1064!       
     1065!--       Finally, define UTM coordinates, which are the x/y-coordinates
     1066!--       plus the origin (lower-left coordinate of the model domain).
     1067          surfaces%es_utm = surfaces%xs + init_model%origin_x
     1068          surfaces%ns_utm = surfaces%ys + init_model%origin_y
     1069!
     1070!--       Initialize NetCDF data output. Please note, local start position for
     1071!--       the surface elements in the NetCDF file is surfaces%s(1), while
     1072!--       the number of surfaces on the subdomain is given by surfaces%ns.
     1073#if defined( __netcdf4_parallel )
     1074
     1075!
     1076!--       Calculate number of time steps to be output
     1077          ntdim_surf(0) = dosurf_time_count(0) + CEILING(                            &
     1078                        ( end_time - MAX(                                            &
     1079                            MERGE(skip_time_dosurf, skip_time_dosurf + spinup_time,  &
     1080                                  data_output_during_spinup ),                       &
     1081                            simulated_time_at_begin )                                &
     1082                        ) / dt_dosurf )
     1083
     1084          ntdim_surf(1) = dosurf_time_count(1) + CEILING(                      &
     1085                        ( end_time - MAX(                                      &
     1086                            MERGE(   skip_time_dosurf_av, skip_time_dosurf_av  &
     1087                                  + spinup_time, data_output_during_spinup ),  &
     1088                            simulated_time_at_begin )                          &
     1089                        ) / dt_dosurf_av )
     1090
     1091!
     1092!--       Create NetCDF4 files for parallel writing
     1093          DO  av = 0, 1
     1094!
     1095!--          If there is no instantaneous data (av=0) or averaged data (av=1)
     1096!--          requested, do not create the corresponding NetCDF file
     1097             IF ( dosurf_no(av) == 0 ) CYCLE
     1098
     1099             IF ( av == 0 )  THEN
     1100                filename = 'SURFACE_DATA_NETCDF' // TRIM( coupling_char )
     1101             ELSE
     1102                filename = 'SURFACE_DATA_AV_NETCDF' // TRIM( coupling_char )
     1103             ENDIF
     1104!
     1105!--          Open file using netCDF4/HDF5 format, parallel
     1106             nc_stat = NF90_CREATE( TRIM(filename),                            &
     1107                                    IOR( NF90_NOCLOBBER,                       &
     1108                                    IOR( NF90_NETCDF4, NF90_MPIIO ) ),         &
     1109                                    id_set_surf(av), COMM = comm2d, INFO = MPI_INFO_NULL )
     1110             CALL netcdf_handle_error( 'surface_data_output_mod', 5550 )
     1111
     1112             !- Write some global attributes
     1113             IF ( av == 0 )  THEN
     1114                CALL netcdf_create_global_atts( id_set_surf(av), 'surface-data', TRIM( run_description_header ), 5551 )
     1115                time_average_text = ' '
     1116             ELSE
     1117                CALL netcdf_create_global_atts( id_set_surf(av), 'surface-data_av', TRIM( run_description_header ), 5552 )
     1118                WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval_surf
     1119                nc_stat = NF90_PUT_ATT( id_set_surf(av), NF90_GLOBAL, 'time_avg',  &
     1120                                        TRIM( time_average_text ) )
     1121                CALL netcdf_handle_error( 'surface_data_output_mod', 5553 )
     1122             ENDIF
     1123
     1124
     1125!
     1126!--          Define time coordinate for surface data.
     1127!--          For parallel output the time dimension has to be limited (ntdim_surf),
     1128!--          otherwise the performance drops significantly.
     1129             CALL netcdf_create_dim( id_set_surf(av), 'time', ntdim_surf(av),  &
     1130                                     id_dim_time_surf(av), 5554 )
     1131
     1132             CALL netcdf_create_var( id_set_surf(av), (/ id_dim_time_surf(av) /),    &
     1133                                     'time', NF90_DOUBLE, id_var_time_surf(av),      &
     1134                                     'seconds since '//TRIM(init_model%origin_time), &
     1135                                     'time', 5555, 5555, 5555 )
     1136             CALL netcdf_create_att( id_set_surf(av), id_var_time_surf(av), 'standard_name', 'time', 5556)
     1137             CALL netcdf_create_att( id_set_surf(av), id_var_time_surf(av), 'axis', 'T', 5557)
     1138!
     1139!--          Define spatial dimensions and coordinates:
     1140!--          Define index of surface element
     1141             CALL netcdf_create_dim( id_set_surf(av), 's', surfaces%ns_total,  &
     1142                                     id_dim_s_surf(av), 5558 )
     1143             CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), &
     1144                                     's', NF90_DOUBLE, id_var_s_surf(av),      &
     1145                                     '1', '', 5559, 5559, 5559 )
     1146!
     1147!--          Define x coordinate
     1148             CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), &
     1149                                     'xs', NF90_DOUBLE, id_var_xs_surf(av),    &
     1150                                     'meters', '', 5561, 5561, 5561 )
     1151!
     1152!--           Define y coordinate
     1153             CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), &
     1154                                     'ys', NF90_DOUBLE, id_var_ys_surf(av),    &
     1155                                     'meters', '', 5562, 5562, 5562 )
     1156!
     1157!--          Define z coordinate
     1158             CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), &
     1159                                     'zs', NF90_DOUBLE, id_var_zs_surf(av),    &
     1160                                     'meters', '', 5560, 5560, 5560 )
     1161
     1162!
     1163!--          Define UTM coordinates
     1164             CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /),    &
     1165                                     'Es_UTM', NF90_DOUBLE, id_var_etum_surf(av), &
     1166                                     'meters', '', 5563, 5563, 5563 )
     1167             CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /),    &
     1168                                     'Ns_UTM', NF90_DOUBLE, id_var_nutm_surf(av), &
     1169                                     'meters', '', 5564, 5564, 5564 )
     1170!
     1171!--          Define the variables
     1172             var_list = ';'
     1173             i = 1
     1174
     1175             DO WHILE ( dosurf(av,i)(1:1) /= ' ' )
     1176
     1177                CALL netcdf_create_var( id_set_surf(av),(/ id_dim_s_surf(av),  &
     1178                                        id_dim_time_surf(av) /), dosurf(av,i), &
     1179                                        NF90_REAL4, id_var_dosurf(av,i),       &
     1180                                        dosurf_unit(av,i), dosurf(av,i), 5565, &
     1181                                        5565, 5565, .TRUE. )
     1182!
     1183!--                Set no fill for every variable to increase performance.
     1184                nc_stat = NF90_DEF_VAR_FILL( id_set_surf(av),     &
     1185                                             id_var_dosurf(av,i), &
     1186                                             1, 0 )
     1187                CALL netcdf_handle_error( 'surface_data_output_init', 5566 )
     1188!
     1189!--                Set collective io operations for parallel io
     1190                nc_stat = NF90_VAR_PAR_ACCESS( id_set_surf(av),     &
     1191                                               id_var_dosurf(av,i), &
     1192                                               NF90_COLLECTIVE )
     1193                CALL netcdf_handle_error( 'surface_data_output_init', 5567 )
     1194                var_list = TRIM( var_list ) // TRIM( dosurf(av,i) ) // ';'
     1195
     1196                i = i + 1
     1197
     1198             ENDDO
     1199!
     1200!--          Write the list of variables as global attribute (this is used by
     1201!--          restart runs and by combine_plot_fields)
     1202             nc_stat = NF90_PUT_ATT( id_set_surf(av), NF90_GLOBAL, 'VAR_LIST', &
     1203                                     var_list )
     1204             CALL netcdf_handle_error( 'surface_data_output_init', 5568 )
     1205
     1206!
     1207!--          Set general no fill, otherwise the performance drops significantly for
     1208!--          parallel output.
     1209             nc_stat = NF90_SET_FILL( id_set_surf(av), NF90_NOFILL, oldmode )
     1210             CALL netcdf_handle_error( 'surface_data_output_init', 5569 )
     1211
     1212!
     1213!--          Leave netCDF define mode
     1214             nc_stat = NF90_ENDDEF( id_set_surf(av) )
     1215             CALL netcdf_handle_error( 'surface_data_output_init', 5570 )
     1216
     1217!
     1218!--          These data are only written by PE0 for parallel output to increase
     1219!--          the performance.
     1220             IF ( myid == 0 )  THEN
     1221!
     1222!--             Write data for surface indices
     1223                ALLOCATE( netcdf_data_1d(1:surfaces%ns_total) )
     1224
     1225                DO  i = 1, surfaces%ns_total
     1226                   netcdf_data_1d(i) = i
     1227                ENDDO
     1228
     1229                nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_s_surf(av),  &
     1230                                        netcdf_data_1d, start = (/ 1 /),     &
     1231                                        count = (/ surfaces%ns_total /) )
     1232                CALL netcdf_handle_error( 'surface_data_output_init', 5571 )
     1233
     1234                DEALLOCATE( netcdf_data_1d )
     1235
     1236             ENDIF
     1237
     1238!
     1239!--          Write surface positions to file
     1240             nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_xs_surf(av),      &
     1241                                     surfaces%xs, start = (/ surfaces%s(1) /), &
     1242                                     count = (/ surfaces%ns /) )
     1243             CALL netcdf_handle_error( 'surface_data_output_init', 5572 )
     1244
     1245             nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_ys_surf(av),      &
     1246                                     surfaces%ys, start = (/ surfaces%s(1) /), &
     1247                                     count = (/ surfaces%ns /) )
     1248             CALL netcdf_handle_error( 'surface_data_output_init', 5573 )
     1249
     1250             nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_zs_surf(av),      &
     1251                                     surfaces%zs, start = (/ surfaces%s(1) /), &
     1252                                     count = (/ surfaces%ns /) )
     1253             CALL netcdf_handle_error( 'surface_data_output_init', 5574 )
     1254
     1255             nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_etum_surf(av),        &
     1256                                     surfaces%es_utm, start = (/ surfaces%s(1) /), &
     1257                                     count = (/ surfaces%ns /) )
     1258             CALL netcdf_handle_error( 'surface_data_output_init', 5575 )
     1259
     1260             nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_nutm_surf(av),        &
     1261                                     surfaces%ns_utm, start = (/ surfaces%s(1) /), &
     1262                                     count = (/ surfaces%ns /) )
     1263             CALL netcdf_handle_error( 'surface_data_output_init', 5576 )
     1264             
     1265          ENDDO
     1266#endif
     1267
     1268       ENDIF
    8621269
    8631270   END SUBROUTINE surface_data_output_init
     
    8921299      IF ( dosurf_no(av) == 0 )  RETURN
    8931300!
    894 !--   Open files
    895       CALL check_open( 25+av )
    896 !
    897 !--   Write coordinates
    898       IF ( .NOT. first_output(av) )  THEN
    899          DO  i = 0, io_blocks-1
    900             IF ( i == io_group )  THEN
    901                WRITE ( 25+av )  surfaces%npoints
    902                WRITE ( 25+av )  surfaces%npoints_total
    903                WRITE ( 25+av )  surfaces%ns
    904                WRITE ( 25+av )  surfaces%ns_total
    905                WRITE ( 25+av )  surfaces%points
    906                WRITE ( 25+av )  surfaces%polygons
    907             ENDIF
     1301!--   In case of VTK output, check if binary files are open and write coordinates.
     1302      IF ( to_vtk )  THEN
     1303
     1304         CALL check_open( 25+av )
     1305
     1306         IF ( .NOT. first_output(av) )  THEN
     1307            DO  i = 0, io_blocks-1
     1308               IF ( i == io_group )  THEN
     1309                  WRITE ( 25+av )  surfaces%npoints
     1310                  WRITE ( 25+av )  surfaces%npoints_total
     1311                  WRITE ( 25+av )  surfaces%ns
     1312                  WRITE ( 25+av )  surfaces%ns_total
     1313                  WRITE ( 25+av )  surfaces%points
     1314                  WRITE ( 25+av )  surfaces%polygons
     1315               ENDIF
    9081316#if defined( __parallel )
    909             CALL MPI_BARRIER( comm2d, ierr )
     1317               CALL MPI_BARRIER( comm2d, ierr )
    9101318#endif
    911             first_output(av) = .TRUE.
    912          ENDDO
     1319               first_output(av) = .TRUE.
     1320            ENDDO
     1321         ENDIF
    9131322      ENDIF
    914 
     1323!
     1324!--   In case of NetCDF output, check if enough time steps are available in file
     1325!--   and update time axis.
     1326      IF ( to_netcdf )  THEN
     1327         IF ( dosurf_time_count(av) + 1 > ntdim_surf(av) )  THEN
     1328            WRITE ( message_string, * ) 'Output of surface data is not given at t=',          &
     1329                                        time_since_reference_point, 's because the maximum ', &
     1330                                        'number of output time levels is exceeded.'
     1331            CALL message( 'surface_data_output', 'PA0???', 0, 1, 0, 6, 0 )
     1332           
     1333            RETURN
     1334
     1335         ENDIF
     1336!
     1337!--      Update the netCDF time axis
     1338!--      In case of parallel output, this is only done by PE0 to increase the
     1339!--      performance.
     1340         dosurf_time_count(av) = dosurf_time_count(av) + 1
     1341         IF ( myid == 0 )  THEN
     1342            nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_time_surf(av),  &
     1343                                    (/ time_since_reference_point /),       &
     1344                                    start = (/ dosurf_time_count(av) /),    &
     1345                                    count = (/ 1 /) )
     1346            CALL netcdf_handle_error( 'surface_data_output', 6666 )
     1347         ENDIF
     1348      ENDIF
     1349
     1350!
     1351!--   Cycle through output quantities and write them to file.
    9151352      n_out = 0
    9161353      DO  WHILE ( dosurf(av,n_out+1)(1:1) /= ' ' )
     
    22692706!--      - Dimension: 1-nsurfaces, 1-npoints - can be ordered consecutively
    22702707!--      - Distinguish between average and non-average data
    2271          DO  i = 0, io_blocks-1
    2272             IF ( i == io_group )  THEN
    2273                WRITE ( 25+av )  LEN_TRIM( 'time' )
    2274                WRITE ( 25+av )  'time'
    2275                WRITE ( 25+av )  time_since_reference_point
    2276                WRITE ( 25+av )  LEN_TRIM( trimvar )
    2277                WRITE ( 25+av )  TRIM( trimvar )
    2278                WRITE ( 25+av )  surfaces%var_out
    2279             ENDIF
     2708         IF ( to_vtk )  THEN
     2709            DO  i = 0, io_blocks-1
     2710               IF ( i == io_group )  THEN
     2711                  WRITE ( 25+av )  LEN_TRIM( 'time' )
     2712                  WRITE ( 25+av )  'time'
     2713                  WRITE ( 25+av )  time_since_reference_point
     2714                  WRITE ( 25+av )  LEN_TRIM( trimvar )
     2715                  WRITE ( 25+av )  TRIM( trimvar )
     2716                  WRITE ( 25+av )  surfaces%var_out
     2717               ENDIF
    22802718#if defined( __parallel )
    2281             CALL MPI_BARRIER( comm2d, ierr )
     2719               CALL MPI_BARRIER( comm2d, ierr )
    22822720#endif
    2283          ENDDO
     2721            ENDDO
     2722         ENDIF
     2723         
     2724         IF ( to_netcdf )  THEN
     2725!
     2726!--         Write output array to file
     2727            nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_dosurf(av,n_out),            &
     2728                                    surfaces%var_out,                                    &
     2729                                    start = (/ surfaces%s(1), dosurf_time_count(av) /),  &
     2730                                    count = (/ surfaces%ns, 1 /) )
     2731            CALL netcdf_handle_error( 'surface_data_output', 6667 )
     2732         
     2733         ENDIF
    22842734
    22852735      ENDDO
     
    35534003                                  averaging_interval_surf, data_output_surf,   &
    35544004                                  dt_dosurf, dt_dosurf_av,                     &
    3555                                   skip_time_dosurf, skip_time_dosurf_av
     4005                                  skip_time_dosurf, skip_time_dosurf_av,       &
     4006                                  to_netcdf, to_vtk
    35564007                                 
    35574008       line = ' '
     
    35994050
    36004051       CHARACTER(LEN=100) ::  trimvar !< dummy for single output variable
     4052       CHARACTER(LEN=100) ::  unit    !< dummy for unit of output variable
    36014053
    36024054       INTEGER(iwp) ::  av    !< id indicating average or non-average data output
     
    36214073                        'PA0087', 1, 2, 0, 6, 0 )
    36224074       ENDIF
    3623 
     4075       
     4076#if ! defined( __netcdf4_parallel )
     4077!
     4078!--    Surface output via NetCDF requires parallel NetCDF
     4079       IF ( to_netcdf )  THEN
     4080          message_string = 'to_netcdf = .True. requires parallel NetCDF'
     4081          CALL message( 'surface_data_output_check_parameters',                &
     4082                        'PA0116', 1, 2, 0, 6, 0 )
     4083       ENDIF
     4084#endif
    36244085!
    36254086!--   Count number of output variables and separate output strings for
     
    36464107
    36474108!
    3648 !--      Check if all output variables are known
     4109!--      Check if all output variables are known and assign a unit
     4110         unit = 'not set'
    36494111         SELECT CASE ( TRIM( trimvar ) )
    36504112
     
    36554117                             'PA0087', 1, 2, 0, 6, 0 )
    36564118
    3657             CASE ( 'us', 'ts', 'qs', 'ss', 'qcs', 'ncs', 'qrs', 'nrs', 'ol' )
    3658 
    3659             CASE ( 'z0', 'z0h', 'z0q' )
    3660 
    3661             CASE ( 'theta1', 'qv1', 'thetav1' )
     4119            CASE ( 'us', 'uvw1' )
     4120               unit = 'm/s'
     4121
     4122            CASE ( 'ss', 'qcs', 'ncs', 'qrs', 'nrs' )
     4123               unit = '1'
     4124
     4125            CASE ( 'z0', 'z0h', 'z0q', 'ol' )
     4126               unit = 'm'
     4127
     4128            CASE ( 'ts', 'theta1', 'thetav1', 'theta_surface', 'thetav_surface' )
     4129               unit = 'K'
    36624130
    36634131            CASE ( 'usws', 'vsws' )
     4132               unit = 'm2/s2'
    36644133           
    3665             CASE ( 'uvw1' )
    3666 
    36674134            CASE ( 'qcsws', 'ncsws', 'qrsws', 'nrsws', 'sasws' )
    36684135
    3669             CASE ( 'shf', 'qsws', 'ssws' )
    3670 
    3671             CASE ( 'q_surface', 'theta_surface', 'thetav_surface' )
     4136            CASE ( 'shf' )
     4137               unit = 'K m/s'
     4138
     4139            CASE ( 'qsws' )
     4140               unit = 'kg/kg m/s'
     4141
     4142            CASE ( 'ssws' )
     4143               unit = 'kg/m2/s'
     4144
     4145            CASE ( 'qs', 'q_surface', 'qv1' )
     4146               unit = 'kg/kg'
    36724147
    36734148            CASE ( 'rad_net' )
     4149               unit = 'W/m2'
    36744150
    36754151            CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_lw_dif', 'rad_lw_ref',     &
    36764152                   'rad_lw_res' )
     4153               unit = 'W/m2'
    36774154
    36784155            CASE ( 'rad_sw_in', 'rad_sw_out', 'rad_sw_dif', 'rad_sw_ref',     &
    36794156                   'rad_sw_res', 'rad_sw_dir' )
     4157               unit = 'W/m2'
    36804158
    36814159            CASE ( 'ghf' )
     4160               unit = 'W/m2'
    36824161
    36834162            CASE ( 'r_a', 'r_canopy', 'r_soil', 'r_s' )
     4163               unit = 's/m'
    36844164
    36854165            CASE DEFAULT
     
    36894169                             'PA0087', 1, 2, 0, 6, 0 )
    36904170         END SELECT
     4171
     4172         dosurf_unit(av,dosurf_no(av)) = unit
    36914173
    36924174       ENDDO
     
    37174199      IF ( dosurf_no(av) == 0 )  RETURN
    37184200!
    3719 !--   Open files
    3720       CALL check_open( 25+av )
    3721 !
    3722 !--   Write time coordinate
    3723       DO  i = 0, io_blocks-1
    3724          IF ( i == io_group )  THEN
    3725             WRITE ( 25+av )  LEN_TRIM( 'END' )
    3726             WRITE ( 25+av )  'END'
    3727          ENDIF
     4201!--   If output to VTK files is enabled, check if files are open and write
     4202!--   an end-of-file statement.
     4203      IF ( to_vtk )  THEN
     4204         CALL check_open( 25+av )
     4205!
     4206!--      Write time coordinate
     4207         DO  i = 0, io_blocks-1
     4208            IF ( i == io_group )  THEN
     4209               WRITE ( 25+av )  LEN_TRIM( 'END' )
     4210               WRITE ( 25+av )  'END'
     4211            ENDIF
    37284212#if defined( __parallel )
    3729          CALL MPI_BARRIER( comm2d, ierr )
     4213            CALL MPI_BARRIER( comm2d, ierr )
    37304214#endif
    3731       ENDDO   
     4215         ENDDO   
     4216      ENDIF
    37324217
    37334218    END SUBROUTINE surface_data_output_last_action   
Note: See TracChangeset for help on using the changeset viewer.