source: palm/trunk/UTIL/surface_output_processing/surface_output_to_vtk.f90 @ 3841

Last change on this file since 3841 was 3841, checked in by suehring, 2 years ago

In VTK postprocessor add suffix to indicate average data in VTK files

  • Property svn:keywords set to Id
File size: 25.8 KB
Line 
1!> @file surface_output_merge.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2018  Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: surface_output_to_vtk.f90 3841 2019-03-29 11:38:44Z suehring $
27! Add suffix to VTK output to indicate average data
28!
29! 3805 2019-03-20 15:26:35Z raasch
30! output format adjusted
31!
32! 3755 2019-02-19 17:03:19Z suehring
33! - Change format description to avoid *** for larger domains
34! - switch back from to non-standard Fortran ftell and fseek since posix
35!   interface do not work any more for some unknown reason
36!
37! 3523 2018-11-13 16:09:31Z suehring
38! Implement interface for posix conform C-systemcalls in order to replace
39! non-standard FORTRAN functions ftell and fseek.
40!
41! 3496 2018-11-06 15:59:50Z suehring
42! Use subroutine call for fseek instead of function call. gfortran has some
43! problems with this.
44!
45! 3494 2018-11-06 14:51:27Z suehring
46! Initial version
47!
48! Authors:
49! --------
50! @author Matthias Suehring and Klaus Ketelsen
51!
52!------------------------------------------------------------------------------!
53! Description:
54! ------------
55!> Interface to posix systemcalls. The module allows to call the following
56!> C system calls from FORTRAN: ftell, lseek.
57!------------------------------------------------------------------------------!
58 MODULE posix_interface
59 
60    INTERFACE
61       INTEGER (C_SIZE_T) FUNCTION C_LSEEK (fd, offset, whence)                &
62                                   BIND(C, NAME='lseek')
63         
64          USE ISO_C_BINDING
65         
66          IMPLICIT NONE
67         
68          INTEGER(KIND=C_INT),    VALUE ::  fd
69          INTEGER(KIND=C_SIZE_T), VALUE ::  offset
70          INTEGER(KIND=C_INT),    VALUE ::  whence
71         
72       END FUNCTION C_LSEEK
73       
74    END INTERFACE
75   
76    INTERFACE
77       INTEGER (C_SIZE_T) FUNCTION C_FTELL ( fd )                              &
78                                   BIND(C, NAME='ftell')
79         
80          USE ISO_C_BINDING
81         
82          IMPLICIT NONE
83         
84          INTEGER(KIND=C_INT),    VALUE ::  fd
85
86       END FUNCTION C_FTELL
87       
88    END INTERFACE
89   
90    PUBLIC posix_lseek, posix_ftell
91   
92    CONTAINS
93!
94!------------------------------------------------------------------------------!
95! Description:
96! ------------
97!> Interface for the C-routine lseek.
98!------------------------------------------------------------------------------!
99    SUBROUTINE posix_lseek( fid, offset )
100   
101       USE ISO_C_BINDING
102   
103       IMPLICIT NONE
104
105       INTEGER,INTENT(IN)                     ::  fid    !< file unit
106       INTEGER(KIND=C_SIZE_T),INTENT(IN)      ::  offset !< file offset from the beginning
107
108       INTEGER(KIND=C_INT)                    ::  my_fid !< file unit
109       INTEGER(KIND=C_SIZE_T)                 ::  retval !< return value
110       INTEGER(KIND=C_INT)                    ::  whence !< mode, here start from the beginning
111
112       my_fid = fid
113       whence = 0
114
115       retval = C_LSEEK( fid, offset, whence )
116       
117    END SUBROUTINE posix_lseek
118!
119!------------------------------------------------------------------------------!
120! Description:
121! ------------
122!> Interface for the C-routine ftell.
123!------------------------------------------------------------------------------!
124    SUBROUTINE posix_ftell( fid, filepos )
125   
126       USE ISO_C_BINDING
127   
128       IMPLICIT NONE
129
130       INTEGER,INTENT(IN)      ::  fid                     !< file unit
131       INTEGER(KIND=C_SIZE_T), INTENT(INOUT)   ::  filepos !< file position
132
133       filepos = C_FTELL( fid )
134       
135    END SUBROUTINE posix_ftell
136
137 END MODULE posix_interface
138!
139!------------------------------------------------------------------------------!
140! Description:
141! ------------
142!> This routine combines surface output from PALM-subdomains into one file.
143!> Output from every processor element is opened and read and output
144!> from all processor elements are written into one file for each timestep
145!> according to Paraview's VTK standard.
146!> Output is distinguished between instantaneous and time-averaged data.
147!------------------------------------------------------------------------------!
148 PROGRAM surface_output_to_vtk_netcdf
149 
150    USE, INTRINSIC ::  ISO_C_BINDING
151   
152#if defined( __netcdf ) 
153    USE NETCDF
154#endif
155 
156    USE posix_interface,                                                       &
157        ONLY:  posix_ftell, posix_lseek
158
159    IMPLICIT NONE
160   
161    CHARACTER(LEN=4)   ::  char_time              !< string indicating simulated time
162    CHARACTER(LEN=4)   ::  file_suffix = '.bin'   !< string which contain the suffix indicating surface data
163   
164    CHARACTER(LEN=10)  ::  char_dum               !< dummy string
165   
166    CHARACTER(LEN=30)  ::  myid_char              !< combined string indicating binary file
167   
168
169    CHARACTER(LEN=100) ::  path                   !< path to the binary data
170    CHARACTER(LEN=100) ::  run                    !< name of the run
171    CHARACTER(LEN=100) ::  variable_name          !< name of the processed output variable   
172
173    INTEGER(4)   ::  ftell                      !< intrinsic function, get current position in file
174    INTEGER(4)   ::  ndum                       !< return parameter of intrinsic function fseek   
175       
176    INTEGER, PARAMETER ::  iwp = 4                !< integer precision
177    INTEGER, PARAMETER ::  wp  = 8                !< float precision
178    INTEGER, PARAMETER ::  OFFSET_KIND = C_SIZE_T !< unsigned integer for the C-interface
179
180    INTEGER(iwp) ::  cycle_number = 0             !< cycle number
181    INTEGER(iwp) ::  f                            !< running index over all binary files
182    INTEGER(iwp) ::  file_id_in = 18              !< file unit for input binaray file   
183    INTEGER(iwp) ::  file_id_out = 20             !< file unit for output VTK file       
184    INTEGER(iwp) ::  file_id_out_header = 19      !< file unit for temporary header file
185    INTEGER(iwp) ::  length                       !< length of string on file
186    INTEGER(iwp) ::  n                            !< running index over all points and polygons
187    INTEGER(iwp) ::  npoints_total                !< total number of points
188    INTEGER(iwp) ::  ns_total                     !< total number of polygons
189    INTEGER(iwp) ::  num_pe                       !< number of processors of the run
190
191!     INTEGER(OFFSET_KIND),DIMENSION(:), ALLOCATABLE ::  filepos !< current fileposition in binary file
192    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  filepos !< current fileposition in binary file
193   
194    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  npoints !< number of points/vertices in a binaray file
195    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ns      !< number of surface elements in a binaray file
196   
197    LOGICAL ::  convert_average_data = .FALSE.          !< namelist parameter to decide whether average or instantaneous data should be converted
198    LOGICAL ::  convert_to_netcdf    = .FALSE.          !< namelist parameter indicating that binary data is also converted to NetCDF files
199    LOGICAL, DIMENSION(:), ALLOCATABLE      ::  eof     !< flag indicating that end of binary file is reached
200   
201    REAL(wp)                              ::  simulated_time !< output time   
202   
203    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  var            !< actual surface data
204   
205    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  points         !< point / vertex data
206    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  polygons       !< polygon data
207   
208    logical :: flag
209
210!
211!-- Read namelist.
212    CALL surface_output_parin
213!
214!-- Allocate array which contains the file position in each output file,
215!-- in order to skip some reading.
216    ALLOCATE( eof(0:num_pe-1)     )
217    ALLOCATE( filepos(0:num_pe-1) )
218    ALLOCATE( npoints(0:num_pe-1) )
219    ALLOCATE( ns(0:num_pe-1)      )
220!
221!-- Initialize file position.
222    filepos = 0
223!
224!-- Open a temporary file which contains all necessary header information for the
225!-- VTK format.
226    OPEN ( file_id_out_header, FILE = 'HEADER', STATUS = 'REPLACE',            &
227           FORM = 'FORMATTED' )
228!
229!-- READ grid setup, i.e. the number and position of vertices and surface elements
230!-- and merge all this information into one file. Further, create all the required
231!-- header information of the VTK file.
232!-- Note, PARAVIEW expects one VTK file for each variable and each timestep.
233!-- Hence, header information needs to be duplicated multiple times.
234!-- Moreover, Paraview expects consecutive vertex and polygon data, which are
235!-- all distributed over the binaray files. Hence, first read vertex data from
236!-- binary file, write this to the HEADER file, close the binary file and read
237!-- data from the next binary file and so on. This requires several openenings
238!-- and closings of the binary files and temporarily storage of the
239!-- file-positions.
240    DO  f = 0, num_pe - 1
241!
242!--    Create filename of the treated binary file.
243       CALL surface_output_create_file_string
244!
245!--    Open file with surface output for processor f.
246       OPEN ( file_id_in, FILE = TRIM( path ) // TRIM( run ) //                &
247              TRIM( myid_char ), FORM = 'UNFORMATTED' )
248!
249!--    Read number of vertices / points and surface elements
250       READ ( file_id_in )  npoints(f)
251       READ ( file_id_in )  npoints_total
252       READ ( file_id_in )  ns(f)
253       READ ( file_id_in )  ns_total
254!
255!--    If binary surface dater shall be converted into NetCDF, create the file.
256!--    Only for first call.
257       IF ( convert_to_netcdf  .AND.  f == 0 )  THEN
258          CALL surface_output_create_netcdf_file
259       ENDIF
260!
261!--    Allocate arrays where all the surfaces and vertices will be stored.
262       ALLOCATE( points(3,1:npoints(f))   )
263!
264!--    Read polygon data and store them in a temporary file.
265       READ ( file_id_in )  points
266!
267!--    Obtain current file position. Will be stored for next file opening.
268!        CALL posix_ftell( file_id_in, filepos(f) )
269       filepos(f) = ftell( file_id_in )
270!
271!--    Write header information. Only one time required.
272       IF ( f == 0 )  THEN
273          WRITE( file_id_out_header,'(A)' ) "# vtk DataFile Version 3.0"
274          WRITE( file_id_out_header,'(A,F8.2,A)' ) "legacy vtk File generated by PALM,  simulation time = xxx sec"
275          WRITE( file_id_out_header,'(A)') "ASCII"
276
277          WRITE( file_id_out_header,'(A)') "DATASET POLYDATA"
278          WRITE( file_id_out_header,'(A,I12,A)') "POINTS ", npoints_total, " float"
279       ENDIF     
280!
281!--    Write the vertex data into header file.
282       DO  n = 1, npoints(f)
283          WRITE( file_id_out_header, '(8F15.1)' )  points(1:3,n)
284       ENDDO
285!
286!--    Deallocate vertex data and close binary file.
287       DEALLOCATE( points )
288       
289       CLOSE ( file_id_in )
290    ENDDO
291!
292!-- Now, treat polygon data.
293    DO  f = 0, num_pe - 1
294!
295!--    Create filename of the treated binary file .   
296       CALL surface_output_create_file_string
297!
298!--    Open file with surface output for processor f.
299       OPEN ( file_id_in, FILE = TRIM( path ) // TRIM( run ) //                &
300              TRIM( myid_char ), FORM = 'UNFORMATTED' )
301!
302!--    Move to last postion.
303!        CALL posix_lseek( file_id_in, filepos(f) )
304       CALL FSEEK( file_id_in, filepos(f), 0, ndum )
305!
306!--    Allocate array for polygon data
307       ALLOCATE( polygons(5,1:ns(f)) )
308!
309!--    Read polygon data and store them in a temporary file.
310       READ ( file_id_in )  polygons
311!
312!--    Obtain current file position after reading the local polygon data.
313!--    Will be used for next file opening.
314!        CALL posix_ftell( file_id_in, filepos(f) )
315       filepos(f) = ftell( file_id_in )
316!
317!--    Write further header information. Only one time required.
318       IF ( f == 0 )                                                           &
319          WRITE ( file_id_out_header, '(A,2I10)') "POLYGONS ",                 &
320                                                  ns_total, 5 * ns_total
321!
322!--    Write the polygons into the header file. Note, polygon data is of type
323!--    integer, as it just connects the point-indices which describe the given
324!--    surface element.
325       DO n = 1, ns(f)
326          WRITE ( file_id_out_header, '(5I18)' )  INT( polygons(1:5,n) )
327       ENDDO
328!
329!--    Deallocate array for polygon data and close the file.
330       DEALLOCATE( polygons )
331       
332       CLOSE ( file_id_in )
333       
334    ENDDO
335   
336    f = 0
337    CALL surface_output_create_file_string
338!
339!-- Write further header information. Only once required.
340    WRITE ( file_id_out_header, '(A,I10)') "CELL_DATA ", ns_total
341    WRITE ( file_id_out_header, '(A,I10)') "SCALARS cell_scalars float 1 "
342    WRITE ( file_id_out_header, '(A,I10)') "LOOKUP_TABLE default "
343!
344!-- Header creation has now been completed. File can be closed.   
345    CLOSE( file_id_out_header )
346!
347!-- Now, read all the actual data from the surface output. Please note, Paraview
348!-- VTK format expects 1 file per variable and time step.
349!-- The output file is only created once and includes the variable and the
350!-- simulated time.
351!-- In the binaray files, several variables and timesteps are stored, data for a
352!-- given variable, however, is distributed over all binary files. Hence, read
353!-- variable data for a given timestep from the binary file, write this data into
354!-- the target output file, remember file position in binary file and close it,
355!-- open nex binary file and read the variable, and so on, until all variables
356!-- for all timesteps are processed.
357    eof = .FALSE.
358    DO
359       DO  f = 0, num_pe - 1
360!
361!--       Clean up strings-
362          char_time = ''
363          variable_name = ''
364!
365!--       Create filename of the treated binary file.           
366          CALL surface_output_create_file_string
367!
368!--       Open binary file with surface output for processor f.
369          OPEN ( file_id_in, FILE = TRIM( path ) // TRIM( run ) //         &
370                 TRIM( myid_char ), FORM = 'UNFORMATTED' )
371!
372!--       Move to last postion.
373!           CALL posix_lseek( file_id_in, filepos(f) )
374          CALL FSEEK( file_id_in, filepos(f), 0, ndum )
375!
376!--       Read string length and string indicating the output time step.
377          READ ( file_id_in ) length
378          READ ( file_id_in ) char_time(1:length)
379!
380!--       If string for the output time indicates the end-of-file, set the eof
381!--       flag and skip the read of the loop.
382          IF ( char_time(1:length) == 'END' )  THEN
383             eof(f) = .TRUE. 
384             CLOSE ( file_id_in )
385             CYCLE
386          ENDIF
387!
388!--       Read output time, and variable name.
389          READ ( file_id_in ) simulated_time
390          READ ( file_id_in ) length
391          READ ( file_id_in ) variable_name(1:length)
392!
393!--       For first loop index, open the target output file. First create the
394!--       filename string. Further, copy HEADER file with the given filename
395!--       string. The header information must be given in each VTK file!
396          IF ( f == 0 )  THEN
397             IF ( simulated_time < 10.0_wp )  THEN
398                WRITE( char_dum, '(I1)' )  INT( simulated_time ) 
399             ELSEIF ( simulated_time < 100.0_wp )  THEN
400                WRITE( char_dum, '(I2)' )  INT( simulated_time )
401             ELSEIF ( simulated_time < 1000.0_wp )  THEN
402                WRITE( char_dum, '(I3)' )  INT( simulated_time )
403             ELSEIF ( simulated_time < 10000.0_wp )  THEN
404                WRITE( char_dum, '(I4)' )  INT( simulated_time )   
405             ELSEIF ( simulated_time < 100000.0_wp )  THEN
406                WRITE( char_dum, '(I5)' )  INT( simulated_time )   
407             ELSEIF ( simulated_time < 1000000.0_wp )  THEN
408                WRITE( char_dum, '(I6)' )  INT( simulated_time )
409             ELSEIF ( simulated_time < 10000000.0_wp )  THEN
410                WRITE( char_dum, '(I7)' )  INT( simulated_time )   
411             ELSEIF ( simulated_time < 100000000.0_wp )  THEN
412                WRITE( char_dum, '(I8)' )  INT( simulated_time )   
413             ELSEIF ( simulated_time < 1000000000.0_wp )  THEN
414                WRITE( char_dum, '(I9)' )  INT( simulated_time )
415             ENDIF             
416!
417!--          Copy HEADER file and open VTK file
418             IF ( convert_average_data )  THEN
419                CALL system('cp HEADER ' //                                    &
420                            TRIM( path ) // TRIM( char_dum ) //                & 
421                            '_AV_' // 's_' // TRIM( variable_name ) // '.vtk' )
422                           
423                OPEN ( file_id_out, FILE = TRIM( path ) // TRIM( char_dum ) // & 
424                       '_AV_' // 's_' // TRIM( variable_name ) // '.vtk',      &
425                       FORM='FORMATTED', POSITION = 'APPEND' ) 
426             ELSE
427                CALL system('cp HEADER ' //                                    &
428                            TRIM( path ) // TRIM( char_dum ) //                & 
429                            's_' // TRIM( variable_name ) // '.vtk' )
430                           
431                OPEN ( file_id_out, FILE = TRIM( path ) // TRIM( char_dum ) // & 
432                       's_' // TRIM( variable_name ) // '.vtk',                &
433                       FORM='FORMATTED', POSITION = 'APPEND' ) 
434             ENDIF
435
436          ENDIF
437!
438!--       Allocate and read array for variable data. 
439          ALLOCATE( var(1:ns(f)) )
440       
441          READ( file_id_in ) var
442!
443!--       Write variable data into output VTK file.
444          DO n = 1, ns(f)
445             WRITE( file_id_out, * ) var(n)
446          ENDDO
447!
448!--       Remember file position in binary file and close it.
449!           CALL posix_ftell( file_id_in, filepos(f) )
450          filepos(f) = ftell( file_id_in )
451       
452          CLOSE ( file_id_in )
453!
454!--       Deallocate temporary array for variable data.
455          DEALLOCATE( var )
456       
457       ENDDO
458!
459!--    After data for a variable for one time step is read, close the output
460!--    VTK file and go to next variable or timestep.
461       CLOSE ( file_id_out )
462!
463!--    If all files reached the end-of-file, exit the loop.
464       IF ( ALL( eof ) )  EXIT
465       
466    ENDDO
467!
468!-- Finally, remove HEADER file
469    CALL system( 'rm HEADER' )
470   
471 CONTAINS
472 
473!------------------------------------------------------------------------------!
474! Description:
475! ------------
476!> This subroutine read the namelist file.
477!------------------------------------------------------------------------------!
478    SUBROUTINE surface_output_parin
479       
480       IMPLICIT NONE
481       
482       INTEGER(iwp) ::  file_id_parin = 90
483       
484       NAMELIST /surface_output/  convert_average_data, convert_to_netcdf,     &
485                                  cycle_number, num_pe, path, run
486
487!
488!--    Open namelist file.
489       OPEN( file_id_parin, FILE='surface_output_parin',                       &
490             STATUS='OLD', FORM='FORMATTED')
491!
492!--    Read namelist.
493       READ ( file_id_parin, surface_output )
494!
495!--    Close namelist file.
496       CLOSE( file_id_parin )
497       
498    END SUBROUTINE surface_output_parin
499     
500!------------------------------------------------------------------------------!
501! Description:
502! ------------
503!> This subroutine creates the filename string of the treated binary file.
504!------------------------------------------------------------------------------!
505    SUBROUTINE surface_output_create_file_string
506       
507       IMPLICIT NONE
508       
509       CHARACTER(LEN=3) ::  char_av = ''
510       CHARACTER(LEN=4) ::  char_cycle = ''
511       
512!
513!--    Create substring for the cycle number.
514       IF ( cycle_number /= 0 )  THEN
515          IF ( cycle_number < 10 )  THEN
516             WRITE( char_cycle, '(I1)')  cycle_number
517             char_cycle = '.00' // TRIM( char_cycle )
518          ELSEIF ( cycle_number < 100 )  THEN
519             WRITE( char_cycle, '(I2)')  cycle_number
520             char_cycle = '.0' // TRIM( char_cycle )
521          ELSEIF ( cycle_number < 1000 )  THEN
522             WRITE( char_cycle, '(I3)')  cycle_number
523             char_cycle = '.' // TRIM( char_cycle )
524          ENDIF
525       ENDIF
526!
527!--    Create substring for averaged data.
528       IF ( convert_average_data )  char_av = '_av'
529!
530!--    Create substring for the processor id and combine all substrings.
531       IF ( f < 10 )  THEN
532          WRITE( char_dum, '(I1)')  f
533          myid_char = TRIM( char_av ) // '_surf_00000' // TRIM( char_dum ) //  &
534                      TRIM( char_cycle ) // file_suffix
535       ELSEIF ( f < 100     )  THEN
536          WRITE( char_dum, '(I2)')  f
537          myid_char = TRIM( char_av ) // '_surf_0000'  // TRIM( char_dum ) //  &
538                      TRIM( char_cycle ) // file_suffix
539       ELSEIF ( f < 1000    )  THEN
540          WRITE( char_dum, '(I3)')  f
541          myid_char = TRIM( char_av ) // '_surf_000'   // TRIM( char_dum ) //  &
542                      TRIM( char_cycle ) // file_suffix
543       ELSEIF ( f < 10000   )  THEN
544          WRITE( char_dum, '(I4)')  f
545          myid_char = TRIM( char_av ) // '_surf_00'    // TRIM( char_dum ) //  &
546                      TRIM( char_cycle ) // file_suffix
547       ELSEIF ( f < 100000  )  THEN
548          WRITE( char_dum, '(I5)')  f
549          myid_char = TRIM( char_av ) // '_surf_0'     // TRIM( char_dum ) //  &
550                      TRIM( char_cycle ) // file_suffix
551       ELSEIF ( f < 1000000 )  THEN
552          WRITE( char_dum, '(I6)')  f
553          myid_char = TRIM( char_av ) // '_surf_'      // TRIM( char_dum ) //  &
554                      TRIM( char_cycle ) // file_suffix
555       ENDIF
556       
557    END SUBROUTINE surface_output_create_file_string
558   
559   
560!------------------------------------------------------------------------------!
561! Description:
562! ------------
563!> This subroutine creates the NetCDF file and defines the dimensions
564!------------------------------------------------------------------------------!
565    SUBROUTINE surface_output_create_netcdf_file
566       
567       IMPLICIT NONE
568       
569       
570       CHARACTER(LEN=5)   ::  char_cycle = ''
571       CHARACTER(LEN=200) ::  nc_filename
572       
573!
574!--    Create substring for the cycle number.
575       IF ( cycle_number /= 0 )  THEN
576          IF ( cycle_number < 10 )  THEN
577             WRITE( char_cycle, '(I1)')  cycle_number
578             char_cycle = '.00' // TRIM( char_cycle ) // '.'
579          ELSEIF ( cycle_number < 100 )  THEN
580             WRITE( char_cycle, '(I2)')  cycle_number
581             char_cycle = '.0' // TRIM( char_cycle ) // '.'
582          ELSEIF ( cycle_number < 1000 )  THEN
583             WRITE( char_cycle, '(I3)')  cycle_number
584             char_cycle = '.' // TRIM( char_cycle ) // '.'
585          ENDIF
586       ELSE
587          char_cycle = '.'
588       ENDIF
589#if defined( __netcdf ) 
590!
591!-     Create file
592       nc_filename = TRIM( path ) // TRIM(run) // '_surf' //                   &
593                     TRIM( char_cycle ) // 'nc'
594                     
595       status_nc = NF90_CREATE( nc_filename, IOR( NF90_CLOBBER, NF90_NETCDF4 ),&
596                                nc_id )
597       CALL handle_error( 'create netcdf file' )
598!
599!--    Define dimensions
600       status_nc = NF90_DEF_DIM( nc_id, 'xs', ns_total, id_xs )
601       CALL handle_error( 'define dimension xs' )
602       status_nc = NF90_DEF_DIM( nc_id, 'ys', ns_total, id_ys )
603       CALL handle_error( 'define dimension ys' )
604       status_nc = NF90_DEF_DIM( nc_id, 'zs', ns_total, id_zs )
605       CALL handle_error( 'define dimension zs' )
606       status_nc = NF90_DEF_DIM( nc_id, 'azimuth', ns_total, id_az )
607       CALL handle_error( 'define dimension azimuth' )
608       status_nc = NF90_DEF_DIM( nc_id, 'zenith',  ns_total, id_ze )
609       CALL handle_error( 'define dimension zenith' )
610       status_nc = NF90_DEF_DIM( nc_id, 'time',  NF90_UNLIMITED, id_time )
611       CALL handle_error( 'define dimension time' )
612!
613!--    Define dimension variables
614       status_nc = NF90_DEF_VAR( nc_id, 'xs', NF90_FLOAT,                      &
615                                 (/ id_xs /), id_var_xs )
616       CALL handle_error( 'define variable xs' )
617       status_nc = NF90_DEF_VAR( nc_id, 'ys', NF90_FLOAT,                      &
618                                 (/ id_ys /), id_var_ys )
619       CALL handle_error( 'define variable ys' )
620       status_nc = NF90_DEF_VAR( nc_id, 'zs', NF90_FLOAT,                      &
621                                 (/ id_zs /), id_var_zs )
622       CALL handle_error( 'define variable zs' ) 
623       status_nc = NF90_DEF_VAR( nc_id, 'azimuth', NF90_FLOAT,                 &
624                                 (/ id_az /), id_var_az ) 
625       CALL handle_error( 'define variable azimuth' )         
626       status_nc = NF90_DEF_VAR( nc_id, 'zenith', NF90_FLOAT,                  &
627                                 (/ id_ze /), id_var_ze ) 
628       CALL handle_error( 'define variable zenith' )                         
629       status_nc = NF90_DEF_VAR( nc_id, 'time', NF90_FLOAT,                    &
630                                 (/ id_time /), id_var_time )
631       CALL handle_error( 'define variable time' )
632#endif       
633    END SUBROUTINE surface_output_create_netcdf_file
634   
635   
636!------------------------------------------------------------------------------!
637! Description:
638! ------------
639!>  Output of NetCDF error code
640!------------------------------------------------------------------------------!   
641    SUBROUTINE handle_error( action )
642
643       IMPLICIT NONE
644
645       CHARACTER(LEN=*) ::  action
646       
647#if defined( __netcdf )       
648       IF ( status_nc /= NF90_NOERR )  THEN
649          PRINT*, TRIM( NF90_STRERROR( status_nc ) ) // ' -- ' // action
650          STOP
651       ENDIF
652#endif
653
654    END SUBROUTINE handle_error
655   
656 END PROGRAM surface_output_to_vtk_netcdf
657
658
659
Note: See TracBrowser for help on using the repository browser.