Changeset 4763


Ignore:
Timestamp:
Oct 30, 2020 12:06:47 PM (3 years ago)
Author:
suehring
Message:

Revision of setting-up virtual measurements. virtual_measurement_mod: Simplified input of coordinates. All coordinate and height arrays are now 1D, independent on featuretype. This allows easier usage also for other campaign data sets independent of UC2; Avoid type conversion from 64 to 32 bit when output the data; Increase dimension size of sample-variable string (sometimes more than 100 variables are listed for heavy measured locations); Remove quantities that cannot be sampled; | Further, revision of the script palm_cvd to convert measurement coordinates from UC2 data standard towards a PALM readable format: Convert trajectory and timeseriesProfile coordinates into 1-D coordinates equivalent to timeseries coordiates. This simplifies processing in PALM and makes the virtual-measurement module also applicable to other campaigns; Check automatically for data organization (stored in subdirectories or not).

Location:
palm/trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SCRIPTS/palm_cvd

    r4758 r4763  
    2626# -----------------
    2727# $Id$
     28# - Check automatically for data organization (stored in subdirectories or not)
     29# - Convert trajectory and timeseriesProfile coordinates into 1-D coordinates
     30#   equivalent to timeseries coordiates. This simplifies processing in PALM and
     31#   makes the virtual-measurement module also applicable to other campaigns.
     32#
     33# 4758 2020-10-26 13:03:52Z suehring
    2834# In order to do not omit observations that are on the same site but have different
    2935# coordinates or feature-types, process all files rather than only one and omit
     
    5460import os
    5561import numpy as np
     62import time
    5663
    5764
     
    187194name_featuretype   = "featureType"
    188195name_ts            = "timeSeries"
     196name_tspr          = "timeSeriesProfile"
    189197name_traj          = "trajectory"
    190198name_ntime         = "ntime"
     
    196204name_eutm          = "E_UTM"
    197205name_nutm          = "N_UTM"
    198 name_hao           = "height_above_origin"
     206name_hao           = "height"
    199207name_station_h     = "station_h"
    200208name_z             = "z"
     
    271279if ( input_from_observations ):
    272280
    273    # Run loop over all subdirectories, detect the files and extract a list of sites.
     281   # Run loop over all listed input data. Depending on the data set, this could be
     282   # a list of files or a list of subdirectories.
    274283   # This is done to reduce the number of virtual measurements in the model. Each
    275284   # virtual measurement has an overhead and consumes memory.
    276285   sites = []
    277286   input_files = []
     287   input_files_orig = []
    278288   for dirname in list_input_data:
    279289       data_file = data_path + dirname
    280290
    281        # Directory may contain various file versions.
    282        # Take the one with highest cycle number.
    283        highest_cycle_nr = 0
    284        for filename in os.listdir(data_file):
    285           start_seq = len( filename ) - 6
    286           end_seq   = len( filename ) - 3
    287           if int( filename[start_seq:end_seq] ) > highest_cycle_nr:
    288              highest_cycle_nr = int(filename[start_seq:end_seq])
    289              latest_file      = filename
     291       if ( os.path.isdir(data_file) == True ):
     292          # Directory may contain various file versions.
     293          # Take the one with highest cycle number.
     294          highest_cycle_nr = 0
     295          for filename in os.listdir(data_file):
     296             start_seq = len( filename ) - 6
     297             end_seq   = len( filename ) - 3
     298             if int( filename[start_seq:end_seq] ) > highest_cycle_nr:
     299                highest_cycle_nr = int(filename[start_seq:end_seq])
     300                latest_file      = filename
     301          input_file = data_file + "/" + latest_file
     302          input_file_orig = latest_file
     303       else:
     304          input_file = data_file
     305          input_file_orig = dirname
    290306
    291307       # Open the NetCDF file
    292        input_file = data_file + "/" + latest_file
    293308       ncfile_in = Dataset( input_file, "r", format="NETCDF4", encoding='ascii')
    294309       input_files.append(input_file)
    295 
    296    # Define a nested list of default variables that shall be measured. Based on this list,
    297    # the final number of measured variables is determined.
    298    measured_variables_all_sites = [ ['u', 'v', 'w', 'theta', 'hus'] for var in range(0, len(input_files))]
    299 
    300    # Run loop over all subdirectories that contain observational data
    301    for counter, file in enumerate(input_files):
    302 
    303       # Open the NetCDF input file
    304       input_file = input_files[counter]
     310       input_files_orig.append(input_file_orig)
     311
     312#  Gather all files according to their feature type and all sites for the respective feature type
     313   files_traj  = []
     314   files_ts    = []
     315   files_tspr  = []
     316   sites_traj  = []
     317   sites_ts    = []
     318   sites_tspr  = []
     319   for input_file in input_files:
    305320      ncfile_in = Dataset( input_file, "r", format="NETCDF4", encoding='ascii' )
    306321
    307       # Determine index for the treated site
    308       num_vmeas = input_files.index( input_file ) + 1
    309       print( counter, num_vmeas )
    310 
    311       # Read global attributes and write them immediately into the output file
    312322      for att in ncfile_in.ncattrs():
    313323         if ( att == name_featuretype ):
    314324            feature = ncfile_in.getncattr(att)
    315          if ( att == name_datacontent ):
    316             content = ncfile_in.getncattr(att)
    317325         if ( att == name_site ):
    318326            site = ncfile_in.getncattr(att)
    319327
    320          if ( att in atts_float ):
    321             ncfile_out.setncattr( att + str(num_vmeas), np.double(ncfile_in.getncattr(att)) )
    322          else:
    323             ncfile_out.setncattr( att + str(num_vmeas), ncfile_in.getncattr(att) )
    324 
    325       #timeSeries
    326       if ( feature == name_ts ):
    327          ntime = len( ncfile_in.dimensions[name_ntime]   )
    328          nstat = len( ncfile_in.dimensions[name_station] )
    329          ncfile_out.createDimension( name_ntime   + str(num_vmeas), ntime )
    330          ncfile_out.createDimension( name_station + str(num_vmeas), nstat )
    331 
    332       #trajectory
    333       elif ( feature == name_traj ):
    334          ntime = len( ncfile_in.dimensions[name_ntime]   )
    335          ntraj = len( ncfile_in.dimensions[name_traj_dim] )
    336          ncfile_out.createDimension( name_ntime    + str(num_vmeas), ntime )
    337          ncfile_out.createDimension( name_traj_dim + str(num_vmeas), ntraj )
    338 
    339       #timeseriesProfile
     328      if ( feature == name_traj ):
     329         files_traj.append(input_file)
     330      elif ( feature == name_ts ):
     331         files_ts.append(input_file)
    340332      else:
    341          ntime = len( ncfile_in.dimensions[name_ntime]   )
    342          nstat = len( ncfile_in.dimensions[name_station] )
    343          nz    = len( ncfile_in.dimensions[name_nz]      )
    344          ncfile_out.createDimension( name_ntime   + str(num_vmeas), ntime )
    345          ncfile_out.createDimension( name_station + str(num_vmeas), nstat )
    346          ncfile_out.createDimension( name_nz      + str(num_vmeas), nz    )
    347 
    348       for var in ncfile_in.variables.keys():
    349          if ( var in dims_out ):
    350             # Create a variable and write it to file after it is read. In order to
    351             # avoid fill values in the dimensions, these are converted to zero
    352             # before written to file. Depending on the featureType of the measurement,
    353             # the array shape is different. For more informations, please see
    354             # [UC]2 data standard.
    355             # Timeseries
    356             if ( feature == name_ts  ):
    357                temp_ts = ncfile_out.createVariable( var + str(num_vmeas), float, \
    358                                                     name_station + str(num_vmeas))
    359                temp_ts[:] = np.nan_to_num( ncfile_in.variables[var][:] )
    360 
    361             # Trajectories
    362             elif ( feature == name_traj ):
    363                # @note: If there are files where station_h is present although featureType is
    364                #        trajectory, station_h must not be read.
    365                if var != name_station_h:
    366                   temp_traj = ncfile_out.createVariable( var + str(num_vmeas), float, \
    367                                                          ( name_traj_dim + str(num_vmeas), \
    368                                                            name_ntime + str(num_vmeas) ) )
    369                   temp_traj[:,:] = np.nan_to_num( ncfile_in.variables[var][:,:] )
    370 
    371             # TimeseriesProfiles
    372             else:
    373                if ( var == 'z' ):
    374                   temp_pr = ncfile_out.createVariable( var + str(num_vmeas), float, \
    375                                                       ( name_station + str(num_vmeas), \
    376                                                         name_nz + str(num_vmeas) ) )
    377                   temp_pr[:] = np.nan_to_num( ncfile_in.variables[var][:,0,:] )
    378                else:
    379                   temp_pr = ncfile_out.createVariable( var + str(num_vmeas), float, \
    380                                                        name_station + str(num_vmeas))
    381                   temp_pr[:] = np.nan_to_num( ncfile_in.variables[var][:] )
    382 
    383       # Search for variables to be measured. In case the variable isn't already defined,
    384       # append the variable to the list.
    385       for var in ncfile_in.variables.keys():
    386          if ( var not in non_measurable_vars  and  \
    387               var not in vars_default         and  \
    388               var not in measured_variables_all_sites[input_files.index( input_file )] ):
    389 
    390             measured_variables_all_sites[input_files.index( input_file )].append(var)
    391 
    392       # Close the NetCDF input file
     333         files_tspr.append(input_file)
     334
     335      if ( feature == name_traj  and  site not in sites_traj ):
     336         sites_traj.append(site)
     337      if ( feature == name_ts  and  site not in sites_ts ):
     338         sites_ts.append(site)
     339      if ( feature == name_tspr  and  site not in sites_tspr ):
     340         sites_tspr.append(site)
     341
    393342      ncfile_in.close()
    394343
    395    # After variables are gathered and dimensions / attributes are already written to file,
    396    # the list of measured variables is written to file.
    397    for site in input_files:
    398 
    399       num_vmeas = input_files.index( site ) + 1
    400 
    401       ncfile_out.createDimension( "nvar"+ str(num_vmeas), \
    402                                   len( measured_variables_all_sites[input_files.index( site )] ) )
    403 
    404       measured = ncfile_out.createVariable( 'measured_variables' + str(num_vmeas), 'S1', \
    405                                             ("nvar" + str(num_vmeas), "string_len")) # must be NC_CHAR
    406 
    407       for counter, meas in enumerate( measured_variables_all_sites[input_files.index( site )] ):
    408          measured[counter] = stringtochar( np.array( meas,"S%s"%(max_string_len) ) )
    409 
    410       # Check if any of the measured variables is a soil variable. Set flag accordingly.
     344   for input_file in files_traj:
     345      print( "traj", input_file )
     346   for site in sites_traj:
     347      print( "traj", site )
     348
     349   for site in sites_ts:
     350      print( "ts", site )
     351   for site in sites_tspr:
     352      print( "tspr", site )
     353     
     354   for file in files_tspr:
     355      print( "tspr", file )
     356   counter_id = 1
     357   for site_traj in sites_traj:
     358      # For the given site already define the featureTpye and site
     359      ncfile_out.setncattr( name_featuretype + str(counter_id), name_traj )
     360      ncfile_out.setncattr( name_site + str(counter_id), site_traj )
     361
     362      # Define the number of coordinates for the site
     363      num_coords = 0
     364
     365      e_utm_traj = np.array([])
     366      n_utm_traj = np.array([])
     367      h_traj     = np.array([])
     368      measured_variables = ['u', 'v', 'w', 'theta', 'hus']
     369      for input_file in files_traj:
     370         ncfile_in = Dataset( input_file, "r", format="NETCDF4", encoding='ascii' )
     371         for att in ncfile_in.ncattrs():
     372            if ( att == name_site ):
     373               site = ncfile_in.getncattr(att)
     374
     375         if ( site == site_traj ):
     376
     377            for att in ncfile_in.ncattrs():
     378               if ( att == name_orig_x ):
     379                  orig_x = ncfile_in.getncattr(att)
     380               if ( att == name_orig_y ):
     381                  orig_y = ncfile_in.getncattr(att)
     382               if ( att == name_orig_z ):
     383                  orig_z = ncfile_in.getncattr(att)
     384
     385            ntime = len( ncfile_in.dimensions[name_ntime]    )
     386            ntraj = len( ncfile_in.dimensions[name_traj_dim] )
     387
     388            num_coords += ntime * ntraj
     389#           Gather UTM and height coordinates and merge them into one array. Further, gather
     390#           the variables that shall be sampled. Coordinates are checked to for NaN values and
     391#           are tranformed to arithmetric numbers. Further, 2D input array is transformed into
     392#           a 1D array.
     393            for var in ncfile_in.variables.keys():
     394               if ( var in dims_out  and  var == name_eutm ):
     395                  e_utm_traj = np.append(e_utm_traj, np.nan_to_num( ncfile_in.variables[var][:,:] ).flatten())
     396                  #e_utm_traj.append( np.nan_to_num( ncfile_in.variables[var][:,:] ).flatten() )
     397               if ( var in dims_out  and  var == name_nutm ):
     398                  n_utm_traj = np.append(n_utm_traj, np.nan_to_num( ncfile_in.variables[var][:,:] ).flatten())
     399               if ( var in dims_out  and  var == name_hao ):
     400                  h_traj = np.append(h_traj, np.nan_to_num( ncfile_in.variables[var][:,:] ).flatten())
     401
     402               if ( var not in non_measurable_vars  and  \
     403                    var not in vars_default         and  \
     404                    var not in measured_variables ):
     405                  measured_variables.append(var)
     406
     407         ncfile_in.close()
     408
     409#     After all files for the current site are processed, write the origin-coordinates for x,y,z
     410      ncfile_out.setncattr( name_orig_x + str(counter_id), orig_x )
     411      ncfile_out.setncattr( name_orig_y + str(counter_id), orig_y )
     412      ncfile_out.setncattr( name_orig_z + str(counter_id), orig_z )
     413#     Create the dimensions
     414      ncfile_out.createDimension( name_station + str(counter_id), num_coords )
     415
     416      temp_traj = ncfile_out.createVariable( name_eutm + str(counter_id), float, name_station + str(counter_id) )
     417      temp_traj[:] = e_utm_traj
     418
     419      temp_traj = ncfile_out.createVariable( name_nutm + str(counter_id), float, name_station + str(counter_id) )
     420      temp_traj[:] = n_utm_traj
     421
     422      temp_traj = ncfile_out.createVariable( name_hao  + str(counter_id), float, name_station + str(counter_id) )
     423      temp_traj[:] = h_traj
     424
     425#     Check if any of the measured variables is a soil variable. Set flag accordingly.
    411426      soil = False
    412       for var in measured_variables_all_sites[input_files.index( site )]:
     427      for var in measured_variables:
    413428         if ( var in soil_vars ):
    414429            soil = True
    415 
    416       # Write soil flag
    417       ncfile_out.setncattr( name_soil_sampling + str( num_vmeas), np.int8(soil) )
    418 
    419 
    420 
    421     # Store the number of observational sites
    422    num_sites += len( input_files )
     430#     Write soil flag
     431      ncfile_out.setncattr( name_soil_sampling + str( counter_id), np.int8(soil) )
     432
     433#     Create dimension for sample-variable string
     434      ncfile_out.createDimension( "nvar"+ str(counter_id), len( measured_variables ) )
     435
     436      measured_var = ncfile_out.createVariable( 'measured_variables' + str(counter_id), 'S1', \
     437                                                ("nvar" + str(counter_id), "string_len") ) # must be NC_CHAR
     438
     439      # Write the variables to the file
     440      for counter, meas in enumerate( measured_variables ):
     441         measured_var[counter] = stringtochar( np.array( meas,"S%s"%(max_string_len) ) )
     442
     443#     Increment the counter
     444      counter_id += 1
     445   
     446   
     447   for site_tspr in sites_tspr:
     448      # For the given site already define the featureTpye and site
     449      ncfile_out.setncattr( name_featuretype + str(counter_id), name_tspr )
     450      ncfile_out.setncattr( name_site + str(counter_id), site_tspr )
     451
     452      # Define the number of coordinates for the site
     453      num_coords = 0
     454      e_utm_tspr     = np.array([])
     455      n_utm_tspr     = np.array([])
     456      station_h_tspr = np.array([])
     457      z_tspr         = np.array([])
     458
     459      measured_variables = ['u', 'v', 'w', 'theta', 'hus']
     460      for input_file in files_tspr:
     461         ncfile_in = Dataset( input_file, "r", format="NETCDF4", encoding='ascii' )
     462         for att in ncfile_in.ncattrs():
     463            if ( att == name_site ):
     464               site = ncfile_in.getncattr(att)
     465
     466         if ( site == site_tspr ):
     467            for att in ncfile_in.ncattrs():
     468               if ( att == name_orig_x ):
     469                  orig_x = ncfile_in.getncattr(att)
     470               if ( att == name_orig_y ):
     471                  orig_y = ncfile_in.getncattr(att)
     472               if ( att == name_orig_z ):
     473                  orig_z = ncfile_in.getncattr(att)
     474
     475            nstation = len( ncfile_in.dimensions[name_station] )
     476            ntime    = len( ncfile_in.dimensions[name_ntime]   )
     477            nz       = len( ncfile_in.dimensions[name_nz]   )
     478
     479            num_coords += nstation * ntime * nz
     480#           Gather UTM and height coordinates and merge them into one array. Further, gather
     481#           the variables that shall be sampled. Coordinates are checked to for NaN values and
     482#           are tranformed to arithmetric numbers. Further, 2D input array is transformed into
     483#           a 1D array.
     484            for var in ncfile_in.variables.keys():
     485               tspr_tmp1 = np.zeros((nstation))
     486               tspr_tmp2 = np.zeros((ntime*nz))
     487               if ( var in dims_out  and  var == name_eutm ):
     488                  tspr_tmp1 = np.nan_to_num( ncfile_in.variables[var][:] )
     489                  for ns in range(0,int(nstation)):
     490                     tspr_tmp2[:] = tspr_tmp1[ns]
     491                     e_utm_tspr = np.append(e_utm_tspr, tspr_tmp2)
     492               if ( var in dims_out  and  var == name_nutm ):
     493                  tspr_tmp1 = np.nan_to_num( ncfile_in.variables[var][:] )
     494                  for ns in range(0,int(nstation)):
     495                     tspr_tmp2[:] = tspr_tmp1[ns]
     496                     n_utm_tspr = np.append(n_utm_tspr, tspr_tmp2)
     497               if ( var in dims_out  and  var == name_z ):
     498                  z_tspr_tmp = np.nan_to_num( ncfile_in.variables[var][:,:,:] )
     499                  z_tspr = np.append(z_tspr, np.concatenate( z_tspr_tmp ))
     500               if ( var in dims_out  and  var == name_station_h ):
     501                  tspr_tmp1 = np.nan_to_num( ncfile_in.variables[var][:] )
     502                  for ns in range(0,int(nstation)):
     503                     tspr_tmp2[:] = tspr_tmp1[ns]
     504                     station_h_tspr = np.append(station_h_tspr, tspr_tmp2)
     505
     506               if ( var not in non_measurable_vars  and  \
     507                    var not in vars_default         and  \
     508                    var not in measured_variables ):
     509                  measured_variables.append(var)
     510
     511         ncfile_in.close()
     512
     513#     After all files for the current site are processed, write the origin-coordinates for x,y,z
     514      ncfile_out.setncattr( name_orig_x + str(counter_id), orig_x )
     515      ncfile_out.setncattr( name_orig_y + str(counter_id), orig_y )
     516      ncfile_out.setncattr( name_orig_z + str(counter_id), orig_z )
     517#     Create the dimensions
     518      ncfile_out.createDimension( name_station + str(counter_id), num_coords )
     519
     520      temp_tspr = ncfile_out.createVariable( name_eutm      + str(counter_id), float, name_station + str(counter_id) )
     521      temp_tspr[:] = e_utm_tspr
     522
     523      temp_tspr = ncfile_out.createVariable( name_nutm      + str(counter_id), float, name_station + str(counter_id) )
     524      temp_tspr[:] = n_utm_tspr
     525
     526      temp_tspr = ncfile_out.createVariable( name_z         + str(counter_id), float, name_station + str(counter_id) )
     527      temp_tspr[:] = z_tspr
     528
     529      temp_tspr = ncfile_out.createVariable( name_station_h + str(counter_id), float, name_station + str(counter_id) )
     530      temp_tspr[:] = station_h_tspr
     531
     532#     Check if any of the measured variables is a soil variable. Set flag accordingly.
     533      soil = False
     534      for var in measured_variables:
     535         if ( var in soil_vars ):
     536            soil = True
     537#     Write soil flag
     538      ncfile_out.setncattr( name_soil_sampling + str( counter_id), np.int8(soil) )
     539
     540#     Create dimension for sample-variable string
     541      ncfile_out.createDimension( "nvar"+ str(counter_id), len( measured_variables ) )
     542
     543      measured_var = ncfile_out.createVariable( 'measured_variables' + str(counter_id), 'S1', \
     544                                                ("nvar" + str(counter_id), "string_len") ) # must be NC_CHAR
     545
     546      # Write the variables to the file
     547      for counter, meas in enumerate( measured_variables ):
     548         measured_var[counter] = stringtochar( np.array( meas,"S%s"%(max_string_len) ) )
     549
     550#     Increment the counter
     551      counter_id += 1
     552   
     553   
     554   for site_ts in sites_ts:
     555      # For the given site already define the featureTpye and site
     556      ncfile_out.setncattr( name_featuretype + str(counter_id), name_ts )
     557      ncfile_out.setncattr( name_site + str(counter_id), site_ts )
     558
     559      # Define the number of coordinates for the site
     560      num_coords = 0
     561      e_utm_ts     = np.array([])
     562      n_utm_ts     = np.array([])
     563      station_h_ts = np.array([])
     564      z_ts         = np.array([])
     565     
     566      measured_variables = ['u', 'v', 'w', 'theta', 'hus']
     567      for input_file in files_ts:
     568         ncfile_in = Dataset( input_file, "r", format="NETCDF4", encoding='ascii' )
     569         for att in ncfile_in.ncattrs():
     570            if ( att == name_site ):
     571               site = ncfile_in.getncattr(att)
     572
     573         if ( site == site_ts ):
     574
     575            for att in ncfile_in.ncattrs():
     576               if ( att == name_orig_x ):
     577                  orig_x = ncfile_in.getncattr(att)
     578               if ( att == name_orig_y ):
     579                  orig_y = ncfile_in.getncattr(att)
     580               if ( att == name_orig_z ):
     581                  orig_z = ncfile_in.getncattr(att)
     582
     583            nstation = len( ncfile_in.dimensions[name_station]    )
     584            num_coords += nstation
     585#           Gather UTM and height coordinates and merge them into one array. Further, gather
     586#           the variables that shall be sampled. Coordinates are checked to for NaN values and
     587#           are tranformed to arithmetric numbers.
     588            for var in ncfile_in.variables.keys():
     589               if ( var in dims_out  and  var == name_eutm ):
     590                  e_utm_ts = np.append(e_utm_ts, np.nan_to_num( ncfile_in.variables[var][:] ))
     591               if ( var in dims_out  and  var == name_nutm ):
     592                  n_utm_ts = np.append(n_utm_ts, np.nan_to_num( ncfile_in.variables[var][:] ))
     593               if ( var in dims_out  and  var == name_z ):
     594                  z_ts = np.append(z_ts, np.nan_to_num( ncfile_in.variables[var][:] ))
     595               if ( var in dims_out  and  var == name_station_h ):
     596                  station_h_ts = np.append(station_h_ts, np.nan_to_num( ncfile_in.variables[var][:] ))
     597
     598               if ( var not in non_measurable_vars  and  \
     599                    var not in vars_default         and  \
     600                    var not in measured_variables ):
     601                  measured_variables.append(var)
     602
     603         ncfile_in.close()
     604
     605#     After all files for the current site are processed, write the origin-coordinates for x,y,z
     606      ncfile_out.setncattr( name_orig_x + str(counter_id), orig_x )
     607      ncfile_out.setncattr( name_orig_y + str(counter_id), orig_y )
     608      ncfile_out.setncattr( name_orig_z + str(counter_id), orig_z )
     609#     Create the dimensions
     610      ncfile_out.createDimension( name_station + str(counter_id), num_coords )
     611
     612      temp_ts = ncfile_out.createVariable( name_eutm      + str(counter_id), float, name_station + str(counter_id) )
     613      temp_ts[:] = e_utm_ts
     614
     615      temp_ts = ncfile_out.createVariable( name_nutm      + str(counter_id), float, name_station + str(counter_id) )
     616      temp_ts[:] = n_utm_ts
     617
     618      temp_ts = ncfile_out.createVariable( name_z         + str(counter_id), float, name_station + str(counter_id) )
     619      temp_ts[:] = z_ts
     620
     621      temp_ts = ncfile_out.createVariable( name_station_h + str(counter_id), float, name_station + str(counter_id) )
     622      temp_ts[:] = station_h_ts
     623
     624#     Check if any of the measured variables is a soil variable. Set flag accordingly.
     625      soil = False
     626      for var in measured_variables:
     627         if ( var in soil_vars ):
     628            soil = True
     629#     Write soil flag
     630      ncfile_out.setncattr( name_soil_sampling + str( counter_id), np.int8(soil) )
     631
     632#     Create dimension for sample-variable string
     633      ncfile_out.createDimension( "nvar"+ str(counter_id), len( measured_variables ) )
     634
     635      measured_var = ncfile_out.createVariable( 'measured_variables' + str(counter_id), 'S1', \
     636                                                ("nvar" + str(counter_id), "string_len") ) # must be NC_CHAR
     637
     638      # Write the variables to the file
     639      for counter, meas in enumerate( measured_variables ):
     640         measured_var[counter] = stringtochar( np.array( meas,"S%s"%(max_string_len) ) )
     641
     642#     Increment the counter
     643      counter_id += 1
     644
     645# Store the number of observational sites
     646   num_sites = len( sites_traj ) + len( sites_ts ) + len( sites_tspr )
    423647
    424648
     
    426650# are possible.
    427651if ( custom_coordinates ):
    428 
    429    count_site = num_sites + 1
     652   num_sites = counter_id - 1
    430653   for coord in coordinates:
    431654      # Define mandatory attributes
    432       ncfile_out.setncattr( name_featuretype + str(count_site),  \
     655      ncfile_out.setncattr( name_featuretype + str(counter_id),  \
    433656                            name_ts )
    434       ncfile_out.setncattr( name_site        + str(count_site),  \
    435                             "custom"         + str(count_site - num_sites) )
    436       ncfile_out.setncattr( name_orig_x      + str(count_site),  \
     657      ncfile_out.setncattr( name_site        + str(counter_id),  \
     658                            "custom"         + str(counter_id - num_sites) )
     659      ncfile_out.setncattr( name_orig_x      + str(counter_id),  \
    437660                            coord[0] )
    438       ncfile_out.setncattr( name_orig_y      + str(count_site),  \
     661      ncfile_out.setncattr( name_orig_y      + str(counter_id),  \
    439662                            coord[1] )
    440       ncfile_out.setncattr( name_orig_z      + str(count_site),  \
     663      ncfile_out.setncattr( name_orig_z      + str(counter_id),  \
    441664                            0.0 )
    442665
     
    444667      ntime = 1
    445668      nstat = 1
    446       ncfile_out.createDimension( name_ntime   + str(count_site), ntime )
    447       ncfile_out.createDimension( name_station + str(count_site), nstat )
     669      ncfile_out.createDimension( name_ntime   + str(counter_id), ntime )
     670      ncfile_out.createDimension( name_station + str(counter_id), nstat )
    448671
    449672      # Define coordinate variables
    450       temp_ts = ncfile_out.createVariable( name_eutm      + str(count_site), \
     673      temp_ts = ncfile_out.createVariable( name_eutm      + str(counter_id), \
    451674                                           float,                            \
    452                                            name_station   + str(count_site) )
     675                                           name_station   + str(counter_id) )
    453676      temp_ts[:] = np.array( coord[0] )
    454677
    455       temp_ts = ncfile_out.createVariable( name_nutm      + str(count_site), \
     678      temp_ts = ncfile_out.createVariable( name_nutm      + str(counter_id), \
    456679                                           float,                            \
    457                                            name_station   + str(count_site) )
     680                                           name_station   + str(counter_id) )
    458681      temp_ts[:] = np.array( coord[1] )
    459682
    460       temp_ts = ncfile_out.createVariable( name_z         + str(count_site), \
     683      temp_ts = ncfile_out.createVariable( name_z         + str(counter_id), \
    461684                                           float,                            \
    462                                            name_station   + str(count_site) )
     685                                           name_station   + str(counter_id) )
    463686      temp_ts[:] = np.array( coord[2] )
    464687
    465       temp_ts = ncfile_out.createVariable( name_station_h + str(count_site), \
     688      temp_ts = ncfile_out.createVariable( name_station_h + str(counter_id), \
    466689                                           float,                            \
    467                                            name_station   + str(count_site) )
     690                                           name_station   + str(counter_id) )
    468691      temp_ts[:] = np.array( 0.0 )
    469692
    470693
    471       count_site += 1
     694      counter_id += 1
    472695
    473696   # Reset counter variable
    474    count_site = num_sites + 1
     697   counter_id = num_sites + 1
    475698
    476699   # check if variables are prescribed. If so, prepare final output string
     
    491714                measured_variables.append(var)
    492715
    493          ncfile_out.createDimension( "nvar"+ str(count_site), \
     716         ncfile_out.createDimension( "nvar"+ str(counter_id), \
    494717                                     len( measured_variables ) )
    495718
    496          measured_var = ncfile_out.createVariable( 'measured_variables' + str(count_site), 'S1', \
    497                                                   ("nvar" + str(count_site), "string_len") ) # must be NC_CHAR
     719         measured_var = ncfile_out.createVariable( 'measured_variables' + str(counter_id), 'S1', \
     720                                                  ("nvar" + str(counter_id), "string_len") ) # must be NC_CHAR
    498721
    499722         # Write the variables to the file
     
    507730
    508731         # Write soil flag
    509          ncfile_out.setncattr( name_soil_sampling + str( count_site), np.int8(soil) )
     732         ncfile_out.setncattr( name_soil_sampling + str( counter_id), np.int8(soil) )
    510733
    511734         # Increment counter variable
    512          count_site += 1
     735         counter_id += 1
    513736
    514737         del ( measured_variables[:] )
    515738
    516739   # Add the number of customized sites.
    517    num_sites += int( number_positions )
     740   num_sites = counter_id - 1
    518741
    519742
  • palm/trunk/SOURCE/virtual_measurement_mod.f90

    r4752 r4763  
    2525! -----------------
    2626! $Id$
     27! - Simplified input of coordinates. All coordinate and height arrays are now 1D, independent on
     28!   featuretype. This allows easier usage also for other campaign data sets independent of UC2.
     29! - Avoid type conversion from 64 to 32 bit when output the data
     30! - Increase dimension size of sample-variable string (sometimes more than 100 variables are listed
     31!   for heavy measured locations)
     32! - Remove quantities that cannot be sampled
     33!
     34! 4752 2020-10-21 13:54:51Z suehring
    2735! Remove unnecessary output and merge output for quantities that are represented by the same
    2836! variable in PALM, e.g. surface temperature and brightness temperature
     
    191199        ONLY:  air_chemistry,                                                                      &
    192200               coupling_char,                                                                      &
     201               debug_output,                                                                       &
    193202               dz,                                                                                 &
    194203               end_time,                                                                           &
     
    305314       INTEGER(iwp) ::  ns              = 0  !< number of observation coordinates on subdomain, for atmospheric measurements
    306315       INTEGER(iwp) ::  ns_tot          = 0  !< total number of observation coordinates, for atmospheric measurements
    307        INTEGER(iwp) ::  n_tr_st              !< number of trajectories / station of a measurement
     316       INTEGER(iwp) ::  nst                  !< number of coordinate points given for a measurement
    308317       INTEGER(iwp) ::  nmeas                !< number of measured variables (atmosphere + soil)
    309318       INTEGER(iwp) ::  ns_soil         = 0  !< number of observation coordinates on subdomain, for soil measurements
     
    677686          output_variable%units         = 'm2 s-2'
    678687
    679        CASE ( 'plev' )
    680           output_variable%long_name     = 'air pressure'
    681           output_variable%standard_name = 'air_pressure'
    682           output_variable%units         = 'Pa'
    683 
    684688       CASE ( 'm_soil' )
    685689          output_variable%long_name     = 'soil moisture volumetric'
     
    731735          output_variable%units         = 'hPa'
    732736
    733        CASE ( 'pswrtg' )
    734           output_variable%long_name     = 'platform speed wrt ground'
    735           output_variable%standard_name = 'platform_speed_wrt_ground'
    736           output_variable%units         = 'm s-1'
    737 
    738        CASE ( 'pswrta' )
    739           output_variable%long_name     = 'platform speed wrt air'
    740           output_variable%standard_name = 'platform_speed_wrt_air'
    741           output_variable%units         = 'm s-1'
    742 
    743737       CASE ( 'pwv' )
    744738          output_variable%long_name     = 'water vapor partial pressure in air'
     
    746740          output_variable%units         = 'hPa'
    747741
    748        CASE ( 'ssdu' )
    749           output_variable%long_name     = 'duration of sunshine'
    750           output_variable%standard_name = 'duration_of_sunshine'
    751           output_variable%units         = 's'
    752 
    753742       CASE ( 't_lw' )
    754743          output_variable%long_name     = 'land water temperature'
     
    771760          output_variable%long_name     = 'upward kinematic latent heat flux in air'
    772761          output_variable%units         = 'g kg-1 m s-1'
    773 
    774        CASE ( 'zcb' )
    775           output_variable%long_name     = 'cloud base altitude'
    776           output_variable%standard_name = 'cloud_base_altitude'
    777           output_variable%units         = 'm'
    778 
    779        CASE ( 'zmla' )
    780           output_variable%long_name     = 'atmosphere boundary layer thickness'
    781           output_variable%standard_name = 'atmosphere_boundary_layer_thickness'
    782           output_variable%units         = 'm'
    783762
    784763       CASE ( 'mcpm1' )
     
    880859
    881860    CHARACTER(LEN=5)                  ::  dum                           !< dummy string indicating station id
    882     CHARACTER(LEN=100), DIMENSION(50) ::  measured_variables_file = ''  !< array with all measured variables read from NetCDF
    883     CHARACTER(LEN=100), DIMENSION(50) ::  measured_variables      = ''  !< dummy array with all measured variables that are allowed
    884 
    885     INTEGER(iwp) ::  dim_ntime  !< dimension size of time coordinate
     861    CHARACTER(LEN=100), DIMENSION(200) ::  measured_variables_file = ''  !< array with all measured variables read from NetCDF
     862    CHARACTER(LEN=100), DIMENSION(200) ::  measured_variables      = ''  !< dummy array with all measured variables that are allowed
     863
    886864    INTEGER(iwp) ::  i          !< grid index of virtual observation point in x-direction
    887865    INTEGER(iwp) ::  is         !< grid index of real observation point of the respective station in x-direction
     
    897875    INTEGER(iwp) ::  ll         !< running index over all measured variables in file
    898876    INTEGER(iwp) ::  m          !< running index for surface elements
    899     INTEGER(iwp) ::  n          !< running index over trajectory coordinates
    900877    INTEGER(iwp) ::  nofill     !< dummy for nofill return value (not used)
    901878    INTEGER(iwp) ::  ns         !< counter variable for number of observation points on subdomain
     
    922899    REAL(wp) ::  fill_zar    !< _FillValue for zar coordinate
    923900
    924     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  e_utm      !< easting UTM coordinate, temporary variable
    925     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  e_utm_tmp  !< EUTM coordinate before rotation
    926     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  height     !< observation height above ground (for trajectories)
    927     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  n_utm      !< northing UTM coordinate, temporary variable
    928     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  n_utm_tmp  !< NUTM coordinate before rotation
    929     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  station_h  !< station height above reference
    930     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  zar        !< observation height above reference
     901    REAL(wp), DIMENSION(:), ALLOCATABLE ::  e_utm      !< easting UTM coordinate, temporary variable
     902    REAL(wp), DIMENSION(:), ALLOCATABLE ::  e_utm_tmp  !< EUTM coordinate before rotation
     903    REAL(wp), DIMENSION(:), ALLOCATABLE ::  height     !< observation height above ground (for trajectories)
     904    REAL(wp), DIMENSION(:), ALLOCATABLE ::  n_utm      !< northing UTM coordinate, temporary variable
     905    REAL(wp), DIMENSION(:), ALLOCATABLE ::  n_utm_tmp  !< NUTM coordinate before rotation
     906    REAL(wp), DIMENSION(:), ALLOCATABLE ::  station_h  !< station height above reference
     907    REAL(wp), DIMENSION(:), ALLOCATABLE ::  zar        !< observation height above reference
     908
     909    IF ( debug_output )  CALL debug_message( 'vm_init', 'start' )
    931910#if defined( __netcdf )
    932911!
     
    10421021       IF ( vmea(l)%nmeas > 0 )  THEN
    10431022!
    1044 !--       For stationary measurements UTM coordinates are just one value and its dimension is
    1045 !--       "station", while for mobile measurements UTM coordinates are arrays depending on the
    1046 !--       number of trajectories and time, according to (UC)2 standard. First, inquire dimension
    1047 !--       length of the UTM coordinates.
    1048           IF ( vmea(l)%trajectory )  THEN
    1049 !
    1050 !--          For non-stationary measurements read the number of trajectories and the number of time
    1051 !--          coordinates.
    1052              CALL get_dimension_length( pids_id, vmea(l)%n_tr_st, "traj" // TRIM( dum ) )
    1053              CALL get_dimension_length( pids_id, dim_ntime, "ntime" // TRIM( dum ) )
    1054 !
    1055 !--       For stationary measurements the dimension for UTM is station and for the time-coordinate
    1056 !--       it is one.
    1057           ELSE
    1058              CALL get_dimension_length( pids_id, vmea(l)%n_tr_st, "station" // TRIM( dum ) )
    1059              dim_ntime = 1
    1060           ENDIF
    1061 !
    1062 !-        Allocate array which defines individual time/space frame for each trajectory or station.
    1063           ALLOCATE( vmea(l)%dim_t(1:vmea(l)%n_tr_st) )
     1023!--       For stationary and mobile measurements UTM coordinates are just one value and its dimension
     1024!--       is "station".
     1025          CALL get_dimension_length( pids_id, vmea(l)%nst, "station" // TRIM( dum ) )
    10641026!
    10651027!--       Allocate temporary arrays for UTM and height coordinates. Note, on file UTM coordinates
    10661028!--       might be 1D or 2D variables
    1067           ALLOCATE( e_utm(1:vmea(l)%n_tr_st,1:dim_ntime)       )
    1068           ALLOCATE( n_utm(1:vmea(l)%n_tr_st,1:dim_ntime)       )
    1069           ALLOCATE( station_h(1:vmea(l)%n_tr_st,1:dim_ntime)   )
    1070           ALLOCATE( zar(1:vmea(l)%n_tr_st,1:dim_ntime)         )
    1071           IF ( vmea(l)%trajectory )  ALLOCATE( height(1:vmea(l)%n_tr_st,1:dim_ntime) )
     1029          ALLOCATE( e_utm(1:vmea(l)%nst)       )
     1030          ALLOCATE( n_utm(1:vmea(l)%nst)       )
     1031          ALLOCATE( station_h(1:vmea(l)%nst)   )
     1032          ALLOCATE( zar(1:vmea(l)%nst)         )
     1033          IF ( vmea(l)%trajectory )  ALLOCATE( height(1:vmea(l)%nst) )
    10721034          e_utm     = 0.0_wp
    10731035          n_utm     = 0.0_wp
     
    10761038          IF ( vmea(l)%trajectory )  height = 0.0_wp
    10771039
    1078           ALLOCATE( e_utm_tmp(1:vmea(l)%n_tr_st,1:dim_ntime) )
    1079           ALLOCATE( n_utm_tmp(1:vmea(l)%n_tr_st,1:dim_ntime) )
     1040          ALLOCATE( e_utm_tmp(1:vmea(l)%nst) )
     1041          ALLOCATE( n_utm_tmp(1:vmea(l)%nst) )
    10801042!
    10811043!--       Read UTM and height coordinates for all trajectories and times. Note, in case
     
    10941056!--       UC2-standard.
    10951057          IF ( vmea(l)%trajectory )  THEN
    1096              CALL get_variable( pids_id, char_eutm   // TRIM( dum ), e_utm,  0, dim_ntime-1, 0,    &
    1097                                 vmea(l)%n_tr_st-1 )
    1098              CALL get_variable( pids_id, char_nutm   // TRIM( dum ), n_utm,  0, dim_ntime-1, 0,    &
    1099                                 vmea(l)%n_tr_st-1 )
    1100              CALL get_variable( pids_id, char_zar    // TRIM( dum ), zar,    0, dim_ntime-1, 0,    &
    1101                                 vmea(l)%n_tr_st-1 )
    1102              CALL get_variable( pids_id, char_height // TRIM( dum ), height, 0, dim_ntime-1, 0,    &
    1103                                 vmea(l)%n_tr_st-1 )
     1058             CALL get_variable( pids_id, char_eutm      // TRIM( dum ), e_utm(:)     )
     1059             CALL get_variable( pids_id, char_nutm      // TRIM( dum ), n_utm(:)     )
     1060             CALL get_variable( pids_id, char_height    // TRIM( dum ), height(:)    )
    11041061          ELSE
    1105              CALL get_variable( pids_id, char_eutm      // TRIM( dum ), e_utm(:,1)     )
    1106              CALL get_variable( pids_id, char_nutm      // TRIM( dum ), n_utm(:,1)     )
    1107              CALL get_variable( pids_id, char_station_h // TRIM( dum ), station_h(:,1) )
    1108              CALL get_variable( pids_id, char_zar       // TRIM( dum ), zar(:,1)       )
     1062             CALL get_variable( pids_id, char_eutm      // TRIM( dum ), e_utm(:)     )
     1063             CALL get_variable( pids_id, char_nutm      // TRIM( dum ), n_utm(:)     )
     1064             CALL get_variable( pids_id, char_station_h // TRIM( dum ), station_h(:) )
     1065             CALL get_variable( pids_id, char_zar       // TRIM( dum ), zar(:)       )
    11091066          ENDIF
    11101067
     
    11161073!
    11171074!--       Compute observation height above ground. Note, for trajectory measurements the height
    1118 !--       above the surface is actually stored in 'height'
     1075!--       above the surface is actually stored in 'height'.
    11191076          IF ( vmea(l)%trajectory )  THEN
    11201077             zar      = height
     
    11271084!--       on subdomain. This case, setup grid index space sample these quantities.
    11281085          meas_flag = 0
    1129           DO  t = 1, vmea(l)%n_tr_st
     1086          DO  t = 1, vmea(l)%nst
    11301087!
    11311088!--          First, compute relative x- and y-coordinates with respect to the lower-left origin of
     
    11331090!--          is not correct, the virtual sites will be misplaced. Further, in case of an rotated
    11341091!--          model domain, the UTM coordinates must also be rotated.
    1135              e_utm_tmp(t,1:dim_ntime) = e_utm(t,1:dim_ntime) - init_model%origin_x
    1136              n_utm_tmp(t,1:dim_ntime) = n_utm(t,1:dim_ntime) - init_model%origin_y
    1137              e_utm(t,1:dim_ntime) = COS( init_model%rotation_angle * pi / 180.0_wp )               &
    1138                                     * e_utm_tmp(t,1:dim_ntime)                                     &
     1092             e_utm_tmp(t) = e_utm(t) - init_model%origin_x
     1093             n_utm_tmp(t) = n_utm(t) - init_model%origin_y
     1094             e_utm(t) = COS( init_model%rotation_angle * pi / 180.0_wp )                           &
     1095                                    * e_utm_tmp(t)                                                 &
    11391096                                  - SIN( init_model%rotation_angle * pi / 180.0_wp )               &
    1140                                     * n_utm_tmp(t,1:dim_ntime)
    1141              n_utm(t,1:dim_ntime) = SIN( init_model%rotation_angle * pi / 180.0_wp )               &
    1142                                     * e_utm_tmp(t,1:dim_ntime)                                     &
     1097                                    * n_utm_tmp(t)
     1098             n_utm(t) = SIN( init_model%rotation_angle * pi / 180.0_wp )                           &
     1099                                    * e_utm_tmp(t)                                                 &
    11431100                                  + COS( init_model%rotation_angle * pi / 180.0_wp )               &
    1144                                     * n_utm_tmp(t,1:dim_ntime)
    1145 !
    1146 !--          Determine the individual time coordinate length for each station and trajectory. This
    1147 !--          is required as several stations and trajectories are merged into one file but they do
    1148 !--          not have the same number of points in time, hence, missing values may occur and cannot
    1149 !--          be processed further. This is actually a work-around for the specific (UC)2 dataset,
    1150 !--          but it won't harm anyway.
    1151              vmea(l)%dim_t(t) = 0
    1152              DO  n = 1, dim_ntime
    1153                 IF ( e_utm(t,n) /= fill_eutm  .AND.  n_utm(t,n) /= fill_nutm  .AND.                &
    1154                      zar(t,n)   /= fill_zar )  vmea(l)%dim_t(t) = n
    1155              ENDDO
     1101                                    * n_utm_tmp(t)
    11561102!
    11571103!--          Compute grid indices relative to origin and check if these are on the subdomain. Note,
     
    11701116             ENDIF
    11711117
    1172              DO  n = 1, vmea(l)%dim_t(t)
    1173                  is = INT( ( e_utm(t,n) + 0.5_wp * dx ) * ddx, KIND = iwp )
    1174                  js = INT( ( n_utm(t,n) + 0.5_wp * dy ) * ddy, KIND = iwp )
    1175 !
    1176 !--             Is the observation point on subdomain?
    1177                 on_pe = ( is >= nxl  .AND.  is <= nxr  .AND.  js >= nys  .AND.  js <= nyn )
    1178 !
    1179 !--             Check if observation coordinate is on subdomain.
    1180                 IF ( on_pe )  THEN
    1181 !
    1182 !--                Determine vertical index which corresponds to the observation height.
    1183                    ksurf = topo_top_ind(js,is,0)
    1184                    ks = MINLOC( ABS( zu - zw(ksurf) - zar(t,n) ), DIM = 1 ) - 1
    1185 !
    1186 !--                Set mask array at the observation coordinates. Also, flag the surrounding
    1187 !--                coordinate points, but first check whether the surrounding coordinate points are
    1188 !--                on the subdomain.
    1189                    kl = MERGE( ks-off_z, ksurf, ks-off_z >= nzb  .AND. ks-off_z >= ksurf )
    1190                    ku = MERGE( ks+off_z, nzt,   ks+off_z < nzt+1 )
    1191 
    1192                    DO  i = is-off, is+off
    1193                       DO  j = js-off, js+off
    1194                          DO  k = kl, ku
    1195                             meas_flag(k,j,i) = MERGE( IBSET( meas_flag(k,j,i), 0 ), 0,             &
    1196                                                       BTEST( wall_flags_total_0(k,j,i), 0 ) )
    1197                          ENDDO
     1118             is = INT( ( e_utm(t) + 0.5_wp * dx ) * ddx, KIND = iwp )
     1119             js = INT( ( n_utm(t) + 0.5_wp * dy ) * ddy, KIND = iwp )
     1120!
     1121!--          Is the observation point on subdomain?
     1122             on_pe = ( is >= nxl  .AND.  is <= nxr  .AND.  js >= nys  .AND.  js <= nyn )
     1123!
     1124!--          Check if observation coordinate is on subdomain.
     1125             IF ( on_pe )  THEN
     1126!
     1127!--             Determine vertical index which corresponds to the observation height.
     1128                ksurf = topo_top_ind(js,is,0)
     1129                ks = MINLOC( ABS( zu - zw(ksurf) - zar(t) ), DIM = 1 ) - 1
     1130!
     1131!--             Set mask array at the observation coordinates. Also, flag the surrounding
     1132!--             coordinate points, but first check whether the surrounding coordinate points are
     1133!--             on the subdomain.
     1134                kl = MERGE( ks-off_z, ksurf, ks-off_z >= nzb  .AND. ks-off_z >= ksurf )
     1135                ku = MERGE( ks+off_z, nzt,   ks+off_z < nzt+1 )
     1136
     1137                DO  i = is-off, is+off
     1138                   DO  j = js-off, js+off
     1139                      DO  k = kl, ku
     1140                         meas_flag(k,j,i) = MERGE( IBSET( meas_flag(k,j,i), 0 ), 0,                &
     1141                                                          BTEST( wall_flags_total_0(k,j,i), 0 ) )
    11981142                      ENDDO
    11991143                   ENDDO
    1200                 ENDIF
    1201              ENDDO
     1144                ENDDO
     1145             ENDIF
    12021146
    12031147          ENDDO
     
    13171261          IF ( ALLOCATED( n_utm_tmp ) )  DEALLOCATE( n_utm_tmp )
    13181262          IF ( ALLOCATED( n_utm )     )  DEALLOCATE( n_utm )
    1319           IF ( ALLOCATED( zar  )      )  DEALLOCATE( vmea(l)%dim_t )
    13201263          IF ( ALLOCATED( zar  )      )  DEALLOCATE( zar  )
    13211264          IF ( ALLOCATED( height )    )  DEALLOCATE( height )
     
    13361279    ns_all = 0
    13371280#if defined( __parallel )
    1338     CALL MPI_ALLREDUCE( vmea(:)%ns, ns_all(:), vmea_general%nvm,                                   &
    1339                         MPI_INTEGER, MPI_SUM, comm2d, ierr )
     1281    CALL MPI_ALLREDUCE( vmea(:)%ns, ns_all(:), vmea_general%nvm, MPI_INTEGER, MPI_SUM, comm2d, ierr )
    13401282#else
    13411283    ns_all(:) = vmea(:)%ns
     
    13461288    ns_all = 0
    13471289#if defined( __parallel )
    1348     CALL MPI_ALLREDUCE( vmea(:)%ns_soil, ns_all(:), vmea_general%nvm,                              &
    1349                         MPI_INTEGER, MPI_SUM, comm2d, ierr )
     1290    CALL MPI_ALLREDUCE( vmea(:)%ns_soil, ns_all(:), vmea_general%nvm, MPI_INTEGER, MPI_SUM, comm2d, ierr )
    13501291#else
    13511292    ns_all(:) = vmea(:)%ns_soil
     
    14031344
    14041345#endif
    1405 
     1346    IF ( debug_output )  CALL debug_message( 'vm_init', 'end' )
    14061347 END SUBROUTINE vm_init
    14071348
     
    19991940    REAL(wp), DIMENSION(:), ALLOCATABLE           ::  dum_lon                   !< transformed geographical coordinate (longitude)
    20001941    REAL(wp), DIMENSION(:), ALLOCATABLE           ::  oro_rel                   !< relative altitude of model surface
    2001     REAL(wp), DIMENSION(:), POINTER               ::  output_values_1d_pointer  !< pointer for 1d output array
    2002     REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET   ::  output_values_1d_target   !< target for 1d output array
    2003     REAL(wp), DIMENSION(:,:), POINTER             ::  output_values_2d_pointer  !< pointer for 2d output array
    2004     REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET ::  output_values_2d_target   !< target for 2d output array
     1942    REAL(sp), DIMENSION(:), POINTER               ::  output_values_1d_pointer  !< pointer for 1d output array
     1943    REAL(sp), DIMENSION(:), ALLOCATABLE, TARGET   ::  output_values_1d_target   !< target for 1d output array
     1944    REAL(sp), DIMENSION(:,:), POINTER             ::  output_values_2d_pointer  !< pointer for 2d output array
     1945    REAL(sp), DIMENSION(:,:), ALLOCATABLE, TARGET ::  output_values_2d_target   !< target for 2d output array
    20051946
    20061947    CALL cpu_log( log_point_s(26), 'VM output', 'start' )
     
    20271968
    20281969          return_value = dom_write_var( vmea(l)%nc_filename, 'E_UTM',                              &
    2029                                         values_realwp_1d = output_values_1d_pointer,               &
     1970                                        values_real32_1d = output_values_1d_pointer,               &
    20301971                                        bounds_start = (/vmea(l)%start_coord_a/),                  &
    20311972                                        bounds_end   = (/vmea(l)%end_coord_a  /) )
     
    20401981          output_values_1d_pointer => output_values_1d_target
    20411982          return_value = dom_write_var( vmea(l)%nc_filename, 'N_UTM',                              &
    2042                                         values_realwp_1d = output_values_1d_pointer,               &
     1983                                        values_real32_1d = output_values_1d_pointer,               &
    20431984                                        bounds_start = (/vmea(l)%start_coord_a/),                  &
    20441985                                        bounds_end   = (/vmea(l)%end_coord_a  /) )
     
    20662007          output_values_1d_pointer => output_values_1d_target
    20672008          return_value = dom_write_var( vmea(l)%nc_filename, 'lat',                                &
    2068                                         values_realwp_1d = output_values_1d_pointer,               &
     2009                                        values_real32_1d = output_values_1d_pointer,               &
    20692010                                        bounds_start = (/vmea(l)%start_coord_a/),                  &
    20702011                                        bounds_end   = (/vmea(l)%end_coord_a  /) )
     
    20732014          output_values_1d_pointer => output_values_1d_target
    20742015          return_value = dom_write_var( vmea(l)%nc_filename, 'lon',                                &
    2075                                         values_realwp_1d = output_values_1d_pointer,               &
     2016                                        values_real32_1d = output_values_1d_pointer,               &
    20762017                                        bounds_start = (/vmea(l)%start_coord_a/),                  &
    20772018                                        bounds_end   = (/vmea(l)%end_coord_a  /) )
     
    20892030          output_values_1d_pointer => output_values_1d_target
    20902031          return_value = dom_write_var( vmea(l)%nc_filename, 'z',                                  &
    2091                                         values_realwp_1d = output_values_1d_pointer,               &
     2032                                        values_real32_1d = output_values_1d_pointer,               &
    20922033                                        bounds_start = (/vmea(l)%start_coord_a/),                  &
    20932034                                        bounds_end   = (/vmea(l)%end_coord_a  /) )
     
    20982039          output_values_1d_pointer => output_values_1d_target
    20992040          return_value = dom_write_var( vmea(l)%nc_filename, 'station_h',                          &
    2100                                         values_realwp_1d = output_values_1d_pointer,               &
     2041                                        values_real32_1d = output_values_1d_pointer,               &
    21012042                                        bounds_start = (/vmea(l)%start_coord_a/),                  &
    21022043                                        bounds_end   = (/vmea(l)%end_coord_a  /) )
     
    21412082             output_values_1d_pointer => output_values_1d_target
    21422083             return_value = dom_write_var( vmea(l)%nc_filename, 'E_UTM_soil',                      &
    2143                                            values_realwp_1d = output_values_1d_pointer,            &
     2084                                           values_real32_1d = output_values_1d_pointer,            &
    21442085                                           bounds_start = (/vmea(l)%start_coord_s/),               &
    21452086                                           bounds_end   = (/vmea(l)%end_coord_s  /) )
     
    21542095             output_values_1d_pointer => output_values_1d_target
    21552096             return_value = dom_write_var( vmea(l)%nc_filename, 'N_UTM_soil',                      &
    2156                                            values_realwp_1d = output_values_1d_pointer,            &
     2097                                           values_real32_1d = output_values_1d_pointer,            &
    21572098                                           bounds_start = (/vmea(l)%start_coord_s/),               &
    21582099                                           bounds_end   = (/vmea(l)%end_coord_s  /) )
     
    21802121             output_values_1d_pointer => output_values_1d_target
    21812122             return_value = dom_write_var( vmea(l)%nc_filename, 'lat_soil',                        &
    2182                                            values_realwp_1d = output_values_1d_pointer,            &
     2123                                           values_real32_1d = output_values_1d_pointer,            &
    21832124                                           bounds_start = (/vmea(l)%start_coord_s/),               &
    21842125                                           bounds_end   = (/vmea(l)%end_coord_s  /) )
     
    21872128             output_values_1d_pointer => output_values_1d_target
    21882129             return_value = dom_write_var( vmea(l)%nc_filename, 'lon_soil',                        &
    2189                                            values_realwp_1d = output_values_1d_pointer,            &
     2130                                           values_real32_1d = output_values_1d_pointer,            &
    21902131                                           bounds_start = (/vmea(l)%start_coord_s/),               &
    21912132                                           bounds_end   = (/vmea(l)%end_coord_s  /) )
     
    22032144             output_values_1d_pointer => output_values_1d_target
    22042145             return_value = dom_write_var( vmea(l)%nc_filename, 'z_soil',                          &
    2205                                            values_realwp_1d = output_values_1d_pointer,            &
     2146                                           values_real32_1d = output_values_1d_pointer,            &
    22062147                                           bounds_start = (/vmea(l)%start_coord_s/),               &
    22072148                                           bounds_end   = (/vmea(l)%end_coord_s  /) )
     
    22122153             output_values_1d_pointer => output_values_1d_target
    22132154             return_value = dom_write_var( vmea(l)%nc_filename, 'station_h_soil',                  &
    2214                                            values_realwp_1d = output_values_1d_pointer,            &
     2155                                           values_real32_1d = output_values_1d_pointer,            &
    22152156                                           bounds_start = (/vmea(l)%start_coord_s/),               &
    22162157                                           bounds_end   = (/vmea(l)%end_coord_s  /) )
     
    22702211
    22712212             return_value = dom_write_var( vmea(l)%nc_filename, variable_name,                     &
    2272                                            values_realwp_2d = output_values_2d_pointer,            &
     2213                                           values_real32_2d = output_values_2d_pointer,            &
    22732214                                           bounds_start = (/vmea(l)%start_coord_s, t_ind/),        &
    22742215                                           bounds_end   = (/vmea(l)%end_coord_s, t_ind /) )
     
    22782219             output_values_2d_pointer => output_values_2d_target
    22792220             return_value = dom_write_var( vmea(l)%nc_filename, variable_name,                     &
    2280                                            values_realwp_2d = output_values_2d_pointer,            &
     2221                                           values_real32_2d = output_values_2d_pointer,            &
    22812222                                           bounds_start = (/vmea(l)%start_coord_s, t_ind/),        &
    22822223                                           bounds_end   = (/vmea(l)%end_coord_s, t_ind  /) )
     
    22912232
    22922233             return_value = dom_write_var( vmea(l)%nc_filename, variable_name,                     &
    2293                                            values_realwp_2d = output_values_2d_pointer,            &
     2234                                           values_real32_2d = output_values_2d_pointer,            &
    22942235                                           bounds_start = (/vmea(l)%start_coord_a, t_ind/),        &
    22952236                                           bounds_end   = (/vmea(l)%end_coord_a, t_ind/) )
     
    23002241             output_values_2d_pointer => output_values_2d_target
    23012242             return_value = dom_write_var( vmea(l)%nc_filename, variable_name,                     &
    2302                                            values_realwp_2d = output_values_2d_pointer,            &
     2243                                           values_real32_2d = output_values_2d_pointer,            &
    23032244                                           bounds_start = (/ vmea(l)%start_coord_a, t_ind /),      &
    23042245                                           bounds_end   = (/ vmea(l)%end_coord_a, t_ind /) )
Note: See TracChangeset for help on using the changeset viewer.