!> @virtual_measurement_mod.f90 !------------------------------------------------------------------------------! ! This file is part of the PALM model system. ! ! PALM is free software: you can redistribute it and/or modify it under the ! terms of the GNU General Public License as published by the Free Software ! Foundation, either version 3 of the License, or (at your option) any later ! version. ! ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR ! A PARTICULAR PURPOSE. See the GNU General Public License for more details. ! ! You should have received a copy of the GNU General Public License along with ! PALM. If not, see . ! ! Copyright 2017 Leibniz Universitaet Hannover !------------------------------------------------------------------------------! ! ! Current revisions: ! ----------------- ! ! ! Former revisions: ! ----------------- ! $Id: virtual_measurement_mod.f90 4182 2019-08-22 15:20:23Z suehring $ ! Corrected "Former revisions" section ! ! 4168 2019-08-16 13:50:17Z suehring ! Replace function get_topography_top_index by topo_top_ind ! ! 3988 2019-05-22 11:32:37Z kanani ! Add variables to enable steering of output interval for virtual measurements ! ! 3913 2019-04-17 15:12:28Z gronemeier ! Bugfix: rotate positions of measurements before writing them into file ! ! 3910 2019-04-17 11:46:56Z suehring ! Bugfix in rotation of UTM coordinates ! ! 3904 2019-04-16 18:22:51Z gronemeier ! Rotate coordinates of stations by given rotation_angle ! ! 3876 2019-04-08 18:41:49Z knoop ! Remove print statement ! ! 3854 2019-04-02 16:59:33Z suehring ! renamed nvar to nmeas, replaced USE chem_modules by USE chem_gasphase_mod and ! nspec by nvar ! ! 3766 2019-02-26 16:23:41Z raasch ! unused variables removed ! ! 3718 2019-02-06 11:08:28Z suehring ! Adjust variable name connections between UC2 and chemistry variables ! ! 3717 2019-02-05 17:21:16Z suehring ! Additional check + error numbers adjusted ! ! 3706 2019-01-29 20:02:26Z suehring ! unused variables removed ! ! 3705 2019-01-29 19:56:39Z suehring ! - initialization revised ! - binary data output ! - list of allowed variables extended ! ! 3704 2019-01-29 19:51:41Z suehring ! Sampling of variables ! ! 3473 2018-10-30 20:50:15Z suehring ! Initial revision ! ! Authors: ! -------- ! @author Matthias Suehring ! ! Description: ! ------------ !> The module acts as an interface between 'real-world' observations and !> model simulations. Virtual measurements will be taken in the model at the !> coordinates representative for the 'real-world' observation coordinates. !> More precisely, coordinates and measured quanties will be read from a !> NetCDF file which contains all required information. In the model, !> the same quantities (as long as all the required components are switched-on) !> will be sampled at the respective positions and output into an extra file, !> which allows for straight-forward comparison of model results with !> observations. !> !> @todo list_of_allowed variables needs careful checking !> @todo Check if sign of surface fluxes for heat, radiation, etc., follows !> the (UC)2 standard !> @note Fluxes are not processed !------------------------------------------------------------------------------! MODULE virtual_measurement_mod USE arrays_3d, & ONLY: q, pt, u, v, w, zu, zw USE basic_constants_and_equations_mod, & ONLY: pi USE chem_gasphase_mod, & ONLY: nvar USE chem_modules, & ONLY: chem_species USE control_parameters, & ONLY: air_chemistry, dz, humidity, io_blocks, io_group, neutral, & message_string, time_since_reference_point, virtual_measurement USE cpulog, & ONLY: cpu_log, log_point USE grid_variables, & ONLY: dx, dy USE indices, & ONLY: nzb, nzt, nxl, nxr, nys, nyn, nx, ny, topo_top_ind, wall_flags_0 USE kinds USE netcdf_data_input_mod, & ONLY: init_model USE pegrid USE surface_mod, & ONLY: surf_lsm_h, surf_usm_h USE land_surface_model_mod, & ONLY: nzb_soil, nzs, nzt_soil, zs, t_soil_h, m_soil_h USE radiation_model_mod USE urban_surface_mod, & ONLY: nzb_wall, nzt_wall, t_wall_h IMPLICIT NONE TYPE virt_general INTEGER(iwp) :: id_vm !< NetCDF file id for virtual measurements INTEGER(iwp) :: nvm = 0 !< number of virtual measurements END TYPE virt_general TYPE virt_mea CHARACTER(LEN=100) :: feature_type !< type of the measurement CHARACTER(LEN=100) :: filename_original !< name of the original file CHARACTER(LEN=100) :: site !< name of the measurement site CHARACTER(LEN=10), DIMENSION(:), ALLOCATABLE :: measured_vars_name !< name of the measured variables INTEGER(iwp) :: ns = 0 !< number of observation coordinates on subdomain, for atmospheric measurements INTEGER(iwp) :: ns_tot = 0 !< total number of observation coordinates, for atmospheric measurements INTEGER(iwp) :: ntraj !< number of trajectories of a measurement INTEGER(iwp) :: nmeas !< number of measured variables (atmosphere + soil) INTEGER(iwp) :: ns_soil = 0 !< number of observation coordinates on subdomain, for soil measurements INTEGER(iwp) :: ns_soil_tot = 0 !< total number of observation coordinates, for soil measurements INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: dim_t !< number observations individual for each trajectory or station that are no _FillValues INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: i !< grid index for measurement position in x-direction INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: j !< grid index for measurement position in y-direction INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k !< grid index for measurement position in k-direction INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: i_soil !< grid index for measurement position in x-direction INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: j_soil !< grid index for measurement position in y-direction INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k_soil !< grid index for measurement position in k-direction LOGICAL :: trajectory = .FALSE. !< flag indicating that the observation is a mobile observation LOGICAL :: timseries = .FALSE. !< flag indicating that the observation is a stationary point measurement LOGICAL :: timseries_profile = .FALSE. !< flag indicating that the observation is a stationary profile measurement LOGICAL :: soil_sampling = .FALSE. !< flag indicating that soil state variables were sampled REAL(wp) :: fill_eutm !< fill value for UTM coordinates in case of missing values REAL(wp) :: fill_nutm !< fill value for UTM coordinates in case of missing values REAL(wp) :: fill_zag !< fill value for heigth coordinates in case of missing values REAL(wp) :: fillout = -999.9 !< fill value for output in case a observation is taken from inside a building REAL(wp) :: origin_x_obs !< origin of the observation in UTM coordiates in x-direction REAL(wp) :: origin_y_obs !< origin of the observation in UTM coordiates in y-direction REAL(wp), DIMENSION(:), ALLOCATABLE :: z_ag !< measurement height above ground level REAL(wp), DIMENSION(:), ALLOCATABLE :: depth !< measurement depth in soil REAL(wp), DIMENSION(:,:), ALLOCATABLE :: measured_vars !< measured variables REAL(wp), DIMENSION(:,:), ALLOCATABLE :: measured_vars_soil !< measured variables END TYPE virt_mea CHARACTER(LEN=5) :: char_eutm = "E_UTM" !< dimension name for UTM coordinate easting CHARACTER(LEN=11) :: char_feature = "featureType" !< attribute name for feature type ! This need to generalized CHARACTER(LEN=8) :: char_filename = "filename" !< attribute name for filename CHARACTER(LEN=11) :: char_soil = "soil_sample" !< attribute name for soil sampling indication CHARACTER(LEN=10) :: char_fillvalue = "_FillValue" !< variable attribute name for _FillValue CHARACTER(LEN=18) :: char_mv = "measured_variables" !< variable name for the array with the measured variable names CHARACTER(LEN=5) :: char_nutm = "N_UTM" !< dimension name for UTM coordinate northing CHARACTER(LEN=18) :: char_numstations = "number_of_stations" !< attribute name for number of stations CHARACTER(LEN=8) :: char_origx = "origin_x" !< attribute name for station coordinate in x CHARACTER(LEN=8) :: char_origy = "origin_y" !< attribute name for station coordinate in y CHARACTER(LEN=4) :: char_site = "site" !< attribute name for site name CHARACTER(LEN=19) :: char_zag = "height_above_ground" !< attribute name for height above ground variable CHARACTER(LEN=10) :: type_ts = 'timeSeries' !< name of stationary point measurements CHARACTER(LEN=10) :: type_traj = 'trajectory' !< name of line measurements CHARACTER(LEN=17) :: type_tspr = 'timeSeriesProfile' !< name of stationary profile measurements CHARACTER(LEN=6), DIMENSION(1:5) :: soil_vars = (/ & !< list of soil variables 't_soil', & 'm_soil', & 'lwc ', & 'lwcs ', & 'smp ' /) CHARACTER(LEN=10), DIMENSION(0:1,1:8) :: chem_vars = RESHAPE( (/ & 'mcpm1 ', 'PM1 ', & 'mcpm2p5 ', 'PM2.5 ', & 'mcpm25 ', 'PM25 ', & 'mcpm10 ', 'PM10 ', & 'mfno2 ', 'NO2 ', & 'mfno ', 'NO ', & 'tro3 ', 'O3 ', & 'mfco ', 'CO ' & /), (/ 2, 8 /) ) ! !-- MS: List requires careful revision! CHARACTER(LEN=10), DIMENSION(1:54), PARAMETER :: list_allowed_variables = & !< variables that can be sampled in PALM (/ 'hfls ', & ! surface latent heat flux (W/m2) 'hfss ', & ! surface sensible heat flux (W/m2) 'hur ', & ! relative humidity (-) 'hus ', & ! specific humidity (g/kg) 'haa ', & ! absolute atmospheric humidity (kg/m3) 'mcpm1 ', & ! mass concentration of PM1 (kg/m3) 'mcpm2p5 ', & ! mass concentration of PM2.5 (kg/m3) 'mcpm10 ', & ! mass concentration of PM10 (kg/m3) 'mcco ', & ! mass concentration of CO (kg/m3) 'mcco2 ', & ! mass concentration of CO2 (kg/m3) 'mcbcda ', & ! mass concentration of black carbon paritcles (kg/m3) 'ncaa ', & ! number concentation of particles (1/m3) 'mfco ', & ! mole fraction of CO (mol/mol) 'mfco2 ', & ! mole fraction of CO2 (mol/mol) 'mfch4 ', & ! mole fraction of methane (mol/mol) 'mfnh3 ', & ! mole fraction of amonia (mol/mol) 'mfno ', & ! mole fraction of nitrogen monoxide (mol/mol) 'mfno2 ', & ! mole fraction of nitrogen dioxide (mol/mol) 'mfso2 ', & ! mole fraction of sulfur dioxide (mol/mol) 'mfh20 ', & ! mole fraction of water (mol/mol) 'plev ', & ! ? air pressure - hydrostaic + perturbation? 'rlds ', & ! surface downward longwave flux (W/m2) 'rlus ', & ! surface upward longwave flux (W/m2) 'rsds ', & ! surface downward shortwave flux (W/m2) 'rsus ', & ! surface upward shortwave flux (W/m2) 'ta ', & ! air temperature (degree C) 't_va ', & ! virtual accoustic temperature (K) 'theta ', & ! potential temperature (K) 'tro3 ', & ! mole fraction of ozone air (mol/mol) 'ts ', & ! scaling parameter of temperature (K) 'wspeed ', & ! ? wind speed - horizontal? 'wdir ', & ! wind direction 'us ', & ! friction velocity 'msoil ', & ! ? soil moisture - which depth? 'tsoil ', & ! ? soil temperature - which depth? 'u ', & ! u-component 'utheta ', & ! total eastward kinematic heat flux 'ua ', & ! eastward wind (is there any difference to u?) 'v ', & ! v-component 'vtheta ', & ! total northward kinematic heat flux 'va ', & ! northward wind (is there any difference to v?) 'w ', & ! w-component 'wtheta ', & ! total vertical kinematic heat flux 'rld ', & ! downward longwave radiative flux (W/m2) 'rlu ', & ! upnward longwave radiative flux (W/m2) 'rsd ', & ! downward shortwave radiative flux (W/m2) 'rsu ', & ! upward shortwave radiative flux (W/m2) 'rsddif ', & ! downward shortwave diffuse radiative flux (W/m2) 'rnds ', & ! surface net downward radiative flux (W/m2) 't_soil ', & 'm_soil ', & 'lwc ', & 'lwcs ', & 'smp ' & /) LOGICAL :: global_attribute = .TRUE. !< flag indicating a global attribute LOGICAL :: init = .TRUE. !< flag indicating initialization of data output LOGICAL :: use_virtual_measurement = .FALSE. !< Namelist parameter REAL(wp) :: dt_virtual_measurement = 0.0_wp !< namelist parameter REAL(wp) :: time_virtual_measurement = 0.0_wp !< time since last 3d output REAL(wp) :: vm_time_start = 0.0 !< time after virtual measurements should start TYPE( virt_general ) :: vmea_general !< data structure which encompass general variables TYPE( virt_mea ), DIMENSION(:), ALLOCATABLE :: vmea !< virtual measurement data structure INTERFACE vm_check_parameters MODULE PROCEDURE vm_check_parameters END INTERFACE vm_check_parameters INTERFACE vm_data_output MODULE PROCEDURE vm_data_output END INTERFACE vm_data_output INTERFACE vm_init MODULE PROCEDURE vm_init END INTERFACE vm_init INTERFACE vm_last_actions MODULE PROCEDURE vm_last_actions END INTERFACE vm_last_actions INTERFACE vm_parin MODULE PROCEDURE vm_parin END INTERFACE vm_parin INTERFACE vm_sampling MODULE PROCEDURE vm_sampling END INTERFACE vm_sampling SAVE PRIVATE ! !-- Public interfaces PUBLIC vm_check_parameters, vm_data_output, vm_init, vm_last_actions, & vm_parin, vm_sampling ! !-- Public variables PUBLIC dt_virtual_measurement, time_virtual_measurement, vmea, vmea_general, vm_time_start CONTAINS !------------------------------------------------------------------------------! ! Description: ! ------------ !> Check parameters for virtual measurement module !------------------------------------------------------------------------------! SUBROUTINE vm_check_parameters USE control_parameters, & ONLY: message_string, virtual_measurement USE netcdf_data_input_mod, & ONLY: input_pids_static, input_pids_vm IMPLICIT NONE ! !-- Virtual measurements require a setup file. IF ( virtual_measurement .AND. .NOT. input_pids_vm ) THEN message_string = 'If virtual measurements are taken, a setup input ' // & 'file for the site locations is mandatory.' CALL message( 'vm_check_parameters', 'PA0533', 1, 2, 0, 6, 0 ) ENDIF ! !-- In case virtual measurements are taken, a static input file is required. !-- This is because UTM coordinates for the PALM domain origin are required !-- for correct mapping of the measurements. !-- ToDo: Revise this later and remove this requirement. IF ( virtual_measurement .AND. .NOT. input_pids_static ) THEN message_string = 'If virtual measurements are taken, a static input ' //& 'file is mandatory.' CALL message( 'vm_check_parameters', 'PA0534', 1, 2, 0, 6, 0 ) ENDIF END SUBROUTINE vm_check_parameters !------------------------------------------------------------------------------! ! Description: ! ------------ !> Read namelist for the virtual measurement module !------------------------------------------------------------------------------! SUBROUTINE vm_parin CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file NAMELIST /virtual_measurement_parameters/ dt_virtual_measurement, & use_virtual_measurement, & vm_time_start line = ' ' ! !-- Try to find stg package REWIND ( 11 ) line = ' ' DO WHILE ( INDEX( line, '&virtual_measurement_parameters' ) == 0 ) READ ( 11, '(A)', END=20 ) line ENDDO BACKSPACE ( 11 ) ! !-- Read namelist READ ( 11, virtual_measurement_parameters, ERR = 10, END = 20 ) ! !-- Set flag that indicates that the virtual measurement module is switched on IF ( use_virtual_measurement ) virtual_measurement = .TRUE. GOTO 20 10 BACKSPACE( 11 ) READ( 11 , '(A)') line CALL parin_fail_message( 'virtual_measurement_parameters', line ) 20 CONTINUE END SUBROUTINE vm_parin !------------------------------------------------------------------------------! ! Description: ! ------------ !> Initialize virtual measurements: read coordiante arrays and measured !> variables, set indicies indicating the measurement points, read further !> attributes, etc.. !------------------------------------------------------------------------------! SUBROUTINE vm_init USE arrays_3d, & ONLY: zu, zw USE grid_variables, & ONLY: ddx, ddy, dx, dy USE indices, & ONLY: nxl, nxr, nyn, nys USE netcdf_data_input_mod, & ONLY: init_model, input_file_vm, & netcdf_data_input_get_dimension_length, & netcdf_data_input_att, netcdf_data_input_var IMPLICIT NONE CHARACTER(LEN=5) :: dum !< dummy string indicate station id CHARACTER(LEN=10), DIMENSION(50) :: measured_variables_file = '' !< array with all measured variables read from NetCDF CHARACTER(LEN=10), DIMENSION(50) :: measured_variables = '' !< dummy array with all measured variables that are allowed INTEGER(iwp) :: dim_ntime !< dimension size of time coordinate INTEGER(iwp) :: i !< grid index of virtual observation point in x-direction INTEGER(iwp) :: is !< grid index of real observation point of the respective station in x-direction INTEGER(iwp) :: j !< grid index of observation point in x-direction INTEGER(iwp) :: js !< grid index of real observation point of the respective station in y-direction INTEGER(iwp) :: k !< grid index of observation point in x-direction INTEGER(iwp) :: kl !< lower vertical index of surrounding grid points of an observation coordinate INTEGER(iwp) :: ks !< grid index of real observation point of the respective station in z-direction INTEGER(iwp) :: ksurf !< topography top index INTEGER(iwp) :: ku !< upper vertical index of surrounding grid points of an observation coordinate INTEGER(iwp) :: l !< running index over all stations INTEGER(iwp) :: len_char !< character length of single measured variables without Null character INTEGER(iwp) :: ll !< running index over all measured variables in file INTEGER(iwp) :: lll !< running index over all allowed variables INTEGER(iwp) :: n !< running index over trajectory coordinates INTEGER(iwp) :: ns !< counter variable for number of observation points on subdomain INTEGER(iwp) :: t !< running index over number of trajectories INTEGER(iwp) :: m INTEGER(KIND=1):: soil_dum INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ns_all !< dummy array used to sum-up the number of observation coordinates INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: meas_flag !< mask array indicating measurement positions ! LOGICAL :: chem_include !< flag indicating that chemical species is considered in modelled mechanism LOGICAL :: on_pe !< flag indicating that the respective measurement coordinate is on subdomain REAL(wp) :: fill_eutm !< _FillValue for coordinate array E_UTM REAL(wp) :: fill_nutm !< _FillValue for coordinate array N_UTM REAL(wp) :: fill_zag !< _FillValue for height coordinate REAL(wp), DIMENSION(:,:), ALLOCATABLE :: e_utm !< easting UTM coordinate, temporary variable REAL(wp), DIMENSION(:,:), ALLOCATABLE :: n_utm !< northing UTM coordinate, temporary variable REAL(wp), DIMENSION(:,:), ALLOCATABLE :: e_utm_tmp !< EUTM coordinate before rotation REAL(wp), DIMENSION(:,:), ALLOCATABLE :: n_utm_tmp !< NUTM coordinate before rotation REAL(wp), DIMENSION(:,:), ALLOCATABLE :: z_ag !< height coordinate relative to origin_z, temporary variable ! !-- Obtain number of sites. Also, pass the 'open' string, in order to initially !-- open the measurement driver. CALL netcdf_data_input_att( vmea_general%nvm, char_numstations, & vmea_general%id_vm, input_file_vm, & global_attribute, 'open', '' ) ! !-- Allocate data structure which encompass all required information, such as !-- grid points indicies, absolute UTM coordinates, the measured quantities, !-- etc. . ALLOCATE( vmea(1:vmea_general%nvm) ) ! !-- Allocate flag array. This dummy array is used to identify grid points !-- where virtual measurements should be taken. Please note, at least one !-- ghost point is required, in order to include also the surrounding !-- grid points of the original coordinate. ALLOCATE( meas_flag(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ) meas_flag = 0 ! !-- Loop over all sites. DO l = 1, vmea_general%nvm ! !-- Determine suffix which contains the ID, ordered according to the number !-- of measurements. IF( l < 10 ) THEN WRITE( dum, '(I1)') l ELSEIF( l < 100 ) THEN WRITE( dum, '(I2)') l ELSEIF( l < 1000 ) THEN WRITE( dum, '(I3)') l ELSEIF( l < 10000 ) THEN WRITE( dum, '(I4)') l ELSEIF( l < 100000 ) THEN WRITE( dum, '(I5)') l ENDIF ! !-- Read site coordinates (UTM). CALL netcdf_data_input_att( vmea(l)%origin_x_obs, char_origx // & TRIM( dum ), vmea_general%id_vm, '', & global_attribute, '', '' ) CALL netcdf_data_input_att( vmea(l)%origin_y_obs, char_origy // & TRIM( dum ), vmea_general%id_vm, '', & global_attribute, '', '' ) ! !-- Read site name CALL netcdf_data_input_att( vmea(l)%site, char_site // TRIM( dum ), & vmea_general%id_vm, '', global_attribute, & '', '' ) ! !-- Read type of the measurement (trajectory, profile, timeseries). CALL netcdf_data_input_att( vmea(l)%feature_type, char_feature // & TRIM( dum ), vmea_general%id_vm, '', & global_attribute, '', '' ) ! !-- Read the name of the original file where observational data is stored. CALL netcdf_data_input_att( vmea(l)%filename_original, char_filename // & TRIM( dum ), vmea_general%id_vm, '', & global_attribute, '', '' ) ! !-- Read a flag which indicates that also soil quantities are take at the !-- respective site (is part of the virtual measurement driver). CALL netcdf_data_input_att( soil_dum, char_soil // TRIM( dum ), & vmea_general%id_vm, '', global_attribute, & '', '' ) ! !-- Set flag for soil-sampling. IF ( soil_dum == 1 ) vmea(l)%soil_sampling = .TRUE. ! !--- Set logicals depending on the type of the measurement IF ( INDEX( vmea(l)%feature_type, type_tspr ) /= 0 ) THEN vmea(l)%timseries_profile = .TRUE. ELSEIF ( INDEX( vmea(l)%feature_type, type_ts ) /= 0 ) THEN vmea(l)%timseries = .TRUE. ELSEIF ( INDEX( vmea(l)%feature_type, type_traj ) /= 0 ) THEN vmea(l)%trajectory = .TRUE. ! !-- Give error message in case the type matches non of the pre-defined types. ELSE message_string = 'Attribue featureType = ' // & TRIM( vmea(l)%feature_type ) // & ' is not allowed.' CALL message( 'vm_init', 'PA0535', 1, 2, 0, 6, 0 ) ENDIF ! !-- Read string with all measured variables at this site measured_variables_file = '' CALL netcdf_data_input_var( measured_variables_file, & char_mv // TRIM( dum ), vmea_general%id_vm ) ! !-- Count the number of measured variables. Only count variables that match !-- with the allowed variables. !-- Please note, for some NetCDF interal reasons characters end with a NULL, !-- i.e. also empty characters contain a NULL. Therefore, check the strings !-- for a NULL to get the correct character length in order to compare !-- them with the list of allowed variables. vmea(l)%nmeas = 0 DO ll = 1, SIZE( measured_variables_file ) IF ( measured_variables_file(ll)(1:1) /= CHAR(0) .AND. & measured_variables_file(ll)(1:1) /= ' ') THEN ! !-- Obtain character length of the character len_char = 1 DO WHILE ( measured_variables_file(ll)(len_char:len_char) /= CHAR(0)& .AND. measured_variables_file(ll)(len_char:len_char) /= ' ' ) len_char = len_char + 1 ENDDO len_char = len_char - 1 ! !-- Now, compare the measured variable with the list of allowed !-- variables. DO lll= 1, SIZE( list_allowed_variables ) IF ( measured_variables_file(ll)(1:len_char) == & TRIM( list_allowed_variables(lll) ) ) THEN vmea(l)%nmeas = vmea(l)%nmeas + 1 measured_variables(vmea(l)%nmeas) = & measured_variables_file(ll)(1:len_char) ENDIF ENDDO ENDIF ENDDO ! !-- Allocate array for the measured variables names for the respective site. ALLOCATE( vmea(l)%measured_vars_name(1:vmea(l)%nmeas) ) DO ll = 1, vmea(l)%nmeas vmea(l)%measured_vars_name(ll) = TRIM( measured_variables(ll) ) ENDDO ! !-- In case of chemistry, check if species is considered in the modelled !-- chemistry mechanism. ! IF ( air_chemistry ) THEN ! DO ll = 1, vmea(l)%nmeas ! chem_include = .FALSE. ! DO n = 1, nvar ! IF ( TRIM( vmea(l)%measured_vars_name(ll) ) == & ! TRIM( chem_species(n)%name ) ) chem_include = .TRUE. ! ENDDO ! ! ! !-- Revise this. It should only check for chemistry variables and not for all! ! IF ( .NOT. chem_include ) THEN ! message_string = TRIM( vmea(l)%measured_vars_name(ll) ) // & ! ' is not considered in the modelled ' // & ! 'chemistry mechanism' ! CALL message( 'vm_init', 'PA0000', 0, 0, 0, 6, 0 ) ! ENDIF ! ENDDO ! ENDIF ! !-- Read the UTM coordinates for the actual site. Based on the coordinates, !-- define the grid-index space on each subdomain where virtual measurements !-- should be taken. Note, the entire coordinate arrays will not be stored !-- as this would exceed memory requirements, particularly for trajectory !-- measurements. IF ( vmea(l)%nmeas > 0 ) THEN ! !-- For stationary measurements UTM coordinates are just one value and !-- its dimension is "station", while for mobile measurements UTM !-- coordinates are arrays depending on the number of trajectories and !-- time, according to (UC)2 standard. First, inquire dimension length !-- of the UTM coordinates. IF ( vmea(l)%trajectory ) THEN ! !-- For non-stationary measurements read the number of trajectories !-- and the number of time coordinates. CALL netcdf_data_input_get_dimension_length( vmea_general%id_vm, & vmea(l)%ntraj, & "traj" // & TRIM( dum ) ) CALL netcdf_data_input_get_dimension_length( vmea_general%id_vm, & dim_ntime, & "ntime" // & TRIM( dum ) ) ! !-- For stationary measurements the dimension for UTM and time !-- coordinates is 1. ELSE vmea(l)%ntraj = 1 dim_ntime = 1 ENDIF ! !- Allocate array which defines individual time frame for each !-- trajectory or station. ALLOCATE( vmea(l)%dim_t(1:vmea(l)%ntraj) ) ! !-- Allocate temporary arrays for UTM and height coordinates. Note, !-- on file UTM coordinates might be 1D or 2D variables ALLOCATE( e_utm(1:vmea(l)%ntraj,1:dim_ntime) ) ALLOCATE( n_utm(1:vmea(l)%ntraj,1:dim_ntime) ) ALLOCATE( z_ag(1:vmea(l)%ntraj,1:dim_ntime) ) ALLOCATE( e_utm_tmp(1:vmea(l)%ntraj,1:dim_ntime) ) ALLOCATE( n_utm_tmp(1:vmea(l)%ntraj,1:dim_ntime) ) ! !-- Read _FillValue attributes of the coordinate dimensions. CALL netcdf_data_input_att( fill_eutm, char_fillvalue, & vmea_general%id_vm, '', & .NOT. global_attribute, '', & char_eutm // TRIM( dum ) ) CALL netcdf_data_input_att( fill_nutm, char_fillvalue, & vmea_general%id_vm, '', & .NOT. global_attribute, '', & char_nutm // TRIM( dum ) ) CALL netcdf_data_input_att( fill_zag, char_fillvalue, & vmea_general%id_vm, '', & .NOT. global_attribute, '', & char_zag // TRIM( dum ) ) ! !-- Read UTM and height coordinates coordinates for all trajectories and !-- times. IF ( vmea(l)%trajectory ) THEN CALL netcdf_data_input_var( e_utm, char_eutm // TRIM( dum ), & vmea_general%id_vm, & 0, dim_ntime-1, 0, vmea(l)%ntraj-1 ) CALL netcdf_data_input_var( n_utm, char_nutm // TRIM( dum ), & vmea_general%id_vm, & 0, dim_ntime-1, 0, vmea(l)%ntraj-1 ) CALL netcdf_data_input_var( z_ag, char_zag // TRIM( dum ), & vmea_general%id_vm, & 0, dim_ntime-1, 0, vmea(l)%ntraj-1 ) ELSE CALL netcdf_data_input_var( e_utm(1,:), char_eutm // TRIM( dum ), & vmea_general%id_vm ) CALL netcdf_data_input_var( n_utm(1,:), char_nutm // TRIM( dum ), & vmea_general%id_vm ) CALL netcdf_data_input_var( z_ag(1,:), char_zag // TRIM( dum ), & vmea_general%id_vm ) ENDIF ! !-- Based on UTM coordinates, check if the measurement station or parts !-- of the trajectory is on subdomain. This case, setup grid index space !-- sample these quantities. meas_flag = 0 DO t = 1, vmea(l)%ntraj ! !-- First, compute relative x- and y-coordinates with respect to the !-- lower-left origin of the model domain, which is the difference !-- between UTM coordinates. Note, if the origin is not correct, the !-- virtual sites will be misplaced. e_utm_tmp(t,1:dim_ntime) = e_utm(t,1:dim_ntime) - init_model%origin_x n_utm_tmp(t,1:dim_ntime) = n_utm(t,1:dim_ntime) - init_model%origin_y e_utm(t,1:dim_ntime) = COS( init_model%rotation_angle * pi / 180.0_wp ) & * e_utm_tmp(t,1:dim_ntime) & - SIN( init_model%rotation_angle * pi / 180.0_wp ) & * n_utm_tmp(t,1:dim_ntime) n_utm(t,1:dim_ntime) = SIN( init_model%rotation_angle * pi / 180.0_wp ) & * e_utm_tmp(t,1:dim_ntime) & + COS( init_model%rotation_angle * pi / 180.0_wp ) & * n_utm_tmp(t,1:dim_ntime) ! !-- Determine the individual time coordinate length for each station and !-- trajectory. This is required as several stations and trajectories !-- are merged into one file but they do not have the same number of !-- points in time, hence, missing values may occur and cannot be !-- processed further. This is actually a work-around for the specific !-- (UC)2 dataset, but it won't harm in anyway. vmea(l)%dim_t(t) = 0 DO n = 1, dim_ntime IF ( e_utm(t,n) /= fill_eutm .AND. & n_utm(t,n) /= fill_nutm .AND. & z_ag(t,n) /= fill_zag ) vmea(l)%dim_t(t) = n ENDDO ! !-- Compute grid indices relative to origin and check if these are !-- on the subdomain. Note, virtual measurements will be taken also !-- at grid points surrounding the station, hence, check also for !-- these grid points. DO n = 1, vmea(l)%dim_t(t) is = INT( ( e_utm(t,n) + 0.5_wp * dx ) * ddx, KIND = iwp ) js = INT( ( n_utm(t,n) + 0.5_wp * dy ) * ddy, KIND = iwp ) ! !-- Is the observation point on subdomain? on_pe = ( is >= nxl .AND. is <= nxr .AND. & js >= nys .AND. js <= nyn ) ! !-- Check if observation coordinate is on subdomain IF ( on_pe ) THEN ! !-- Determine vertical index which correspond to the observation !-- height. ksurf = topo_top_ind(js,is,0) ks = MINLOC( ABS( zu - zw(ksurf) - z_ag(t,n) ), DIM = 1 ) - 1 ! !-- Set mask array at the observation coordinates. Also, flag the !-- surrounding coordinate points, but first check whether the !-- surrounding coordinate points are on the subdomain. kl = MERGE( ks-1, ks, ks-1 >= nzb .AND. ks-1 >= ksurf ) ku = MERGE( ks+1, ks, ks+1 < nzt+1 ) DO i = is-1, is+1 DO j = js-1, js+1 DO k = kl, ku meas_flag(k,j,i) = MERGE( & IBSET( meas_flag(k,j,i), 0 ), & 0, & BTEST( wall_flags_0(k,j,i), 0 ) & ) ENDDO ENDDO ENDDO ENDIF ENDDO ENDDO ! !-- Based on the flag array count the number of of sampling coordinates. !-- Please note, sampling coordinates in atmosphere and soil may be !-- different, as within the soil all levels will be measured. !-- Hence, count individually. Start with atmoshere. ns = 0 DO i = nxl-1, nxr+1 DO j = nys-1, nyn+1 DO k = nzb, nzt+1 ns = ns + MERGE( 1, 0, BTEST( meas_flag(k,j,i), 0 ) ) ENDDO ENDDO ENDDO ! !-- Store number of observation points on subdomain and allocate index !-- arrays as well as array containing height information. vmea(l)%ns = ns ALLOCATE( vmea(l)%i(1:vmea(l)%ns) ) ALLOCATE( vmea(l)%j(1:vmea(l)%ns) ) ALLOCATE( vmea(l)%k(1:vmea(l)%ns) ) ALLOCATE( vmea(l)%z_ag(1:vmea(l)%ns) ) ! !-- Based on the flag array store the grid indices which correspond to !-- the observation coordinates. ns = 0 DO i = nxl-1, nxr+1 DO j = nys-1, nyn+1 DO k = nzb, nzt+1 IF ( BTEST( meas_flag(k,j,i), 0 ) ) THEN ns = ns + 1 vmea(l)%i(ns) = i vmea(l)%j(ns) = j vmea(l)%k(ns) = k vmea(l)%z_ag(ns) = zu(k) - zw(topo_top_ind(j,i,0)) ENDIF ENDDO ENDDO ENDDO ! !-- Same for the soil. Based on the flag array, count the number of !-- sampling coordinates in soil. Sample at all soil levels in this case. IF ( vmea(l)%soil_sampling ) THEN DO i = nxl, nxr DO j = nys, nyn IF ( ANY( BTEST( meas_flag(:,j,i), 0 ) ) ) THEN IF ( surf_lsm_h%start_index(j,i) <= & surf_lsm_h%end_index(j,i) ) THEN vmea(l)%ns_soil = vmea(l)%ns_soil + & nzt_soil - nzb_soil + 1 ENDIF IF ( surf_usm_h%start_index(j,i) <= & surf_usm_h%end_index(j,i) ) THEN vmea(l)%ns_soil = vmea(l)%ns_soil + & nzt_wall - nzb_wall + 1 ENDIF ENDIF ENDDO ENDDO ENDIF ! !-- Allocate index arrays as well as array containing height information !-- for soil. IF ( vmea(l)%soil_sampling ) THEN ALLOCATE( vmea(l)%i_soil(1:vmea(l)%ns_soil) ) ALLOCATE( vmea(l)%j_soil(1:vmea(l)%ns_soil) ) ALLOCATE( vmea(l)%k_soil(1:vmea(l)%ns_soil) ) ALLOCATE( vmea(l)%depth(1:vmea(l)%ns_soil) ) ENDIF ! !-- For soil, store the grid indices. ns = 0 IF ( vmea(l)%soil_sampling ) THEN DO i = nxl, nxr DO j = nys, nyn IF ( ANY( BTEST( meas_flag(:,j,i), 0 ) ) ) THEN IF ( surf_lsm_h%start_index(j,i) <= & surf_lsm_h%end_index(j,i) ) THEN m = surf_lsm_h%start_index(j,i) DO k = nzb_soil, nzt_soil ns = ns + 1 vmea(l)%i_soil(ns) = i vmea(l)%j_soil(ns) = j vmea(l)%k_soil(ns) = k vmea(l)%depth(ns) = zs(k) ENDDO ENDIF IF ( surf_usm_h%start_index(j,i) <= & surf_usm_h%end_index(j,i) ) THEN m = surf_usm_h%start_index(j,i) DO k = nzb_wall, nzt_wall ns = ns + 1 vmea(l)%i_soil(ns) = i vmea(l)%j_soil(ns) = j vmea(l)%k_soil(ns) = k vmea(l)%depth(ns) = surf_usm_h%zw(k,m) ENDDO ENDIF ENDIF ENDDO ENDDO ENDIF ! !-- Allocate array to save the sampled values. ALLOCATE( vmea(l)%measured_vars(1:vmea(l)%ns,1:vmea(l)%nmeas) ) IF ( vmea(l)%soil_sampling ) & ALLOCATE( vmea(l)%measured_vars_soil(1:vmea(l)%ns_soil, & 1:vmea(l)%nmeas) ) ! !-- Initialize with _FillValues vmea(l)%measured_vars(1:vmea(l)%ns,1:vmea(l)%nmeas) = vmea(l)%fillout IF ( vmea(l)%soil_sampling ) & vmea(l)%measured_vars_soil(1:vmea(l)%ns_soil,1:vmea(l)%nmeas) = & vmea(l)%fillout ! !-- Deallocate temporary coordinate arrays IF ( ALLOCATED( e_utm ) ) DEALLOCATE( e_utm ) IF ( ALLOCATED( n_utm ) ) DEALLOCATE( n_utm ) IF ( ALLOCATED( e_utm_tmp ) ) DEALLOCATE( e_utm_tmp ) IF ( ALLOCATED( n_utm_tmp ) ) DEALLOCATE( n_utm_tmp ) IF ( ALLOCATED( n_utm ) ) DEALLOCATE( n_utm ) IF ( ALLOCATED( z_ag ) ) DEALLOCATE( z_ag ) IF ( ALLOCATED( z_ag ) ) DEALLOCATE( vmea(l)%dim_t ) ENDIF ENDDO ! !-- Close input file for virtual measurements. Therefore, just call !-- the read attribute routine with the "close" option. CALL netcdf_data_input_att( vmea_general%nvm, char_numstations, & vmea_general%id_vm, '', & global_attribute, 'close', '' ) ! !-- Sum-up the number of observation coordiates, for atmosphere first. !-- This is actually only required for data output. ALLOCATE( ns_all(1:vmea_general%nvm) ) ns_all = 0 #if defined( __parallel ) CALL MPI_ALLREDUCE( vmea(:)%ns, ns_all(:), vmea_general%nvm, MPI_INTEGER, & MPI_SUM, comm2d, ierr ) #else ns_all(:) = vmea(:)%ns #endif vmea(:)%ns_tot = ns_all(:) ! !-- Now for soil ns_all = 0 #if defined( __parallel ) CALL MPI_ALLREDUCE( vmea(:)%ns_soil, ns_all(:), vmea_general%nvm, & MPI_INTEGER, MPI_SUM, comm2d, ierr ) #else ns_all(:) = vmea(:)%ns_soil #endif vmea(:)%ns_soil_tot = ns_all(:) DEALLOCATE( ns_all ) ! !-- Dellocate flag array DEALLOCATE( meas_flag ) ! !-- Initialize binary data output of virtual measurements. !-- Open binary output file. CALL check_open( 27 ) ! !-- Output header information. CALL vm_data_output END SUBROUTINE vm_init !------------------------------------------------------------------------------! ! Description: ! ------------ !> Binary data output. !------------------------------------------------------------------------------! SUBROUTINE vm_data_output USE pegrid IMPLICIT NONE INTEGER(iwp) :: i !< running index over IO blocks INTEGER(iwp) :: l !< running index over all stations INTEGER(iwp) :: n !< running index over all measured variables at a station ! !-- Header output on each PE IF ( init ) THEN DO i = 0, io_blocks-1 IF ( i == io_group ) THEN WRITE ( 27 ) 'number of measurements ' WRITE ( 27 ) vmea_general%nvm DO l = 1, vmea_general%nvm WRITE ( 27 ) 'site ' WRITE ( 27 ) vmea(l)%site WRITE ( 27 ) 'file ' WRITE ( 27 ) vmea(l)%filename_original WRITE ( 27 ) 'feature_type ' WRITE ( 27 ) vmea(l)%feature_type WRITE ( 27 ) 'origin_x_obs ' WRITE ( 27 ) vmea(l)%origin_x_obs WRITE ( 27 ) 'origin_y_obs ' WRITE ( 27 ) vmea(l)%origin_y_obs WRITE ( 27 ) 'total number of observation points' WRITE ( 27 ) vmea(l)%ns_tot WRITE ( 27 ) 'number of measured variables ' WRITE ( 27 ) vmea(l)%nmeas WRITE ( 27 ) 'variables ' WRITE ( 27 ) vmea(l)%measured_vars_name(:) WRITE ( 27 ) 'number of observation points ' WRITE ( 27 ) vmea(l)%ns WRITE ( 27 ) 'E_UTM ' WRITE ( 27 ) init_model%origin_x & + REAL( vmea(l)%i(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dx & * COS( init_model%rotation_angle * pi / 180.0_wp ) & + REAL( vmea(l)%j(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dy & * SIN( init_model%rotation_angle * pi / 180.0_wp ) WRITE ( 27 ) 'N_UTM ' WRITE ( 27 ) init_model%origin_y & - REAL( vmea(l)%i(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dx & * SIN( init_model%rotation_angle * pi / 180.0_wp ) & + REAL( vmea(l)%j(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dy & * COS( init_model%rotation_angle * pi / 180.0_wp ) WRITE ( 27 ) 'Z_AG ' WRITE ( 27 ) vmea(l)%z_ag(1:vmea(l)%ns) WRITE ( 27 ) 'soil sampling ' WRITE ( 27 ) MERGE( 'yes ', & 'no ', & vmea(l)%soil_sampling ) IF ( vmea(l)%soil_sampling ) THEN WRITE ( 27 ) 'total number of soil points ' WRITE ( 27 ) vmea(l)%ns_soil_tot WRITE ( 27 ) 'number of soil points ' WRITE ( 27 ) vmea(l)%ns_soil WRITE ( 27 ) 'E_UTM soil ' WRITE ( 27 ) init_model%origin_x & + REAL( vmea(l)%i_soil(1:vmea(l)%ns_soil) + 0.5_wp, & KIND = wp ) * dx & * COS( init_model%rotation_angle * pi / 180.0_wp ) & + REAL( vmea(l)%j_soil(1:vmea(l)%ns_soil) + 0.5_wp, & KIND = wp ) * dy & * SIN( init_model%rotation_angle * pi / 180.0_wp ) WRITE ( 27 ) 'N_UTM soil ' WRITE ( 27 ) init_model%origin_y & - REAL( vmea(l)%i_soil(1:vmea(l)%ns_soil) + 0.5_wp, & KIND = wp ) * dx & * SIN( init_model%rotation_angle * pi / 180.0_wp ) & + REAL( vmea(l)%j_soil(1:vmea(l)%ns_soil) + 0.5_wp, & KIND = wp ) * dy & * COS( init_model%rotation_angle * pi / 180.0_wp ) WRITE ( 27 ) 'DEPTH ' WRITE ( 27 ) vmea(l)%depth(1:vmea(l)%ns_soil) ENDIF ENDDO ENDIF ENDDO #if defined( __parallel ) CALL MPI_BARRIER( comm2d, ierr ) #endif ! !-- After header information is written, set control flag to .FALSE. init = .FALSE. ! !-- Data output at each measurement timestep on each PE ELSE DO i = 0, io_blocks-1 IF ( i == io_group ) THEN WRITE( 27 ) 'output time ' WRITE( 27 ) time_since_reference_point DO l = 1, vmea_general%nvm ! !-- Skip binary writing if no observation points are defined on PE IF ( vmea(l)%ns < 1 .AND. vmea(l)%ns_soil < 1) CYCLE DO n = 1, vmea(l)%nmeas WRITE( 27 ) vmea(l)%measured_vars_name(n) IF ( vmea(l)%soil_sampling .AND. & ANY( TRIM( vmea(l)%measured_vars_name(n)) == & soil_vars ) ) THEN WRITE( 27 ) vmea(l)%measured_vars_soil(:,n) ELSE WRITE( 27 ) vmea(l)%measured_vars(:,n) ENDIF ENDDO ENDDO ENDIF ENDDO #if defined( __parallel ) CALL MPI_BARRIER( comm2d, ierr ) #endif ENDIF END SUBROUTINE vm_data_output !------------------------------------------------------------------------------! ! Description: ! ------------ !> Write end-of-file statement as last action. !------------------------------------------------------------------------------! SUBROUTINE vm_last_actions USE pegrid IMPLICIT NONE INTEGER(iwp) :: i !< running index over IO blocks DO i = 0, io_blocks-1 IF ( i == io_group ) THEN WRITE( 27 ) 'EOF ' ENDIF ENDDO #if defined( __parallel ) CALL MPI_BARRIER( comm2d, ierr ) #endif ! !-- Close binary file CALL close_file( 27 ) END SUBROUTINE vm_last_actions !------------------------------------------------------------------------------! ! Description: ! ------------ !> Sampling of the actual quantities along the observation coordinates !------------------------------------------------------------------------------! SUBROUTINE vm_sampling USE arrays_3d, & ONLY: exner, pt, q, u, v, w USE radiation_model_mod, & ONLY: radiation USE surface_mod, & ONLY: surf_def_h, surf_lsm_h, surf_usm_h IMPLICIT NONE INTEGER(iwp) :: i !< grid index in x-direction INTEGER(iwp) :: j !< grid index in y-direction INTEGER(iwp) :: k !< grid index in z-direction INTEGER(iwp) :: ind_chem !< dummy index to identify chemistry variable and translate it from (UC)2 standard to interal naming INTEGER(iwp) :: l !< running index over the number of stations INTEGER(iwp) :: m !< running index over all virtual observation coordinates INTEGER(iwp) :: mm !< index of surface element which corresponds to the virtual observation coordinate INTEGER(iwp) :: n !< running index over all measured variables at a station INTEGER(iwp) :: nn !< running index over the number of chemcal species LOGICAL :: match_lsm !< flag indicating natural-type surface LOGICAL :: match_usm !< flag indicating urban-type surface ! !-- Loop over all sites. DO l = 1, vmea_general%nvm ! !-- At the beginning, set _FillValues IF ( ALLOCATED( vmea(l)%measured_vars ) ) & vmea(l)%measured_vars = vmea(l)%fillout IF ( ALLOCATED( vmea(l)%measured_vars_soil ) ) & vmea(l)%measured_vars_soil = vmea(l)%fillout ! !-- Loop over all variables measured at this site. DO n = 1, vmea(l)%nmeas SELECT CASE ( TRIM( vmea(l)%measured_vars_name(n) ) ) CASE ( 'theta' ) IF ( .NOT. neutral ) THEN DO m = 1, vmea(l)%ns k = vmea(l)%k(m) j = vmea(l)%j(m) i = vmea(l)%i(m) vmea(l)%measured_vars(m,n) = pt(k,j,i) ENDDO ENDIF CASE ( 'ta' ) IF ( .NOT. neutral ) THEN DO m = 1, vmea(l)%ns k = vmea(l)%k(m) j = vmea(l)%j(m) i = vmea(l)%i(m) vmea(l)%measured_vars(m,n) = pt(k,j,i) * exner( k ) ENDDO ENDIF CASE ( 't_va' ) CASE ( 'hus', 'haa' ) IF ( humidity ) THEN DO m = 1, vmea(l)%ns k = vmea(l)%k(m) j = vmea(l)%j(m) i = vmea(l)%i(m) vmea(l)%measured_vars(m,n) = q(k,j,i) ENDDO ENDIF CASE ( 'u', 'ua' ) DO m = 1, vmea(l)%ns k = vmea(l)%k(m) j = vmea(l)%j(m) i = vmea(l)%i(m) vmea(l)%measured_vars(m,n) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) ENDDO CASE ( 'v', 'va' ) DO m = 1, vmea(l)%ns k = vmea(l)%k(m) j = vmea(l)%j(m) i = vmea(l)%i(m) vmea(l)%measured_vars(m,n) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) ENDDO CASE ( 'w' ) DO m = 1, vmea(l)%ns k = MAX ( 1, vmea(l)%k(m) ) j = vmea(l)%j(m) i = vmea(l)%i(m) vmea(l)%measured_vars(m,n) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) ) ENDDO CASE ( 'wspeed' ) DO m = 1, vmea(l)%ns k = vmea(l)%k(m) j = vmea(l)%j(m) i = vmea(l)%i(m) vmea(l)%measured_vars(m,n) = SQRT( & ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2 + & ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2 & ) ENDDO CASE ( 'wdir' ) DO m = 1, vmea(l)%ns k = vmea(l)%k(m) j = vmea(l)%j(m) i = vmea(l)%i(m) vmea(l)%measured_vars(m,n) = ATAN2( & - 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ), & - 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) & ) * 180.0_wp / pi ENDDO CASE ( 'utheta' ) DO m = 1, vmea(l)%ns k = vmea(l)%k(m) j = vmea(l)%j(m) i = vmea(l)%i(m) vmea(l)%measured_vars(m,n) = 0.5_wp * & ( u(k,j,i) + u(k,j,i+1) ) * & pt(k,j,i) ENDDO CASE ( 'vtheta' ) DO m = 1, vmea(l)%ns k = vmea(l)%k(m) j = vmea(l)%j(m) i = vmea(l)%i(m) vmea(l)%measured_vars(m,n) = 0.5_wp * & ( v(k,j,i) + v(k,j+1,i) ) * & pt(k,j,i) ENDDO CASE ( 'wtheta' ) DO m = 1, vmea(l)%ns k = MAX ( 1, vmea(l)%k(m) ) j = vmea(l)%j(m) i = vmea(l)%i(m) vmea(l)%measured_vars(m,n) = 0.5_wp * & ( w(k-1,j,i) + w(k,j,i) ) * & pt(k,j,i) ENDDO CASE ( 'uw' ) DO m = 1, vmea(l)%ns k = MAX ( 1, vmea(l)%k(m) ) j = vmea(l)%j(m) i = vmea(l)%i(m) vmea(l)%measured_vars(m,n) = 0.25_wp * & ( w(k-1,j,i) + w(k,j,i) ) * & ( u(k,j,i) + u(k,j,i+1) ) ENDDO CASE ( 'vw' ) DO m = 1, vmea(l)%ns k = MAX ( 1, vmea(l)%k(m) ) j = vmea(l)%j(m) i = vmea(l)%i(m) vmea(l)%measured_vars(m,n) = 0.25_wp * & ( w(k-1,j,i) + w(k,j,i) ) * & ( v(k,j,i) + v(k,j+1,i) ) ENDDO CASE ( 'uv' ) DO m = 1, vmea(l)%ns k = MAX ( 1, vmea(l)%k(m) ) j = vmea(l)%j(m) i = vmea(l)%i(m) vmea(l)%measured_vars(m,n) = 0.25_wp * & ( u(k,j,i) + u(k,j,i+1) ) * & ( v(k,j,i) + v(k,j+1,i) ) ENDDO ! !-- List of variables may need extension. CASE ( 'mcpm1', 'mcpm2p5', 'mcpm10', 'mfco', 'mfno', 'mfno2', & 'tro3' ) IF ( air_chemistry ) THEN ! !-- First, search for the measured variable in the chem_vars !-- list, in order to get the internal name of the variable. DO nn = 1, UBOUND( chem_vars, 2 ) IF ( TRIM( vmea(l)%measured_vars_name(m) ) == & TRIM( chem_vars(0,nn) ) ) ind_chem = nn ENDDO ! !-- Run loop over all chemical species, if the measured !-- variable matches the interal name, sample the variable. DO nn = 1, nvar IF ( TRIM( chem_vars(1,ind_chem) ) == & TRIM( chem_species(nn)%name ) ) THEN DO m = 1, vmea(l)%ns k = vmea(l)%k(m) j = vmea(l)%j(m) i = vmea(l)%i(m) vmea(l)%measured_vars(m,n) = & chem_species(nn)%conc(k,j,i) ENDDO ENDIF ENDDO ENDIF CASE ( 'us' ) DO m = 1, vmea(l)%ns ! !-- Surface data is only available on inner subdomains, not !-- on ghost points. Hence, limit the indices. j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys ) j = MERGE( j , nyn, j > nyn ) i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl ) i = MERGE( i , nxr, i > nxr ) DO mm = surf_def_h(0)%start_index(j,i), & surf_def_h(0)%end_index(j,i) vmea(l)%measured_vars(m,n) = surf_def_h(0)%us(mm) ENDDO DO mm = surf_lsm_h%start_index(j,i), & surf_lsm_h%end_index(j,i) vmea(l)%measured_vars(m,n) = surf_lsm_h%us(mm) ENDDO DO mm = surf_usm_h%start_index(j,i), & surf_usm_h%end_index(j,i) vmea(l)%measured_vars(m,n) = surf_usm_h%us(mm) ENDDO ENDDO CASE ( 'ts' ) DO m = 1, vmea(l)%ns ! !-- Surface data is only available on inner subdomains, not !-- on ghost points. Hence, limit the indices. j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys ) j = MERGE( j , nyn, j > nyn ) i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl ) i = MERGE( i , nxr, i > nxr ) DO mm = surf_def_h(0)%start_index(j,i), & surf_def_h(0)%end_index(j,i) vmea(l)%measured_vars(m,n) = surf_def_h(0)%ts(mm) ENDDO DO mm = surf_lsm_h%start_index(j,i), & surf_lsm_h%end_index(j,i) vmea(l)%measured_vars(m,n) = surf_lsm_h%ts(mm) ENDDO DO mm = surf_usm_h%start_index(j,i), & surf_usm_h%end_index(j,i) vmea(l)%measured_vars(m,n) = surf_usm_h%ts(mm) ENDDO ENDDO CASE ( 'hfls' ) DO m = 1, vmea(l)%ns ! !-- Surface data is only available on inner subdomains, not !-- on ghost points. Hence, limit the indices. j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys ) j = MERGE( j , nyn, j > nyn ) i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl ) i = MERGE( i , nxr, i > nxr ) DO mm = surf_def_h(0)%start_index(j,i), & surf_def_h(0)%end_index(j,i) vmea(l)%measured_vars(m,n) = surf_def_h(0)%qsws(mm) ENDDO DO mm = surf_lsm_h%start_index(j,i), & surf_lsm_h%end_index(j,i) vmea(l)%measured_vars(m,n) = surf_lsm_h%qsws(mm) ENDDO DO mm = surf_usm_h%start_index(j,i), & surf_usm_h%end_index(j,i) vmea(l)%measured_vars(m,n) = surf_usm_h%qsws(mm) ENDDO ENDDO CASE ( 'hfss' ) DO m = 1, vmea(l)%ns ! !-- Surface data is only available on inner subdomains, not !-- on ghost points. Hence, limit the indices. j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys ) j = MERGE( j , nyn, j > nyn ) i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl ) i = MERGE( i , nxr, i > nxr ) DO mm = surf_def_h(0)%start_index(j,i), & surf_def_h(0)%end_index(j,i) vmea(l)%measured_vars(m,n) = surf_def_h(0)%shf(mm) ENDDO DO mm = surf_lsm_h%start_index(j,i), & surf_lsm_h%end_index(j,i) vmea(l)%measured_vars(m,n) = surf_lsm_h%shf(mm) ENDDO DO mm = surf_usm_h%start_index(j,i), & surf_usm_h%end_index(j,i) vmea(l)%measured_vars(m,n) = surf_usm_h%shf(mm) ENDDO ENDDO CASE ( 'rnds' ) IF ( radiation ) THEN DO m = 1, vmea(l)%ns ! !-- Surface data is only available on inner subdomains, not !-- on ghost points. Hence, limit the indices. j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys ) j = MERGE( j , nyn, j > nyn ) i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl ) i = MERGE( i , nxr, i > nxr ) DO mm = surf_lsm_h%start_index(j,i), & surf_lsm_h%end_index(j,i) vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_net(mm) ENDDO DO mm = surf_usm_h%start_index(j,i), & surf_usm_h%end_index(j,i) vmea(l)%measured_vars(m,n) = surf_usm_h%rad_net(mm) ENDDO ENDDO ENDIF CASE ( 'rsus' ) IF ( radiation ) THEN DO m = 1, vmea(l)%ns ! !-- Surface data is only available on inner subdomains, not !-- on ghost points. Hence, limit the indices. j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys ) j = MERGE( j , nyn, j > nyn ) i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl ) i = MERGE( i , nxr, i > nxr ) DO mm = surf_lsm_h%start_index(j,i), & surf_lsm_h%end_index(j,i) vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_sw_out(mm) ENDDO DO mm = surf_usm_h%start_index(j,i), & surf_usm_h%end_index(j,i) vmea(l)%measured_vars(m,n) = surf_usm_h%rad_sw_out(mm) ENDDO ENDDO ENDIF CASE ( 'rsds' ) IF ( radiation ) THEN DO m = 1, vmea(l)%ns ! !-- Surface data is only available on inner subdomains, not !-- on ghost points. Hence, limit the indices. j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys ) j = MERGE( j , nyn, j > nyn ) i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl ) i = MERGE( i , nxr, i > nxr ) DO mm = surf_lsm_h%start_index(j,i), & surf_lsm_h%end_index(j,i) vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_sw_in(mm) ENDDO DO mm = surf_usm_h%start_index(j,i), & surf_usm_h%end_index(j,i) vmea(l)%measured_vars(m,n) = surf_usm_h%rad_sw_in(mm) ENDDO ENDDO ENDIF CASE ( 'rlus' ) IF ( radiation ) THEN DO m = 1, vmea(l)%ns ! !-- Surface data is only available on inner subdomains, not !-- on ghost points. Hence, limit the indices. j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys ) j = MERGE( j , nyn, j > nyn ) i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl ) i = MERGE( i , nxr, i > nxr ) DO mm = surf_lsm_h%start_index(j,i), & surf_lsm_h%end_index(j,i) vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_lw_out(mm) ENDDO DO mm = surf_usm_h%start_index(j,i), & surf_usm_h%end_index(j,i) vmea(l)%measured_vars(m,n) = surf_usm_h%rad_lw_out(mm) ENDDO ENDDO ENDIF CASE ( 'rlds' ) IF ( radiation ) THEN DO m = 1, vmea(l)%ns ! !-- Surface data is only available on inner subdomains, not !-- on ghost points. Hence, limit the indices. j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys ) j = MERGE( j , nyn, j > nyn ) i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl ) i = MERGE( i , nxr, i > nxr ) DO mm = surf_lsm_h%start_index(j,i), & surf_lsm_h%end_index(j,i) vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_lw_in(mm) ENDDO DO mm = surf_usm_h%start_index(j,i), & surf_usm_h%end_index(j,i) vmea(l)%measured_vars(m,n) = surf_usm_h%rad_lw_in(mm) ENDDO ENDDO ENDIF CASE ( 'rsd' ) IF ( radiation ) THEN DO m = 1, vmea(l)%ns k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' ) j = vmea(l)%j(m) i = vmea(l)%i(m) vmea(l)%measured_vars(m,n) = rad_sw_in(k,j,i) ENDDO ENDIF CASE ( 'rsu' ) IF ( radiation ) THEN DO m = 1, vmea(l)%ns k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' ) j = vmea(l)%j(m) i = vmea(l)%i(m) vmea(l)%measured_vars(m,n) = rad_sw_out(k,j,i) ENDDO ENDIF CASE ( 'rlu' ) IF ( radiation ) THEN DO m = 1, vmea(l)%ns k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' ) j = vmea(l)%j(m) i = vmea(l)%i(m) vmea(l)%measured_vars(m,n) = rad_lw_out(k,j,i) ENDDO ENDIF CASE ( 'rld' ) IF ( radiation ) THEN DO m = 1, vmea(l)%ns k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' ) j = vmea(l)%j(m) i = vmea(l)%i(m) vmea(l)%measured_vars(m,n) = rad_lw_in(k,j,i) ENDDO ENDIF CASE ( 'rsddif' ) IF ( radiation ) THEN DO m = 1, vmea(l)%ns j = vmea(l)%j(m) i = vmea(l)%i(m) vmea(l)%measured_vars(m,n) = rad_sw_in_diff(j,i) ENDDO ENDIF CASE ( 't_soil' ) DO m = 1, vmea(l)%ns_soil i = vmea(l)%i_soil(m) j = vmea(l)%j_soil(m) k = vmea(l)%k_soil(m) match_lsm = surf_lsm_h%start_index(j,i) <= & surf_lsm_h%end_index(j,i) match_usm = surf_usm_h%start_index(j,i) <= & surf_usm_h%end_index(j,i) IF ( match_lsm ) THEN mm = surf_lsm_h%start_index(j,i) vmea(l)%measured_vars_soil(m,n) = t_soil_h%var_2d(k,mm) ENDIF IF ( match_usm ) THEN mm = surf_usm_h%start_index(j,i) vmea(l)%measured_vars_soil(m,n) = t_wall_h(k,mm) ENDIF ENDDO CASE ( 'm_soil' ) DO m = 1, vmea(l)%ns_soil i = vmea(l)%i_soil(m) j = vmea(l)%j_soil(m) k = vmea(l)%k_soil(m) match_lsm = surf_lsm_h%start_index(j,i) <= & surf_lsm_h%end_index(j,i) IF ( match_lsm ) THEN mm = surf_lsm_h%start_index(j,i) vmea(l)%measured_vars_soil(m,n) = m_soil_h%var_2d(k,mm) ENDIF ENDDO ! !-- More will follow ... ! !-- No match found - just set a fill value CASE DEFAULT vmea(l)%measured_vars(:,n) = vmea(l)%fillout END SELECT ENDDO ENDDO END SUBROUTINE vm_sampling END MODULE virtual_measurement_mod