Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

Changeset 4400


Ignore:
Timestamp:
Feb 10, 2020 8:32:41 PM (5 years ago)
Author:
suehring
Message:

Revision of the virtual-measurement module: data input from NetCDF file; removed binary output - instead parallel NetCDF output using the new data-output module; variable attributes added; further variables added that can be sampled, file connections added; Functions for coordinate transformation moved to basic_constants_and_equations; netcdf_data_input_mod: unused routines netcdf_data_input_att and netcdf_data_input_var removed; new routines to inquire fill values added; Preprocessing script (palm_cvd) to setup virtual-measurement input files provided; postprocessor combine_virtual_measurements deleted

Location:
palm/trunk
Files:
2 added
1 deleted
10 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SCRIPTS/.palm.iofiles

    r3967 r4400  
    6363SURFACE_DATA_AV_NETCDF*    out:tr   *         $base_data/$run_identifier/OUTPUT         _av_surf   nc
    6464#
    65 VIRTUAL_MEAS_BIN*          out:lnpe *         $base_data/$run_identifier/OUTPUT         _vmeas
    66 #
    6765DATA_1D_FL_NETCDF*         out:tr   *         $base_data/$run_identifier/OUTPUT         _fl        nc
    6866DATA_1D_PR_NETCDF*         out:tr   *         $base_data/$run_identifier/OUTPUT         _pr        nc
    … …  
    8179DATA_MASK_AV_NETCDF*       out:tr   *         $base_data/$run_identifier/OUTPUT         _av_masked nc
    8280DATA_AGT_NETCDF*           out:tr   *         $base_data/$run_identifier/OUTPUT         _agt       nc
     81VM_OUTPUT*                 out:pe   *         $base_data/$run_identifier/OUTPUT         _vm
    8382#
    8483DATA_PRT_NETCDF*           out:pe   *         $base_data/$run_identifier/OUTPUT         _prt

  • palm/trunk/SOURCE/Makefile

    r4392 r4400  
    2525# -----------------
    2626# $Id$
     27# Add dependency for data-output module
     28#
     29# 4392 2020-01-31 16:14:57Z pavelkrc
    2730# add dependency on fft for transpose
    2831#
    … …  
    750753        bulk_cloud_model_mod.o \
    751754        chemistry_model_mod.o \
     755        data_output_module.o \
    752756        diagnostic_output_quantities_mod.o \
    753757        dynamics_mod.o \
    … …  
    12731277        user_init_flight.o
    12741278virtual_measurement_mod.o: \
     1279        basic_constants_and_equations_mod.o \
    12751280        cpulog_mod.o \
    12761281        chem_modules.o \
    12771282        chem_gasphase_mod.o \
    1278         mod_kinds.o \
    1279         modules.o \
    1280         netcdf_data_input_mod.o \
    1281         netcdf_interface_mod.o \
    1282         radiation_model_mod.o
     1283        data_output_module.o \
     1284        land_surface_model_mod.o \
     1285        mod_kinds.o \
     1286        modules.o \
     1287        netcdf_data_input_mod.o \
     1288        radiation_model_mod.o \
     1289        surface_mod.o \
     1290        urban_surface_mod.o
    12831291wind_turbine_model_mod.o: \
    12841292        basic_constants_and_equations_mod.o \

  • palm/trunk/SOURCE/basic_constants_and_equations_mod.f90

    r4360 r4400  
    2525! -----------------
    2626! $Id$
     27! Move routine to transform coordinates from netcdf_interface_mod to
     28! basic_constants_and_equations_mod
     29!
     30! 4360 2020-01-07 11:25:50Z suehring
    2731! Corrected "Former revisions" section
    2832!
    … …  
    113117            barometric_formula_1d
    114118
     119
     120    INTERFACE convert_utm_to_geographic
     121       MODULE PROCEDURE convert_utm_to_geographic
     122    END INTERFACE convert_utm_to_geographic
     123
    115124    INTERFACE magnus
    116125       MODULE PROCEDURE magnus_0d
    … …  
    142151       MODULE PROCEDURE barometric_formula_1d
    143152    END INTERFACE barometric_formula
     153!
     154!-- Public routines
     155    PUBLIC convert_utm_to_geographic
    144156
    145157 CONTAINS
     158
     159
     160!------------------------------------------------------------------------------!
     161! Description:
     162! ------------
     163!> Convert UTM coordinates into geographic latitude and longitude. Conversion
     164!> is based on the work of KrÃŒger (1912) DOI: 10.2312/GFZ.b103-krueger28
     165!> and Karney (2013) DOI: 10.1007/s00190-012-0578-z
     166!> Based on a JavaScript of the geodesy function library written by chrisveness
     167!> https://github.com/chrisveness/geodesy
     168!------------------------------------------------------------------------------!
     169    SUBROUTINE convert_utm_to_geographic( crs, eutm, nutm, lon, lat )
     170
     171       INTEGER(iwp) ::  j   !< loop index
     172
     173       REAL(wp), INTENT(in)  ::  eutm !< easting (UTM)
     174       REAL(wp), INTENT(out) ::  lat  !< geographic latitude in degree
     175       REAL(wp), INTENT(out) ::  lon  !< geographic longitude in degree
     176       REAL(wp), INTENT(in)  ::  nutm !< northing (UTM)
     177
     178       REAL(wp) ::  a           !< 2*pi*a is the circumference of a meridian
     179       REAL(wp) ::  cos_eta_s   !< cos(eta_s)
     180       REAL(wp) ::  delta_i     !<
     181       REAL(wp) ::  delta_tau_i !<
     182       REAL(wp) ::  e           !< eccentricity
     183       REAL(wp) ::  eta         !<
     184       REAL(wp) ::  eta_s       !<
     185       REAL(wp) ::  n           !< 3rd flattening
     186       REAL(wp) ::  n2          !< n^2
     187       REAL(wp) ::  n3          !< n^3
     188       REAL(wp) ::  n4          !< n^4
     189       REAL(wp) ::  n5          !< n^5
     190       REAL(wp) ::  n6          !< n^6
     191       REAL(wp) ::  nu          !<
     192       REAL(wp) ::  nu_s        !<
     193       REAL(wp) ::  sin_eta_s   !< sin(eta_s)
     194       REAL(wp) ::  sinh_nu_s   !< sinush(nu_s)
     195       REAL(wp) ::  tau_i       !<
     196       REAL(wp) ::  tau_i_s     !<
     197       REAL(wp) ::  tau_s       !<
     198       REAL(wp) ::  x           !< adjusted easting
     199       REAL(wp) ::  y           !< adjusted northing
     200
     201       REAL(wp), DIMENSION(6) ::  beta !< 6th order KrÃŒger expressions
     202
     203       REAL(wp), DIMENSION(8), INTENT(in) ::  crs !< coordinate reference system, consists of
     204                                                  !< (/semi_major_axis,
     205                                                  !<   inverse_flattening,
     206                                                  !<   longitude_of_prime_meridian,
     207                                                  !<   longitude_of_central_meridian,
     208                                                  !<   scale_factor_at_central_meridian,
     209                                                  !<   latitude_of_projection_origin,
     210                                                  !<   false_easting,
     211                                                  !<   false_northing /)
     212
     213       x = eutm - crs(7)  ! remove false easting
     214       y = nutm - crs(8)  ! remove false northing
     215!
     216!--    from Karney 2011 Eq 15-22, 36:
     217       e = SQRT( 1.0_wp / crs(2) * ( 2.0_wp - 1.0_wp / crs(2) ) )
     218       n = 1.0_wp / crs(2) / ( 2.0_wp - 1.0_wp / crs(2) )
     219       n2 = n * n
     220       n3 = n * n2
     221       n4 = n * n3
     222       n5 = n * n4
     223       n6 = n * n5
     224
     225       a = crs(1) / ( 1.0_wp + n ) * ( 1.0_wp + 0.25_wp * n2       &
     226                                              + 0.015625_wp * n4   &
     227                                              + 3.90625E-3_wp * n6 )
     228
     229       nu  = x / ( crs(5) * a )
     230       eta = y / ( crs(5) * a )
     231
     232!--    According to KrÃŒger (1912), eq. 26*
     233       beta(1) =        0.5_wp                  * n  &
     234               -        2.0_wp /         3.0_wp * n2 &
     235               +       37.0_wp /        96.0_wp * n3 &
     236               -        1.0_wp /       360.0_wp * n4 &
     237               -       81.0_wp /       512.0_wp * n5 &
     238               +    96199.0_wp /    604800.0_wp * n6
     239
     240       beta(2) =        1.0_wp /        48.0_wp * n2 &
     241               +        1.0_wp /        15.0_wp * n3 &
     242               -      437.0_wp /      1440.0_wp * n4 &
     243               +       46.0_wp /       105.0_wp * n5 &
     244               -  1118711.0_wp /   3870720.0_wp * n6
     245
     246       beta(3) =       17.0_wp /       480.0_wp * n3 &
     247               -       37.0_wp /       840.0_wp * n4 &
     248               -      209.0_wp /      4480.0_wp * n5 &
     249               +     5569.0_wp /     90720.0_wp * n6
     250
     251       beta(4) =     4397.0_wp /    161280.0_wp * n4 &
     252               -       11.0_wp /       504.0_wp * n5 &
     253               -   830251.0_wp /   7257600.0_wp * n6
     254
     255       beta(5) =     4583.0_wp /    161280.0_wp * n5 &
     256               -   108847.0_wp /   3991680.0_wp * n6
     257
     258       beta(6) = 20648693.0_wp / 638668800.0_wp * n6
     259
     260       eta_s = eta
     261       nu_s  = nu
     262       DO  j = 1, 6
     263         eta_s = eta_s - beta(j) * SIN(2.0_wp * j * eta) * COSH(2.0_wp * j * nu)
     264         nu_s  = nu_s  - beta(j) * COS(2.0_wp * j * eta) * SINH(2.0_wp * j * nu)
     265       ENDDO
     266
     267       sinh_nu_s = SINH( nu_s )
     268       sin_eta_s = SIN( eta_s )
     269       cos_eta_s = COS( eta_s )
     270
     271       tau_s = sin_eta_s / SQRT( sinh_nu_s**2 + cos_eta_s**2 )
     272
     273       tau_i = tau_s
     274       delta_tau_i = 1.0_wp
     275
     276       DO WHILE ( ABS( delta_tau_i ) > 1.0E-12_wp )
     277
     278         delta_i = SINH( e * ATANH( e * tau_i / SQRT( 1.0_wp + tau_i**2 ) ) )
     279
     280         tau_i_s = tau_i   * SQRT( 1.0_wp + delta_i**2 )  &
     281                  - delta_i * SQRT( 1.0_wp + tau_i**2 )
     282
     283         delta_tau_i = ( tau_s - tau_i_s ) / SQRT( 1.0_wp + tau_i_s**2 )  &
     284                      * ( 1.0_wp + ( 1.0_wp - e**2 ) * tau_i**2 )          &
     285                      / ( ( 1.0_wp - e**2 ) * SQRT( 1.0_wp + tau_i**2 ) )
     286
     287         tau_i = tau_i + delta_tau_i
     288
     289       ENDDO
     290
     291       lat = ATAN( tau_i ) / pi * 180.0_wp
     292       lon = ATAN2( sinh_nu_s, cos_eta_s ) / pi * 180.0_wp + crs(4)
     293
     294    END SUBROUTINE convert_utm_to_geographic
    146295
    147296!------------------------------------------------------------------------------!
    … …  
    161310!--    Saturation vapor pressure for a specific temperature:
    162311       magnus_0d =  611.2_wp * EXP( 17.62_wp * ( t - degc_to_k ) /             &
    163                                                    ( t - 29.65_wp  ) )
     312                                               ( t - 29.65_wp  ) )
    164313
    165314    END FUNCTION magnus_0d

  • palm/trunk/SOURCE/check_open.f90

    r4360 r4400  
    2525! -----------------
    2626! $Id$
     27! Remove binary output for virtual measurements
     28!
     29! 4360 2020-01-07 11:25:50Z suehring
    2730! Corrected "Former revisions" section
    2831!
    … …  
    444447          ENDIF
    445448
    446        CASE ( 27 )
    447 !
    448 !--       Binary files for virtual measurement data
    449           IF ( myid_char == '' )  THEN
    450              OPEN ( 27, FILE='VIRTUAL_MEAS_BIN'//TRIM( coupling_char )//       &
    451                              myid_char, FORM='UNFORMATTED', POSITION='APPEND' )
    452           ELSE
    453 
    454              IF ( myid == 0  .AND. .NOT. openfile(file_id)%opened_before )  THEN
    455                 CALL local_system( 'mkdir  VIRTUAL_MEAS_BIN' //                &
    456                                    TRIM( coupling_char ) )
    457              ENDIF
    458 #if defined( __parallel )
    459 !
    460 !--          Set a barrier in order to allow that all other processors in the
    461 !--          directory created by PE0 can open their file
    462              CALL MPI_BARRIER( comm2d, ierr )
    463 #endif
    464              ioerr = 1
    465              DO WHILE ( ioerr /= 0 )
    466                 OPEN ( 27, FILE='VIRTUAL_MEAS_BIN'//TRIM(coupling_char)//      &
    467                                 '/'//myid_char,                                &
    468                            FORM='UNFORMATTED', IOSTAT=ioerr )
    469                 IF ( ioerr /= 0 )  THEN
    470                    WRITE( 9, * )  '*** could not open "VIRTUAL_MEAS_BIN'//     &
    471                                   TRIM(coupling_char)//'/'//myid_char//        &
    472                                   '"! Trying again in 1 sec.'
    473                    CALL fortran_sleep( 1 )
    474                 ENDIF
    475              ENDDO
    476 
    477           ENDIF
    478 
    479449       CASE ( 30 )
    480450

  • palm/trunk/SOURCE/module_interface.f90

    r4361 r4400  
    2525! -----------------
    2626! $Id$
     27! - Use data-output module for virtual measurement output
     28! - Remove deprecated routines for virtual measurement module
     29!
     30! 4361 2020-01-07 12:22:38Z suehring
    2731! Remove unused arrays in pmc_rrd_local
    2832!
    … …  
    172176    USE kinds
    173177
     178    USE pegrid,                                                                &
     179        ONLY:  comm2d
     180
    174181!
    175182!-- load module-specific control parameters.
    176183!-- ToDo: move all of them to respective module or a dedicated central module
     184    USE data_output_module,                                                    &
     185        ONLY:  dom_def_end,                                                    &
     186               dom_finalize_output,                                            &
     187               dom_init
    177188
    178189    USE dynamics_mod, &
    … …  
    220231        ONLY:  air_chemistry,                                                  &
    221232               biometeorology,                                                 &
     233               coupling_char,                                                  &
    222234               debug_output,                                                   &
    223235               debug_output_timestep,                                          &
    … …  
    493505        ONLY:  vm_check_parameters,                                            &
    494506               vm_init,                                                        &
    495                vm_last_actions,                                                &
     507               vm_init_output,                                                 &
    496508               vm_parin
    497509
    … …  
    543555       module_interface_init,                                                  &
    544556       module_interface_init_checks,                                           &
     557       module_interface_init_output,                                           &
    545558       module_interface_header,                                                &
    546559       module_interface_actions,                                               &
    … …  
    600613       MODULE PROCEDURE module_interface_init_checks
    601614    END INTERFACE module_interface_init_checks
     615
     616    INTERFACE module_interface_init_output
     617       MODULE PROCEDURE module_interface_init_output
     618    END INTERFACE module_interface_init_output
    602619
    603620    INTERFACE module_interface_header
    … …  
    10811098    IF ( debug_output )  CALL debug_message( 'module-specific initialization', 'end' )
    10821099
    1083 
    10841100 END SUBROUTINE module_interface_init
    10851101
     1102!------------------------------------------------------------------------------!
     1103! Description:
     1104! ------------
     1105!> Initialize data output
     1106!------------------------------------------------------------------------------!
     1107 SUBROUTINE module_interface_init_output
     1108
     1109    INTEGER(iwp) ::  return_value  !< returned status value of called function
     1110
     1111!
     1112!-- Initialize data-output module
     1113    CALL dom_init( file_suffix_of_output_group=coupling_char,                  &
     1114                   mpi_comm_of_output_group=comm2d,                            &
     1115                   program_debug_output_unit=6,                                &
     1116                   debug_output=debug_output )
     1117!
     1118!-- Define module-specific output quantities
     1119    IF ( virtual_measurement )  CALL vm_init_output
     1120!
     1121!-- Leave output-definition state
     1122    return_value = dom_def_end()
     1123
     1124 END SUBROUTINE module_interface_init_output
    10861125
    10871126!------------------------------------------------------------------------------!
    … …  
    18901929 SUBROUTINE module_interface_last_actions
    18911930
     1931    INTEGER ::  return_value  !< returned status value of a called function
     1932
    18921933
    18931934    IF ( debug_output )  CALL debug_message( 'module-specific last actions', 'start' )
    18941935
     1936    return_value = dom_finalize_output()
     1937
    18951938    CALL dynamics_last_actions
    18961939
    1897     IF ( virtual_measurement )  CALL vm_last_actions
    1898 
    18991940    IF ( user_module_enabled )  CALL user_last_actions
    19001941

  • palm/trunk/SOURCE/netcdf_data_input_mod.f90

    r4392 r4400  
    2525! -----------------
    2626! $Id$
     27! - Routine to inquire default fill values added
     28! - netcdf_data_input_att and netcdf_data_input_var routines removed
     29!
     30! 4392 2020-01-31 16:14:57Z pavelkrc
    2731! (resler) Decrease length of reading buffer (fix problem of ifort/icc compilers)
    2832!
    … …  
    542546       INTEGER(iwp) ::  version              !< version of data set
    543547
     548       REAL(wp) ::  fillvalue = -9999.0      !< default fill value
    544549       REAL(wp) ::  origin_lat               !< latitude of lower left corner
    545550       REAL(wp) ::  origin_lon               !< longitude of lower left corner
    … …  
    652657    LOGICAL ::  collective_read = .FALSE.      !< Enable NetCDF collective read
    653658
     659    REAL(wp), DIMENSION(8) ::  crs_list        !< list of coord_ref_sys values
     660
    654661    TYPE(global_atts_type) ::  input_file_atts !< global attributes of input file
    655662
    … …  
    666673    END INTERFACE netcdf_data_input_check_static
    667674
    668     INTERFACE netcdf_data_input_chemistry_data                       
     675    INTERFACE netcdf_data_input_chemistry_data
    669676       MODULE PROCEDURE netcdf_data_input_chemistry_data
    670677    END INTERFACE netcdf_data_input_chemistry_data
    671    
    672     INTERFACE get_dimension_length                       
     678
     679    INTERFACE get_dimension_length
    673680       MODULE PROCEDURE get_dimension_length
    674681    END INTERFACE get_dimension_length
     682
     683    INTERFACE inquire_fill_value
     684       MODULE PROCEDURE inquire_fill_value_int
     685       MODULE PROCEDURE inquire_fill_value_real
     686    END INTERFACE inquire_fill_value
    675687
    676688    INTERFACE netcdf_data_input_inquire_file
    … …  
    681693       MODULE PROCEDURE netcdf_data_input_init
    682694    END INTERFACE netcdf_data_input_init
    683    
    684     INTERFACE netcdf_data_input_att
    685        MODULE PROCEDURE netcdf_data_input_att_int8
    686        MODULE PROCEDURE netcdf_data_input_att_int32
    687        MODULE PROCEDURE netcdf_data_input_att_real
    688        MODULE PROCEDURE netcdf_data_input_att_string
    689     END INTERFACE netcdf_data_input_att
    690695
    691696    INTERFACE netcdf_data_input_init_3d
    … …  
    696701       MODULE PROCEDURE netcdf_data_input_surface_data
    697702    END INTERFACE netcdf_data_input_surface_data
    698 
    699     INTERFACE netcdf_data_input_var
    700        MODULE PROCEDURE netcdf_data_input_var_char
    701        MODULE PROCEDURE netcdf_data_input_var_real_1d
    702        MODULE PROCEDURE netcdf_data_input_var_real_2d
    703     END INTERFACE netcdf_data_input_var
    704703
    705704    INTERFACE netcdf_data_input_uvem
    … …  
    752751           chem_emis, chem_emis_att, chem_emis_att_type, chem_emis_val_type,   &
    753752           coord_ref_sys,                                                      &
     753           crs_list,                                                           &
    754754           init_3d, init_model, input_file_atts,                               &
    755755           input_file_dynamic,                                                 &
    … …  
    781781           netcdf_data_input_init,                                             &
    782782           netcdf_data_input_init_3d,                                          &
    783            netcdf_data_input_att,                                              &
    784783           netcdf_data_input_surface_data,                                     &
    785784           netcdf_data_input_topo,                                             &
    786            netcdf_data_input_var,                                              &
    787785           get_attribute,                                                      &
    788786           get_variable,                                                       &
    … …  
    790788           open_read_file,                                                     &
    791789           check_existence,                                                    &
     790           inquire_fill_value,                                                 &
    792791           inquire_num_variables,                                              &
    793792           inquire_variable_names,                                             &
    … …  
    817816       INQUIRE( FILE = TRIM( input_file_chem )    // TRIM( coupling_char ),    &
    818817                EXIST = input_pids_chem )
    819        INQUIRE( FILE = TRIM( input_file_uvem ) // TRIM( coupling_char ),       &
     818       INQUIRE( FILE = TRIM( input_file_uvem )    // TRIM( coupling_char ),    &
    820819                EXIST = input_pids_uvem  )
    821820       INQUIRE( FILE = TRIM( input_file_vm )      // TRIM( coupling_char ),    &
    … …  
    962961       init_model%origin_z        = input_file_atts%origin_z 
    963962       init_model%rotation_angle  = input_file_atts%rotation_angle 
    964            
     963
     964!
     965!--    Define list of crs attributes. This is required for coordinate
     966!--    transformation.
     967       crs_list = (/ coord_ref_sys%semi_major_axis,                            &
     968                     coord_ref_sys%inverse_flattening,                         &
     969                     coord_ref_sys%longitude_of_prime_meridian,                &
     970                     coord_ref_sys%longitude_of_central_meridian,              &
     971                     coord_ref_sys%scale_factor_at_central_meridian,           &
     972                     coord_ref_sys%latitude_of_projection_origin,              &
     973                     coord_ref_sys%false_easting,                              &
     974                     coord_ref_sys%false_northing /)
    965975!
    966976!--    In case of nested runs, each model domain might have different longitude
    … …  
    977987
    978988    END SUBROUTINE netcdf_data_input_init
    979    
    980 !------------------------------------------------------------------------------!
    981 ! Description:
    982 ! ------------
    983 !> Read an array of characters.
    984 !------------------------------------------------------------------------------!
    985     SUBROUTINE netcdf_data_input_var_char( val, search_string, id_mod )
    986 
    987        IMPLICIT NONE
    988 
    989        CHARACTER(LEN=*) ::  search_string     !< name of the variable
    990        CHARACTER(LEN=*), DIMENSION(:) ::  val !< variable which should be read
    991        
    992        INTEGER(iwp) ::  id_mod   !< NetCDF id of input file
    993 
    994 #if defined ( __netcdf )
    995 !
    996 !--    Read variable
    997        CALL get_variable( id_mod, search_string, val )
    998 #endif           
    999 
    1000     END SUBROUTINE netcdf_data_input_var_char
    1001    
    1002 !------------------------------------------------------------------------------!
    1003 ! Description:
    1004 ! ------------
    1005 !> Read an 1D array of REAL values.
    1006 !------------------------------------------------------------------------------!
    1007     SUBROUTINE netcdf_data_input_var_real_1d( val, search_string, id_mod )
    1008 
    1009        IMPLICIT NONE
    1010 
    1011        CHARACTER(LEN=*) ::  search_string     !< name of the variable     
    1012        
    1013        INTEGER(iwp) ::  id_mod        !< NetCDF id of input file
    1014        
    1015        REAL(wp), DIMENSION(:) ::  val !< variable which should be read
    1016 
    1017 #if defined ( __netcdf )
    1018 !
    1019 !--    Read variable
    1020        CALL get_variable( id_mod, search_string, val )
    1021 #endif           
    1022 
    1023     END SUBROUTINE netcdf_data_input_var_real_1d
    1024    
    1025 !------------------------------------------------------------------------------!
    1026 ! Description:
    1027 ! ------------
    1028 !> Read an 1D array of REAL values.
    1029 !------------------------------------------------------------------------------!
    1030     SUBROUTINE netcdf_data_input_var_real_2d( val, search_string,              &
    1031                                               id_mod, d1s, d1e, d2s, d2e )
    1032 
    1033        IMPLICIT NONE
    1034 
    1035        CHARACTER(LEN=*) ::  search_string     !< name of the variable     
    1036        
    1037        INTEGER(iwp) ::  id_mod  !< NetCDF id of input file
    1038        INTEGER(iwp) ::  d1e     !< end index of first dimension to be read
    1039        INTEGER(iwp) ::  d2e     !< end index of second dimension to be read
    1040        INTEGER(iwp) ::  d1s     !< start index of first dimension to be read
    1041        INTEGER(iwp) ::  d2s     !< start index of second dimension to be read
    1042        
    1043        REAL(wp), DIMENSION(:,:) ::  val !< variable which should be read
    1044 
    1045 #if defined ( __netcdf )
    1046 !
    1047 !--    Read character variable
    1048        CALL get_variable( id_mod, search_string, val, d1s, d1e, d2s, d2e )
    1049 #endif           
    1050 
    1051     END SUBROUTINE netcdf_data_input_var_real_2d
    1052    
    1053 !------------------------------------------------------------------------------!
    1054 ! Description:
    1055 ! ------------
    1056 !> Read a global string attribute
    1057 !------------------------------------------------------------------------------!
    1058     SUBROUTINE netcdf_data_input_att_string( val, search_string, id_mod,       &
    1059                                              input_file, global, openclose,    &
    1060                                              variable_name )
    1061 
    1062        IMPLICIT NONE
    1063 
    1064        CHARACTER(LEN=*) ::  search_string !< name of the attribue
    1065        CHARACTER(LEN=*) ::  val           !< attribute
    1066        
    1067        CHARACTER(LEN=*) ::  input_file    !< name of input file
    1068        CHARACTER(LEN=*) ::  openclose     !< string indicating whether NetCDF needs to be opend or closed
    1069        CHARACTER(LEN=*) ::  variable_name !< string indicating whether NetCDF needs to be opend or closed 
    1070        
    1071        INTEGER(iwp) ::  id_mod   !< NetCDF id of input file
    1072        
    1073        LOGICAL ::  global                 !< flag indicating a global or a variable's attribute
    1074 
    1075 #if defined ( __netcdf )
    1076 !
    1077 !--    Open file in read-only mode if necessary
    1078        IF ( openclose == 'open' )  THEN
    1079           CALL open_read_file( TRIM( input_file ) // TRIM( coupling_char ), &
    1080                                   id_mod )
    1081        ENDIF
    1082 !
    1083 !--    Read global attribute
    1084        IF ( global )  THEN
    1085           CALL get_attribute( id_mod, search_string, val, global )
    1086 !
    1087 !--    Read variable attribute
    1088        ELSE
    1089           CALL get_attribute( id_mod, search_string, val, global, variable_name )
    1090        ENDIF
    1091 !
    1092 !--    Close input file
    1093        IF ( openclose == 'close' )  CALL close_input_file( id_mod )
    1094 #endif           
    1095 
    1096     END SUBROUTINE netcdf_data_input_att_string
    1097    
    1098 !------------------------------------------------------------------------------!
    1099 ! Description:
    1100 ! ------------
    1101 !> Read a global 8-bit integer attribute
    1102 !------------------------------------------------------------------------------!
    1103     SUBROUTINE netcdf_data_input_att_int8( val, search_string, id_mod,         &
    1104                                            input_file, global, openclose,      &
    1105                                            variable_name )
    1106 
    1107        IMPLICIT NONE
    1108 
    1109        CHARACTER(LEN=*) ::  search_string !< name of the attribue
    1110        
    1111        CHARACTER(LEN=*) ::  input_file    !< name of input file
    1112        CHARACTER(LEN=*) ::  openclose     !< string indicating whether NetCDF needs to be opend or closed
    1113        CHARACTER(LEN=*) ::  variable_name !< string indicating whether NetCDF needs to be opend or closed
    1114        
    1115        INTEGER(iwp) ::  id_mod   !< NetCDF id of input file
    1116        INTEGER(KIND=1) ::  val      !< value of the attribute
    1117        
    1118        LOGICAL ::  global        !< flag indicating a global or a variable's attribute
    1119 
    1120 #if defined ( __netcdf )
    1121 !
    1122 !--    Open file in read-only mode
    1123        IF ( openclose == 'open' )  THEN
    1124           CALL open_read_file( TRIM( input_file ) // TRIM( coupling_char ), &
    1125                                   id_mod )
    1126        ENDIF
    1127 !
    1128 !--    Read global attribute
    1129        IF ( global )  THEN
    1130           CALL get_attribute( id_mod, search_string, val, global )
    1131 !
    1132 !--    Read variable attribute
    1133        ELSE
    1134           CALL get_attribute( id_mod, search_string, val, global, variable_name )
    1135        ENDIF
    1136 !
    1137 !--    Finally, close input file
    1138        IF ( openclose == 'close' )  CALL close_input_file( id_mod )
    1139 #endif           
    1140 
    1141     END SUBROUTINE netcdf_data_input_att_int8
    1142    
    1143 !------------------------------------------------------------------------------!
    1144 ! Description:
    1145 ! ------------
    1146 !> Read a global 32-bit integer attribute
    1147 !------------------------------------------------------------------------------!
    1148     SUBROUTINE netcdf_data_input_att_int32( val, search_string, id_mod,        &
    1149                                             input_file, global, openclose,     &
    1150                                             variable_name )
    1151 
    1152        IMPLICIT NONE
    1153 
    1154        CHARACTER(LEN=*) ::  search_string !< name of the attribue
    1155        
    1156        CHARACTER(LEN=*) ::  input_file    !< name of input file
    1157        CHARACTER(LEN=*) ::  openclose     !< string indicating whether NetCDF needs to be opend or closed
    1158        CHARACTER(LEN=*) ::  variable_name !< string indicating whether NetCDF needs to be opend or closed
    1159        
    1160        INTEGER(iwp) ::  id_mod   !< NetCDF id of input file
    1161        INTEGER(iwp) ::  val      !< value of the attribute
    1162        
    1163        LOGICAL ::  global        !< flag indicating a global or a variable's attribute
    1164 
    1165 #if defined ( __netcdf )
    1166 !
    1167 !--    Open file in read-only mode
    1168        IF ( openclose == 'open' )  THEN
    1169           CALL open_read_file( TRIM( input_file ) // TRIM( coupling_char ), &
    1170                                   id_mod )
    1171        ENDIF
    1172 !
    1173 !--    Read global attribute
    1174        IF ( global )  THEN
    1175           CALL get_attribute( id_mod, search_string, val, global )
    1176 !
    1177 !--    Read variable attribute
    1178        ELSE
    1179           CALL get_attribute( id_mod, search_string, val, global, variable_name )
    1180        ENDIF
    1181 !
    1182 !--    Finally, close input file
    1183        IF ( openclose == 'close' )  CALL close_input_file( id_mod )
    1184 #endif           
    1185 
    1186     END SUBROUTINE netcdf_data_input_att_int32
    1187    
    1188 !------------------------------------------------------------------------------!
    1189 ! Description:
    1190 ! ------------
    1191 !> Read a global real attribute
    1192 !------------------------------------------------------------------------------!
    1193     SUBROUTINE netcdf_data_input_att_real( val, search_string, id_mod,         &
    1194                                            input_file, global, openclose,      &
    1195                                            variable_name )
    1196 
    1197        IMPLICIT NONE
    1198 
    1199        CHARACTER(LEN=*) ::  search_string !< name of the attribue
    1200        
    1201        CHARACTER(LEN=*) ::  input_file    !< name of input file
    1202        CHARACTER(LEN=*) ::  openclose     !< string indicating whether NetCDF needs to be opend or closed
    1203        CHARACTER(LEN=*) ::  variable_name !< string indicating whether NetCDF needs to be opend or closed
    1204        
    1205        INTEGER(iwp) ::  id_mod            !< NetCDF id of input file
    1206        
    1207        LOGICAL ::  global                 !< flag indicating a global or a variable's attribute
    1208        
    1209        REAL(wp) ::  val                   !< value of the attribute
    1210 
    1211 #if defined ( __netcdf )
    1212 !
    1213 !--    Open file in read-only mode
    1214        IF ( openclose == 'open' )  THEN
    1215           CALL open_read_file( TRIM( input_file ) // TRIM( coupling_char ), &
    1216                                   id_mod )
    1217        ENDIF
    1218 !
    1219 !--    Read global attribute
    1220        IF ( global )  THEN
    1221           CALL get_attribute( id_mod, search_string, val, global )
    1222 !
    1223 !--    Read variable attribute
    1224        ELSE
    1225           CALL get_attribute( id_mod, search_string, val, global, variable_name )
    1226        ENDIF
    1227 !
    1228 !--    Finally, close input file
    1229        IF ( openclose == 'close' )  CALL close_input_file( id_mod )
    1230 #endif           
    1231 
    1232     END SUBROUTINE netcdf_data_input_att_real
     989
    1233990
    1234991!------------------------------------------------------------------------------!
    … …  
    58235580! Description:
    58245581! ------------
     5582!> Inquires the _FillValue settings of an integer variable.
     5583!------------------------------------------------------------------------------!
     5584    SUBROUTINE inquire_fill_value_int( id, var_name, nofill, fill_value )
     5585#if defined( __netcdf )
     5586       CHARACTER(LEN=*), INTENT(IN) ::  var_name    !< variable name
     5587
     5588       INTEGER(iwp), INTENT(IN)  ::  id          !< file id
     5589       INTEGER(iwp)              ::  nofill      !< flag indicating whether fill values are set or not
     5590       INTEGER(iwp)              ::  fill_value  !< fill value
     5591       INTEGER(iwp)              ::  id_var      !< netCDF variable id (varid)
     5592
     5593       nc_stat = NF90_INQ_VARID( id, TRIM( var_name ), id_var )
     5594       nc_stat = NF90_INQ_VAR_FILL(id, id_var, no_fill, fill_value )
     5595#endif
     5596!
     5597!--    Further line is just to avoid compiler warnings. nofill might be used
     5598!--    in future.
     5599       IF ( nofill == 0  .OR.  nofill /= 0 )  CONTINUE
     5600
     5601    END SUBROUTINE inquire_fill_value_int
     5602
     5603!------------------------------------------------------------------------------!
     5604! Description:
     5605! ------------
     5606!> Inquires the _FillValue settings of a real variable.
     5607!------------------------------------------------------------------------------!
     5608    SUBROUTINE inquire_fill_value_real( id, var_name, nofill, fill_value )
     5609#if defined( __netcdf )
     5610       CHARACTER(LEN=*), INTENT(IN) ::  var_name    !< variable name
     5611
     5612       INTEGER(iwp), INTENT(IN)  ::  id          !< file id
     5613       INTEGER(iwp)              ::  nofill      !< flag indicating whether fill values are set or not
     5614       INTEGER(iwp)              ::  id_var      !< netCDF variable id (varid)
     5615
     5616       REAL(wp), INTENT(OUT)     ::  fill_value  !< fill value
     5617
     5618       nc_stat = NF90_INQ_VARID( id, TRIM( var_name ), id_var )
     5619       nc_stat = NF90_INQ_VAR_FILL(id, id_var, no_fill, fill_value )
     5620#endif
     5621!
     5622!--    Further line is just to avoid compiler warnings. nofill might be used
     5623!--    in future.
     5624       IF ( nofill == 0  .OR.  nofill /= 0 )  CONTINUE
     5625
     5626    END SUBROUTINE inquire_fill_value_real
     5627
     5628!------------------------------------------------------------------------------!
     5629! Description:
     5630! ------------
    58255631!> Prints out a text message corresponding to the current status.
    58265632!------------------------------------------------------------------------------!

  • palm/trunk/SOURCE/netcdf_interface_mod.f90

    r4360 r4400  
    2525! -----------------
    2626! $Id$
     27! Move routine to transform coordinates from netcdf_interface_mod to
     28! basic_constants_and_equations_mod
     29!
     30! 4360 2020-01-07 11:25:50Z suehring
    2731! Adjusted output of multi-agent system for biometeorology
    2832!
    … …  
    135139
    136140    USE netcdf_data_input_mod,                                                 &
    137         ONLY: coord_ref_sys, init_model
     141        ONLY: coord_ref_sys,                                                   &
     142              crs_list,                                                        &
     143              init_model
    138144
    139145    PRIVATE
    … …  
    417423
    418424    USE basic_constants_and_equations_mod,                                     &
    419         ONLY:  pi
     425        ONLY:  convert_utm_to_geographic,                                      &
     426               pi
    420427
    421428    USE control_parameters,                                                    &
    … …  
    564571
    565572    REAL(wp), DIMENSION(1) ::  last_time_coordinate          !< last time value in file
    566     REAL(wp), DIMENSION(8) ::  crs_list                      !< list of coord_ref_sys values
    567573
    568574    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  netcdf_data    !<
    … …  
    662668
    663669    ENDIF
    664 !
    665 !-- Convert coord_ref_sys into vector (used for lat/lon calculation)
    666     crs_list = (/ coord_ref_sys%semi_major_axis,                  &
    667                   coord_ref_sys%inverse_flattening,               &
    668                   coord_ref_sys%longitude_of_prime_meridian,      &
    669                   coord_ref_sys%longitude_of_central_meridian,    &
    670                   coord_ref_sys%scale_factor_at_central_meridian, &
    671                   coord_ref_sys%latitude_of_projection_origin,    &
    672                   coord_ref_sys%false_easting,                    &
    673                   coord_ref_sys%false_northing /)
    674670
    675671!
    … …  
    71087104
    71097105
    7110 !------------------------------------------------------------------------------!
    7111 ! Description:
    7112 ! ------------
    7113 !> Convert UTM coordinates into geographic latitude and longitude. Conversion
    7114 !> is based on the work of KrÃŒger (1912) DOI: 10.2312/GFZ.b103-krueger28
    7115 !> and Karney (2013) DOI: 10.1007/s00190-012-0578-z
    7116 !> Based on a JavaScript of the geodesy function library written by chrisveness
    7117 !> https://github.com/chrisveness/geodesy
    7118 !------------------------------------------------------------------------------!
    7119  SUBROUTINE convert_utm_to_geographic( crs, eutm, nutm, lon, lat )
    7120 
    7121     USE basic_constants_and_equations_mod,                                     &
    7122         ONLY:  pi
    7123 
    7124     IMPLICIT NONE
    7125 
    7126     INTEGER(iwp) ::  j   !< loop index
    7127 
    7128     REAL(wp), INTENT(in)  ::  eutm !< easting (UTM)
    7129     REAL(wp), INTENT(out) ::  lat  !< geographic latitude in degree
    7130     REAL(wp), INTENT(out) ::  lon  !< geographic longitude in degree
    7131     REAL(wp), INTENT(in)  ::  nutm !< northing (UTM)
    7132 
    7133     REAL(wp) ::  a           !< 2*pi*a is the circumference of a meridian
    7134     REAL(wp) ::  cos_eta_s   !< cos(eta_s)
    7135     REAL(wp) ::  delta_i     !<
    7136     REAL(wp) ::  delta_tau_i !<
    7137     REAL(wp) ::  e           !< eccentricity
    7138     REAL(wp) ::  eta         !<
    7139     REAL(wp) ::  eta_s       !<
    7140     REAL(wp) ::  n           !< 3rd flattening
    7141     REAL(wp) ::  n2          !< n^2
    7142     REAL(wp) ::  n3          !< n^3
    7143     REAL(wp) ::  n4          !< n^4
    7144     REAL(wp) ::  n5          !< n^5
    7145     REAL(wp) ::  n6          !< n^6
    7146     REAL(wp) ::  nu          !<
    7147     REAL(wp) ::  nu_s        !<
    7148     REAL(wp) ::  sin_eta_s   !< sin(eta_s)
    7149     REAL(wp) ::  sinh_nu_s   !< sinush(nu_s)
    7150     REAL(wp) ::  tau_i       !<
    7151     REAL(wp) ::  tau_i_s     !<
    7152     REAL(wp) ::  tau_s       !<
    7153     REAL(wp) ::  x           !< adjusted easting
    7154     REAL(wp) ::  y           !< adjusted northing
    7155 
    7156     REAL(wp), DIMENSION(6) ::  beta !< 6th order KrÃŒger expressions
    7157 
    7158     REAL(wp), DIMENSION(8), INTENT(in) ::  crs !< coordinate reference system, consists of
    7159                                                !< (/semi_major_axis,
    7160                                                !<   inverse_flattening,
    7161                                                !<   longitude_of_prime_meridian,
    7162                                                !<   longitude_of_central_meridian,
    7163                                                !<   scale_factor_at_central_meridian,
    7164                                                !<   latitude_of_projection_origin,
    7165                                                !<   false_easting,
    7166                                                !<   false_northing /)
    7167 
    7168     x = eutm - crs(7)  ! remove false easting
    7169     y = nutm - crs(8)  ! remove false northing
    7170 
    7171 !-- from Karney 2011 Eq 15-22, 36:
    7172     e = SQRT( 1.0_wp / crs(2) * ( 2.0_wp - 1.0_wp / crs(2) ) )
    7173     n = 1.0_wp / crs(2) / ( 2.0_wp - 1.0_wp / crs(2) )
    7174     n2 = n * n
    7175     n3 = n * n2
    7176     n4 = n * n3
    7177     n5 = n * n4
    7178     n6 = n * n5
    7179 
    7180     a = crs(1) / ( 1.0_wp + n ) * ( 1.0_wp + 0.25_wp * n2       &
    7181                                            + 0.015625_wp * n4   &
    7182                                            + 3.90625E-3_wp * n6 )
    7183 
    7184     nu  = x / ( crs(5) * a )
    7185     eta = y / ( crs(5) * a )
    7186 
    7187 !-- According to KrÃŒger (1912), eq. 26*
    7188     beta(1) =        0.5_wp                  * n  &
    7189             -        2.0_wp /         3.0_wp * n2 &
    7190             +       37.0_wp /        96.0_wp * n3 &
    7191             -        1.0_wp /       360.0_wp * n4 &
    7192             -       81.0_wp /       512.0_wp * n5 &
    7193             +    96199.0_wp /    604800.0_wp * n6
    7194 
    7195     beta(2) =        1.0_wp /        48.0_wp * n2 &
    7196             +        1.0_wp /        15.0_wp * n3 &
    7197             -      437.0_wp /      1440.0_wp * n4 &
    7198             +       46.0_wp /       105.0_wp * n5 &
    7199             -  1118711.0_wp /   3870720.0_wp * n6
    7200 
    7201     beta(3) =       17.0_wp /       480.0_wp * n3 &
    7202             -       37.0_wp /       840.0_wp * n4 &
    7203             -      209.0_wp /      4480.0_wp * n5 &
    7204             +     5569.0_wp /     90720.0_wp * n6
    7205 
    7206     beta(4) =     4397.0_wp /    161280.0_wp * n4 &
    7207             -       11.0_wp /       504.0_wp * n5 &
    7208             -   830251.0_wp /   7257600.0_wp * n6
    7209 
    7210     beta(5) =     4583.0_wp /    161280.0_wp * n5 &
    7211             -   108847.0_wp /   3991680.0_wp * n6
    7212 
    7213     beta(6) = 20648693.0_wp / 638668800.0_wp * n6
    7214 
    7215     eta_s = eta
    7216     nu_s  = nu
    7217     DO  j = 1, 6
    7218       eta_s = eta_s - beta(j) * SIN(2.0_wp * j * eta) * COSH(2.0_wp * j * nu)
    7219       nu_s  = nu_s  - beta(j) * COS(2.0_wp * j * eta) * SINH(2.0_wp * j * nu)
    7220     ENDDO
    7221 
    7222     sinh_nu_s = SINH( nu_s )
    7223     sin_eta_s = SIN( eta_s )
    7224     cos_eta_s = COS( eta_s )
    7225 
    7226     tau_s = sin_eta_s / SQRT( sinh_nu_s**2 + cos_eta_s**2 )
    7227 
    7228     tau_i = tau_s
    7229     delta_tau_i = 1.0_wp
    7230 
    7231     DO WHILE ( ABS( delta_tau_i ) > 1.0E-12_wp )
    7232 
    7233       delta_i = SINH( e * ATANH( e * tau_i / SQRT( 1.0_wp + tau_i**2 ) ) )
    7234 
    7235       tau_i_s = tau_i   * SQRT( 1.0_wp + delta_i**2 )  &
    7236                - delta_i * SQRT( 1.0_wp + tau_i**2 )
    7237 
    7238       delta_tau_i = ( tau_s - tau_i_s ) / SQRT( 1.0_wp + tau_i_s**2 )  &
    7239                    * ( 1.0_wp + ( 1.0_wp - e**2 ) * tau_i**2 )          &
    7240                    / ( ( 1.0_wp - e**2 ) * SQRT( 1.0_wp + tau_i**2 ) )
    7241 
    7242       tau_i = tau_i + delta_tau_i
    7243 
    7244     ENDDO
    7245 
    7246     lat = ATAN( tau_i ) / pi * 180.0_wp
    7247     lon = ATAN2( sinh_nu_s, cos_eta_s ) / pi * 180.0_wp + crs(4)
    7248 
    7249  END SUBROUTINE convert_utm_to_geographic
    7250 
    72517106 END MODULE netcdf_interface

  • palm/trunk/SOURCE/palm.f90

    r4360 r4400  
    2525! -----------------
    2626! $Id$
     27! Add interface to initialize data output with dom
     28!
     29! 4360 2020-01-07 11:25:50Z suehring
    2730! implement new palm_date_time_mod
    2831!
    … …  
    9093
    9194    USE module_interface,                                                      &
    92         ONLY:  module_interface_last_actions
     95        ONLY:  module_interface_init_output,                                   &
     96               module_interface_last_actions
     97
    9398
    9499    USE multi_agent_system_mod,                                                &
    … …  
    269274
    270275    CALL init_3d_model
     276
     277    CALL module_interface_init_output
    271278
    272279!

  • palm/trunk/SOURCE/radiation_model_mod.f90

    r4392 r4400  
    2828! -----------------
    2929! $Id$
     30! Initialize radiation arrays with zero
     31!
     32! 4392 2020-01-31 16:14:57Z pavelkrc
    3033! - Add debug tracing of large radiative fluxes (option trace_fluxes_above)
    3134! - Print exact counts of SVF and CSF if debut_output is enabled
    … …  
    17821785          IF ( .NOT. ALLOCATED ( rad_sw_in ) )  THEN
    17831786             ALLOCATE ( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
     1787             rad_sw_in = 0.0_wp
    17841788          ENDIF
    17851789          IF ( .NOT. ALLOCATED ( rad_sw_out ) )  THEN
    17861790             ALLOCATE ( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
     1791             rad_sw_out = 0.0_wp
    17871792          ENDIF
    17881793
    17891794          IF ( .NOT. ALLOCATED ( rad_lw_in ) )  THEN
    17901795             ALLOCATE ( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
     1796             rad_lw_in = 0.0_wp
    17911797          ENDIF
    17921798          IF ( .NOT. ALLOCATED ( rad_lw_out ) )  THEN
    17931799             ALLOCATE ( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
     1800             rad_lw_out = 0.0_wp
    17941801          ENDIF
    17951802

  • palm/trunk/SOURCE/virtual_measurement_mod.f90

    r4346 r4400  
    2525! -----------------
    2626! $Id$
     27! Revision of the module:
     28! - revised input from NetCDF setup file
     29! - parallel NetCDF output via data-output module ( Tobias Gronemeier )
     30! - variable attributes added
     31! - further variables defined
     32!
     33! 4346 2019-12-18 11:55:56Z motisi
    2734! Introduction of wall_flags_total_0, which currently sets bits based on static
    2835! topography information used in wall_flags_static_0
    29 ! 
     36!
    3037! 4329 2019-12-10 15:46:36Z motisi
    3138! Renamed wall_flags_0 to wall_flags_static_0
    32 ! 
     39!
    3340! 4226 2019-09-10 17:03:24Z suehring
    3441! Netcdf input routine for dimension length renamed
    35 ! 
     42!
    3643! 4182 2019-08-22 15:20:23Z scharf
    3744! Corrected "Former revisions" section
    38 ! 
     45!
    3946! 4168 2019-08-16 13:50:17Z suehring
    4047! Replace function get_topography_top_index by topo_top_ind
    41 ! 
     48!
    4249! 3988 2019-05-22 11:32:37Z kanani
    4350! Add variables to enable steering of output interval for virtual measurements
    44 ! 
     51!
    4552! 3913 2019-04-17 15:12:28Z gronemeier
    4653! Bugfix: rotate positions of measurements before writing them into file
    47 ! 
     54!
    4855! 3910 2019-04-17 11:46:56Z suehring
    4956! Bugfix in rotation of UTM coordinates
    50 ! 
     57!
    5158! 3904 2019-04-16 18:22:51Z gronemeier
    5259! Rotate coordinates of stations by given rotation_angle
    53 ! 
     60!
    5461! 3876 2019-04-08 18:41:49Z knoop
    5562! Remove print statement
    56 ! 
     63!
    5764! 3854 2019-04-02 16:59:33Z suehring
    5865! renamed nvar to nmeas, replaced USE chem_modules by USE chem_gasphase_mod and
    59 ! nspec by nvar 
    60 ! 
     66! nspec by nvar
     67!
    6168! 3766 2019-02-26 16:23:41Z raasch
    6269! unused variables removed
    63 ! 
     70!
    6471! 3718 2019-02-06 11:08:28Z suehring
    6572! Adjust variable name connections between UC2 and chemistry variables
    66 ! 
     73!
    6774! 3717 2019-02-05 17:21:16Z suehring
    6875! Additional check + error numbers adjusted
    69 ! 
     76!
    7077! 3706 2019-01-29 20:02:26Z suehring
    7178! unused variables removed
    72 ! 
     79!
    7380! 3705 2019-01-29 19:56:39Z suehring
    7481! - initialization revised
    7582! - binary data output
    7683! - list of allowed variables extended
    77 ! 
     84!
    7885! 3704 2019-01-29 19:51:41Z suehring
    7986! Sampling of variables
    80 ! 
     87!
    8188! 3473 2018-10-30 20:50:15Z suehring
    8289! Initial revision
    … …  
    8592! --------
    8693! @author Matthias Suehring
     94! @author Tobias Gronemeier
    8795!
    8896! Description:
    8997! ------------
    90 !> The module acts as an interface between 'real-world' observations and 
     98!> The module acts as an interface between 'real-world' observations and
    9199!> model simulations. Virtual measurements will be taken in the model at the
    92 !> coordinates representative for the 'real-world' observation coordinates. 
     100!> coordinates representative for the 'real-world' observation coordinates.
    93101!> More precisely, coordinates and measured quanties will be read from a
    94 !> NetCDF file which contains all required information. In the model, 
     102!> NetCDF file which contains all required information. In the model,
    95103!> the same quantities (as long as all the required components are switched-on)
    96104!> will be sampled at the respective positions and output into an extra file,
    97 !> which allows for straight-forward comparison of model results with 
    98 !> observations. 
     105!> which allows for straight-forward comparison of model results with
     106!> observations.
    99107!>
    100 !> @todo list_of_allowed variables needs careful checking
    101 !> @todo Check if sign of surface fluxes for heat, radiation, etc., follows
    102 !>       the (UC)2 standard
    103 !> @note Fluxes are not processed
     108!> @todo Check why there is an error when _FillValue attributes are added via
     109!>       dom.
     110!> @todo Output of character variable station_name (dom hasn't this feature
     111!>       yet implemented).
    104112!------------------------------------------------------------------------------!
    105113 MODULE virtual_measurement_mod
    106114
    107115    USE arrays_3d,                                                             &
    108         ONLY:  q, pt, u, v, w, zu, zw
     116        ONLY:  dzw,                                                            &
     117               exner,                                                          &
     118               hyp,                                                            &
     119               q,                                                              &
     120               ql,                                                             &
     121               pt,                                                             &
     122               rho_air,                                                        &
     123               u,                                                              &
     124               v,                                                              &
     125               w,                                                              &
     126               zu,                                                             &
     127               zw
    109128
    110129    USE basic_constants_and_equations_mod,                                     &
    111         ONLY:  pi
     130        ONLY:  convert_utm_to_geographic,                                      &
     131               degc_to_k,                                                      &
     132               magnus,                                                         &
     133               pi,                                                             &
     134               rd_d_rv
    112135
    113136    USE chem_gasphase_mod,                                                     &
    … …  
    116139    USE chem_modules,                                                          &
    117140        ONLY:  chem_species
    118        
     141
    119142    USE control_parameters,                                                    &
    120         ONLY:  air_chemistry, dz, humidity, io_blocks, io_group, neutral,      &
    121                message_string, time_since_reference_point, virtual_measurement
     143        ONLY:  air_chemistry,                                                  &
     144               child_domain,                                                   &
     145               coupling_char,                                                  &
     146               dz,                                                             &
     147               end_time,                                                       &
     148               humidity,                                                       &
     149               message_string,                                                 &
     150               neutral,                                                        &
     151               origin_date_time,                                               &
     152               rho_surface,                                                    &
     153               surface_pressure,                                               &
     154               time_since_reference_point,                                     &
     155               virtual_measurement,                                            &
     156               write_binary
    122157
    123158    USE cpulog,                                                                &
    124         ONLY:  cpu_log, log_point
    125        
     159        ONLY:  cpu_log,                                                        &
     160               log_point
     161
     162    USE data_output_module
     163
    126164    USE grid_variables,                                                        &
    127         ONLY:  dx, dy
     165        ONLY:  ddx,                                                            &
     166               ddy,                                                            &
     167               dx,                                                             &
     168               dy
    128169
    129170    USE indices,                                                               &
    130         ONLY:  nzb, nzt, nxl, nxr, nys, nyn, nx, ny, topo_top_ind,             &
     171        ONLY:  nbgp,                                                           &
     172               nzb,                                                            &
     173               nzt,                                                            &
     174               nxl,                                                            &
     175               nxlg,                                                           &
     176               nxr,                                                            &
     177               nxrg,                                                           &
     178               nys,                                                            &
     179               nysg,                                                           &
     180               nyn,                                                            &
     181               nyng,                                                           &
     182               topo_top_ind,                                                   &
    131183               wall_flags_total_0
    132184
    133185    USE kinds
    134    
     186
    135187    USE netcdf_data_input_mod,                                                 &
    136         ONLY:  init_model
    137        
     188        ONLY:  close_input_file,                                               &
     189               coord_ref_sys,                                                  &
     190               crs_list,                                                       &
     191               get_attribute,                                                  &
     192               get_dimension_length,                                           &
     193               get_variable,                                                   &
     194               init_model,                                                     &
     195               input_file_atts,                                                &
     196               input_file_vm,                                                  &
     197               input_pids_static,                                              &
     198               input_pids_vm,                                                  &
     199               inquire_fill_value,                                             &
     200               open_read_file,                                                 &
     201               pids_id
     202
    138203    USE pegrid
    139    
     204
    140205    USE surface_mod,                                                           &
    141         ONLY:  surf_lsm_h, surf_usm_h
    142        
     206        ONLY:  surf_lsm_h,                                                     &
     207               surf_usm_h
     208
    143209    USE land_surface_model_mod,                                                &
    144         ONLY:  nzb_soil, nzs, nzt_soil, zs, t_soil_h, m_soil_h
    145        
    146     USE radiation_model_mod
    147        
     210        ONLY:  m_soil_h,                                                       &
     211               nzb_soil,                                                       &
     212               nzt_soil,                                                       &
     213               t_soil_h,                                                       &
     214               zs
     215
     216    USE radiation_model_mod,                                                   &
     217        ONLY:  rad_lw_in,                                                      &
     218               rad_lw_out,                                                     &
     219               rad_sw_in,                                                      &
     220               rad_sw_in_diff,                                                 &
     221               rad_sw_out,                                                     &
     222               radiation_scheme
     223
    148224    USE urban_surface_mod,                                                     &
    149         ONLY:  nzb_wall, nzt_wall, t_wall_h
     225        ONLY:  nzb_wall,                                                       &
     226               nzt_wall,                                                       &
     227               t_wall_h
    150228
    151229
    152230    IMPLICIT NONE
    153    
     231
    154232    TYPE virt_general
    155        INTEGER(iwp) ::  id_vm     !< NetCDF file id for virtual measurements
    156233       INTEGER(iwp) ::  nvm = 0   !< number of virtual measurements
    157234    END TYPE virt_general
    158235
     236    TYPE virt_var_atts
     237       CHARACTER(LEN=100) ::  coordinates          !< defined longname of the variable
     238       CHARACTER(LEN=100) ::  grid_mapping         !< defined longname of the variable
     239       CHARACTER(LEN=100) ::  long_name            !< defined longname of the variable
     240       CHARACTER(LEN=100) ::  name                 !< variable name
     241       CHARACTER(LEN=100) ::  standard_name        !< defined standard name of the variable
     242       CHARACTER(LEN=100) ::  units                !< unit of the output variable
     243
     244       REAL(wp)           ::  fill_value = -9999.0 !< _FillValue attribute
     245    END TYPE virt_var_atts
     246
    159247    TYPE virt_mea
    160    
    161        CHARACTER(LEN=100)  ::  feature_type      !< type of the measurement
    162        CHARACTER(LEN=100)  ::  filename_original !< name of the original file
    163        CHARACTER(LEN=100)  ::  site              !< name of the measurement site
    164    
    165        CHARACTER(LEN=10), DIMENSION(:), ALLOCATABLE ::  measured_vars_name !< name of the measured variables
    166    
    167        INTEGER(iwp) ::  ns = 0          !< number of observation coordinates on subdomain, for atmospheric measurements
    168        INTEGER(iwp) ::  ns_tot = 0      !< total number of observation coordinates, for atmospheric measurements
    169        INTEGER(iwp) ::  ntraj           !< number of trajectories of a measurement
    170        INTEGER(iwp) ::  nmeas           !< number of measured variables (atmosphere + soil)
    171        
    172        INTEGER(iwp) ::  ns_soil = 0     !< number of observation coordinates on subdomain, for soil measurements
    173        INTEGER(iwp) ::  ns_soil_tot = 0 !< total number of observation coordinates, for soil measurements
    174        
    175        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  dim_t !< number observations individual for each trajectory or station that are no _FillValues
    176        
     248
     249       CHARACTER(LEN=100)  ::  feature_type                     !< type of the real-world measurement
     250       CHARACTER(LEN=100)  ::  feature_type_out = 'timeSeries'  !< type of the virtual measurement
     251                                                                !< (all will be timeSeries, even trajectories)
     252       CHARACTER(LEN=100)  ::  nc_filename                      !< name of the NetCDF output file for the station
     253       CHARACTER(LEN=100)  ::  site                             !< name of the measurement site
     254
     255       CHARACTER(LEN=1000) ::  data_content = REPEAT(' ', 1000) !< string of measured variables (data output only)
     256
     257       INTEGER(iwp) ::  end_coord_a = 0     !< end coordinate in NetCDF file for local atmosphere observations
     258       INTEGER(iwp) ::  end_coord_s = 0     !< end coordinate in NetCDF file for local soil observations
     259       INTEGER(iwp) ::  file_time_index = 0 !< time index in NetCDF output file
     260       INTEGER(iwp) ::  ns = 0              !< number of observation coordinates on subdomain, for atmospheric measurements
     261       INTEGER(iwp) ::  ns_tot = 0          !< total number of observation coordinates, for atmospheric measurements
     262       INTEGER(iwp) ::  n_tr_st             !< number of trajectories / station of a measurement
     263       INTEGER(iwp) ::  nmeas               !< number of measured variables (atmosphere + soil)
     264       INTEGER(iwp) ::  ns_soil = 0         !< number of observation coordinates on subdomain, for soil measurements
     265       INTEGER(iwp) ::  ns_soil_tot = 0     !< total number of observation coordinates, for soil measurements
     266       INTEGER(iwp) ::  start_coord_a = 0   !< start coordinate in NetCDF file for local atmosphere observations
     267       INTEGER(iwp) ::  start_coord_s = 0   !< start coordinate in NetCDF file for local soil observations
     268
     269       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  dim_t !< number observations individual for each trajectory
     270                                                         !< or station that are no _FillValues
     271
    177272       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i       !< grid index for measurement position in x-direction
    178273       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  j       !< grid index for measurement position in y-direction
    179274       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k       !< grid index for measurement position in k-direction
    180        
     275
    181276       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i_soil  !< grid index for measurement position in x-direction
    182277       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  j_soil  !< grid index for measurement position in y-direction
    183278       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k_soil  !< grid index for measurement position in k-direction
    184            
     279
    185280       LOGICAL ::  trajectory         = .FALSE. !< flag indicating that the observation is a mobile observation
    186281       LOGICAL ::  timseries          = .FALSE. !< flag indicating that the observation is a stationary point measurement
    187282       LOGICAL ::  timseries_profile  = .FALSE. !< flag indicating that the observation is a stationary profile measurement
    188283       LOGICAL ::  soil_sampling      = .FALSE. !< flag indicating that soil state variables were sampled
    189        
    190        REAL(wp) ::  fill_eutm          !< fill value for UTM coordinates in case of missing values
    191        REAL(wp) ::  fill_nutm          !< fill value for UTM coordinates in case of missing values
    192        REAL(wp) ::  fill_zag           !< fill value for heigth coordinates in case of missing values
    193        REAL(wp) ::  fillout = -999.9   !< fill value for output in case a observation is taken from inside a building
    194        REAL(wp) ::  origin_x_obs       !< origin of the observation in UTM coordiates in x-direction
    195        REAL(wp) ::  origin_y_obs       !< origin of the observation in UTM coordiates in y-direction
    196        
    197        REAL(wp), DIMENSION(:), ALLOCATABLE   ::  z_ag           !< measurement height above ground level
    198        REAL(wp), DIMENSION(:), ALLOCATABLE   ::  depth          !< measurement depth in soil
    199              
    200        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  measured_vars       !< measured variables
    201        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  measured_vars_soil  !< measured variables
    202        
     284
     285       REAL(wp) ::  fill_eutm                            !< fill value for UTM coordinates in case of missing values
     286       REAL(wp) ::  fill_nutm                            !< fill value for UTM coordinates in case of missing values
     287       REAL(wp) ::  fill_zar                             !< fill value for heigth coordinates in case of missing values
     288       REAL(wp) ::  fillout = -9999.0                    !< fill value for output in case a observation is taken
     289                                                         !< e.g. from inside a building
     290       REAL(wp) ::  origin_x_obs                         !< origin of the observation in UTM coordiates in x-direction
     291       REAL(wp) ::  origin_y_obs                         !< origin of the observation in UTM coordiates in y-direction
     292
     293       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  zar           !< measurement height above ground level
     294       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  depth         !< measurement depth in soil
     295
     296       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  measured_vars       !< measured variables
     297       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  measured_vars_soil  !< measured variables
     298
     299       TYPE( virt_var_atts ), DIMENSION(:), ALLOCATABLE ::  var_atts !< variable attributes
     300
    203301    END TYPE virt_mea
    204302
    205303    CHARACTER(LEN=5)  ::  char_eutm = "E_UTM"                      !< dimension name for UTM coordinate easting
    206304    CHARACTER(LEN=11) ::  char_feature = "featureType"             !< attribute name for feature type
    207    
    208     ! This need to generalized
    209     CHARACTER(LEN=8)  ::  char_filename = "filename"               !< attribute name for filename
     305
     306    ! This need to be generalized
     307!     CHARACTER(LEN=10) ::  char_fill = '_FillValue'
     308    CHARACTER(LEN=9)  ::  char_long = 'long_name'                  !< attribute name for long_name
     309    CHARACTER(LEN=13) ::  char_standard = 'standard_name'          !< attribute name for standard_name
     310    CHARACTER(LEN=5)  ::  char_unit = 'units'                      !< attribute name for standard_name
    210311    CHARACTER(LEN=11) ::  char_soil = "soil_sample"                !< attribute name for soil sampling indication
    211     CHARACTER(LEN=10) ::  char_fillvalue = "_FillValue"            !< variable attribute name for _FillValue
    212312    CHARACTER(LEN=18) ::  char_mv = "measured_variables"           !< variable name for the array with the measured variable names
    213313    CHARACTER(LEN=5)  ::  char_nutm = "N_UTM"                      !< dimension name for UTM coordinate northing
    … …  
    216316    CHARACTER(LEN=8)  ::  char_origy = "origin_y"                  !< attribute name for station coordinate in y
    217317    CHARACTER(LEN=4)  ::  char_site = "site"                       !< attribute name for site name
    218     CHARACTER(LEN=19) ::  char_zag = "height_above_ground"         !< attribute name for height above ground variable
     318    CHARACTER(LEN=9)  ::  char_station_h = "station_h"             !< variable name indicating height of the site
     319    CHARACTER(LEN=1)  ::  char_zar = "z"                           !< attribute name indicating height above reference level
    219320    CHARACTER(LEN=10) ::  type_ts   = 'timeSeries'                 !< name of stationary point measurements
    220321    CHARACTER(LEN=10) ::  type_traj = 'trajectory'                 !< name of line measurements
    221322    CHARACTER(LEN=17) ::  type_tspr = 'timeSeriesProfile'          !< name of stationary profile measurements
    222    
     323
    223324    CHARACTER(LEN=6), DIMENSION(1:5) ::  soil_vars       = (/                  & !< list of soil variables
    224325                            't_soil',                                          &
    … …  
    227328                            'lwcs  ',                                          &
    228329                            'smp   '                       /)
    229                            
     330
    230331    CHARACTER(LEN=10), DIMENSION(0:1,1:8) ::  chem_vars = RESHAPE( (/          &
    231332                                              'mcpm1     ', 'PM1       ',      &
    232333                                              'mcpm2p5   ', 'PM2.5     ',      &
    233                                               'mcpm25    ', 'PM25      ',      &
    234334                                              'mcpm10    ', 'PM10      ',      &
    235335                                              'mfno2     ', 'NO2       ',      &
    236336                                              'mfno      ', 'NO        ',      &
    237                                               'tro3      ', 'O3        ',      &
    238                                               'mfco      ', 'CO        '       &
     337                                              'mcno2     ', 'NO2       ',      &
     338                                              'mcno      ', 'NO        ',      &
     339                                              'tro3      ', 'O3        '       &
    239340                                                                   /), (/ 2, 8 /) )
    240 !
    241 !-- MS: List requires careful revision!
    242     CHARACTER(LEN=10), DIMENSION(1:54), PARAMETER ::  list_allowed_variables = & !< variables that can be sampled in PALM
    243        (/ 'hfls      ',  & ! surface latent heat flux (W/m2)
    244           'hfss      ',  & ! surface sensible heat flux (W/m2)
    245           'hur       ',  & ! relative humidity (-)
    246           'hus       ',  & ! specific humidity (g/kg)
    247           'haa       ',  & ! absolute atmospheric humidity (kg/m3)
    248           'mcpm1     ',  & ! mass concentration of PM1 (kg/m3)
    249           'mcpm2p5   ',  & ! mass concentration of PM2.5 (kg/m3)
    250           'mcpm10    ',  & ! mass concentration of PM10 (kg/m3)
    251           'mcco      ',  & ! mass concentration of CO (kg/m3)
    252           'mcco2     ',  & ! mass concentration of CO2 (kg/m3)
    253           'mcbcda    ',  & ! mass concentration of black carbon paritcles (kg/m3)
    254           'ncaa      ',  & ! number concentation of particles (1/m3)
    255           'mfco      ',  & ! mole fraction of CO (mol/mol)
    256           'mfco2     ',  & ! mole fraction of CO2 (mol/mol)
    257           'mfch4     ',  & ! mole fraction of methane (mol/mol)
    258           'mfnh3     ',  & ! mole fraction of amonia (mol/mol)
    259           'mfno      ',  & ! mole fraction of nitrogen monoxide (mol/mol)
    260           'mfno2     ',  & ! mole fraction of nitrogen dioxide (mol/mol)
    261           'mfso2     ',  & ! mole fraction of sulfur dioxide (mol/mol)
    262           'mfh20     ',  & ! mole fraction of water (mol/mol)
    263           'plev      ',  & ! ? air pressure - hydrostaic + perturbation?
    264           'rlds      ',  & ! surface downward longwave flux  (W/m2)
    265           'rlus      ',  & ! surface upward longwave flux (W/m2)
    266           'rsds      ',  & ! surface downward shortwave flux (W/m2)
    267           'rsus      ',  & ! surface upward shortwave flux (W/m2)
    268           'ta        ',  & ! air temperature (degree C)
    269           't_va      ',  & ! virtual accoustic temperature (K)
    270           'theta     ',  & ! potential temperature (K)
    271           'tro3      ',  & ! mole fraction of ozone air (mol/mol)
    272           'ts        ',  & ! scaling parameter of temperature (K)
    273           'wspeed    ',  & ! ? wind speed - horizontal?
    274           'wdir      ',  & ! wind direction
    275           'us        ',  & ! friction velocity
    276           'msoil     ',  & ! ? soil moisture - which depth? 
    277           'tsoil     ',  & ! ? soil temperature - which depth?                                                               
    278           'u         ',  & ! u-component
    279           'utheta    ',  & ! total eastward kinematic heat flux
    280           'ua        ',  & ! eastward wind (is there any difference to u?)
    281           'v         ',  & ! v-component
    282           'vtheta    ',  & ! total northward kinematic heat flux
    283           'va        ',  & ! northward wind (is there any difference to v?)
    284           'w         ',  & ! w-component
    285           'wtheta    ',  & ! total vertical kinematic heat flux
    286           'rld       ',  & ! downward longwave radiative flux (W/m2)
    287           'rlu       ',  & ! upnward longwave radiative flux (W/m2)
    288           'rsd       ',  & ! downward shortwave radiative flux (W/m2)
    289           'rsu       ',  & ! upward shortwave radiative flux (W/m2)
    290           'rsddif    ',  & ! downward shortwave diffuse radiative flux (W/m2)
    291           'rnds      ',  & ! surface net downward radiative flux (W/m2)
    292           't_soil    ',  &
    293           'm_soil    ',  &
    294           'lwc       ',  &
    295           'lwcs      ',  &
    296           'smp       '   &
    297        /)
    298                                                            
    299    
    300     LOGICAL ::  global_attribute = .TRUE.         !< flag indicating a global attribute
    301     LOGICAL ::  init = .TRUE.                     !< flag indicating initialization of data output
    302     LOGICAL ::  use_virtual_measurement = .FALSE. !< Namelist parameter
    303 
    304     REAL(wp) ::  dt_virtual_measurement = 0.0_wp    !< namelist parameter
    305     REAL(wp) ::  time_virtual_measurement = 0.0_wp  !< time since last 3d output
    306     REAL(wp) ::  vm_time_start = 0.0                !< time after virtual measurements should start
    307 
    308     TYPE( virt_general )                        ::  vmea_general !< data structure which encompass general variables
    309     TYPE( virt_mea ), DIMENSION(:), ALLOCATABLE ::  vmea !< virtual measurement data structure
    310    
     341
     342    LOGICAL ::  global_attribute = .TRUE.           !< flag indicating a global attribute
     343    LOGICAL ::  initial_write_coordinates = .FALSE. !< flag indicating a global attribute
     344    LOGICAL ::  use_virtual_measurement = .FALSE.   !< Namelist parameter
     345
     346    INTEGER(iwp) ::  maximum_name_length = 32 !< maximum name length of station names
     347    INTEGER(iwp) ::  ntimesteps               !< number of timesteps defined in NetCDF output file
     348    INTEGER(iwp) ::  ntimesteps_max = 80000   !< number of maximum timesteps defined in NetCDF output file
     349    INTEGER(iwp) ::  off_pr = 1               !< number neighboring grid points (in each direction) where virtual profile
     350                                              !< measurements shall be taken, in addition to the given coordinates in the driver
     351    INTEGER(iwp) ::  off_ts = 1               !< number neighboring grid points (in each direction) where virtual timeseries
     352                                              !< measurements shall be taken, in addition to the given coordinates in the driver
     353    INTEGER(iwp) ::  off_tr = 1               !< number neighboring grid points (in each direction) where virtual trajectory
     354                                              !< measurements shall be taken, in addition to the given coordinates in the driver
     355
     356    REAL(wp) ::  dt_virtual_measurement = 0.0_wp    !< sampling interval
     357    REAL(wp) ::  time_virtual_measurement = 0.0_wp  !< time since last sampling
     358    REAL(wp) ::  vm_time_start = 0.0                !< time after which sampling shall start
     359
     360    TYPE( virt_general )                        ::  vmea_general !< data structure which encompass global variables
     361    TYPE( virt_mea ), DIMENSION(:), ALLOCATABLE ::  vmea         !< data structure contain station-specific variables
     362
    311363    INTERFACE vm_check_parameters
    312364       MODULE PROCEDURE vm_check_parameters
    313365    END INTERFACE vm_check_parameters
    314    
     366
    315367    INTERFACE vm_data_output
    316368       MODULE PROCEDURE vm_data_output
    317369    END INTERFACE vm_data_output
    318    
     370
    319371    INTERFACE vm_init
    320372       MODULE PROCEDURE vm_init
    321373    END INTERFACE vm_init
    322    
    323     INTERFACE vm_last_actions
    324        MODULE PROCEDURE vm_last_actions
    325     END INTERFACE vm_last_actions
    326    
     374
     375    INTERFACE vm_init_output
     376       MODULE PROCEDURE vm_init_output
     377    END INTERFACE vm_init_output
     378
    327379    INTERFACE vm_parin
    328380       MODULE PROCEDURE vm_parin
    329381    END INTERFACE vm_parin
    330    
     382
    331383    INTERFACE vm_sampling
    332384       MODULE PROCEDURE vm_sampling
    … …  
    339391!
    340392!-- Public interfaces
    341     PUBLIC  vm_check_parameters, vm_data_output, vm_init, vm_last_actions,     &
    342             vm_parin, vm_sampling
     393    PUBLIC  vm_check_parameters,                                               &
     394            vm_data_output,                                                    &
     395            vm_init,                                                           &
     396            vm_init_output,                                                    &
     397            vm_parin,                                                          &
     398            vm_sampling
    343399
    344400!
    345401!-- Public variables
    346     PUBLIC  dt_virtual_measurement, time_virtual_measurement, vmea, vmea_general, vm_time_start
     402    PUBLIC  dt_virtual_measurement,                                            &
     403            time_virtual_measurement,                                          &
     404            vmea,                                                              &
     405            vmea_general,                                                      &
     406            vm_time_start
    347407
    348408 CONTAINS
    … …  
    356416 SUBROUTINE vm_check_parameters
    357417
    358     USE control_parameters,                                                    &
    359         ONLY:  message_string, virtual_measurement
    360  
    361     USE netcdf_data_input_mod,                                                 &
    362         ONLY:  input_pids_static, input_pids_vm
    363        
    364     IMPLICIT NONE
    365 
    366 !
    367 !-- Virtual measurements require a setup file.
    368     IF ( virtual_measurement  .AND.  .NOT. input_pids_vm )  THEN
     418    IF ( .NOT. virtual_measurement )  RETURN
     419!
     420!-- Virtual measurements require a setup file.
     421    IF ( .NOT. input_pids_vm )  THEN
    369422       message_string = 'If virtual measurements are taken, a setup input ' // &
    370423                        'file for the site locations is mandatory.'
    371424       CALL message( 'vm_check_parameters', 'PA0533', 1, 2, 0, 6, 0 )
    372     ENDIF   
     425    ENDIF
    373426!
    374427!-- In case virtual measurements are taken, a static input file is required.
    375428!-- This is because UTM coordinates for the PALM domain origin are required
    376 !-- for correct mapping of the measurements. 
     429!-- for correct mapping of the measurements.
    377430!-- ToDo: Revise this later and remove this requirement.
    378     IF ( virtual_measurement  .AND.  .NOT. input_pids_static )  THEN
     431    IF ( .NOT. input_pids_static )  THEN
    379432       message_string = 'If virtual measurements are taken, a static input ' //&
    380433                        'file is mandatory.'
    381434       CALL message( 'vm_check_parameters', 'PA0534', 1, 2, 0, 6, 0 )
    382435    ENDIF
    383  
     436
     437#if !defined( __netcdf4_parallel )
     438!
     439!-- In case of non-parallel NetCDF the virtual measurement output is not
     440!-- working. This is only designed for parallel NetCDF.
     441    message_string = 'If virtual measurements are taken, parallel ' //         &
     442                     'NetCDF is required.'
     443    CALL message( 'vm_check_parameters', 'PA0708', 1, 2, 0, 6, 0 )
     444#endif
     445!
     446!-- Check if the given number of neighboring grid points do not exceeds the number
     447!-- of ghost points.
     448    IF ( off_pr > nbgp - 1  .OR.  off_ts > nbgp - 1  .OR.  off_tr > nbgp - 1 ) &
     449    THEN
     450       WRITE(message_string,*)                                                 &
     451                        'If virtual measurements are taken, the number ' //    &
     452                        'of surrounding grid points must not be larger ' //    &
     453                        'than the number of ghost points - 1, which is: ',     &
     454                        nbgp - 1
     455       CALL message( 'vm_check_parameters', 'PA0705', 1, 2, 0, 6, 0 )
     456    ENDIF
     457!
     458!-- Check if restart activation string is set. Virtual measurements may exceed
     459!-- the maximum number of allowed timesteps in a NetCDF file. In order to handle
     460!-- this, the module will trigger a programm abortion with writing the restart
     461!-- data.
     462    IF ( .NOT. write_binary )  THEN
     463       message_string = 'virtual measurements require file activation ' //     &
     464                        'string "restart"'
     465       CALL message( 'check_parameters', 'PA0706', 1, 2, 0, 6, 0 )
     466    ENDIF
     467
    384468 END SUBROUTINE vm_check_parameters
    385469 
     470!------------------------------------------------------------------------------!
     471! Description:
     472! ------------
     473!> Subroutine defines variable attributes according to UC2 standard. Note, later
     474!> this list can be moved to the data-output module where it can be re-used also
     475!> for other output.
     476!------------------------------------------------------------------------------!
     477  SUBROUTINE vm_set_attributes( output_variable )
     478
     479     TYPE( virt_var_atts ), INTENT(INOUT) ::  output_variable !< data structure with attributes that need to be set
     480
     481     output_variable%long_name     = 'none'
     482     output_variable%standard_name = 'none'
     483     output_variable%units         = 'none'
     484     output_variable%coordinates   = 'lon lat E_UTM N_UTM x y z time station_name'
     485     output_variable%grid_mapping  = 'crs'
     486
     487     SELECT CASE ( TRIM( output_variable%name ) )
     488
     489        CASE ( 'u' )
     490           output_variable%long_name     = 'u wind component'
     491           output_variable%units         = 'm s-1'
     492
     493        CASE ( 'ua' )
     494           output_variable%long_name     = 'eastward wind'
     495           output_variable%standard_name = 'eastward_wind'
     496           output_variable%units         = 'm s-1'
     497
     498        CASE ( 'v' )
     499           output_variable%long_name     = 'v wind component'
     500           output_variable%units         = 'm s-1'
     501
     502        CASE ( 'va' )
     503           output_variable%long_name     = 'northward wind'
     504           output_variable%standard_name = 'northward_wind'
     505           output_variable%units         = 'm s-1'
     506
     507        CASE ( 'w' )
     508           output_variable%long_name     = 'w wind component'
     509           output_variable%standard_name = 'upward_air_velocity'
     510           output_variable%units         = 'm s-1'
     511
     512        CASE ( 'wspeed' )
     513           output_variable%long_name     = 'wind speed'
     514           output_variable%standard_name = 'wind_speed'
     515           output_variable%units         = 'm s-1'
     516
     517        CASE ( 'wdir' )
     518           output_variable%long_name     = 'wind from direction'
     519           output_variable%standard_name = 'wind_from_direction'
     520           output_variable%units         = 'degrees'
     521
     522        CASE ( 'theta' )
     523           output_variable%long_name     = 'air potential temperature'
     524           output_variable%standard_name = 'air_potential_temperature'
     525           output_variable%units         = 'K'
     526
     527        CASE ( 'utheta' )
     528           output_variable%long_name     = 'eastward kinematic sensible heat flux in air'
     529           output_variable%units         = 'K m s-1'
     530
     531        CASE ( 'vtheta' )
     532           output_variable%long_name     = 'northward kinematic sensible heat flux in air'
     533           output_variable%units         = 'K m s-1'
     534
     535        CASE ( 'wtheta' )
     536           output_variable%long_name     = 'upward kinematic sensible heat flux in air'
     537           output_variable%units         = 'K m s-1'
     538
     539        CASE ( 'ta' )
     540           output_variable%long_name     = 'air temperature'
     541           output_variable%standard_name = 'air_temperature'
     542           output_variable%units         = 'degree_C'
     543
     544        CASE ( 'tva' )
     545           output_variable%long_name     = 'virtual acoustic temperature'
     546           output_variable%units         = 'K'
     547
     548        CASE ( 'haa' )
     549           output_variable%long_name     = 'absolute atmospheric humidity'
     550           output_variable%units         = 'kg m-3'
     551
     552        CASE ( 'hus' )
     553           output_variable%long_name     = 'specific humidity'
     554           output_variable%standard_name = 'specific_humidity'
     555           output_variable%units         = 'kg kg-1'
     556
     557        CASE ( 'hur' )
     558           output_variable%long_name     = 'relative humidity'
     559           output_variable%standard_name = 'relative_humidity'
     560           output_variable%units         = '1'
     561
     562        CASE ( 'rlu' )
     563           output_variable%long_name     = 'upwelling longwave flux in air'
     564           output_variable%standard_name = 'upwelling_longwave_flux_in_air'
     565           output_variable%units         = 'W m-2'
     566
     567        CASE ( 'rlus' )
     568           output_variable%long_name     = 'surface upwelling longwave flux in air'
     569           output_variable%standard_name = 'surface_upwelling_longwave_flux_in_air'
     570           output_variable%units         = 'W m-2'
     571
     572        CASE ( 'rld' )
     573           output_variable%long_name     = 'downwelling longwave flux in air'
     574           output_variable%standard_name = 'downwelling_longwave_flux_in_air'
     575           output_variable%units         = 'W m-2'
     576
     577        CASE ( 'rsddif' )
     578           output_variable%long_name     = 'diffuse downwelling shortwave flux in air'
     579           output_variable%standard_name = 'diffuse_downwelling_shortwave_flux_in_air'
     580           output_variable%units         = 'W m-2'
     581
     582        CASE ( 'rsd' )
     583           output_variable%long_name     = 'downwelling shortwave flux in air'
     584           output_variable%standard_name = 'downwelling_shortwave_flux_in_air'
     585           output_variable%units         = 'W m-2'
     586
     587        CASE ( 'rnds' )
     588           output_variable%long_name     = 'surface net downward radiative flux'
     589           output_variable%standard_name = 'surface_net_downward_radiative_flux'
     590           output_variable%units         = 'W m-2'
     591
     592        CASE ( 'rsu' )
     593           output_variable%long_name     = 'upwelling shortwave flux in air'
     594           output_variable%standard_name = 'upwelling_shortwave_flux_in_air'
     595           output_variable%units         = 'W m-2'
     596
     597        CASE ( 'rsus' )
     598           output_variable%long_name     = 'surface upwelling shortwave flux in air'
     599           output_variable%standard_name = 'surface_upwelling_shortwave_flux_in_air'
     600           output_variable%units         = 'W m-2'
     601
     602        CASE ( 'rsds' )
     603           output_variable%long_name     = 'surface downwelling shortwave flux in air'
     604           output_variable%standard_name = 'surface_downwelling_shortwave_flux_in_air'
     605           output_variable%units         = 'W m-2'
     606
     607        CASE ( 'hfss' )
     608           output_variable%long_name     = 'surface upward sensible heat flux'
     609           output_variable%standard_name = 'surface_upward_sensible_heat_flux'
     610           output_variable%units         = 'W m-2'
     611
     612        CASE ( 'hfls' )
     613           output_variable%long_name     = 'surface upward latent heat flux'
     614           output_variable%standard_name = 'surface_upward_latent_heat_flux'
     615           output_variable%units         = 'W m-2'
     616
     617        CASE ( 'ts' )
     618           output_variable%long_name     = 'surface temperature'
     619           output_variable%standard_name = 'surface_temperature'
     620           output_variable%units         = 'K'
     621
     622        CASE ( 'thetas' )
     623           output_variable%long_name     = 'surface layer temperature scale'
     624           output_variable%units         = 'K'
     625
     626        CASE ( 'us' )
     627           output_variable%long_name     = 'friction velocity'
     628           output_variable%units         = 'm s-1'
     629
     630        CASE ( 'uw' )
     631           output_variable%long_name     = 'upward eastward kinematic momentum flux in air'
     632           output_variable%units         = 'm2 s-2'
     633
     634        CASE ( 'vw' )
     635           output_variable%long_name     = 'upward northward kinematic momentum flux in air'
     636           output_variable%units         = 'm2 s-2'
     637
     638        CASE ( 'uv' )
     639           output_variable%long_name     = 'eastward northward kinematic momentum flux in air'
     640           output_variable%units         = 'm2 s-2'
     641
     642        CASE ( 'plev' )
     643           output_variable%long_name     = 'air pressure'
     644           output_variable%standard_name = 'air_pressure'
     645           output_variable%units         = 'Pa'
     646
     647        CASE ( 'm_soil' )
     648           output_variable%long_name     = 'soil moisture volumetric'
     649           output_variable%units         = 'm3 m-3'
     650
     651        CASE ( 't_soil' )
     652           output_variable%long_name     = 'soil temperature'
     653           output_variable%standard_name = 'soil_temperature'
     654           output_variable%units         = 'degree_C'
     655
     656        CASE ( 'hfdg' )
     657           output_variable%long_name     = 'downward heat flux at ground level in soil'
     658           output_variable%standard_name = 'downward_heat_flux_at_ground_level_in_soil'
     659           output_variable%units         = 'W m-2'
     660
     661        CASE ( 'hfds' )
     662           output_variable%long_name     = 'downward heat flux in soil'
     663           output_variable%standard_name = 'downward_heat_flux_in_soil'
     664           output_variable%units         = 'W m-2'
     665
     666        CASE ( 'hfla' )
     667           output_variable%long_name     = 'upward latent heat flux in air'
     668           output_variable%standard_name = 'upward_latent_heat_flux_in_air'
     669           output_variable%units         = 'W m-2'
     670
     671        CASE ( 'hfsa' )
     672           output_variable%long_name     = 'upward latent heat flux in air'
     673           output_variable%standard_name = 'upward_sensible_heat_flux_in_air'
     674           output_variable%units         = 'W m-2'
     675
     676        CASE ( 'jno2' )
     677           output_variable%long_name     = 'photolysis rate of nitrogen dioxide'
     678           output_variable%standard_name = 'photolysis_rate_of_nitrogen_dioxide'
     679           output_variable%units         = 's-1'
     680
     681        CASE ( 'lwcs' )
     682           output_variable%long_name     = 'liquid water content of soil layer'
     683           output_variable%standard_name = 'liquid_water_content_of_soil_layer'
     684           output_variable%units         = 'kg m-2'
     685
     686        CASE ( 'lwp' )
     687           output_variable%long_name     = 'liquid water path'
     688           output_variable%standard_name = 'atmosphere_mass_content_of_cloud_liquid_water'
     689           output_variable%units         = 'kg m-2'
     690
     691        CASE ( 'ps' )
     692           output_variable%long_name     = 'surface air pressure'
     693           output_variable%standard_name = 'surface_air_pressure'
     694           output_variable%units         = 'hPa'
     695
     696        CASE ( 'pswrtg' )
     697           output_variable%long_name     = 'platform speed wrt ground'
     698           output_variable%standard_name = 'platform_speed_wrt_ground'
     699           output_variable%units         = 'm s-1'
     700
     701        CASE ( 'pswrta' )
     702           output_variable%long_name     = 'platform speed wrt air'
     703           output_variable%standard_name = 'platform_speed_wrt_air'
     704           output_variable%units         = 'm s-1'
     705
     706        CASE ( 'pwv' )
     707           output_variable%long_name     = 'water vapor partial pressure in air'
     708           output_variable%standard_name = 'water_vapor_partial_pressure_in_air'
     709           output_variable%units         = 'hPa'
     710
     711        CASE ( 'ssdu' )
     712           output_variable%long_name     = 'duration of sunshine'
     713           output_variable%standard_name = 'duration_of_sunshine'
     714           output_variable%units         = 's'
     715
     716        CASE ( 't_lw' )
     717           output_variable%long_name     = 'land water temperature'
     718           output_variable%units         = 'degree_C'
     719
     720        CASE ( 'tb' )
     721           output_variable%long_name     = 'brightness temperature'
     722           output_variable%standard_name = 'brightness_temperature'
     723           output_variable%units         = 'K'
     724
     725        CASE ( 'uqv' )
     726           output_variable%long_name     = 'eastward kinematic latent heat flux in air'
     727           output_variable%units         = 'g kg-1 m s-1'
     728
     729        CASE ( 'vqv' )
     730           output_variable%long_name     = 'northward kinematic latent heat flux in air'
     731           output_variable%units         = 'g kg-1 m s-1'
     732
     733        CASE ( 'wqv' )
     734           output_variable%long_name     = 'upward kinematic latent heat flux in air'
     735           output_variable%units         = 'g kg-1 m s-1'
     736
     737        CASE ( 'zcb' )
     738           output_variable%long_name     = 'cloud base altitude'
     739           output_variable%standard_name = 'cloud_base_altitude'
     740           output_variable%units         = 'm'
     741
     742        CASE ( 'zmla' )
     743           output_variable%long_name     = 'atmosphere boundary layer thickness'
     744           output_variable%standard_name = 'atmosphere_boundary_layer_thickness'
     745           output_variable%units         = 'm'
     746
     747        CASE ( 'mcpm1' )
     748           output_variable%long_name     = 'mass concentration of pm1 ambient aerosol particles in air'
     749           output_variable%standard_name = 'mass_concentration_of_pm1_ambient_aerosol_particles_in_air'
     750           output_variable%units         = 'kg m-3'
     751
     752        CASE ( 'mcpm10' )
     753           output_variable%long_name     = 'mass concentration of pm10 ambient aerosol particles in air'
     754           output_variable%standard_name = 'mass_concentration_of_pm10_ambient_aerosol_particles_in_air'
     755           output_variable%units         = 'kg m-3'
     756
     757        CASE ( 'mcpm2p5' )
     758           output_variable%long_name     = 'mass concentration of pm2p5 ambient aerosol particles in air'
     759           output_variable%standard_name = 'mass_concentration_of_pm2p5_ambient_aerosol_particles_in_air'
     760           output_variable%units         = 'kg m-3'
     761
     762        CASE ( 'mfno', 'mcno'  )
     763           output_variable%long_name     = 'mole fraction of nitrogen monoxide in air'
     764           output_variable%standard_name = 'mole_fraction_of_nitrogen_monoxide_in_air'
     765           output_variable%units         = 'ppm' !'mol mol-1'
     766
     767        CASE ( 'mfno2', 'mcno2'  )
     768           output_variable%long_name     = 'mole fraction of nitrogen dioxide in air'
     769           output_variable%standard_name = 'mole_fraction_of_nitrogen_dioxide_in_air'
     770           output_variable%units         = 'ppm' !'mol mol-1'
     771
     772        CASE ( 'tro3'  )
     773           output_variable%long_name     = 'mole fraction of ozone in air'
     774           output_variable%standard_name = 'mole_fraction_of_ozone_in_air'
     775           output_variable%units         = 'ppm' !'mol mol-1'
     776
     777        CASE DEFAULT
     778
     779     END SELECT
     780
     781  END SUBROUTINE vm_set_attributes
     782
     783
    386784!------------------------------------------------------------------------------!
    387785! Description:
    … …  
    390788!------------------------------------------------------------------------------!
    391789 SUBROUTINE vm_parin
    392  
     790
    393791    CHARACTER (LEN=80) ::  line   !< dummy string that contains the current line of the parameter file
    394  
    395     NAMELIST /virtual_measurement_parameters/  dt_virtual_measurement,                             &
    396                                                use_virtual_measurement,                            &
     792
     793    NAMELIST /virtual_measurement_parameters/  dt_virtual_measurement,         &
     794                                               off_ts,                         &
     795                                               off_pr,                         &
     796                                               off_tr,                         &
     797                                               use_virtual_measurement,        &
    397798                                               vm_time_start
    398799
    399800    line = ' '
    400 
    401801!
    402802!-- Try to find stg package
    … …  
    415815!-- Set flag that indicates that the virtual measurement module is switched on
    416816    IF ( use_virtual_measurement )  virtual_measurement = .TRUE.
    417    
     817
    418818    GOTO 20
    419819
    … …  
    423823
    424824 20 CONTINUE
    425  
     825
    426826 END SUBROUTINE vm_parin
    427827
    … …  
    430830! Description:
    431831! ------------
    432 !> Initialize virtual measurements: read coordiante arrays and measured 
     832!> Initialize virtual measurements: read coordiante arrays and measured
    433833!> variables, set indicies indicating the measurement points, read further
    434834!> attributes, etc..
    … …  
    436836 SUBROUTINE vm_init
    437837
    438     USE arrays_3d,                                                             &
    439         ONLY:  zu, zw
    440        
    441     USE grid_variables,                                                        &
    442         ONLY:  ddx, ddy, dx, dy
    443        
    444     USE indices,                                                               &
    445         ONLY:  nxl, nxr, nyn, nys
    446  
    447     USE netcdf_data_input_mod,                                                 &
    448         ONLY:  get_dimension_length,                                           &
    449                init_model,                                                     &
    450                input_file_vm,                                                  &
    451                netcdf_data_input_att,                                          &
    452                netcdf_data_input_var
    453        
    454     IMPLICIT NONE
    455    
    456     CHARACTER(LEN=5)    ::  dum                !< dummy string indicate station id
    457     CHARACTER(LEN=10), DIMENSION(50) ::  measured_variables_file = '' !< array with all measured variables read from NetCDF
    458     CHARACTER(LEN=10), DIMENSION(50) ::  measured_variables      = '' !< dummy array with all measured variables that are allowed   
    459    
    460     INTEGER(iwp) ::  dim_ntime !< dimension size of time coordinate
     838    CHARACTER(LEN=5)                  ::  dum                          !< dummy string indicating station id
     839    CHARACTER(LEN=100), DIMENSION(50) ::  measured_variables_file = '' !< array with all measured variables read from NetCDF
     840    CHARACTER(LEN=100), DIMENSION(50) ::  measured_variables      = '' !< dummy array with all measured variables that are allowed
     841
     842    INTEGER(iwp) ::  dim_ntime !< dimension size of time coordinate
    461843    INTEGER(iwp) ::  i         !< grid index of virtual observation point in x-direction
    462844    INTEGER(iwp) ::  is        !< grid index of real observation point of the respective station in x-direction
    … …  
    471853    INTEGER(iwp) ::  len_char  !< character length of single measured variables without Null character
    472854    INTEGER(iwp) ::  ll        !< running index over all measured variables in file
    473     INTEGER(iwp) ::  lll       !< running index over all allowed variables
     855    INTEGER(iwp) ::  m         !< running index for surface elements
    474856    INTEGER(iwp) ::  n         !< running index over trajectory coordinates
     857    INTEGER(iwp) ::  nofill    !< dummy for nofill return value (not used)
    475858    INTEGER(iwp) ::  ns        !< counter variable for number of observation points on subdomain
     859    INTEGER(iwp) ::  off       !< number of surrounding grid points to be sampled
    476860    INTEGER(iwp) ::  t         !< running index over number of trajectories
    477     INTEGER(iwp) ::  m
    478    
    479     INTEGER(KIND=1)::  soil_dum
    480    
    481     INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ns_all !< dummy array used to sum-up the number of observation coordinates
    482    
     861
     862    INTEGER(KIND=1)                             ::  soil_dum !< dummy variable to input a soil flag
     863
     864    INTEGER(iwp), DIMENSION(:), ALLOCATABLE     ::  ns_all !< dummy array used to sum-up the number of observation coordinates
     865
     866    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE   ::  ns_atmos !< number of observation points for each station on each mpi rank
     867    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE   ::  ns_soil  !< number of observation points for each station on each mpi rank
     868
    483869    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  meas_flag !< mask array indicating measurement positions
    484    
    485 !    LOGICAL ::  chem_include !< flag indicating that chemical species is considered in modelled mechanism
    486     LOGICAL ::  on_pe        !< flag indicating that the respective measurement coordinate is on subdomain
    487    
    488     REAL(wp)     ::  fill_eutm !< _FillValue for coordinate array E_UTM
    489     REAL(wp)     ::  fill_nutm !< _FillValue for coordinate array N_UTM
    490     REAL(wp)     ::  fill_zag  !< _FillValue for height coordinate
    491    
     870
     871    LOGICAL  ::  on_pe        !< flag indicating that the respective measurement coordinate is on subdomain
     872
     873    REAL(wp) ::  fill_eutm !< _FillValue for coordinate array E_UTM
     874    REAL(wp) ::  fill_nutm !< _FillValue for coordinate array N_UTM
     875    REAL(wp) ::  fill_zar  !< _FillValue for height coordinate
     876
    492877    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  e_utm     !< easting UTM coordinate, temporary variable
    493878    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  n_utm     !< northing UTM coordinate, temporary variable
    494879    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  e_utm_tmp !< EUTM coordinate before rotation
    495880    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  n_utm_tmp !< NUTM coordinate before rotation
    496     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  z_ag      !< height coordinate relative to origin_z, temporary variable
    497 !
    498 !-- Obtain number of sites. Also, pass the 'open' string, in order to initially
    499 !-- open the measurement driver.
    500     CALL netcdf_data_input_att( vmea_general%nvm, char_numstations,            &
    501                                 vmea_general%id_vm, input_file_vm,             &
    502                                 global_attribute, 'open', '' )
     881    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  station_h !< station height above reference
     882    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  zar       !< observation height above reference
     883#if defined( __netcdf )
     884!
     885!-- Open the input file.
     886    CALL open_read_file( input_file_vm, pids_id )
     887!
     888!-- Obtain number of sites.
     889    CALL get_attribute( pids_id,                                               &
     890                        char_numstations,                                      &
     891                        vmea_general%nvm,                                      &
     892                        global_attribute )
    503893!
    504894!-- Allocate data structure which encompass all required information, such as
    505 !-- grid points indicies, absolute UTM coordinates, the measured quantities, 
     895!-- grid points indicies, absolute UTM coordinates, the measured quantities,
    506896!-- etc. .
    507897    ALLOCATE( vmea(1:vmea_general%nvm) )
    508898!
    509899!-- Allocate flag array. This dummy array is used to identify grid points
    510 !-- where virtual measurements should be taken. Please note, at least one
    511 !-- ghost point is required, in order to include also the surrounding
    512 !-- grid points of the original coordinate. 
    513     ALLOCATE( meas_flag(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
     900!-- where virtual measurements should be taken. Please note, in order to
     901!-- include also the surrounding grid points of the original coordinate
     902!-- ghost points are required.
     903    ALLOCATE( meas_flag(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    514904    meas_flag = 0
    515905!
    516 !-- Loop over all sites.
     906!-- Loop over all sites in the setup file.
    517907    DO  l = 1, vmea_general%nvm
    518908!
    519909!--    Determine suffix which contains the ID, ordered according to the number
    520 !--    of measurements. 
     910!--    of measurements.
    521911       IF( l < 10 )  THEN
    522912          WRITE( dum, '(I1)')  l
    … …  
    531921       ENDIF
    532922!
    533 !--    Read site coordinates (UTM).
    534        CALL netcdf_data_input_att( vmea(l)%origin_x_obs, char_origx //         &
    535                                    TRIM( dum ), vmea_general%id_vm, '',        &
    536                                    global_attribute, '', '' )
    537        CALL netcdf_data_input_att( vmea(l)%origin_y_obs, char_origy //         &
    538                                    TRIM( dum ), vmea_general%id_vm, '',        &
    539                                    global_attribute, '', '' )
    540 !
    541 !--    Read site name                 
    542        CALL netcdf_data_input_att( vmea(l)%site, char_site // TRIM( dum ),     &
    543                                    vmea_general%id_vm, '', global_attribute,   &
    544                                    '', '' )
     923!--    Read the origin site coordinates (UTM).
     924       CALL get_attribute( pids_id,                                            &
     925                           char_origx // TRIM( dum ),                          &
     926                           vmea(l)%origin_x_obs,                               &
     927                           global_attribute )
     928       CALL get_attribute( pids_id,                                            &
     929                           char_origy // TRIM( dum ),                          &
     930                           vmea(l)%origin_y_obs,                               &
     931                           global_attribute )
     932!
     933!--    Read site name.
     934       CALL get_attribute( pids_id,                                            &
     935                           char_site // TRIM( dum ),                           &
     936                           vmea(l)%site,                                       &
     937                           global_attribute )
     938!
     939!--    Read a flag which indicates that also soil quantities are take at the
     940!--    respective site (is part of the virtual measurement driver).
     941       CALL get_attribute( pids_id,                                            &
     942                           char_soil // TRIM( dum ),                           &
     943                           soil_dum,                                           &
     944                           global_attribute )
     945!
     946!--    Set flag indicating soil-sampling.
     947       IF ( soil_dum == 1 )  vmea(l)%soil_sampling = .TRUE.
    545948!
    546949!--    Read type of the measurement (trajectory, profile, timeseries).
    547        CALL netcdf_data_input_att( vmea(l)%feature_type, char_feature //       &
    548                                    TRIM( dum ), vmea_general%id_vm, '',        &
    549                                    global_attribute, '', '' )
    550 !
    551 !--    Read the name of the original file where observational data is stored.
    552        CALL netcdf_data_input_att( vmea(l)%filename_original, char_filename // &
    553                                    TRIM( dum ), vmea_general%id_vm, '',        &
    554                                    global_attribute, '', '' )
    555 !
    556 !--    Read a flag which indicates that also soil quantities are take at the
    557 !--    respective site (is part of the virtual measurement driver). 
    558        CALL netcdf_data_input_att( soil_dum, char_soil // TRIM( dum ),         &
    559                                    vmea_general%id_vm, '', global_attribute,   &
    560                                    '', '' )
    561 !
    562 !--    Set flag for soil-sampling.
    563        IF ( soil_dum == 1 )  vmea(l)%soil_sampling = .TRUE.
     950       CALL get_attribute( pids_id,                                            &
     951                           char_feature // TRIM( dum ),                        &
     952                           vmea(l)%feature_type,                               &
     953                           global_attribute )
    564954!
    565955!---   Set logicals depending on the type of the measurement
    … …  
    571961          vmea(l)%trajectory        = .TRUE.
    572962!
    573 !--   Give error message in case the type matches non of the pre-defined types.
     963!--    Give error message in case the type matches non of the pre-defined types.
    574964       ELSE
    575965          message_string = 'Attribue featureType = ' //                        &
    576966                           TRIM( vmea(l)%feature_type ) //                     &
    577                            ' is not allowed.' 
     967                           ' is not allowed.'
    578968          CALL message( 'vm_init', 'PA0535', 1, 2, 0, 6, 0 )
    579969       ENDIF
    580970!
    581 !--    Read string with all measured variables at this site
     971!--    Read string with all measured variables at this site.
    582972       measured_variables_file = ''
    583        CALL netcdf_data_input_var( measured_variables_file,                    &
    584                                    char_mv // TRIM( dum ), vmea_general%id_vm )
    585 !
    586 !--    Count the number of measured variables. Only count variables that match
    587 !--    with the allowed variables.
     973       CALL get_variable( pids_id,                                             &
     974                          char_mv // TRIM( dum ),                              &
     975                          measured_variables_file )
     976!
     977!--    Count the number of measured variables.
    588978!--    Please note, for some NetCDF interal reasons characters end with a NULL,
    589979!--    i.e. also empty characters contain a NULL. Therefore, check the strings
    590 !--    for a NULL to get the correct character length in order to compare 
    591 !--    them with the list of allowed variables. 
    592        vmea(l)%nmeas  = 0
     980!--    for a NULL to get the correct character length in order to compare
     981!--    them with the list of allowed variables.
     982       vmea(l)%nmeas  = 1
    593983       DO ll = 1, SIZE( measured_variables_file )
    594984          IF ( measured_variables_file(ll)(1:1) /= CHAR(0)  .AND.              &
    … …  
    602992             ENDDO
    603993             len_char = len_char - 1
    604 !
    605 !--          Now, compare the measured variable with the list of allowed
    606 !--          variables.
    607              DO  lll= 1, SIZE( list_allowed_variables )
    608                 IF ( measured_variables_file(ll)(1:len_char) ==                &
    609                      TRIM( list_allowed_variables(lll) ) )  THEN
    610                    vmea(l)%nmeas = vmea(l)%nmeas + 1
    611                    measured_variables(vmea(l)%nmeas) =                         &
     994
     995             measured_variables(vmea(l)%nmeas) =                               &
    612996                                       measured_variables_file(ll)(1:len_char)
    613                 ENDIF
    614              ENDDO
     997             vmea(l)%nmeas = vmea(l)%nmeas + 1
     998
    615999          ENDIF
    6161000       ENDDO
    617        
    618        
    619 !
    620 !--    Allocate array for the measured variables names for the respective site.
    621        ALLOCATE( vmea(l)%measured_vars_name(1:vmea(l)%nmeas) )
    622 
     1001       vmea(l)%nmeas = vmea(l)%nmeas - 1
     1002!
     1003!--    Allocate data-type array for the measured variables names and attributes
     1004!--    at the respective site.
     1005       ALLOCATE( vmea(l)%var_atts(1:vmea(l)%nmeas) )
     1006!
     1007!--    Store the variable names in a data structures, which assigns further
     1008!--    attributes to this name. Further, for data output reasons, create a
     1009!--    string of output variables, which will be written into the attribute
     1010!--    data_content.
    6231011       DO  ll = 1, vmea(l)%nmeas
    624           vmea(l)%measured_vars_name(ll) = TRIM( measured_variables(ll) )
     1012          vmea(l)%var_atts(ll)%name = TRIM( measured_variables(ll) )
     1013
     1014          vmea(l)%data_content = TRIM( vmea(l)%data_content ) // " " //        &
     1015                                 TRIM( vmea(l)%var_atts(ll)%name )
    6251016       ENDDO
    6261017!
    627 !--    In case of chemistry, check if species is considered in the modelled
    628 !--    chemistry mechanism.
    629 !        IF ( air_chemistry )  THEN
    630 !           DO  ll = 1, vmea(l)%nmeas
    631 !              chem_include = .FALSE.
    632 !              DO  n = 1, nvar
    633 !                 IF ( TRIM( vmea(l)%measured_vars_name(ll) ) ==                 &
    634 !                      TRIM( chem_species(n)%name ) )  chem_include = .TRUE.
    635 !              ENDDO
    636 ! !
    637 ! !--  Revise this. It should only check for chemistry variables and not for all!
    638 !              IF ( .NOT. chem_include )  THEN
    639 !                 message_string = TRIM( vmea(l)%measured_vars_name(ll) ) //     &
    640 !                                  ' is not considered in the modelled '  //     &
    641 !                                  'chemistry mechanism'
    642 !                 CALL message( 'vm_init', 'PA0000', 0, 0, 0, 6, 0 )
    643 !              ENDIF
    644 !           ENDDO
    645 !        ENDIF
    646 !
    647 !--    Read the UTM coordinates for the actual site. Based on the coordinates,
     1018!--    Read all the UTM coordinates for the site. Based on the coordinates,
    6481019!--    define the grid-index space on each subdomain where virtual measurements
    649 !--    should be taken. Note, the entire coordinate arrays will not be stored
    650 !--    as this would exceed memory requirements, particularly for trajectory
    651 !--    measurements.
     1020!--    should be taken. Note, the entire coordinate array (on the entire model
     1021!--    domain) won't be stored as this would exceed memory requirements,
     1022!--    particularly for trajectories.
    6521023       IF ( vmea(l)%nmeas > 0 )  THEN
    6531024!
    654 !--       For stationary measurements UTM coordinates are just one value and 
    655 !--       its dimension is "station", while for mobile measurements UTM 
     1025!--       For stationary measurements UTM coordinates are just one value and
     1026!--       its dimension is "station", while for mobile measurements UTM
    6561027!--       coordinates are arrays depending on the number of trajectories and
    657 !--       time, according to (UC)2 standard. First, inquire dimension length 
     1028!--       time, according to (UC)2 standard. First, inquire dimension length
    6581029!--       of the UTM coordinates.
    6591030          IF ( vmea(l)%trajectory )  THEN
    … …  
    6611032!--          For non-stationary measurements read the number of trajectories
    6621033!--          and the number of time coordinates.
    663              CALL get_dimension_length( vmea_general%id_vm,                    &
    664                                         vmea(l)%ntraj,                         &
    665                                         "traj" //                              &
    666                                         TRIM( dum ) )
    667              CALL get_dimension_length( vmea_general%id_vm,                    &
     1034             CALL get_dimension_length( pids_id,                               &
     1035                                        vmea(l)%n_tr_st,                       &
     1036                                        "traj" // TRIM( dum ) )
     1037             CALL get_dimension_length( pids_id,                               &
    6681038                                        dim_ntime,                             &
    669                                         "ntime" //                             &
    670                                         TRIM( dum ) )
    671 !
    672 !--       For stationary measurements the dimension for UTM and time
    673 !--       coordinates is 1.
     1039                                        "ntime" // TRIM( dum ) )
     1040!
     1041!--       For stationary measurements the dimension for UTM is station
     1042!--       and for the time-coordinate it is one.
    6741043          ELSE
    675              vmea(l)%ntraj  = 1
     1044             CALL get_dimension_length( pids_id,                               &
     1045                                        vmea(l)%n_tr_st,                       &
     1046                                        "station" // TRIM( dum ) )
    6761047             dim_ntime = 1
    6771048          ENDIF
    6781049!
    679 !-        Allocate array which defines individual time frame for each
     1050!-        Allocate array which defines individual time/space frame for each
    6801051!--       trajectory or station.
    681           ALLOCATE( vmea(l)%dim_t(1:vmea(l)%ntraj) )
    682 !
    683 !--       Allocate temporary arrays for UTM and height coordinates. Note, 
     1052          ALLOCATE( vmea(l)%dim_t(1:vmea(l)%n_tr_st) )
     1053!
     1054!--       Allocate temporary arrays for UTM and height coordinates. Note,
    6841055!--       on file UTM coordinates might be 1D or 2D variables
    685           ALLOCATE( e_utm(1:vmea(l)%ntraj,1:dim_ntime) )
    686           ALLOCATE( n_utm(1:vmea(l)%ntraj,1:dim_ntime) )
    687           ALLOCATE( z_ag(1:vmea(l)%ntraj,1:dim_ntime)  )
    688          
    689           ALLOCATE( e_utm_tmp(1:vmea(l)%ntraj,1:dim_ntime) )
    690           ALLOCATE( n_utm_tmp(1:vmea(l)%ntraj,1:dim_ntime) )
    691 !
    692 !--       Read _FillValue attributes of the coordinate dimensions.
    693           CALL netcdf_data_input_att( fill_eutm, char_fillvalue,               &
    694                                       vmea_general%id_vm, '',                  &
    695                                       .NOT. global_attribute, '',              &
    696                                       char_eutm // TRIM( dum ) )
    697           CALL netcdf_data_input_att( fill_nutm, char_fillvalue,               &
    698                                       vmea_general%id_vm, '',                  &
    699                                       .NOT. global_attribute, '',              &
    700                                       char_nutm // TRIM( dum ) )
    701           CALL netcdf_data_input_att( fill_zag, char_fillvalue,                &
    702                                       vmea_general%id_vm, '',                  &
    703                                       .NOT. global_attribute, '',              &
    704                                       char_zag  // TRIM( dum ) )
    705 !
    706 !--       Read UTM and height coordinates coordinates for all trajectories and
    707 !--       times.
     1056          ALLOCATE( e_utm(1:vmea(l)%n_tr_st,1:dim_ntime)       )
     1057          ALLOCATE( n_utm(1:vmea(l)%n_tr_st,1:dim_ntime)       )
     1058          ALLOCATE( station_h(1:vmea(l)%n_tr_st,1:dim_ntime)   )
     1059          ALLOCATE( zar(1:vmea(l)%n_tr_st,1:dim_ntime)         )
     1060          e_utm     = 0.0_wp
     1061          n_utm     = 0.0_wp
     1062          station_h = 0.0_wp
     1063          zar       = 0.0_wp
     1064
     1065          ALLOCATE( e_utm_tmp(1:vmea(l)%n_tr_st,1:dim_ntime) )
     1066          ALLOCATE( n_utm_tmp(1:vmea(l)%n_tr_st,1:dim_ntime) )
     1067!
     1068!--       Read UTM and height coordinates coordinates for all trajectories and
     1069!--       times. Note, in case these obtain any missing values, replace them
     1070!--       with default _FillValues.
     1071          CALL inquire_fill_value( pids_id,                                    &
     1072                                   char_eutm // TRIM( dum ),                   &
     1073                                   nofill,                                     &
     1074                                   fill_eutm )
     1075          CALL inquire_fill_value( pids_id,                                    &
     1076                                   char_nutm // TRIM( dum ),                   &
     1077                                   nofill,                                     &
     1078                                   fill_nutm )
     1079          CALL inquire_fill_value( pids_id,                                    &
     1080                                   char_zar // TRIM( dum ),                    &
     1081                                   nofill,                                     &
     1082                                   fill_zar )
     1083!
     1084!--       Further line is just to avoid compiler warnings. nofill might be used
     1085!--       in future.
     1086          IF ( nofill == 0  .OR.  nofill /= 0 )  CONTINUE
     1087!
     1088!--       Read observation coordinates. Please note, for trajectories the
     1089!--       observation height is stored directly in z, while for timeSeries
     1090!--       it is stored in z - station_h, according to UC2-standard.
    7081091          IF ( vmea(l)%trajectory )  THEN
    709              CALL netcdf_data_input_var( e_utm, char_eutm // TRIM( dum ),      &
    710                                          vmea_general%id_vm,                   &
    711                                          0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
    712              CALL netcdf_data_input_var( n_utm, char_nutm // TRIM( dum ),      &
    713                                          vmea_general%id_vm,                   &
    714                                          0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
    715              CALL netcdf_data_input_var( z_ag, char_zag // TRIM( dum ),        &
    716                                          vmea_general%id_vm,                   &
    717                                          0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
     1092             CALL get_variable( pids_id,                                       &
     1093                                char_eutm // TRIM( dum ),                      &
     1094                                e_utm,                                         &
     1095                                0, dim_ntime-1,                                &
     1096                                0, vmea(l)%n_tr_st-1 )
     1097             CALL get_variable( pids_id,                                       &
     1098                                char_nutm // TRIM( dum ),                      &
     1099                                n_utm,                                         &
     1100                                0, dim_ntime-1,                                &
     1101                                0, vmea(l)%n_tr_st-1 )
     1102             CALL get_variable( pids_id,                                       &
     1103                                char_zar // TRIM( dum ),                       &
     1104                                zar,                                           &
     1105                                0, dim_ntime-1,                                &
     1106                                0, vmea(l)%n_tr_st-1 )
    7181107          ELSE
    719              CALL netcdf_data_input_var( e_utm(1,:), char_eutm // TRIM( dum ), &
    720                                          vmea_general%id_vm )
    721              CALL netcdf_data_input_var( n_utm(1,:), char_nutm // TRIM( dum ), &
    722                                          vmea_general%id_vm )
    723              CALL netcdf_data_input_var( z_ag(1,:),  char_zag  // TRIM( dum ), &
    724                                          vmea_general%id_vm )
    725           ENDIF         
     1108             CALL get_variable( pids_id,                                       &
     1109                                char_eutm // TRIM( dum ),                      &
     1110                                e_utm(:,1) )
     1111             CALL get_variable( pids_id,                                       &
     1112                                char_nutm // TRIM( dum ),                      &
     1113                                n_utm(:,1) )
     1114             CALL get_variable( pids_id,                                       &
     1115                                char_station_h // TRIM( dum ),                 &
     1116                                station_h(:,1) )
     1117             CALL get_variable( pids_id,                                       &
     1118                                char_zar // TRIM( dum ),                       &
     1119                                zar(:,1) )
     1120          ENDIF
     1121
     1122          e_utm = MERGE( e_utm, vmea(l)%fillout, e_utm /= fill_eutm )
     1123          n_utm = MERGE( n_utm, vmea(l)%fillout, n_utm /= fill_nutm )
     1124          zar   = MERGE( zar,   vmea(l)%fillout, zar   /= fill_zar  )
     1125!
     1126!--       Compute observation height above ground.
     1127          zar  = zar - station_h
    7261128!
    7271129!--       Based on UTM coordinates, check if the measurement station or parts
    728 !--       of the trajectory is on subdomain. This case, setup grid index space
    729 !--       sample these quantities. 
     1130!--       of the trajectory are on subdomain. This case, setup grid index space
     1131!--       sample these quantities.
    7301132          meas_flag = 0
    731           DO  t = 1, vmea(l)%ntraj
    732 !             
    733 !--          First, compute relative x- and y-coordinates with respect to the 
    734 !--          lower-left origin of the model domain, which is the difference 
     1133          DO  t = 1, vmea(l)%n_tr_st
     1134!
     1135!--          First, compute relative x- and y-coordinates with respect to the
     1136!--          lower-left origin of the model domain, which is the difference
    7351137!--          between UTM coordinates. Note, if the origin is not correct, the
    736 !--          virtual sites will be misplaced.
     1138!--          virtual sites will be misplaced. Further, in case of an rotated
     1139!--          model domain, the UTM coordinates must be also rotated.
    7371140             e_utm_tmp(t,1:dim_ntime) = e_utm(t,1:dim_ntime) - init_model%origin_x
    7381141             n_utm_tmp(t,1:dim_ntime) = n_utm(t,1:dim_ntime) - init_model%origin_y
    … …  
    7491152!--          trajectory. This is required as several stations and trajectories
    7501153!--          are merged into one file but they do not have the same number of
    751 !--          points in time, hence, missing values may occur and cannot be 
     1154!--          points in time, hence, missing values may occur and cannot be
    7521155!--          processed further. This is actually a work-around for the specific
    753 !--          (UC)2 dataset, but it won't harm in anyway.
     1156!--          (UC)2 dataset, but it won't harm anyway.
    7541157             vmea(l)%dim_t(t) = 0
    7551158             DO  n = 1, dim_ntime
    7561159                IF ( e_utm(t,n) /= fill_eutm  .AND.                            &
    7571160                     n_utm(t,n) /= fill_nutm  .AND.                            &
    758                      z_ag(t,n)  /= fill_zag )  vmea(l)%dim_t(t) = n
     1161                     zar(t,n)   /= fill_zar )  vmea(l)%dim_t(t) = n
    7591162             ENDDO
    7601163!
    7611164!--          Compute grid indices relative to origin and check if these are
    762 !--          on the subdomain. Note, virtual measurements will be taken also 
    763 !--          at grid points surrounding the station, hence, check also for 
     1165!--          on the subdomain. Note, virtual measurements will be taken also
     1166!--          at grid points surrounding the station, hence, check also for
    7641167!--          these grid points.
     1168!--          The number of surrounding grid points is set according to the
     1169!--          featureType.
     1170             IF ( vmea(l)%timseries_profile )  THEN
     1171                off = off_pr
     1172             ELSEIF ( vmea(l)%timseries     )  THEN
     1173                off = off_ts
     1174             ELSEIF ( vmea(l)%trajectory    )  THEN
     1175                off = off_tr
     1176             ENDIF
     1177
    7651178             DO  n = 1, vmea(l)%dim_t(t)
    7661179                is = INT( ( e_utm(t,n) + 0.5_wp * dx ) * ddx, KIND = iwp )
    767                 js = INT( ( n_utm(t,n) + 0.5_wp * dy ) * ddy, KIND = iwp )             
     1180                js = INT( ( n_utm(t,n) + 0.5_wp * dy ) * ddy, KIND = iwp )
    7681181!
    7691182!--             Is the observation point on subdomain?
    … …  
    7771190!--                height.
    7781191                   ksurf = topo_top_ind(js,is,0)
    779                    ks = MINLOC( ABS( zu - zw(ksurf) - z_ag(t,n) ), DIM = 1 ) - 1
     1192                   ks = MINLOC( ABS( zu - zw(ksurf) - zar(t,n) ), DIM = 1 ) - 1
    7801193!
    7811194!--                Set mask array at the observation coordinates. Also, flag the
    7821195!--                surrounding coordinate points, but first check whether the
    783 !--                surrounding coordinate points are on the subdomain. 
    784                    kl = MERGE( ks-1, ks, ks-1 >= nzb  .AND. ks-1 >= ksurf )
    785                    ku = MERGE( ks+1, ks, ks+1 < nzt+1 )
    786                  
    787                    DO  i = is-1, is+1
    788                       DO  j = js-1, js+1
     1196!--                surrounding coordinate points are on the subdomain.
     1197                   kl = MERGE( ks-off, ksurf, ks-off >= nzb  .AND. ks-off >= ksurf )
     1198                   ku = MERGE( ks+off, nzt,   ks+off < nzt+1 )
     1199
     1200                   DO  i = is-off, is+off
     1201                      DO  j = js-off, js+off
    7891202                         DO  k = kl, ku
    790                             meas_flag(k,j,i) = MERGE(                              &
    791                                              IBSET( meas_flag(k,j,i), 0 ),         &
    792                                              0,                                    &
    793                                              BTEST( wall_flags_total_0(k,j,i), 0 ) &
     1203                            meas_flag(k,j,i) = MERGE(                           &
     1204                                          IBSET( meas_flag(k,j,i), 0 ),         &
     1205                                          0,                                    &
     1206                                          BTEST( wall_flags_total_0(k,j,i), 0 ) &
    7941207                                                    )
    7951208                         ENDDO
    … …  
    7981211                ENDIF
    7991212             ENDDO
    800              
     1213
    8011214          ENDDO
    8021215!
    803 !--       Based on the flag array count the number of of sampling coordinates.
    804 !--       Please note, sampling coordinates in atmosphere and soil may be 
    805 !--       different, as within the soil all levels will be measured.           
     1216!--       Based on the flag array count the number of sampling coordinates.
     1217!--       Please note, sampling coordinates in atmosphere and soil may be
     1218!--       different, as within the soil all levels will be measured.
    8061219!--       Hence, count individually. Start with atmoshere.
    8071220          ns = 0
    808           DO  i = nxl-1, nxr+1
    809              DO  j = nys-1, nyn+1
     1221          DO  i = nxl-off, nxr+off
     1222             DO  j = nys-off, nyn+off
    8101223                DO  k = nzb, nzt+1
    8111224                   ns = ns + MERGE( 1, 0, BTEST( meas_flag(k,j,i), 0 ) )
    … …  
    8131226             ENDDO
    8141227          ENDDO
    815          
     1228
    8161229!
    8171230!--       Store number of observation points on subdomain and allocate index
    8181231!--       arrays as well as array containing height information.
    8191232          vmea(l)%ns = ns
    820          
     1233
    8211234          ALLOCATE( vmea(l)%i(1:vmea(l)%ns) )
    8221235          ALLOCATE( vmea(l)%j(1:vmea(l)%ns) )
    8231236          ALLOCATE( vmea(l)%k(1:vmea(l)%ns) )
    824           ALLOCATE( vmea(l)%z_ag(1:vmea(l)%ns) )     
    825 !
    826 !--       Based on the flag array store the grid indices which correspond to 
    827 !--       the observation coordinates. 
     1237          ALLOCATE( vmea(l)%zar(1:vmea(l)%ns) )
     1238!
     1239!--       Based on the flag array store the grid indices which correspond to
     1240!--       the observation coordinates.
    8281241          ns = 0
    829           DO  i = nxl-1, nxr+1
    830              DO  j = nys-1, nyn+1
     1242          DO  i = nxl-off, nxr+off
     1243             DO  j = nys-off, nyn+off
    8311244                DO  k = nzb, nzt+1
    8321245                   IF ( BTEST( meas_flag(k,j,i), 0 ) )  THEN
    … …  
    8351248                      vmea(l)%j(ns) = j
    8361249                      vmea(l)%k(ns) = k
    837                       vmea(l)%z_ag(ns)  = zu(k) - zw(topo_top_ind(j,i,0))
     1250                      vmea(l)%zar(ns)  = zu(k) - zw(topo_top_ind(j,i,0))
    8381251                   ENDIF
    8391252                ENDDO
    … …  
    8411254          ENDDO
    8421255!
    843 !--       Same for the soil. Based on the flag array, count the number of
    844 !--       sampling coordinates in soil. Sample at all soil levels in this case.
     1256!--       Same for the soil. Based on the flag array, count the number of
     1257!--       sampling coordinates in soil. Sample at all soil levels in this case.
     1258!--       Please note, soil variables can only be sampled on subdomains, not
     1259!--       on ghost layers.
    8451260          IF ( vmea(l)%soil_sampling )  THEN
    8461261             DO  i = nxl, nxr
    … …  
    8501265                           surf_lsm_h%end_index(j,i) )  THEN
    8511266                         vmea(l)%ns_soil = vmea(l)%ns_soil +                   &
    852                                                       nzt_soil - nzb_soil + 1 
     1267                                                      nzt_soil - nzb_soil + 1
    8531268                      ENDIF
    8541269                      IF ( surf_usm_h%start_index(j,i) <=                      &
    8551270                           surf_usm_h%end_index(j,i) )  THEN
    8561271                         vmea(l)%ns_soil = vmea(l)%ns_soil +                   &
    857                                                       nzt_wall - nzb_wall + 1 
     1272                                                      nzt_wall - nzb_wall + 1
    8581273                      ENDIF
    8591274                   ENDIF
    8601275                ENDDO
    8611276             ENDDO
    862           ENDIF         
    863 !
    864 !--       Allocate index arrays as well as array containing height information 
     1277          ENDIF
     1278!
     1279!--       Allocate index arrays as well as array containing height information
    8651280!--       for soil.
    8661281          IF ( vmea(l)%soil_sampling )  THEN
    … …  
    8681283             ALLOCATE( vmea(l)%j_soil(1:vmea(l)%ns_soil) )
    8691284             ALLOCATE( vmea(l)%k_soil(1:vmea(l)%ns_soil) )
    870              ALLOCATE( vmea(l)%depth(1:vmea(l)%ns_soil) )
    871           ENDIF     
     1285             ALLOCATE( vmea(l)%depth(1:vmea(l)%ns_soil)  )
     1286          ENDIF
    8721287!
    8731288!--       For soil, store the grid indices.
    … …  
    8851300                            vmea(l)%j_soil(ns) = j
    8861301                            vmea(l)%k_soil(ns) = k
    887                             vmea(l)%depth(ns)  = zs(k)
     1302                            vmea(l)%depth(ns)  = - zs(k)
    8881303                         ENDDO
    8891304                      ENDIF
    890                      
     1305
    8911306                      IF ( surf_usm_h%start_index(j,i) <=                      &
    8921307                           surf_usm_h%end_index(j,i) )  THEN
    … …  
    8971312                            vmea(l)%j_soil(ns) = j
    8981313                            vmea(l)%k_soil(ns) = k
    899                             vmea(l)%depth(ns)  = surf_usm_h%zw(k,m)
     1314                            vmea(l)%depth(ns)  = - surf_usm_h%zw(k,m)
    9001315                         ENDDO
    9011316                      ENDIF
    … …  
    9071322!--       Allocate array to save the sampled values.
    9081323          ALLOCATE( vmea(l)%measured_vars(1:vmea(l)%ns,1:vmea(l)%nmeas) )
    909          
     1324
    9101325          IF ( vmea(l)%soil_sampling )                                         &
    9111326             ALLOCATE( vmea(l)%measured_vars_soil(1:vmea(l)%ns_soil,           &
    … …  
    9241339          IF ( ALLOCATED( n_utm_tmp ) )  DEALLOCATE( n_utm_tmp )
    9251340          IF ( ALLOCATED( n_utm )     )  DEALLOCATE( n_utm )
    926           IF ( ALLOCATED( z_ag  )     )  DEALLOCATE( z_ag  )
    927           IF ( ALLOCATED( z_ag  )     )  DEALLOCATE( vmea(l)%dim_t )
     1341          IF ( ALLOCATED( zar  )      )  DEALLOCATE( vmea(l)%dim_t )
     1342          IF ( ALLOCATED( zar  )      )  DEALLOCATE( zar  )
     1343          IF ( ALLOCATED( station_h ) )  DEALLOCATE( station_h )
     1344
    9281345       ENDIF
    9291346    ENDDO
    9301347!
    931 !-- Close input file for virtual measurements. Therefore, just call
    932 !-- the read attribute routine with the "close" option.
    933     CALL netcdf_data_input_att( vmea_general%nvm, char_numstations,            &
    934                                 vmea_general%id_vm, '',                        &
    935                                 global_attribute, 'close', '' )
    936 !
    937 !-- Sum-up the number of observation coordiates, for atmosphere first. 
     1348!-- Dellocate flag array
     1349    DEALLOCATE( meas_flag )
     1350!
     1351!-- Close input file for virtual measurements.
     1352    CALL close_input_file( pids_id )
     1353!
     1354!-- Sum-up the number of observation coordiates, for atmosphere first.
    9381355!-- This is actually only required for data output.
    9391356    ALLOCATE( ns_all(1:vmea_general%nvm) )
    940     ns_all = 0   
     1357    ns_all = 0
    9411358#if defined( __parallel )
    9421359    CALL MPI_ALLREDUCE( vmea(:)%ns, ns_all(:), vmea_general%nvm, MPI_INTEGER,  &
    … …  
    9481365!
    9491366!-- Now for soil
    950     ns_all = 0   
     1367    ns_all = 0
    9511368#if defined( __parallel )
    9521369    CALL MPI_ALLREDUCE( vmea(:)%ns_soil, ns_all(:), vmea_general%nvm,          &
    … …  
    9561373#endif
    9571374    vmea(:)%ns_soil_tot = ns_all(:)
    958    
     1375
    9591376    DEALLOCATE( ns_all )
    960 !                               
    961 !-- Dellocate flag array
    962     DEALLOCATE( meas_flag )
    963 !
    964 !-- Initialize binary data output of virtual measurements.
    965 !-- Open binary output file.
    966     CALL check_open( 27 )
    967 !
    968 !-- Output header information.
    969     CALL vm_data_output
    970        
     1377!
     1378!-- In case of parallel NetCDF the start coordinate for each mpi rank needs to
     1379!-- be defined, so that each processor knows where to write the data.
     1380#if defined( __netcdf4_parallel )
     1381    ALLOCATE( ns_atmos(0:numprocs-1,1:vmea_general%nvm) )
     1382    ALLOCATE( ns_soil(0:numprocs-1,1:vmea_general%nvm)  )
     1383    ns_atmos = 0
     1384    ns_soil  = 0
     1385
     1386    DO  l = 1, vmea_general%nvm
     1387       ns_atmos(myid,l) = vmea(l)%ns
     1388       ns_soil(myid,l)  = vmea(l)%ns_soil
     1389    ENDDO
     1390
     1391#if defined( __parallel )
     1392    CALL MPI_ALLREDUCE( MPI_IN_PLACE, ns_atmos, numprocs * vmea_general%nvm,   &
     1393                        MPI_INTEGER, MPI_SUM, comm2d, ierr )
     1394    CALL MPI_ALLREDUCE( MPI_IN_PLACE, ns_soil, numprocs * vmea_general%nvm,    &
     1395                        MPI_INTEGER, MPI_SUM, comm2d, ierr )
     1396#else
     1397    ns_atmos(0,:) = vmea(:)%ns
     1398    ns_soil(0,:)  = vmea(:)%ns_soil
     1399#endif
     1400
     1401!
     1402!-- Determine the start coordinate in NetCDF file for the local arrays.
     1403!-- Note, start coordinates are initialized with zero for sake of simplicity
     1404!-- in summation. However, in NetCDF the start coordinates must be >= 1,
     1405!-- so that a one needs to be added at the end.
     1406    DO  l = 1, vmea_general%nvm
     1407       DO  n  = 0, myid - 1
     1408          vmea(l)%start_coord_a = vmea(l)%start_coord_a + ns_atmos(n,l)
     1409          vmea(l)%start_coord_s = vmea(l)%start_coord_s + ns_soil(n,l)
     1410       ENDDO
     1411!
     1412!--    Start coordinate in NetCDF starts always at one not at 0.
     1413       vmea(l)%start_coord_a = vmea(l)%start_coord_a + 1
     1414       vmea(l)%start_coord_s = vmea(l)%start_coord_s + 1
     1415!
     1416!--    Determine the local end coordinate
     1417       vmea(l)%end_coord_a = vmea(l)%start_coord_a + vmea(l)%ns - 1
     1418       vmea(l)%end_coord_s = vmea(l)%start_coord_s + vmea(l)%ns_soil - 1
     1419    ENDDO
     1420
     1421    DEALLOCATE( ns_atmos )
     1422    DEALLOCATE( ns_soil  )
     1423
     1424#endif
     1425
     1426#endif
     1427
    9711428  END SUBROUTINE vm_init
    972  
    973  
     1429
     1430
    9741431!------------------------------------------------------------------------------!
    9751432! Description:
    9761433! ------------
    977 !> Binary data output.
     1434!> Initialize output using data-output module
    9781435!------------------------------------------------------------------------------!
    979   SUBROUTINE vm_data_output
    980    
    981      USE pegrid
    982    
    983      IMPLICIT NONE
    984          
    985      INTEGER(iwp) ::  i  !< running index over IO blocks   
    986      INTEGER(iwp) ::  l  !< running index over all stations
    987      INTEGER(iwp) ::  n  !< running index over all measured variables at a station
    988 !
    989 !--  Header output on each PE
    990      IF ( init )  THEN
    991 
    992         DO  i = 0, io_blocks-1
    993            IF ( i == io_group )  THEN
    994               WRITE ( 27 )  'number of measurements            '
    995               WRITE ( 27 )  vmea_general%nvm
    996 
    997               DO  l = 1, vmea_general%nvm
    998                  WRITE ( 27 )  'site                              '
    999                  WRITE ( 27 )  vmea(l)%site
    1000                  WRITE ( 27 )  'file                              '
    1001                  WRITE ( 27 )  vmea(l)%filename_original
    1002                  WRITE ( 27 )  'feature_type                      '
    1003                  WRITE ( 27 )  vmea(l)%feature_type
    1004                  WRITE ( 27 )  'origin_x_obs                      '
    1005                  WRITE ( 27 )  vmea(l)%origin_x_obs
    1006                  WRITE ( 27 )  'origin_y_obs                      '
    1007                  WRITE ( 27 )  vmea(l)%origin_y_obs
    1008                  WRITE ( 27 )  'total number of observation points'                               
    1009                  WRITE ( 27 )  vmea(l)%ns_tot
    1010                  WRITE ( 27 )  'number of measured variables      '
    1011                  WRITE ( 27 )  vmea(l)%nmeas
    1012                  WRITE ( 27 )  'variables                         '
    1013                  WRITE ( 27 )  vmea(l)%measured_vars_name(:)
    1014                  WRITE ( 27 )  'number of observation points      '
    1015                  WRITE ( 27 )  vmea(l)%ns
    1016                  WRITE ( 27 )  'E_UTM                             '
    1017                  WRITE ( 27 )  init_model%origin_x                             &
     1436 SUBROUTINE vm_init_output
     1437
     1438    CHARACTER(LEN=100) ::  variable_name  !< name of output variable
     1439
     1440    INTEGER(iwp) ::  l              !< loop index
     1441    INTEGER(iwp) ::  n              !< loop index
     1442    INTEGER      ::  return_value   !< returned status value of called function
     1443
     1444    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ndim !< dummy to write dimension
     1445
     1446    REAL(wp) ::  dum_lat !< transformed geographical coordinate (latitude)
     1447    REAL(wp) ::  dum_lon !< transformed geographical coordinate (longitude)
     1448
     1449!
     1450!-- Determine the number of output timesteps. Set a maximum value of 80000
     1451!-- timesteps.
     1452    ntimesteps = MIN( CEILING(                                                 &
     1453                  ( end_time - MAX( vm_time_start, time_since_reference_point )&
     1454                  ) / MAX( 0.1_wp, dt_virtual_measurement ) ), ntimesteps_max )
     1455!
     1456!-- Create directory where output files will be stored.
     1457    CALL local_system( 'mkdir -p VM_OUTPUT' // TRIM( coupling_char ) )
     1458!
     1459!-- Loop over all sites.
     1460    DO  l = 1, vmea_general%nvm
     1461!
     1462!--    Skip if no observations will be taken for this site.
     1463       IF ( vmea(l)%ns_tot == 0  .AND.  vmea(l)%ns_soil_tot == 0 )  CYCLE
     1464!
     1465!--    Define output file.
     1466       WRITE( vmea(l)%nc_filename, '(A,I4.4)')  'VM_OUTPUT' //                 &
     1467                                                TRIM( coupling_char ) // '/' //&
     1468                                                'site', l
     1469
     1470
     1471       return_value = dom_def_file( vmea(l)%nc_filename, 'netcdf4-parallel' )
     1472!
     1473!--    Define global attributes.
     1474!--    Before, transform UTM into geographical coordinates.
     1475       CALL convert_utm_to_geographic( crs_list,                               &
     1476                                       vmea(l)%origin_x_obs,                   &
     1477                                       vmea(l)%origin_y_obs,                   &
     1478                                       dum_lon,                                &
     1479                                       dum_lat )
     1480
     1481       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1482                                   attribute_name = 'site',                    &
     1483                                   value = TRIM( vmea(l)%site ) )
     1484       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1485                                   attribute_name = 'title',                   &
     1486                                   value = 'Virtual measurement output')
     1487       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1488                                   attribute_name = 'source',                  &
     1489                                   value = 'PALM-4U')
     1490       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1491                                   attribute_name = 'institution',             &
     1492                                   value = input_file_atts%institution )
     1493       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1494                                   attribute_name = 'acronym',                 &
     1495                                   value = input_file_atts%acronym )
     1496       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1497                                   attribute_name = 'author',                  &
     1498                                   value = input_file_atts%author )
     1499       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1500                                   attribute_name = 'contact_person',          &
     1501                                   value = input_file_atts%contact_person )
     1502       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1503                                   attribute_name = 'iop',                     &
     1504                                   value = input_file_atts%campaign )
     1505       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1506                                   attribute_name = 'campaign',                &
     1507                                   value = 'PALM-4U' )
     1508       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1509                                   attribute_name = 'origin_time ',            &
     1510                                   value = origin_date_time)
     1511       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1512                                   attribute_name = 'location',                &
     1513                                   value = input_file_atts%location )
     1514       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1515                                   attribute_name = 'origin_x',                &
     1516                                   value = vmea(l)%origin_x_obs )
     1517       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1518                                   attribute_name = 'origin_y',                &
     1519                                   value = vmea(l)%origin_y_obs )
     1520       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1521                                   attribute_name = 'origin_lon',              &
     1522                                   value = dum_lon )
     1523       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1524                                   attribute_name = 'origin_lat',              &
     1525                                   value = dum_lat )
     1526       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1527                                   attribute_name = 'origin_z',                &
     1528                                   value = 0.0 )
     1529       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1530                                   attribute_name = 'rotation_angle',          &
     1531                                   value = input_file_atts%rotation_angle )
     1532       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1533                                   attribute_name = 'featureType',             &
     1534                                   value = TRIM( vmea(l)%feature_type_out ) )
     1535       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1536                                   attribute_name = 'data_content',            &
     1537                                   value = TRIM( vmea(l)%data_content ) )
     1538       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1539                                   attribute_name = 'creation_time',           &
     1540                                   value = input_file_atts%creation_time )
     1541       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1542                                   attribute_name = 'version',                 &
     1543                                   value = 1 ) !input_file_atts%version )
     1544       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1545                                   attribute_name = 'creation_time',           &
     1546                                   value = TRIM( vmea(l)%site ) )
     1547       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1548                                   attribute_name = 'Conventions',             &
     1549                                   value = input_file_atts%conventions )
     1550       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1551                                   attribute_name = 'dependencies',            &
     1552                                   value = input_file_atts%dependencies )
     1553       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1554                                   attribute_name = 'history',                 &
     1555                                   value = input_file_atts%history )
     1556       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1557                                   attribute_name = 'references',              &
     1558                                   value = input_file_atts%references )
     1559       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1560                                   attribute_name = 'comment',                 &
     1561                                   value = input_file_atts%comment )
     1562       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1563                                   attribute_name = 'keywords',                &
     1564                                   value = input_file_atts%keywords )
     1565       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1566                                   attribute_name = 'licence',                 &
     1567                                   value = '[UC]2 Open Licence; see [UC]2 ' // &
     1568                                           'data policy available at ' //      &
     1569                                           'www.uc2-program.org/uc2_data_policy.pdf' )
     1570!
     1571!--    Define dimensions.
     1572!--    station
     1573       ALLOCATE( ndim(1:vmea(l)%ns_tot) )
     1574       DO  n = 1, vmea(l)%ns_tot
     1575          ndim(n) = n
     1576       ENDDO
     1577       return_value = dom_def_dim( vmea(l)%nc_filename,                        &
     1578                                   dimension_name = 'station',                 &
     1579                                   output_type = 'int32',                      &
     1580                                   bounds = (/1_iwp, vmea(l)%ns_tot/),         &
     1581                                   values_int32 = ndim )
     1582       DEALLOCATE( ndim )
     1583!
     1584!--    ntime
     1585       ALLOCATE( ndim(1:ntimesteps) )
     1586       DO  n = 1, ntimesteps
     1587          ndim(n) = n
     1588       ENDDO
     1589
     1590       return_value = dom_def_dim( vmea(l)%nc_filename,                        &
     1591                                   dimension_name = 'ntime',                   &
     1592                                   output_type = 'int32',                      &
     1593                                   bounds = (/1_iwp, ntimesteps/),             &
     1594                                   values_int32 = ndim )
     1595       DEALLOCATE( ndim )
     1596!
     1597!--    nv
     1598       ALLOCATE( ndim(1:2) )
     1599       DO  n = 1, 2
     1600          ndim(n) = n
     1601       ENDDO
     1602
     1603       return_value = dom_def_dim( vmea(l)%nc_filename,                        &
     1604                                   dimension_name = 'nv',                      &
     1605                                   output_type = 'int32',                      &
     1606                                   bounds = (/1_iwp, 2_iwp/),                  &
     1607                                   values_int32 = ndim )
     1608       DEALLOCATE( ndim )
     1609!
     1610!--    maximum name length
     1611       ALLOCATE( ndim(1:maximum_name_length) )
     1612       DO  n = 1, maximum_name_length
     1613          ndim(n) = n
     1614       ENDDO
     1615
     1616       return_value = dom_def_dim( vmea(l)%nc_filename,                        &
     1617                                   dimension_name = 'max_name_len',            &
     1618                                   output_type = 'int32',                      &
     1619                                   bounds = (/1_iwp, 32_iwp/),                 &
     1620                                   values_int32 = ndim )
     1621       DEALLOCATE( ndim )
     1622!
     1623!--    Define coordinate variables.
     1624!--    time
     1625       variable_name = 'time'
     1626       return_value = dom_def_var( vmea(l)%nc_filename,                        &
     1627                                   variable_name = variable_name,              &
     1628                                   dimension_names = (/ 'station  ',           &
     1629                                                        'ntime    '/),         &
     1630                                   output_type = 'real32' )
     1631!
     1632!--    station_name. DOM needs to be enabled to define CHARACTER variables.
     1633!        variable_name = 'station_name'
     1634!        return_value = dom_def_var( vmea(l)%nc_filename,                        &
     1635!                                    variable_name = variable_name,              &
     1636!                                    dimension_names = (/ 'max_name_len',        &
     1637!                                                         'station     '/),      &
     1638!                                    output_type = 'char' )
     1639!
     1640!--    vrs (vertical reference system)
     1641       variable_name = 'vrs'
     1642       return_value = dom_def_var( vmea(l)%nc_filename,                        &
     1643                                   variable_name = variable_name,              &
     1644                                   dimension_names = (/ 'station' /),          &
     1645                                   output_type = 'int8' )
     1646!
     1647!--    crs (coordinate reference system)
     1648       variable_name = 'crs'
     1649       return_value = dom_def_var( vmea(l)%nc_filename,                        &
     1650                                   variable_name = variable_name,              &
     1651                                   dimension_names = (/ 'station' /),          &
     1652                                   output_type = 'int8' )
     1653!
     1654!--    z
     1655       variable_name = 'z'
     1656       return_value = dom_def_var( vmea(l)%nc_filename,                        &
     1657                                   variable_name = variable_name,              &
     1658                                   dimension_names = (/'station'/),            &
     1659                                   output_type = 'real32' )
     1660!
     1661!--    station_h
     1662       variable_name = 'station_h'
     1663       return_value = dom_def_var( vmea(l)%nc_filename,                        &
     1664                                   variable_name = variable_name,              &
     1665                                   dimension_names = (/'station'/),            &
     1666                                   output_type = 'real32' )
     1667!
     1668!--    x
     1669       variable_name = 'x'
     1670       return_value = dom_def_var( vmea(l)%nc_filename,                        &
     1671                                   variable_name = variable_name,              &
     1672                                   dimension_names = (/'station'/),            &
     1673                                   output_type = 'real32' )
     1674!
     1675!--    y
     1676       variable_name = 'y'
     1677       return_value = dom_def_var( vmea(l)%nc_filename,                        &
     1678                                   variable_name = variable_name,              &
     1679                                   dimension_names = (/'station'/),            &
     1680                                   output_type = 'real32' )
     1681!
     1682!--    E-UTM
     1683       variable_name = 'E_UTM'
     1684       return_value = dom_def_var( vmea(l)%nc_filename,                        &
     1685                                   variable_name = variable_name,              &
     1686                                   dimension_names = (/'station'/),            &
     1687                                   output_type = 'real32' )
     1688!
     1689!--    N-UTM
     1690       variable_name = 'N_UTM'
     1691       return_value = dom_def_var( vmea(l)%nc_filename,                        &
     1692                                   variable_name = variable_name,              &
     1693                                   dimension_names = (/'station'/),            &
     1694                                   output_type = 'real32' )
     1695!
     1696!--    latitude
     1697       variable_name = 'lat'
     1698       return_value = dom_def_var( vmea(l)%nc_filename,                        &
     1699                                   variable_name = variable_name,              &
     1700                                   dimension_names = (/'station'/),            &
     1701                                   output_type = 'real32' )
     1702!
     1703!--    longitude
     1704       variable_name = 'lon'
     1705       return_value = dom_def_var( vmea(l)%nc_filename,                        &
     1706                                   variable_name = variable_name,              &
     1707                                   dimension_names = (/'station'/),            &
     1708                                   output_type = 'real32' )
     1709!
     1710!--    Set attributes for the coordinate variables. Note, not all coordinates
     1711!--    have the same number of attributes.
     1712!--    Units
     1713       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1714                                   variable_name = 'time',                     &
     1715                                   attribute_name = char_unit,                 &
     1716                                   value = 'seconds since ' // origin_date_time )
     1717       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1718                                   variable_name = 'z',                        &
     1719                                   attribute_name = char_unit,                 &
     1720                                   value = 'm' )
     1721       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1722                                   variable_name = 'station_h',                &
     1723                                   attribute_name = char_unit,                 &
     1724                                   value = 'm' )
     1725       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1726                                   variable_name = 'x',                        &
     1727                                   attribute_name = char_unit,                 &
     1728                                   value = 'm' )
     1729       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1730                                   variable_name = 'y',                        &
     1731                                   attribute_name = char_unit,                 &
     1732                                   value = 'm' )
     1733       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1734                                   variable_name = 'E_UTM',                    &
     1735                                   attribute_name = char_unit,                 &
     1736                                   value = 'm' )
     1737       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1738                                   variable_name = 'N_UTM',                    &
     1739                                   attribute_name = char_unit,                 &
     1740                                   value = 'm' )
     1741       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1742                                   variable_name = 'lat',                      &
     1743                                   attribute_name = char_unit,                 &
     1744                                   value = 'degrees_north' )
     1745       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1746                                   variable_name = 'lon',                      &
     1747                                   attribute_name = char_unit,                 &
     1748                                   value = 'degrees_east' )
     1749!
     1750!--    Long name
     1751       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1752                                   variable_name = 'station_name',             &
     1753                                   attribute_name = char_long,                 &
     1754                                   value = 'station name')
     1755       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1756                                   variable_name = 'time',                     &
     1757                                   attribute_name = char_long,                 &
     1758                                   value = 'time')
     1759       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1760                                   variable_name = 'z',                        &
     1761                                   attribute_name = char_long,                 &
     1762                                   value = 'height above origin' )
     1763       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1764                                   variable_name = 'station_h',                &
     1765                                   attribute_name = char_long,                 &
     1766                                   value = 'surface altitude' )
     1767       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1768                                   variable_name = 'x',                        &
     1769                                   attribute_name = char_long,                 &
     1770                                   value = 'distance to origin in x-direction' )
     1771       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1772                                   variable_name = 'y',                        &
     1773                                   attribute_name = char_long,                 &
     1774                                   value = 'distance to origin in y-direction' )
     1775       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1776                                   variable_name = 'E_UTM',                    &
     1777                                   attribute_name = char_long,                 &
     1778                                   value = 'easting' )
     1779       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1780                                   variable_name = 'N_UTM',                    &
     1781                                   attribute_name = char_long,                 &
     1782                                   value = 'northing' )
     1783       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1784                                   variable_name = 'lat',                      &
     1785                                   attribute_name = char_long,                 &
     1786                                   value = 'latitude' )
     1787       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1788                                   variable_name = 'lon',                      &
     1789                                   attribute_name = char_long,                 &
     1790                                   value = 'longitude' )
     1791!
     1792!--    Standard name
     1793       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1794                                   variable_name = 'station_name',             &
     1795                                   attribute_name = char_standard,             &
     1796                                   value = 'platform_name')
     1797       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1798                                   variable_name = 'time',                     &
     1799                                   attribute_name = char_standard,             &
     1800                                   value = 'time')
     1801       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1802                                   variable_name = 'z',                        &
     1803                                   attribute_name = char_standard,             &
     1804                                   value = 'height_above_mean_sea_level' )
     1805       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1806                                   variable_name = 'station_h',                &
     1807                                   attribute_name = char_standard,             &
     1808                                   value = 'surface_altitude' )
     1809       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1810                                   variable_name = 'E_UTM',                    &
     1811                                   attribute_name = char_standard,             &
     1812                                   value = 'projection_x_coordinate' )
     1813       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1814                                   variable_name = 'N_UTM',                    &
     1815                                   attribute_name = char_standard,             &
     1816                                   value = 'projection_y_coordinate' )
     1817       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1818                                   variable_name = 'lat',                      &
     1819                                   attribute_name = char_standard,             &
     1820                                   value = 'latitude' )
     1821       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1822                                   variable_name = 'lon',                      &
     1823                                   attribute_name = char_standard,             &
     1824                                   value = 'longitude' )
     1825!
     1826!--    Axis
     1827       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1828                                   variable_name = 'time',                     &
     1829                                   attribute_name = 'axis',                    &
     1830                                   value = 'T')
     1831       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1832                                   variable_name = 'z',                        &
     1833                                   attribute_name = 'axis',                    &
     1834                                   value = 'Z' )
     1835       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1836                                   variable_name = 'x',                        &
     1837                                   attribute_name = 'axis',                    &
     1838                                   value = 'X' )
     1839       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1840                                   variable_name = 'y',                        &
     1841                                   attribute_name = 'axis',                    &
     1842                                   value = 'Y' )
     1843!
     1844!--    Set further individual attributes for the coordinate variables.
     1845!--    For station name
     1846       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1847                                   variable_name = 'station_name',             &
     1848                                   attribute_name = 'cf_role',                 &
     1849                                   value = 'timeseries_id' )
     1850!
     1851!--    For time
     1852       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1853                                   variable_name = 'time',                     &
     1854                                   attribute_name = 'calendar',                &
     1855                                   value = 'proleptic_gregorian' )
     1856       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1857                                   variable_name = 'time',                     &
     1858                                   attribute_name = 'bounds',                  &
     1859                                   value = 'time_bounds' )
     1860!
     1861!--    For vertical reference system
     1862       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1863                                   variable_name = 'vrs',                      &
     1864                                   attribute_name = char_long,                 &
     1865                                   value = 'vertical reference system' )
     1866       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1867                                   variable_name = 'vrs',                      &
     1868                                   attribute_name = 'system_name',             &
     1869                                   value = 'DHHN2016' )
     1870!
     1871!--    For z
     1872       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1873                                   variable_name = 'z',                        &
     1874                                   attribute_name = 'positive',                &
     1875                                   value = 'up' )
     1876!
     1877!--    For coordinate reference system
     1878       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1879                                   variable_name = 'crs',                      &
     1880                                   attribute_name = 'epsg_code',               &
     1881                                   value = coord_ref_sys%epsg_code )
     1882       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1883                                   variable_name = 'crs',                      &
     1884                                   attribute_name = 'false_easting',           &
     1885                                   value = coord_ref_sys%false_easting )
     1886       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1887                                   variable_name = 'crs',                      &
     1888                                   attribute_name = 'false_northing',          &
     1889                                   value = coord_ref_sys%false_northing )
     1890       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1891                                   variable_name = 'crs',                      &
     1892                                   attribute_name = 'grid_mapping_name',       &
     1893                                   value = coord_ref_sys%grid_mapping_name )
     1894       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1895                                   variable_name = 'crs',                      &
     1896                                   attribute_name = 'inverse_flattening',      &
     1897                                   value = coord_ref_sys%inverse_flattening )
     1898       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1899                                   variable_name = 'crs',                      &
     1900                                   attribute_name = 'latitude_of_projection_origin',&
     1901                                   value = coord_ref_sys%latitude_of_projection_origin )
     1902       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1903                                   variable_name = 'crs',                      &
     1904                                   attribute_name = char_long,                 &
     1905                                   value = coord_ref_sys%long_name )
     1906       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1907                                   variable_name = 'crs',                      &
     1908                                   attribute_name = 'longitude_of_central_meridian', &
     1909                                   value = coord_ref_sys%longitude_of_central_meridian )
     1910       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1911                                   variable_name = 'crs',                      &
     1912                                   attribute_name = 'longitude_of_prime_meridian', &
     1913                                   value = coord_ref_sys%longitude_of_prime_meridian )
     1914       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1915                                   variable_name = 'crs',                      &
     1916                                   attribute_name = 'scale_factor_at_central_meridian', &
     1917                                   value = coord_ref_sys%scale_factor_at_central_meridian )
     1918       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1919                                   variable_name = 'crs',                      &
     1920                                   attribute_name = 'semi_major_axis',         &
     1921                                   value = coord_ref_sys%semi_major_axis )
     1922       return_value = dom_def_att( vmea(l)%nc_filename,                        &
     1923                                   variable_name = 'crs',                      &
     1924                                   attribute_name = char_unit,                 &
     1925                                   value = coord_ref_sys%units )
     1926!
     1927!--    In case of sampled soil quantities, define further dimensions and
     1928!--    coordinates.
     1929       IF ( vmea(l)%soil_sampling )  THEN
     1930!
     1931!--       station for soil
     1932          ALLOCATE( ndim(1:vmea(l)%ns_soil_tot) )
     1933          DO  n = 1, vmea(l)%ns_soil_tot
     1934             ndim(n) = n
     1935          ENDDO
     1936
     1937          return_value = dom_def_dim( vmea(l)%nc_filename,                     &
     1938                                      dimension_name = 'station_soil',         &
     1939                                      output_type = 'int32',                   &
     1940                                      bounds = (/1_iwp,vmea(l)%ns_soil_tot/),  &
     1941                                      values_int32 = ndim )
     1942          DEALLOCATE( ndim )
     1943!
     1944!--       ntime for soil
     1945          ALLOCATE( ndim(1:ntimesteps) )
     1946          DO  n = 1, ntimesteps
     1947             ndim(n) = n
     1948          ENDDO
     1949
     1950          return_value = dom_def_dim( vmea(l)%nc_filename,                     &
     1951                                      dimension_name = 'ntime_soil',           &
     1952                                      output_type = 'int32',                   &
     1953                                      bounds = (/1_iwp,ntimesteps/),           &
     1954                                      values_int32 = ndim )
     1955          DEALLOCATE( ndim )
     1956!
     1957!--       time for soil
     1958          variable_name = 'time_soil'
     1959          return_value = dom_def_var( vmea(l)%nc_filename,                     &
     1960                                      variable_name = variable_name,           &
     1961                                      dimension_names = (/'station_soil',      &
     1962                                                          'ntime_soil  '/),    &
     1963                                      output_type = 'real32' )
     1964!
     1965!--       station_name for soil
     1966          variable_name = 'station_name_soil'
     1967          return_value = dom_def_var( vmea(l)%nc_filename,                     &
     1968                                      variable_name = variable_name,           &
     1969                                      dimension_names = (/ 'max_name_len',     &
     1970                                                           'station_soil'/),   &
     1971                                      output_type = 'char' )
     1972!
     1973!--       z
     1974          variable_name = 'z_soil'
     1975          return_value = dom_def_var( vmea(l)%nc_filename,                     &
     1976                                      variable_name = variable_name,           &
     1977                                      dimension_names = (/'station_soil'/),    &
     1978                                      output_type = 'real32' )
     1979!
     1980!--       station_h for soil
     1981          variable_name = 'station_h_soil'
     1982          return_value = dom_def_var( vmea(l)%nc_filename,                     &
     1983                                      variable_name = variable_name,           &
     1984                                      dimension_names = (/'station_soil'/),    &
     1985                                      output_type = 'real32' )
     1986!
     1987!--       x soil
     1988          variable_name = 'x_soil'
     1989          return_value = dom_def_var( vmea(l)%nc_filename,                     &
     1990                                      variable_name = variable_name,           &
     1991                                      dimension_names = (/'station_soil'/),    &
     1992                                      output_type = 'real32' )
     1993!
     1994!-        y soil
     1995          variable_name = 'y_soil'
     1996          return_value = dom_def_var( vmea(l)%nc_filename,                     &
     1997                                      variable_name = variable_name,           &
     1998                                      dimension_names = (/'station_soil'/),    &
     1999                                      output_type = 'real32' )
     2000!
     2001!--       E-UTM soil
     2002          variable_name = 'E_UTM_soil'
     2003          return_value = dom_def_var( vmea(l)%nc_filename,                     &
     2004                                      variable_name = variable_name,           &
     2005                                      dimension_names = (/'station_soil'/),    &
     2006                                      output_type = 'real32' )
     2007!
     2008!--       N-UTM soil
     2009          variable_name = 'N_UTM_soil'
     2010          return_value = dom_def_var( vmea(l)%nc_filename,                     &
     2011                                      variable_name = variable_name,           &
     2012                                      dimension_names = (/'station_soil'/),    &
     2013                                      output_type = 'real32' )
     2014!
     2015!--       latitude soil
     2016          variable_name = 'lat_soil'
     2017          return_value = dom_def_var( vmea(l)%nc_filename,                     &
     2018                                      variable_name = variable_name,           &
     2019                                      dimension_names = (/'station_soil'/),    &
     2020                                      output_type = 'real32' )
     2021!
     2022!--       longitude soil
     2023          variable_name = 'lon_soil'
     2024          return_value = dom_def_var( vmea(l)%nc_filename,                     &
     2025                                      variable_name = variable_name,           &
     2026                                      dimension_names = (/'station_soil'/),    &
     2027                                      output_type = 'real32' )
     2028!
     2029!--       Set attributes for the coordinate variables. Note, not all coordinates
     2030!--       have the same number of attributes.
     2031!--       Units
     2032          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2033                                      variable_name = 'time_soil',             &
     2034                                      attribute_name = char_unit,              &
     2035                                      value = 'seconds since ' // origin_date_time )
     2036          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2037                                      variable_name = 'z_soil',                &
     2038                                      attribute_name = char_unit,              &
     2039                                      value = 'm' )
     2040          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2041                                      variable_name = 'station_h_soil',        &
     2042                                      attribute_name = char_unit,              &
     2043                                      value = 'm' )
     2044          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2045                                      variable_name = 'x_soil',                &
     2046                                      attribute_name = char_unit,              &
     2047                                      value = 'm' )
     2048          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2049                                      variable_name = 'y_soil',                &
     2050                                      attribute_name = char_unit,              &
     2051                                      value = 'm' )
     2052          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2053                                      variable_name = 'E_UTM_soil',            &
     2054                                      attribute_name = char_unit,              &
     2055                                      value = 'm' )
     2056          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2057                                      variable_name = 'N_UTM_soil',            &
     2058                                      attribute_name = char_unit,              &
     2059                                      value = 'm' )
     2060          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2061                                      variable_name = 'lat_soil',              &
     2062                                      attribute_name = char_unit,              &
     2063                                      value = 'degrees_north' )
     2064          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2065                                      variable_name = 'lon_soil',              &
     2066                                      attribute_name = char_unit,              &
     2067                                      value = 'degrees_east' )
     2068!
     2069!--       Long name
     2070          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2071                                      variable_name = 'station_name_soil',     &
     2072                                      attribute_name = char_long,              &
     2073                                      value = 'station name')
     2074          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2075                                      variable_name = 'time_soil',             &
     2076                                      attribute_name = char_long,              &
     2077                                      value = 'time')
     2078          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2079                                      variable_name = 'z_soil',                &
     2080                                      attribute_name = char_long,              &
     2081                                      value = 'height above origin' )
     2082          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2083                                      variable_name = 'station_h_soil',        &
     2084                                      attribute_name = char_long,              &
     2085                                      value = 'surface altitude' )
     2086          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2087                                      variable_name = 'x_soil',                &
     2088                                      attribute_name = char_long,              &
     2089                                      value = 'distance to origin in x-direction' )
     2090          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2091                                      variable_name = 'y_soil',                &
     2092                                      attribute_name = char_long,              &
     2093                                      value = 'distance to origin in y-direction' )
     2094          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2095                                      variable_name = 'E_UTM_soil',            &
     2096                                      attribute_name = char_long,              &
     2097                                      value = 'easting' )
     2098          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2099                                      variable_name = 'N_UTM_soil',            &
     2100                                      attribute_name = char_long,              &
     2101                                      value = 'northing' )
     2102          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2103                                      variable_name = 'lat_soil',              &
     2104                                      attribute_name = char_long,              &
     2105                                      value = 'latitude' )
     2106          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2107                                      variable_name = 'lon_soil',              &
     2108                                      attribute_name = char_long,              &
     2109                                      value = 'longitude' )
     2110!
     2111!--       Standard name
     2112          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2113                                      variable_name = 'station_name_soil',     &
     2114                                      attribute_name = char_standard,          &
     2115                                      value = 'platform_name')
     2116          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2117                                      variable_name = 'time_soil',             &
     2118                                      attribute_name = char_standard,          &
     2119                                      value = 'time')
     2120          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2121                                      variable_name = 'z_soil',                &
     2122                                      attribute_name = char_standard,          &
     2123                                      value = 'height_above_mean_sea_level' )
     2124          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2125                                      variable_name = 'station_h_soil',        &
     2126                                      attribute_name = char_standard,          &
     2127                                      value = 'surface_altitude' )
     2128          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2129                                      variable_name = 'E_UTM_soil',            &
     2130                                      attribute_name = char_standard,          &
     2131                                      value = 'projection_x_coordinate' )
     2132          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2133                                      variable_name = 'N_UTM_soil',            &
     2134                                      attribute_name = char_standard,          &
     2135                                      value = 'projection_y_coordinate' )
     2136          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2137                                      variable_name = 'lat_soil',              &
     2138                                      attribute_name = char_standard,          &
     2139                                      value = 'latitude' )
     2140          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2141                                      variable_name = 'lon_soil',              &
     2142                                      attribute_name = char_standard,          &
     2143                                      value = 'longitude' )
     2144!
     2145!--       Axis
     2146          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2147                                      variable_name = 'time_soil',             &
     2148                                      attribute_name = 'axis',                 &
     2149                                      value = 'T')
     2150          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2151                                      variable_name = 'z_soil',                &
     2152                                      attribute_name = 'axis',                 &
     2153                                      value = 'Z' )
     2154          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2155                                      variable_name = 'x_soil',                &
     2156                                      attribute_name = 'axis',                 &
     2157                                      value = 'X' )
     2158          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2159                                      variable_name = 'y_soil',                &
     2160                                      attribute_name = 'axis',                 &
     2161                                      value = 'Y' )
     2162!
     2163!--       Set further individual attributes for the coordinate variables.
     2164!--       For station name soil
     2165          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2166                                      variable_name = 'station_name_soil',     &
     2167                                      attribute_name = 'cf_role',              &
     2168                                      value = 'timeseries_id' )
     2169!
     2170!--       For time soil
     2171          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2172                                      variable_name = 'time_soil',             &
     2173                                      attribute_name = 'calendar',             &
     2174                                      value = 'proleptic_gregorian' )
     2175          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2176                                      variable_name = 'time_soil',             &
     2177                                      attribute_name = 'bounds',               &
     2178                                      value = 'time_bounds' )
     2179!
     2180!--       For z soil
     2181          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2182                                      variable_name = 'z_soil',                &
     2183                                      attribute_name = 'positive',             &
     2184                                      value = 'up' )
     2185       ENDIF
     2186!
     2187!--    Define variables that shall be sampled.
     2188       DO  n = 1, vmea(l)%nmeas
     2189          variable_name = TRIM( vmea(l)%var_atts(n)%name )
     2190!
     2191!--       In order to link the correct dimension names, atmosphere and soil
     2192!--       variables need to be distinguished.
     2193          IF ( vmea(l)%soil_sampling  .AND.                                    &
     2194               ANY( TRIM( vmea(l)%var_atts(n)%name) == soil_vars ) )  THEN
     2195
     2196             return_value = dom_def_var( vmea(l)%nc_filename,                  &
     2197                                         variable_name = variable_name,        &
     2198                                         dimension_names = (/'station_soil',   &
     2199                                                             'ntime_soil  '/), &
     2200                                         output_type = 'real32' )
     2201          ELSE
     2202
     2203             return_value = dom_def_var( vmea(l)%nc_filename,                  &
     2204                                         variable_name = variable_name,        &
     2205                                         dimension_names = (/'station',        &
     2206                                                             'ntime  '/),      &
     2207                                         output_type = 'real32' )
     2208          ENDIF
     2209!
     2210!--       Set variable attributes. Please note, for some variables not all
     2211!--       attributes are defined, e.g. standard_name for the horizontal wind
     2212!--       components.
     2213          CALL vm_set_attributes( vmea(l)%var_atts(n) )
     2214
     2215          IF ( vmea(l)%var_atts(n)%long_name /= 'none' )  THEN
     2216             return_value = dom_def_att( vmea(l)%nc_filename,                  &
     2217                                         variable_name = variable_name,        &
     2218                                         attribute_name = char_long,           &
     2219                                         value = TRIM( vmea(l)%var_atts(n)%long_name ) )
     2220          ENDIF
     2221          IF ( vmea(l)%var_atts(n)%standard_name /= 'none' )  THEN
     2222             return_value = dom_def_att( vmea(l)%nc_filename,                  &
     2223                                         variable_name = variable_name,        &
     2224                                         attribute_name = char_standard,       &
     2225                                         value = TRIM( vmea(l)%var_atts(n)%standard_name ) )
     2226          ENDIF
     2227          IF ( vmea(l)%var_atts(n)%units /= 'none' )  THEN
     2228             return_value = dom_def_att( vmea(l)%nc_filename,                  &
     2229                                         variable_name = variable_name,        &
     2230                                         attribute_name = char_unit,           &
     2231                                         value = TRIM( vmea(l)%var_atts(n)%units ) )
     2232          ENDIF
     2233
     2234          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2235                                      variable_name = variable_name,           &
     2236                                      attribute_name = 'grid_mapping',         &
     2237                                      value = TRIM( vmea(l)%var_atts(n)%grid_mapping ) )
     2238
     2239          return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2240                                      variable_name = variable_name,           &
     2241                                      attribute_name = 'coordinates',          &
     2242                                      value = TRIM( vmea(l)%var_atts(n)%coordinates ) )
     2243
     2244!           return_value = dom_def_att( vmea(l)%nc_filename,                     &
     2245!                                       variable_name = variable_name,           &
     2246!                                       attribute_name = char_fill,              &
     2247!                                       value = vmea(l)%var_atts(n)%fill_value )
     2248
     2249       ENDDO  ! loop over variables per site
     2250
     2251    ENDDO  ! loop over sites
     2252
     2253
     2254 END SUBROUTINE vm_init_output
     2255
     2256!------------------------------------------------------------------------------!
     2257! Description:
     2258! ------------
     2259!> Parallel NetCDF output via data-output module.
     2260!------------------------------------------------------------------------------!
     2261 SUBROUTINE vm_data_output
     2262
     2263    CHARACTER(LEN=100) ::  variable_name !< name of output variable
     2264
     2265    INTEGER(iwp)       ::  l             !< loop index
     2266    INTEGER(iwp)       ::  n             !< loop index
     2267    INTEGER            ::  return_value  !< returned status value of called function
     2268
     2269    INTEGER(iwp)       ::  t_ind         !< time index
     2270
     2271    LOGICAL            ::  proceed_run   !< flag to indicate that the maximum number of
     2272                                         !< timesteps in the file will be exceeded
     2273
     2274    REAL(wp), DIMENSION(:), ALLOCATABLE           ::  oro_rel                  !< relative altitude of model surface
     2275    REAL(wp), DIMENSION(:), POINTER               ::  output_values_1d_pointer !< pointer for 1d output array
     2276    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET   ::  output_values_1d_target  !< target for 1d output array
     2277    REAL(wp), DIMENSION(:,:), POINTER             ::  output_values_2d_pointer !< pointer for 2d output array
     2278    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET ::  output_values_2d_target  !< target for 2d output array
     2279
     2280!
     2281!-- At the first call of this routine write the spatial coordinates.
     2282    IF ( .NOT. initial_write_coordinates )  THEN
     2283!
     2284!--    Write spatial coordinates.
     2285       DO  l = 1, vmea_general%nvm
     2286!
     2287!--       Skip if no observations were taken.
     2288          IF ( vmea(l)%ns_tot == 0  .AND.  vmea(l)%ns_soil_tot == 0 )  CYCLE
     2289
     2290          ALLOCATE( output_values_1d_target(vmea(l)%start_coord_a:vmea(l)%end_coord_a) )
     2291!
     2292!--       Output of Easting coordinate. Before output, recalculate EUTM.
     2293          output_values_1d_target = init_model%origin_x                        &
    10182294                    + REAL( vmea(l)%i(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dx &
    10192295                    * COS( init_model%rotation_angle * pi / 180.0_wp )         &
    10202296                    + REAL( vmea(l)%j(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dy &
    10212297                    * SIN( init_model%rotation_angle * pi / 180.0_wp )
    1022                  WRITE ( 27 )  'N_UTM                             '
    1023                  WRITE ( 27 )  init_model%origin_y                             &
     2298
     2299          output_values_1d_pointer => output_values_1d_target
     2300
     2301          return_value =                                                       &
     2302                  dom_write_var( vmea(l)%nc_filename,                          &
     2303                                 'E_UTM',                                      &
     2304                                 values_realwp_1d = output_values_1d_pointer,  &
     2305                                 bounds_start = (/vmea(l)%start_coord_a/),     &
     2306                                 bounds_end   = (/vmea(l)%end_coord_a  /) )
     2307!
     2308!--       Output of Northing coordinate. Before output, recalculate NUTM.
     2309          output_values_1d_target = init_model%origin_y                        &
    10242310                    - REAL( vmea(l)%i(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dx &
    10252311                    * SIN( init_model%rotation_angle * pi / 180.0_wp )         &
    10262312                    + REAL( vmea(l)%j(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dy &
    10272313                    * COS( init_model%rotation_angle * pi / 180.0_wp )
    1028                  WRITE ( 27 )  'Z_AG                              '
    1029                  WRITE ( 27 )  vmea(l)%z_ag(1:vmea(l)%ns)
    1030                  WRITE ( 27 )  'soil sampling                     '
    1031                  WRITE ( 27 )  MERGE( 'yes                               ',    &
    1032                                       'no                                ',    &
    1033                                       vmea(l)%soil_sampling )
    1034  
    1035                  IF ( vmea(l)%soil_sampling )  THEN                 
    1036                     WRITE ( 27 )  'total number of soil points       '                               
    1037                     WRITE ( 27 )  vmea(l)%ns_soil_tot
    1038                     WRITE ( 27 )  'number of soil points             '
    1039                     WRITE ( 27 )  vmea(l)%ns_soil
    1040                     WRITE ( 27 )  'E_UTM soil                        '
    1041                     WRITE ( 27 )  init_model%origin_x                          &
    1042                        + REAL( vmea(l)%i_soil(1:vmea(l)%ns_soil) + 0.5_wp,     &
    1043                                KIND = wp ) * dx                                &
    1044                        * COS( init_model%rotation_angle * pi / 180.0_wp )      &
    1045                        + REAL( vmea(l)%j_soil(1:vmea(l)%ns_soil) + 0.5_wp,     &
    1046                                KIND = wp ) * dy                                &
    1047                        * SIN( init_model%rotation_angle * pi / 180.0_wp )
    1048                     WRITE ( 27 )  'N_UTM soil                        '
    1049                     WRITE ( 27 )  init_model%origin_y                          &
    1050                        - REAL( vmea(l)%i_soil(1:vmea(l)%ns_soil) + 0.5_wp,     &
    1051                                KIND = wp ) * dx                                &
    1052                        * SIN( init_model%rotation_angle * pi / 180.0_wp )      &
    1053                        + REAL( vmea(l)%j_soil(1:vmea(l)%ns_soil) + 0.5_wp,     &
    1054                                KIND = wp ) * dy                                &
    1055                        * COS( init_model%rotation_angle * pi / 180.0_wp )
    1056                     WRITE ( 27 )  'DEPTH                             '
    1057                     WRITE ( 27 )  vmea(l)%depth(1:vmea(l)%ns_soil)
    1058                  ENDIF
    1059               ENDDO
    1060 
    1061            ENDIF
    1062         ENDDO
    1063        
     2314
     2315          output_values_1d_pointer => output_values_1d_target
     2316          return_value =                                                       &
     2317                  dom_write_var( vmea(l)%nc_filename,                          &
     2318                                 'N_UTM',                                      &
     2319                                 values_realwp_1d = output_values_1d_pointer,  &
     2320                                 bounds_start = (/vmea(l)%start_coord_a/),     &
     2321                                 bounds_end   = (/vmea(l)%end_coord_a  /) )
     2322!
     2323!--       Output of relative height coordinate.
     2324!--       Before this is output, first define the relative orographie height
     2325!--       and add this to z.
     2326          ALLOCATE( oro_rel(1:vmea(l)%ns) )
     2327          DO  n = 1, vmea(l)%ns
     2328             oro_rel(n) = zw(topo_top_ind(vmea(l)%j(n),vmea(l)%i(n),3))
     2329          ENDDO
     2330
     2331          output_values_1d_target = vmea(l)%zar(1:vmea(l)%ns) + oro_rel(:)
     2332          output_values_1d_pointer => output_values_1d_target
     2333          return_value =                                                       &
     2334                  dom_write_var( vmea(l)%nc_filename,                          &
     2335                                 'z',                                          &
     2336                                 values_realwp_1d = output_values_1d_pointer,  &
     2337                                 bounds_start = (/vmea(l)%start_coord_a/),     &
     2338                                 bounds_end   = (/vmea(l)%end_coord_a  /) )
     2339!
     2340!--       Write surface altitude for the station. Note, since z is already
     2341!--       a relative observation height, station_h must be zero, in order
     2342!--       to obtain the observation level.
     2343          output_values_1d_target = oro_rel(:)
     2344          output_values_1d_pointer => output_values_1d_target
     2345          return_value =                                                       &
     2346                  dom_write_var( vmea(l)%nc_filename,                          &
     2347                                 'station_h',                                  &
     2348                                 values_realwp_1d = output_values_1d_pointer,  &
     2349                                 bounds_start = (/vmea(l)%start_coord_a/),     &
     2350                                 bounds_end   = (/vmea(l)%end_coord_a  /) )
     2351
     2352          DEALLOCATE( oro_rel )
     2353          DEALLOCATE( output_values_1d_target )
     2354!
     2355!--       In case of sampled soil quantities, output also the respective
     2356!--       coordinate arrays.
     2357          IF ( vmea(l)%soil_sampling )  THEN
     2358             ALLOCATE( output_values_1d_target(vmea(l)%start_coord_s:vmea(l)%end_coord_s) )
     2359!
     2360!--          Output of Easting coordinate. Before output, recalculate EUTM.
     2361             output_values_1d_target = init_model%origin_x                     &
     2362               + REAL( vmea(l)%i(1:vmea(l)%ns_soil) + 0.5_wp, KIND = wp ) * dx &
     2363               * COS( init_model%rotation_angle * pi / 180.0_wp )              &
     2364               + REAL( vmea(l)%j(1:vmea(l)%ns_soil) + 0.5_wp, KIND = wp ) * dy &
     2365               * SIN( init_model%rotation_angle * pi / 180.0_wp )
     2366             output_values_1d_pointer => output_values_1d_target
     2367             return_value =                                                    &
     2368                  dom_write_var( vmea(l)%nc_filename,                          &
     2369                                 'E_UTM_soil',                                 &
     2370                                 values_realwp_1d = output_values_1d_pointer,  &
     2371                                 bounds_start = (/vmea(l)%start_coord_s/),     &
     2372                                 bounds_end   = (/vmea(l)%end_coord_s  /) )
     2373!
     2374!--          Output of Northing coordinate. Before output, recalculate NUTM.
     2375             output_values_1d_target = init_model%origin_y                     &
     2376               - REAL( vmea(l)%i(1:vmea(l)%ns_soil) + 0.5_wp, KIND = wp ) * dx &
     2377               * SIN( init_model%rotation_angle * pi / 180.0_wp )              &
     2378               + REAL( vmea(l)%j(1:vmea(l)%ns_soil) + 0.5_wp, KIND = wp ) * dy &
     2379               * COS( init_model%rotation_angle * pi / 180.0_wp )
     2380
     2381             output_values_1d_pointer => output_values_1d_target
     2382             return_value =                                                    &
     2383                  dom_write_var( vmea(l)%nc_filename,                          &
     2384                                 'N_UTM_soil',                                 &
     2385                                 values_realwp_1d = output_values_1d_pointer,  &
     2386                                 bounds_start = (/vmea(l)%start_coord_s/),     &
     2387                                 bounds_end   = (/vmea(l)%end_coord_s  /) )
     2388!
     2389!--          Output of relative height coordinate.
     2390!--          Before this is output, first define the relative orographie height
     2391!--          and add this to z.
     2392             ALLOCATE( oro_rel(1:vmea(l)%ns_soil) )
     2393             DO  n = 1, vmea(l)%ns
     2394                oro_rel(n) = zw(topo_top_ind(vmea(l)%j_soil(n),vmea(l)%i_soil(n),3))
     2395             ENDDO
     2396
     2397             output_values_1d_target = vmea(l)%depth(1:vmea(l)%ns_soil) + oro_rel(:)
     2398             output_values_1d_pointer => output_values_1d_target
     2399             return_value =                                                    &
     2400                  dom_write_var( vmea(l)%nc_filename,                          &
     2401                                 'z_soil',                                     &
     2402                                 values_realwp_1d = output_values_1d_pointer,  &
     2403                                 bounds_start = (/vmea(l)%start_coord_s/),     &
     2404                                 bounds_end   = (/vmea(l)%end_coord_s  /) )
     2405!
     2406!--          Write surface altitude for the station. Note, since z is already
     2407!--          a relative observation height, station_h must be zero, in order
     2408!--          to obtain the observation level.
     2409             output_values_1d_target = oro_rel(:)
     2410             output_values_1d_pointer => output_values_1d_target
     2411             return_value =                                                    &
     2412                  dom_write_var( vmea(l)%nc_filename,                          &
     2413                                 'station_h_soil',                             &
     2414                                 values_realwp_1d = output_values_1d_pointer,  &
     2415                                 bounds_start = (/vmea(l)%start_coord_s/),     &
     2416                                 bounds_end   = (/vmea(l)%end_coord_s  /) )
     2417
     2418             DEALLOCATE( oro_rel )
     2419             DEALLOCATE( output_values_1d_target )
     2420!
     2421!--          Write the stations name
     2422             
     2423          ENDIF
     2424
     2425       ENDDO  ! loop over sites
     2426
     2427       initial_write_coordinates = .TRUE.
     2428    ENDIF
     2429!
     2430!-- Loop over all sites.
     2431    DO  l = 1, vmea_general%nvm
     2432!
     2433!--    Skip if no observations were taken.
     2434       IF ( vmea(l)%ns_tot == 0  .AND.  vmea(l)%ns_soil_tot == 0 )  CYCLE
     2435!
     2436!--    Determine time index in file.
     2437       t_ind = vmea(l)%file_time_index + 1
     2438!
     2439!--    Write output variables. Distinguish between atmosphere and soil variables.
     2440       DO  n = 1, vmea(l)%nmeas
     2441          IF ( vmea(l)%soil_sampling  .AND.                                    &
     2442            ANY( TRIM( vmea(l)%var_atts(n)%name) == soil_vars ) )  THEN
     2443!
     2444!--          Write time coordinate to file
     2445             variable_name = 'time_soil'
     2446             ALLOCATE( output_values_2d_target(t_ind:t_ind,vmea(l)%start_coord_s:vmea(l)%end_coord_s) )
     2447             output_values_2d_target(t_ind,:) = time_since_reference_point
     2448             output_values_2d_pointer => output_values_2d_target
     2449
     2450             return_value = dom_write_var( vmea(l)%nc_filename,                &
     2451                                           variable_name,                      &
     2452                                           values_realwp_2d = output_values_2d_pointer, &
     2453                                           bounds_start = (/vmea(l)%start_coord_s, t_ind/), &
     2454                                           bounds_end   = (/vmea(l)%end_coord_s, t_ind /) )
     2455
     2456             variable_name = TRIM( vmea(l)%var_atts(n)%name )
     2457             output_values_2d_target(t_ind,:) = vmea(l)%measured_vars_soil(:,n)
     2458             output_values_2d_pointer => output_values_2d_target
     2459             return_value =                                                    &
     2460                      dom_write_var( vmea(l)%nc_filename,                      &
     2461                                     variable_name,                            &
     2462                                     values_realwp_2d = output_values_2d_pointer, &
     2463                                     bounds_start = (/vmea(l)%start_coord_s, t_ind/), &
     2464                                     bounds_end   = (/vmea(l)%end_coord_s, t_ind  /) )
     2465             DEALLOCATE( output_values_2d_target )
     2466          ELSE
     2467!
     2468!--          Write time coordinate to file
     2469             variable_name = 'time'
     2470             ALLOCATE( output_values_2d_target(t_ind:t_ind,vmea(l)%start_coord_a:vmea(l)%end_coord_a) )
     2471             output_values_2d_target(t_ind,:) = time_since_reference_point
     2472             output_values_2d_pointer => output_values_2d_target
     2473
     2474             return_value = dom_write_var( vmea(l)%nc_filename,                &
     2475                                           variable_name,                      &
     2476                                           values_realwp_2d = output_values_2d_pointer, &
     2477                                           bounds_start = (/vmea(l)%start_coord_a, t_ind/), &
     2478                                           bounds_end   = (/vmea(l)%end_coord_a, t_ind/) )
     2479
     2480             variable_name = TRIM( vmea(l)%var_atts(n)%name )
     2481
     2482             output_values_2d_target(t_ind,:) = vmea(l)%measured_vars(:,n)
     2483             output_values_2d_pointer => output_values_2d_target
     2484             return_value =                                                    &
     2485                      dom_write_var( vmea(l)%nc_filename,                      &
     2486                                     variable_name,                            &
     2487                                     values_realwp_2d = output_values_2d_pointer, &
     2488                                     bounds_start = (/ vmea(l)%start_coord_a, t_ind /), &
     2489                                     bounds_end   = (/ vmea(l)%end_coord_a, t_ind /) )
     2490
     2491             DEALLOCATE( output_values_2d_target )
     2492          ENDIF
     2493       ENDDO
     2494!
     2495!--    Update number of written time indices
     2496       vmea(l)%file_time_index = t_ind
     2497
     2498    ENDDO  ! loop over sites
     2499!
     2500!-- Check if the time index exceeds its maximum values. This case, terminate the
     2501!-- run and force a restart. This is simply done by creating the file
     2502!-- DO_RESTART_NOW.
     2503    proceed_run = .TRUE.
     2504    IF ( ANY( vmea(:)%file_time_index >= ntimesteps - 1) )  proceed_run = .FALSE.
     2505
    10642506#if defined( __parallel )
    1065         CALL MPI_BARRIER( comm2d, ierr )
     2507    CALL MPI_ALLREDUCE( MPI_IN_PLACE,                                          &
     2508                        proceed_run,                                           &
     2509                        1,                                                     &
     2510                        MPI_LOGICAL,                                           &
     2511                        MPI_LAND,                                              &
     2512                        MPI_COMM_WORLD,                                        &
     2513                        ierr )
    10662514#endif
    10672515!
    1068 !--     After header information is written, set control flag to .FALSE.
    1069         init = .FALSE.
    1070 !
    1071 !--  Data output at each measurement timestep on each PE
    1072      ELSE
    1073         DO  i = 0, io_blocks-1
    1074 
    1075            IF ( i == io_group )  THEN
    1076               WRITE( 27 )  'output time                       '
    1077               WRITE( 27 )  time_since_reference_point
    1078               DO  l = 1, vmea_general%nvm
    1079 !
    1080 !--              Skip binary writing if no observation points are defined on PE
    1081                  IF ( vmea(l)%ns < 1  .AND.  vmea(l)%ns_soil < 1)  CYCLE                 
    1082                  DO  n = 1, vmea(l)%nmeas
    1083                     WRITE( 27 )  vmea(l)%measured_vars_name(n)
    1084                     IF ( vmea(l)%soil_sampling  .AND.                           &
    1085                          ANY( TRIM( vmea(l)%measured_vars_name(n))  ==          &
    1086                               soil_vars ) )  THEN                   
    1087                        WRITE( 27 )  vmea(l)%measured_vars_soil(:,n)
    1088                     ELSE
    1089                        WRITE( 27 )  vmea(l)%measured_vars(:,n)
    1090                     ENDIF
    1091                  ENDDO
    1092            
    1093               ENDDO
    1094            ENDIF
    1095         ENDDO
    1096 #if defined( __parallel )
    1097         CALL MPI_BARRIER( comm2d, ierr )
    1098 #endif
    1099      ENDIF
    1100  
    1101   END SUBROUTINE vm_data_output 
    1102  
    1103  
    1104 !------------------------------------------------------------------------------!
    1105 ! Description:
    1106 ! ------------
    1107 !> Write end-of-file statement as last action.
    1108 !------------------------------------------------------------------------------!
    1109   SUBROUTINE vm_last_actions
    1110    
    1111      USE pegrid
    1112    
    1113      IMPLICIT NONE
    1114          
    1115      INTEGER(iwp) ::  i         !< running index over IO blocks   
    1116  
    1117      DO  i = 0, io_blocks-1
    1118         IF ( i == io_group )  THEN
    1119            WRITE( 27 )  'EOF                               '
    1120         ENDIF
    1121      ENDDO
    1122 #if defined( __parallel )
    1123         CALL MPI_BARRIER( comm2d, ierr )
    1124 #endif
    1125 !
    1126 !--  Close binary file
    1127      CALL close_file( 27 )
    1128  
    1129   END SUBROUTINE vm_last_actions
    1130  
     2516!-- Create file DO_RESTART_NOW to force a restart. This file is only created
     2517!-- by rank 0 of the root domain.
     2518    IF ( myid == 0  .AND.  .NOT. child_domain  .AND.  .NOT. proceed_run )  THEN
     2519       OPEN( 999, FILE='DO_RESTART_NOW' )
     2520       CLOSE( 999 )
     2521
     2522       message_string = 'Run will be terminated because virtual ' //           &
     2523                        'measurement times index exceeds its maximum value ' //&
     2524                        'in the output files. A restart run is forced.'
     2525       CALL message( 'vm_data_output', 'PA0707', 0, 0, 0, 6, 0 )
     2526    ENDIF
     2527
     2528  END SUBROUTINE vm_data_output
     2529
    11312530!------------------------------------------------------------------------------!
    11322531! Description:
    … …  
    11362535  SUBROUTINE vm_sampling
    11372536
    1138     USE arrays_3d,                                                             &
    1139         ONLY:  exner, pt, q, u, v, w
    1140 
    11412537    USE radiation_model_mod,                                                   &
    1142         ONLY:  radiation 
     2538        ONLY:  radiation
    11432539
    11442540    USE surface_mod,                                                           &
    1145         ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
    1146    
    1147      IMPLICIT NONE
    1148      
     2541        ONLY:  surf_def_h,                                                     &
     2542               surf_lsm_h,                                                     &
     2543               surf_usm_h
     2544
    11492545     INTEGER(iwp) ::  i         !< grid index in x-direction
    11502546     INTEGER(iwp) ::  j         !< grid index in y-direction
    … …  
    11562552     INTEGER(iwp) ::  n         !< running index over all measured variables at a station
    11572553     INTEGER(iwp) ::  nn        !< running index over the number of chemcal species
    1158      
     2554
    11592555     LOGICAL ::  match_lsm !< flag indicating natural-type surface
    11602556     LOGICAL ::  match_usm !< flag indicating urban-type surface
    1161 !
    1162 !--  Loop over all sites.
     2557
     2558     REAL(wp) ::  e_s      !< saturation water vapor pressure
     2559     REAL(wp) ::  q_s      !< saturation mixing ratio
     2560     REAL(wp) ::  q_wv     !< mixing ratio
     2561!
     2562!--  Loop over all sites.
    11632563     DO  l = 1, vmea_general%nvm
    11642564!
    11652565!--     At the beginning, set _FillValues
    11662566        IF ( ALLOCATED( vmea(l)%measured_vars      ) )                         &
    1167            vmea(l)%measured_vars      = vmea(l)%fillout 
     2567           vmea(l)%measured_vars      = vmea(l)%fillout
    11682568        IF ( ALLOCATED( vmea(l)%measured_vars_soil ) )                         &
    1169            vmea(l)%measured_vars_soil = vmea(l)%fillout 
    1170 !
    1171 !--     Loop over all variables measured at this site. 
     2569           vmea(l)%measured_vars_soil = vmea(l)%fillout
     2570!
     2571!--     Loop over all variables measured at this site.
    11722572        DO  n = 1, vmea(l)%nmeas
    1173        
    1174            SELECT CASE ( TRIM( vmea(l)%measured_vars_name(n) ) )
    1175            
    1176               CASE ( 'theta' )
     2573
     2574           SELECT CASE ( TRIM( vmea(l)%var_atts(n)%name ) )
     2575
     2576              CASE ( 'theta' ) ! potential temperature
    11772577                 IF ( .NOT. neutral )  THEN
    11782578                    DO  m = 1, vmea(l)%ns
    … …  
    11832583                    ENDDO
    11842584                 ENDIF
    1185                  
    1186               CASE ( 'ta' )
     2585
     2586              CASE ( 'ta' ) ! absolute temperature
    11872587                 IF ( .NOT. neutral )  THEN
    11882588                    DO  m = 1, vmea(l)%ns
    … …  
    11902590                       j = vmea(l)%j(m)
    11912591                       i = vmea(l)%i(m)
    1192                        vmea(l)%measured_vars(m,n) = pt(k,j,i) * exner( k )
     2592                       vmea(l)%measured_vars(m,n) = pt(k,j,i) * exner( k )     &
     2593                                                  - degc_to_k
    11932594                    ENDDO
    11942595                 ENDIF
    1195              
     2596
    11962597              CASE ( 't_va' )
    1197                  
    1198               CASE ( 'hus', 'haa' )
     2598
     2599              CASE ( 'hus' ) ! mixing ratio
    11992600                 IF ( humidity )  THEN
    12002601                    DO  m = 1, vmea(l)%ns
    … …  
    12052606                    ENDDO
    12062607                 ENDIF
    1207                  
    1208               CASE ( 'u', 'ua' )
     2608
     2609              CASE ( 'haa' ) ! absolute humidity
     2610                 IF ( humidity )  THEN
     2611                    DO  m = 1, vmea(l)%ns
     2612                       k = vmea(l)%k(m)
     2613                       j = vmea(l)%j(m)
     2614                       i = vmea(l)%i(m)
     2615                       vmea(l)%measured_vars(m,n) = ( q(k,j,i)                 &
     2616                                                    / ( 1.0_wp - q(k,j,i) ) )  &
     2617                                                  * rho_air(k)
     2618                    ENDDO
     2619                 ENDIF
     2620
     2621              CASE ( 'pwv' ) ! water vapor partial pressure
     2622                 IF ( humidity )  THEN
     2623!                     DO  m = 1, vmea(l)%ns
     2624!                        k = vmea(l)%k(m)
     2625!                        j = vmea(l)%j(m)
     2626!                        i = vmea(l)%i(m)
     2627!                        vmea(l)%measured_vars(m,n) = ( q(k,j,i)                 &
     2628!                                                     / ( 1.0_wp - q(k,j,i) ) )  &
     2629!                                                   * rho_air(k)
     2630!                     ENDDO
     2631                 ENDIF
     2632
     2633              CASE ( 'hur' ) ! relative humidity
     2634                 IF ( humidity )  THEN
     2635                    DO  m = 1, vmea(l)%ns
     2636                       k = vmea(l)%k(m)
     2637                       j = vmea(l)%j(m)
     2638                       i = vmea(l)%i(m)
     2639!
     2640!--                    Calculate actual temperature, water vapor saturation
     2641!--                    pressure, and based on this  the saturation mixing ratio.
     2642                       e_s  = magnus( exner(k) * pt(k,j,i) )
     2643                       q_s  = rd_d_rv * e_s / ( hyp(k) - e_s )
     2644                       q_wv = ( q(k,j,i) / ( 1.0_wp - q(k,j,i) ) ) * rho_air(k)
     2645
     2646                       vmea(l)%measured_vars(m,n) = q_wv / ( q_s + 1E-10_wp )
     2647                    ENDDO
     2648                 ENDIF
     2649
     2650              CASE ( 'u', 'ua' ) ! u-component
    12092651                 DO  m = 1, vmea(l)%ns
    12102652                    k = vmea(l)%k(m)
    … …  
    12132655                    vmea(l)%measured_vars(m,n) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
    12142656                 ENDDO
    1215                  
    1216               CASE ( 'v', 'va' )
     2657
     2658              CASE ( 'v', 'va' ) ! v-component
    12172659                 DO  m = 1, vmea(l)%ns
    12182660                    k = vmea(l)%k(m)
    … …  
    12212663                    vmea(l)%measured_vars(m,n) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
    12222664                 ENDDO
    1223                  
    1224               CASE ( 'w' )
     2665
     2666              CASE ( 'w' ) ! w-component
    12252667                 DO  m = 1, vmea(l)%ns
    1226                     k = MAX ( 1, vmea(l)%k(m) ) 
     2668                    k = MAX ( 1, vmea(l)%k(m) )
    12272669                    j = vmea(l)%j(m)
    12282670                    i = vmea(l)%i(m)
    12292671                    vmea(l)%measured_vars(m,n) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
    12302672                 ENDDO
    1231                  
    1232               CASE ( 'wspeed' )
     2673
     2674              CASE ( 'wspeed' ) ! horizontal wind speed
    12332675                 DO  m = 1, vmea(l)%ns
    12342676                    k = vmea(l)%k(m)
    … …  
    12402682                                                     )
    12412683                 ENDDO
    1242                  
    1243               CASE ( 'wdir' )
     2684
     2685              CASE ( 'wdir' ) ! wind direction
    12442686                 DO  m = 1, vmea(l)%ns
    12452687                    k = vmea(l)%k(m)
    12462688                    j = vmea(l)%j(m)
    12472689                    i = vmea(l)%i(m)
    1248                    
    1249                     vmea(l)%measured_vars(m,n) = ATAN2(                        &
    1250                                        - 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ),   &
    1251                                        - 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )    &
    1252                                                       ) * 180.0_wp / pi
     2690
     2691                    vmea(l)%measured_vars(m,n) = 180.0_wp + 180.0_wp / pi      &
     2692                                               * ATAN2(                        &
     2693                                          0.5_wp * ( v(k,j,i) + v(k,j+1,i) ),  &
     2694                                          0.5_wp * ( u(k,j,i) + u(k,j,i+1) )   &
     2695                                                      )
    12532696                 ENDDO
    1254                  
     2697
    12552698              CASE ( 'utheta' )
    12562699                 DO  m = 1, vmea(l)%ns
    … …  
    12622705                                                   pt(k,j,i)
    12632706                 ENDDO
    1264                  
     2707
    12652708              CASE ( 'vtheta' )
    12662709                 DO  m = 1, vmea(l)%ns
    … …  
    12722715                                                   pt(k,j,i)
    12732716                 ENDDO
    1274                  
     2717
    12752718              CASE ( 'wtheta' )
    12762719                 DO  m = 1, vmea(l)%ns
    … …  
    12822725                                                   pt(k,j,i)
    12832726                 ENDDO
    1284                  
     2727
     2728              CASE ( 'uqv' )
     2729                 IF ( humidity )  THEN
     2730                    DO  m = 1, vmea(l)%ns
     2731                       k = vmea(l)%k(m)
     2732                       j = vmea(l)%j(m)
     2733                       i = vmea(l)%i(m)
     2734                       vmea(l)%measured_vars(m,n) = 0.5_wp *                   &
     2735                                                    ( u(k,j,i) + u(k,j,i+1) ) *&
     2736                                                      q(k,j,i)
     2737                    ENDDO
     2738                 ENDIF
     2739
     2740              CASE ( 'vqv' )
     2741                 IF ( humidity )  THEN
     2742                    DO  m = 1, vmea(l)%ns
     2743                       k = vmea(l)%k(m)
     2744                       j = vmea(l)%j(m)
     2745                       i = vmea(l)%i(m)
     2746                       vmea(l)%measured_vars(m,n) = 0.5_wp *                   &
     2747                                                    ( v(k,j,i) + v(k,j+1,i) ) *&
     2748                                                      q(k,j,i)
     2749                    ENDDO
     2750                 ENDIF
     2751
     2752              CASE ( 'wqv' )
     2753                 IF ( humidity )  THEN
     2754                    DO  m = 1, vmea(l)%ns
     2755                       k = MAX ( 1, vmea(l)%k(m) )
     2756                       j = vmea(l)%j(m)
     2757                       i = vmea(l)%i(m)
     2758                       vmea(l)%measured_vars(m,n) = 0.5_wp *                   &
     2759                                                    ( w(k-1,j,i) + w(k,j,i) ) *&
     2760                                                      q(k,j,i)
     2761                    ENDDO
     2762                 ENDIF
     2763
    12852764              CASE ( 'uw' )
    12862765                 DO  m = 1, vmea(l)%ns
    … …  
    12922771                                                 ( u(k,j,i)   + u(k,j,i+1) )
    12932772                 ENDDO
    1294                  
     2773
    12952774              CASE ( 'vw' )
    12962775                 DO  m = 1, vmea(l)%ns
    … …  
    13022781                                                 ( v(k,j,i)   + v(k,j+1,i) )
    13032782                 ENDDO
    1304                  
     2783
    13052784              CASE ( 'uv' )
    13062785                 DO  m = 1, vmea(l)%ns
    1307                     k = MAX ( 1, vmea(l)%k(m) )
     2786                    k = vmea(l)%k(m)
    13082787                    j = vmea(l)%j(m)
    13092788                    i = vmea(l)%i(m)
    … …  
    13132792                 ENDDO
    13142793!
    1315 !--           List of variables may need extension.
    1316               CASE ( 'mcpm1', 'mcpm2p5', 'mcpm10', 'mfco', 'mfno', 'mfno2',    &
    1317                      'tro3' )                     
     2794!--           Chemistry variables. List of variables may need extension.
     2795!--           Note, gas species in PALM are in ppm and no distinction is made
     2796!--           between mole-fraction and concentration quantities (all are
     2797!--           output in ppm so far).
     2798              CASE ( 'mcpm1', 'mcpm2p5', 'mcpm10', 'mfno', 'mfno2',            &
     2799                     'mcno', 'mcno2', 'tro3' )
    13182800                 IF ( air_chemistry )  THEN
    13192801!
    1320 !--                 First, search for the measured variable in the chem_vars 
    1321 !--                 list, in order to get the internal name of the variable. 
     2802!--                 First, search for the measured variable in the chem_vars
     2803!--                 list, in order to get the internal name of the variable.
    13222804                    DO  nn = 1, UBOUND( chem_vars, 2 )
    1323                        IF ( TRIM( vmea(l)%measured_vars_name(m) ) ==           &
     2805                       IF ( TRIM( vmea(l)%var_atts(n)%name ) ==                &
    13242806                            TRIM( chem_vars(0,nn) ) )  ind_chem = nn
    13252807                    ENDDO
    13262808!
    1327 !--                 Run loop over all chemical species, if the measured 
     2809!--                 Run loop over all chemical species, if the measured
    13282810!--                 variable matches the interal name, sample the variable.
    1329                     DO  nn = 1, nvar                   
     2811!--                 Note, nvar as a chemistry-module variable.
     2812                    DO  nn = 1, nvar
    13302813                       IF ( TRIM( chem_vars(1,ind_chem) ) ==                   &
    1331                             TRIM( chem_species(nn)%name ) )  THEN                           
    1332                           DO  m = 1, vmea(l)%ns             
     2814                            TRIM( chem_species(nn)%name ) )  THEN
     2815                          DO  m = 1, vmea(l)%ns
    13332816                             k = vmea(l)%k(m)
    13342817                             j = vmea(l)%j(m)
    1335                              i = vmea(l)%i(m)                   
     2818                             i = vmea(l)%i(m)
    13362819                             vmea(l)%measured_vars(m,n) =                      &
    13372820                                                   chem_species(nn)%conc(k,j,i)
    … …  
    13402823                    ENDDO
    13412824                 ENDIF
    1342                  
    1343               CASE ( 'us' )
     2825
     2826              CASE ( 'us' ) ! friction velocity
    13442827                 DO  m = 1, vmea(l)%ns
    13452828!
    1346 !--                 Surface data is only available on inner subdomains, not 
    1347 !--                 on ghost points. Hence, limit the indices. 
    1348                     j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
    1349                     j = MERGE( j           , nyn, j            > nyn )
    1350                     i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
    1351                     i = MERGE( i           , nxr, i            > nxr )
    1352                    
     2829!--                 Surface data is only available on inner subdomains, not
     2830!--                 on ghost points. Hence, limit the indices.
     2831                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
     2832                    j = MERGE( j           , nyn, j            < nyn )
     2833                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
     2834                    i = MERGE( i           , nxr, i            < nxr )
     2835
    13532836                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
    13542837                             surf_def_h(0)%end_index(j,i)
    … …  
    13642847                    ENDDO
    13652848                 ENDDO
    1366                  
    1367               CASE ( 'ts' )
     2849
     2850              CASE ( 'thetas' ) ! scaling parameter temperature
    13682851                 DO  m = 1, vmea(l)%ns
    13692852!
    1370 !--                 Surface data is only available on inner subdomains, not 
    1371 !--                 on ghost points. Hence, limit the indices. 
    1372                     j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
    1373                     j = MERGE( j           , nyn, j            > nyn )
    1374                     i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
    1375                     i = MERGE( i           , nxr, i            > nxr )
    1376                    
     2853!--                 Surface data is only available on inner subdomains, not
     2854!--                 on ghost points. Hence, limit the indices.
     2855                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
     2856                    j = MERGE( j           , nyn, j            < nyn )
     2857                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
     2858                    i = MERGE( i           , nxr, i            < nxr )
     2859
    13772860                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
    13782861                             surf_def_h(0)%end_index(j,i)
    … …  
    13882871                    ENDDO
    13892872                 ENDDO
    1390                  
    1391               CASE ( 'hfls' )
     2873
     2874              CASE ( 'hfls' ) ! surface latent heat flux
    13922875                 DO  m = 1, vmea(l)%ns
    13932876!
    1394 !--                 Surface data is only available on inner subdomains, not 
    1395 !--                 on ghost points. Hence, limit the indices. 
    1396                     j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
    1397                     j = MERGE( j           , nyn, j            > nyn )
    1398                     i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
    1399                     i = MERGE( i           , nxr, i            > nxr )
    1400                    
     2877!--                 Surface data is only available on inner subdomains, not
     2878!--                 on ghost points. Hence, limit the indices.
     2879                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
     2880                    j = MERGE( j           , nyn, j            < nyn )
     2881                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
     2882                    i = MERGE( i           , nxr, i            < nxr )
     2883
    14012884                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
    14022885                             surf_def_h(0)%end_index(j,i)
    … …  
    14122895                    ENDDO
    14132896                 ENDDO
    1414                  
    1415               CASE ( 'hfss' )
     2897
     2898              CASE ( 'hfss' ) ! surface sensible heat flux
    14162899                 DO  m = 1, vmea(l)%ns
    14172900!
    1418 !--                 Surface data is only available on inner subdomains, not 
    1419 !--                 on ghost points. Hence, limit the indices. 
    1420                     j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
    1421                     j = MERGE( j           , nyn, j            > nyn )
    1422                     i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
    1423                     i = MERGE( i           , nxr, i            > nxr )
    1424                    
     2901!--                 Surface data is only available on inner subdomains, not
     2902!--                 on ghost points. Hence, limit the indices.
     2903                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
     2904                    j = MERGE( j           , nyn, j            < nyn )
     2905                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
     2906                    i = MERGE( i           , nxr, i            < nxr )
     2907
    14252908                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
    14262909                             surf_def_h(0)%end_index(j,i)
    … …  
    14362919                    ENDDO
    14372920                 ENDDO
    1438                  
    1439               CASE ( 'rnds' )
     2921
     2922              CASE ( 'hfdg' ) ! ground heat flux
     2923                 DO  m = 1, vmea(l)%ns
     2924!
     2925!--                 Surface data is only available on inner subdomains, not
     2926!--                 on ghost points. Hence, limit the indices.
     2927                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
     2928                    j = MERGE( j           , nyn, j            < nyn )
     2929                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
     2930                    i = MERGE( i           , nxr, i            < nxr )
     2931
     2932                    DO  mm = surf_lsm_h%start_index(j,i),                      &
     2933                             surf_lsm_h%end_index(j,i)
     2934                       vmea(l)%measured_vars(m,n) = surf_lsm_h%ghf(mm)
     2935                    ENDDO
     2936                 ENDDO
     2937
     2938              CASE ( 'lwcs' )  ! liquid water of soil layer
     2939!                  DO  m = 1, vmea(l)%ns
     2940! !
     2941! !--                 Surface data is only available on inner subdomains, not
     2942! !--                 on ghost points. Hence, limit the indices.
     2943!                     j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
     2944!                     j = MERGE( j           , nyn, j            < nyn )
     2945!                     i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
     2946!                     i = MERGE( i           , nxr, i            < nxr )
     2947!
     2948!                     DO  mm = surf_lsm_h%start_index(j,i),                      &
     2949!                              surf_lsm_h%end_index(j,i)
     2950!                        vmea(l)%measured_vars(m,n) = ?
     2951!                     ENDDO
     2952!                  ENDDO
     2953
     2954              CASE ( 'rnds' ) ! surface net radiation
    14402955                 IF ( radiation )  THEN
    14412956                    DO  m = 1, vmea(l)%ns
    14422957!
    1443 !--                    Surface data is only available on inner subdomains, not 
    1444 !--                    on ghost points. Hence, limit the indices. 
    1445                        j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
    1446                        j = MERGE( j           , nyn, j            > nyn )
    1447                        i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
    1448                        i = MERGE( i           , nxr, i            > nxr )
    1449                    
     2958!--                    Surface data is only available on inner subdomains, not
     2959!--                    on ghost points. Hence, limit the indices.
     2960                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
     2961                       j = MERGE( j           , nyn, j            < nyn )
     2962                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
     2963                       i = MERGE( i           , nxr, i            < nxr )
     2964
    14502965                       DO  mm = surf_lsm_h%start_index(j,i),                   &
    14512966                                surf_lsm_h%end_index(j,i)
    … …  
    14582973                    ENDDO
    14592974                 ENDIF
    1460                  
    1461               CASE ( 'rsus' )
     2975
     2976              CASE ( 'rsus' ) ! surface shortwave out
    14622977                 IF ( radiation )  THEN
    14632978                    DO  m = 1, vmea(l)%ns
    14642979!
    1465 !--                    Surface data is only available on inner subdomains, not 
    1466 !--                    on ghost points. Hence, limit the indices. 
    1467                        j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
    1468                        j = MERGE( j           , nyn, j            > nyn )
    1469                        i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
    1470                        i = MERGE( i           , nxr, i            > nxr )
    1471                    
     2980!--                    Surface data is only available on inner subdomains, not
     2981!--                    on ghost points. Hence, limit the indices.
     2982                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
     2983                       j = MERGE( j           , nyn, j            < nyn )
     2984                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
     2985                       i = MERGE( i           , nxr, i            < nxr )
     2986
    14722987                       DO  mm = surf_lsm_h%start_index(j,i),                   &
    14732988                                surf_lsm_h%end_index(j,i)
    … …  
    14802995                    ENDDO
    14812996                 ENDIF
    1482                  
    1483               CASE ( 'rsds' )
     2997
     2998              CASE ( 'rsds' ) ! surface shortwave in
    14842999                 IF ( radiation )  THEN
    14853000                    DO  m = 1, vmea(l)%ns
    14863001!
    1487 !--                    Surface data is only available on inner subdomains, not 
    1488 !--                    on ghost points. Hence, limit the indices. 
    1489                        j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
    1490                        j = MERGE( j           , nyn, j            > nyn )
    1491                        i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
    1492                        i = MERGE( i           , nxr, i            > nxr )
    1493                    
     3002!--                    Surface data is only available on inner subdomains, not
     3003!--                    on ghost points. Hence, limit the indices.
     3004                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
     3005                       j = MERGE( j           , nyn, j            < nyn )
     3006                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
     3007                       i = MERGE( i           , nxr, i            < nxr )
     3008
    14943009                       DO  mm = surf_lsm_h%start_index(j,i),                   &
    14953010                                surf_lsm_h%end_index(j,i)
    … …  
    15023017                    ENDDO
    15033018                 ENDIF
    1504                  
    1505               CASE ( 'rlus' )
     3019
     3020              CASE ( 'rlus' ) ! surface longwave out
    15063021                 IF ( radiation )  THEN
    15073022                    DO  m = 1, vmea(l)%ns
    15083023!
    1509 !--                    Surface data is only available on inner subdomains, not 
    1510 !--                    on ghost points. Hence, limit the indices. 
    1511                        j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
    1512                        j = MERGE( j           , nyn, j            > nyn )
    1513                        i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
    1514                        i = MERGE( i           , nxr, i            > nxr )
    1515                    
     3024!--                    Surface data is only available on inner subdomains, not
     3025!--                    on ghost points. Hence, limit the indices.
     3026                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
     3027                       j = MERGE( j           , nyn, j            < nyn )
     3028                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
     3029                       i = MERGE( i           , nxr, i            < nxr )
     3030
    15163031                       DO  mm = surf_lsm_h%start_index(j,i),                   &
    15173032                                surf_lsm_h%end_index(j,i)
    … …  
    15243039                    ENDDO
    15253040                 ENDIF
    1526                  
    1527               CASE ( 'rlds' )
     3041
     3042              CASE ( 'rlds' ) ! surface longwave in
    15283043                 IF ( radiation )  THEN
    15293044                    DO  m = 1, vmea(l)%ns
    15303045!
    1531 !--                    Surface data is only available on inner subdomains, not 
    1532 !--                    on ghost points. Hence, limit the indices. 
    1533                        j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
    1534                        j = MERGE( j           , nyn, j            > nyn )
    1535                        i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
    1536                        i = MERGE( i           , nxr, i            > nxr )
    1537                    
     3046!--                    Surface data is only available on inner subdomains, not
     3047!--                    on ghost points. Hence, limit the indices.
     3048                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
     3049                       j = MERGE( j           , nyn, j            < nyn )
     3050                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
     3051                       i = MERGE( i           , nxr, i            < nxr )
     3052
    15383053                       DO  mm = surf_lsm_h%start_index(j,i),                   &
    15393054                                surf_lsm_h%end_index(j,i)
    … …  
    15463061                    ENDDO
    15473062                 ENDIF
    1548                  
    1549               CASE ( 'rsd' )
     3063
     3064              CASE ( 'rsd' ) ! shortwave in
    15503065                 IF ( radiation )  THEN
    1551                     DO  m = 1, vmea(l)%ns
    1552                        k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' )
    1553                        j = vmea(l)%j(m)
    1554                        i = vmea(l)%i(m)
    1555                    
    1556                        vmea(l)%measured_vars(m,n) = rad_sw_in(k,j,i)
    1557                     ENDDO
     3066                    IF ( radiation_scheme /= 'rrtmg' )  THEN
     3067                       DO  m = 1, vmea(l)%ns
     3068                          k = 0
     3069                          j = vmea(l)%j(m)
     3070                          i = vmea(l)%i(m)
     3071                          vmea(l)%measured_vars(m,n) = rad_sw_in(k,j,i)
     3072                       ENDDO
     3073                    ELSE
     3074                       DO  m = 1, vmea(l)%ns
     3075                          k = vmea(l)%k(m)
     3076                          j = vmea(l)%j(m)
     3077                          i = vmea(l)%i(m)
     3078                          vmea(l)%measured_vars(m,n) = rad_sw_in(k,j,i)
     3079                       ENDDO
     3080                    ENDIF
    15583081                 ENDIF
    1559                  
    1560               CASE ( 'rsu' )
     3082
     3083              CASE ( 'rsu' ) ! shortwave out
    15613084                 IF ( radiation )  THEN
    1562                     DO  m = 1, vmea(l)%ns
    1563                        k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' )
    1564                        j = vmea(l)%j(m)
    1565                        i = vmea(l)%i(m)
    1566                    
    1567                        vmea(l)%measured_vars(m,n) = rad_sw_out(k,j,i)
    1568                     ENDDO
     3085                    IF ( radiation_scheme /= 'rrtmg' )  THEN
     3086                       DO  m = 1, vmea(l)%ns
     3087                          k = 0
     3088                          j = vmea(l)%j(m)
     3089                          i = vmea(l)%i(m)
     3090                          vmea(l)%measured_vars(m,n) = rad_sw_out(k,j,i)
     3091                       ENDDO
     3092                    ELSE
     3093                       DO  m = 1, vmea(l)%ns
     3094                          k = vmea(l)%k(m)
     3095                          j = vmea(l)%j(m)
     3096                          i = vmea(l)%i(m)
     3097                          vmea(l)%measured_vars(m,n) = rad_sw_out(k,j,i)
     3098                       ENDDO
     3099                    ENDIF
    15693100                 ENDIF
    1570                  
    1571               CASE ( 'rlu' )
     3101
     3102              CASE ( 'rlu' ) ! longwave out
    15723103                 IF ( radiation )  THEN
    1573                     DO  m = 1, vmea(l)%ns
    1574                        k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' )
    1575                        j = vmea(l)%j(m)
    1576                        i = vmea(l)%i(m)
    1577                    
    1578                        vmea(l)%measured_vars(m,n) = rad_lw_out(k,j,i)
    1579                     ENDDO
     3104                    IF ( radiation_scheme /= 'rrtmg' )  THEN
     3105                       DO  m = 1, vmea(l)%ns
     3106                          k = 0
     3107                          j = vmea(l)%j(m)
     3108                          i = vmea(l)%i(m)
     3109                          vmea(l)%measured_vars(m,n) = rad_lw_out(k,j,i)
     3110                       ENDDO
     3111                    ELSE
     3112                       DO  m = 1, vmea(l)%ns
     3113                          k = vmea(l)%k(m)
     3114                          j = vmea(l)%j(m)
     3115                          i = vmea(l)%i(m)
     3116                          vmea(l)%measured_vars(m,n) = rad_lw_out(k,j,i)
     3117                       ENDDO
     3118                    ENDIF
    15803119                 ENDIF
    1581                  
    1582               CASE ( 'rld' )
     3120
     3121              CASE ( 'rld' ) ! longwave in
    15833122                 IF ( radiation )  THEN
    1584                     DO  m = 1, vmea(l)%ns
    1585                        k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' )
    1586                        j = vmea(l)%j(m)
    1587                        i = vmea(l)%i(m)
    1588                    
    1589                        vmea(l)%measured_vars(m,n) = rad_lw_in(k,j,i)
    1590                     ENDDO
     3123                    IF ( radiation_scheme /= 'rrtmg' )  THEN
     3124                       DO  m = 1, vmea(l)%ns
     3125                          k = 0
     3126!
     3127!--                       Surface data is only available on inner subdomains,
     3128!--                       not on ghost points. Hence, limit the indices.
     3129                          j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
     3130                          j = MERGE( j           , nyn, j            < nyn )
     3131                          i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
     3132                          i = MERGE( i           , nxr, i            < nxr )
     3133
     3134                          vmea(l)%measured_vars(m,n) = rad_lw_in(k,j,i)
     3135                       ENDDO
     3136                    ELSE
     3137                       DO  m = 1, vmea(l)%ns
     3138                          k = vmea(l)%k(m)
     3139                          j = vmea(l)%j(m)
     3140                          i = vmea(l)%i(m)
     3141                          vmea(l)%measured_vars(m,n) = rad_lw_in(k,j,i)
     3142                       ENDDO
     3143                    ENDIF
    15913144                 ENDIF
    1592                  
    1593               CASE ( 'rsddif' )
     3145
     3146              CASE ( 'rsddif' ) ! shortwave in, diffuse part
    15943147                 IF ( radiation )  THEN
    15953148                    DO  m = 1, vmea(l)%ns
    15963149                       j = vmea(l)%j(m)
    15973150                       i = vmea(l)%i(m)
    1598                    
     3151
    15993152                       vmea(l)%measured_vars(m,n) = rad_sw_in_diff(j,i)
    16003153                    ENDDO
    16013154                 ENDIF
    1602                  
    1603               CASE ( 't_soil' )
     3155
     3156              CASE ( 't_soil' ) ! soil and wall temperature
    16043157                 DO  m = 1, vmea(l)%ns_soil
    1605                     i = vmea(l)%i_soil(m)
    1606                     j = vmea(l)%j_soil(m)
     3158                    j = MERGE( vmea(l)%j_soil(m), nys, vmea(l)%j_soil(m) > nys )
     3159                    j = MERGE( j                , nyn, j                 < nyn )
     3160                    i = MERGE( vmea(l)%i_soil(m), nxl, vmea(l)%i_soil(m) > nxl )
     3161                    i = MERGE( i                , nxr, i                 < nxr )
    16073162                    k = vmea(l)%k_soil(m)
    1608                    
     3163
    16093164                    match_lsm = surf_lsm_h%start_index(j,i) <=                 &
    16103165                                surf_lsm_h%end_index(j,i)
    16113166                    match_usm = surf_usm_h%start_index(j,i) <=                 &
    16123167                                surf_usm_h%end_index(j,i)
    1613                                
     3168
    16143169                    IF ( match_lsm )  THEN
    16153170                       mm = surf_lsm_h%start_index(j,i)
    16163171                       vmea(l)%measured_vars_soil(m,n) = t_soil_h%var_2d(k,mm)
    16173172                    ENDIF
    1618                    
     3173
    16193174                    IF ( match_usm )  THEN
    16203175                       mm = surf_usm_h%start_index(j,i)
    … …  
    16223177                    ENDIF
    16233178                 ENDDO
    1624                  
    1625               CASE ( 'm_soil' )
     3179
     3180              CASE ( 'm_soil' ) ! soil moisture
    16263181                 DO  m = 1, vmea(l)%ns_soil
    1627                     i = vmea(l)%i_soil(m)
    1628                     j = vmea(l)%j_soil(m)
     3182                    j = MERGE( vmea(l)%j_soil(m), nys, vmea(l)%j_soil(m) > nys )
     3183                    j = MERGE( j                , nyn, j                 < nyn )
     3184                    i = MERGE( vmea(l)%i_soil(m), nxl, vmea(l)%i_soil(m) > nxl )
     3185                    i = MERGE( i                , nxr, i                 < nxr )
    16293186                    k = vmea(l)%k_soil(m)
    1630                    
     3187
    16313188                    match_lsm = surf_lsm_h%start_index(j,i) <=                 &
    16323189                                surf_lsm_h%end_index(j,i)
    1633                                
     3190
    16343191                    IF ( match_lsm )  THEN
    16353192                       mm = surf_lsm_h%start_index(j,i)
    16363193                       vmea(l)%measured_vars_soil(m,n) = m_soil_h%var_2d(k,mm)
    16373194                    ENDIF
    1638                    
     3195
     3196                 ENDDO
     3197
     3198              CASE ( 'ts' ) ! surface temperature
     3199                 DO  m = 1, vmea(l)%ns
     3200!
     3201!--                 Surface data is only available on inner subdomains, not
     3202!--                 on ghost points. Hence, limit the indices.
     3203                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
     3204                    j = MERGE( j           , nyn, j            < nyn )
     3205                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
     3206                    i = MERGE( i           , nxr, i            < nxr )
     3207
     3208                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
     3209                             surf_def_h(0)%end_index(j,i)
     3210                       vmea(l)%measured_vars(m,n) = surf_def_h(0)%pt_surface(mm)
     3211                    ENDDO
     3212                    DO  mm = surf_lsm_h%start_index(j,i),                      &
     3213                             surf_lsm_h%end_index(j,i)
     3214                       vmea(l)%measured_vars(m,n) = surf_lsm_h%pt_surface(mm)
     3215                    ENDDO
     3216                    DO  mm = surf_usm_h%start_index(j,i),                      &
     3217                             surf_usm_h%end_index(j,i)
     3218                       vmea(l)%measured_vars(m,n) = surf_usm_h%pt_surface(mm)
     3219                    ENDDO
     3220                 ENDDO
     3221
     3222              CASE ( 'lwp' ) ! liquid water path
     3223                 IF ( ASSOCIATED( ql ) )  THEN
     3224                    DO  m = 1, vmea(l)%ns
     3225                       j = vmea(l)%j(m)
     3226                       i = vmea(l)%i(m)
     3227
     3228                       vmea(l)%measured_vars(m,n) = SUM( ql(nzb:nzt,j,i)       &
     3229                                                       * dzw(1:nzt+1) )        &
     3230                                                  * rho_surface
     3231                    ENDDO
     3232                 ENDIF
     3233
     3234              CASE ( 'ps' ) ! surface pressure
     3235                 vmea(l)%measured_vars(:,n) = surface_pressure
     3236
     3237              CASE ( 'pswrtg' ) ! platform speed above ground
     3238                 vmea(l)%measured_vars(:,n) = 0.0_wp
     3239
     3240              CASE ( 'pswrta' ) ! platform speed in air
     3241                 vmea(l)%measured_vars(:,n) = 0.0_wp
     3242
     3243              CASE ( 't_lw' ) ! water temperature
     3244                 DO  m = 1, vmea(l)%ns
     3245!
     3246!--                 Surface data is only available on inner subdomains, not
     3247!--                 on ghost points. Hence, limit the indices.
     3248                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
     3249                    j = MERGE( j           , nyn, j            < nyn )
     3250                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
     3251                    i = MERGE( i           , nxr, i            < nxr )
     3252
     3253                    DO  mm = surf_lsm_h%start_index(j,i),                      &
     3254                             surf_lsm_h%end_index(j,i)
     3255                       IF ( surf_lsm_h%water_surface(m) )                      &
     3256                          vmea(l)%measured_vars(m,n) = t_soil_h%var_2d(nzt,m)
     3257                    ENDDO
     3258
    16393259                 ENDDO
    16403260!
    … …  
    16503270
    16513271     ENDDO
    1652          
     3272
    16533273  END SUBROUTINE vm_sampling
    1654  
     3274
    16553275
    16563276 END MODULE virtual_measurement_mod
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |