Changeset 3704
- Timestamp:
- Jan 29, 2019 7:51:41 PM (6 years ago)
- Location:
- palm/trunk
- Files:
-
- 3 added
- 9 edited
-
SCRIPTS/.palm.iofiles (modified) (1 diff)
-
SOURCE/check_open.f90 (modified) (2 diffs)
-
SOURCE/module_interface.f90 (modified) (3 diffs)
-
SOURCE/nesting_offl_mod.f90 (modified) (2 diffs)
-
SOURCE/netcdf_data_input_mod.f90 (modified) (4 diffs)
-
SOURCE/radiation_model_mod.f90 (modified) (2 diffs)
-
SOURCE/time_integration.f90 (modified) (3 diffs)
-
SOURCE/urban_surface_mod.f90 (modified) (3 diffs)
-
SOURCE/virtual_measurement_mod.f90 (modified) (50 diffs)
-
UTIL/combine_virtual_measurements (added)
-
UTIL/combine_virtual_measurements/combine_virtual_measurements.f90 (added)
-
UTIL/combine_virtual_measurements/vm_parin (added)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SCRIPTS/.palm.iofiles
r3619 r3704 60 60 SURFACE_DATA_AV_BIN* out:tr * $base_data/$run_identifier/OUTPUT _av_surf bin 61 61 # 62 VIRTUAL_MEAS_BIN* out * $base_data_output/$run_identifier/OUTPUT _vmeas bin 63 # 62 64 DATA_1D_PR_NETCDF* out:tr * $base_data/$run_identifier/OUTPUT _pr nc 63 65 DATA_1D_SP_NETCDF* out:tr * $base_data/$run_identifier/OUTPUT _sp nc -
palm/trunk/SOURCE/check_open.f90
r3655 r3704 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Open binary files for virtual measurements 23 23 ! 24 24 ! Former revisions: … … 470 470 OPEN ( 26, FILE='SURFACE_DATA_AV_BIN'//TRIM( coupling_char )//myid_char, & 471 471 FORM='UNFORMATTED', POSITION='APPEND' ) 472 472 473 CASE ( 27 ) 474 ! 475 !-- Binary files for virtual measurement data 476 OPEN ( 27, FILE='VIRTUAL_MEAS_BIN'//TRIM( coupling_char )//myid_char,& 477 FORM='UNFORMATTED', POSITION='APPEND' ) 473 478 474 479 CASE ( 30 ) -
palm/trunk/SOURCE/module_interface.f90
r3701 r3704 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Add last_actions for virtual measurements 23 23 ! 24 24 ! Former revisions: … … 290 290 291 291 USE virtual_measurement_mod, & 292 ONLY: vm_parin, & 293 vm_init 292 ONLY: vm_init, & 293 vm_last_actions, & 294 vm_parin 294 295 295 296 USE wind_turbine_model_mod, & … … 1227 1228 SUBROUTINE module_interface_last_actions 1228 1229 1229 1230 IF ( virtual_measurement ) CALL vm_last_actions 1230 1231 IF ( user_module_enabled ) CALL user_last_actions 1231 1232 -
palm/trunk/SOURCE/nesting_offl_mod.f90
r3655 r3704 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Formatting adjustments 23 23 ! 24 24 ! Former revisions: … … 623 623 ENDIF 624 624 #else 625 u_ref = u_ref_l626 v_ref = v_ref_l627 IF ( humidity ) q_ref = q_ref_l628 IF ( .NOT. neutral ) pt_ref = pt_ref_l625 u_ref = u_ref_l 626 v_ref = v_ref_l 627 IF ( humidity ) q_ref = q_ref_l 628 IF ( .NOT. neutral ) pt_ref = pt_ref_l 629 629 #endif 630 630 ! 631 !-- Average data. Note, reference profiles up to nzt are derived from lateral 632 !-- boundaries, at the model top it is derived from the top boundary. Thus, 633 !-- number of input data is different from nzb:nzt compared to nzt+1. 634 !-- Derived from lateral boundaries. 635 u_ref(nzb:nzt) = u_ref(nzb:nzt) / REAL( 2.0_wp * ( ny + 1 + nx ), & 636 KIND = wp ) 637 v_ref(nzb:nzt) = v_ref(nzb:nzt) / REAL( 2.0_wp * ( ny + nx + 1 ), & 638 KIND = wp ) 639 IF ( humidity ) & 640 q_ref(nzb:nzt) = q_ref(nzb:nzt) / REAL( 2.0_wp * ( ny + 1 + nx + 1 ), & 641 KIND = wp ) 642 IF ( .NOT. neutral ) & 643 pt_ref(nzb:nzt) = pt_ref(nzb:nzt) / REAL( 2.0_wp * ( ny + 1 + nx + 1 ), & 644 KIND = wp ) 645 ! 646 !-- Derived from top boundary. 647 u_ref(nzt+1) = u_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx ), KIND = wp ) 648 v_ref(nzt+1) = v_ref(nzt+1) / REAL( ( ny ) * ( nx + 1 ), KIND = wp ) 649 IF ( humidity ) & 650 q_ref(nzt+1) = q_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx + 1 ), KIND = wp ) 651 IF ( .NOT. neutral ) & 652 pt_ref(nzt+1) = pt_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx + 1 ), KIND = wp ) 653 ! 654 !-- Write onto init profiles, which are used for damping 655 u_init = u_ref 656 v_init = v_ref 657 IF ( humidity ) q_init = q_ref 658 IF ( .NOT. neutral ) pt_init = pt_ref 659 ! 660 !-- Set bottom boundary condition 661 IF ( humidity ) q_init(nzb) = q_init(nzb+1) 662 IF ( .NOT. neutral ) pt_init(nzb) = pt_init(nzb+1) 663 ! 664 !-- Further, adjust Rayleigh damping height in case of time-changing conditions. 665 !-- Therefore, calculate boundary-layer depth first. 666 CALL calc_zi 667 CALL adjust_sponge_layer 668 669 ! 670 !-- Update geostrophic wind components from dynamic input file. 671 DO k = nzb+1, nzt 672 ug(k) = interpolate_in_time( nest_offl%ug(0,k), nest_offl%ug(1,k), & 673 fac_dt ) 674 vg(k) = interpolate_in_time( nest_offl%vg(0,k), nest_offl%vg(1,k), & 675 fac_dt ) 676 ENDDO 677 ug(nzt+1) = ug(nzt) 678 vg(nzt+1) = vg(nzt) 679 680 CALL cpu_log( log_point(58), 'offline nesting', 'stop' ) 631 !-- Average data. Note, reference profiles up to nzt are derived from lateral 632 !-- boundaries, at the model top it is derived from the top boundary. Thus, 633 !-- number of input data is different from nzb:nzt compared to nzt+1. 634 !-- Derived from lateral boundaries. 635 u_ref(nzb:nzt) = u_ref(nzb:nzt) / REAL( 2.0_wp * ( ny + 1 + nx ), & 636 KIND = wp ) 637 v_ref(nzb:nzt) = v_ref(nzb:nzt) / REAL( 2.0_wp * ( ny + nx + 1 ), & 638 KIND = wp ) 639 IF ( humidity ) & 640 q_ref(nzb:nzt) = q_ref(nzb:nzt) / REAL( 2.0_wp * & 641 ( ny + 1 + nx + 1 ), & 642 KIND = wp ) 643 IF ( .NOT. neutral ) & 644 pt_ref(nzb:nzt) = pt_ref(nzb:nzt) / REAL( 2.0_wp * & 645 ( ny + 1 + nx + 1 ), & 646 KIND = wp ) 647 ! 648 !-- Derived from top boundary. 649 u_ref(nzt+1) = u_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx ), KIND = wp ) 650 v_ref(nzt+1) = v_ref(nzt+1) / REAL( ( ny ) * ( nx + 1 ), KIND = wp ) 651 IF ( humidity ) & 652 q_ref(nzt+1) = q_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx + 1 ), & 653 KIND = wp ) 654 IF ( .NOT. neutral ) & 655 pt_ref(nzt+1) = pt_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx + 1 ), & 656 KIND = wp ) 657 ! 658 !-- Write onto init profiles, which are used for damping 659 u_init = u_ref 660 v_init = v_ref 661 IF ( humidity ) q_init = q_ref 662 IF ( .NOT. neutral ) pt_init = pt_ref 663 ! 664 !-- Set bottom boundary condition 665 IF ( humidity ) q_init(nzb) = q_init(nzb+1) 666 IF ( .NOT. neutral ) pt_init(nzb) = pt_init(nzb+1) 667 ! 668 !-- Further, adjust Rayleigh damping height in case of time-changing conditions. 669 !-- Therefore, calculate boundary-layer depth first. 670 CALL calc_zi 671 CALL adjust_sponge_layer 672 673 ! 674 !-- Update geostrophic wind components from dynamic input file. 675 DO k = nzb+1, nzt 676 ug(k) = interpolate_in_time( nest_offl%ug(0,k), nest_offl%ug(1,k), & 677 fac_dt ) 678 vg(k) = interpolate_in_time( nest_offl%vg(0,k), nest_offl%vg(1,k), & 679 fac_dt ) 680 ENDDO 681 ug(nzt+1) = ug(nzt) 682 vg(nzt+1) = vg(nzt) 683 684 CALL cpu_log( log_point(58), 'offline nesting', 'stop' ) 681 685 682 686 END SUBROUTINE nesting_offl_bc -
palm/trunk/SOURCE/netcdf_data_input_mod.f90
r3655 r3704 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Interface for attribute input of 8-bit and 32-bit integer 23 23 ! 24 24 ! Former revisions: … … 777 777 778 778 INTERFACE netcdf_data_input_att 779 MODULE PROCEDURE netcdf_data_input_att_int 779 MODULE PROCEDURE netcdf_data_input_att_int8 780 MODULE PROCEDURE netcdf_data_input_att_int32 780 781 MODULE PROCEDURE netcdf_data_input_att_real 781 782 MODULE PROCEDURE netcdf_data_input_att_string … … 1174 1175 ! Description: 1175 1176 ! ------------ 1176 !> Read a global integer attribute 1177 !------------------------------------------------------------------------------! 1178 SUBROUTINE netcdf_data_input_att_int( val, search_string, id_mod, & 1179 input_file, global, openclose, & 1180 variable_name ) 1177 !> Read a global 8-bit integer attribute 1178 !------------------------------------------------------------------------------! 1179 SUBROUTINE netcdf_data_input_att_int8( val, search_string, id_mod, & 1180 input_file, global, openclose, & 1181 variable_name ) 1182 1183 IMPLICIT NONE 1184 1185 CHARACTER(LEN=*) :: search_string !< name of the attribue 1186 1187 CHARACTER(LEN=*) :: input_file !< name of input file 1188 CHARACTER(LEN=*) :: openclose !< string indicating whether NetCDF needs to be opend or closed 1189 CHARACTER(LEN=*) :: variable_name !< string indicating whether NetCDF needs to be opend or closed 1190 1191 INTEGER(iwp) :: id_mod !< NetCDF id of input file 1192 INTEGER(KIND=1) :: val !< value of the attribute 1193 1194 LOGICAL :: global !< flag indicating a global or a variable's attribute 1195 1196 #if defined ( __netcdf ) 1197 ! 1198 !-- Open file in read-only mode 1199 IF ( openclose == 'open' ) THEN 1200 CALL open_read_file( TRIM( input_file ) // TRIM( coupling_char ), & 1201 id_mod ) 1202 ENDIF 1203 ! 1204 !-- Read global attribute 1205 IF ( global ) THEN 1206 CALL get_attribute( id_mod, search_string, val, global ) 1207 ! 1208 !-- Read variable attribute 1209 ELSE 1210 CALL get_attribute( id_mod, search_string, val, global, variable_name ) 1211 ENDIF 1212 ! 1213 !-- Finally, close input file 1214 IF ( openclose == 'close' ) CALL close_input_file( id_mod ) 1215 #endif 1216 1217 END SUBROUTINE netcdf_data_input_att_int8 1218 1219 !------------------------------------------------------------------------------! 1220 ! Description: 1221 ! ------------ 1222 !> Read a global 32-bit integer attribute 1223 !------------------------------------------------------------------------------! 1224 SUBROUTINE netcdf_data_input_att_int32( val, search_string, id_mod, & 1225 input_file, global, openclose, & 1226 variable_name ) 1181 1227 1182 1228 IMPLICIT NONE … … 1214 1260 #endif 1215 1261 1216 END SUBROUTINE netcdf_data_input_att_int 1262 END SUBROUTINE netcdf_data_input_att_int32 1217 1263 1218 1264 !------------------------------------------------------------------------------! -
palm/trunk/SOURCE/radiation_model_mod.f90
r3700 r3704 23 23 ! Current revisions: 24 24 ! ------------------ 25 ! 25 ! Make variables that are sampled in virtual measurement module public 26 26 ! 27 27 ! Former revisions: … … 1276 1276 idsvf, ndsvf, idcsf, ndcsf, kdcsf, pct, & 1277 1277 radiation_interactions, startwall, startland, endland, endwall, & 1278 skyvf, skyvft, radiation_interactions_on, average_radiation 1278 skyvf, skyvft, radiation_interactions_on, average_radiation, & 1279 rad_sw_in_diff, rad_sw_in_dir 1279 1280 1280 1281 -
palm/trunk/SOURCE/time_integration.f90
r3684 r3704 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Data output for virtual measurements added 23 23 ! 24 24 ! Former revisions: … … 632 632 633 633 USE virtual_measurement_mod, & 634 ONLY: vm_ sampling, vm_time_start634 ONLY: vm_data_output, vm_sampling, vm_time_start 635 635 636 636 IMPLICIT NONE … … 1612 1612 !-- Take virtual measurements 1613 1613 IF ( virtual_measurement .AND. & 1614 vm_time_start <= time_since_reference_point ) CALL vm_sampling 1614 vm_time_start <= time_since_reference_point ) THEN 1615 CALL vm_sampling 1616 CALL vm_data_output 1617 ENDIF 1615 1618 1616 1619 ! -
palm/trunk/SOURCE/urban_surface_mod.f90
r3685 r3704 23 23 ! Current revisions: 24 24 ! ------------------ 25 ! 25 ! make nzb_wall public, required for virtual-measurements 26 26 ! 27 27 ! Former revisions: … … 1150 1150 !-- Public parameters, constants and initial values 1151 1151 PUBLIC usm_anthropogenic_heat, usm_material_model, usm_wall_mod, & 1152 usm_green_heat_model, building_pars, &1153 nz t_wall, t_wall_h, t_wall_v,&1152 usm_green_heat_model, building_pars, & 1153 nzb_wall, nzt_wall, t_wall_h, t_wall_v, & 1154 1154 t_window_h, t_window_v, building_type 1155 1155 … … 7469 7469 IMPLICIT NONE 7470 7470 7471 INTEGER(iwp) 7471 INTEGER(iwp) :: i, j, k, l, m !< running indices 7472 7472 7473 7473 INTEGER(iwp) :: i_off !< offset to determine index of surface element, seen from atmospheric grid point, for x -
palm/trunk/SOURCE/virtual_measurement_mod.f90
r3665 r3704 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! - initialization revised 23 ! - binary data output 24 ! - list of allowed variables extended 23 25 ! 24 26 ! Former revisions: 25 27 ! ----------------- 26 28 ! $Id$ 27 ! unused variables removed28 !29 ! 3522 2018-11-13 12:14:36Z suehring30 29 ! Sampling of variables 31 30 ! … … 44 43 !> The module acts as an interface between 'real-world' observations and 45 44 !> model simulations. Virtual measurements will be taken in the model at the 46 !> coordinates representative for the 'real-world' measurement positions.45 !> coordinates representative for the 'real-world' observation coordinates. 47 46 !> More precisely, coordinates and measured quanties will be read from a 48 47 !> NetCDF file which contains all required information. In the model, … … 53 52 !> 54 53 !> @todo list_of_allowed variables needs careful checking 55 !> @todo output (binary or NetCDF) needs to be implemented56 !> @todo clean-up anything from current test modus57 54 !> @todo Check if sign of surface fluxes for heat, radiation, etc., follows 58 55 !> the (UC)2 standard … … 71 68 72 69 USE control_parameters, & 73 ONLY: air_chemistry, dz, humidity, neutral, message_string,&74 virtual_measurement70 ONLY: air_chemistry, dz, humidity, io_blocks, io_group, neutral, & 71 message_string, time_since_reference_point, virtual_measurement 75 72 76 73 USE cpulog, & … … 81 78 82 79 USE indices, & 83 ONLY: nzb, nzt, nxl, nxr, nys, nyn, nx, ny 80 ONLY: nzb, nzt, nxl, nxr, nys, nyn, nx, ny, wall_flags_0 84 81 85 82 USE kinds 83 84 USE netcdf_data_input_mod, & 85 ONLY: init_model 86 87 USE pegrid 88 89 USE surface_mod, & 90 ONLY: surf_lsm_h, surf_usm_h 91 92 USE land_surface_model_mod, & 93 ONLY: nzb_soil, nzs, nzt_soil, zs, t_soil_h, m_soil_h 94 95 USE radiation_model_mod 96 97 USE urban_surface_mod, & 98 ONLY: nzb_wall, nzt_wall, t_wall_h 86 99 87 100 88 101 IMPLICIT NONE 102 103 TYPE virt_general 104 INTEGER(iwp) :: id_vm !< NetCDF file id for virtual measurements 105 INTEGER(iwp) :: nvm = 0 !< number of virtual measurements 106 END TYPE virt_general 89 107 90 108 TYPE virt_mea 91 109 92 CHARACTER(LEN=100) :: feature_type !< type of the measurement 93 CHARACTER(LEN=100) :: site !< name of the measurement site 110 CHARACTER(LEN=100) :: feature_type !< type of the measurement 111 CHARACTER(LEN=100) :: filename_original !< name of the original file 112 CHARACTER(LEN=100) :: site !< name of the measurement site 94 113 95 114 CHARACTER(LEN=10), DIMENSION(:), ALLOCATABLE :: measured_vars_name !< name of the measured variables 96 115 97 INTEGER(iwp) :: ns = 0 !< total number of observation points for a site on subdomain, i.e. sum of all trajectories 98 INTEGER(iwp) :: ntraj !< number of trajectories of a measurement 99 INTEGER(iwp) :: nvar !< number of measured variables 116 INTEGER(iwp) :: ns = 0 !< number of observation coordinates on subdomain, for atmospheric measurements 117 INTEGER(iwp) :: ns_tot = 0 !< total number of observation coordinates, for atmospheric measurements 118 INTEGER(iwp) :: ntraj !< number of trajectories of a measurement 119 INTEGER(iwp) :: nvar !< number of measured variables (atmosphere + soil) 120 121 INTEGER(iwp) :: ns_soil = 0 !< number of observation coordinates on subdomain, for soil measurements 122 INTEGER(iwp) :: ns_soil_tot = 0 !< total number of observation coordinates, for soil measurements 100 123 101 124 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: dim_t !< number observations individual for each trajectory or station that are no _FillValues 102 125 103 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: i !< grid index for measurement position in x-direction 104 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: j !< grid index for measurement position in y-direction 105 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k !< grid index for measurement position in k-direction 126 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: i !< grid index for measurement position in x-direction 127 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: j !< grid index for measurement position in y-direction 128 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k !< grid index for measurement position in k-direction 129 130 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: i_soil !< grid index for measurement position in x-direction 131 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: j_soil !< grid index for measurement position in y-direction 132 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k_soil !< grid index for measurement position in k-direction 106 133 107 134 LOGICAL :: trajectory = .FALSE. !< flag indicating that the observation is a mobile observation 108 135 LOGICAL :: timseries = .FALSE. !< flag indicating that the observation is a stationary point measurement 109 136 LOGICAL :: timseries_profile = .FALSE. !< flag indicating that the observation is a stationary profile measurement 137 LOGICAL :: soil_sampling = .FALSE. !< flag indicating that soil state variables were sampled 110 138 111 REAL(wp) :: fill_eutm !< fill value for UTM coordinates in case of missing values 112 REAL(wp) :: fill_nutm !< fill value for UTM coordinates in case of missing values 113 REAL(wp) :: fill_zag !< fill value for heigth coordinates in case of missing values 114 REAL(wp) :: fillout = -999.9 !< fill value for output in case a observation is taken from inside a building 115 REAL(wp) :: origin_x_obs !< origin of the observation in UTM coordiates in x-direction 116 REAL(wp) :: origin_y_obs !< origin of the observation in UTM coordiates in y-direction 139 REAL(wp) :: fill_eutm !< fill value for UTM coordinates in case of missing values 140 REAL(wp) :: fill_nutm !< fill value for UTM coordinates in case of missing values 141 REAL(wp) :: fill_zag !< fill value for heigth coordinates in case of missing values 142 REAL(wp) :: fillout = -999.9 !< fill value for output in case a observation is taken from inside a building 143 REAL(wp) :: origin_x_obs !< origin of the observation in UTM coordiates in x-direction 144 REAL(wp) :: origin_y_obs !< origin of the observation in UTM coordiates in y-direction 145 146 REAL(wp), DIMENSION(:), ALLOCATABLE :: z_ag !< measurement height above ground level 147 REAL(wp), DIMENSION(:), ALLOCATABLE :: depth !< measurement depth in soil 117 148 118 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: measured_vars !< measured variables 149 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: measured_vars !< measured variables 150 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: measured_vars_soil !< measured variables 119 151 120 152 END TYPE virt_mea … … 122 154 CHARACTER(LEN=5) :: char_eutm = "E_UTM" !< dimension name for UTM coordinate easting 123 155 CHARACTER(LEN=11) :: char_feature = "featureType" !< attribute name for feature type 156 157 ! This need to generalized 158 CHARACTER(LEN=8) :: char_filename = "filename" !< attribute name for filename 159 CHARACTER(LEN=11) :: char_soil = "soil_sample" !< attribute name for soil sampling indication 124 160 CHARACTER(LEN=10) :: char_fillvalue = "_FillValue" !< variable attribute name for _FillValue 125 161 CHARACTER(LEN=18) :: char_mv = "measured_variables" !< variable name for the array with the measured variable names … … 133 169 CHARACTER(LEN=10) :: type_traj = 'trajectory' !< name of line measurements 134 170 CHARACTER(LEN=17) :: type_tspr = 'timeSeriesProfile' !< name of stationary profile measurements 171 172 CHARACTER(LEN=6), DIMENSION(1:5) :: soil_vars = (/ & !< list of soil variables 173 't_soil', & 174 'm_soil', & 175 'lwc ', & 176 'lwcs ', & 177 'smp ' /) 178 179 CHARACTER(LEN=10), DIMENSION(0:1,1:8) :: chem_vars = RESHAPE( (/ & 180 'mcpm1 ', 'kc_PM1 ', & 181 'mcpm2p5 ', 'kc_PM2.5 ', & 182 'mcpm25 ', 'kc_PM25 ', & 183 'mcpm10 ', 'kc_PM10 ', & 184 'mfno2 ', 'kc_NO2 ', & 185 'mfno ', 'kc_NO ', & 186 'tro3 ', 'kc_O3 ', & 187 'mfco ', 'kc_CO ' & 188 /), (/ 2, 8 /) ) 135 189 ! 136 190 !-- MS: List requires careful revision! 137 CHARACTER(LEN=10), DIMENSION(1: 47), PARAMETER :: list_allowed_variables = & !< variables that can be sampled in PALM191 CHARACTER(LEN=10), DIMENSION(1:54), PARAMETER :: list_allowed_variables = & !< variables that can be sampled in PALM 138 192 (/ 'hfls ', & ! surface latent heat flux (W/m2) 139 193 'hfss ', & ! surface sensible heat flux (W/m2) … … 143 197 'mcpm1 ', & ! mass concentration of PM1 (kg/m3) 144 198 'mcpm2p5 ', & ! mass concentration of PM2.5 (kg/m3) 145 'mcpm10 ', & ! mass concentration of PM10 (kg/m3)146 199 'mcpm10 ', & ! mass concentration of PM10 (kg/m3) 147 200 'mcco ', & ! mass concentration of CO (kg/m3) … … 173 226 'tsoil ', & ! ? soil temperature - which depth? 174 227 'u ', & ! u-component 228 'utheta ', & ! total eastward kinematic heat flux 175 229 'ua ', & ! eastward wind (is there any difference to u?) 176 230 'v ', & ! v-component 231 'vtheta ', & ! total northward kinematic heat flux 177 232 'va ', & ! northward wind (is there any difference to v?) 178 233 'w ', & ! w-component 234 'wtheta ', & ! total vertical kinematic heat flux 179 235 'rld ', & ! downward longwave radiative flux (W/m2) 180 236 'rlu ', & ! upnward longwave radiative flux (W/m2) … … 182 238 'rsu ', & ! upward shortwave radiative flux (W/m2) 183 239 'rsddif ', & ! downward shortwave diffuse radiative flux (W/m2) 184 'rnds ' & ! surface net downward radiative flux (W/m2) 240 'rnds ', & ! surface net downward radiative flux (W/m2) 241 't_soil ', & 242 'm_soil ', & 243 'lwc ', & 244 'lwcs ', & 245 'smp ' & 185 246 /) 186 187 INTEGER(iwp) :: id_vm !< NetCDF file id for virtual measurements 188 INTEGER(iwp) :: nvm = 0 !< number of virtual measurements 189 ! INTEGER(iwp) :: observation_coverage_xy = 0 !< horizontal distance from the measurement point where observations should be taken in the surrounding 190 ! INTEGER(iwp) :: observation_coverage_z = 0 !< vertical distance from the measurement point where observations should be taken in the surrounding 191 247 248 249 LOGICAL :: global_attribute = .TRUE. !< flag indicating a global attribute 250 LOGICAL :: init = .TRUE. !< flag indicating initialization of data output 192 251 LOGICAL :: use_virtual_measurement = .FALSE. !< Namelist parameter 193 LOGICAL :: global_attribute = .TRUE. !< flag indicating a global attribute194 252 195 253 REAL(wp) :: vm_time_start = 0.0 !< time after virtual measurements should start 196 197 254 255 TYPE( virt_general ) :: vmea_general !< data structure which encompass general variables 198 256 TYPE( virt_mea ), DIMENSION(:), ALLOCATABLE :: vmea !< virtual measurement data structure 199 200 257 201 258 INTERFACE vm_check_parameters … … 203 260 END INTERFACE vm_check_parameters 204 261 262 INTERFACE vm_data_output 263 MODULE PROCEDURE vm_data_output 264 END INTERFACE vm_data_output 265 205 266 INTERFACE vm_init 206 267 MODULE PROCEDURE vm_init 207 268 END INTERFACE vm_init 208 269 270 INTERFACE vm_last_actions 271 MODULE PROCEDURE vm_last_actions 272 END INTERFACE vm_last_actions 273 209 274 INTERFACE vm_parin 210 275 MODULE PROCEDURE vm_parin … … 221 286 ! 222 287 !-- Public interfaces 223 PUBLIC vm_check_parameters, vm_init, vm_parin, vm_sampling 288 PUBLIC vm_check_parameters, vm_data_output, vm_init, vm_last_actions, & 289 vm_parin, vm_sampling 224 290 225 291 ! 226 292 !-- Public variables 227 PUBLIC vmea, vm _time_start293 PUBLIC vmea, vmea_general, vm_time_start 228 294 229 295 CONTAINS … … 251 317 !-- ToDo: Revise this later and remove this requirement. 252 318 IF ( virtual_measurement .AND. .NOT. input_pids_static ) THEN 253 message_string = 'If virtual measurements are taken a static input ' //&319 message_string = 'If virtual measurements are taken, a static input ' //& 254 320 'file is mandatory.' 255 321 CALL message( 'vm_check_parameters', 'PA0000', 1, 2, 0, 6, 0 ) … … 319 385 320 386 USE netcdf_data_input_mod, & 321 ONLY: input_file_vm, netcdf_data_input_get_dimension_length, & 387 ONLY: init_model, input_file_vm, & 388 netcdf_data_input_get_dimension_length, & 322 389 netcdf_data_input_att, netcdf_data_input_var 323 390 … … 328 395 329 396 CHARACTER(LEN=5) :: dum !< dummy string indicate station id 397 CHARACTER(LEN=5) :: dummy_read !< dummy string indicate station id 330 398 CHARACTER(LEN=10), DIMENSION(50) :: measured_variables_file = '' !< array with all measured variables read from NetCDF 331 399 CHARACTER(LEN=10), DIMENSION(50) :: measured_variables = '' !< dummy array with all measured variables that are allowed 332 400 333 !INTEGER(iwp) :: dim_eutm !< dimension size of UTM easting coordinate334 !INTEGER(iwp) :: dim_nutm !< dimension size of UTM northing coordinate401 INTEGER(iwp) :: dim_eutm !< dimension size of UTM easting coordinate 402 INTEGER(iwp) :: dim_nutm !< dimension size of UTM northing coordinate 335 403 INTEGER(iwp) :: dim_ntime !< dimension size of time coordinate 336 !INTEGER(iwp) :: dim_zag !< dimension size of height coordinate337 !INTEGER(iwp) :: i !< grid index of virtual observation point in x-direction338 !INTEGER(iwp) :: ii !< running index over all coordinate points of a measurement404 INTEGER(iwp) :: dim_zag !< dimension size of height coordinate 405 INTEGER(iwp) :: i !< grid index of virtual observation point in x-direction 406 INTEGER(iwp) :: ii !< running index over all coordinate points of a measurement 339 407 INTEGER(iwp) :: is !< grid index of real observation point of the respective station in x-direction 340 !INTEGER(iwp) :: j !< grid index of observation point in x-direction408 INTEGER(iwp) :: j !< grid index of observation point in x-direction 341 409 INTEGER(iwp) :: js !< grid index of real observation point of the respective station in y-direction 342 !INTEGER(iwp) :: k !< grid index of observation point in x-direction410 INTEGER(iwp) :: k !< grid index of observation point in x-direction 343 411 INTEGER(iwp) :: kl !< lower vertical index of surrounding grid points of an observation coordinate 344 412 INTEGER(iwp) :: ks !< grid index of real observation point of the respective station in z-direction … … 352 420 INTEGER(iwp) :: ns !< counter variable for number of observation points on subdomain 353 421 INTEGER(iwp) :: t !< running index over number of trajectories 422 INTEGER(iwp) :: m 423 424 INTEGER(KIND=1):: soil_dum 425 426 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ns_all !< dummy array used to sum-up the number of observation coordinates 354 427 355 428 INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: meas_flag !< mask array indicating measurement positions … … 366 439 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: z_ag !< height coordinate relative to origin_z, temporary variable 367 440 ! 368 !-- Obtain number of virtual measurement stations 369 CALL netcdf_data_input_att( nvm, char_numstations, id_vm, input_file_vm, & 441 !-- Obtain number of sites. Also, pass the 'open' string, in order to initially 442 !-- open the measurement driver. 443 CALL netcdf_data_input_att( vmea_general%nvm, char_numstations, & 444 vmea_general%id_vm, input_file_vm, & 370 445 global_attribute, 'open', '' ) 371 446 372 ! write(9,*) "num stationi", nvm 373 ! flush(9) 374 ! 375 !-- ALLOCATE data structure 376 ALLOCATE( vmea(1:nvm) ) 377 ! 378 !-- Allocate flag array 447 ! 448 !-- Allocate data structure which encompass all required information, such as 449 !-- grid points indicies, absolute UTM coordinates, the measured quantities, 450 !-- etc. . 451 ALLOCATE( vmea(1:vmea_general%nvm) ) 452 ! 453 !-- Allocate flag array. This dummy array is used to identify grid points 454 !-- where virtual measurements should be taken. Please note, at least one 455 !-- ghost point is required, in order to include also the surrounding 456 !-- grid points of the original coordinate. 379 457 ALLOCATE( meas_flag(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ) 380 458 meas_flag = 0 381 382 ! 383 !-- Read station coordinates and further attributes. 384 !-- Note all coordinates are in UTM coordinates. 385 DO l = 1, nvm 386 ! 387 !-- Determine suffix which contains the ID 459 ! 460 !-- Loop over all sites. 461 DO l = 1, vmea_general%nvm 462 ! 463 !-- Determine suffix which contains the ID, ordered according to the number 464 !-- of measurements. 388 465 IF( l < 10 ) THEN 389 466 WRITE( dum, '(I1)') l … … 397 474 WRITE( dum, '(I5)') l 398 475 ENDIF 399 400 CALL netcdf_data_input_att( vmea(l)%origin_x_obs, char_origx & 401 // TRIM( dum ), id_vm, '', global_attribute,& 476 ! 477 !-- Read site coordinates (UTM). 478 CALL netcdf_data_input_att( vmea(l)%origin_x_obs, char_origx // & 479 TRIM( dum ), vmea_general%id_vm, '', & 480 global_attribute, '', '' ) 481 CALL netcdf_data_input_att( vmea(l)%origin_y_obs, char_origy // & 482 TRIM( dum ), vmea_general%id_vm, '', & 483 global_attribute, '', '' ) 484 ! 485 !-- Read site name 486 CALL netcdf_data_input_att( vmea(l)%site, char_site // TRIM( dum ), & 487 vmea_general%id_vm, '', global_attribute, & 402 488 '', '' ) 403 CALL netcdf_data_input_att( vmea(l)%origin_y_obs, char_origy & 404 // TRIM( dum ), id_vm, '', global_attribute,& 489 ! 490 !-- Read type of the measurement (trajectory, profile, timeseries). 491 CALL netcdf_data_input_att( vmea(l)%feature_type, char_feature // & 492 TRIM( dum ), vmea_general%id_vm, '', & 493 global_attribute, '', '' ) 494 ! 495 !-- Read the name of the original file where observational data is stored. 496 CALL netcdf_data_input_att( vmea(l)%filename_original, char_filename // & 497 TRIM( dum ), vmea_general%id_vm, '', & 498 global_attribute, '', '' ) 499 ! 500 !-- Read a flag which indicates that also soil quantities are take at the 501 !-- respective site (is part of the virtual measurement driver). 502 CALL netcdf_data_input_att( soil_dum, char_soil // TRIM( dum ), & 503 vmea_general%id_vm, '', global_attribute, & 405 504 '', '' ) 406 CALL netcdf_data_input_att( vmea(l)%site, char_site & 407 // TRIM( dum ), id_vm, '', global_attribute,& 408 '', '' ) 409 CALL netcdf_data_input_att( vmea(l)%feature_type, char_feature & 410 // TRIM( dum ), id_vm, '', global_attribute,& 411 '', '' ) 412 505 ! 506 !-- Set flag for soil-sampling. 507 IF ( soil_dum == 1 ) vmea(l)%soil_sampling = .TRUE. 413 508 ! 414 509 !--- Set logicals depending on the type of the measurement … … 419 514 ELSEIF ( INDEX( vmea(l)%feature_type, type_traj ) /= 0 ) THEN 420 515 vmea(l)%trajectory = .TRUE. 516 ! 517 !-- Give error message in case the type matches non of the pre-defined types. 421 518 ELSE 422 !423 !-- Give error message424 519 message_string = 'Attribue featureType = ' // & 425 520 TRIM( vmea(l)%feature_type ) // & … … 428 523 ENDIF 429 524 ! 430 !-- Read string with all measured variables at this s tation525 !-- Read string with all measured variables at this site 431 526 measured_variables_file = '' 432 527 CALL netcdf_data_input_var( measured_variables_file, & 433 char_mv // TRIM( dum ), id_vm )434 ! 435 !-- Count the number of measured variables which match with the variables436 !-- w hich are allowed to be measured in PALM. Please note, for some437 !-- NetCDF interal reasons characters end with a NULL, i.e. also empty438 !-- characters contain a NULL. Therefore, check the strings for a Null to439 !-- get the correct character length in order to compare them with the list440 !-- of allowed variables.441 vmea(l)%nvar = 0528 char_mv // TRIM( dum ), vmea_general%id_vm ) 529 ! 530 !-- Count the number of measured variables. Only count variables that match 531 !-- with the allowed variables. 532 !-- Please note, for some NetCDF interal reasons characters end with a NULL, 533 !-- i.e. also empty characters contain a NULL. Therefore, check the strings 534 !-- for a NULL to get the correct character length in order to compare 535 !-- them with the list of allowed variables. 536 vmea(l)%nvar = 0 442 537 DO ll = 1, SIZE( measured_variables_file ) 443 538 IF ( measured_variables_file(ll)(1:1) /= CHAR(0) .AND. & … … 465 560 ENDDO 466 561 ! 467 !-- Allocate array for the measured variables names for the station l.562 !-- Allocate array for the measured variables names for the respective site. 468 563 ALLOCATE( vmea(l)%measured_vars_name(1:vmea(l)%nvar) ) 469 564 … … 474 569 !-- In case of chemistry, check if species is considered in the modelled 475 570 !-- chemistry mechanism. 476 IF ( air_chemistry ) THEN 477 DO ll = 1, vmea(l)%nvar 478 chem_include = .FALSE. 479 DO n = 1, nspec 480 IF ( TRIM( vmea(l)%measured_vars_name(ll) ) == & 481 TRIM( chem_species(n)%name ) ) chem_include = .TRUE. 482 ENDDO 483 IF ( .NOT. chem_include ) THEN 484 message_string = TRIM( vmea(l)%measured_vars_name(ll) ) // & 485 ' is not considered in the modelled ' // & 486 'chemistry mechanism' 487 CALL message( 'vm_init', 'PA0000', 0, 0, 0, 6, 0 ) 488 ENDIF 489 ENDDO 490 ENDIF 491 ! 492 !-- For the actual measurement ID read the UTM coordinates. Based on these, 493 !-- define the index space on each subdomain where measurements should be 494 !-- taken. Note, the entire coordinate arrays will not be stored on data 495 !-- type as this would exceed memory requirements, particularly for 496 !-- trajectory measurements. If no variable will be virtually measured, 497 !-- skip the reading. 571 ! IF ( air_chemistry ) THEN 572 ! DO ll = 1, vmea(l)%nvar 573 ! chem_include = .FALSE. 574 ! DO n = 1, nspec 575 ! IF ( TRIM( vmea(l)%measured_vars_name(ll) ) == & 576 ! TRIM( chem_species(n)%name ) ) chem_include = .TRUE. 577 ! ENDDO 578 ! ! 579 ! !-- Revise this. It should only check for chemistry variables and not for all! 580 ! IF ( .NOT. chem_include ) THEN 581 ! message_string = TRIM( vmea(l)%measured_vars_name(ll) ) // & 582 ! ' is not considered in the modelled ' // & 583 ! 'chemistry mechanism' 584 ! CALL message( 'vm_init', 'PA0000', 0, 0, 0, 6, 0 ) 585 ! ENDIF 586 ! ENDDO 587 ! ENDIF 588 ! 589 !-- Read the UTM coordinates for the actual site. Based on the coordinates, 590 !-- define the grid-index space on each subdomain where virtual measurements 591 !-- should be taken. Note, the entire coordinate arrays will not be stored 592 !-- as this would exceed memory requirements, particularly for trajectory 593 !-- measurements. 498 594 IF ( vmea(l)%nvar > 0 ) THEN 499 595 ! 500 596 !-- For stationary measurements UTM coordinates are just one value and 501 597 !-- its dimension is "station", while for mobile measurements UTM 502 !-- coordinates are arrays. First, inquire dimension length for 503 !-- UTM coordinates. 598 !-- coordinates are arrays depending on the number of trajectories and 599 !-- time, according to (UC)2 standard. First, inquire dimension length 600 !-- of the UTM coordinates. 504 601 IF ( vmea(l)%trajectory ) THEN 505 602 ! 506 603 !-- For non-stationary measurements read the number of trajectories 507 CALL netcdf_data_input_get_dimension_length( id_vm, & 604 !-- and the number of time coordinates. 605 CALL netcdf_data_input_get_dimension_length( vmea_general%id_vm, & 508 606 vmea(l)%ntraj, & 509 607 "traj" // & 510 608 TRIM( dum ) ) 511 CALL netcdf_data_input_get_dimension_length( id_vm, dim_ntime, & 609 CALL netcdf_data_input_get_dimension_length( vmea_general%id_vm, & 610 dim_ntime, & 512 611 "ntime" // & 513 612 TRIM( dum ) ) 514 613 ! 515 !-- For stationary measurements the dimension for UTM coordinates is 1 614 !-- For stationary measurements the dimension for UTM and time 615 !-- coordinates is 1. 516 616 ELSE 517 617 vmea(l)%ntraj = 1 518 618 dim_ntime = 1 519 619 ENDIF 520 521 620 ! 522 621 !- Allocate array which defines individual time frame for each 523 !-- trajectory or station 622 !-- trajectory or station. 524 623 ALLOCATE( vmea(l)%dim_t(1:vmea(l)%ntraj) ) 525 624 ! … … 530 629 ALLOCATE( z_ag(1:vmea(l)%ntraj,1:dim_ntime) ) 531 630 ! 532 !-- Read _FillValue attributes 631 !-- Read _FillValue attributes of the coordinate dimensions. 533 632 CALL netcdf_data_input_att( fill_eutm, char_fillvalue, & 534 id_vm, '', .NOT. global_attribute, '', & 633 vmea_general%id_vm, '', & 634 .NOT. global_attribute, '', & 535 635 char_eutm // TRIM( dum ) ) 536 636 CALL netcdf_data_input_att( fill_nutm, char_fillvalue, & 537 id_vm, '', .NOT. global_attribute, '', & 637 vmea_general%id_vm, '', & 638 .NOT. global_attribute, '', & 538 639 char_nutm // TRIM( dum ) ) 539 640 CALL netcdf_data_input_att( fill_zag, char_fillvalue, & 540 id_vm, '', .NOT. global_attribute, '', & 641 vmea_general%id_vm, '', & 642 .NOT. global_attribute, '', & 541 643 char_zag // TRIM( dum ) ) 542 644 ! … … 544 646 !-- times. 545 647 IF ( vmea(l)%trajectory ) THEN 546 CALL netcdf_data_input_var( e_utm, char_eutm // TRIM( dum ), id_vm, & 648 CALL netcdf_data_input_var( e_utm, char_eutm // TRIM( dum ), & 649 vmea_general%id_vm, & 547 650 0, dim_ntime-1, 0, vmea(l)%ntraj-1 ) 548 CALL netcdf_data_input_var( n_utm, char_nutm // TRIM( dum ), id_vm, & 651 CALL netcdf_data_input_var( n_utm, char_nutm // TRIM( dum ), & 652 vmea_general%id_vm, & 549 653 0, dim_ntime-1, 0, vmea(l)%ntraj-1 ) 550 CALL netcdf_data_input_var( z_ag, char_zag // TRIM( dum ), id_vm, & 654 CALL netcdf_data_input_var( z_ag, char_zag // TRIM( dum ), & 655 vmea_general%id_vm, & 551 656 0, dim_ntime-1, 0, vmea(l)%ntraj-1 ) 552 657 ELSE 553 CALL netcdf_data_input_var( e_utm(1,:), char_eutm // TRIM( dum ), id_vm ) 554 CALL netcdf_data_input_var( n_utm(1,:), char_nutm // TRIM( dum ), id_vm ) 555 CALL netcdf_data_input_var( z_ag(1,:), char_zag // TRIM( dum ), id_vm ) 556 ENDIF 557 ! 558 !-- For testing: 559 e_utm = e_utm - e_utm(1,1) 560 n_utm = n_utm - n_utm(1,1) 561 562 ! 563 !-- First, compute relative x- and y-coordinates with respect to the 564 !-- lower-left origin of the model domain, which is the difference 565 !-- betwen UTM coordinates. 566 ! e_utm(t,1:vmea(l)%dim_t(t)) = e_utm(t,1:vmea(l)%dim_t(t)) & 567 ! - init_model%origin_x 568 ! n_utm(t,1:vmea(l)%dim_t(t)) = n_utm(t,1:vmea(l)%dim_t(t)) & 569 ! - init_model%origin_y 570 658 CALL netcdf_data_input_var( e_utm(1,:), char_eutm // TRIM( dum ), & 659 vmea_general%id_vm ) 660 CALL netcdf_data_input_var( n_utm(1,:), char_nutm // TRIM( dum ), & 661 vmea_general%id_vm ) 662 CALL netcdf_data_input_var( z_ag(1,:), char_zag // TRIM( dum ), & 663 vmea_general%id_vm ) 664 ENDIF 571 665 ! 572 666 !-- Based on UTM coordinates, check if the measurement station or parts … … 575 669 meas_flag = 0 576 670 DO t = 1, vmea(l)%ntraj 671 ! 672 !-- First, compute relative x- and y-coordinates with respect to the 673 !-- lower-left origin of the model domain, which is the difference 674 !-- betwen UTM coordinates. Note, if the origin is not correct, the 675 !-- virtual sites will be misplaced. 676 e_utm(t,1:dim_ntime) = e_utm(t,1:dim_ntime) - init_model%origin_x 677 n_utm(t,1:dim_ntime) = n_utm(t,1:dim_ntime) - init_model%origin_y 577 678 ! 578 679 !-- Determine the individual time coordinate length for each station and … … 580 681 !-- are merged into one file but they do not have the same number of 581 682 !-- points in time, hence, missing values may occur and cannot be 582 !-- processed further. 683 !-- processed further. This is actually a work-around for the specific 684 !-- (UC)2 dataset, but it won't harm in anyway. 583 685 vmea(l)%dim_t(t) = 0 584 686 DO n = 1, dim_ntime … … 587 689 z_ag(t,n) /= fill_zag ) vmea(l)%dim_t(t) = n 588 690 ENDDO 589 590 691 ! 591 692 !-- Compute grid indices relative to origin and check if these are … … 612 713 !-- surrounding coordinate points, but first check whether the 613 714 !-- surrounding coordinate points are on the subdomain. 614 kl = MERGE( ks-1, ks, ks-1 >= nzb .AND. ks-1 > ksurf ) 615 ku = MERGE( ks+1, ks, ks+1 <= nzt+1 ) 616 617 meas_flag(kl:ku,js-1:js+1,is-1:is+1) = & 618 IBSET( meas_flag(kl:ku,js-1:js+1,is-1:is+1), 0 ) 715 kl = MERGE( ks-1, ks, ks-1 >= nzb .AND. ks-1 >= ksurf ) 716 ku = MERGE( ks+1, ks, ks+1 < nzt+1 ) 717 718 DO i = is-1, is+1 719 DO j = js-1, js+1 720 DO k = kl, ku 721 meas_flag(k,j,i) = MERGE( & 722 IBSET( meas_flag(k,j,i), 0 ), & 723 0, & 724 BTEST( wall_flags_0(k,j,i), 0 ) & 725 ) 726 ENDDO 727 ENDDO 728 ENDDO 619 729 ENDIF 620 730 ENDDO … … 622 732 ENDDO 623 733 ! 624 !-- Based on the flag array count the number of observation points. 734 !-- Based on the flag array count the number of of sampling coordinates. 735 !-- Please note, sampling coordinates in atmosphere and soil may be 736 !-- different, as within the soil all levels will be measured. 737 !-- Hence, count individually. Start with atmoshere. 625 738 ns = 0 626 DO i s= nxl-1, nxr+1627 DO j s= nys-1, nyn+1628 DO k s= nzb, nzt+1629 ns = ns + MERGE( 1, 0, BTEST( meas_flag(k s,js,is), 0 ) )739 DO i = nxl-1, nxr+1 740 DO j = nys-1, nyn+1 741 DO k = nzb, nzt+1 742 ns = ns + MERGE( 1, 0, BTEST( meas_flag(k,j,i), 0 ) ) 630 743 ENDDO 631 744 ENDDO 632 745 ENDDO 746 633 747 ! 634 748 !-- Store number of observation points on subdomain and allocate index 635 !-- arrays .749 !-- arrays as well as array containing height information. 636 750 vmea(l)%ns = ns 637 ns = 0638 751 639 752 ALLOCATE( vmea(l)%i(1:vmea(l)%ns) ) 640 753 ALLOCATE( vmea(l)%j(1:vmea(l)%ns) ) 641 754 ALLOCATE( vmea(l)%k(1:vmea(l)%ns) ) 755 ALLOCATE( vmea(l)%z_ag(1:vmea(l)%ns) ) 642 756 ! 643 757 !-- Based on the flag array store the grid indices which correspond to 644 758 !-- the observation coordinates. 645 DO is = nxl-1, nxr+1 646 DO js = nys-1, nyn+1 647 DO ks = nzb, nzt+1 648 IF ( BTEST( meas_flag(ks,js,is), 0 ) ) THEN 759 ns = 0 760 DO i = nxl-1, nxr+1 761 DO j = nys-1, nyn+1 762 DO k = nzb, nzt+1 763 IF ( BTEST( meas_flag(k,j,i), 0 ) ) THEN 649 764 ns = ns + 1 650 vmea(l)%i(ns) = is 651 vmea(l)%j(ns) = js 652 vmea(l)%k(ns) = ks 765 vmea(l)%i(ns) = i 766 vmea(l)%j(ns) = j 767 vmea(l)%k(ns) = k 768 vmea(l)%z_ag(ns) = zu(k) - & 769 zw(get_topography_top_index_ji( j, i, 's' )) 653 770 ENDIF 654 771 ENDDO 655 772 ENDDO 656 773 ENDDO 774 ! 775 !-- Same for the soil. Based on the flag array, count the number of 776 !-- sampling coordinates in soil. Sample at all soil levels in this case. 777 IF ( vmea(l)%soil_sampling ) THEN 778 DO i = nxl, nxr 779 DO j = nys, nyn 780 IF ( ANY( BTEST( meas_flag(:,j,i), 0 ) ) ) THEN 781 IF ( surf_lsm_h%start_index(j,i) <= & 782 surf_lsm_h%end_index(j,i) ) THEN 783 vmea(l)%ns_soil = vmea(l)%ns_soil + & 784 nzt_soil - nzb_soil + 1 785 ENDIF 786 IF ( surf_usm_h%start_index(j,i) <= & 787 surf_usm_h%end_index(j,i) ) THEN 788 vmea(l)%ns_soil = vmea(l)%ns_soil + & 789 nzt_wall - nzb_wall + 1 790 ENDIF 791 ENDIF 792 ENDDO 793 ENDDO 794 ENDIF 795 ! 796 !-- Allocate index arrays as well as array containing height information 797 !-- for soil. 798 IF ( vmea(l)%soil_sampling ) THEN 799 ALLOCATE( vmea(l)%i_soil(1:vmea(l)%ns_soil) ) 800 ALLOCATE( vmea(l)%j_soil(1:vmea(l)%ns_soil) ) 801 ALLOCATE( vmea(l)%k_soil(1:vmea(l)%ns_soil) ) 802 ALLOCATE( vmea(l)%depth(1:vmea(l)%ns_soil) ) 803 ENDIF 804 ! 805 !-- For soil, store the grid indices. 806 ns = 0 807 IF ( vmea(l)%soil_sampling ) THEN 808 DO i = nxl, nxr 809 DO j = nys, nyn 810 IF ( ANY( BTEST( meas_flag(:,j,i), 0 ) ) ) THEN 811 IF ( surf_lsm_h%start_index(j,i) <= & 812 surf_lsm_h%end_index(j,i) ) THEN 813 m = surf_lsm_h%start_index(j,i) 814 DO k = nzb_soil, nzt_soil 815 ns = ns + 1 816 vmea(l)%i_soil(ns) = i 817 vmea(l)%j_soil(ns) = j 818 vmea(l)%k_soil(ns) = k 819 vmea(l)%depth(ns) = zs(k) 820 ENDDO 821 ENDIF 822 823 IF ( surf_usm_h%start_index(j,i) <= & 824 surf_usm_h%end_index(j,i) ) THEN 825 m = surf_usm_h%start_index(j,i) 826 DO k = nzb_wall, nzt_wall 827 ns = ns + 1 828 vmea(l)%i_soil(ns) = i 829 vmea(l)%j_soil(ns) = j 830 vmea(l)%k_soil(ns) = k 831 vmea(l)%depth(ns) = surf_usm_h%zw(k,m) 832 ENDDO 833 ENDIF 834 ENDIF 835 ENDDO 836 ENDDO 837 ENDIF 838 ! 839 !-- Allocate array to save the sampled values. 840 ALLOCATE( vmea(l)%measured_vars(1:vmea(l)%ns,1:vmea(l)%nvar) ) 657 841 658 ! write(9,*) l, "NS: ", ns 659 ! IF ( ns > 0 ) THEN 660 ! write(9,*) "i ", vmea(l)%i 661 ! write(9,*) "j ", vmea(l)%j 662 ! write(9,*) "k ", vmea(l)%k 663 ! write(9,*) 664 ! ENDIF 665 ! 666 !-- Allocate array to save the sampled values. 667 !-- Todo: Is it better to allocate for all variables at a station 668 !-- and store all the values before writing, or sample the variables 669 !-- directly in the data output? 670 ALLOCATE( vmea(l)%measured_vars(1:vmea(l)%nvar,1:vmea(l)%ns) ) 671 ! 672 !-- Initialize with _FillValue 673 vmea(l)%measured_vars(1:vmea(l)%nvar,1:vmea(l)%ns) = vmea(l)%fillout 842 IF ( vmea(l)%soil_sampling ) & 843 ALLOCATE( vmea(l)%measured_vars_soil(1:vmea(l)%ns_soil, & 844 1:vmea(l)%nvar) ) 845 ! 846 !-- Initialize with _FillValues 847 vmea(l)%measured_vars(1:vmea(l)%ns,1:vmea(l)%nvar) = vmea(l)%fillout 848 IF ( vmea(l)%soil_sampling ) & 849 vmea(l)%measured_vars_soil(1:vmea(l)%ns_soil,1:vmea(l)%nvar) = & 850 vmea(l)%fillout 674 851 ! 675 852 !-- Deallocate temporary coordinate arrays … … 677 854 IF ( ALLOCATED( n_utm ) ) DEALLOCATE( n_utm ) 678 855 IF ( ALLOCATED( z_ag ) ) DEALLOCATE( z_ag ) 856 IF ( ALLOCATED( z_ag ) ) DEALLOCATE( vmea(l)%dim_t ) 679 857 ENDIF 680 858 ENDDO 681 ! flush(9)682 683 859 ! 684 860 !-- Close input file for virtual measurements. Therefore, just call 685 861 !-- the read attribute routine with the "close" option. 686 CALL netcdf_data_input_att( nvm, char_numstations, id_vm, '', & 862 CALL netcdf_data_input_att( vmea_general%nvm, char_numstations, & 863 vmea_general%id_vm, '', & 687 864 global_attribute, 'close', '' ) 865 ! 866 !-- Sum-up the number of observation coordiates, for atmosphere first. 867 !-- This is actually only required for data output. 868 ALLOCATE( ns_all(1:vmea_general%nvm) ) 869 ns_all = 0 870 #if defined( __parallel ) 871 CALL MPI_ALLREDUCE( vmea(:)%ns, ns_all(:), vmea_general%nvm, MPI_INTEGER, & 872 MPI_SUM, comm2d, ierr ) 873 #else 874 ns_all(:) = vmea(:)%ns 875 #endif 876 vmea(:)%ns_tot = ns_all(:) 877 ! 878 !-- Now for soil 879 ns_all = 0 880 #if defined( __parallel ) 881 CALL MPI_ALLREDUCE( vmea(:)%ns_soil, ns_all(:), vmea_general%nvm, & 882 MPI_INTEGER, MPI_SUM, comm2d, ierr ) 883 #else 884 ns_all(:) = vmea(:)%ns_soil 885 #endif 886 vmea(:)%ns_soil_tot = ns_all(:) 887 888 DEALLOCATE( ns_all ) 688 889 ! 689 890 !-- Dellocate flag array 690 891 DEALLOCATE( meas_flag ) 892 ! 893 !-- Initialize binary data output of virtual measurements. 894 !-- Open binary output file. 895 CALL check_open( 27 ) 896 ! 897 !-- Output header information. 898 CALL vm_data_output 691 899 692 900 END SUBROUTINE vm_init … … 696 904 ! Description: 697 905 ! ------------ 906 !> Binary data output. 907 !------------------------------------------------------------------------------! 908 SUBROUTINE vm_data_output 909 910 USE pegrid 911 912 IMPLICIT NONE 913 914 INTEGER(iwp) :: i !< running index over IO blocks 915 INTEGER(iwp) :: l !< running index over all stations 916 INTEGER(iwp) :: n !< running index over all measured variables at a station 917 ! 918 !-- Header output on each PE 919 IF ( init ) THEN 920 921 DO i = 0, io_blocks-1 922 IF ( i == io_group ) THEN 923 WRITE ( 27 ) 'number of measurements ' 924 WRITE ( 27 ) vmea_general%nvm 925 926 DO l = 1, vmea_general%nvm 927 WRITE ( 27 ) 'site ' 928 WRITE ( 27 ) vmea(l)%site 929 WRITE ( 27 ) 'file ' 930 WRITE ( 27 ) vmea(l)%filename_original 931 WRITE ( 27 ) 'feature_type ' 932 WRITE ( 27 ) vmea(l)%feature_type 933 WRITE ( 27 ) 'origin_x_obs ' 934 WRITE ( 27 ) vmea(l)%origin_x_obs 935 WRITE ( 27 ) 'origin_y_obs ' 936 WRITE ( 27 ) vmea(l)%origin_y_obs 937 WRITE ( 27 ) 'total number of observation points' 938 WRITE ( 27 ) vmea(l)%ns_tot 939 WRITE ( 27 ) 'number of measured variables ' 940 WRITE ( 27 ) vmea(l)%nvar 941 WRITE ( 27 ) 'variables ' 942 WRITE ( 27 ) vmea(l)%measured_vars_name(:) 943 WRITE ( 27 ) 'number of observation points ' 944 WRITE ( 27 ) vmea(l)%ns 945 WRITE ( 27 ) 'E_UTM ' 946 WRITE ( 27 ) init_model%origin_x + & 947 REAL( vmea(l)%i(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dx 948 WRITE ( 27 ) 'N_UTM ' 949 WRITE ( 27 ) init_model%origin_y + & 950 REAL( vmea(l)%j(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dy 951 WRITE ( 27 ) 'Z_AG ' 952 WRITE ( 27 ) vmea(l)%z_ag(1:vmea(l)%ns) 953 WRITE ( 27 ) 'soil sampling ' 954 WRITE ( 27 ) MERGE( 'yes ', & 955 'no ', & 956 vmea(l)%soil_sampling ) 957 958 IF ( vmea(l)%soil_sampling ) THEN 959 WRITE ( 27 ) 'total number of soil points ' 960 WRITE ( 27 ) vmea(l)%ns_soil_tot 961 print*, "vmea(l)%ns_soil_tot", vmea(l)%ns_soil_tot 962 WRITE ( 27 ) 'number of soil points ' 963 WRITE ( 27 ) vmea(l)%ns_soil 964 WRITE ( 27 ) 'E_UTM soil ' 965 WRITE ( 27 ) init_model%origin_x + & 966 REAL( vmea(l)%i_soil(1:vmea(l)%ns_soil) + 0.5_wp, & 967 KIND = wp ) * dx 968 WRITE ( 27 ) 'N_UTM soil ' 969 WRITE ( 27 ) init_model%origin_y + & 970 REAL( vmea(l)%j_soil(1:vmea(l)%ns_soil) + 0.5_wp, & 971 KIND = wp ) * dy 972 WRITE ( 27 ) 'DEPTH ' 973 WRITE ( 27 ) vmea(l)%depth(1:vmea(l)%ns_soil) 974 ENDIF 975 ENDDO 976 977 ENDIF 978 ENDDO 979 980 #if defined( __parallel ) 981 CALL MPI_BARRIER( comm2d, ierr ) 982 #endif 983 ! 984 !-- After header information is written, set control flag to .FALSE. 985 init = .FALSE. 986 ! 987 !-- Data output at each measurement timestep on each PE 988 ELSE 989 DO i = 0, io_blocks-1 990 991 IF ( i == io_group ) THEN 992 WRITE( 27 ) 'output time ' 993 WRITE( 27 ) time_since_reference_point 994 DO l = 1, vmea_general%nvm 995 ! 996 !-- Skip binary writing if no observation points are defined on PE 997 IF ( vmea(l)%ns < 1 .AND. vmea(l)%ns_soil < 1) CYCLE 998 DO n = 1, vmea(l)%nvar 999 WRITE( 27 ) vmea(l)%measured_vars_name(n) 1000 IF ( vmea(l)%soil_sampling .AND. & 1001 ANY( TRIM( vmea(l)%measured_vars_name(n)) == & 1002 soil_vars ) ) THEN 1003 WRITE( 27 ) vmea(l)%measured_vars_soil(:,n) 1004 ELSE 1005 WRITE( 27 ) vmea(l)%measured_vars(:,n) 1006 ENDIF 1007 ENDDO 1008 1009 ENDDO 1010 ENDIF 1011 ENDDO 1012 #if defined( __parallel ) 1013 CALL MPI_BARRIER( comm2d, ierr ) 1014 #endif 1015 ENDIF 1016 1017 END SUBROUTINE vm_data_output 1018 1019 1020 !------------------------------------------------------------------------------! 1021 ! Description: 1022 ! ------------ 1023 !> Write end-of-file statement as last action. 1024 !------------------------------------------------------------------------------! 1025 SUBROUTINE vm_last_actions 1026 1027 USE pegrid 1028 1029 IMPLICIT NONE 1030 1031 INTEGER(iwp) :: i !< running index over IO blocks 1032 INTEGER(iwp) :: l !< running index over all stations 1033 INTEGER(iwp) :: n !< running index over all measured variables at a station 1034 1035 DO i = 0, io_blocks-1 1036 IF ( i == io_group ) THEN 1037 WRITE( 27 ) 'EOF ' 1038 ENDIF 1039 ENDDO 1040 #if defined( __parallel ) 1041 CALL MPI_BARRIER( comm2d, ierr ) 1042 #endif 1043 ! 1044 !-- Close binary file 1045 CALL close_file( 27 ) 1046 1047 END SUBROUTINE vm_last_actions 1048 1049 !------------------------------------------------------------------------------! 1050 ! Description: 1051 ! ------------ 698 1052 !> Sampling of the actual quantities along the observation coordinates 699 1053 !------------------------------------------------------------------------------! … … 714 1068 IMPLICIT NONE 715 1069 716 !CHARACTER(LEN=10) :: trimvar !< dummy for the measured variable name1070 CHARACTER(LEN=10) :: trimvar !< dummy for the measured variable name 717 1071 718 INTEGER(iwp) :: i !< grid index in x-direction 719 INTEGER(iwp) :: j !< grid index in y-direction 720 INTEGER(iwp) :: k !< grid index in z-direction 721 INTEGER(iwp) :: l !< running index over the number of stations 722 INTEGER(iwp) :: m !< running index over all virtual observation coordinates 723 INTEGER(iwp) :: mm !< index of surface element which corresponds to the virtual observation coordinate 724 INTEGER(iwp) :: n !< running index over all measured variables at a station 725 INTEGER(iwp) :: nn !< running index over the number of chemcal species 726 ! 727 !-- Loop over all stations. For each possible variable loop over all 728 !-- observation points 729 DO l = 1, nvm 730 ! 731 !-- Loop over all measured variables at this station. Please note, 732 !-- velocity components are interpolated onto scalar grid. 1072 INTEGER(iwp) :: i !< grid index in x-direction 1073 INTEGER(iwp) :: j !< grid index in y-direction 1074 INTEGER(iwp) :: k !< grid index in z-direction 1075 INTEGER(iwp) :: ind_chem !< dummy index to identify chemistry variable and translate it from (UC)2 standard to interal naming 1076 INTEGER(iwp) :: l !< running index over the number of stations 1077 INTEGER(iwp) :: m !< running index over all virtual observation coordinates 1078 INTEGER(iwp) :: mm !< index of surface element which corresponds to the virtual observation coordinate 1079 INTEGER(iwp) :: n !< running index over all measured variables at a station 1080 INTEGER(iwp) :: nn !< running index over the number of chemcal species 1081 1082 LOGICAL :: match_lsm !< flag indicating natural-type surface 1083 LOGICAL :: match_usm !< flag indicating urban-type surface 1084 ! 1085 !-- Loop over all sites. 1086 DO l = 1, vmea_general%nvm 1087 ! 1088 !-- At the beginning, set _FillValues 1089 IF ( ALLOCATED( vmea(l)%measured_vars ) ) & 1090 vmea(l)%measured_vars = vmea(l)%fillout 1091 IF ( ALLOCATED( vmea(l)%measured_vars_soil ) ) & 1092 vmea(l)%measured_vars_soil = vmea(l)%fillout 1093 ! 1094 !-- Loop over all variables measured at this site. 733 1095 DO n = 1, vmea(l)%nvar 734 1096 … … 741 1103 j = vmea(l)%j(m) 742 1104 i = vmea(l)%i(m) 743 vmea(l)%measured_vars( n,m) = pt(k,j,i)1105 vmea(l)%measured_vars(m,n) = pt(k,j,i) 744 1106 ENDDO 745 1107 ENDIF 746 1108 747 CASE ( 'ta' , 't_va')1109 CASE ( 'ta' ) 748 1110 IF ( .NOT. neutral ) THEN 749 1111 DO m = 1, vmea(l)%ns … … 751 1113 j = vmea(l)%j(m) 752 1114 i = vmea(l)%i(m) 753 vmea(l)%measured_vars( n,m) = pt(k,j,i) * exner( k )1115 vmea(l)%measured_vars(m,n) = pt(k,j,i) * exner( k ) 754 1116 ENDDO 755 1117 ENDIF 1118 1119 CASE ( 't_va' ) 756 1120 757 1121 CASE ( 'hus', 'haa' ) … … 761 1125 j = vmea(l)%j(m) 762 1126 i = vmea(l)%i(m) 763 vmea(l)%measured_vars( n,m) = q(k,j,i)1127 vmea(l)%measured_vars(m,n) = q(k,j,i) 764 1128 ENDDO 765 1129 ENDIF … … 770 1134 j = vmea(l)%j(m) 771 1135 i = vmea(l)%i(m) 772 vmea(l)%measured_vars( n,m) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )1136 vmea(l)%measured_vars(m,n) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) 773 1137 ENDDO 774 1138 … … 778 1142 j = vmea(l)%j(m) 779 1143 i = vmea(l)%i(m) 780 vmea(l)%measured_vars( n,m) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )1144 vmea(l)%measured_vars(m,n) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) 781 1145 ENDDO 782 1146 783 1147 CASE ( 'w' ) 784 1148 DO m = 1, vmea(l)%ns 785 k = vmea(l)%k(m)1149 k = MAX ( 1, vmea(l)%k(m) ) 786 1150 j = vmea(l)%j(m) 787 1151 i = vmea(l)%i(m) 788 vmea(l)%measured_vars( n,m) = w(k,j,i) !0.5_wp * ( w(k,j,i) + w(k-1,j,i) )1152 vmea(l)%measured_vars(m,n) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) ) 789 1153 ENDDO 790 1154 … … 794 1158 j = vmea(l)%j(m) 795 1159 i = vmea(l)%i(m) 796 vmea(l)%measured_vars( n,m) = SQRT( &1160 vmea(l)%measured_vars(m,n) = SQRT( & 797 1161 ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2 + & 798 1162 ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2 & … … 806 1170 i = vmea(l)%i(m) 807 1171 808 vmea(l)%measured_vars( n,m) = ATAN2( &1172 vmea(l)%measured_vars(m,n) = ATAN2( & 809 1173 - 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ), & 810 1174 - 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) & 811 1175 ) * 180.0_wp / pi 812 1176 ENDDO 813 ! 814 !-- MS: list of variables needs extension. 815 CASE ( 'mcpm1', 'mcpm2p5', 'mcpm10', 'mcco', 'mcco2', 'mcbcda', & 816 'ncaa', 'mfco2', 'mfco', 'mfch4', 'mfnh3', 'mfno', & 817 'mfno2', 'mfso2', 'mfh20', 'tr03' ) 818 IF ( air_chemistry ) THEN 1177 1178 CASE ( 'utheta' ) 1179 DO m = 1, vmea(l)%ns 1180 k = vmea(l)%k(m) 1181 j = vmea(l)%j(m) 1182 i = vmea(l)%i(m) 1183 vmea(l)%measured_vars(m,n) = 0.5_wp * & 1184 ( u(k,j,i) + u(k,j,i+1) ) * & 1185 pt(k,j,i) 1186 ENDDO 1187 1188 CASE ( 'vtheta' ) 1189 DO m = 1, vmea(l)%ns 1190 k = vmea(l)%k(m) 1191 j = vmea(l)%j(m) 1192 i = vmea(l)%i(m) 1193 vmea(l)%measured_vars(m,n) = 0.5_wp * & 1194 ( v(k,j,i) + v(k,j+1,i) ) * & 1195 pt(k,j,i) 1196 ENDDO 1197 1198 CASE ( 'wtheta' ) 1199 DO m = 1, vmea(l)%ns 1200 k = MAX ( 1, vmea(l)%k(m) ) 1201 j = vmea(l)%j(m) 1202 i = vmea(l)%i(m) 1203 vmea(l)%measured_vars(m,n) = 0.5_wp * & 1204 ( w(k-1,j,i) + w(k,j,i) ) * & 1205 pt(k,j,i) 1206 ENDDO 1207 1208 CASE ( 'uw' ) 1209 DO m = 1, vmea(l)%ns 1210 k = MAX ( 1, vmea(l)%k(m) ) 1211 j = vmea(l)%j(m) 1212 i = vmea(l)%i(m) 1213 vmea(l)%measured_vars(m,n) = 0.25_wp * & 1214 ( w(k-1,j,i) + w(k,j,i) ) * & 1215 ( u(k,j,i) + u(k,j,i+1) ) 1216 ENDDO 1217 1218 CASE ( 'vw' ) 1219 DO m = 1, vmea(l)%ns 1220 k = MAX ( 1, vmea(l)%k(m) ) 1221 j = vmea(l)%j(m) 1222 i = vmea(l)%i(m) 1223 vmea(l)%measured_vars(m,n) = 0.25_wp * & 1224 ( w(k-1,j,i) + w(k,j,i) ) * & 1225 ( v(k,j,i) + v(k,j+1,i) ) 1226 ENDDO 1227 1228 CASE ( 'uv' ) 1229 DO m = 1, vmea(l)%ns 1230 k = MAX ( 1, vmea(l)%k(m) ) 1231 j = vmea(l)%j(m) 1232 i = vmea(l)%i(m) 1233 vmea(l)%measured_vars(m,n) = 0.25_wp * & 1234 ( u(k,j,i) + u(k,j,i+1) ) * & 1235 ( v(k,j,i) + v(k,j+1,i) ) 1236 ENDDO 1237 ! 1238 !-- List of variables may need extension. 1239 CASE ( 'mcpm1', 'mcpm2p5', 'mcpm10', 'mfco', 'mfno', 'mfno2', & 1240 'tro3' ) 1241 IF ( air_chemistry ) THEN 1242 ! 1243 !-- First, search for the measured variable in the chem_vars 1244 !-- list, in order to get the internal name of the variable. 1245 DO nn = 1, UBOUND( chem_vars, 2 ) 1246 IF ( TRIM( vmea(l)%measured_vars_name(m) ) == & 1247 TRIM( chem_vars(0,nn) ) ) ind_chem = nn 1248 ENDDO 1249 ! 1250 !-- Run loop over all chemical species, if the measured 1251 !-- variable matches the interal name, sample the variable. 819 1252 DO nn = 1, nspec 820 IF ( TRIM( vmea(l)%measured_vars_name(m) ) ==&1253 IF ( TRIM( chem_vars(1,ind_chem) ) == & 821 1254 TRIM( chem_species(nn)%name ) ) THEN 822 1255 DO m = 1, vmea(l)%ns … … 824 1257 j = vmea(l)%j(m) 825 1258 i = vmea(l)%i(m) 826 vmea(l)%measured_vars( n,m) = &1259 vmea(l)%measured_vars(m,n) = & 827 1260 chem_species(nn)%conc(k,j,i) 828 1261 ENDDO … … 843 1276 DO mm = surf_def_h(0)%start_index(j,i), & 844 1277 surf_def_h(0)%end_index(j,i) 845 vmea(l)%measured_vars( n,m) = surf_def_h(0)%us(mm)1278 vmea(l)%measured_vars(m,n) = surf_def_h(0)%us(mm) 846 1279 ENDDO 847 1280 DO mm = surf_lsm_h%start_index(j,i), & 848 1281 surf_lsm_h%end_index(j,i) 849 vmea(l)%measured_vars( n,m) = surf_lsm_h%us(mm)1282 vmea(l)%measured_vars(m,n) = surf_lsm_h%us(mm) 850 1283 ENDDO 851 1284 DO mm = surf_usm_h%start_index(j,i), & 852 1285 surf_usm_h%end_index(j,i) 853 vmea(l)%measured_vars( n,m) = surf_usm_h%us(mm)1286 vmea(l)%measured_vars(m,n) = surf_usm_h%us(mm) 854 1287 ENDDO 855 1288 ENDDO … … 867 1300 DO mm = surf_def_h(0)%start_index(j,i), & 868 1301 surf_def_h(0)%end_index(j,i) 869 vmea(l)%measured_vars( n,m) = surf_def_h(0)%ts(mm)1302 vmea(l)%measured_vars(m,n) = surf_def_h(0)%ts(mm) 870 1303 ENDDO 871 1304 DO mm = surf_lsm_h%start_index(j,i), & 872 1305 surf_lsm_h%end_index(j,i) 873 vmea(l)%measured_vars( n,m) = surf_lsm_h%ts(mm)1306 vmea(l)%measured_vars(m,n) = surf_lsm_h%ts(mm) 874 1307 ENDDO 875 1308 DO mm = surf_usm_h%start_index(j,i), & 876 1309 surf_usm_h%end_index(j,i) 877 vmea(l)%measured_vars( n,m) = surf_usm_h%ts(mm)1310 vmea(l)%measured_vars(m,n) = surf_usm_h%ts(mm) 878 1311 ENDDO 879 1312 ENDDO … … 891 1324 DO mm = surf_def_h(0)%start_index(j,i), & 892 1325 surf_def_h(0)%end_index(j,i) 893 vmea(l)%measured_vars( n,m) = surf_def_h(0)%qsws(mm)1326 vmea(l)%measured_vars(m,n) = surf_def_h(0)%qsws(mm) 894 1327 ENDDO 895 1328 DO mm = surf_lsm_h%start_index(j,i), & 896 1329 surf_lsm_h%end_index(j,i) 897 vmea(l)%measured_vars( n,m) = surf_lsm_h%qsws(mm)1330 vmea(l)%measured_vars(m,n) = surf_lsm_h%qsws(mm) 898 1331 ENDDO 899 1332 DO mm = surf_usm_h%start_index(j,i), & 900 1333 surf_usm_h%end_index(j,i) 901 vmea(l)%measured_vars( n,m) = surf_usm_h%qsws(mm)1334 vmea(l)%measured_vars(m,n) = surf_usm_h%qsws(mm) 902 1335 ENDDO 903 1336 ENDDO … … 915 1348 DO mm = surf_def_h(0)%start_index(j,i), & 916 1349 surf_def_h(0)%end_index(j,i) 917 vmea(l)%measured_vars( n,m) = surf_def_h(0)%shf(mm)1350 vmea(l)%measured_vars(m,n) = surf_def_h(0)%shf(mm) 918 1351 ENDDO 919 1352 DO mm = surf_lsm_h%start_index(j,i), & 920 1353 surf_lsm_h%end_index(j,i) 921 vmea(l)%measured_vars( n,m) = surf_lsm_h%shf(mm)1354 vmea(l)%measured_vars(m,n) = surf_lsm_h%shf(mm) 922 1355 ENDDO 923 1356 DO mm = surf_usm_h%start_index(j,i), & 924 1357 surf_usm_h%end_index(j,i) 925 vmea(l)%measured_vars( n,m) = surf_usm_h%shf(mm)1358 vmea(l)%measured_vars(m,n) = surf_usm_h%shf(mm) 926 1359 ENDDO 927 1360 ENDDO … … 940 1373 DO mm = surf_lsm_h%start_index(j,i), & 941 1374 surf_lsm_h%end_index(j,i) 942 vmea(l)%measured_vars( n,m) = surf_lsm_h%rad_net(mm)1375 vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_net(mm) 943 1376 ENDDO 944 1377 DO mm = surf_usm_h%start_index(j,i), & 945 1378 surf_usm_h%end_index(j,i) 946 vmea(l)%measured_vars( n,m) = surf_usm_h%rad_net(mm)1379 vmea(l)%measured_vars(m,n) = surf_usm_h%rad_net(mm) 947 1380 ENDDO 948 1381 ENDDO 949 1382 ENDIF 950 1383 951 CASE ( 'rsus' , 'rsu')1384 CASE ( 'rsus' ) 952 1385 IF ( radiation ) THEN 953 1386 DO m = 1, vmea(l)%ns … … 962 1395 DO mm = surf_lsm_h%start_index(j,i), & 963 1396 surf_lsm_h%end_index(j,i) 964 vmea(l)%measured_vars( n,m) = surf_lsm_h%rad_sw_out(mm)1397 vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_sw_out(mm) 965 1398 ENDDO 966 1399 DO mm = surf_usm_h%start_index(j,i), & 967 1400 surf_usm_h%end_index(j,i) 968 vmea(l)%measured_vars( n,m) = surf_usm_h%rad_sw_out(mm)1401 vmea(l)%measured_vars(m,n) = surf_usm_h%rad_sw_out(mm) 969 1402 ENDDO 970 1403 ENDDO 971 1404 ENDIF 972 1405 973 CASE ( 'rsds' , 'rsd')1406 CASE ( 'rsds' ) 974 1407 IF ( radiation ) THEN 975 1408 DO m = 1, vmea(l)%ns … … 984 1417 DO mm = surf_lsm_h%start_index(j,i), & 985 1418 surf_lsm_h%end_index(j,i) 986 vmea(l)%measured_vars( n,m) = surf_lsm_h%rad_sw_in(mm)1419 vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_sw_in(mm) 987 1420 ENDDO 988 1421 DO mm = surf_usm_h%start_index(j,i), & 989 1422 surf_usm_h%end_index(j,i) 990 vmea(l)%measured_vars( n,m) = surf_usm_h%rad_sw_in(mm)1423 vmea(l)%measured_vars(m,n) = surf_usm_h%rad_sw_in(mm) 991 1424 ENDDO 992 1425 ENDDO 993 1426 ENDIF 994 1427 995 CASE ( 'rlus' , 'rlu')1428 CASE ( 'rlus' ) 996 1429 IF ( radiation ) THEN 997 1430 DO m = 1, vmea(l)%ns … … 1006 1439 DO mm = surf_lsm_h%start_index(j,i), & 1007 1440 surf_lsm_h%end_index(j,i) 1008 vmea(l)%measured_vars( n,m) = surf_lsm_h%rad_lw_out(mm)1441 vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_lw_out(mm) 1009 1442 ENDDO 1010 1443 DO mm = surf_usm_h%start_index(j,i), & 1011 1444 surf_usm_h%end_index(j,i) 1012 vmea(l)%measured_vars( n,m) = surf_usm_h%rad_lw_out(mm)1445 vmea(l)%measured_vars(m,n) = surf_usm_h%rad_lw_out(mm) 1013 1446 ENDDO 1014 1447 ENDDO 1015 1448 ENDIF 1016 1449 1017 CASE ( 'rlds' , 'rld')1450 CASE ( 'rlds' ) 1018 1451 IF ( radiation ) THEN 1019 1452 DO m = 1, vmea(l)%ns … … 1028 1461 DO mm = surf_lsm_h%start_index(j,i), & 1029 1462 surf_lsm_h%end_index(j,i) 1030 vmea(l)%measured_vars( n,m) = surf_lsm_h%rad_lw_in(mm)1463 vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_lw_in(mm) 1031 1464 ENDDO 1032 1465 DO mm = surf_usm_h%start_index(j,i), & 1033 1466 surf_usm_h%end_index(j,i) 1034 vmea(l)%measured_vars( n,m) = surf_usm_h%rad_lw_in(mm)1467 vmea(l)%measured_vars(m,n) = surf_usm_h%rad_lw_in(mm) 1035 1468 ENDDO 1036 1469 ENDDO 1037 1470 ENDIF 1471 1472 CASE ( 'rsd' ) 1473 IF ( radiation ) THEN 1474 DO m = 1, vmea(l)%ns 1475 k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' ) 1476 j = vmea(l)%j(m) 1477 i = vmea(l)%i(m) 1478 1479 vmea(l)%measured_vars(m,n) = rad_sw_in(k,j,i) 1480 ENDDO 1481 ENDIF 1482 1483 CASE ( 'rsu' ) 1484 IF ( radiation ) THEN 1485 DO m = 1, vmea(l)%ns 1486 k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' ) 1487 j = vmea(l)%j(m) 1488 i = vmea(l)%i(m) 1489 1490 vmea(l)%measured_vars(m,n) = rad_sw_out(k,j,i) 1491 ENDDO 1492 ENDIF 1493 1494 CASE ( 'rlu' ) 1495 IF ( radiation ) THEN 1496 DO m = 1, vmea(l)%ns 1497 k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' ) 1498 j = vmea(l)%j(m) 1499 i = vmea(l)%i(m) 1500 1501 vmea(l)%measured_vars(m,n) = rad_lw_out(k,j,i) 1502 ENDDO 1503 ENDIF 1504 1505 CASE ( 'rld' ) 1506 IF ( radiation ) THEN 1507 DO m = 1, vmea(l)%ns 1508 k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' ) 1509 j = vmea(l)%j(m) 1510 i = vmea(l)%i(m) 1511 1512 vmea(l)%measured_vars(m,n) = rad_lw_in(k,j,i) 1513 ENDDO 1514 ENDIF 1515 1516 CASE ( 'rsddif' ) 1517 IF ( radiation ) THEN 1518 DO m = 1, vmea(l)%ns 1519 j = vmea(l)%j(m) 1520 i = vmea(l)%i(m) 1521 1522 vmea(l)%measured_vars(m,n) = rad_sw_in_diff(j,i) 1523 ENDDO 1524 ENDIF 1525 1526 CASE ( 't_soil' ) 1527 DO m = 1, vmea(l)%ns_soil 1528 i = vmea(l)%i_soil(m) 1529 j = vmea(l)%j_soil(m) 1530 k = vmea(l)%k_soil(m) 1531 1532 match_lsm = surf_lsm_h%start_index(j,i) <= & 1533 surf_lsm_h%end_index(j,i) 1534 match_usm = surf_usm_h%start_index(j,i) <= & 1535 surf_usm_h%end_index(j,i) 1536 1537 IF ( match_lsm ) THEN 1538 mm = surf_lsm_h%start_index(j,i) 1539 vmea(l)%measured_vars_soil(m,n) = t_soil_h%var_2d(k,mm) 1540 ENDIF 1541 1542 IF ( match_usm ) THEN 1543 mm = surf_usm_h%start_index(j,i) 1544 vmea(l)%measured_vars_soil(m,n) = t_wall_h(k,mm) 1545 ENDIF 1546 ENDDO 1547 1548 CASE ( 'm_soil' ) 1549 DO m = 1, vmea(l)%ns_soil 1550 i = vmea(l)%i_soil(m) 1551 j = vmea(l)%j_soil(m) 1552 k = vmea(l)%k_soil(m) 1553 1554 match_lsm = surf_lsm_h%start_index(j,i) <= & 1555 surf_lsm_h%end_index(j,i) 1556 1557 IF ( match_lsm ) THEN 1558 mm = surf_lsm_h%start_index(j,i) 1559 vmea(l)%measured_vars_soil(m,n) = m_soil_h%var_2d(k,mm) 1560 ENDIF 1561 1562 ENDDO 1038 1563 ! 1039 1564 !-- More will follow ... 1040 1565 1566 ! 1567 !-- No match found - just set a fill value 1568 CASE DEFAULT 1569 vmea(l)%measured_vars(:,n) = vmea(l)%fillout 1041 1570 END SELECT 1042 1571 … … 1044 1573 1045 1574 ENDDO 1046 1575 1047 1576 END SUBROUTINE vm_sampling 1048 1577
Note: See TracChangeset
for help on using the changeset viewer.
