Changeset 3704


Ignore:
Timestamp:
Jan 29, 2019 7:51:41 PM (3 years ago)
Author:
suehring
Message:

Revision of virtual-measurement module and data output enabled. Further, post-processing tool added to merge distributed virtually sampled data and to output it into NetCDF files.

Location:
palm/trunk
Files:
3 added
9 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SCRIPTS/.palm.iofiles

    r3619 r3704  
    6060SURFACE_DATA_AV_BIN*       out:tr   *         $base_data/$run_identifier/OUTPUT         _av_surf   bin
    6161#
     62VIRTUAL_MEAS_BIN*          out   *            $base_data_output/$run_identifier/OUTPUT  _vmeas     bin
     63#
    6264DATA_1D_PR_NETCDF*         out:tr   *         $base_data/$run_identifier/OUTPUT         _pr        nc
    6365DATA_1D_SP_NETCDF*         out:tr   *         $base_data/$run_identifier/OUTPUT         _sp        nc
  • palm/trunk/SOURCE/check_open.f90

    r3655 r3704  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Open binary files for virtual measurements
    2323!
    2424! Former revisions:
     
    470470          OPEN ( 26, FILE='SURFACE_DATA_AV_BIN'//TRIM( coupling_char )//myid_char,  &
    471471                 FORM='UNFORMATTED', POSITION='APPEND' )
    472        
     472                 
     473       CASE ( 27 )
     474!
     475!--       Binary files for virtual measurement data
     476          OPEN ( 27, FILE='VIRTUAL_MEAS_BIN'//TRIM( coupling_char )//myid_char,&
     477                 FORM='UNFORMATTED', POSITION='APPEND' )
    473478
    474479       CASE ( 30 )
  • palm/trunk/SOURCE/module_interface.f90

    r3701 r3704  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Add last_actions for virtual measurements
    2323!
    2424! Former revisions:
     
    290290
    291291   USE virtual_measurement_mod,                                                &
    292        ONLY:  vm_parin,                                                        &
    293               vm_init
     292       ONLY:  vm_init,                                                         &
     293              vm_last_actions,                                                 &
     294              vm_parin
    294295
    295296   USE wind_turbine_model_mod,                                                 &
     
    12271228SUBROUTINE module_interface_last_actions
    12281229
    1229 
     1230   IF ( virtual_measurement )  CALL vm_last_actions
    12301231   IF ( user_module_enabled )  CALL user_last_actions
    12311232
  • palm/trunk/SOURCE/nesting_offl_mod.f90

    r3655 r3704  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Formatting adjustments
    2323!
    2424! Former revisions:
     
    623623       ENDIF
    624624#else
    625     u_ref  = u_ref_l
    626     v_ref  = v_ref_l
    627     IF ( humidity )       q_ref  = q_ref_l
    628     IF ( .NOT. neutral )  pt_ref = pt_ref_l
     625       u_ref  = u_ref_l
     626       v_ref  = v_ref_l
     627       IF ( humidity )       q_ref  = q_ref_l
     628       IF ( .NOT. neutral )  pt_ref = pt_ref_l
    629629#endif
    630630!
    631 !-- Average data. Note, reference profiles up to nzt are derived from lateral
    632 !-- boundaries, at the model top it is derived from the top boundary. Thus,
    633 !-- number of input data is different from nzb:nzt compared to nzt+1.
    634 !-- Derived from lateral boundaries.
    635     u_ref(nzb:nzt) = u_ref(nzb:nzt) / REAL( 2.0_wp * ( ny + 1 + nx     ),      &
    636                                             KIND = wp ) 
    637     v_ref(nzb:nzt) = v_ref(nzb:nzt) / REAL( 2.0_wp * ( ny   + nx + 1   ),      &
    638                                             KIND = wp )
    639     IF ( humidity )                                                            &
    640        q_ref(nzb:nzt) = q_ref(nzb:nzt)   / REAL( 2.0_wp * ( ny + 1 + nx + 1 ), &
    641                                                  KIND = wp )
    642     IF ( .NOT. neutral )                                                       &
    643        pt_ref(nzb:nzt) = pt_ref(nzb:nzt) / REAL( 2.0_wp * ( ny + 1 + nx + 1 ), &
    644                                            KIND = wp )
    645 !
    646 !-- Derived from top boundary.   
    647     u_ref(nzt+1) = u_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx     ), KIND = wp )
    648     v_ref(nzt+1) = v_ref(nzt+1) / REAL( ( ny     ) * ( nx + 1 ), KIND = wp )
    649     IF ( humidity )                                                            &
    650        q_ref(nzt+1) = q_ref(nzt+1)   / REAL( ( ny + 1 ) * ( nx + 1 ), KIND = wp )
    651     IF ( .NOT. neutral )                                                       &
    652        pt_ref(nzt+1) = pt_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx + 1 ), KIND = wp )
    653 !
    654 !-- Write onto init profiles, which are used for damping
    655     u_init = u_ref
    656     v_init = v_ref
    657     IF ( humidity      )  q_init  = q_ref
    658     IF ( .NOT. neutral )  pt_init = pt_ref
    659 !
    660 !-- Set bottom boundary condition
    661     IF ( humidity      )  q_init(nzb)  = q_init(nzb+1)
    662     IF ( .NOT. neutral )  pt_init(nzb) = pt_init(nzb+1)
    663 !
    664 !-- Further, adjust Rayleigh damping height in case of time-changing conditions.
    665 !-- Therefore, calculate boundary-layer depth first.
    666     CALL calc_zi
    667     CALL adjust_sponge_layer
    668    
    669 !
    670 !-- Update geostrophic wind components from dynamic input file.
    671     DO  k = nzb+1, nzt
    672        ug(k) = interpolate_in_time( nest_offl%ug(0,k), nest_offl%ug(1,k),      &
    673                                     fac_dt )
    674        vg(k) = interpolate_in_time( nest_offl%vg(0,k), nest_offl%vg(1,k),      &
    675                                     fac_dt )
    676     ENDDO
    677     ug(nzt+1) = ug(nzt)
    678     vg(nzt+1) = vg(nzt)
    679    
    680     CALL  cpu_log( log_point(58), 'offline nesting', 'stop' )
     631!--    Average data. Note, reference profiles up to nzt are derived from lateral
     632!--    boundaries, at the model top it is derived from the top boundary. Thus,
     633!--    number of input data is different from nzb:nzt compared to nzt+1.
     634!--    Derived from lateral boundaries.
     635       u_ref(nzb:nzt) = u_ref(nzb:nzt) / REAL( 2.0_wp * ( ny + 1 + nx     ),   &
     636                                               KIND = wp ) 
     637       v_ref(nzb:nzt) = v_ref(nzb:nzt) / REAL( 2.0_wp * ( ny   + nx + 1   ),   &
     638                                               KIND = wp )
     639       IF ( humidity )                                                         &
     640          q_ref(nzb:nzt) = q_ref(nzb:nzt)   / REAL( 2.0_wp *                   &
     641                                                          ( ny + 1 + nx + 1 ), &
     642                                                    KIND = wp )
     643       IF ( .NOT. neutral )                                                    &
     644          pt_ref(nzb:nzt) = pt_ref(nzb:nzt) / REAL( 2.0_wp *                   &
     645                                                          ( ny + 1 + nx + 1 ), &
     646                                              KIND = wp )
     647!
     648!--    Derived from top boundary.   
     649       u_ref(nzt+1) = u_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx     ), KIND = wp )
     650       v_ref(nzt+1) = v_ref(nzt+1) / REAL( ( ny     ) * ( nx + 1 ), KIND = wp )
     651       IF ( humidity )                                                         &
     652          q_ref(nzt+1) = q_ref(nzt+1)   / REAL( ( ny + 1 ) * ( nx + 1 ),       &
     653                                                KIND = wp )
     654       IF ( .NOT. neutral )                                                    &
     655          pt_ref(nzt+1) = pt_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx + 1 ),       &
     656                                                KIND = wp )
     657!
     658!--    Write onto init profiles, which are used for damping
     659       u_init = u_ref
     660       v_init = v_ref
     661       IF ( humidity      )  q_init  = q_ref
     662       IF ( .NOT. neutral )  pt_init = pt_ref
     663!
     664!--    Set bottom boundary condition
     665       IF ( humidity      )  q_init(nzb)  = q_init(nzb+1)
     666       IF ( .NOT. neutral )  pt_init(nzb) = pt_init(nzb+1)
     667!
     668!--    Further, adjust Rayleigh damping height in case of time-changing conditions.
     669!--    Therefore, calculate boundary-layer depth first.
     670       CALL calc_zi
     671       CALL adjust_sponge_layer
     672   
     673!
     674!--    Update geostrophic wind components from dynamic input file.
     675       DO  k = nzb+1, nzt
     676          ug(k) = interpolate_in_time( nest_offl%ug(0,k), nest_offl%ug(1,k),   &
     677                                       fac_dt )
     678          vg(k) = interpolate_in_time( nest_offl%vg(0,k), nest_offl%vg(1,k),   &
     679                                       fac_dt )
     680       ENDDO
     681       ug(nzt+1) = ug(nzt)
     682       vg(nzt+1) = vg(nzt)
     683   
     684       CALL  cpu_log( log_point(58), 'offline nesting', 'stop' )
    681685
    682686    END SUBROUTINE nesting_offl_bc
  • palm/trunk/SOURCE/netcdf_data_input_mod.f90

    r3655 r3704  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Interface for attribute input of 8-bit and 32-bit integer
    2323!
    2424! Former revisions:
     
    777777   
    778778    INTERFACE netcdf_data_input_att
    779        MODULE PROCEDURE netcdf_data_input_att_int
     779       MODULE PROCEDURE netcdf_data_input_att_int8
     780       MODULE PROCEDURE netcdf_data_input_att_int32
    780781       MODULE PROCEDURE netcdf_data_input_att_real
    781782       MODULE PROCEDURE netcdf_data_input_att_string
     
    11741175! Description:
    11751176! ------------
    1176 !> Read a global integer attribute
    1177 !------------------------------------------------------------------------------!
    1178     SUBROUTINE netcdf_data_input_att_int( val, search_string, id_mod,          &
    1179                                           input_file, global, openclose,       &
    1180                                           variable_name )
     1177!> Read a global 8-bit integer attribute
     1178!------------------------------------------------------------------------------!
     1179    SUBROUTINE netcdf_data_input_att_int8( val, search_string, id_mod,         &
     1180                                           input_file, global, openclose,      &
     1181                                           variable_name )
     1182
     1183       IMPLICIT NONE
     1184
     1185       CHARACTER(LEN=*) ::  search_string !< name of the attribue
     1186       
     1187       CHARACTER(LEN=*) ::  input_file    !< name of input file
     1188       CHARACTER(LEN=*) ::  openclose     !< string indicating whether NetCDF needs to be opend or closed
     1189       CHARACTER(LEN=*) ::  variable_name !< string indicating whether NetCDF needs to be opend or closed
     1190       
     1191       INTEGER(iwp) ::  id_mod   !< NetCDF id of input file
     1192       INTEGER(KIND=1) ::  val      !< value of the attribute
     1193       
     1194       LOGICAL ::  global        !< flag indicating a global or a variable's attribute
     1195
     1196#if defined ( __netcdf )
     1197!
     1198!--    Open file in read-only mode
     1199       IF ( openclose == 'open' )  THEN
     1200          CALL open_read_file( TRIM( input_file ) // TRIM( coupling_char ), &
     1201                                  id_mod )
     1202       ENDIF
     1203!
     1204!--    Read global attribute
     1205       IF ( global )  THEN
     1206          CALL get_attribute( id_mod, search_string, val, global )
     1207!
     1208!--    Read variable attribute
     1209       ELSE
     1210          CALL get_attribute( id_mod, search_string, val, global, variable_name )
     1211       ENDIF
     1212!
     1213!--    Finally, close input file
     1214       IF ( openclose == 'close' )  CALL close_input_file( id_mod )
     1215#endif           
     1216
     1217    END SUBROUTINE netcdf_data_input_att_int8
     1218   
     1219!------------------------------------------------------------------------------!
     1220! Description:
     1221! ------------
     1222!> Read a global 32-bit integer attribute
     1223!------------------------------------------------------------------------------!
     1224    SUBROUTINE netcdf_data_input_att_int32( val, search_string, id_mod,        &
     1225                                            input_file, global, openclose,     &
     1226                                            variable_name )
    11811227
    11821228       IMPLICIT NONE
     
    12141260#endif           
    12151261
    1216     END SUBROUTINE netcdf_data_input_att_int
     1262    END SUBROUTINE netcdf_data_input_att_int32
    12171263   
    12181264!------------------------------------------------------------------------------!
  • palm/trunk/SOURCE/radiation_model_mod.f90

    r3700 r3704  
    2323! Current revisions:
    2424! ------------------
    25 !
     25! Make variables that are sampled in virtual measurement module public
    2626!
    2727! Former revisions:
     
    12761276           idsvf, ndsvf, idcsf, ndcsf, kdcsf, pct,                             &
    12771277           radiation_interactions, startwall, startland, endland, endwall,     &
    1278            skyvf, skyvft, radiation_interactions_on, average_radiation
     1278           skyvf, skyvft, radiation_interactions_on, average_radiation,        &
     1279           rad_sw_in_diff, rad_sw_in_dir
    12791280
    12801281
  • palm/trunk/SOURCE/time_integration.f90

    r3684 r3704  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Data output for virtual measurements added
    2323!
    2424! Former revisions:
     
    632632               
    633633    USE virtual_measurement_mod,                                               &
    634         ONLY:  vm_sampling, vm_time_start
     634        ONLY:  vm_data_output, vm_sampling, vm_time_start
    635635
    636636    IMPLICIT NONE
     
    16121612!--    Take virtual measurements
    16131613       IF ( virtual_measurement  .AND.                                         &
    1614             vm_time_start <= time_since_reference_point )  CALL vm_sampling
     1614            vm_time_start <= time_since_reference_point )  THEN
     1615          CALL vm_sampling
     1616          CALL vm_data_output
     1617       ENDIF
    16151618
    16161619!
  • palm/trunk/SOURCE/urban_surface_mod.f90

    r3685 r3704  
    2323! Current revisions:
    2424! ------------------
    25 !
     25! make nzb_wall public, required for virtual-measurements
    2626!
    2727! Former revisions:
     
    11501150!-- Public parameters, constants and initial values
    11511151    PUBLIC usm_anthropogenic_heat, usm_material_model, usm_wall_mod,           &
    1152            usm_green_heat_model, building_pars,  &
    1153            nzt_wall, t_wall_h, t_wall_v,         &
     1152           usm_green_heat_model, building_pars,                                &
     1153           nzb_wall, nzt_wall, t_wall_h, t_wall_v,                             &
    11541154           t_window_h, t_window_v, building_type
    11551155
     
    74697469        IMPLICIT NONE
    74707470
    7471         INTEGER(iwp)                          :: i, j, k, l, m   !< running indices
     7471        INTEGER(iwp) :: i, j, k, l, m   !< running indices
    74727472       
    74737473        INTEGER(iwp) ::  i_off     !< offset to determine index of surface element, seen from atmospheric grid point, for x
  • palm/trunk/SOURCE/virtual_measurement_mod.f90

    r3665 r3704  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! - initialization revised
     23! - binary data output
     24! - list of allowed variables extended
    2325!
    2426! Former revisions:
    2527! -----------------
    2628! $Id$
    27 ! unused variables removed
    28 !
    29 ! 3522 2018-11-13 12:14:36Z suehring
    3029! Sampling of variables
    3130!
     
    4443!> The module acts as an interface between 'real-world' observations and
    4544!> model simulations. Virtual measurements will be taken in the model at the
    46 !> coordinates representative for the 'real-world' measurement positions.
     45!> coordinates representative for the 'real-world' observation coordinates.
    4746!> More precisely, coordinates and measured quanties will be read from a
    4847!> NetCDF file which contains all required information. In the model,
     
    5352!>
    5453!> @todo list_of_allowed variables needs careful checking
    55 !> @todo output (binary or NetCDF) needs to be implemented
    56 !> @todo clean-up anything from current test modus
    5754!> @todo Check if sign of surface fluxes for heat, radiation, etc., follows
    5855!>       the (UC)2 standard
     
    7168       
    7269    USE control_parameters,                                                    &
    73         ONLY:  air_chemistry, dz, humidity, neutral, message_string,           &
    74                virtual_measurement
     70        ONLY:  air_chemistry, dz, humidity, io_blocks, io_group, neutral,      &
     71               message_string, time_since_reference_point, virtual_measurement
    7572
    7673    USE cpulog,                                                                &
     
    8178
    8279    USE indices,                                                               &
    83         ONLY:  nzb, nzt, nxl, nxr, nys, nyn, nx, ny
     80        ONLY:  nzb, nzt, nxl, nxr, nys, nyn, nx, ny, wall_flags_0
    8481
    8582    USE kinds
     83   
     84    USE netcdf_data_input_mod,                                                 &
     85        ONLY:  init_model
     86       
     87    USE pegrid
     88   
     89    USE surface_mod,                                                           &
     90        ONLY:  surf_lsm_h, surf_usm_h
     91       
     92    USE land_surface_model_mod,                                                &
     93        ONLY:  nzb_soil, nzs, nzt_soil, zs, t_soil_h, m_soil_h
     94       
     95    USE radiation_model_mod
     96       
     97    USE urban_surface_mod,                                                     &
     98        ONLY:  nzb_wall, nzt_wall, t_wall_h
    8699
    87100
    88101    IMPLICIT NONE
     102   
     103    TYPE virt_general
     104       INTEGER(iwp) ::  id_vm     !< NetCDF file id for virtual measurements
     105       INTEGER(iwp) ::  nvm = 0   !< number of virtual measurements
     106    END TYPE virt_general
    89107
    90108    TYPE virt_mea
    91109   
    92        CHARACTER(LEN=100)  ::  feature_type  !< type of the measurement
    93        CHARACTER(LEN=100)  ::  site          !< name of the measurement site
     110       CHARACTER(LEN=100)  ::  feature_type      !< type of the measurement
     111       CHARACTER(LEN=100)  ::  filename_original !< name of the original file
     112       CHARACTER(LEN=100)  ::  site              !< name of the measurement site
    94113   
    95114       CHARACTER(LEN=10), DIMENSION(:), ALLOCATABLE ::  measured_vars_name !< name of the measured variables
    96115   
    97        INTEGER(iwp) ::  ns = 0 !< total number of observation points for a site on subdomain, i.e. sum of all trajectories
    98        INTEGER(iwp) ::  ntraj  !< number of trajectories of a measurement
    99        INTEGER(iwp) ::  nvar   !< number of measured variables
     116       INTEGER(iwp) ::  ns = 0          !< number of observation coordinates on subdomain, for atmospheric measurements
     117       INTEGER(iwp) ::  ns_tot = 0      !< total number of observation coordinates, for atmospheric measurements
     118       INTEGER(iwp) ::  ntraj           !< number of trajectories of a measurement
     119       INTEGER(iwp) ::  nvar            !< number of measured variables (atmosphere + soil)
     120       
     121       INTEGER(iwp) ::  ns_soil = 0     !< number of observation coordinates on subdomain, for soil measurements
     122       INTEGER(iwp) ::  ns_soil_tot = 0 !< total number of observation coordinates, for soil measurements
    100123       
    101124       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  dim_t !< number observations individual for each trajectory or station that are no _FillValues
    102125       
    103        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i  !< grid index for measurement position in x-direction
    104        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  j  !< grid index for measurement position in y-direction
    105        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k  !< grid index for measurement position in k-direction
     126       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i       !< grid index for measurement position in x-direction
     127       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  j       !< grid index for measurement position in y-direction
     128       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k       !< grid index for measurement position in k-direction
     129       
     130       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i_soil  !< grid index for measurement position in x-direction
     131       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  j_soil  !< grid index for measurement position in y-direction
     132       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k_soil  !< grid index for measurement position in k-direction
    106133           
    107134       LOGICAL ::  trajectory         = .FALSE. !< flag indicating that the observation is a mobile observation
    108135       LOGICAL ::  timseries          = .FALSE. !< flag indicating that the observation is a stationary point measurement
    109136       LOGICAL ::  timseries_profile  = .FALSE. !< flag indicating that the observation is a stationary profile measurement
     137       LOGICAL ::  soil_sampling      = .FALSE. !< flag indicating that soil state variables were sampled
    110138       
    111        REAL(wp) ::  fill_eutm                         !< fill value for UTM coordinates in case of missing values
    112        REAL(wp) ::  fill_nutm                         !< fill value for UTM coordinates in case of missing values
    113        REAL(wp) ::  fill_zag                          !< fill value for heigth coordinates in case of missing values
    114        REAL(wp) ::  fillout = -999.9                  !< fill value for output in case a observation is taken from inside a building
    115        REAL(wp) ::  origin_x_obs                      !< origin of the observation in UTM coordiates in x-direction
    116        REAL(wp) ::  origin_y_obs                      !< origin of the observation in UTM coordiates in y-direction
     139       REAL(wp) ::  fill_eutm          !< fill value for UTM coordinates in case of missing values
     140       REAL(wp) ::  fill_nutm          !< fill value for UTM coordinates in case of missing values
     141       REAL(wp) ::  fill_zag           !< fill value for heigth coordinates in case of missing values
     142       REAL(wp) ::  fillout = -999.9   !< fill value for output in case a observation is taken from inside a building
     143       REAL(wp) ::  origin_x_obs       !< origin of the observation in UTM coordiates in x-direction
     144       REAL(wp) ::  origin_y_obs       !< origin of the observation in UTM coordiates in y-direction
     145       
     146       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  z_ag           !< measurement height above ground level
     147       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  depth          !< measurement depth in soil
    117148             
    118        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  measured_vars  !< measured variables
     149       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  measured_vars       !< measured variables
     150       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  measured_vars_soil  !< measured variables
    119151       
    120152    END TYPE virt_mea
     
    122154    CHARACTER(LEN=5)  ::  char_eutm = "E_UTM"                      !< dimension name for UTM coordinate easting
    123155    CHARACTER(LEN=11) ::  char_feature = "featureType"             !< attribute name for feature type
     156   
     157    ! This need to generalized
     158    CHARACTER(LEN=8)  ::  char_filename = "filename"               !< attribute name for filename
     159    CHARACTER(LEN=11) ::  char_soil = "soil_sample"                !< attribute name for soil sampling indication
    124160    CHARACTER(LEN=10) ::  char_fillvalue = "_FillValue"            !< variable attribute name for _FillValue
    125161    CHARACTER(LEN=18) ::  char_mv = "measured_variables"           !< variable name for the array with the measured variable names
     
    133169    CHARACTER(LEN=10) ::  type_traj = 'trajectory'                 !< name of line measurements
    134170    CHARACTER(LEN=17) ::  type_tspr = 'timeSeriesProfile'          !< name of stationary profile measurements
     171   
     172    CHARACTER(LEN=6), DIMENSION(1:5) ::  soil_vars       = (/                  & !< list of soil variables
     173                            't_soil',                                          &
     174                            'm_soil',                                          &
     175                            'lwc   ',                                          &
     176                            'lwcs  ',                                          &
     177                            'smp   '                       /)
     178                           
     179    CHARACTER(LEN=10), DIMENSION(0:1,1:8) ::  chem_vars = RESHAPE( (/          &
     180                                              'mcpm1     ', 'kc_PM1    ',      &
     181                                              'mcpm2p5   ', 'kc_PM2.5  ',      &
     182                                              'mcpm25    ', 'kc_PM25   ',      &
     183                                              'mcpm10    ', 'kc_PM10   ',      &
     184                                              'mfno2     ', 'kc_NO2    ',      &
     185                                              'mfno      ', 'kc_NO     ',      &
     186                                              'tro3      ', 'kc_O3     ',      &
     187                                              'mfco      ', 'kc_CO     '       &
     188                                                                   /), (/ 2, 8 /) )
    135189!
    136190!-- MS: List requires careful revision!
    137     CHARACTER(LEN=10), DIMENSION(1:47), PARAMETER ::  list_allowed_variables = & !< variables that can be sampled in PALM
     191    CHARACTER(LEN=10), DIMENSION(1:54), PARAMETER ::  list_allowed_variables = & !< variables that can be sampled in PALM
    138192       (/ 'hfls      ',  & ! surface latent heat flux (W/m2)
    139193          'hfss      ',  & ! surface sensible heat flux (W/m2)
     
    143197          'mcpm1     ',  & ! mass concentration of PM1 (kg/m3)
    144198          'mcpm2p5   ',  & ! mass concentration of PM2.5 (kg/m3)
    145           'mcpm10    ',  & ! mass concentration of PM10 (kg/m3)
    146199          'mcpm10    ',  & ! mass concentration of PM10 (kg/m3)
    147200          'mcco      ',  & ! mass concentration of CO (kg/m3)
     
    173226          'tsoil     ',  & ! ? soil temperature - which depth?                                                               
    174227          'u         ',  & ! u-component
     228          'utheta    ',  & ! total eastward kinematic heat flux
    175229          'ua        ',  & ! eastward wind (is there any difference to u?)
    176230          'v         ',  & ! v-component
     231          'vtheta    ',  & ! total northward kinematic heat flux
    177232          'va        ',  & ! northward wind (is there any difference to v?)
    178233          'w         ',  & ! w-component
     234          'wtheta    ',  & ! total vertical kinematic heat flux
    179235          'rld       ',  & ! downward longwave radiative flux (W/m2)
    180236          'rlu       ',  & ! upnward longwave radiative flux (W/m2)
     
    182238          'rsu       ',  & ! upward shortwave radiative flux (W/m2)
    183239          'rsddif    ',  & ! downward shortwave diffuse radiative flux (W/m2)
    184           'rnds      '   & ! surface net downward radiative flux (W/m2)
     240          'rnds      ',  & ! surface net downward radiative flux (W/m2)
     241          't_soil    ',  &
     242          'm_soil    ',  &
     243          'lwc       ',  &
     244          'lwcs      ',  &
     245          'smp       '   &
    185246       /)
    186                                                              
    187     INTEGER(iwp) ::  id_vm                           !< NetCDF file id for virtual measurements
    188     INTEGER(iwp) ::  nvm = 0                         !< number of virtual measurements
    189 !    INTEGER(iwp) ::  observation_coverage_xy = 0     !< horizontal distance from the measurement point where observations should be taken in the surrounding
    190 !    INTEGER(iwp) ::  observation_coverage_z  = 0     !< vertical distance from the measurement point where observations should be taken in the surrounding
    191    
     247                                                           
     248   
     249    LOGICAL ::  global_attribute = .TRUE.         !< flag indicating a global attribute
     250    LOGICAL ::  init = .TRUE.                     !< flag indicating initialization of data output
    192251    LOGICAL ::  use_virtual_measurement = .FALSE. !< Namelist parameter
    193     LOGICAL ::  global_attribute = .TRUE.         !< flag indicating a global attribute
    194252   
    195253    REAL(wp) ::  vm_time_start = 0.0              !< time after virtual measurements should start
    196    
    197 
     254
     255    TYPE( virt_general )                        ::  vmea_general !< data structure which encompass general variables
    198256    TYPE( virt_mea ), DIMENSION(:), ALLOCATABLE ::  vmea !< virtual measurement data structure
    199    
    200257   
    201258    INTERFACE vm_check_parameters
     
    203260    END INTERFACE vm_check_parameters
    204261   
     262    INTERFACE vm_data_output
     263       MODULE PROCEDURE vm_data_output
     264    END INTERFACE vm_data_output
     265   
    205266    INTERFACE vm_init
    206267       MODULE PROCEDURE vm_init
    207268    END INTERFACE vm_init
    208269   
     270    INTERFACE vm_last_actions
     271       MODULE PROCEDURE vm_last_actions
     272    END INTERFACE vm_last_actions
     273   
    209274    INTERFACE vm_parin
    210275       MODULE PROCEDURE vm_parin
     
    221286!
    222287!-- Public interfaces
    223     PUBLIC  vm_check_parameters, vm_init, vm_parin, vm_sampling
     288    PUBLIC  vm_check_parameters, vm_data_output, vm_init, vm_last_actions,     &
     289            vm_parin, vm_sampling
    224290
    225291!
    226292!-- Public variables
    227     PUBLIC  vmea, vm_time_start
     293    PUBLIC  vmea, vmea_general, vm_time_start
    228294
    229295 CONTAINS
     
    251317!-- ToDo: Revise this later and remove this requirement.
    252318    IF ( virtual_measurement  .AND.  .NOT. input_pids_static )  THEN
    253        message_string = 'If virtual measurements are taken a static input ' // &
     319       message_string = 'If virtual measurements are taken, a static input ' //&
    254320                        'file is mandatory.'
    255321       CALL message( 'vm_check_parameters', 'PA0000', 1, 2, 0, 6, 0 )
     
    319385 
    320386    USE netcdf_data_input_mod,                                                 &
    321         ONLY:  input_file_vm, netcdf_data_input_get_dimension_length,          &
     387        ONLY:  init_model, input_file_vm,                                      &
     388               netcdf_data_input_get_dimension_length,                         &
    322389               netcdf_data_input_att, netcdf_data_input_var
    323390               
     
    328395   
    329396    CHARACTER(LEN=5)    ::  dum                !< dummy string indicate station id
     397    CHARACTER(LEN=5)    ::  dummy_read                !< dummy string indicate station id
    330398    CHARACTER(LEN=10), DIMENSION(50) ::  measured_variables_file = '' !< array with all measured variables read from NetCDF
    331399    CHARACTER(LEN=10), DIMENSION(50) ::  measured_variables      = '' !< dummy array with all measured variables that are allowed   
    332400   
    333 !    INTEGER(iwp) ::  dim_eutm  !< dimension size of UTM easting coordinate
    334 !    INTEGER(iwp) ::  dim_nutm  !< dimension size of UTM northing coordinate
     401    INTEGER(iwp) ::  dim_eutm  !< dimension size of UTM easting coordinate
     402    INTEGER(iwp) ::  dim_nutm  !< dimension size of UTM northing coordinate
    335403    INTEGER(iwp) ::  dim_ntime !< dimension size of time coordinate
    336 !    INTEGER(iwp) ::  dim_zag   !< dimension size of height coordinate
    337 !    INTEGER(iwp) ::  i         !< grid index of virtual observation point in x-direction
    338 !    INTEGER(iwp) ::  ii        !< running index over all coordinate points of a measurement
     404    INTEGER(iwp) ::  dim_zag   !< dimension size of height coordinate
     405    INTEGER(iwp) ::  i         !< grid index of virtual observation point in x-direction
     406    INTEGER(iwp) ::  ii        !< running index over all coordinate points of a measurement
    339407    INTEGER(iwp) ::  is        !< grid index of real observation point of the respective station in x-direction
    340 !    INTEGER(iwp) ::  j         !< grid index of observation point in x-direction
     408    INTEGER(iwp) ::  j         !< grid index of observation point in x-direction
    341409    INTEGER(iwp) ::  js        !< grid index of real observation point of the respective station in y-direction
    342 !    INTEGER(iwp) ::  k         !< grid index of observation point in x-direction
     410    INTEGER(iwp) ::  k         !< grid index of observation point in x-direction
    343411    INTEGER(iwp) ::  kl        !< lower vertical index of surrounding grid points of an observation coordinate
    344412    INTEGER(iwp) ::  ks        !< grid index of real observation point of the respective station in z-direction
     
    352420    INTEGER(iwp) ::  ns        !< counter variable for number of observation points on subdomain
    353421    INTEGER(iwp) ::  t         !< running index over number of trajectories
     422    INTEGER(iwp) ::  m
     423   
     424    INTEGER(KIND=1)::  soil_dum
     425   
     426    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ns_all !< dummy array used to sum-up the number of observation coordinates
    354427   
    355428    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  meas_flag !< mask array indicating measurement positions
     
    366439    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  z_ag  !< height coordinate relative to origin_z, temporary variable
    367440!
    368 !-- Obtain number of virtual measurement stations
    369     CALL netcdf_data_input_att( nvm, char_numstations, id_vm, input_file_vm,   &
     441!-- Obtain number of sites. Also, pass the 'open' string, in order to initially
     442!-- open the measurement driver.
     443    CALL netcdf_data_input_att( vmea_general%nvm, char_numstations,            &
     444                                vmea_general%id_vm, input_file_vm,             &
    370445                                global_attribute, 'open', '' )
    371446                               
    372 !     write(9,*) "num stationi", nvm
    373 !     flush(9)
    374 !
    375 !-- ALLOCATE data structure
    376     ALLOCATE( vmea(1:nvm) )
    377 !
    378 !-- Allocate flag array
     447!
     448!-- Allocate data structure which encompass all required information, such as
     449!-- grid points indicies, absolute UTM coordinates, the measured quantities,
     450!-- etc. .
     451    ALLOCATE( vmea(1:vmea_general%nvm) )
     452!
     453!-- Allocate flag array. This dummy array is used to identify grid points
     454!-- where virtual measurements should be taken. Please note, at least one
     455!-- ghost point is required, in order to include also the surrounding
     456!-- grid points of the original coordinate. 
    379457    ALLOCATE( meas_flag(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
    380458    meas_flag = 0
    381  
    382 !
    383 !-- Read station coordinates and further attributes.
    384 !-- Note all coordinates are in UTM coordinates.
    385     DO  l = 1, nvm
    386 !
    387 !--    Determine suffix which contains the ID
     459!
     460!-- Loop over all sites.
     461    DO  l = 1, vmea_general%nvm
     462!
     463!--    Determine suffix which contains the ID, ordered according to the number
     464!--    of measurements.
    388465       IF( l < 10 )  THEN
    389466          WRITE( dum, '(I1)')  l
     
    397474          WRITE( dum, '(I5)')  l
    398475       ENDIF
    399        
    400        CALL netcdf_data_input_att( vmea(l)%origin_x_obs, char_origx            &
    401                                    // TRIM( dum ), id_vm, '', global_attribute,&
     476!
     477!--    Read site coordinates (UTM).
     478       CALL netcdf_data_input_att( vmea(l)%origin_x_obs, char_origx //         &
     479                                   TRIM( dum ), vmea_general%id_vm, '',        &
     480                                   global_attribute, '', '' )
     481       CALL netcdf_data_input_att( vmea(l)%origin_y_obs, char_origy //         &
     482                                   TRIM( dum ), vmea_general%id_vm, '',        &
     483                                   global_attribute, '', '' )
     484!
     485!--    Read site name                 
     486       CALL netcdf_data_input_att( vmea(l)%site, char_site // TRIM( dum ),     &
     487                                   vmea_general%id_vm, '', global_attribute,   &
    402488                                   '', '' )
    403        CALL netcdf_data_input_att( vmea(l)%origin_y_obs, char_origy            &
    404                                    // TRIM( dum ), id_vm, '', global_attribute,&
     489!
     490!--    Read type of the measurement (trajectory, profile, timeseries).
     491       CALL netcdf_data_input_att( vmea(l)%feature_type, char_feature //       &
     492                                   TRIM( dum ), vmea_general%id_vm, '',        &
     493                                   global_attribute, '', '' )
     494!
     495!--    Read the name of the original file where observational data is stored.
     496       CALL netcdf_data_input_att( vmea(l)%filename_original, char_filename // &
     497                                   TRIM( dum ), vmea_general%id_vm, '',        &
     498                                   global_attribute, '', '' )
     499!
     500!--    Read a flag which indicates that also soil quantities are take at the
     501!--    respective site (is part of the virtual measurement driver). 
     502       CALL netcdf_data_input_att( soil_dum, char_soil // TRIM( dum ),         &
     503                                   vmea_general%id_vm, '', global_attribute,   &
    405504                                   '', '' )
    406        CALL netcdf_data_input_att( vmea(l)%site,         char_site             &
    407                                    // TRIM( dum ), id_vm, '', global_attribute,&
    408                                    '', '' )
    409        CALL netcdf_data_input_att( vmea(l)%feature_type, char_feature          &
    410                                    // TRIM( dum ), id_vm, '', global_attribute,&
    411                                    '', '' )
    412        
     505!
     506!--    Set flag for soil-sampling.
     507       IF ( soil_dum == 1 )  vmea(l)%soil_sampling = .TRUE.
    413508!
    414509!---   Set logicals depending on the type of the measurement
     
    419514       ELSEIF ( INDEX( vmea(l)%feature_type, type_traj ) /= 0 )  THEN
    420515          vmea(l)%trajectory        = .TRUE.
     516!
     517!--   Give error message in case the type matches non of the pre-defined types.
    421518       ELSE
    422 !
    423 !--       Give error message
    424519          message_string = 'Attribue featureType = ' //                        &
    425520                           TRIM( vmea(l)%feature_type ) //                     &
     
    428523       ENDIF
    429524!
    430 !--    Read string with all measured variables at this station
     525!--    Read string with all measured variables at this site
    431526       measured_variables_file = ''
    432527       CALL netcdf_data_input_var( measured_variables_file,                    &
    433                                    char_mv // TRIM( dum ), id_vm )
    434 !
    435 !--    Count the number of measured variables which match with the variables
    436 !--    which are allowed to be measured in PALM. Please note, for some
    437 !--    NetCDF interal reasons characters end with a NULL, i.e. also empty
    438 !--    characters contain a NULL. Therefore, check the strings for a Null to
    439 !--    get the correct character length in order to compare them with the list
    440 !--    of allowed variables.
    441        vmea(l)%nvar = 0
     528                                   char_mv // TRIM( dum ), vmea_general%id_vm )
     529!
     530!--    Count the number of measured variables. Only count variables that match
     531!--    with the allowed variables.
     532!--    Please note, for some NetCDF interal reasons characters end with a NULL,
     533!--    i.e. also empty characters contain a NULL. Therefore, check the strings
     534!--    for a NULL to get the correct character length in order to compare
     535!--    them with the list of allowed variables.
     536       vmea(l)%nvar   = 0
    442537       DO ll = 1, SIZE( measured_variables_file )
    443538          IF ( measured_variables_file(ll)(1:1) /= CHAR(0)  .AND.              &
     
    465560       ENDDO
    466561!
    467 !--    Allocate array for the measured variables names for the station l.
     562!--    Allocate array for the measured variables names for the respective site.
    468563       ALLOCATE( vmea(l)%measured_vars_name(1:vmea(l)%nvar) )
    469564
     
    474569!--    In case of chemistry, check if species is considered in the modelled
    475570!--    chemistry mechanism.
    476        IF ( air_chemistry )  THEN
    477           DO  ll = 1, vmea(l)%nvar
    478              chem_include = .FALSE.
    479              DO  n = 1, nspec
    480                 IF ( TRIM( vmea(l)%measured_vars_name(ll) ) ==                 &
    481                      TRIM( chem_species(n)%name ) )  chem_include = .TRUE.
    482              ENDDO
    483              IF ( .NOT. chem_include )  THEN
    484                 message_string = TRIM( vmea(l)%measured_vars_name(ll) ) //     &
    485                                  ' is not considered in the modelled '  //     &
    486                                  'chemistry mechanism'
    487                 CALL message( 'vm_init', 'PA0000', 0, 0, 0, 6, 0 )
    488              ENDIF
    489           ENDDO
    490        ENDIF
    491 !
    492 !--    For the actual measurement ID read the UTM coordinates. Based on these,
    493 !--    define the index space on each subdomain where measurements should be
    494 !--    taken. Note, the entire coordinate arrays will not be stored on data
    495 !--    type as this would exceed memory requirements, particularly for
    496 !--    trajectory measurements. If no variable will be virtually measured,
    497 !--    skip the reading.
     571!        IF ( air_chemistry )  THEN
     572!           DO  ll = 1, vmea(l)%nvar
     573!              chem_include = .FALSE.
     574!              DO  n = 1, nspec
     575!                 IF ( TRIM( vmea(l)%measured_vars_name(ll) ) ==                 &
     576!                      TRIM( chem_species(n)%name ) )  chem_include = .TRUE.
     577!              ENDDO
     578! !
     579! !--  Revise this. It should only check for chemistry variables and not for all!
     580!              IF ( .NOT. chem_include )  THEN
     581!                 message_string = TRIM( vmea(l)%measured_vars_name(ll) ) //     &
     582!                                  ' is not considered in the modelled '  //     &
     583!                                  'chemistry mechanism'
     584!                 CALL message( 'vm_init', 'PA0000', 0, 0, 0, 6, 0 )
     585!              ENDIF
     586!           ENDDO
     587!        ENDIF
     588!
     589!--    Read the UTM coordinates for the actual site. Based on the coordinates,
     590!--    define the grid-index space on each subdomain where virtual measurements
     591!--    should be taken. Note, the entire coordinate arrays will not be stored
     592!--    as this would exceed memory requirements, particularly for trajectory
     593!--    measurements.
    498594       IF ( vmea(l)%nvar > 0 )  THEN
    499595!
    500596!--       For stationary measurements UTM coordinates are just one value and
    501597!--       its dimension is "station", while for mobile measurements UTM
    502 !--       coordinates are arrays. First, inquire dimension length for
    503 !--       UTM coordinates.
     598!--       coordinates are arrays depending on the number of trajectories and
     599!--       time, according to (UC)2 standard. First, inquire dimension length
     600!--       of the UTM coordinates.
    504601          IF ( vmea(l)%trajectory )  THEN
    505602!
    506603!--          For non-stationary measurements read the number of trajectories
    507              CALL netcdf_data_input_get_dimension_length( id_vm,              &
     604!--          and the number of time coordinates.
     605             CALL netcdf_data_input_get_dimension_length( vmea_general%id_vm, &
    508606                                                          vmea(l)%ntraj,      &
    509607                                                          "traj" //           &
    510608                                                          TRIM( dum ) )
    511              CALL netcdf_data_input_get_dimension_length( id_vm, dim_ntime,   &
     609             CALL netcdf_data_input_get_dimension_length( vmea_general%id_vm, &
     610                                                          dim_ntime,          &
    512611                                                          "ntime" //          &
    513612                                                          TRIM( dum ) )
    514613!
    515 !--       For stationary measurements the dimension for UTM coordinates is 1
     614!--       For stationary measurements the dimension for UTM and time
     615!--       coordinates is 1.
    516616          ELSE
    517617             vmea(l)%ntraj  = 1
    518618             dim_ntime = 1
    519619          ENDIF
    520          
    521620!
    522621!-        Allocate array which defines individual time frame for each
    523 !--       trajectory or station
     622!--       trajectory or station.
    524623          ALLOCATE( vmea(l)%dim_t(1:vmea(l)%ntraj) )
    525624!
     
    530629          ALLOCATE( z_ag(1:vmea(l)%ntraj,1:dim_ntime)  )
    531630!
    532 !--       Read _FillValue attributes
     631!--       Read _FillValue attributes of the coordinate dimensions.
    533632          CALL netcdf_data_input_att( fill_eutm, char_fillvalue,               &
    534                                       id_vm, '', .NOT. global_attribute, '',   &
     633                                      vmea_general%id_vm, '',                  &
     634                                      .NOT. global_attribute, '',              &
    535635                                      char_eutm // TRIM( dum ) )
    536636          CALL netcdf_data_input_att( fill_nutm, char_fillvalue,               &
    537                                       id_vm, '', .NOT. global_attribute, '',   &
     637                                      vmea_general%id_vm, '',                  &
     638                                      .NOT. global_attribute, '',              &
    538639                                      char_nutm // TRIM( dum ) )
    539640          CALL netcdf_data_input_att( fill_zag, char_fillvalue,                &
    540                                       id_vm, '', .NOT. global_attribute, '',   &
     641                                      vmea_general%id_vm, '',                  &
     642                                      .NOT. global_attribute, '',              &
    541643                                      char_zag  // TRIM( dum ) )
    542644!
     
    544646!--       times.
    545647          IF ( vmea(l)%trajectory )  THEN
    546              CALL netcdf_data_input_var( e_utm, char_eutm // TRIM( dum ), id_vm,  &
     648             CALL netcdf_data_input_var( e_utm, char_eutm // TRIM( dum ),      &
     649                                         vmea_general%id_vm,                   &
    547650                                         0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
    548              CALL netcdf_data_input_var( n_utm, char_nutm // TRIM( dum ), id_vm,  &
     651             CALL netcdf_data_input_var( n_utm, char_nutm // TRIM( dum ),      &
     652                                         vmea_general%id_vm,                   &
    549653                                         0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
    550              CALL netcdf_data_input_var( z_ag, char_zag // TRIM( dum ), id_vm,  &
     654             CALL netcdf_data_input_var( z_ag, char_zag // TRIM( dum ),        &
     655                                         vmea_general%id_vm,                   &
    551656                                         0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
    552657          ELSE
    553              CALL netcdf_data_input_var( e_utm(1,:), char_eutm // TRIM( dum ), id_vm )
    554              CALL netcdf_data_input_var( n_utm(1,:), char_nutm // TRIM( dum ), id_vm )
    555              CALL netcdf_data_input_var( z_ag(1,:),  char_zag  // TRIM( dum ), id_vm )
    556           ENDIF
    557 !
    558 !-- For testing:
    559           e_utm = e_utm - e_utm(1,1)
    560           n_utm = n_utm - n_utm(1,1)
    561          
    562 !             
    563 !--       First, compute relative x- and y-coordinates with respect to the
    564 !--       lower-left origin of the model domain, which is the difference
    565 !--       betwen UTM coordinates. 
    566 !           e_utm(t,1:vmea(l)%dim_t(t)) = e_utm(t,1:vmea(l)%dim_t(t))            &
    567 !                                       - init_model%origin_x
    568 !           n_utm(t,1:vmea(l)%dim_t(t)) = n_utm(t,1:vmea(l)%dim_t(t))            &
    569 !                                       - init_model%origin_y
    570 
     658             CALL netcdf_data_input_var( e_utm(1,:), char_eutm // TRIM( dum ), &
     659                                         vmea_general%id_vm )
     660             CALL netcdf_data_input_var( n_utm(1,:), char_nutm // TRIM( dum ), &
     661                                         vmea_general%id_vm )
     662             CALL netcdf_data_input_var( z_ag(1,:),  char_zag  // TRIM( dum ), &
     663                                         vmea_general%id_vm )
     664          ENDIF         
    571665!
    572666!--       Based on UTM coordinates, check if the measurement station or parts
     
    575669          meas_flag = 0
    576670          DO  t = 1, vmea(l)%ntraj
     671!             
     672!--          First, compute relative x- and y-coordinates with respect to the
     673!--          lower-left origin of the model domain, which is the difference
     674!--          betwen UTM coordinates. Note, if the origin is not correct, the
     675!--          virtual sites will be misplaced.
     676             e_utm(t,1:dim_ntime) = e_utm(t,1:dim_ntime) - init_model%origin_x
     677             n_utm(t,1:dim_ntime) = n_utm(t,1:dim_ntime) - init_model%origin_y
    577678!
    578679!--          Determine the individual time coordinate length for each station and
     
    580681!--          are merged into one file but they do not have the same number of
    581682!--          points in time, hence, missing values may occur and cannot be
    582 !--          processed further.
     683!--          processed further. This is actually a work-around for the specific
     684!--          (UC)2 dataset, but it won't harm in anyway.
    583685             vmea(l)%dim_t(t) = 0
    584686             DO  n = 1, dim_ntime
     
    587689                     z_ag(t,n)  /= fill_zag )  vmea(l)%dim_t(t) = n
    588690             ENDDO
    589 
    590691!
    591692!--          Compute grid indices relative to origin and check if these are
     
    612713!--                surrounding coordinate points, but first check whether the
    613714!--                surrounding coordinate points are on the subdomain.
    614                    kl = MERGE( ks-1, ks, ks-1 >= nzb  .AND. ks-1 > ksurf )
    615                    ku = MERGE( ks+1, ks, ks+1 <= nzt+1 )
    616                  
    617                    meas_flag(kl:ku,js-1:js+1,is-1:is+1) =                      &
    618                                IBSET( meas_flag(kl:ku,js-1:js+1,is-1:is+1), 0 )
     715                   kl = MERGE( ks-1, ks, ks-1 >= nzb  .AND. ks-1 >= ksurf )
     716                   ku = MERGE( ks+1, ks, ks+1 < nzt+1 )
     717                 
     718                   DO  i = is-1, is+1
     719                      DO  j = js-1, js+1
     720                         DO  k = kl, ku
     721                            meas_flag(k,j,i) = MERGE(                          &
     722                                             IBSET( meas_flag(k,j,i), 0 ),     &
     723                                             0,                                &
     724                                             BTEST( wall_flags_0(k,j,i), 0 )   &
     725                                                    )
     726                         ENDDO
     727                      ENDDO
     728                   ENDDO
    619729                ENDIF
    620730             ENDDO
     
    622732          ENDDO
    623733!
    624 !--       Based on the flag array count the number of observation points.
     734!--       Based on the flag array count the number of of sampling coordinates.
     735!--       Please note, sampling coordinates in atmosphere and soil may be
     736!--       different, as within the soil all levels will be measured.           
     737!--       Hence, count individually. Start with atmoshere.
    625738          ns = 0
    626           DO  is = nxl-1, nxr+1
    627              DO  js = nys-1, nyn+1
    628                 DO  ks = nzb, nzt+1
    629                    ns = ns + MERGE( 1, 0, BTEST( meas_flag(ks,js,is), 0 ) )
     739          DO  i = nxl-1, nxr+1
     740             DO  j = nys-1, nyn+1
     741                DO  k = nzb, nzt+1
     742                   ns = ns + MERGE( 1, 0, BTEST( meas_flag(k,j,i), 0 ) )
    630743                ENDDO
    631744             ENDDO
    632745          ENDDO
     746         
    633747!
    634748!--       Store number of observation points on subdomain and allocate index
    635 !--       arrays.
     749!--       arrays as well as array containing height information.
    636750          vmea(l)%ns = ns
    637           ns = 0
    638751         
    639752          ALLOCATE( vmea(l)%i(1:vmea(l)%ns) )
    640753          ALLOCATE( vmea(l)%j(1:vmea(l)%ns) )
    641754          ALLOCATE( vmea(l)%k(1:vmea(l)%ns) )
     755          ALLOCATE( vmea(l)%z_ag(1:vmea(l)%ns) )         
    642756!
    643757!--       Based on the flag array store the grid indices which correspond to
    644758!--       the observation coordinates.
    645           DO  is = nxl-1, nxr+1
    646              DO  js = nys-1, nyn+1
    647                 DO  ks = nzb, nzt+1
    648                    IF ( BTEST( meas_flag(ks,js,is), 0 ) )  THEN
     759          ns = 0
     760          DO  i = nxl-1, nxr+1
     761             DO  j = nys-1, nyn+1
     762                DO  k = nzb, nzt+1
     763                   IF ( BTEST( meas_flag(k,j,i), 0 ) )  THEN
    649764                      ns = ns + 1
    650                       vmea(l)%i(ns) = is
    651                       vmea(l)%j(ns) = js
    652                       vmea(l)%k(ns) = ks
     765                      vmea(l)%i(ns) = i
     766                      vmea(l)%j(ns) = j
     767                      vmea(l)%k(ns) = k
     768                      vmea(l)%z_ag(ns)  = zu(k) -                              &
     769                                   zw(get_topography_top_index_ji( j, i, 's' ))
    653770                   ENDIF
    654771                ENDDO
    655772             ENDDO
    656773          ENDDO
     774!
     775!--       Same for the soil. Based on the flag array, count the number of
     776!--       sampling coordinates in soil. Sample at all soil levels in this case.
     777          IF ( vmea(l)%soil_sampling )  THEN
     778             DO  i = nxl, nxr
     779                DO  j = nys, nyn
     780                   IF ( ANY( BTEST( meas_flag(:,j,i), 0 ) ) )  THEN
     781                      IF ( surf_lsm_h%start_index(j,i) <=                      &
     782                           surf_lsm_h%end_index(j,i) )  THEN
     783                         vmea(l)%ns_soil = vmea(l)%ns_soil +                   &
     784                                                      nzt_soil - nzb_soil + 1
     785                      ENDIF
     786                      IF ( surf_usm_h%start_index(j,i) <=                      &
     787                           surf_usm_h%end_index(j,i) )  THEN
     788                         vmea(l)%ns_soil = vmea(l)%ns_soil +                   &
     789                                                      nzt_wall - nzb_wall + 1
     790                      ENDIF
     791                   ENDIF
     792                ENDDO
     793             ENDDO
     794          ENDIF         
     795!
     796!--       Allocate index arrays as well as array containing height information
     797!--       for soil.
     798          IF ( vmea(l)%soil_sampling )  THEN
     799             ALLOCATE( vmea(l)%i_soil(1:vmea(l)%ns_soil) )
     800             ALLOCATE( vmea(l)%j_soil(1:vmea(l)%ns_soil) )
     801             ALLOCATE( vmea(l)%k_soil(1:vmea(l)%ns_soil) )
     802             ALLOCATE( vmea(l)%depth(1:vmea(l)%ns_soil) )
     803          ENDIF     
     804!
     805!--       For soil, store the grid indices.
     806          ns = 0
     807          IF ( vmea(l)%soil_sampling )  THEN
     808             DO  i = nxl, nxr
     809                DO  j = nys, nyn
     810                   IF ( ANY( BTEST( meas_flag(:,j,i), 0 ) ) )  THEN
     811                      IF ( surf_lsm_h%start_index(j,i) <=                      &
     812                           surf_lsm_h%end_index(j,i) )  THEN
     813                         m = surf_lsm_h%start_index(j,i)
     814                         DO  k = nzb_soil, nzt_soil
     815                            ns = ns + 1
     816                            vmea(l)%i_soil(ns) = i
     817                            vmea(l)%j_soil(ns) = j
     818                            vmea(l)%k_soil(ns) = k
     819                            vmea(l)%depth(ns)  = zs(k)
     820                         ENDDO
     821                      ENDIF
     822                     
     823                      IF ( surf_usm_h%start_index(j,i) <=                      &
     824                           surf_usm_h%end_index(j,i) )  THEN
     825                         m = surf_usm_h%start_index(j,i)
     826                         DO  k = nzb_wall, nzt_wall
     827                            ns = ns + 1
     828                            vmea(l)%i_soil(ns) = i
     829                            vmea(l)%j_soil(ns) = j
     830                            vmea(l)%k_soil(ns) = k
     831                            vmea(l)%depth(ns)  = surf_usm_h%zw(k,m)
     832                         ENDDO
     833                      ENDIF
     834                   ENDIF
     835                ENDDO
     836             ENDDO
     837          ENDIF
     838!
     839!--       Allocate array to save the sampled values.
     840          ALLOCATE( vmea(l)%measured_vars(1:vmea(l)%ns,1:vmea(l)%nvar) )
    657841         
    658 !           write(9,*) l, "NS: ", ns
    659 !           IF ( ns > 0 )   THEN
    660 !              write(9,*) "i ", vmea(l)%i
    661 !              write(9,*) "j ", vmea(l)%j
    662 !              write(9,*) "k ", vmea(l)%k
    663 !              write(9,*)
    664 !           ENDIF
    665 !
    666 !--       Allocate array to save the sampled values.
    667 !--       Todo: Is it better to allocate for all variables at a station
    668 !--       and store all the values before writing, or sample the variables
    669 !--       directly in the data output?
    670           ALLOCATE( vmea(l)%measured_vars(1:vmea(l)%nvar,1:vmea(l)%ns) )
    671 !
    672 !--       Initialize with _FillValue
    673           vmea(l)%measured_vars(1:vmea(l)%nvar,1:vmea(l)%ns) = vmea(l)%fillout
     842          IF ( vmea(l)%soil_sampling )                                         &
     843             ALLOCATE( vmea(l)%measured_vars_soil(1:vmea(l)%ns_soil,           &
     844                                                  1:vmea(l)%nvar) )
     845!
     846!--       Initialize with _FillValues
     847          vmea(l)%measured_vars(1:vmea(l)%ns,1:vmea(l)%nvar) = vmea(l)%fillout
     848          IF ( vmea(l)%soil_sampling )                                         &
     849             vmea(l)%measured_vars_soil(1:vmea(l)%ns_soil,1:vmea(l)%nvar) =    &
     850                                                                vmea(l)%fillout
    674851!
    675852!--       Deallocate temporary coordinate arrays
     
    677854          IF ( ALLOCATED( n_utm ) )  DEALLOCATE( n_utm )
    678855          IF ( ALLOCATED( z_ag  ) )  DEALLOCATE( z_ag  )
     856          IF ( ALLOCATED( z_ag  ) )  DEALLOCATE( vmea(l)%dim_t )
    679857       ENDIF
    680858    ENDDO
    681 !     flush(9)
    682    
    683859!
    684860!-- Close input file for virtual measurements. Therefore, just call
    685861!-- the read attribute routine with the "close" option.
    686     CALL netcdf_data_input_att( nvm, char_numstations, id_vm, '',              &
     862    CALL netcdf_data_input_att( vmea_general%nvm, char_numstations,            &
     863                                vmea_general%id_vm, '',                        &
    687864                                global_attribute, 'close', '' )
     865!
     866!-- Sum-up the number of observation coordiates, for atmosphere first.
     867!-- This is actually only required for data output.
     868    ALLOCATE( ns_all(1:vmea_general%nvm) )
     869    ns_all = 0   
     870#if defined( __parallel )
     871    CALL MPI_ALLREDUCE( vmea(:)%ns, ns_all(:), vmea_general%nvm, MPI_INTEGER,  &
     872                        MPI_SUM, comm2d, ierr )
     873#else
     874    ns_all(:) = vmea(:)%ns
     875#endif
     876    vmea(:)%ns_tot = ns_all(:)
     877!
     878!-- Now for soil
     879    ns_all = 0   
     880#if defined( __parallel )
     881    CALL MPI_ALLREDUCE( vmea(:)%ns_soil, ns_all(:), vmea_general%nvm,          &
     882                        MPI_INTEGER, MPI_SUM, comm2d, ierr )
     883#else
     884    ns_all(:) = vmea(:)%ns_soil
     885#endif
     886    vmea(:)%ns_soil_tot = ns_all(:)
     887   
     888    DEALLOCATE( ns_all )
    688889!                               
    689890!-- Dellocate flag array
    690891    DEALLOCATE( meas_flag )
     892!
     893!-- Initialize binary data output of virtual measurements.
     894!-- Open binary output file.
     895    CALL check_open( 27 )
     896!
     897!-- Output header information.
     898    CALL vm_data_output
    691899       
    692900  END SUBROUTINE vm_init
     
    696904! Description:
    697905! ------------
     906!> Binary data output.
     907!------------------------------------------------------------------------------!
     908  SUBROUTINE vm_data_output
     909   
     910     USE pegrid
     911   
     912     IMPLICIT NONE
     913         
     914     INTEGER(iwp) ::  i         !< running index over IO blocks   
     915     INTEGER(iwp) ::  l         !< running index over all stations
     916     INTEGER(iwp) ::  n         !< running index over all measured variables at a station
     917!
     918!--  Header output on each PE
     919     IF ( init )  THEN
     920
     921        DO  i = 0, io_blocks-1
     922           IF ( i == io_group )  THEN
     923              WRITE ( 27 )  'number of measurements            '
     924              WRITE ( 27 )  vmea_general%nvm
     925
     926              DO  l = 1, vmea_general%nvm
     927                 WRITE ( 27 )  'site                              '
     928                 WRITE ( 27 )  vmea(l)%site
     929                 WRITE ( 27 )  'file                              '
     930                 WRITE ( 27 )  vmea(l)%filename_original
     931                 WRITE ( 27 )  'feature_type                      '
     932                 WRITE ( 27 )  vmea(l)%feature_type
     933                 WRITE ( 27 )  'origin_x_obs                      '
     934                 WRITE ( 27 )  vmea(l)%origin_x_obs
     935                 WRITE ( 27 )  'origin_y_obs                      '
     936                 WRITE ( 27 )  vmea(l)%origin_y_obs
     937                 WRITE ( 27 )  'total number of observation points'                               
     938                 WRITE ( 27 )  vmea(l)%ns_tot
     939                 WRITE ( 27 )  'number of measured variables      '
     940                 WRITE ( 27 )  vmea(l)%nvar
     941                 WRITE ( 27 )  'variables                         '
     942                 WRITE ( 27 )  vmea(l)%measured_vars_name(:)
     943                 WRITE ( 27 )  'number of observation points      '
     944                 WRITE ( 27 )  vmea(l)%ns
     945                 WRITE ( 27 )  'E_UTM                             '
     946                 WRITE ( 27 )  init_model%origin_x +                           &
     947                        REAL( vmea(l)%i(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dx
     948                 WRITE ( 27 )  'N_UTM                             '
     949                 WRITE ( 27 )  init_model%origin_y +                           &
     950                        REAL( vmea(l)%j(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dy
     951                 WRITE ( 27 )  'Z_AG                              '
     952                 WRITE ( 27 )  vmea(l)%z_ag(1:vmea(l)%ns)
     953                 WRITE ( 27 )  'soil sampling                     '
     954                 WRITE ( 27 )  MERGE( 'yes                               ',    &
     955                                      'no                                ',    &
     956                                      vmea(l)%soil_sampling )
     957 
     958                 IF ( vmea(l)%soil_sampling )  THEN                 
     959                    WRITE ( 27 )  'total number of soil points       '                               
     960                    WRITE ( 27 )  vmea(l)%ns_soil_tot
     961                    print*, "vmea(l)%ns_soil_tot", vmea(l)%ns_soil_tot
     962                    WRITE ( 27 )  'number of soil points             '
     963                    WRITE ( 27 )  vmea(l)%ns_soil
     964                    WRITE ( 27 )  'E_UTM soil                        '
     965                    WRITE ( 27 )  init_model%origin_x +                        &
     966                           REAL( vmea(l)%i_soil(1:vmea(l)%ns_soil) + 0.5_wp,   &
     967                                 KIND = wp ) * dx
     968                    WRITE ( 27 )  'N_UTM soil                        '
     969                    WRITE ( 27 )  init_model%origin_y +                        &
     970                           REAL( vmea(l)%j_soil(1:vmea(l)%ns_soil) + 0.5_wp,   &
     971                                 KIND = wp ) * dy
     972                    WRITE ( 27 )  'DEPTH                             '
     973                    WRITE ( 27 )  vmea(l)%depth(1:vmea(l)%ns_soil)
     974                 ENDIF
     975              ENDDO
     976
     977           ENDIF
     978        ENDDO
     979       
     980#if defined( __parallel )
     981        CALL MPI_BARRIER( comm2d, ierr )
     982#endif
     983!
     984!--     After header information is written, set control flag to .FALSE.
     985        init = .FALSE.
     986!
     987!--  Data output at each measurement timestep on each PE
     988     ELSE
     989        DO  i = 0, io_blocks-1
     990
     991           IF ( i == io_group )  THEN
     992              WRITE( 27 )  'output time                       '
     993              WRITE( 27 )  time_since_reference_point
     994              DO  l = 1, vmea_general%nvm
     995!
     996!--              Skip binary writing if no observation points are defined on PE
     997                 IF ( vmea(l)%ns < 1  .AND.  vmea(l)%ns_soil < 1)  CYCLE                 
     998                 DO  n = 1, vmea(l)%nvar
     999                    WRITE( 27 )  vmea(l)%measured_vars_name(n)
     1000                    IF ( vmea(l)%soil_sampling  .AND.                           &
     1001                         ANY( TRIM( vmea(l)%measured_vars_name(n))  ==          &
     1002                              soil_vars ) )  THEN                   
     1003                       WRITE( 27 )  vmea(l)%measured_vars_soil(:,n)
     1004                    ELSE
     1005                       WRITE( 27 )  vmea(l)%measured_vars(:,n)
     1006                    ENDIF
     1007                 ENDDO
     1008           
     1009              ENDDO
     1010           ENDIF
     1011        ENDDO
     1012#if defined( __parallel )
     1013        CALL MPI_BARRIER( comm2d, ierr )
     1014#endif
     1015     ENDIF
     1016 
     1017  END SUBROUTINE vm_data_output 
     1018 
     1019 
     1020!------------------------------------------------------------------------------!
     1021! Description:
     1022! ------------
     1023!> Write end-of-file statement as last action.
     1024!------------------------------------------------------------------------------!
     1025  SUBROUTINE vm_last_actions
     1026   
     1027     USE pegrid
     1028   
     1029     IMPLICIT NONE
     1030         
     1031     INTEGER(iwp) ::  i         !< running index over IO blocks   
     1032     INTEGER(iwp) ::  l         !< running index over all stations
     1033     INTEGER(iwp) ::  n         !< running index over all measured variables at a station
     1034 
     1035     DO  i = 0, io_blocks-1
     1036        IF ( i == io_group )  THEN
     1037           WRITE( 27 )  'EOF                               '
     1038        ENDIF
     1039     ENDDO
     1040#if defined( __parallel )
     1041        CALL MPI_BARRIER( comm2d, ierr )
     1042#endif
     1043!
     1044!--  Close binary file
     1045     CALL close_file( 27 )
     1046 
     1047  END SUBROUTINE vm_last_actions
     1048 
     1049!------------------------------------------------------------------------------!
     1050! Description:
     1051! ------------
    6981052!> Sampling of the actual quantities along the observation coordinates
    6991053!------------------------------------------------------------------------------!
     
    7141068     IMPLICIT NONE
    7151069     
    716 !     CHARACTER(LEN=10) ::  trimvar !< dummy for the measured variable name
     1070     CHARACTER(LEN=10) ::  trimvar !< dummy for the measured variable name
    7171071     
    718      INTEGER(iwp) ::  i  !< grid index in x-direction
    719      INTEGER(iwp) ::  j  !< grid index in y-direction
    720      INTEGER(iwp) ::  k  !< grid index in z-direction
    721      INTEGER(iwp) ::  l  !< running index over the number of stations
    722      INTEGER(iwp) ::  m  !< running index over all virtual observation coordinates
    723      INTEGER(iwp) ::  mm !< index of surface element which corresponds to the virtual observation coordinate
    724      INTEGER(iwp) ::  n  !< running index over all measured variables at a station
    725      INTEGER(iwp) ::  nn !< running index over the number of chemcal species
    726 !
    727 !--  Loop over all stations. For each possible variable loop over all
    728 !--  observation points
    729      DO  l = 1, nvm
    730 !
    731 !--     Loop over all measured variables at this station. Please note,
    732 !--     velocity components are interpolated onto scalar grid. 
     1072     INTEGER(iwp) ::  i         !< grid index in x-direction
     1073     INTEGER(iwp) ::  j         !< grid index in y-direction
     1074     INTEGER(iwp) ::  k         !< grid index in z-direction
     1075     INTEGER(iwp) ::  ind_chem  !< dummy index to identify chemistry variable and translate it from (UC)2 standard to interal naming
     1076     INTEGER(iwp) ::  l         !< running index over the number of stations
     1077     INTEGER(iwp) ::  m         !< running index over all virtual observation coordinates
     1078     INTEGER(iwp) ::  mm        !< index of surface element which corresponds to the virtual observation coordinate
     1079     INTEGER(iwp) ::  n         !< running index over all measured variables at a station
     1080     INTEGER(iwp) ::  nn        !< running index over the number of chemcal species
     1081     
     1082     LOGICAL ::  match_lsm !< flag indicating natural-type surface
     1083     LOGICAL ::  match_usm !< flag indicating urban-type surface
     1084!
     1085!--  Loop over all sites.
     1086     DO  l = 1, vmea_general%nvm
     1087!
     1088!--     At the beginning, set _FillValues
     1089        IF ( ALLOCATED( vmea(l)%measured_vars      ) )                         &
     1090           vmea(l)%measured_vars      = vmea(l)%fillout
     1091        IF ( ALLOCATED( vmea(l)%measured_vars_soil ) )                         &
     1092           vmea(l)%measured_vars_soil = vmea(l)%fillout
     1093!
     1094!--     Loop over all variables measured at this site. 
    7331095        DO  n = 1, vmea(l)%nvar
    7341096       
     
    7411103                       j = vmea(l)%j(m)
    7421104                       i = vmea(l)%i(m)
    743                        vmea(l)%measured_vars(n,m) = pt(k,j,i)
     1105                       vmea(l)%measured_vars(m,n) = pt(k,j,i)
    7441106                    ENDDO
    7451107                 ENDIF
    7461108                 
    747               CASE ( 'ta', 't_va' )
     1109              CASE ( 'ta' )
    7481110                 IF ( .NOT. neutral )  THEN
    7491111                    DO  m = 1, vmea(l)%ns
     
    7511113                       j = vmea(l)%j(m)
    7521114                       i = vmea(l)%i(m)
    753                        vmea(l)%measured_vars(n,m) = pt(k,j,i) * exner( k )
     1115                       vmea(l)%measured_vars(m,n) = pt(k,j,i) * exner( k )
    7541116                    ENDDO
    7551117                 ENDIF
     1118             
     1119              CASE ( 't_va' )
    7561120                 
    7571121              CASE ( 'hus', 'haa' )
     
    7611125                       j = vmea(l)%j(m)
    7621126                       i = vmea(l)%i(m)
    763                        vmea(l)%measured_vars(n,m) = q(k,j,i)
     1127                       vmea(l)%measured_vars(m,n) = q(k,j,i)
    7641128                    ENDDO
    7651129                 ENDIF
     
    7701134                    j = vmea(l)%j(m)
    7711135                    i = vmea(l)%i(m)
    772                     vmea(l)%measured_vars(n,m) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
     1136                    vmea(l)%measured_vars(m,n) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
    7731137                 ENDDO
    7741138                 
     
    7781142                    j = vmea(l)%j(m)
    7791143                    i = vmea(l)%i(m)
    780                     vmea(l)%measured_vars(n,m) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
     1144                    vmea(l)%measured_vars(m,n) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
    7811145                 ENDDO
    7821146                 
    7831147              CASE ( 'w' )
    7841148                 DO  m = 1, vmea(l)%ns
    785                     k = vmea(l)%k(m)
     1149                    k = MAX ( 1, vmea(l)%k(m) ) 
    7861150                    j = vmea(l)%j(m)
    7871151                    i = vmea(l)%i(m)
    788                     vmea(l)%measured_vars(n,m) = w(k,j,i) !0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
     1152                    vmea(l)%measured_vars(m,n) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
    7891153                 ENDDO
    7901154                 
     
    7941158                    j = vmea(l)%j(m)
    7951159                    i = vmea(l)%i(m)
    796                     vmea(l)%measured_vars(n,m) = SQRT(                         &
     1160                    vmea(l)%measured_vars(m,n) = SQRT(                         &
    7971161                                   ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2 + &
    7981162                                   ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2   &
     
    8061170                    i = vmea(l)%i(m)
    8071171                   
    808                     vmea(l)%measured_vars(n,m) = ATAN2(                        &
     1172                    vmea(l)%measured_vars(m,n) = ATAN2(                        &
    8091173                                       - 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ),   &
    8101174                                       - 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )    &
    8111175                                                      ) * 180.0_wp / pi
    8121176                 ENDDO
    813 !
    814 !--           MS: list of variables needs extension.
    815               CASE ( 'mcpm1', 'mcpm2p5', 'mcpm10', 'mcco', 'mcco2', 'mcbcda',  &
    816                      'ncaa', 'mfco2', 'mfco', 'mfch4', 'mfnh3', 'mfno',        &
    817                      'mfno2', 'mfso2', 'mfh20', 'tr03' )
    818                  IF ( air_chemistry )  THEN                 
     1177                 
     1178              CASE ( 'utheta' )
     1179                 DO  m = 1, vmea(l)%ns
     1180                    k = vmea(l)%k(m)
     1181                    j = vmea(l)%j(m)
     1182                    i = vmea(l)%i(m)
     1183                    vmea(l)%measured_vars(m,n) = 0.5_wp *                      &
     1184                                                 ( u(k,j,i) + u(k,j,i+1) ) *   &
     1185                                                   pt(k,j,i)
     1186                 ENDDO
     1187                 
     1188              CASE ( 'vtheta' )
     1189                 DO  m = 1, vmea(l)%ns
     1190                    k = vmea(l)%k(m)
     1191                    j = vmea(l)%j(m)
     1192                    i = vmea(l)%i(m)
     1193                    vmea(l)%measured_vars(m,n) = 0.5_wp *                      &
     1194                                                 ( v(k,j,i) + v(k,j+1,i) ) *   &
     1195                                                   pt(k,j,i)
     1196                 ENDDO
     1197                 
     1198              CASE ( 'wtheta' )
     1199                 DO  m = 1, vmea(l)%ns
     1200                    k = MAX ( 1, vmea(l)%k(m) )
     1201                    j = vmea(l)%j(m)
     1202                    i = vmea(l)%i(m)
     1203                    vmea(l)%measured_vars(m,n) = 0.5_wp *                      &
     1204                                                 ( w(k-1,j,i) + w(k,j,i) ) *   &
     1205                                                   pt(k,j,i)
     1206                 ENDDO
     1207                 
     1208              CASE ( 'uw' )
     1209                 DO  m = 1, vmea(l)%ns
     1210                    k = MAX ( 1, vmea(l)%k(m) )
     1211                    j = vmea(l)%j(m)
     1212                    i = vmea(l)%i(m)
     1213                    vmea(l)%measured_vars(m,n) = 0.25_wp *                     &
     1214                                                 ( w(k-1,j,i) + w(k,j,i) ) *   &
     1215                                                 ( u(k,j,i)   + u(k,j,i+1) )
     1216                 ENDDO
     1217                 
     1218              CASE ( 'vw' )
     1219                 DO  m = 1, vmea(l)%ns
     1220                    k = MAX ( 1, vmea(l)%k(m) )
     1221                    j = vmea(l)%j(m)
     1222                    i = vmea(l)%i(m)
     1223                    vmea(l)%measured_vars(m,n) = 0.25_wp *                     &
     1224                                                 ( w(k-1,j,i) + w(k,j,i) ) *   &
     1225                                                 ( v(k,j,i)   + v(k,j+1,i) )
     1226                 ENDDO
     1227                 
     1228              CASE ( 'uv' )
     1229                 DO  m = 1, vmea(l)%ns
     1230                    k = MAX ( 1, vmea(l)%k(m) )
     1231                    j = vmea(l)%j(m)
     1232                    i = vmea(l)%i(m)
     1233                    vmea(l)%measured_vars(m,n) = 0.25_wp *                     &
     1234                                                 ( u(k,j,i)   + u(k,j,i+1) ) * &
     1235                                                 ( v(k,j,i)   + v(k,j+1,i) )
     1236                 ENDDO
     1237!
     1238!--           List of variables may need extension.
     1239              CASE ( 'mcpm1', 'mcpm2p5', 'mcpm10', 'mfco', 'mfno', 'mfno2',    &
     1240                     'tro3' )                     
     1241                 IF ( air_chemistry )  THEN
     1242!
     1243!--                 First, search for the measured variable in the chem_vars
     1244!--                 list, in order to get the internal name of the variable.
     1245                    DO  nn = 1, UBOUND( chem_vars, 2 )
     1246                       IF ( TRIM( vmea(l)%measured_vars_name(m) ) ==           &
     1247                            TRIM( chem_vars(0,nn) ) )  ind_chem = nn
     1248                    ENDDO
     1249!
     1250!--                 Run loop over all chemical species, if the measured
     1251!--                 variable matches the interal name, sample the variable.
    8191252                    DO  nn = 1, nspec                   
    820                        IF ( TRIM( vmea(l)%measured_vars_name(m) ) ==           &
     1253                       IF ( TRIM( chem_vars(1,ind_chem) ) ==                   &
    8211254                            TRIM( chem_species(nn)%name ) )  THEN                           
    8221255                          DO  m = 1, vmea(l)%ns             
     
    8241257                             j = vmea(l)%j(m)
    8251258                             i = vmea(l)%i(m)                   
    826                              vmea(l)%measured_vars(n,m) =                      &
     1259                             vmea(l)%measured_vars(m,n) =                      &
    8271260                                                   chem_species(nn)%conc(k,j,i)
    8281261                          ENDDO
     
    8431276                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
    8441277                             surf_def_h(0)%end_index(j,i)
    845                        vmea(l)%measured_vars(n,m) = surf_def_h(0)%us(mm)
     1278                       vmea(l)%measured_vars(m,n) = surf_def_h(0)%us(mm)
    8461279                    ENDDO
    8471280                    DO  mm = surf_lsm_h%start_index(j,i),                      &
    8481281                             surf_lsm_h%end_index(j,i)
    849                        vmea(l)%measured_vars(n,m) = surf_lsm_h%us(mm)
     1282                       vmea(l)%measured_vars(m,n) = surf_lsm_h%us(mm)
    8501283                    ENDDO
    8511284                    DO  mm = surf_usm_h%start_index(j,i),                      &
    8521285                             surf_usm_h%end_index(j,i)
    853                        vmea(l)%measured_vars(n,m) = surf_usm_h%us(mm)
     1286                       vmea(l)%measured_vars(m,n) = surf_usm_h%us(mm)
    8541287                    ENDDO
    8551288                 ENDDO
     
    8671300                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
    8681301                             surf_def_h(0)%end_index(j,i)
    869                        vmea(l)%measured_vars(n,m) = surf_def_h(0)%ts(mm)
     1302                       vmea(l)%measured_vars(m,n) = surf_def_h(0)%ts(mm)
    8701303                    ENDDO
    8711304                    DO  mm = surf_lsm_h%start_index(j,i),                      &
    8721305                             surf_lsm_h%end_index(j,i)
    873                        vmea(l)%measured_vars(n,m) = surf_lsm_h%ts(mm)
     1306                       vmea(l)%measured_vars(m,n) = surf_lsm_h%ts(mm)
    8741307                    ENDDO
    8751308                    DO  mm = surf_usm_h%start_index(j,i),                      &
    8761309                             surf_usm_h%end_index(j,i)
    877                        vmea(l)%measured_vars(n,m) = surf_usm_h%ts(mm)
     1310                       vmea(l)%measured_vars(m,n) = surf_usm_h%ts(mm)
    8781311                    ENDDO
    8791312                 ENDDO
     
    8911324                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
    8921325                             surf_def_h(0)%end_index(j,i)
    893                        vmea(l)%measured_vars(n,m) = surf_def_h(0)%qsws(mm)
     1326                       vmea(l)%measured_vars(m,n) = surf_def_h(0)%qsws(mm)
    8941327                    ENDDO
    8951328                    DO  mm = surf_lsm_h%start_index(j,i),                      &
    8961329                             surf_lsm_h%end_index(j,i)
    897                        vmea(l)%measured_vars(n,m) = surf_lsm_h%qsws(mm)
     1330                       vmea(l)%measured_vars(m,n) = surf_lsm_h%qsws(mm)
    8981331                    ENDDO
    8991332                    DO  mm = surf_usm_h%start_index(j,i),                      &
    9001333                             surf_usm_h%end_index(j,i)
    901                        vmea(l)%measured_vars(n,m) = surf_usm_h%qsws(mm)
     1334                       vmea(l)%measured_vars(m,n) = surf_usm_h%qsws(mm)
    9021335                    ENDDO
    9031336                 ENDDO
     
    9151348                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
    9161349                             surf_def_h(0)%end_index(j,i)
    917                        vmea(l)%measured_vars(n,m) = surf_def_h(0)%shf(mm)
     1350                       vmea(l)%measured_vars(m,n) = surf_def_h(0)%shf(mm)
    9181351                    ENDDO
    9191352                    DO  mm = surf_lsm_h%start_index(j,i),                      &
    9201353                             surf_lsm_h%end_index(j,i)
    921                        vmea(l)%measured_vars(n,m) = surf_lsm_h%shf(mm)
     1354                       vmea(l)%measured_vars(m,n) = surf_lsm_h%shf(mm)
    9221355                    ENDDO
    9231356                    DO  mm = surf_usm_h%start_index(j,i),                      &
    9241357                             surf_usm_h%end_index(j,i)
    925                        vmea(l)%measured_vars(n,m) = surf_usm_h%shf(mm)
     1358                       vmea(l)%measured_vars(m,n) = surf_usm_h%shf(mm)
    9261359                    ENDDO
    9271360                 ENDDO
     
    9401373                       DO  mm = surf_lsm_h%start_index(j,i),                   &
    9411374                                surf_lsm_h%end_index(j,i)
    942                           vmea(l)%measured_vars(n,m) = surf_lsm_h%rad_net(mm)
     1375                          vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_net(mm)
    9431376                       ENDDO
    9441377                       DO  mm = surf_usm_h%start_index(j,i),                   &
    9451378                                surf_usm_h%end_index(j,i)
    946                           vmea(l)%measured_vars(n,m) = surf_usm_h%rad_net(mm)
     1379                          vmea(l)%measured_vars(m,n) = surf_usm_h%rad_net(mm)
    9471380                       ENDDO
    9481381                    ENDDO
    9491382                 ENDIF
    9501383                 
    951               CASE ( 'rsus', 'rsu' )
     1384              CASE ( 'rsus' )
    9521385                 IF ( radiation )  THEN
    9531386                    DO  m = 1, vmea(l)%ns
     
    9621395                       DO  mm = surf_lsm_h%start_index(j,i),                   &
    9631396                                surf_lsm_h%end_index(j,i)
    964                           vmea(l)%measured_vars(n,m) = surf_lsm_h%rad_sw_out(mm)
     1397                          vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_sw_out(mm)
    9651398                       ENDDO
    9661399                       DO  mm = surf_usm_h%start_index(j,i),                   &
    9671400                                surf_usm_h%end_index(j,i)
    968                           vmea(l)%measured_vars(n,m) = surf_usm_h%rad_sw_out(mm)
     1401                          vmea(l)%measured_vars(m,n) = surf_usm_h%rad_sw_out(mm)
    9691402                       ENDDO
    9701403                    ENDDO
    9711404                 ENDIF
    9721405                 
    973               CASE ( 'rsds', 'rsd' )
     1406              CASE ( 'rsds' )
    9741407                 IF ( radiation )  THEN
    9751408                    DO  m = 1, vmea(l)%ns
     
    9841417                       DO  mm = surf_lsm_h%start_index(j,i),                   &
    9851418                                surf_lsm_h%end_index(j,i)
    986                           vmea(l)%measured_vars(n,m) = surf_lsm_h%rad_sw_in(mm)
     1419                          vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_sw_in(mm)
    9871420                       ENDDO
    9881421                       DO  mm = surf_usm_h%start_index(j,i),                   &
    9891422                                surf_usm_h%end_index(j,i)
    990                           vmea(l)%measured_vars(n,m) = surf_usm_h%rad_sw_in(mm)
     1423                          vmea(l)%measured_vars(m,n) = surf_usm_h%rad_sw_in(mm)
    9911424                       ENDDO
    9921425                    ENDDO
    9931426                 ENDIF
    9941427                 
    995               CASE ( 'rlus', 'rlu' )
     1428              CASE ( 'rlus' )
    9961429                 IF ( radiation )  THEN
    9971430                    DO  m = 1, vmea(l)%ns
     
    10061439                       DO  mm = surf_lsm_h%start_index(j,i),                   &
    10071440                                surf_lsm_h%end_index(j,i)
    1008                           vmea(l)%measured_vars(n,m) = surf_lsm_h%rad_lw_out(mm)
     1441                          vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_lw_out(mm)
    10091442                       ENDDO
    10101443                       DO  mm = surf_usm_h%start_index(j,i),                   &
    10111444                                surf_usm_h%end_index(j,i)
    1012                           vmea(l)%measured_vars(n,m) = surf_usm_h%rad_lw_out(mm)
     1445                          vmea(l)%measured_vars(m,n) = surf_usm_h%rad_lw_out(mm)
    10131446                       ENDDO
    10141447                    ENDDO
    10151448                 ENDIF
    10161449                 
    1017               CASE ( 'rlds', 'rld' )
     1450              CASE ( 'rlds' )
    10181451                 IF ( radiation )  THEN
    10191452                    DO  m = 1, vmea(l)%ns
     
    10281461                       DO  mm = surf_lsm_h%start_index(j,i),                   &
    10291462                                surf_lsm_h%end_index(j,i)
    1030                           vmea(l)%measured_vars(n,m) = surf_lsm_h%rad_lw_in(mm)
     1463                          vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_lw_in(mm)
    10311464                       ENDDO
    10321465                       DO  mm = surf_usm_h%start_index(j,i),                   &
    10331466                                surf_usm_h%end_index(j,i)
    1034                           vmea(l)%measured_vars(n,m) = surf_usm_h%rad_lw_in(mm)
     1467                          vmea(l)%measured_vars(m,n) = surf_usm_h%rad_lw_in(mm)
    10351468                       ENDDO
    10361469                    ENDDO
    10371470                 ENDIF
     1471                 
     1472              CASE ( 'rsd' )
     1473                 IF ( radiation )  THEN
     1474                    DO  m = 1, vmea(l)%ns
     1475                       k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' )
     1476                       j = vmea(l)%j(m)
     1477                       i = vmea(l)%i(m)
     1478                   
     1479                       vmea(l)%measured_vars(m,n) = rad_sw_in(k,j,i)
     1480                    ENDDO
     1481                 ENDIF
     1482                 
     1483              CASE ( 'rsu' )
     1484                 IF ( radiation )  THEN
     1485                    DO  m = 1, vmea(l)%ns
     1486                       k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' )
     1487                       j = vmea(l)%j(m)
     1488                       i = vmea(l)%i(m)
     1489                   
     1490                       vmea(l)%measured_vars(m,n) = rad_sw_out(k,j,i)
     1491                    ENDDO
     1492                 ENDIF
     1493                 
     1494              CASE ( 'rlu' )
     1495                 IF ( radiation )  THEN
     1496                    DO  m = 1, vmea(l)%ns
     1497                       k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' )
     1498                       j = vmea(l)%j(m)
     1499                       i = vmea(l)%i(m)
     1500                   
     1501                       vmea(l)%measured_vars(m,n) = rad_lw_out(k,j,i)
     1502                    ENDDO
     1503                 ENDIF
     1504                 
     1505              CASE ( 'rld' )
     1506                 IF ( radiation )  THEN
     1507                    DO  m = 1, vmea(l)%ns
     1508                       k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' )
     1509                       j = vmea(l)%j(m)
     1510                       i = vmea(l)%i(m)
     1511                   
     1512                       vmea(l)%measured_vars(m,n) = rad_lw_in(k,j,i)
     1513                    ENDDO
     1514                 ENDIF
     1515                 
     1516              CASE ( 'rsddif' )
     1517                 IF ( radiation )  THEN
     1518                    DO  m = 1, vmea(l)%ns
     1519                       j = vmea(l)%j(m)
     1520                       i = vmea(l)%i(m)
     1521                   
     1522                       vmea(l)%measured_vars(m,n) = rad_sw_in_diff(j,i)
     1523                    ENDDO
     1524                 ENDIF
     1525                 
     1526              CASE ( 't_soil' )
     1527                 DO  m = 1, vmea(l)%ns_soil
     1528                    i = vmea(l)%i_soil(m)
     1529                    j = vmea(l)%j_soil(m)
     1530                    k = vmea(l)%k_soil(m)
     1531                   
     1532                    match_lsm = surf_lsm_h%start_index(j,i) <=                 &
     1533                                surf_lsm_h%end_index(j,i)
     1534                    match_usm = surf_usm_h%start_index(j,i) <=                 &
     1535                                surf_usm_h%end_index(j,i)
     1536                               
     1537                    IF ( match_lsm )  THEN
     1538                       mm = surf_lsm_h%start_index(j,i)
     1539                       vmea(l)%measured_vars_soil(m,n) = t_soil_h%var_2d(k,mm)
     1540                    ENDIF
     1541                   
     1542                    IF ( match_usm )  THEN
     1543                       mm = surf_usm_h%start_index(j,i)
     1544                       vmea(l)%measured_vars_soil(m,n) = t_wall_h(k,mm)
     1545                    ENDIF
     1546                 ENDDO
     1547                 
     1548              CASE ( 'm_soil' )
     1549                 DO  m = 1, vmea(l)%ns_soil
     1550                    i = vmea(l)%i_soil(m)
     1551                    j = vmea(l)%j_soil(m)
     1552                    k = vmea(l)%k_soil(m)
     1553                   
     1554                    match_lsm = surf_lsm_h%start_index(j,i) <=                 &
     1555                                surf_lsm_h%end_index(j,i)
     1556                               
     1557                    IF ( match_lsm )  THEN
     1558                       mm = surf_lsm_h%start_index(j,i)
     1559                       vmea(l)%measured_vars_soil(m,n) = m_soil_h%var_2d(k,mm)
     1560                    ENDIF
     1561                   
     1562                 ENDDO
    10381563!
    10391564!--           More will follow ...
    1040                  
     1565
     1566!
     1567!--           No match found - just set a fill value
     1568              CASE DEFAULT
     1569                 vmea(l)%measured_vars(:,n) = vmea(l)%fillout
    10411570           END SELECT
    10421571
     
    10441573
    10451574     ENDDO
    1046      
     1575          
    10471576  END SUBROUTINE vm_sampling
    10481577 
Note: See TracChangeset for help on using the changeset viewer.