diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_io/unc_write_his.F90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_io/unc_write_his.F90 index cf553deb20..0cd5a7a0a1 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_io/unc_write_his.F90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_io/unc_write_his.F90 @@ -256,9 +256,10 @@ subroutine unc_write_his(tim) ! wrihis call check_netcdf_error(nf90_def_dim(ihisfile, 'name_len', strlen_netcdf, id_strlendim)) - if (kmx > 0) then - call check_netcdf_error(nf90_def_dim(ihisfile, 'laydim', kmx, id_laydim)) - call check_netcdf_error(nf90_def_dim(ihisfile, 'laydimw', kmx + 1, id_laydimw)) + ! TK_Temp: Allow for writing in case of kmx = 0 + if (kmx >= 0) then + call check_netcdf_error(nf90_def_dim(ihisfile, 'laydim', max(kmx,1) , id_laydim)) + call check_netcdf_error(nf90_def_dim(ihisfile, 'laydimw',max(kmx,1) + 1, id_laydimw)) end if if (stm_included .and. jahissed > 0) then @@ -1070,8 +1071,9 @@ function unc_def_his_station_coord_vars_z(ihisfile, id_laydim, id_laydimw, id_st type(ug_nc_attribute) :: extra_attributes(1) ierr = DFM_NOERR + ! TK_Temp: Also write vertical coordinates to his file for kmx == 0! if (.not. model_is_3D()) then - return +! return end if call ncu_set_att(extra_attributes(1), 'positive', 'up') ! If so specified, add the zcoordinate_c @@ -1228,19 +1230,21 @@ function unc_put_his_station_coord_vars_z(ihisfile, numobs, nummovobs, jawrizc, integer :: layer ierr = DFM_NOERR - + ! TK_Temp: Write vertical coordinates for kmx = 0 if (.not. model_is_3D()) then - return +! return end if if (jawrizc == 1) then - do layer = 1, kmx + ! TK_Temp: for KMX = 0 + do layer = 1, max(kmx,1) call check_netcdf_error(nf90_put_var(ihisfile, id_zcs, valobs(:, IPNT_ZCS + layer - 1), start=[layer, 1, it_his], count=[1, numobs + nummovobs, 1])) end do end if if (jawrizw == 1) then - do layer = 1, kmx + 1 + ! TK_Temp: for KMX = 0 + do layer = 1, max(kmx,1) + 1 call check_netcdf_error(nf90_put_var(ihisfile, id_zws, valobs(:, IPNT_ZWS + layer - 1), start=[layer, 1, it_his], count=[1, numobs + nummovobs, 1])) call check_netcdf_error(nf90_put_var(ihisfile, id_zwu, valobs(:, IPNT_ZWU + layer - 1), start=[layer, 1, it_his], count=[1, numobs + nummovobs, 1])) end do diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/fill_valobs.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/fill_valobs.f90 index badc2de48d..87c69b195c 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/fill_valobs.f90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/fill_valobs.f90 @@ -56,22 +56,22 @@ subroutine fill_valobs() use m_flowtimes, only: handle_extra use m_transport, only: constituents, isalt, itemp, itra1, ised1 use m_flowgeom, only: ndx, lnx, bl, nd, ln, wcl, bob, ba - use m_observations_data, only: valobs, numobs, nummovobs, kobs, lobs, ipnt_s1, ipnt_hs, ipnt_bl, ipnt_cmx, cmxobs, & - ipnt_wx, ipnt_wy, ipnt_patm, ipnt_waver, ipnt_waveh, ipnt_wavet, ipnt_waved, ipnt_wavel, ipnt_waveu, ipnt_taux, & - ipnt_tauy, ival_sbcx1, ival_sbcxn, ipnt_sbcx1, ival_sbcy1, ival_sbcyn, ipnt_sbcy1, ival_sscx1, ival_sscxn, & - ipnt_sscx1, ival_sscy1, ival_sscyn, ipnt_sscy1, ival_sbwx1, ival_sbwxn, ipnt_sbwx1, ival_sbwy1, ival_sbwyn, & - ipnt_sbwy1, ival_sswx1, ival_sswxn, ipnt_sswx1, ival_sswy1, ival_sswyn, ipnt_sswy1, ipnt_taub, ival_bodsed1, & - ival_bodsedn, ipnt_bodsed1, ipnt_dpsed, ival_msed1, ival_msedn, ipnt_msed1, ival_lyrfrac1, ival_lyrfracn, & - ipnt_lyrfrac1, ipnt_poros, ipnt_thlyr, ival_frac1, ival_fracn, ipnt_frac1, ipnt_mudfrac, ipnt_sandfrac, & - ival_mfluff1, ival_mfluffn, ipnt_mfluff1, ival_fixfac1, ival_fixfacn, ipnt_fixfac1, ival_hidexp1, ival_hidexpn, & - ipnt_hidexp1, ival_sour1, ival_sourn, ipnt_sour1, ival_sink1, ival_sinkn, ipnt_sink1, ival_wqb1, ival_wqbn, & - ipnt_wqb1, ipnt_ucxq, ipnt_ucyq, ipnt_zcs, ipnt_ucx, ipnt_ucy, ipnt_ucxst, ipnt_ucyst, ipnt_ucz, ipnt_sa1, & - ipnt_tem1, ipnt_viu, ipnt_rhop, ipnt_rho, ipnt_umag, ipnt_qmag, ival_tra1, ival_tran, ipnt_tra1, ival_hwq1, & - ival_hwqn, ipnt_hwq1, ival_wqb3d1, ival_wqb3dn, ipnt_wqb3d1, ival_sf1, ival_sfn, ipnt_sf1, ival_ws1, ival_wsn, & - ipnt_ws1, ipnt_sed, ipnt_smx, smxobs, ipnt_zws, ipnt_vicwws, ipnt_difwws, ipnt_bruv, ipnt_richs, ival_seddif1, & - ival_seddifn, ipnt_seddif1, ipnt_zwu, ipnt_vicwwu, ipnt_tkin, ipnt_teps, ipnt_rich, ipnt_rain, ipnt_airdensity, & - ipnt_infiltcap, ipnt_infiltact, ipnt_wind, ipnt_tair, ipnt_rhum, ipnt_clou, ipnt_qsun, ipnt_qeva, ipnt_qcon, & - ipnt_qlon, ipnt_qfre, ipnt_qfrc, ipnt_qtot + use m_observations_data, only : valobs, numobs, nummovobs, kobs, lobs, ipnt_s1, ipnt_hs, ipnt_bl, ipnt_cmx, cmxobs, & + ipnt_wx, ipnt_wy, ipnt_patm, ipnt_waver, ipnt_waveh, ipnt_wavet, ipnt_waved, ipnt_wavel, ipnt_waveu, ipnt_taux, & + ipnt_tauy, ival_sbcx1, ival_sbcxn, ipnt_sbcx1, ival_sbcy1, ival_sbcyn, ipnt_sbcy1, ival_sscx1, ival_sscxn, & + ipnt_sscx1, ival_sscy1, ival_sscyn, ipnt_sscy1, ival_sbwx1, ival_sbwxn, ipnt_sbwx1, ival_sbwy1, ival_sbwyn, & + ipnt_sbwy1, ival_sswx1, ival_sswxn, ipnt_sswx1, ival_sswy1, ival_sswyn, ipnt_sswy1, ipnt_taub, ival_bodsed1, & + ival_bodsedn, ipnt_bodsed1, ipnt_dpsed, ival_msed1, ival_msedn, ipnt_msed1, ival_lyrfrac1, ival_lyrfracn, & + ipnt_lyrfrac1, ipnt_poros, ipnt_thlyr, ival_frac1, ival_fracn, ipnt_frac1, ipnt_mudfrac, ipnt_sandfrac, & + ival_mfluff1, ival_mfluffn, ipnt_mfluff1, ival_fixfac1, ival_fixfacn, ipnt_fixfac1, ival_hidexp1, ival_hidexpn, & + ipnt_hidexp1, ival_sour1, ival_sourn, ipnt_sour1, ival_sink1, ival_sinkn, ipnt_sink1, ival_wqb1, ival_wqbn, & + ipnt_wqb1, ipnt_ucxq, ipnt_ucyq, ipnt_zcs, ipnt_ucx, ipnt_ucy, ipnt_ucxst, ipnt_ucyst, ipnt_ucz, ipnt_sa1, & + ipnt_tem1, ipnt_viu, ipnt_rhop, ipnt_rho, ipnt_umag, ipnt_qmag, ival_tra1, ival_tran, ipnt_tra1, ival_hwq1, & + ival_hwqn, ipnt_hwq1, ival_wqb3d1, ival_wqb3dn, ipnt_wqb3d1, ival_sf1, ival_sfn, ipnt_sf1, ival_ws1, ival_wsn, & + ipnt_ws1, ipnt_sed, ipnt_smx, smxobs, ipnt_zws, ipnt_vicwws, ipnt_difwws, ipnt_bruv, ipnt_richs, ival_seddif1, & + ival_seddifn, ipnt_seddif1, ipnt_zwu, ipnt_vicwwu, ipnt_tkin, ipnt_teps, ipnt_rich, ipnt_rain, ipnt_airdensity, & + ipnt_infiltcap, ipnt_infiltact, ipnt_wind, ipnt_tair, ipnt_rhum, ipnt_clou, ipnt_qsun, ipnt_qeva, ipnt_qcon, & + ipnt_qlon, ipnt_qfre, ipnt_qfrc, ipnt_qtot, neighbour_nodes_obs,neighbour_weights_obs, intobs use m_sediment, only: jahissigwav, stm_included, stmpar, ustokes, hwav, twav, phiwav, rlabda, uorb, sedtra, fp, mtd, sed use Timers, only: timon, timstrt, timstop use m_gettaus, only: gettaus @@ -89,6 +89,7 @@ subroutine fill_valobs() use fm_statistical_output, only: model_is_3d use m_links_to_centers, only: links_to_centers use m_wind, only: wx, wy, jawind, air_pressure_available, air_pressure, jarain, rain, air_density, air_temperature, relative_humidity, cloudiness + use fm_location_types implicit none @@ -103,9 +104,13 @@ subroutine fill_valobs() real(kind=dp), allocatable :: poros(:) real(kind=dp), allocatable :: ueux(:) real(kind=dp), allocatable :: ueuy(:) + real(kind=dp), allocatable :: tmp_interp(:) real(kind=dp), allocatable :: vius(:) !< Flowlink-averaged horizontal viscosity (viu) at s-point - + kmx_const = kmx + if (kmx == 0) then + kmx_const = 1 ! to make numbering work + end if nlyrs = 0 if (timon) then @@ -116,6 +121,10 @@ subroutine fill_valobs() call realloc(ueux, ndkx, keepExisting=.false., fill=0.0_dp) call realloc(ueuy, ndkx, keepExisting=.false., fill=0.0_dp) end if + if (.not. allocated(tmp_interp)) then + ! Allocate as 2D arry for water levels etc. + call realloc(tmp_interp, ndx, keepExisting=.false., fill=0.0_dp) + end if ! if (jawave > NO_WAVES) then if (jahissigwav == 0) then @@ -202,6 +211,16 @@ subroutine fill_valobs() do i = 1, numobs + nummovobs k = max(kobs(i), 1) link_id_nearest = lobs(i) + if (intobs(i) == 0) then + ! Treat snapped stations as interpolated ones! + neighbour_nodes_obs(1,i) = k + neighbour_nodes_obs(2,i) = k + neighbour_nodes_obs(3,i) = k + neighbour_weights_obs(1,i) = 1.0 + neighbour_weights_obs(2,i) = 0.0 + neighbour_weights_obs(3,i) = 0.0 + end if + if (kobs(i) > 0) then ! rely on reduce_kobs to have selected the right global flow nodes if (model_is_3D()) then @@ -219,21 +238,88 @@ subroutine fill_valobs() call linkstocentercartcomp(k, ustokes, wa) ! wa now 2*1 value or 2*1 vertical slice end if -! store values in valobs work array valobs(i, :) = dmiss ! Intended to have dmiss on inactive layers for output. - ! It is taken care of in subroutine reduce_valobs for parallel computation. - valobs(i, IPNT_S1) = s1(k) - if (nshiptxy > 0) then + ! Fill valobs: Start with interpolateds values of hydrodynamic quantities needed for off-line nesting + ! (water levells, velocities, sality and temperature). Treat other quanitities (water quality, morpholgy, turbulence) as before (snapped) + ! + ! Water levels + + call interpolate_horizontal (s1,i,IPNT_S1,UNC_LOC_S) + + if (nshiptxy > 0) then if (allocated(zsp)) then - valobs(i, IPNT_S1) = valobs(i, IPNT_S1) + zsp(k) + tmp_interp = s1 + zsp + call interpolate_horizontal (tmp_interp,i,IPNT_S1,UNC_LOC_S) end if end if - valobs(i, IPNT_HS) = s1(k) - bl(k) - - valobs(i, IPNT_BL) = bl(k) + ! Water Depth + tmp_interp = s1 - bl + call interpolate_horizontal (tmp_interp,i,IPNT_HS,UNC_LOC_S) + ! Bed level + call interpolate_horizontal (bl ,i,IPNT_BL,UNC_LOC_S) valobs(i, IPNT_CMX) = cmxobs(i) + + ! For now here: interpolate velocities, salinity and temperature (not within loop from kb to ke, taken care of in interpolate horizontal) + ! First : allocate tmp_interp for 3D quantities + call realloc(tmp_interp, ndkx, keepExisting=.false., fill=0.0_dp) + + ! Horizontal velocities (3D) + if (jahisvelocity > 0 .or. jahisvelvec > 0) then + call interpolate_horizontal (ueux,i,IPNT_UCX,UNC_LOC_S3D) + call interpolate_horizontal (ueuy,i,IPNT_UCY,UNC_LOC_S3D) + end if + + ! Vertical velocities (3D) + if (model_is_3D()) then + call interpolate_horizontal (ucz,i,IPNT_UCZ,UNC_LOC_S3D) + end if + + ! Velocity magnitude (3D) + if (jahisvelocity > 0) then + call interpolate_horizontal (ucmag,i,IPNT_UMAG,UNC_LOC_S3D) + end if + + ! Depth averaged velocities (first ndx points of ucx/ucy array) + if (model_is_3D()) then + call interpolate_horizontal (ucx,i,IPNT_UCXQ,UNC_LOC_S) + call interpolate_horizontal (ucy,i,IPNT_UCYQ,UNC_LOC_S) + end if + + ! Salinity (interpolated) + if (jasal > 0) then + tmp_interp = constituents(isalt,:) + call interpolate_horizontal (tmp_interp,i,IPNT_SA1,UNC_LOC_S3D) + end if + + ! Temperature + ! if (jatem > 0) then + if (temperature_model /= TEMPERATURE_MODEL_NONE) then + tmp_interp = constituents(itemp,:) + call interpolate_horizontal (tmp_interp,i,IPNT_TEM1,UNC_LOC_S3D) + end if + + ! Finally; vertical positions + ! Maybe not the right place to do this, fille interfaces with bed level and water surface in case of 2d model + + + if (model_is_3D()) then + ! interface + call interpolate_horizontal (zws,i,IPNT_ZWS,UNC_LOC_W) + ! centre: make temporary array with cellcentres + do j = 2, ndkx + tmp_interp(j) = 0.5_dp * (zws(j) + zws(j - 1)) + end do + call interpolate_horizontal (tmp_interp,i,IPNT_ZCS,UNC_LOC_S3D) + else + ! TK_Temp: Fill interfaces with surface and bed, centre with average (woulde be nicer to fill zws with correct values) + valobs(i,IPNT_ZWS) = valobs(i,IPNT_BL) + valobs(i,IPNT_ZWS + 1) = valobs(i,IPNT_S1) + valobs(i,IPNT_ZCS) = 0.5_dp * (valobs(i,IPNT_BL) + valobs(i,IPNT_S1)) + end if + + ! Frome here: everything as snapped!!! if (jawind > 0) then valobs(i, IPNT_wx) = 0.0_dp valobs(i, IPNT_wy) = 0.0_dp @@ -399,22 +485,24 @@ subroutine fill_valobs() end do end if - if (model_is_3D()) then - valobs(i, IPNT_UCXQ) = ucx(k) - valobs(i, IPNT_UCYQ) = ucy(k) - end if + ! Taken care off by interpolate_horizontal +! if (model_is_3D()) then +! valobs(i, IPNT_UCXQ) = ucx(k) +! valobs(i, IPNT_UCYQ) = ucy(k) +! end if do kk = kb, kt klay = kk - kb + nlayb + ! Taken care off by interpolate_horizontal +! if (model_is_3D()) then +! valobs(i, IPNT_ZCS + klay - 1) = 0.5_dp * (zws(kk) + zws(kk - 1)) +! end if - if (model_is_3D()) then - valobs(i, IPNT_ZCS + klay - 1) = 0.5_dp * (zws(kk) + zws(kk - 1)) - end if - - if (jahisvelocity > 0 .or. jahisvelvec > 0) then - valobs(i, IPNT_UCX + klay - 1) = ueux(kk) - valobs(i, IPNT_UCY + klay - 1) = ueuy(kk) - end if + ! Taken care off by interpolate_horizontal +! if (jahisvelocity > 0 .or. jahisvelvec > 0) then +! valobs(i, IPNT_UCX + klay - 1) = ueux(kk) +! valobs(i, IPNT_UCY + klay - 1) = ueuy(kk) +! end if if (jawave > NO_WAVES .and. .not. flow_without_waves) then if (hs(k) > epshu) then @@ -431,12 +519,12 @@ subroutine fill_valobs() if (model_is_3D()) then valobs(i, IPNT_UCZ + klay - 1) = ucz(kk) end if - if (jasal > 0) then - valobs(i, IPNT_SA1 + klay - 1) = constituents(isalt, kk) - end if - if (temperature_model /= TEMPERATURE_MODEL_NONE) then - valobs(i, IPNT_TEM1 + klay - 1) = constituents(itemp, kk) - end if +! if (jasal > 0) then +! valobs(i, IPNT_SA1 + klay - 1) = constituents(isalt, kk) +! end if +! if (jatem > 0) then +! valobs(i, IPNT_TEM1 + klay - 1) = constituents(itemp, kk) +! end if if (jahistur > 0) then valobs(i, IPNT_VIU + klay - 1) = vius(kk) end if @@ -446,9 +534,11 @@ subroutine fill_valobs() valobs(i, IPNT_RHO + klay - 1) = in_situ_density(kk) end if end if - if (jahisvelocity > 0) then - valobs(i, IPNT_UMAG + klay - 1) = ucmag(kk) - end if + + ! Taken care of by interpolate_horizontal + !if (jahisvelocity > 0) then + ! valobs(i, IPNT_UMAG + klay - 1) = ucmag(kk) + !end if valobs(i, IPNT_QMAG + klay - 1) = 0.5_dp * (squ(kk) + sqi(kk)) if (kmx == 0) then @@ -502,7 +592,8 @@ subroutine fill_valobs() call getlayerindices(k, nlayb, nrlay) do kk = kb - 1, kt klay = kk - kb + nlayb + 1 - valobs(i, IPNT_ZWS + klay - 1) = zws(kk) + ! Taken care of by interpolate_horizontal +! valobs(i, IPNT_ZWS + klay - 1) = zws(kk) if (iturbulencemodel >= 2 .and. jahistur > 0) then valobs(i, IPNT_VICWWS + klay - 1) = vicwws(kk) valobs(i, IPNT_DIFWWS + klay - 1) = difwws(kk) @@ -653,5 +744,66 @@ subroutine conditional_assign(valobs, i, ipnt, array, k) valobs(i, ipnt) = real(array(k), dp) end if end subroutine conditional_assign + + subroutine interpolate_horizontal (rarray,istat,IPNT,locType) + + ! Interpolate (horizontally, within a computational layer) to a position from 3 surrounding snapped points + ! rarray can be constituents (3D) or water levels (2D) + ! + use precision, only: dp + use fm_statistical_output, only: model_is_3d + use m_observations_data + use m_get_kbot_ktop + use m_get_layer_indices + use fm_location_types + + integer , intent(in) :: istat, IPNT,locType + real(kind=dp), intent(in), allocatable :: rarray (:) + + real(kind=dp) :: value + real(kind=dp) :: weighttot + + integer :: kb_tmp(3), kt_tmp(3), nlayb_tmp(3), nrlay_tmp(3), iwght, kstart, kstop, pntnr, klay, oneDown + + oneDown = 0 + + do iwght = 1, 3 + if (model_is_3D() .and. (locType == UNC_LOC_S3D .or. LocType == UNC_LOC_W)) then + call getkbotktop (neighbour_nodes_obs(iwght,istat), kb_tmp(iwght), kt_tmp(iwght)) + call getlayerindices(neighbour_nodes_obs(iwght,istat), nlayb_tmp(iwght), nrlay_tmp(iwght)) + + else + kb_tmp (iwght) = neighbour_nodes_obs(iwght,istat) + kt_tmp (iwght) = neighbour_nodes_obs(iwght,istat) + nlayb_tmp(iwght) = 1 + nrlay_tmp(iwght) = 1 + end if + end do + + ! Values ar interfaces stored 1 below (interface 1 effectively corresponds with layer 0) + if (LocType == UNC_LOC_W) then + nrlay_tmp = nrlay_tmp + 1 + oneDown = 1 + end if + + ! Determine start and stop layer nr for interpolated values + kstart = minval(nlayb_tmp) + kstop = maxval(nlayb_tmp + nrlay_tmp - 1) + + do klay = kstart, kstop + + value = 0.0_dp + weighttot = 0.0_dp + + do iwght = 1, 3 + if ((klay >= nlayb_tmp(iwght)) .and. (klay <= nlayb_tmp(iwght) + nrlay_tmp(iwght) - 1)) then + pntnr = kb_tmp(iwght) - nlayb_tmp(iwght) + klay - oneDown + value = value + rarray(pntnr)*neighbour_weights_obs(iwght,istat) + weighttot = weighttot + neighbour_weights_obs(iwght,istat) + end if + end do + valobs(istat, IPNT + klay - 1) = value/weighttot + end do + end subroutine interpolate_horizontal end module m_fill_valobs diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_obs_on_flowgeom.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_obs_on_flowgeom.f90 index 72c22d579e..9f7d672af1 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_obs_on_flowgeom.f90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_obs_on_flowgeom.f90 @@ -41,7 +41,7 @@ module m_obs_on_flowgeom !! Results are stored in `m_observations` state arrays. subroutine obs_on_flowgeom(iobstype) - use m_observations_data, only: numobs, nummovobs, kobs, namobs + use m_observations_data, only: numobs, nummovobs, kobs, namobs, neighbour_nodes_obs, neighbour_weights_obs use unstruc_messages, only: loglevel_StdOut use messagehandling, only: LEVEL_DEBUG, LEVEL_INFO, msgbuf, mess use m_flowgeom, only: ndx2D, ndxi @@ -83,7 +83,11 @@ subroutine obs_on_flowgeom(iobstype) end if else ! No cache, so process all requested observation points. - call find_flownodes_and_links_for_all_observation_stations(n1, n2) + if (n2 - n1 >= 0) then + call find_flownodes_and_links_for_all_observation_stations(n1, n2) + + call init_interpolation_data_for_all_observation_stations(n1, n2, neighbour_nodes_obs, neighbour_weights_obs) + end if end if if (loglevel_StdOut == LEVEL_DEBUG) then @@ -255,4 +259,54 @@ subroutine find_flownodes_and_links_for_all_observation_stations(nstart, nend) end subroutine find_flownodes_and_links_for_all_observation_stations + subroutine init_interpolation_data_for_all_observation_stations(n_start, n_end,neighbour_nodes_obs,neighbour_weights_obs) + + use m_observations_data , only: xobs, yobs, numobs, nummovobs + use m_flowgeom , only: xz, yz, ndx2d + use m_missing , only: dmiss + use m_sferic , only: jsferic, jasfer3D + use fm_external_forcings_data, only: transformcoef + use m_polygon , only: npl, xpl, ypl, zpl + use m_samples , only: mxsam, mysam + use precision , only: dp + use m_alloc + + use m_ec_basic_interpolation, only: triinterp2 + ! use m_ec_triangle, only: jagetwf + use m_ec_triangle, only: jagetwf, indxx, wfxx + + + integer , intent(in) :: n_start !< Starting index of obs + integer , intent(in) :: n_end !< Ending index of obs + integer , dimension(3,n_start:n_end),intent(inout) :: neighbour_nodes_obs ! Table of nearby flow node numbers for each station + real(kind=dp), dimension(3,n_start:n_end),intent(inout) :: neighbour_weights_obs ! Table of interpolation weights for nearby flow node numbers for each station + + integer :: jdla, i, jagetwf_org + real(kind=dp), allocatable :: dummyZ (:) + real(kind=dp), allocatable :: dumout (:) + + ! Create arrays with dummy z values(samples and output, needto deallocate?) + call realloc(dummyZ,ndx2d , keepexisting=.false., fill=dmiss) + call realloc(dumout,numobs + nummovobs, keepexisting=.false., fill=dmiss) + + ! interpolate + jagetwf_org = jagetwf + jdla = 1 + jagetwf = 1 + + call realloc(indxx, [3, numobs + nummovobs], keepexisting=.false., fill=0) + call realloc(wfxx, [3, numobs + nummovobs], keepexisting=.false., fill=0.0_dp) + + call triinterp2(xobs, yobs,dumout, numobs + nummovobs, jdla ,xz(1:ndx2d), yz(1:ndx2d), dummyZ, ndx2d, dmiss, jsferic, 1 , & + jasfer3D, NPL, MXSAM, MYSAM, XPL, YPL, ZPL, transformcoef) + + do i = 1, numobs + nummovobs + neighbour_nodes_obs (:, i) = indxx(:, i) + neighbour_weights_obs(:, i) = wfxx (:, i) + end do + + jagetwf = jagetwf_org + + end subroutine init_interpolation_data_for_all_observation_stations + end module m_obs_on_flowgeom diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_observations.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_observations.f90 index c4f2214a11..06fd31fb9c 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_observations.f90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_observations.f90 @@ -361,7 +361,8 @@ subroutine init_valobs_pointers() if (jased > 0 .and. .not. stm_included) then IVAL_SED = next_index(i) end if - if (kmx > 0) then + !TK_Temp: Also for KMX = 0 + if (kmx >= 0) then IVAL_ZCS = next_index(i) end if if (use_density()) then @@ -377,9 +378,12 @@ subroutine init_valobs_pointers() ! 3D, layer interfaces i0 = i - if (kmx > 0) then + ! TK_Temp: also for kmx == 0 + if (kmx >= 0) then IVAL_ZWS = next_index(i) IVAL_ZWU = next_index(i) + end if + if (kmx > 0) then IVAL_BRUV = next_index(i) if (iturbulencemodel > 0 .and. jahistur > 0) then IVAL_TKIN = next_index(i) @@ -735,11 +739,14 @@ subroutine addObservation(x, y, name, isMoving, loctype, iOP) call realloc(xyobs, 2 * (nummovobs + capacity_)) call realloc(kobs, numobs + nummovobs + capacity_) call realloc(lobs, numobs + nummovobs + capacity_) + call realloc(neighbour_nodes_obs, [3, numobs + nummovobs + capacity_]) + call realloc(neighbour_weights_obs, [3, numobs + nummovobs + capacity_]) call realloc(namobs, numobs + nummovobs + capacity_) call realloc(smxobs, numobs + nummovobs + capacity_) call realloc(cmxobs, numobs + nummovobs + capacity_) call realloc(locTpObs, numobs + nummovobs + capacity_) call realloc(obs2OP, numobs + nummovobs + capacity_) + call realloc(intobs, numobs + nummovobs + capacity_) end if ! Before adding new normal observation station: @@ -755,6 +762,7 @@ subroutine addObservation(x, y, name, isMoving, loctype, iOP) cmxobs(i + 1) = cmxobs(i) locTpObs(i + 1) = locTpObs(i) obs2OP(i + 1) = obs2OP(i) + intobs(i + 1) = intobs(i) end do numobs = numobs + 1 inew = numobs @@ -766,6 +774,7 @@ subroutine addObservation(x, y, name, isMoving, loctype, iOP) ! Add the actual station (moving or static) xobs(inew) = x yobs(inew) = y + intobs(inew) = 0 namobs(inew) = name_ kobs(inew) = -999 ! Cell number is set elsewhere lobs(inew) = -999 ! Flow link number is set elsewhere @@ -880,6 +889,7 @@ subroutine purgeObservations() k = k + 1 xobs(k) = xobs(i) yobs(k) = yobs(i) + intobs(k) = intobs(i) kobs(k) = kobs(i) lobs(k) = lobs(i) namobs(k) = namobs(i) @@ -910,6 +920,9 @@ subroutine deleteObservations() deallocate (cmxobs) deallocate (locTpObs) deallocate (obs2OP) + deallocate (intobs) + deallocate (neighbour_nodes_obs) + deallocate (neighbour_weights_obs) end if call dealloc(network%obs) ! deallocate obs (defined in *.ini file) @@ -924,6 +937,9 @@ subroutine deleteObservations() allocate (cmxobs(capacity_)) allocate (locTpObs(capacity_)) allocate (obs2OP(capacity_)) + allocate (intobs(capacity_)) + allocate (neighbour_nodes_obs(3, capacity_)) + allocate (neighbour_weights_obs(3, capacity_)) kobs = -999 lobs = -999 @@ -957,12 +973,15 @@ subroutine loadObservations(filename, jadoorladen) tok = index(filename, '.xy') if (tok > 0) then call loadObservations_from_xyn(filename) - else - tok = index(filename, '.ini') - if (tok > 0) then - call readObservationPoints(network, filename) - call addObservation_from_ini(network, filename) - end if + end if + tok = index(filename, '.pli') + if (tok > 0) then + call loadObservations_from_pli(filename) + end if + tok = index(filename, '.ini') + if (tok > 0) then + call readObservationPoints(network, filename) + call addObservation_from_ini(network, filename) end if else call mess(LEVEL_ERROR, "Observation file '"//trim(filename)//"' not found!") @@ -1032,4 +1051,28 @@ subroutine saveObservations(filename) end subroutine saveObservations + subroutine loadObservations_from_pli(filename) + + use m_filez, only: oldfil + use m_polygon + use m_reapol_nampli, only: reapol_nampli + + + implicit none + character(len=*), intent(in) :: filename + + ! locals + integer :: minp,istat, ipli + character(5) :: numstr + + call oldfil(minp, filename) + ipli = 0 + call reapol_nampli(minp, 0, 1, ipli) + + do istat = 1, npl + write(numstr,'(i4.4)') istat + call addObservation(xpl(istat), ypl(istat), trim(nampli(1))//'_'//numstr) + intobs(numobs) = 1 + end do + end subroutine loadObservations_from_pli end module m_observations diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_observations_data.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_observations_data.f90 index 39b6041815..ff8aa90d9b 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_observations_data.f90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_observations_data.f90 @@ -45,6 +45,12 @@ module m_observations_data real(kind=dp), allocatable :: cmxobs(:) !< maximum 2D flow velocity of observation points, 3D: maximum over all layers and time integer, allocatable :: kobs(:) !< node nrs of ACTIVE observation points integer, allocatable :: lobs(:) !< flowlink nrs of active observation points + integer, allocatable :: intobs(:) !< interpolated station or not +! For storing number and weights of stations wher you want to get interpolated values + + integer, allocatable :: neighbour_nodes_obs(:,:) !< [3, numobs+nummovobs] List of nearby flow node numbers for each observation point + real(kind=dp), allocatable :: neighbour_weights_obs(:,:) !< [3, numobs+nummovobs] List of weights for the nearby flow node numbers for each observation point + ! NOTE: kobs is not maintained here (so also not after deleteObservation, etc.) All done once by obs_on_flowgrid. character(len=IdLen), allocatable :: namobs(:) ! names of observation points integer, allocatable :: locTpObs(:) !< location type of observation points, determining to which flownodes to snap to (0=1d2d, 1=1d, 2=2d, 3=1d defined by branchID+chainage)