!> @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: ! ----------------- ! Initial revision ! ! Former revisions: ! ----------------- ! $Id: virtual_measurement_mod.f90 3472 2018-10-30 20:43:50Z suehring $ ! ! 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' measurement positions. !> 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. !------------------------------------------------------------------------------! MODULE virtual_measurement_mod USE arrays_3d, & ONLY: q, pt, u, v, w, zu, zw USE control_parameters, & ONLY: dz, message_string, virtual_measurement USE cpulog, & ONLY: cpu_log, log_point USE grid_variables, & ONLY: dx, dy USE indices, & ONLY: nzb, nzt, nxl, nxr, nys, nyn USE kinds IMPLICIT NONE TYPE virt_mea CHARACTER(LEN=100) :: feature_type !< type of the measurement 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 !< total number of observation points for a site on subdomain, i.e. sum of all trajectories INTEGER(iwp) :: ntraj !< number of trajectories of a measurement INTEGER(iwp) :: nvar !< number of measured variables INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: dim_t !< number observations individual for each trajectory or station that are no _FillValues INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ngp !< number of grid points where observations for a site took place, ! 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 IMPLICIT NONE ! !-- 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', 'PA0000', 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/ 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 control_parameters, & ONLY: message_string 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 USE surface_mod, & ONLY: get_topography_top_index_ji 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 LOGICAL :: on_pe !< flag indicating that the respective measurement coordinate is on subdomain INTEGER(iwp) :: dim_eutm !< dimension size of UTM easting coordinate INTEGER(iwp) :: dim_nutm !< dimension size of UTM northing coordinate INTEGER(iwp) :: dim_ntime !< dimension size of time coordinate INTEGER(iwp) :: dim_zag !< dimension size of height coordinate INTEGER(iwp) :: i !< grid index of virtual observation point in x-direction INTEGER(iwp) :: icov !< index range where observations should be taken in x-direction INTEGER(iwp) :: ii !< running index over all coordinate points of a measurement INTEGER(iwp) :: i_prev !< grid index along x for UTM coordinate at previous observation time step 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) :: jcov !< index range where observations should be taken in y-direction INTEGER(iwp) :: j_prev !< grid index along y for UTM coordinate at previous observation time step 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) :: kcov !< index range where observations should be taken in z-direction INTEGER(iwp) :: ks !< grid index of real observation point of the respective station in z-direction INTEGER(iwp) :: k_prev !< grid index along z for UTM coordinate at previous observation time step INTEGER(iwp) :: ksurf !< topography top index 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 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 :: z_ag !< height coordinate relative to origin_z, temporary variable ! !-- Obtain number of virtual measurement stations CALL netcdf_data_input_att( nvm, char_numstations, id_vm, input_file_vm, & global_attribute, 'open', '' ) ! !-- ALLOCATE data structure ALLOCATE( vmea(1:nvm) ) ! print*, "nvm", nvm ! !-- Read station coordinates and further attributes. !-- Note all coordinates are in UTM coordinates. DO l = 1, nvm ! !-- Determine suffix which contains the ID 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 CALL netcdf_data_input_att( vmea(l)%origin_x_obs, char_origx & // TRIM( dum ), id_vm, '', global_attribute,& '', '' ) CALL netcdf_data_input_att( vmea(l)%origin_y_obs, char_origy & // TRIM( dum ), id_vm, '', global_attribute,& '', '' ) CALL netcdf_data_input_att( vmea(l)%site, char_site & // TRIM( dum ), id_vm, '', global_attribute,& '', '' ) CALL netcdf_data_input_att( vmea(l)%feature_type, char_feature & // TRIM( dum ), id_vm, '', global_attribute,& '', '' ) ! !--- 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. ELSE ! !-- Give error message message_string = 'Attribue featureType = ' // & TRIM( vmea(l)%feature_type ) // & ' is not allowed.' CALL message( 'vm_init', 'PA0000', 1, 2, 0, 6, 0 ) ENDIF ! !-- Read string with all measured variables at this station measured_variables_file = '' CALL netcdf_data_input_var( measured_variables_file, & char_mv // TRIM( dum ), id_vm ) ! !-- Count the number of measured variables which match with the variables !-- which are allowed to be measured in PALM. 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)%nvar = 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)%nvar = vmea(l)%nvar + 1 measured_variables(vmea(l)%nvar) = & measured_variables_file(ll)(1:len_char) ENDIF ENDDO ENDIF ENDDO ! !-- Allocate array for the measured variables names for the station l. ALLOCATE( vmea(l)%measured_vars_name(1:vmea(l)%nvar) ) DO ll = 1, vmea(l)%nvar vmea(l)%measured_vars_name(ll) = TRIM( measured_variables(ll) ) ENDDO ! print*, "numvars", vmea(l)%nvar, vmea(l)%measured_vars_name(1:vmea(l)%nvar) ! !-- For the actual measurement ID read the UTM coordinates. Based on these, !-- define the index space on each subdomain where measurements should be !-- taken. Note, the entire coordinate arrays will not be stored on data !-- type as this would exceed memory requirements, particularly for !-- trajectory measurements. If no variable will be virtually measured, !-- skip the reading. IF ( vmea(l)%nvar > 0 ) THEN ! !-- For stationary measurements UTM coordinates are just one value and !-- its dimension is "station", while for mobile measurements UTM !-- coordinates are arrays. First, inquire dimension length for !-- UTM coordinates. IF ( vmea(l)%trajectory ) THEN ! !-- For non-stationary measurements read the number of trajectories CALL netcdf_data_input_get_dimension_length( id_vm, & vmea(l)%ntraj, & "traj" // & TRIM( dum ) ) CALL netcdf_data_input_get_dimension_length( id_vm, dim_ntime, & "ntime" // & TRIM( dum ) ) ! !-- For stationary measurements the dimension for UTM 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( vmea(l)%ngp(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) ) ! !-- Read _FillValue attributes CALL netcdf_data_input_att( fill_eutm, char_fillvalue, & id_vm, '', .NOT. global_attribute, '', & char_eutm // TRIM( dum ) ) CALL netcdf_data_input_att( fill_nutm, char_fillvalue, & id_vm, '', .NOT. global_attribute, '', & char_nutm // TRIM( dum ) ) CALL netcdf_data_input_att( fill_zag, char_fillvalue, & 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 ), id_vm, & 0, dim_ntime-1, 0, vmea(l)%ntraj-1 ) CALL netcdf_data_input_var( n_utm, char_nutm // TRIM( dum ), id_vm, & 0, dim_ntime-1, 0, vmea(l)%ntraj-1 ) CALL netcdf_data_input_var( z_ag, char_zag // TRIM( dum ), id_vm, & 0, dim_ntime-1, 0, vmea(l)%ntraj-1 ) ELSE CALL netcdf_data_input_var( e_utm(1,:), char_eutm // TRIM( dum ), id_vm ) CALL netcdf_data_input_var( n_utm(1,:), char_nutm // TRIM( dum ), id_vm ) CALL netcdf_data_input_var( z_ag(1,:), char_zag // TRIM( dum ), 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. ns = 0 DO t = 1, vmea(l)%ntraj ! !-- 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. 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 ! !-- First, compute relative x- and y-coordinates with respect to the !-- lower-left origin of the model domain, which is the difference !-- betwen UTM coordinates. ! write(9,*) l, "before Eutm", e_utm(t,1:vmea(l)%dim_t(t)), "orig", init_model%origin_x ! write(9,*) l, "before Nutm", n_utm(t,1:vmea(l)%dim_t(t)), "orig", init_model%origin_y e_utm(t,1:vmea(l)%dim_t(t)) = e_utm(t,1:vmea(l)%dim_t(t)) & - init_model%origin_x n_utm(t,1:vmea(l)%dim_t(t)) = n_utm(t,1:vmea(l)%dim_t(t)) & - init_model%origin_y ! write(9,*) l, "Eutm", e_utm(t,1:vmea(l)%dim_t(t)) ! write(9,*) l, "Nutm", n_utm(t,1:vmea(l)%dim_t(t)) ! write(9,*) ! !-- 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. vmea(l)%ngp(t) = 0 k_prev = -999 j_prev = -999 i_prev = -999 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 ) ! !-- If the measurement is on subdomain, determine vertical index !-- which refers to the observation height above ground level. ks = k_prev IF ( on_pe ) THEN ! write(9,*) " I am onpe" ksurf = get_topography_top_index_ji( js, is, 's' ) ks = MINLOC( ABS( zu - zw(ksurf) - z_ag(t,n) ), DIM = 1 ) - 1 ENDIF ! !-- Count the number of observation points in index space on !-- subdomain. Only increment if grid indices are different from !-- the previous one. IF ( on_pe .AND. is /= i_prev .AND. js /= j_prev .AND. & ks /= k_prev ) THEN ns = ns + 1 vmea(l)%ngp(t) = vmea(l)%ngp(t) + 1 ENDIF !-- Store arrays for next iteration - avoid double counting i_prev = is j_prev = js k_prev = ks ENDDO ENDDO ! !-- Store number of observation points on subdomain and allocate index !-- arrays. 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) ) ! print*, "Num ns: ", vmea(l)%ns, "per traj", vmea(l)%ngp(:) ! !-- Repeat the prior loop and save the grid indices relevant for !-- sampling. ns = 0 DO t = 1, vmea(l)%ntraj ! !-- 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. k_prev = -999 j_prev = -999 i_prev = -999 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 ) ! !-- If the measurement is on subdomain, determine vertical index !-- which refers to the observation height above ground level. ks = k_prev IF ( on_pe ) THEN ksurf = get_topography_top_index_ji( js, is, 's' ) ks = MINLOC( ABS( zu - zw(ksurf) - z_ag(t,n) ), DIM = 1 ) - 1 ENDIF ! !-- Count the number of observation points in index space on !-- subdomain. Only increment if grid indices are different from !-- the previous one. IF ( on_pe .AND. is /= i_prev .AND. js /= j_prev .AND. & ks /= k_prev ) THEN ns = ns + 1 vmea(l)%i(ns) = is vmea(l)%j(ns) = js vmea(l)%k(ns) = ks ! write(9,*) TRIM( vmea(l)%feature_type), l, "ns", ns, "ijk", vmea(l)%i(ns), vmea(l)%j(ns), vmea(l)%k(ns) ENDIF ! !-- Store arrays for next iteration - avoid double counting i_prev = is j_prev = js k_prev = ks ENDDO ENDDO ! !-- Allocate array to save the sampled values. !-- Todo: Is it better to allocate for all variables at a station !-- and store all the values before writing, or sample the variables !-- directly in the data output? ALLOCATE( vmea(l)%measured_vars(1:vmea(l)%nvar,1:vmea(l)%ns) ) ! !-- Initialize with _FillValue vmea(l)%measured_vars(1:vmea(l)%nvar,1:vmea(l)%ns) = vmea(l)%fillout ! !-- Deallocate temporary coordinate arrays IF ( ALLOCATED( e_utm ) ) DEALLOCATE( e_utm ) IF ( ALLOCATED( n_utm ) ) DEALLOCATE( n_utm ) IF ( ALLOCATED( z_ag ) ) DEALLOCATE( z_ag ) ENDIF ENDDO flush(9) ! !-- Close input file for virtual measurements. Therefore, just call !-- the read attribute routine with the "close" option. CALL netcdf_data_input_att( nvm, char_numstations, id_vm, '', & global_attribute, 'close', '' ) END SUBROUTINE vm_init !------------------------------------------------------------------------------! ! Description: ! ------------ !> Sampling of the actual quantities along the observation coordinates !------------------------------------------------------------------------------! SUBROUTINE vm_sampling USE arrays_3d !, & ! ONLY: pt USE surface_mod IMPLICIT NONE CHARACTER(LEN=10) :: trimvar !< dummy for the measured variable name INTEGER(iwp) :: l !< INTEGER(iwp) :: m !< INTEGER(iwp) :: var !< INTEGER(iwp) :: mm, j, i ! write(9,*) "sampling" ! flush(9) ! !-- Loop over all stations. For each possible variable loop over all !-- observation points DO l = 1, nvm ! !-- Loop over all measured variables. Please note, for the moment !-- the same indices for scalar and velocity components are used. !-- ToDo: Revise this later. ! DO m = 1, vmea(l)%ns ! j = vmea(l)%j(m) ! i = vmea(l)%i(m) ! ! IF ( i >= nxl .AND. i <= nxr .AND. & ! j >= nys .AND. j <= nyn ) THEN ! IF ( surf_def_h(0)%start_index(j,i) <= & ! surf_def_h(0)%end_index(j,i) ) THEN ! ! write(9,*) "sampled" ! flush(9) ! mm = surf_def_h(0)%end_index(j,i) ! ! surf_def_h(0)%pt_surface(mm) = -99.0 ! ENDIF ! ENDIF ! ENDDO ! DO var = 1, vmea(l)%nvar ! trimvar = TRIM( vmea(l)%measured_vars_name(var) ) ! ! IF ( TRIM( trimvar ) == 'theta' ) THEN ! DO m = 1, vmea(l)%ns ! vmea(l)%measured_vars(var,m) = pt(vmea(l)%k(m),vmea(l)%j(m),vmea(l)%i(m)) ! ENDDO ! ENDIF ! IF ( TRIM( trimvar ) == 'w' ) THEN ! DO m = 1, vmea(l)%ns ! vmea(l)%measured_vars(var,m) = w(vmea(l)%k(m),vmea(l)%j(m),vmea(l)%i(m)) ! ENDDO ! ENDIF ! IF ( TRIM( trimvar ) == 'v' ) THEN ! DO m = 1, vmea(l)%ns ! vmea(l)%measured_vars(var,m) = v(vmea(l)%k(m),vmea(l)%j(m),vmea(l)%i(m)) ! ENDDO ! ENDIF ! IF ( TRIM( trimvar ) == 'u' ) THEN ! DO m = 1, vmea(l)%ns ! vmea(l)%measured_vars(var,m) = u(vmea(l)%k(m),vmea(l)%j(m),vmea(l)%i(m)) ! ENDDO ! ENDIF ! ENDDO ENDDO END SUBROUTINE vm_sampling END MODULE virtual_measurement_mod