[3494] | 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 $ |
---|
[3841] | 27 | ! Add suffix to VTK output to indicate average data |
---|
| 28 | ! |
---|
| 29 | ! 3805 2019-03-20 15:26:35Z raasch |
---|
[3805] | 30 | ! output format adjusted |
---|
| 31 | ! |
---|
| 32 | ! 3755 2019-02-19 17:03:19Z suehring |
---|
[3755] | 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 |
---|
[3523] | 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 |
---|
[3496] | 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 |
---|
[3494] | 46 | ! Initial version |
---|
| 47 | ! |
---|
| 48 | ! Authors: |
---|
| 49 | ! -------- |
---|
| 50 | ! @author Matthias Suehring and Klaus Ketelsen |
---|
| 51 | ! |
---|
| 52 | !------------------------------------------------------------------------------! |
---|
| 53 | ! Description: |
---|
| 54 | ! ------------ |
---|
[3523] | 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 | ! ------------ |
---|
[3494] | 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 | !------------------------------------------------------------------------------! |
---|
[3841] | 148 | PROGRAM surface_output_to_vtk_netcdf |
---|
[3523] | 149 | |
---|
| 150 | USE, INTRINSIC :: ISO_C_BINDING |
---|
[3841] | 151 | |
---|
| 152 | #if defined( __netcdf ) |
---|
| 153 | USE NETCDF |
---|
| 154 | #endif |
---|
[3523] | 155 | |
---|
| 156 | USE posix_interface, & |
---|
| 157 | ONLY: posix_ftell, posix_lseek |
---|
[3494] | 158 | |
---|
| 159 | IMPLICIT NONE |
---|
| 160 | |
---|
[3523] | 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 |
---|
[3494] | 163 | |
---|
[3523] | 164 | CHARACTER(LEN=10) :: char_dum !< dummy string |
---|
[3494] | 165 | |
---|
[3523] | 166 | CHARACTER(LEN=30) :: myid_char !< combined string indicating binary file |
---|
[3494] | 167 | |
---|
[3755] | 168 | |
---|
[3523] | 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 |
---|
[3755] | 172 | |
---|
| 173 | INTEGER(4) :: ftell !< intrinsic function, get current position in file |
---|
| 174 | INTEGER(4) :: ndum !< return parameter of intrinsic function fseek |
---|
[3523] | 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 |
---|
[3755] | 179 | |
---|
[3841] | 180 | INTEGER(iwp) :: cycle_number = 0 !< cycle number |
---|
[3523] | 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 |
---|
[3755] | 190 | |
---|
| 191 | ! INTEGER(OFFSET_KIND),DIMENSION(:), ALLOCATABLE :: filepos !< current fileposition in binary file |
---|
| 192 | INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: filepos !< current fileposition in binary file |
---|
[3494] | 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 |
---|
[3841] | 198 | LOGICAL :: convert_to_netcdf = .FALSE. !< namelist parameter indicating that binary data is also converted to NetCDF files |
---|
[3494] | 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 |
---|
[3496] | 207 | |
---|
| 208 | logical :: flag |
---|
[3494] | 209 | |
---|
| 210 | ! |
---|
| 211 | !-- Read namelist. |
---|
[3496] | 212 | CALL surface_output_parin |
---|
[3494] | 213 | ! |
---|
| 214 | !-- Allocate array which contains the file position in each output file, |
---|
| 215 | !-- in order to skip some reading. |
---|
[3755] | 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) ) |
---|
[3494] | 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', & |
---|
[3841] | 227 | FORM = 'FORMATTED' ) |
---|
[3494] | 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. |
---|
[3841] | 233 | !-- Hence, header information needs to be duplicated multiple times. |
---|
[3494] | 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. |
---|
[3496] | 243 | CALL surface_output_create_file_string |
---|
[3494] | 244 | ! |
---|
| 245 | !-- Open file with surface output for processor f. |
---|
[3496] | 246 | OPEN ( file_id_in, FILE = TRIM( path ) // TRIM( run ) // & |
---|
| 247 | TRIM( myid_char ), FORM = 'UNFORMATTED' ) |
---|
[3494] | 248 | ! |
---|
| 249 | !-- Read number of vertices / points and surface elements |
---|
[3496] | 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 |
---|
[3494] | 254 | ! |
---|
[3841] | 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 | ! |
---|
[3494] | 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. |
---|
[3496] | 265 | READ ( file_id_in ) points |
---|
[3494] | 266 | ! |
---|
| 267 | !-- Obtain current file position. Will be stored for next file opening. |
---|
[3755] | 268 | ! CALL posix_ftell( file_id_in, filepos(f) ) |
---|
| 269 | filepos(f) = ftell( file_id_in ) |
---|
[3494] | 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" |
---|
[3805] | 278 | WRITE( file_id_out_header,'(A,I12,A)') "POINTS ", npoints_total, " float" |
---|
[3494] | 279 | ENDIF |
---|
| 280 | ! |
---|
| 281 | !-- Write the vertex data into header file. |
---|
| 282 | DO n = 1, npoints(f) |
---|
[3755] | 283 | WRITE( file_id_out_header, '(8F15.1)' ) points(1:3,n) |
---|
[3494] | 284 | ENDDO |
---|
| 285 | ! |
---|
| 286 | !-- Deallocate vertex data and close binary file. |
---|
| 287 | DEALLOCATE( points ) |
---|
| 288 | |
---|
[3496] | 289 | CLOSE ( file_id_in ) |
---|
[3494] | 290 | ENDDO |
---|
| 291 | ! |
---|
| 292 | !-- Now, treat polygon data. |
---|
| 293 | DO f = 0, num_pe - 1 |
---|
| 294 | ! |
---|
| 295 | !-- Create filename of the treated binary file . |
---|
[3496] | 296 | CALL surface_output_create_file_string |
---|
[3494] | 297 | ! |
---|
| 298 | !-- Open file with surface output for processor f. |
---|
[3496] | 299 | OPEN ( file_id_in, FILE = TRIM( path ) // TRIM( run ) // & |
---|
| 300 | TRIM( myid_char ), FORM = 'UNFORMATTED' ) |
---|
[3494] | 301 | ! |
---|
| 302 | !-- Move to last postion. |
---|
[3755] | 303 | ! CALL posix_lseek( file_id_in, filepos(f) ) |
---|
| 304 | CALL FSEEK( file_id_in, filepos(f), 0, ndum ) |
---|
[3494] | 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. |
---|
[3496] | 310 | READ ( file_id_in ) polygons |
---|
[3494] | 311 | ! |
---|
| 312 | !-- Obtain current file position after reading the local polygon data. |
---|
| 313 | !-- Will be used for next file opening. |
---|
[3755] | 314 | ! CALL posix_ftell( file_id_in, filepos(f) ) |
---|
| 315 | filepos(f) = ftell( file_id_in ) |
---|
[3494] | 316 | ! |
---|
| 317 | !-- Write further header information. Only one time required. |
---|
| 318 | IF ( f == 0 ) & |
---|
[3755] | 319 | WRITE ( file_id_out_header, '(A,2I10)') "POLYGONS ", & |
---|
[3494] | 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) |
---|
[3755] | 326 | WRITE ( file_id_out_header, '(5I18)' ) INT( polygons(1:5,n) ) |
---|
[3494] | 327 | ENDDO |
---|
| 328 | ! |
---|
| 329 | !-- Deallocate array for polygon data and close the file. |
---|
| 330 | DEALLOCATE( polygons ) |
---|
| 331 | |
---|
[3496] | 332 | CLOSE ( file_id_in ) |
---|
[3494] | 333 | |
---|
| 334 | ENDDO |
---|
| 335 | |
---|
| 336 | f = 0 |
---|
[3496] | 337 | CALL surface_output_create_file_string |
---|
[3494] | 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. |
---|
[3496] | 366 | CALL surface_output_create_file_string |
---|
[3494] | 367 | ! |
---|
| 368 | !-- Open binary file with surface output for processor f. |
---|
[3496] | 369 | OPEN ( file_id_in, FILE = TRIM( path ) // TRIM( run ) // & |
---|
| 370 | TRIM( myid_char ), FORM = 'UNFORMATTED' ) |
---|
[3494] | 371 | ! |
---|
| 372 | !-- Move to last postion. |
---|
[3755] | 373 | ! CALL posix_lseek( file_id_in, filepos(f) ) |
---|
| 374 | CALL FSEEK( file_id_in, filepos(f), 0, ndum ) |
---|
[3494] | 375 | ! |
---|
| 376 | !-- Read string length and string indicating the output time step. |
---|
[3496] | 377 | READ ( file_id_in ) length |
---|
| 378 | READ ( file_id_in ) char_time(1:length) |
---|
[3494] | 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. |
---|
[3496] | 384 | CLOSE ( file_id_in ) |
---|
[3494] | 385 | CYCLE |
---|
| 386 | ENDIF |
---|
| 387 | ! |
---|
| 388 | !-- Read output time, and variable name. |
---|
[3496] | 389 | READ ( file_id_in ) simulated_time |
---|
| 390 | READ ( file_id_in ) length |
---|
| 391 | READ ( file_id_in ) variable_name(1:length) |
---|
[3494] | 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 | ! |
---|
[3841] | 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 | |
---|
[3494] | 436 | ENDIF |
---|
| 437 | ! |
---|
| 438 | !-- Allocate and read array for variable data. |
---|
| 439 | ALLOCATE( var(1:ns(f)) ) |
---|
| 440 | |
---|
[3496] | 441 | READ( file_id_in ) var |
---|
[3494] | 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. |
---|
[3755] | 449 | ! CALL posix_ftell( file_id_in, filepos(f) ) |
---|
| 450 | filepos(f) = ftell( file_id_in ) |
---|
[3494] | 451 | |
---|
[3496] | 452 | CLOSE ( file_id_in ) |
---|
[3494] | 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 | !------------------------------------------------------------------------------! |
---|
[3496] | 478 | SUBROUTINE surface_output_parin |
---|
[3494] | 479 | |
---|
| 480 | IMPLICIT NONE |
---|
| 481 | |
---|
| 482 | INTEGER(iwp) :: file_id_parin = 90 |
---|
| 483 | |
---|
[3841] | 484 | NAMELIST /surface_output/ convert_average_data, convert_to_netcdf, & |
---|
| 485 | cycle_number, num_pe, path, run |
---|
[3494] | 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 | |
---|
[3496] | 498 | END SUBROUTINE surface_output_parin |
---|
[3494] | 499 | |
---|
| 500 | !------------------------------------------------------------------------------! |
---|
| 501 | ! Description: |
---|
| 502 | ! ------------ |
---|
| 503 | !> This subroutine creates the filename string of the treated binary file. |
---|
| 504 | !------------------------------------------------------------------------------! |
---|
[3496] | 505 | SUBROUTINE surface_output_create_file_string |
---|
[3494] | 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 | |
---|
[3496] | 557 | END SUBROUTINE surface_output_create_file_string |
---|
[3494] | 558 | |
---|
[3841] | 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 ) |
---|
[3494] | 642 | |
---|
[3841] | 643 | IMPLICIT NONE |
---|
[3494] | 644 | |
---|
[3841] | 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 |
---|
[3494] | 653 | |
---|
[3841] | 654 | END SUBROUTINE handle_error |
---|
| 655 | |
---|
| 656 | END PROGRAM surface_output_to_vtk_netcdf |
---|
| 657 | |
---|
| 658 | |
---|
| 659 | |
---|